Cloud patterns and mixing properties in shallow moist Rayleigh–Bénard convection

Three-dimensional direct numerical simulations of idealized moist turbulent Rayleigh–Bénard convection are presented. The thermodynamics of moist air is linearized close to the phase boundary between water vapor and liquid water. This formulation allows for a simplified saturation condition for the cloud formation, but omits supersaturation and rain. The sensitivity of this problem to changes of the Rayleigh number, the aspect ratio of the convection layer and the water vapor concentration is studied. The Rayleigh number is found to impact the behavior of the system in multiple ways. First, the relaxation time toward a well-mixed turbulent state increases with the Rayleigh number. Similarly, the flow exhibits a higher spatial and temporal intermittency at higher Rayleigh number. This is in line with an enhanced intermittency of the upward buoyancy flux, which we quantify by a multifractal analysis. In addition, phase transition introduces an asymmetry in the distribution of the thermodynamic properties of the well-mixed state. This asymmetry is most pronounced in layers where clouds are partially present. Furthermore, the geometrical properties of the cloud formations averaged with respect to the height of the layer are studied. Similar to isocontours in scalar mixing, the boundaries of isolated clouds show no strict (mono-)fractal behavior. The results of the perimeter–area analysis of the largest isolated clouds agree well with those of large eddy simulations of cumulus convection. This perimeter–area scaling is also similar to that of percolation processes in a plane.

The correct version of figure 3 is shown here. The caption of the figure and the rest of the discussions and results in the article remain unaffected.

Introduction
The bottom of the Earth's atmosphere is heated by sensible and latent heat flux from the ground surface, while its interior cools as it emits more infrared radiation than it absorbs. This differential heating continuously destabilizes the atmospheric column and, as a result, convective motions develop. Convection is omnipresent in the Earth's atmosphere, but its scale can vary widely, from a few hundred meters for shallow cumulus clouds over the subtropical ocean to 16 km for deep cumulonimbus clouds in the equatorial regions [14]. Most of these convective processes, although not all, involve clouds. Indeed, as moist air rises, it expands and cools adiabatically. An air parcel can thus become saturated after a sufficient ascent and a cloud is formed. The term 'cumulus'-from the Latin word for 'heap'-refers in meteorology to the type of clouds most directly associated with convection.
Clouds are not solely a passive indicator of convective processes, but also have a direct impact on the dynamics of atmospheric convection [7,41]. Indeed, the condensation of water vapor is associated with a release of latent heat, which compensates in part for the adiabatic cooling. The temperature of an unsaturated air parcel rising adiabatically drops by approximately 10 K km −1 (which is known as the dry adiabatic lapse rate d = g/c p with g being the gravity acceleration and c p the specific heat at constant pressure). In contrast, the temperature drop in a saturated air parcel can be as low as 4 K km −1 for very moist air in the tropics (which is known as the moist adiabatic lapse rate s ) [32]. The difference in thermodynamic behavior between saturated and unsaturated air parcels has fundamental implications not solely for convection, but also for many other aspects of the atmospheric circulation such as hurricanes, midlatitude weather systems or the planetary scale circulation.
In this paper, we continue the discussion of a simple mathematical formulation to account for the effect of phase transitions on the dynamical behavior of turbulent convection. This idea was originally proposed by Kuo [16] and then significantly extended by Bretherton [3,4]. Bretherton's model has recently been revisited and further extended in [28,36]. The main simplification arises from replacing the equation of state for moist air by a piecewise linear approximation about the phase boundary. On the one hand, the approximation on which this 3 model is based limits its accuracy to the study of relatively shallow convection, with a depth of at most a few kilometers. On the other hand, this approximation greatly simplifies the mathematical expression for the equation of state while still accounting for the thermodynamic effects of the latent heat release.
This simple equation of state, which is reviewed in section 2, can be directly implemented in a Boussinesq model to simulate convection in a shallow atmospheric slab. In this idealized moist Rayleigh-Bénard convection problem, the bottom of the atmosphere is maintained at a higher temperature and water content, while the top of the atmosphere has smaller fixed values of temperature and water content. Once this is done, the problem can be characterized by five non-dimensional parameters. Thus, while the moist Rayleigh-Bénard problem omits a range of physical processes, including radiative transfer and precipitation, its small number of free parameters should make it more amenable to a systematic exploration of the parameter space.
The direct numerical simulations (DNS) of moist Rayleigh-Bénard convection for Rayleigh number up to 10 8 are analyzed in section 3. Here, we focus primarily on the sensitivity of the simulated cloud field impact to the Rayleigh number and to the water vapor supply. Among other benefits, these DNS provide a set of benchmark simulations to assess whether other models relying on turbulent closure, such as large eddy simulations (LES) [40] or the single column model, can accurately reproduce the cloud field. We summarize our findings and give a brief outlook on future research directions in section 4.

