A flexible three-dimensional stratocumulus, cumulus and cirrus cloud generator (3DCLOUD) based on drastically simplified atmospheric equations and the Fourier transform framework

The 3DCLOUD algorithm for generating stochastic three-dimensional (3-D) cloud fields is described in this paper. The generated outputs are 3-D optical depth ( τ ) for stratocumulus and cumulus fields and 3-D ice water content (IWC) for cirrus clouds. This model is designed to generate cloud fields that share some statistical properties observed in real clouds such as the inhomogeneity parameter ρ (standard deviation normalized by the mean of the studied quantity), the Fourier spectral slope β close to−5/3 between the smallest scale of the simulation to the outer Lout (where the spectrum becomes flat). Firstly, 3DCLOUD assimilates meteorological profiles (humidity, pressure, temperature and wind velocity). The cloud coverage C, defined by the user, can also be assimilated, but only for stratocumulus and cumulus regime. 3DCLOUD solves drastically simplified basic atmospheric equations, in order to simulate 3-D cloud structures of liquid or ice water content. Secondly, the Fourier filtering method is used to constrain the intensity of ρ, β, Lout and the mean ofτ or IWC of these 3-D cloud structures. The 3DCLOUD model was developed to run on a personal computer under Matlab environment with the Matlab statistics toolbox. It is used to study 3-D interactions between cloudy atmosphere and radiation.


