In this work, the reconstruction of sub-grid scales in large eddy simulation (LES) of turbulent flows in stratocumulus clouds is addressed. The approach is based on the fractality assumption of turbulent velocity field. The fractal model reconstructs sub-grid velocity field from known filtered values on LES grid, by means of fractal interpolation, proposed by Scotti and Meneveau (Physica D 127, 198–232 1999). The characteristics of the reconstructed signal depend on the stretching parameter d, which is related to the fractal dimension of the signal. In many previous studies, the stretching parameter values were assumed to be constant in space and time. To improve the fractal interpolation approach, we account for the stretching parameter variability. The local stretching parameter is calculated from direct numerical simulation (DNS) data with an algorithm proposed by Mazel and Hayes (IEEE Trans. Signal Process 40(7), 1724–1734, 1992), and its probability density function (PDF) is determined. It is found that the PDFs of d have a universal form when the velocity field is filtered to wave-numbers within the inertial range. In order to investigate Reynolds number (Re) dependence, we compare the inertial-range PDFs of d in DNS and large eddy simulation (LES) of stratocumulus cloud-top and experimental airborne data from Physics of Stratocumulus Top (POST) research campaign. Next, fractal reconstruction of the subgrid velocity is performed and energy spectra and statistics of velocity increments are compared with DNS data. It is assumed that the stretching parameter d is a random variable with the prescribed PDF. We show that the agreement with the DNS is in such case better and the error in mass conservation is smaller compared to the use of constant values of d. The motivation of this study is to reproduce effect of sub-grid scales on a motion of Lagrangian particles (e.g. droplets) in clouds.