Towards the shallow moist convection model
Modeling atmospheric convection requires a careful treatment of phase transition of water vapor and of the behavior of the large number of water droplets and ice crystals. While the governing equation can be derived from fundamental physical principles (see [24] for a recent discussion), our approach here is to obtain a simplified model that still captures the impact of phase transition on the dynamics of convection. In the following, we review in more detail the assumptions that lead to the moist Rayleigh-Bénard model. As we will see in section 2.2, fluid motion is caused by a buoyancy term that is added to the momentum balance equation in the case of thermal convection. The buoyancy B in atmospheric convection is given by the following relation (see e.g. the textbook by Emanuel [7]): with g being the gravity acceleration,ρ a mean density, p the pressure, S the entropy and q v and q l the mixing ratios of water vapor and liquid water. We exclude the existence of ice in our considerations and limit the discussion to the so-called warm clouds. The sequence of simplifications that result in the present model for shallow moist Rayleigh-Bénard convection reads as follows: • Local thermodynamic equilibrium. It is assumed that the air parcels are in a local thermodynamic equilibrium, which means for example that precipitation or super-cooled water are absent. As soon as the saturation threshold is exceeded, liquid water will be formed. It is well known that supersaturation in the atmosphere does not exceed a few per cent due to the presence of aerosols (∼10 6 m −3 air). Heterogeneous nucleation is thus the 4 dominant process that leads to the formation of cloud water droplets [32,37]. The two mixing ratios are combined with the total water mixing ratio, a new variable of state that remains constant for adiabatic non-precipitating processes: This reduces the number of variables to three and thus B(S, q T , p).
• Boussinesq approximation. We consider small pressure variations about a mean hydrostatic profile. This is known as the Boussinesq approximation that leaves us with a buoyancy B(S, q T , z) and an incompressible fluid (for more details see also [27]).
• Linearization at the phase boundary. The functional dependence of the buoyancy field on the two remaining variables of state and height is highly nonlinear and contains the full thermodynamics of phase changes (see e.g. [2,7]). The next step is to approximate B(S, q T , z) as a piecewise linear function of S and q T around the phase boundary between unsaturated (u) gaseous and saturated (s) liquid phases. The step preserves an important physical property, the discontinuity of partial derivatives of B(S, q T , z) at the phase boundary, which are given by This means that the constants B S,u = B S,s and B q T ,u = B q T,s . These last two steps will limit our considerations to the shallow convection case where the height variations remain moderate. In the atmosphere, this means that heights up to 2-3 km can be considered.
• Dry and moist buoyancy fields. Since B is a linear function of the two variables of state S and q T , we can introduce two new prognostic buoyancy fields, a dry buoyancy field D and a moist buoyancy field M 4 as new linear combinations of S and q T : An implicit relation F(M, D, z) = 0 determines the phase boundary and thus whether an air parcel is saturated or not. In the case of the linearized thermodynamics, an air parcel is thus saturated (and forms in the present model immediately liquid water) when The unknown function f (z) can be shown to be f (z) = −N 2 s z [28], where N s is the Brunt-Vaisala frequency, which is defined as Here, d and s are the dry and moist adiabatic lapse rates and T ref is a reference temperature, e.g. the temperature at the bottom plane. Consequently, the buoyancy field B(M, D, z) is given by

5
The final saturation condition (8) is an explicit and linear relation that can be carried out straightforwardly in the DNS. It is a macroscopic condition that does not resolve cloud microphysics, such as local fluctuations of the supersaturation and the activation of cloud condensation nuclei [37]. On the basis of (8), we can define the liquid water content field as q l = M − D + N 2 s z. Clouds are consequently defined by and the cloud boundary is given by q l = 0.