Introduction
Clouds have a significant effect on the Earth radiation budget.They reflect the solar radiation and reduce the warming of the Earth (albedo effect).They also create a greenhouse effect by trapping the thermal radiation emitted from the Earth's surface, reducing the radiative cooling of the Earth (Collins and Satoh, 2009).Cloud feedback has remained, however, the largest uncertainty in the study of climate sensitivity for almost 20 years (Bony et al., 2006).In almost all climate models, clouds are assumed plane-parallel with homogeneous optical properties (PPH), and radiation codes use a one-dimensional (1-D) scheme.Therefore, improving parameterisations of clouds in large-scale model, especially their interaction with radiation, is a challenge in order to reduce uncertainty in model projections of the future climate (Illingworth and Bony, 2009).Improved global characterization of the three dimensional (3-D) spatial distribution of clouds is, thus, necessary (Clothiaux et al., 2004).
Moreover, satellite passive sensors, such as multi-spectral and multi-angular radiometers, and satellite active sensor, such as LIDAR and RADAR in the A-train mission, allow the retrieval of cloud horizontal and vertical optical properties with an adequate spatial and temporal coverage.For practical and computational cost purposes, interpretation of such measurements generally also assumes 1-D radiative algorithm and PPH cloud.This assumption can be far from being realistic and leads to biases on the retrieved properties from passive sensors (Barker and Liu, 1995;Várnai and F. Szczap et al.: A flexible three-dimensional cloud generator (3DCLOUD) Marshak, 2002Marshak, , 2007;;Lafont and Guillemet, 2004;Cornet et al., 2005Cornet et al., , 2013) ) and active sensors (Battaglia and Tanelli, 2011).These biases depend at least, on the cloud coverage and on the variability of cloud optical depth or water content.This variability is quantified by an inhomogeneity parameter, often defined as the standard deviation normalized by the mean of the studied quantity (Szczap et al., 2000;Carlin et al., 2002;Oreopoulos and Cahalan, 2005;Sassen et al., 2007;Hill et al., 2012).
Determining the significance of the 3-D inhomogeneity of clouds for climate and remote sensing applications requires the measurement and the simulation of the full range of actual cloud structure.Apart from the computational time, accurate 3-D cloudy radiative transfer problem is not an issue, per se (Evans and Wiscombe, 2004).Monte Carlo transfer models can indeed accurately and efficiently compute radiative properties for arbitrary cloud fields (Battaglia and Mantovani, 2005;Pincus and Evans, 2009;Mayer, 2009;Cornet et al., 2010Cornet et al., , 2013;;Battaglia and Tanelli, 2011;Fauchez et al., 2013).The difficulty is to generate cloud property fields that are statistically representative of cloud fields in nature.
Cloud fields generated by dynamic cloud models, such as the cloud resolving model (CRM) or the large-eddy simulation model (LES), are very attractive, as they contain the state of the art of physical processes (resolution of atmospheric equations, detailed microphysics, radiation, etc.).The goal of the LES approach is to simulate the three-dimensional atmospheric turbulent flows.There are different scales of turbulent eddies; large eddies (from 100 to 1000 m and more) that are produced directly by the instability of the mean flow and small eddies (from a few centimetres to 100 m) as well as by the energy-cascade process from the larger eddies (Moeng, 1984).LES seeks to capture accurately the larger eddies, while only modelling the smaller ones.Instead of reproducing all the scales of turbulence flow, they can integrate a flow in which small scale details are removed from the solution.The spatial filtered equations can, therefore, be integrated with available resources (Bryan et al., 2003).Nevertheless, they are still very expensive to run in a 3-D domain.
Stochastic models have the capability to simulate quickly realistic 2-D and 3-D cloud structures with just a few parameters.Examples of these types of cloud models are: the bounded cascade model (Cahalan et al., 1994;Marshak et al., 1998), the iterative amplitude adapted Fourier transform (IAAFT) algorithm (Venema et al., 2006), the SITCOM model (Di Guiseppe and Thompkins, 2003), the tree-driven mass accumulation process (tdMAP) model (Benassi et al., 2004), the model developed by Evans and Wiscombe (2004) for low liquid clouds (stratocumulus and cumulus) or by Alexandrov et al. (2010) and the Cloudgen model (Hogan and Kew, 2005) for high ice clouds (cirrus).These stochastic models are based on fractal or Fourier framework.The scale invariant properties observed in real clouds can be controlled.The power spectra of the logarithm of their optical properties (optical depth, liquid water content or liquid water path for low clouds and ice water content for high clouds) typically exhibits a spectral slope of around −5/3 (Davis et al., 1994(Davis et al., , 1996(Davis et al., , 1997(Davis et al., , 1999;;Cahalan et al., 1994;Benassi et al., 2004;Hogan and Kew, 2005;Hill et al., 2012;Fauchez et al., 2014) from small scale (a few metres) to the "integral scale" or the outer scale (few tenths of a kilometre to one-hundred kilometres), where the spectrum becomes flat (i.e.decorrelation occurs).The disadvantage of such models arises from the fact that effects of meteorological processes are not always considered and dominant scales of organization related to turbulent eddy due, for example, to wind shear, convection, and entrainment are not directly modelled.At the same time, it should be noted that Cloudgen does consider the effect of wind shear on cirrus cloud.
The aim of the 3DCLOUD algorithm is to reconcile these two approaches.In Sect.2, we describe the 3DCLOUD generator.In Sect.3, 3DCLOUD outputs are compared to LES outputs to check the validity of the chosen basic atmospheric equations.In Sect.4, stratocumulus, cumulus and cirrus examples provided by 3DCLOUD are presented.
2 The 3DCloud generator 3DCLOUD generates, in two distinct steps (see Fig. 1), a 3-D optical depth field for stratocumulus and cumulus or a 3-D ice water content field for cirrus clouds.These cloud fields were chosen as most of the papers dealing with scale invariant properties focus on liquid water path and optical depth for stratocumulus and cumulus and on ice water content for cirrus.During the first step, meteorological vertical profiles (temperature, pressure, wind, humidity), defined by the user, are assimilated and basic atmospheric equations are resolved.During the second step, cloud scale invariant properties are constrained in a Fourier framework.At the same time, a gamma distribution of local optical depth or 3-D ice water content (IWC) is mapped onto the liquid water content (LWC), or IWC, generated during the first step.This gamma distribution is iteratively computed in such a way that the mean optical depth or IWC and the inhomogeneity parameter satisfy the values imposed by the user.Details of these two steps are presented below.

Step 1: the 3-D LWC/IWC generator
The essential basic quantities to generate cloud fields are the condensed water mixing ratio q c = q l +q i where q l is the liquid water mixing ratio and q i is the ice water mixing ratio, the wind velocity vector u, air pressure p, temperature T , and vapour water mixing ratio q v .Mixing ratios are the mass of vapour or condensed water per unit of dry air mass.We describe in this section the equations used to generate clouds with the associated simplifications used.

The simplification of basic atmospheric equations
The continuity and momentum equations of the atmosphere can be written as follows (Houze, 1993): where t is time, ρ the air density, f the Coriolis parameter, g the acceleration due to gravity and F the acceleration due to other forces (frictional acceleration for example).D dt = ∂ ∂t + u • ∇ is the Lagrangian derivative operator following a parcel of air, ∂ ∂t is the Eulerian derivative operator and ∇ is the three-dimensional gradient operator.u = ui + vj + wk is the wind velocity vector with horizontal components u, v and vertical component w projected in the Cartesian geometry system, where i, j and k are the unit vectors in the x, y and z directions.The continuity and momentum equations of the atmosphere under the anelastic and Boussinesq approximation, assuming shallow motion, neglecting Coriolis parameter, neglecting frictional forces, and neglecting the molecular viscosity can be written as follows (Holton, 2004, p. 117;Houze, 1993, p. 35;Emanuel, 1994, p. 11): where B is the buoyancy acceleration, ρ 0 is the constant mean value of air density and p * is the pressure perturbations.The above differential operators are valid only in the limit when δt, δx, δy and δz approach 0 (Pielke, 2002, p. 41).Nevertheless, turbulent motions (shear induced eddies, convection eddies) have spatial and temporal variations at scales much smaller than those resolved by LES and 3DCLOUD.
If we assume field variables can be separated in slowly varying mean field and rapidly varying turbulent component, and if we apply the Reynolds decomposition, we can rewrite the above equation set as where is the three dimensional convergence of the eddy flux of moment (Houze, 1993, p. 42), the turbulent flux (Holton, 2002, p. 119) or the sub-grid correlation term (Pielke, 2002, p. 44).The Reynolds decomposition is not used in LES.The atmospheric equations are derived by spatial filtering, where a special function is applied.Thus, the filtering operation acts on atmospheric quantities and separates them in two categories: the resolved one (large eddy) and unresolved one (subgrid scale).An unknown term remains in the filtered equations of LES, often called the subgrid-scale stress, which needs to be parameterized or estimated with the help of subgrid-scale modelling.This subgrid-scale stress for LES equations is analogous to the term for Reynolds decomposition.In 3DCLOUD, the term is voluntarily neglected.Indeed, the guiding idea of 3DCLOUD is to simulate, in the fastest way, 3-D fluctuations of LWC/IWC of a cloud showing turbulent properties (or invariant scale properties).
When water phase changes are only associated with condensation and evaporation (or sublimation), the first law of thermodynamics can be written (Houze, 1993) as follows: www.where L = 2500 kJ kg −1 and L = 2800 kJ kg −1 are the usual latent heat of vaporization of water and ice, respectively.C p = 1.004 kJ kg −1 K −1 is the usual specific heat of dry air at constant pressure, θ is the potential temperature and = p p 0 κ = T θ is the Exner function where p 0 = 1000 hPa and κ = 0.286.In addition to the equation of motion and the first law of thermodynamics, air parcels follow the water continuity equation: where S i are the sum of the sources and sinks for a particular category (among n categories) of water indicated by i (vapour, solid, liquid water category for example).
As the horizontal extension of the simulated cloud fields is around a few km, horizontal pressure is assumed to be constant.Therefore, the current version of 3DCLOUD does not have a large enough domain to contain power in the mesoscale.All these considerations lead to a dramatic simplification of the dynamic equations.The simplified equations of 3DCLOUD governing the formation of 3-D cloud structures are where the reference state is denoted by subscript 0 and the deviation from the reference state by an asterisk, θ v = θ (1 + 0.61q v ) is the virtual potential temperature.For stratocumulus and cumulus fields ξ is estimated as follows: ξ = min (q vs − q v , q c ) t (7), where t is the simulation time step and q vs (T , p) is the saturation mixing ratio derived from Thetens and Magnus formula: q vs (T , p) = 0.622P sat (p/100−0.378Psat ) , where P sat = 6.107 exp 4028(T −273.15)234.82(T −38.33) for water, P sat = 6.107 × 10 9.5(T −273.15)265.5+(T −273.15) for ice.Computation of ξ at each simulation step is based on the work of Asai (1965).For cirrus clouds, condensation, evaporation and ice crystals sedimentation processes are very complex and still not well understood (Kärcher and Spichtinger, 2009).In order to take into account super-saturation and subsaturation regions in cirrus clouds in a very simple way, we used the parameterisation of Starr and Cox (1985) to compute the values of ξ every 2.5 min.Sedimentation processes are taken into account in Eq. (5.5) by adding ice fall speed v fall taken from Starr and Cox (1985): where v fall is in m s −1 and the ice water content IWC in g m −3 .

Assimilation of meteorological profiles and cloud coverage
In order to control the structure and nature of clouds and especially vertical position and extension, it is necessary to impose a large-scale environment.Practically, forcing terms are added to the 3DCLOUD equations to nudge the solutions towards observations.Our state observations are the initial meteorological profiles (provided by the user for example) and do not change during the simulation.The technique used is based on the initialization integration method (Pielke, 2002).Consequently, 3DCLOUD equations become where for variables X, X (z) is the mean of X at height z and quantities G X (z) are adjusted during the simulation in such a way that X or X (z) do not diverge far from the initial conditions X ini (z).In a general way, G is the inverse of timescale but, because the contribution of G is artificial, it must not be a dominant term in the governing equations and should be scaled by the slowest physical adjustment processes in the model (Cheng et al., 2001).This timescale was first set to 1 h but this value was found to be too large and must be adjusted as a function of altitude, especially at height where large vertical gradients of X appear (e.g. at the top of a stratocumulus cloud, for example).Therefore, we developed a very fast and simple numerical method to adjust the values of G X (z) during the simulation.At each level, we compute the relative difference α X (z) = X ini (z)−X(z) X(z) • 100.G X (z) is assumed to be proportional to α X (z) and is estimated as Values of α X,max were estimated during our numerous numerical experiments.For stratocumulus and cumulus, α X,max values for horizontal wind, temperature and humidity are 20 %, 2 % and 20 % respectively and for cirrus, α X,max values are 20 %, 10 % and 10 %, respectively.
The cloud coverage C is defined as the fraction of the number of cloudy pixels to the total number of pixels in the 2-D horizontal plan.The value of C is chosen with the assimilation of initial meteorological profiles.At each time step, the initial profile of vapour mixing ratio q v ini (z) is modified between cloud base and top height if C ≥ 50 % or between ground and cloud top height if C < 50 % until C agrees with the desired value within few percents.The new "initial" profile of vapour mixing ratio q new v ini (z) is computed from the currently simulated (old) profile of vapour mixing ratio q old v ini (z) as q new v ini (z) = q old v ini (z) ± n z −n base n top −n base q old v ini (z) × 0.1 100 where z is height, and n z , n top and n base are the levels indexes (in z direction) corresponding to cloud top height and to cloud base height (or ground), respectively.This method gives satisfactory results for stratocumulus and cumulus cloud fields (see Sect. 4.1.2),but not for cirrus fields.This is because condensation/evaporation and dynamic processes are different for stratocumulus/cumulus and cirrus regimes.Indeed, for liquid and warm stratocumulus/cumulus regime, liquid super or sub-saturation regions are not allowed in 3DCLOUD.Therefore, the distinction between cloudy and free cloud voxels is sharp.Moreover, as stratocumulus/cumulus fields are often driven by convection processes in a well-mixed planetary boundary layer, vertical correlation occurs between cloudy voxels (free cloud voxels) and updrafts (downdraughts).Thus, the fractional cloud coverage is easily controlled by adjusting the vertical profile of vapour mixing ratio during the simulation.By contrast, in ice cirrus regimes, (large) ice crystals can survive even if ice relative humidity is less than 100 %.Ice super or subsaturation regions are often observed in cirrus and are taken into account in the Starr and Cox parameterization used in 3DCLOUD.Therefore, many cloudy voxels still exist in our cirrus simulations, even if the ice water content is very small.The distinction between cloudy and free cloud voxels is, thus, very tenuous.Moreover, cirrus dynamics are often driven by wind shear; small fractional cloud coverage can exist at the top of the cirrus field due to convection or radiative cooling coexisting with large fractional cloud coverage and can also exist at the bottom of cirrus field due to wind shear.Finally, the total cloud coverage could be large.If we adjust the vertical profile of vapour mixing ratio during the simulation in the same way as for the stratocumulus/cumulus field, the total cloud coverage will be difficult to control.Further investigations are thus needed to perfectly control the cloud coverage of cirrus simulated by 3DCLOUD.

Implementation of 3DCLOUD algorithm
To implement the previously described equations, space is divided in (N x + 2) × N y + 2 × (N z + 2) cells or voxels where N x , N y and N z are the voxel numbers in each direction.A voxel is characterized by its spatial resolution with x = y = z.Horizontal extensions are L x = L y and can be different from the vertical extension L z .In order to take into account the boundary conditions, one layer of voxel is added around the simulation domain.
A semi-Lagrangian scheme was chosen to solve the equation: where X is a scalar advected by the wind velocity u.X can be the potential temperature θ, the condensed water q c or the vapour mixing ratios q v , and also the three components of wind velocity u, v and w.Two steps are needed in order to compute the value of X (x, t + t) at a fixed position x and at time t + t.X (x, t) and u (x, t) are known values and t is the time step.First, we compute the previous position p (X, t − t) = x − u (x, t) t of X at time t − t.In a second step, we compute the value of X (p, t) at the position p and at the time t by an interpolation scheme.This interpolated value X (p, t) is the desired value X (x, t + t).
The main advantage of this approach is that the time step is not restricted by the Courant-Friedrichs-Lewy stability limit, but by the less restrictive condition that parcels do not overtake each other during the time step (Riddaway, 2001).Therefore, at each iteration, the maximum value of time step t max can be roughly estimated as t max = 1 max u x + max v y + max w z .The accuracy of this approach depends on the accuracy of the interpolation scheme.Due to CPU time, we chose a linear interpolation, which, unfortunately, provides numerical dissipations.However, this drawback is overcome using the Fourier transform performed in the second step of the 3DCLOUD algorithm (see Sect. 2.2.2).
As the Fourier transform is easy to implement, this method was chosen to solve the equation ∇ • u = 0.In the Fourier domain, the gradient operator ∇ is equivalent to the multiplication by ik, where i ≡ √ −1 and k is the wave number vector.Thus, the following equation ik.û (k) = 0, where û is the transform of wind velocity u in the Fourier domain, has to be solved.This implies that the Fourier transform of the velocity of a divergent free field is always perpendicular to its wave numbers.Therefore, the quantity 1 k 2 k. û (k) k is removed from û. Keeping the real part of inverse transform of û provides the new wind velocity u with the desired free divergent property.
Lateral periodic conditions and continuity conditions to bottom and top are applied.For wind velocity, free slip boundary conditions are applied at the bottom and top of the domain, which are assumed to be a solid wall (i.e.w = 0).But, as the Fourier transform (which is needed to solve the equation ∇ • u = 0) requires periodic conditions, it provides spurious oscillations during the simulations.In order to limit this effect, extra levels with wind velocity set to zero are added under and above the model domain.
The 3DCLOUD algorithm to simulate 3-D structures of liquid water content (LWC) or IWC is, in summary: 1. Definition of initial meteorological profiles u ini , v ini , θ ini , q v ini from idealized cloud conceptual models or from the user.The vertical pressure profile is generally computed from the hydrostatic law, but can be provided by the user.
2. Initial perturbations u are added to wind velocity u. u is free-divergent and turbulent with a spectral slope of −5/3 (see more explanations in Sect.2.2.2).
5. Computation of ice fall speed v fall (only for cirrus cloud, see Eq. 8).
6. Advection of u, v, w, θ, q v and q c by wind velocity u (see Eq. 10).
7. Modification of θ, q v and q c due to evaporation/condensation processes (see Eq. 7).
8. Modification of the vertical velocity due to buoyancy (see Eq. 9). 9. Modification of q v ini in order to assimilate cloud coverage C (optional, only for stratocumulus and cumulus, see Sect.2.1.2).
10.Return to (3) until maximum iteration number is reached.

Step 2: statistical adjustment
Hereafter, we present the second step of the 3DCLOUD algorithm that is the methodology to adjust, according to user requirements, the mean optical depth τ (or the mean ice water content IWC) and the inhomogeneity parameter of the optical depth ρ τ (or the ice water content ρ IWC ) from the LWC (or from the IWC) simulated at the step 1.The distribution of τ or IWC is assumed to follow a gamma distribution.Indeed, distribution of τ and IWC are usually well represented by a lognormal or gamma distribution (Cahalan et al., 1994;Barker et al., 1996;Carlin et al., 2002;Hogan and Illingworth, 2003;Hogan and Kew, 2005).The scale invariant cloud properties, controlled at each level, are characterised by the spectral exponent β 1-D close to −5/3 (slope β of the one dimension wave number spectrum in log-log axes of the Fourier space).This spectral slope is computed from the outer scale L out (defined by the user) to the smaller scale (voxel horizontal size).

Control of the mean and of the inhomogeneity parameter
The relation between local optical depth τ (x, y, z), liquid water path (LWP) and density of water ρ l in each voxel is given by (Liou, 2002) where x, y, z are the spatial positions inside the simulation domain, LWP (x, y, z) is the local liquid water path and z is the vertical resolution.Local quantity means that the quantity is estimated inside a voxel.The optical depth τ (x, y) for each pixel is the sum of local optical depths along the z axis: The mean optical depth τ is then defined as follows: In the same way, for ice cloud, the mean IWC is obtained with where N * z is the number of layers between cloud top and cloud bottom.
To describe the amplitude of the optical depth for 1-D and 2-D overcast cloud, Szczap et al. (2000) defined the inhomogeneity parameter of optical depth ρ τ .For 3-D broken fields, this parameter is defined according to where σ τ >0 (x, y) and τ >0 (x, y) are the standard deviation and the mean of the strictly positive optical depth.
Due to the flexibility of the mathematical formulation of the gamma distribution and to its ability to mimic the attributes of other positive-value distributions, such as lognormal and exponential distributions, we choose to control τ or IWC and ρ τ or ρ IWC by mapping theoretical gammadistributed properties onto the simulated properties.This mapping technique is analog to the "amplitude adaptation" technique explained in Venema et al. (2006), where amplitudes are adjusted based on their ranking.The gamma distribution is a two-parameter family of continuous probability distribution.It has a shape parameter a and scale parameter b.The equation defining the probability density as a function of a gamma-distributed random variable Y is as follows: where (.) is the gamma function.We develop a simple iterative algorithm, where values of a and b are adjusted until mean and inhomogeneity parameters reach the required user values within few percents.

Control of invariant scale properties by adjustment of spectral exponent in Fourier space
The spectral slope value β 1-D of the horizontally 2-D field is adjusted according to the following methodology.As proposed by Hogan and Kew (2005), we choose to manipulate the 2-D plan of Fourier amplitudes of local optical depth τ 2-D (or IWC 2-D ) with a 2-D Fourier transform performed at each height of cloudy layer.Suppose a 2-D isotropic field g (x, y) characterized by a Gaussian probability density function (PDF) and a 1-D power spectrum E 1 (k) with a spectral slope β 1-D at all scales defined as follows: where k is the wave number in any direction and E 1 is the spectral energy density at k = 1 −1 m .Following Hogan and Kew (2005), for the idealized case where g (x, y) is continuous at small scales and infinite in extent, its 2-D spectral density matrix E 2 k x , k y can be written as where k = k 2 x + k 2 y and κ a constant.In general, a 2-D cloud layer of τ 2-D (or IWC 2-D ) is anisotropic and, in our case, the optical depth (or IWC) is gamma-distributed.Therefore, the 1-D power spectrum E 1 (k) seldom has the required spectral slope β 1-D .In this context, a numerical method has to be developed to perform our objectives.
Let set Y 2-D be the 2-D Fourier transform of Y 2-D , where Y 2-D can be τ 2-D (or IWC 2-D ) at a given cloudy layer.This quantity, estimated with the help of a direct 2-D fast Fourier transform algorithm can be written as follows: where E 2-D = Y 2-D is the magnitude or spectral energy, φ 2-D (k) are the phase angles and k = k 2 x + k 2 y is the absolute wave number.The cloud field domain is defined to measure L x and L y and they have spatial resolutions of x and y.The resulting wave number for E 2-D ranges from −K x to +K x with a resolution of k x = 1 L x , where It is similar for k y direction.
Our objective is to modify spectral energy E 2-D (k) in such a way that the 1-D spectral slope value µ 1-D estimated in one dimension from Y 2-D for k ≥ k out (k out = 1 L out ) satisfies the desired value β 1-D required by the user.Practically, µ 1-D is defined as follows: where β x and β y are the 1-D spectral slope values of Y 2-D estimated in the x and y directions respectively.In order to conserve the spatial repartition of Y 2-D , we keep φ 2-D (k) phase angles unchanged for all values of k.We also 18) and is defined as follows: whereas E * * 2-D (k) is defined as follows: where X is the mean of variable X.If the degree of anisotropy of Y 2-D is small, such as for stratocumulus and cumulus, we use E * 2-D (k) and if not, such as for cirrus clouds, E * * 2-D (k).Nonetheless, the user can also choose one of these methods.
Finally, the new 2-D local optical depth (or the 2-D new ice water content Y new 2-D ) at the given cloud layer is computed by keeping the real part of the inverse 2-D fast Fourier transform of the new quantity: But as a result, the distribution of Y new 2-D is not the same Y 2-D at the given cloud layer, and the equality between the estimated spectral slope µ 1-D of Y new 2-D and the required valuesβ 1-D is not always guaranteed.Therefore, we have to redo an "amplitude adaptation", as explained in Venema et al. (2006), and iterate the process explained in this section by changing the value of β 1-D in Eqs. ( 21) or ( 22), until the estimated value µ 1-D reaches the required value within a few percent.

Implementation
We describe here the part of the 3DCLOUD algorithm that establishes the cloud field mean optical depth τ (IWC), the inhomogeneity parameter ρ τ (or ρ IWC ) and the spectral exponent β 1-D .
For stratocumulus and cumulus clouds, the algorithm is the following.

F. Szczap et al.: A flexible three-dimensional cloud generator (3DCLOUD)
For cirrus clouds, the algorithm is as follows.
1. Transformation of IWC (x, y, z) to IWC (x, y, z) with the algorithm explained in Sect.2.2.1 in order to constrain IWC and ρ IWC values.
2. Application of the algorithm explained in Sect.2.2.2 in order to constrain β 1-D of each cloud layer of IWC (x, y, z).We obtain IWC (x, y, z).
3. Transformation of IWC (x, y, z) to IWC 3-D (x, y, z) with the algorithm explained in Sect.2.2.1 in order to constrain IWC and ρ IWC values.

Differences between 3DCLOUD, IAFFT method and Cloudgen models
Both IAAFT (Venema et al., 2006) and Cloudgen (Hogan and Kew, 2005) models are purely stochastic Fourier based approaches that are able to generate synthetic or surrogate cloud.On the contrary, 3DCLOUD solves, in a first step, basic atmospheric equations, in order to generate an intermediate cloud field.In its second step, as for both IAAFT and Cloudgen models, it uses Fourier tools (manipulation of energy and phase in frequency space) and amplitude adaptation (manipulation of distributions) in order to generate the final cloud field.IAAFT and Cloudgen are designed to simulate stratocumulus/cumulus fields for the first and cirrus fields for the second, when 3DCLOUD is able to simulate stratocumulus, cumulus and cirrus field within the same framework.More specifically, the IAAFT method is designed to generate surrogate clouds having both the amplitude distribution and power of the original cloud (2-D LWC from 1-D LWP measurement, 3-D LWC from 2-D LWC fields or 3-D LWC from 3-D fields generated by LES).It needs LES inputs or measurements.As explained in Venema et al. (2006), stratocumulus often display beautiful cell structures, similar to Bénard convection, and LES clouds show such features.But their 3-D IAFFT surrogates show these much less and do not show fallstreak or a filamentous structure.Due to the specific manipulations of Fourier coefficients presented in the paper, we show that 3DCLOUD is able to simulate the cell structure of stratocumulus (see Figs. 7c and 8), the filamentous structure of cirrus (see Fig. 13) and the cirrus fallstreaks (see Figs. 14 and 15) relatively well.Moreover, the objective of 3DCLOUD is not to provide many surrogate clouds with the same amplitude distribution and power spectrum from an LES original cloud, but to provide 3-D LWC (or optical depth) with the required cloud coverage, the −5/3 spectral slope (often observed in real clouds), the mean value of the gamma distribution of the optical depth and the inhomogeneity parameter, all these parameters being very pertinent for radiative transfer.
Cloudgen is designed to simulate surrogate cirrus with the cirrus specific structural properties: fallstreak geometry and shear-induced mixing.It first generates a 3-D fractal field by performing an inverse 3-D Fourier transform on a matrix of simulated Fourier coefficients with amplitude consistent with observed 1-D spectra.Then random phases are generated for the coefficient allowing multiple cloud realizations with the same statistical properties.Horizontal slices from the domain are manipulated in turn to simulate horizontal displacement and to change the spectra with height.The final field is scaled to produce the observed mean and fractional standard deviation of ice water content.3DCLOUD does not use a 3-D fractal field, but a 3-D IWC field simulated by the simplified atmospheric equation set.Therefore, cloud structures due to wind shear are physically obtained by taking into account the advection (a nonlinear term in momentum equation) rather than by a linear horizontal displacement of phase.Afterwards, in the current version of 3DCLOUD, for each level, 2-D horizontal slices of this 3-D IWC are manipulated in 2-D Fourier domain in such a way that the Fourier coefficient amplitude is consistent with the 2-D spectra of the simulated IWC, with the constraint that the 1-D spectral slope is equal to −5/3 (this value can be change easily in future version of 3DCLOUD).At each level, the 2-D phase for the coefficient is kept unchanged.Finally, the mean value of the 3-D IWC and the inhomogeneity parameter are adjusted.As explained by Hogan and Kew (2005), it is difficult with Cloudgen to generated anisotropic cirrus structure such as roll-like structure near cloud top.3DCLOUD, using physically based equations, allows simulating such kinds of anisotropy, as for example, 3DCLOUD Kelvin-Helmholtz wave breaking (see Fig. 14).

Comparison between 3DCLOUD and large-eddy simulation (LES) outputs
The objective of this section is to check that the basic atmospheric equations used in 3DCLOUD (see Sect. 2.1.1)are solved correctly.We compare 3DCLOUD and LES outputs found in the scientific literature for marine stratocumulus, cumulus and cirrus regimes.Note that assimilation techniques of meteorological profiles and of cloud coverage (see Sect. 2.1.2) are not used here, except in Sect.3.4.
The test cases come from the output LES numerical files provided by the Global Water and Energy Experiment (GEWEX) Cloud System Studies (GCSS) Working Group 1 (WG1) and Working Group 2 (WG2), easily downloadable from the web.They are often used as a benchmark.We choose the DYCOMS2-RF01 case (the first Research Flight of the second Dynamics and Chemistry of Marine Stratocumulus) for the marine stratocumulus regime (Stevens et al., 2005), the BOMEX case (Barbados Oceanographic and Meteorological Experiment) for the shallow cumulus regime (Siebesma et al., 2003), and the ICMCP case (Idealized Cirrus Model Comparison Project) for cirrus regimes (Starr, 2000).

DYCOMS2-RF01 (GCSS-WG1) case
We remember briefly the conditions of simulations and configurations, explained in detail in Stevens et al. (2005).A 4 h simulation on a horizontal grid of 96 by 96 points with 35 m spacing between grid nodes was required.Vertical spacing was required to be 5 m or less.In 3DCLOUD, we thus set N x = N y = 96 and N z = 240, L x = L y = 3.5 km and L z = 1200 m, so that x = y ≈ 36.5 m and z = 5 m.Initial profiles of the liquid water potential temperature θ l and of the total water mixing ratio q t are θ l = 289.0K and q t = 9.0 g kg −1 if z ≤ z i and θ l = 297.5 + (z − z i ) 1/3 K and q t = 1. Figure 2 shows the evolution of the mean cloud top height, the mean cloud base height, the cloud coverage and the liquid water path from the master ensemble and for 3DCLOUD during the 4 h simulations.Even though we can notice slight discrepancies between 3DCLOUD and master ensemble results in the first 2 h ("spin-up" period), 3DCLOUD results are quite consistent with master ensemble results, especially at the end of the simulation.Nevertheless, 3DCLOUD tends to generate a lower cloud height than the mean results with a higher liquid water path.
Figure 3 shows the mean profiles averaged over the fourth hour of the long wave net flux, the liquid water potential temperature, the total water mixing ratio, the liquid water mixing ratio, the horizontal velocity components, and the air density.Even though the 3DCLOUD long wave net flux is smaller compared to master ensemble, again 3DCLOUD results are quite consistent with other results for all the variables.

BOMEX (GCSS-WG1) case
For the BOMEX case (Siebesma et al., 2003), a 6 h simulation on a horizontal grid of 64 by 64 points with 100 m spacing between grid node was required.Vertical spacing was required to be 40 m.In 3DCLOUD, we set thus N x = N y = 64 and N z = 76 with L x = L y = 6.4 km and L z = 2980 m, so x = y = 100 m and z ≈ 39.2 m.Initial profiles of the liquid water potential temperature θ l and the total water mixing ratio q t and the other requirement including geostrophic winds, divergence due to the subsidence, surface sensible heat flux, surface latent heat flux, momentum surface fluxes, moisture large scale horizontal advection term and long wave radiative cooling (radiative effects due to the presence of clouds are neglected) are presented in Appendix B in Siebesma et al., 2003).
Figure 4 shows the evolution of the cloud coverage and the liquid water path from the master ensemble and for 3DCLOUD during the 6 h simulation.We can notice   the small value of the cloud coverage (less than 10 %).3DCLOUD results are quite consistent with the master ensemble results, even if the simulated 3DCLOUD liquid water path (LWP) may be too low at the end of the simulation.
Figure 5 shows mean profiles, averaged over the fifth hour of the cloud coverage, potential temperature, water vapour mixing ratio, liquid water mixing ratio, horizontal velocity components, and air density.The 3DCLOUD results are again quite consistent with the master ensemble results.We note, however, that the 3DCLOUD cloud coverage and liquid water mixing ratio are smaller at all heights.We also see small differences (less than 1 m s −1 ) for the wind velocity below 500 m, and for the potential temperature (less than 1 K) and water vapour mixing ratio (less than 1 g kg −1 ) for altitudes 1800 m.

ICMCP (GCSS-WG2) case
For the cirrus case detailed in Starr et al. (2000), the baseline simulations include night-time "warm" cirrus and "cold" cirrus cases where cloud top initially occurs at −47 • C and −66 • C, respectively.The cloud is generated in an ice super-saturated layer with a geometric thickness around 1 km (120 % in 0.5 km layer) and with a neutral ice pseudoadiabatic thermal stratification.Cloud formation is forced via an imposed diabatic cooling over a 4 h time span followed by a 2 h dissipation stage without cooling.All models simulate radiative transfer, contrary to 3DCLOUD.In 3DCLOUD, we set N x = N y = 60 and N z = 140 with L x = L y = 6.3 km and L z = 14 km, so that x = y = 105 m and z ≈ 100 m.
Figure 6 shows the evolution of the ice water path (IWP) from the master ensemble and for 3DCLOUD during the 6 h simulation.In a general way, most of the tested models and 3DCLOUD have similar behaviour: indeed, the IWP increases during the first 4 h simulation (cirrus formation due to imposed cooling) and decreases after (cirrus dissipation due to non-imposed cooling).The IWP range of the tested models is very large (factor of 10), but we can notice that 3DCLOUD behaviour is closer to bulk microphysics models behaviours, especially for "warm cirrus".
For "cold" cirrus, the 3DCLOUD IWP is smaller than most participating models during all the simulation duration.It is probably because 3DCLOUD does not account for the radiative transfer, as opposed to the participating models.Indeed, neglecting cirrus top cooling due to radiative processes restricts the formation of thin "cold" cirrus.This radiative diabatic effect is probably less important for the "warm" cirrus because the latent heat diabatic effect is larger.

Comparison between 3DCLOUD and BRAMS for the DYCOM2-RF01 case
In order to underscore differences between 3DCLOUD and LES for comparable scenes, we choose again the well documented DYCOMS2-RF01 case.Snapshots can be found, for example, in Stevens et al. (2005) and in Yamaguchi and Feingold (2012).We performed the 4 h simulations of the DYCOMS2-RF01 case with 3DCLOUD and with the Brazilian Regional Atmospheric Modelling System (BRAMS v4) model (Pielke at al., 1992;Cotton et al., 2003).BRAMS simulations were provided by G. Penide (Penide et al., 2010).The BRAMS model is constructed around the full set of nonhydrostatic, compressible equations.The cloud microphysics parameterization is based on a two-moment scheme (Meyers et al., 1997).Subgrid scale fluxes are modelled following Deardroff (1980).The base calculations are performed on a 100 × 100 × 100 point mesh with a step time of 0.3 s.
Figure 7 shows the instantaneous cloud-field snapshots of the pseudo albedo (see definition in Sect.performed in the second step of the 3DCLOUD algorithm.By contrast, the BRAMS optical depth spectral slope is close to −5/3 only in the 2 × 10 −3 : 5 × 10 −3 ≈ 1/ (5 x) m −1 wave number range.Depending on their degree of sophistication, LES do not always guarantee cloud invariant scale properties at the larger wave numbers.Indeed, Bryan et al. (2003) have shown, that for the finite-difference model, the vertical wind velocity spectral slope is steeper than −5/3 for scales shorter than 6 x.Table 1 shows the computation performance of 3DCLOUD and BRAMS.For this specific case, 3DCLOUD simulation is 30 times faster than BRAMS simulation.

Examples of 3DCLOUD possibilities
In this section, we present cloud fields generated by 3DCLOUD with the assimilation of idealized meteorological profiles and fractional cloud coverage defined by the user.We also show the effect of the outer scale L out and the inhomogeneity parameter of optical depth ρ τ on the generated optical depth field.We also give an example of cirrus clouds with fallstreaks.In order to have a spatial representation of the clouds as seen from above, we choose to show the socalled pseudo-albedo α defined as follows: where the asymmetry parameter g is set to 0.86 and τ is the optical depth.

Stratocumulus and cumulus fields with assimilation of meteorological profiles based on DYCOMS2-RF01 and BOMEX cases
We choose to simulate stratocumulus and cumulus LWC in the context of DYCOMS2-RF01 and BOMEX cases.With this aim, we assimilate temperature and humidity initial profiles for stratocumulus and cumulus given by Stevens et al. (2005) and Siebesma et al. (2003), respectively.However, in order to mimic the sensible and latent heat, these profiles have to be slightly modified.At sea surface (z = 0), for stratocumulus (DYCOMS2-RF01 case), the liquid potential temperature and total mixing ratio are set to 290 K and 10 g kg −1 , respectively (instead of 289 K and 9 g kg −1 , respectively).For cumulus (BOMEX case), the liquid potential temperature is set to 299.7 K instead of 298.7 K.In addition, wind profiles assimilated by 3DCLOUD are those computed by the master ensemble at the end of simulation.

Effects of numerical spatial resolution
The effects of the numerical spatial resolution on 3DCLOUD simulations are presented herein.Figure 8 shows pseudoalbedo and cross sections of the vertical velocity and cloud water at the end of the simulation for the stratocumulus case based on the DYCOMS2-RF01 experiment.It also shows the mean profiles of potential temperature, liquid water mixing ratio, horizontal velocity components, and vapour water mixing ratio for different numerical spatial resolutions x = y = 200 m, 100 m, 50 m and 25 m.Horizontal extensions L x = L y are set to 10 km and vertical resolution z to 24 m for all simulations.Figure 9 is the same as Fig. 8, but for the cumulus case with assimilation of meteorological profiles based on the BOMEX case.The vertical resolution z is set to 38.5 m in this last case.
It is obvious that change in the horizontal mesh leads to a more pleasant and detailed flow visualization but there is no significant impact on the mean statistics of the simulated temperature vertical profile, water vapour mixing ratio and wind velocity.The water mixing ratio simulated by 3DCLOUD for the DYCOMS2-RF01 case is very close to the mean profile averaged over the fourth hour and provided by the master ensemble, even if the vertical resolution used in this section is only z = 38.5 m compared to z = 5 m in Sect.3.1.For the BOMEX case, the water mixing ratio simulated by 3DCLOUD changes as a function of the numerical spatial resolution.This behaviour is quite understandable as results drawn on Figs. 8 and 9 are snapshots at the end of the 3DCLOUD simulation and not average results over 1 h as done on Figs. 3 and 5.Moreover, BOMEX meteorological conditions cause time dependent cumulus fields, contrary to DYCOMS2 meteorological conditions that cause more stationary stratocumulus fields.
In addition, it is expected for the BOMEX case, that cloud spacing converges at high spatial resolution.In order to investigate it, we defined an estimator of the cloud spacing called the mean distance D mean .To compute D mean , the 3-D LWC is vertically projected on the 2-D x − y plan in order to obtain the 2-D binary image of the cloud coverage with free cloud areas set to 0 and cloudy areas set to 1. Then we compute the mean distance between the cloud cell for the x and y directions to obtain D mean .Figure 10 shows time series of D mean for different horizontal spatial resolution ( x = y = 192, 100, 50, 25, 12.5 and 8.3 m) with a constant vertical resolution ( z = 38.5 m), for cumulus cloud fields simulated by 3DCLOUD after assimilation of the BOMEX case meteorological profiles.The main difference between these simulations and the BOMEX case simulation is the smaller horizontal extension L x = L y = 5 km instead of 10 km in order to access high numerical spatial resolution x = 8.3 m (N x = N y = 600, N Z = 70).Cumulus clouds appear 10 to 20 min after the beginning of the simulation.After 1 h of simulation, D mean is relatively constant with time, meaning that 3DCLOUD has converged.The mean distance averaged over the last half-hour of the 2 h simulation D mean , is also presented in Fig. 10 as a function of the numerical spatial resolution x.D mean is relatively constant for a spatial resolution x smaller than 20 m, showing that BOMEX cloud spacing converges for spatial resolution close to x = 10 m, a value smaller than x = 25 m used in Fig. 9d.
Table 1.Time step, process time for one time step and process time for 2 h-simulation with 3DCLOUD model, as a function of the numerical resolution.DYCOMS2-RF01 and BOMEX cases are presented.A comparison between 3DCLOUD and BRAMS LES computation time for a specific DYCOMS2-RF01 case is added.3DCLOUD (Matlab code) runs on a personal computer with Intel Xeon E5520 (2.26 GHz) and BRAMS (Fortran code) runs on a PowerEdge R720 with Intel Xeon E5-2670 (2.60 GHz), both of them having a single-processor configuration.

Study case
Point mesh  Table 1 shows the time step, process time for one time step and process time for 2 h-simulation with 3DCLOUD model, as a function of the numerical resolution.DYCOMS2-RF01 and BOMEX cases are presented.The process time for 2 hsimulation is indicated because 3DCLOUD algorithm convergence is achieved after 2 h (or less) of simulation for stratocumulus, cumulus and cirrus regimes (see Fig. 10 for cumulus case).For both cases, the smaller the spatial resolution, the smaller the step time and the larger the process time.A comparison between 3DCLOUD and BRAMS LES computation time for a specific DYCOMS2-RF01 case is added (see Sect. 3.4).For this specific case, 3DCLOUD simulation is 30 times faster than BRAMS simulation.Note that 3DCLOUD (Matlab code) runs on a personal computer with Intel Xeon E5520 (2.26 GHz) and BRAMS (Fortran code) runs on a PowerEdge R720 with Intel Xeon E5-2670 (2.60 GHz), both of them having a single-processor configuration.

Assimilation of the fractional cloud coverage C
Results shown in Fig. 11 are the same as Fig. 8 but with the addition of the cloud coverage assimilation C = 99 %, C = 80 %, C = 50 % and C = 20 %.Horizontal extensions L x and L y are set to 10 km and vertical resolution z is set to 24 m.They show that 3DCLOUD is able to assimilate correctly fractional cloud coverage of stratocumulus for very different values of C, even though the extreme example with C = 20 % is a fair weather cumulus field rather than a stratocumulus.For each value of C assimilated, it is interesting to note that cloud base and cloud top heights are still localised around 600 m and 800 m, respectively.Temperature vertical profiles are almost unchanged.The water mixing ratio vertical profiles decrease with the assimilated C value.

Effect of the outer scale L out and inhomogeneity parameter ρ τ on the optical depth field
We saw that 3DCLOUD can, at the end of step 1, simulate stratocumulus and cumulus fields with enough coherent statistics profiles.However, optical depth (for stratocumulus and cumulus) or IWC (for cirrus) generated during step 1 of 3DCLOUD does not show scale invariant properties observed in real cloud and often characterised by the spectral exponent β 1-D close to −5/3.As described in Sect.2.2, it is the main task of the step 2 of 3DCLOUD, in addition to the adjustment of the mean and standard deviation of optical thickness or IWC.We focus on the DYCOMS2-RF01 case at the spatial resolution x = 50 m (Fig. 8c).The effective radius R eff is set to 10 µm to compute optical depth from liquid water content.The mean optical depth of this initial field is set to 10 and we change the inhomogeneity parameter ρ τ and the outer scale L out to 0. case 1 to 0.7 and 1 km for case 2 and to 0.7 and 10 km for case 3. Figure 12 shows pseudo-albedo, mean power spectra, probability density function of optical depth fields, mean vertical profiles of the horizontally averaged optical depth for the three cases and a volume rendering of optical depth.
First, we notice that the pseudo-albedo of the initial optical depth field (see Fig. 12a) is smoother than the pseudoalbedo of case 1, 2, and 3. Between cases 1 and 2, we clearly see an increase in heterogeneity as case 1 is a quasihomogenous stratocumulus with a small value of ρ τ = 0.2 and case 2 is more inhomogeneous with a larger value of ρ τ = 0.7.Between cases 2 and 3, we can see the effect of the outer scale.In accordance with smooth variations, the spectral slope of the initial optical depth is close to −3 for the 10 −3 : 10 −2 m −1 wave number range (Fig. 12e).Cases 1, 2 and 3 present the proper spectral slope value of −5/3.For cases 1 and 2, this slope is obtained for the 10 −3 : 10 −2 m −1 wave number range, which is coherent with the imposed value of outer scale L out = 1 km.For case 3, L out = 10 km, so the spectral slope should be −5/3 on the 10 −4 : 10 −2 m −1 wave number range.However, we note that this spectral slope value is achieved only for the 5 × 10 −3 : 10 −2 m −1 wave number range, because we keep the phase angles unchanged in the 3DCLOUD algorithm.
In Fig. 12f, we represent the optical depth distributions.The initial optical depth distribution does not follow a common distribution, whereas the optical depth distribution for cases 1 and 2 are log-normal.Indeed, in the 3DCLOUD algorithm, a gamma distribution for the optical depth is imposed.
For case 2 and 3, optical depth distributions are very close, even if the outer scales are different.Thus, changing the L out value does not affect significantly the shape of optical depth distribution.In Fig. 12g, we can see that the horizontal mean optical depth profiles are quasi identical for all cases.
These results show undeniably the flexibility of the 3DCLOUD algorithm.Indeed, in step 2, 3DCLOUD is able, by mapping a theoretical gamma-distributed optical depth onto the optical depth field simulated at step1, to adjust, quasi-independently, the optical depth mean value, the inhomogenenity parameter value of optical depth and the spectral slope value of optical depth for 1 L out : 1 2 x m −1 wave number range.

Cirrus fields with assimilation of idealized meteorological profiles
Including some modifications presented in Sect.2.1.1,3DCLOUD is also able to generate cirrus cloud.We briefly present in this section an example of ice water content (IWC) of cirrus with fallstreaks.For cirrus, we chose to generate IWC field instead of optical depth field as for stratocumulus or cumulus.
Figure 13 shows idealized vertical profiles of potential temperature, relative humidity, and horizontal velocity components assimilated by 3DCLOUD as well as the ice water path (IWP) simulated at step 1 by 3DCLOUD.It also shows   the IWP simulated by 3DCLOUD during step 2, the initial and corrected mean power spectra, the initial and corrected probability density functions and the IWC volume rendering.Horizontal extensions L x = L y and vertical extension L z are set to 10 km and 12.5 km, respectively.Horizontal resolutions x = y and vertical resolution z are set to 24 m and 83.3 m, respectively.IWC is set to the value obtained at the end of step 1 (0.54 mg m −3 ).The inhomogeneity parameter ρ IWC is set to 1.0 and the outer scale L out to 1 km.
Initial meteorological profiles assimilated by 3DCLOUD have been constructed in such a way that thin cirrus is generated between 9.5 km and 10.5 km with fallstreaks.Vertical profiles of potential temperature, and especially their vertical gradients under and above the cirrus are based on those proposed by Liu et al. (2003).In order to generate instabilities due to radiative cooling (not simulated by 3DCLOUD), we imposed a null vertical gradient of the potential temperature near the cirrus top height.We imposed a mean relative humidity with respect to ice (RHI) of 104 % between 9.5 km and 10.5 km.Just above the cloud, RHI is set to 50 % then 20 % near 12 km in altitude.Under the cirrus, RHI decreases with height to 85 % near 8 km.To generate fallstreaks, we imposed larger wind shear inside the cirrus than under the cirrus.
In Fig. 13, we note that IWP obtained after step 1 is smoother than IWP obtained after step 2, and that the initial IWC spectral slope value after step 1 is much smaller (around −5.5) than the corrected IWC spectral slope after step 2 (around −1.6) in the 10 −3 : 2 × 10 −2 m −1 wave number range.For wave number smaller than 1 L out = 10 −3 m −1 , the power spectra are constant.The corrected IWC probability distribution is exponential-like distribution after step 2. This is due to the larger value of ρ IWC = 1.0 used in this example, compared to ρ τ = 0.7 used for stratocumulus in Sect.4.1.3.

Cirrus field and wind shear
We investigate briefly the aspect of cloud organization due to wind shear with 3DCLOUD model and with other stochastic models.We focus on the work of Marsham and Dobbie (2005) and of Hogan and Kew (2005).These two studies are very pertinent together.Indeed, based on RADAR retrievals of IWC from the Chilbolton 94 GHz RADAR on 27 December 1999, which shows a strongly sheared ice cloud (named hereafter RC99 case), Marsham and Dobbie (2005) investigated shear effects by simulating the RC99 case with the UK Met office LES.In contrast, Hogan and Kew (2005)  To configure 3DCLOUD in order to simulate the RC99 case, we assimilate meteorological profiles (potential temperature, horizontal wind velocity) based on those drawn in Fig. 2 in Marsham and Dobbie (2005).We run also the RC99a case with no wind (and therefore no wind shear), and the RC99b case where the potential temperature profile (drawn in Fig. 15  Figure 14 shows a 2-D vertical slice of 3DCLOUD IWC at an angle parallel to the wind of the RC99a case, the RC99 case and the RC99b case.Note that the 3DCLOUD fields are smooth because they are obtained at the first step of the algorithm.These three snapshots are very similar to those presented in Marsham and Dobbie (2005), allowing us to confirm that our basic atmospheric equations are correctly solved.In Fig. 14a, we can see small structures (few km) at above 7 km, due to radiative cooling at the cloud top and latent heat release in the updraughts.Below, we can observe fallstreaks advected (or not if there is no wind) relative to their source in the convective layer by the shear.The shear homogenizes the fallstreaks.Figure 14b clearly shows that 3DCLOUD simulations at the first step of the algorithm homogenize the fallstreaks a lot, certainly too much compared to the RADAR retrievals of IWC from the Chilbolton 94 GHz RADAR on 27 December 1999 (see Fig. 1 in Hogan and Kew, 2005).Figure 14c shows the RC99b case where 3DCLOUD model is able to simulate Kelvin-Helmholtz wave breaking, a dynamic aspect difficult to simulate with purely stochastic models.
Figure 15a  vertical inhomogeneity parameter to match the one estimated crudely from the RADAR retrievals, we obtain Fig. 15b.Details in the layer under 7 km of this snapshot are quite similar to those obtained by Hogan and Kew (2005).
Finally, Fig. 16, which is similar to Fig. 2 in Hogan and Kew (2005), shows vertical profiles of the cloud fraction, the mean of logarithm of IWC for the cloudy voxels and the standard deviation of logarithm of IWC for cloudy voxels computed from two 3DCLOUD fields obtained at the second step of the algorithm for RC99 case.As for Fig. 15, inhomogeneity parameter is ρ IWC = 0.4 for the first cloud field and depends on height and it is derived from RADAR retrievals for the second case.For both cloud fields, the cloud coverage is equal to 1 between 5 km and 10 km.RADAR retrievals of IWC shows that the cloud coverage is equal to 1 only between 5.5 km and 7 km, and decreases to 0 at 10 km.However, as discussed in Sect.2.1.2,the current version of 3DCLOUD is not able to readily simulate fractional cloud coverage in the cirrus regime.For both cloud fields, vertical profiles of the mean IWC are quite similar and consistent with those retrieved from RADAR.The inhomogeneity parameter vertical profile simulated by the current version of 3DCLOUD is too small, leading to a smoothing of the layer under 7 km that can be improved if the vertical profile of the inhomogeneity parameter is known.

Conclusions
3DCLOUD is a flexible three-dimensional cloud generator developed to simulate with a personal computer and under Matlab environment, synthetic but realistic stratocumulus, cumulus and cirrus cloud fields.Simplified dynamic and thermodynamic laws allow the generation of realistic liquid or ice water content from meteorological profiles.The stochastic process with the Fourier framework allow us to provide ice water content or optical depth sharing similar statistical properties to those observed in real clouds such as the inhomogeneity parameter (set by the user) and the invariant scale properties characterised by a spectral slope close to −5/3 from the smaller scale (set by spatial resolution of grid computation) to the outer scale (set by the user).In order to simulate cloud structures, 3DCLOUD solves simplified basic atmospheric equations and assimilates the cloud coverage set by the user (only for the stratocumulus and cumulus regimes) and meteorological profiles (pressure, humidity, wind velocity) defined by the user.
The 3DCLOUD outputs were compared to LES ones for three classical test cases.We chose the case of DYCOMS2-RF01, the case of BOMEX, and the case of ICMCP.For these cases, results show that 3DCLOUD outputs are relatively consistent with LES outputs, and confirm that the chosen basic atmospheric equations of 3DCLOUD are solved correctly.We also show that, under the condition that the user provides coherent meteorological profiles, 3DCLOUD algorithm is able to assimilate them and generate realistic cloud structures.
3DCLOUD is a very interesting research tool to better understand 3-D interactions between cloudy atmosphere and atmospheric radiation, which is of primary importance in order to make progress in the direct radiative problem (global climate models context) and in the inverse radiative problem (remote sensing context, development of the next generation of atmospheric sensors).For example, 3DCLOUD was used to quantify the impact of stratocumulus heterogeneities on polarized radiation measurements performed by POLDER/PARASOL (Cornet et al., 2013) as well as the influence of cirrus heterogeneities on brightness temperature measured by IIR/CALIPSO (Fauchez et al., 2013(Fauchez et al., , 2014)).
We still have to develop a stochastic process to generate 3-D field of cloud effective radius.In a longer term, investigations will focus on the generation of 3-D mixed phase cloud and eventually on the simulation of 3-D rain rate.Another task will be to provide a FORTRAN code of 3DCLOUD, assumed to be faster than the current Matlab code.

Figure 1 .
Figure 1.General flow chart of the stratocumulus, cumulus and cirrus generator 3DCloud.Note that 3DCloud algorithm is divided in two distinct steps.

Figure 2 .Figure 2 .
Figure 2. Time series of (a) the mean cloud top height, (b) the mean cloud base height, (c) the 2 cloud coverage, and (d) the liquid water path.The DYCOM2-RF01 case is displayed.The 3 solid lines indicate 3DCLOUD results.The dotted lines indicate a mean over all LES results.4 The light shading around this mean delimits the maximum and minimum values within the 5 master ensemble at any given time.6 7 8 9 10 11 12 13 5 g kg −1 if z > z i .Other required forcings include geostrophic winds (U g = 7 m s −1 and V g = −5.5 m s −1 ), divergence of the large-scale winds (D = 3.75 × 10 −6 s −1 ), surface sensible heat flux (15 W m −2 ) and surface latent heat flux (115 W m −2 ).The momentum surface fluxes where the total momentum is specified by setting u * = 0.25 m s −1 and the radiation schemes is based on a simple model of the net long wave radiative flux (see Eqs.3 and 4 in Stevens et al., 2005).

1 Figure 3 .Figure 3 .
Figure 3. Mean profiles averaged over the fourth hour of (a) the longwave net flux, (b) the 2 liquid water potential temperature, (c) the total water mixing ratio, (d) the liquid water mixing 3 ratio, (e) the horizontal velocity components, and (f) the air density.The DYCOM2-RF01 4 case is displayed.The solid lines indicate 3DCLOUD results.The dotted lines indicate a mean 5 over all LES results.The light shading around this mean delimits by the maximum and 6 minimum values within the master ensemble at any given height.7 8 9 10 11 12 13 14 15

Figure 4 .
Figure 4. Time series of (a) the cloud coverage and (b) the liquid water path.The BOMEX case is displayed.The solid lines indicate 3DCLOUD results.The light shading delimits the maximum and minimum values within the master ensemble at any given time.

Figure 4 .
Figure 4. Time series of (a) the cloud coverage and (b) the liquid water path.The BOMEX case is displayed.The solid lines indicate 3DCLOUD results.The light shading delimits the maximum and minimum values within the master ensemble at any given time.

Figure 5 .
Figure 5. Mean profiles averaged over the fifth hour of (a) cloud coverage, (b) the potential temperature, (c) the water vapour mixing ratio, (d) the liquid water mixing ratio, (e) the horizontal velocity components, and (f) the air density.The BOMEX case is displayed.The solid lines indicate 3DCLOUD results.The light shading delimits the maximum and minimum values within the master ensemble at any given height.
4) at 4 hours simulated by (a) the UCLA-0 model (picture taken from Stevens et al., 2005), (b) the BRAMS model, both configured following the DYCOMS2-RF01 case (Stevens et al., 2005) and (c) from 3DCLOUD with assimilation of meteorological profiles based on the DYCOMS-RF01 case.Both BRAMS and 3DCLOUD cases are drawn from simulations where x = y = 40 m and z = 12 m.These three snapshots of cloud fields are characterized by closed cellular convection with large cloud cover, as argued in Yamaguchi and Feingold (2012), who did simulation of DYCOMS-RF01 case with the LES mode of the Advanced Research weather research and forcasting (WRF) model.Figure 7 also shows the power spectra computed following the x and y directions and then averaged, for BRAMS and 3DCLOUD optical depth fields.The 3DCLOUD optical depth spectral slope is close to −5/3 in the L out : 1/ (2 x) m −1 wave number range, as expected, because of the statistical adjustment

Figure 6 .
Figure 6.Time series of vertically-integrated ice water path (IWP) from different cirrus models, which participated in the Idealized Cirrus Model Comparison Project and from the 3DCLOUD model (thick blue lines).The upper panel is for the cold cirrus case and the bottom one is for the warm cirrus case.Cyan line represents models with bin microphysics, red line models with bulk microphysics, green line single column models, and thin black lines models with heritage in the study of deep convection or boundary layer clouds.This figure is made from the one taken from Starr et al. (2000) and Yang et al. (2012).

Figure 6 .
Figure 6.Time series of vertically integrated ice water path (IWP) from different cirrus models, which participated in the Idealized Cirrus Model Comparison Project and from the 3DCLOUD model (thick blue lines).The upper panel is for the cold cirrus case and the bottom one is for the warm cirrus case.Cyan line represents models with bin microphysics, red line models with bulk microphysics, green line single column models, and thin black lines models with heritage in the study of deep convection or boundary layer clouds.This figure is made from the one taken from Starr et al. (2000) and Yang et al. (2012).

Figure 7 .
Figure 7.The instantaneous cloud-field snapshots of the pseudo albedo at 4 hours simulated by (a) the UCLA-0 model (picture taken from Stevens et al., 2005), (b) the BRAMS model, both configured following the DYCOMS2-RF01 case (Stevens et al., 2005) and (c) from 3DCLOUD with assimilation of meteorological profiles based on the DYCOMS-RF01 case.The UCLA-0 field is drawn from simulation where N x = N y = 192 and x = y = 20 m.Both BRAMS and 3DCLOUD are drawn from simulations where N x = N y = N z = 100, x = y = 40 m and z = 12 m.Note that the 3DCLOUD field is obtained at the second step of the algorithm, with the inhomogeneity parameter ρ τ = 0.3, mean optical depth τ = 10 and L out = 2 km.(d) is the optical depth power spectra computed following the x and the y directions and then averaged, for BRAMS (points) and 3DCLOUD (circles).A theoretical power spectrum with spectral slope β = −5/3 is added (black line).
Figure 8. (a), (b), (c) and (d) pseudo albedo and (e), (f), (g), and (h) cross sections of the 2 vertical velocity (shaded) and the cloud water (contoured), at the end of simulation, for the 3 stratocumulus simulated by 3DCLOUD with assimilation of meteorological profiles based on 4 the DYCOMS2-RF01 case.Different numerical spatial resolutions are presented with 5 x y ∆ = ∆ : (a) and (e) 200 x ∆ = m, (b) and (f) 100 x ∆ = m, (c) and (g) 50 x ∆ = m and (d) and 6 (h) 25 x ∆ = m.(i), (j), (k) and (l) mean profiles of the potential temperature, the liquid water 7 mixing ratio, the horizontal velocity components, and the vapour water mixing ratio.The 8 solid lines indicate meteorological profiles based on DYCOMS2-RF01 case and assimilated 9

Figure 8 . 1 Figure 9 .Figure 9 .
Figure 8. (a), (b), (c) and (d) pseudo albedo and (e), (f), (g), and (h) cross sections of the vertical velocity (shaded) and the cloud water (contoured), at the end of simulation, for the stratocumulus simulated by 3DCLOUD with assimilation of meteorological profiles based on the DYCOMS2-RF01 case.Different numerical spatial resolutions are presented with x = y: (a) and (e) x = 200 m, (b) and (f) x = 100 m, (c) and (g) x = 50 m and (d) and (h) x = 25 m.(i), (j), (k) and (l) mean profiles of the potential temperature, the liquid water mixing ratio, the horizontal velocity components, and the vapour water mixing ratio.The solid lines indicate meteorological profiles based on DYCOMS2-RF01 case and assimilated by 3DCLOUD.Points, dotted lines, dashed lines and dash-dot lines indicate 3DCLOUD results at the end of simulation for the different numerical spatial resolution x = 200 m, x = 100 m, x = 50 m and x = 25 m, respectively.Number of iterations is 700.

1 Figure 10 .BOMEXFigure 10 .
Figure 10.Time series of the mean distance between cloud areas for different horizontal 2 numerical spatial resolutions (colored lines) with a constant vertical numerical spatial 3 resolution (∆ = 38.5 m) and mean distance averaged over the last-half hour of a 2h 4 simulation as a function of numerical spatial resolution (black line with circles).The cumulus 5 cloud is simulated by 3DCLOUD with assimilation of meteorological profiles based on the 6 BOMEX case.Horizontal extensions are = = 5 km and vertical extension is = 7 2700 m. 8 9 10 11 12 13 14 15 16

Figure 12 .
Figure 12.(a) Pseudo albedo estimated from optical depth (initial field) simulated in the step 1 by 3DCLOUD for the DYCOMS2-RF01 case (see Fig. 8c), (b), (c) and (d) pseudo albedo adjusted in the step 2 of 3DCLOUD for different values of the inhomogeneity parameter ρ τ and of the outer scale L out .(e) mean power spectra of optical depth along x and y directions.The power spectra are scaled for better visualization.(f)probability density function of optical depth, (g) mean vertical profiles of horizontally averaged optical depth and (h) volume rendering of optical depth for the case 3. ρ τ and L out are 0.2 and 1 km for case 1, 0.7 and 1 km for case 2 and 0.7 and 10 km for case 3, respectively.Solid lines, dotted lines, dashed lines and dash-dot lines indicate initial field, case 1, case 2 and case 3 fields, respectively.

Figure 13 .
Figure 13.Idealized vertical profiles assimilated (dashed lines) and simulated (solid lines) by 3DCLOUD during step 1 of (a) the potential temperature and relative humidity and of (b) the horizontal velocity components and the ice water content (IWC), (c) ice water path (IWP) simulated by 3DCloud in step 1, (d) IWP simulated by 3DCloud in step 2, (e) mean power spectra of IWC along x and y directions after the step 1 and the step 2. (f) IWC probability density functions after step 1 and step 2. (g) IWC volume rendering after step 2. ρ IWC is set to 1 and L out is set to 1 km.Number of iterations is 1000.
Figure 14.2-D vertical slice of 3DCLOUD ice water content (IWC gm 3 ) through a 3-D simulation at an angle parallel to the wind (a) RC99a case, (b) RC99 case and (c) RC99b case.Fields are obtained from simulations where N x = N y = 200 and N z = 66 and x = y = 250 m and z = 120 m.Horizontal extensions are L x = L y = 50 km and vertical extension is L z = 8 km between 4 km and 12 km.Note that the 3DCLOUD fields are smooth because obtained at the first step of the algorithm.
in Marsham and Dobbie, 2005) reduces atmospheric stability in order to give more extensive Kelvin-Helmholtz wave braking.All our simulations are done with N x = N y = 200 and N z = 66 and x = y = 250 m and z = 120 m.Horizontal extensions are L x = L y = 50 km and vertical extension is L z = 8 km between 4 km and 12 km.Note that 3DCLOUD, Marsham and Dobbie (2005) and Hogan and Kew (2005) numerical resolution are x = y = 250 m, x = 100 m and x ≈ 780 m, respectively.Note also that 3DCLOUD, Marsham and Dobbie (2005) and Hogan and Kew (2005) horizontal extensions are L x = L y = 50 km, L x = 50 km and L x = L y = 200 km.