Moist Boussinesq equations
The dry and moist buoyancy fields can be decomposed into Here Note that the first term on the right-hand side is horizontally uniform. This implies that it can be balanced by a horizontally uniform pressure field. We can thus remove the mean contribution from the buoyancy field without any loss of generality. The moist Boussinesq equations together with the decompositions (10) and (11) and the saturation condition (8) are given by ∂u ∂t Here u is the velocity field, p is the kinematic pressure field, ν is the kinematic viscosity and κ is the diffusivity of the buoyancy fields.
Under most circumstances, the amount of water vapor in the atmosphere decreases with height. This implies that the moist Rayleigh number should be larger than the dry Rayleigh 6 number, Ra M Ra D . 5 In addition to the three parameters explicitly present in the nondimensional versions of equations (13)- (16), two more parameters are hidden implicitly within the definition (12) of the buoyancy B that is given in dimensionless form by These are the so-called surface saturation deficit SSD and the condensation in saturated ascent CSA. The parameters are defined as These two new non-dimensional parameters quantify the amount of saturation that can take place within the domain. The parameter SSD determines the degree of subsaturation of the air parcels at the bottom z = 0. For SSD > 0 the air is unsaturated, and D 0 − M 0 is proportional to the 'water deficit', i.e. the amount of water vapor that must be added for the parcel to become saturated at z = 0. For SSD = 0, moist air is in thermodynamic equilibrium with liquid water at the lower boundary. The second parameter CSA is proportional to the amount of cloud water formed in a saturated parcel rising from the bottom to the top. An alternative choice for the non-dimensional parameters can be obtained based on the degree of saturation of the air parcels at z = 0 and H . From (8), it follows that M 0 − D 0 and M H − D H + N 2 s z determine the amount of cloud water at z = 0 and H , respectively. We can thus define relative amounts of cloud water at the bottom and top by A positive value is proportional to the amount of cloud water present at either the lower or upper boundary. A negative value for either non-dimensional parameters corresponds to an unsaturated parcel and is proportional to the amount of water vapor that must be added to reach saturation. The amount of water here is normalized by N 2 s H , i.e. by the amount of condensation that results from lifting a saturated parcel from the lower boundary to the upper boundary. While the two parameters CSA and SSD arise directly from the non-dimensionalization procedure, we will use here the alternative choice CW 0 and CW H , which can be more easily interpreted in terms of the thermodynamic properties of the lower and the upper boundary. The full moist Rayleigh-Bénard problem is uniquely defined by the five non-dimensional parameters Ra D , Ra M , Pr , CW 0 and CW H . The present study is restricted to the Prandtl number of air at common temperatures, Pr = 0.7, and to a lower boundary in thermodynamic equilibrium with liquid water, CW 0 = SSD = 0.

Numerical scheme and initial configurations
The equations of motion are solved by a pseudospectral scheme with volumetric fast Fourier transformations and 2/3 de-aliasing in a Cartesian slab with side lengths H × H × H . Here is the aspect ratio of the slab and L = L x = L y = H . In the lateral directions x and y, we apply periodic boundary conditions. In the vertical z-direction, we apply free-slip boundary conditions at planes z = 0 and z = H , Here, D and M are the fluctuating fields arising from decompositions of both total buoyancy fields (10) and (11), respectively. These are the same boundary conditions that have been used in the moist convection studies by Kuo [16] and Bretherton [3,4]. The free-slip boundary conditions can be considered as a rough approximation of the ocean surface at the bottom and a temperature inversion at the top. A further very important and interesting case is constant flux boundary conditions, which will not be discussed here. Time stepping is done by a second-order Runge-Kutta scheme. Table 1 summarizes the grid resolutions and dimensionless parameter sets that are taken in the DNS. The spectral resolution does not go below k max η K = 2.45 for all DNS, where k max = 2 √ 2π N x /(3L) is the maximum resolved wavenumber and η K = ν 3/4 / 1/4 the Kolmogorov length. The quantity is the statistical mean of the energy dissipation rate that is given by (x, t) = (ν/2)(∂ i u j + ∂ j u i ) 2 with i, j = 1, 2, 3. We apply a uniform viscosity throughout the simulations, even for the largest aspect ratios.The effects of different viscosities in the lateral and vertical directions for convection have been discussed in [31]. Technically, we use B in the momentum equation (13) instead of B since the mean contribution isB(z), which can be added to the kinematic pressure, i.e. −∂ z p + B = −∂ zp + B . In our simulations, we distinguish two classes of initial equilibrium configurations. These are fully saturated slabs for positive values of CW H and fully unsaturated slabs for negative ones (see figure 1). This means that we start with a layer either completely filled with clouds or free of clouds.

Evolution into the statistically stationary turbulent state
The numerical simulations always started with a static equilibrium solution. As sketched in figure 1, we applied two different classes of initial equilibria. The infinitesimal perturbation is here always implemented by a few three-dimensional small-wavenumber modes. They trigger the primary instability and cause the growth of a regular two-dimensional roll pattern. The pattern is linearly unstable for our Ra M , Ra D Ra c (where Ra c is the critical Rayleigh number for the onset of convection) and crosses over into a fully three-dimensional flow state, via a zigzag-like instability. For the present Rayleigh numbers, the system relaxes eventually into a fully turbulent state after about 100 free-fall time units T f (see figure 2). In the case of dry convection, this sequence of bifurcations to a turbulent state has been studied in much detail by Hartlep et al [13]. Other perturbations, e.g. by white noise, are also possible and lead to the same relaxation times of ∼100T f . This relaxation to a turbulent and well-mixed state can also be investigated in the plane that is spanned by the joint distribution of the two buoyancy fields, D and M. In [3,28], it is demonstrated that the two fields can be combined to a new conserved scalar field which is freely decaying. Over long time scales, the joint distribution for D and M converges toward the diagonal mixing line.
The collapse of the distribution on the mixing line is not instantaneous. This is relevant for atmospheric problems such as the evolution of convection through the diurnal cycle in which the boundary conditions evolve over time. The mixing time scales can be derived from the numerical simulation. For this purpose, we determine the joint probability density function (JPDF) of M and D and follow their evolution in time. In figure 3, this relaxation is shown for run 7. The evolution for run 10 is similar. In order to monitor this relaxation over a longer time period, we took special initial conditions, a fully relaxed turbulent state with SSD > 0 and different values of Ra D and Ra M . To be consistent with the rest of the parametric study, the values of Ra D and Ra M in these initial configurations differed by a factor of ten as in the subsequent relaxation runs. Such an initial condition is a 'badly mixed state' as can be seen in figure 3. The support of the JPDF takes a large fraction of the plane initially (figure 3(a)); it steadily shrinks to an eventually diagonal support. This relaxation can be studied more quantitatively as follows. The mean values of both buoyancy fields are given by The covariance matrixĈ(t) of the JPDF is a 2 × 2 matrix with the following time-dependent elements: where C M D (t) = C D M (t). The evolution of the first eigenvalue of the singular value decomposition ofĈ(t), denoted as σ 1 (t), is shown in figure 4(a). This value that measures the characteristic length of the support around the mixing line remains nearly unchanged after an initial variation. The second eigenvalue, σ 2 (t), characterizes the width of the JPDF perpendicular to the diagonal mixing line. The time evolution of this singular value is shown in figure 4(b) for two Rayleigh numbers. The larger Ra M obeys a slower relaxation into a wellmixed state. We determined the characteristic decay time scale, which is defined via exp(−t/τ ), from the initial slope of the curve. This decay time scale grows from 1.96T f at Ra M = 1.9 × 10 7 to 4.35T f at Ra M = 1.9 × 10 8 .

Statistics of the mixture fraction χ
The description of turbulence in a two-component fluid with a single quantity, the mixture fraction χ, is a concept that has been developed in non-premixed combustion [29]. It turned out that entrainment processes in the cloud top can be described within this formalism as well [22,23,38]. Here, we can define a mixture fraction χ based on the buoyancy fields in the present model. It is given by Note that the standard definition of the mixture fraction χ is based on the ratio between clear and cloudy air [22,23,38]. After a sufficiently long time for the joint distribution to have collapsed on the mixing line, the same field χ would arise when we use the moist buoyancy field M(x, t) instead. In figure 5, we compare the probability density functions (PDF) of χ for different values of CW H at a given value of Ra D in panel (a) and as a function of Ra M at a fixed CW H in panel (b). For Rayleigh-Bénard convection, the first moment of the mixing line must be 1/2 due to the up-down symmetry of the problem. However, this symmetry is broken in the moist Rayleigh-Bénard problem due to the formulation of the buoyancy. As a result, the first moment χ of the mixing fraction can be different from 1/2, as shown in figure 5. The skewness and flatness for the mixture fraction χ are defined as follows: where χ is the first moment of the PDF. The left panel indicates a small dependence on CW H , which manifests in an asymmetry of the PDF. For the run at CW H = −0.43, one of the two largest skewness values of the PDF is found as displayed in the inset of figure 5(a). This observation is in line with the recent findings in [36], where the asymmetry of the vertical profile of the root mean square of the vertical velocity component, u 2 z (z) A,t , is strongest at exactly the same parameters. The inset shows that S(χ) decreases for smaller and larger CW H parameters. Figure 5(b) indicates that the distribution gets narrower with increasing Rayleigh number. We compared therefore runs 3, 7 and 10 at CW H = −1.43. The flatness F(χ) increases from 3.7 at Ra M = 1.9 × 10 6 , 4.4 at Ra M = 1.9 × 10 7 to 8.5 at Ra M = 1.9 × 10 8 , which deviates already significantly from the Gaussian value of 3. The skewness S(χ) varies from 0.04 at Ra M = 1.9 × 10 6 , 0.02 at Ra M = 1.9 × 10 7 to −0.13 at Ra M = 1.9 × 10 8 . The values of F indicate that the statistics of the buoyancy fields (and therefore of the mixture fraction) becomes increasingly intermittent as is known from dry convection [8,19].

Connection between the cloud base and the mixture fraction
In meteorology, the lifting condensation level (LCL) is defined as the first level at which the parcel becomes saturated. In the case of the piecewise linear equation of state, it is given by and a parcel is saturated if z z LCL (M, D). Note that if z LCL is negative, the parcel will always be saturated. In addition, the two non-dimensional cloud water amount parameters CW 0 and CW H defined in (20) are tied to the LCL for air parcels originating from the lower and upper boundaries, as we have In this case, the LCL is a function of the mixture fraction This implies that the LCL of any parcel will always be between the LCL of parcels originating from the upper and lower boundaries. In particular, in our experiments, the lower boundary is saturated, corresponding to CW 0 = 0, so that the LCL at the lower boundary is zero. If the upper boundary is also saturated, i.e. CW H > 0, then all parcels will be saturated for z H (1 − CW H ). This situation guarantees the presence of a stratocumulus layer in the upper part of the domain. Hence, a pure cumulus regime, defined by the absence of a fully saturated layer, can only occur when the upper boundary is unsaturated, i.e. for CW H 0. Equation (32) also indicates that the PDF for the LCL is tied to the PDF of the mixture fraction. Figure 5(b) shows that the distribution for the mixture fraction χ becomes narrower at higher Rayleigh number, and so would the distribution for the LCL. As the distribution of χ becomes narrower, the probability of finding an unsaturated parcel above the mean value of the LCL, or a saturated parcel below the same level, diminishes gradually. One would thus expect a much sharper transition between the saturated and unsaturated regions at high Rayleigh number, a point that is discussed further in section 3.4.4. It remains to be seen whether the sharpening of the distribution of χ shown in figure 5(b) continues for Rayleigh numbers larger than 1.9 × 10 8 , but if this were the case, this would preclude the presence of a cumulus regime in the current configuration of the moist Rayleigh-Bénard convection at very high Rayleigh number.

Cloud pattern analysis
We turn now to the investigation of the cloud patterns and recall first that clouds are defined in our moist convection by inequality (9). The analysis has been performed in parts for runs at different aspect ratios but the same set of Ra D , Ra M , CW 0 and CW H (see table 1). This is to confirm that no artefacts are introduced by the periodic boundary conditions in x and y. Note that the resolution remains unchanged, which causes significant numerical efforts for the largest . Furthermore, the larger the aspect ratio, the longer the relaxation period into the fully developed turbulent state. This is illustrated in figure 6(a) where we plot the turbulent kinetic energy versus time for aspect ratios = 4, 8, 16 and 32. It can be seen in figure 6(b) that the cloud cover-which is defined as the probability to find clouds at a height z (see [36])-is not affected significantly by .
Our geometric analysis will be based on vertically averaged quantities. The motivation for this kind of analysis comes from LES of mesoscale convection and satellite image analysis of cloud patterns in atmospheric mesoscale systems with extensions L ∼ 10 2 km (see e.g. [1] for a review). We define the vertically averaged cloud cover as where the Heaviside function (q l ) = 1 for q l > 0 and is 0 otherwise. Similarly, we define the upward buoyancy flux as (see [26] for a similar definition for a convective-radiative flow) into a plane-time average and the rest. Figure 7 shows two instantaneous snapshots of the distribution of clouds (gray isosurfaces) and the related strong upward motion of fluid. Clouds are formed where the air rises up and can form a condensate. While the data for the smaller value of CW H in the upper panel display isolated clouds, the data in the lower panel display an almost closed cloud layer interrupted by fjord-like break-ups. Both cloud distributions are connected with a skeleton of upward motion that displays an arrangement in circulation cells with a diameter ∼4H and larger. It implies that the aspect ratio should be at least four.

Cloud size distribution.
The calculation of the cloud size distribution requires first identifying the subsets of grid points that form a cloud, a typical task of image segmentation. This can be done by complex multiscale algorithms that have been applied e.g. in [17] to study the maxima of the scalar dissipation rate. Here, we apply a so-called flood-fill algorithm. The result is shown in figure 8, where we compare the snapshot of Q(x, y) in the left panel with the separated clouds in the right panel. The left panel of figure 9 reports the analysis of the relative cloud size A as a function of the aspect ratio in a plot with double-logarithmic axes. The relative size A is defined as the ratio of the number of grid points that belong to a cloud to N x × N y . The cumulated density function (CDF) is plotted. We see in figure 9 (a) that the CDF is basically independent of , which confirms our results of figure 6(b). Figure 9(b) displays the dependence of the cloud size on the moist Rayleigh number. We found (not shown) that the algebraic decay of the PDF p(A) ∼ A −γ for  The horizontal bar in the picture has a length of L/2 6 and will be discussed in the geometric analysis of the height-averaged cloud patterns.

Box counting analysis of cloud fields.
Contours of Q can be considered as level sets and inspected by the box counting procedure in order to deduce a monofractal behavior of the turbulent mixing or not. We define the set of clouds as The biggest length scale that is associated with the plateau in the local slope is also indicated in the right panel of figure 8. Note that for the data in (a) and (b) we excluded clouds with an area smaller than 10 3 2 , where = L/N x is the grid spacing of the DNS. and the level set that corresponds to the cloud boundary as Note that in the case of the largest aspect ratio, we can monitor geometric properties over more than three orders of magnitude. We decompose the plane A = L x × L y = L 2 into boxes of size δ × δ with side length δ = 2 m L/N x and m = 0, . . . , log 2 N x . The box counting dimension is defined as where N (δ) is the number of boxes of side length δ that are necessary to cover the sets L 0 or L 1 . It is clear that this limit in definition (39) cannot be carried out for a finite grid resolution. In figure 10, we have summarized our findings. Instead of plotting N (δ), we show the local slope that is given by A plateau of the local slope would indicate an algebraic scaling of N (δ) and thus a monofractal behavior. In figures 10(a) and (b), we compare the data for the two runs at = 32; figure 10(a) is for L 0 and figure 10(b) for L 1 . A closed cloud layer would result in a constant local slope D loc = 2. The fjord-like break-ups of the cloud layer cause a slight decrease of the local slope, D loc < 2, for the smallest scales at CW H = −0.43. The isolated cloud case at CW H = −1.43 displays a significant decrease to a minimal local slope of 1.4. The minimum is found at scales of δ ≈ 2 −5 L. The level set analysis for the cloud boundary in figure 10(b) reveals a plateau of the local slope in the same range of scales, between δ = 2 −7 L and 2 −5 L, and at the same magnitude of 1.4, indicating a short range of scales with monofractal behavior. At smaller scales the data decrease to a local slope 1, which can be interpreted as (interrupted) lines. At larger scales, the local slopes converge to 2 since the level set L 1 becomes area-filling. We have also indicated the scale of the plateau in the right panel of figure 8. Similar to studies of high-Schmidt-number mixing of passive scalars in turbulence, we do not detect a strict monofractal behavior of level sets [9,33,34] for the box counting method.

Perimeter-area analysis of cloud fields.
A second method was performed in [6] with satellite image data and in [39] with data from LES of cumulus convection. The perimeter of a cloud c is therefore related to the square root of its area content √ A, which can be considered as a 'cloud diameter'. Figure 11(b) illustrates how the perimeter (black cells) is counted for a cloud (blue cells) following [11]. Our results for isolated cloud patterns at Ra M = 1.9 × 10 6 and 1.9 × 10 8 are shown in figures 11(c) and (d) in a double-logarithmic scatter plot. We can identify two subgroups in the data sets following different algebraic power laws. For the smallest clouds, the power-law exponent of c ∼ √ A D p is found to be D p ≈ 0.94-0.95. The largest clouds seem to follow a 4/3-scaling for all the analyzed data, which agrees well with D p = 1.32 as found by Siebesma and Jonkers [39]. Our result for D p is also within the range of scaling dimensions that have been found in [6]. Note also that the perimeter dimension can vary with the cloud type as has been discussed in [21].
In the light of percolation processes in the plane, the 4/3-scaling is considered as a universal result that follows independently of the underlying physics [11,30]. The crossover to the 4/3-scaling for the largest clouds can thus be interpreted here as independent of the specifics of the underlying turbulent flow that advects and forms the liquid water content. Interestingly, this crossover coincides quite well with the Taylor microscale of the turbulence, which we calculated as λ = √ 5ν/ u 2 i 1/2 V,t (see the dashed vertical line in figures 11(c) and (d)). This scale is of the order of the mean curvature of the vortex filaments in the flow. It is also considered as the Markov-Einstein coherence length above which velocity increment statistics can be described well as a Markov process [20]. This interpretation would also rationalize why in [39] a 4/3-scaling is detected over the whole data set: LES will probably not resolve λ.
While the cloud size distribution depends on the Rayleigh number, this geometric relation between perimeter and area is found to be rather insensitive. It can be rationalized as follows. When the Rayleigh number (and thus the Reynolds number) increase, similar flow patterns on smaller spatial scales will be generated that stir and generate local level sets of the liquid water content in the same way, but on smaller scales. This dynamical process is shifted consequently to smaller scales as figures 11(c) and (d) suggest.
We can go one step further and define a quantity that is the ratio of both quantities, the so-called shape complexity c/ √ A [5]. For any circle, it gives 2 √ π. Recent analysis in dry convection [42] suggested that the shape complexity follows a log-normal distribution. The present results for the shape complexity of the clouds do not support this finding as can be seen in the right panel of figure 11(a). Recall that the cloud fields are structures that arise from the competition of two scalar fields via equation (8), which is probably one reason for the difference.

Intermittency of the upward buoyancy flux.
Finally, we will discuss the trends of the buoyancy flux with Rayleigh number. In [36], we showed that the buoyancy flux is bounded between the fluxes of a completely unsaturated layer and a completely saturated layer, This suggests to define a dimensionless saturation fraction which is bounded between 0 and 1. In figure 12(a), we plot F(z) and demonstrate that crossover from a fully unsaturated to a fully saturated slab is sharper with increasing moist Rayleigh number. This is a manifestation of the reduced diffusivity present in the layer. Figure 12(b) shows the same quantity for the case of isolated clouds. The increase from unsaturated to partially saturated is again sharper, but superposed with an overall drier layer, and thus a confirmation of the results from figures 5(b) and 9(b).
The intermittency of the upward buoyancy flux B can be quantified by a multifractal analysis that generalizes the box counting. While the box counting analysis (see section 3.4.2) counts the number of cubes (or squares) of a particular side length δ (or r/L) that form a cover of a (level) set, the multifractal analysis assigns a weight-in mathematical terms, a measure µ i -with each box i of the cover. This is, for example, important if one wants to count how frequently a trajectory passes the boxes of a cover of an attractor or, as in the following, to quantify the spatial intermittency of a scalar field over a plane. Since the support of the field is the whole plane here, D 0 = 2. The spectrum of generalized dimensions D q for q real is defined as [12,15] Definition (39) follows for q = 0. In the present case of the upward heat flux, the measure that is accumulated in squares i of side lengthr = r/L is given by the following normalized integral: A (i) r is for the present case the area content of the ith square segment of side lengthr centered aroundx i and it holds A = ∪ N (r ) i=1 A (i) r . Similar analysis is conducted for kinetic energy [10] or scalar dissipation rate fields [33] or for the heat flux in radiatively driven convection [26]. The trend with the moist Rayleigh number is shown in figure 13. The three snapshots in figures 13(a)-(c), which are taken at three different Rayleigh numbers and the same CW H , illustrate the increasingly intermittent nature of the upward buoyancy flux. The contours of B(x, y) are increasingly filamented and the spatial regions of strongly enhanced flux narrower. Quantitatively, this manifests in the spectra D q as follows: it was discussed in [33] that increasing spatial intermittency is in line with a spectrum of generalized dimensions ever closer to the monofractal behaviour, namely that D q → D 0 = 2 for all powers q. The fourth data set at Ra M = 1.1 × 10 8 (run 9 at CW H = 0.17) has been added in figure 13(d). Since more liquid water is present in the case of a closed cloud layer, more pronounced maxima of B are obtained in comparison to figure 13(c). However, the multifractal spectrum collapses almost perfectly with the data of Ra M = 1.9 × 10 8 (run 10 at CW H = −1.43) as can be seen in figure 14.
It implies that the present measure of the intermittency is insensitive to whether the cloud layer is closed or not.

Summary and outlook
We have conducted three-dimensional DNS of moist turbulent convection in the Boussinesq approximation with a piecewise linear thermodynamics of phase change at the gaseous-liquid phase boundary. The simulation volume is a Cartesian cell with periodic side walls and freeslip top and bottom plates for aspect ratios up to 32. The assumption of local thermodynamic equilibrium prevents the fallout of rain as well as supersaturation. The latter is small, but present in real clouds and might be a central element connecting microphysics and turbulence in clouds [37]. The Boussinesq approximation limits the present studies to the case of shallow convection. We cover a subdomain of the five-dimensional parameter space only: both the dry and moist buoyancy fields are unstably stratified and all air parcels at the bottom plane z = 0 are neither subsaturated nor saturated. The focus is on the mixing and geometric properties of the turbulence. The main results that complement earlier studies [28,36] can be summarized as follows: • The relaxation to a well-mixed turbulent state, in which the dry and moist buoyancy fields are highly correlated, is slower when the Rayleigh number grows.
• The statistics of the mixture fraction χ , which is defined here by means of the buoyancy fields, becomes increasingly non-Gaussian with growing Ra M . The PDF of χ is most asymmetric for the case in which the largest fraction of the layer is neither unsaturated nor saturated, for CW H = −0.43. This is the case where the strongest asymmetry of the vertical root-mean-square fluctuations has been found [28].
• Similar to scalar mixing, the boundaries of isolated clouds show no strict monofractal behaviour although a local algebraic scaling with an exponent of about 1.4 for an intermediate range of scales is detected.
• The perimeter-area analysis in broken cloud layers agrees well with the results in large eddy simulations of cumulus convection.
• The multifractal analysis of the upward buoyancy flux confirms increasingly intermittent transport as the Rayleigh number grows.
Although the idealized moist Rayleigh-Bénard problem does not include several physical processes such as radiation, precipitation or the diurnal cycle, it still exhibits many features of atmospheric convection. Despite the many simplifications in the derivations of our model, it does produce a transition from the stratocumulus regime to the cumulus regime, qualitatively similar to those noted in [40] in a more realistic setting, which is of great interest given its impact on the Earth's energy budget. Our analysis indicates that such a transition is an intrinsic feature of moist turbulence and can be quantified in terms of a few non-dimensional parameters.
Our study also found that the simulated cloud field exhibits a sensitivity to the Rayleigh number, at least up to Ra ∼ 10 8 , which raises the question of whether the regime transition can be adequately captured with LES that rely on turbulent closure. A better understanding of the regime changes in the moist Rayleigh-Bénard convection, in particular for very high Rayleigh numbers, is necessary and could provide the much needed insights into the behavior of atmospheric convection.
Our future efforts will be focused on conditional instability, characterized by a stable dry Rayleigh number Ra D <0 and an unstable moist Rayleigh number Ra M > 0, and on obtaining a Lagrangian perspective on moist convection, as recently applied to the case of dry convection [35]. This could shed new light on the local buoyancy flux transport as well as on the path of air parcels through the clouds. Finally, further extension of the present model to include some of the missing physical processes, such as radiative cooling effects, is necessary to fully address the question of whether the regime transitions in the idealized moist Rayleigh-Bénard problem can be robustly extended to the Earth's atmosphere.