Probability density functions in the cloud-top mixing layer

The cloud-top mixing layer is an idealized configuration often employed in the literature to study local aspects (over length scales of the order of 10 m) of the top of stratocumulus-topped mixed layers. Latent heat effects are further investigated here by means of direct numerical simulations, discussing the probability density functions of the horizontal and vertical velocities, as well as the mixture fraction (equal to a normalized enthalpy and total water-specific humidity). The focus is on the turbulent convection layer that develops from the buoyancy reversal instability as a consequence of the evaporative cooling at the upper cloud boundary. An approximately self-similar behavior is found, where the convection scales based on the molecular buoyancy flux at the cloud top characterize the distributions at different times, at least to leading order and within the statistical convergence achieved in the simulations. However, a very strong vertical variation in the density functions across the turbulent convection layer is found, which is of relevance to possible models. Non-Gaussian behavior is often observed, even in the horizontal component of the velocity vector. In particular, large values of skewness and flatness are measured at the lower end of the turbulent zone, where external intermittency is very strong.


Introduction
Turbulent free convection is a state of a fluid system very common in nature [1]. Well-known examples are the diurnal cycle of the planetary boundary layer caused by surface heating or the upper ocean mixing layer due to nighttime surface cooling-both are radiatively driven. Another example is due to water phase change at cloud boundaries, which leads to evaporatively driven configurations. The latter is the topic of this paper, namely, the turbulent convection layer that develops on the top of stratocumulus clouds when buoyancy reversing conditions are encountered.
The understanding of stratocumulus, particularly in marine environments where the surface albedo is otherwise low, is important because marine stratocumuli play a fundamental role in setting the radiative balance of the earth system. Different models have been proposed since the original work by Lilly [2], but all hinge on the treatment of the entrainment interface: a narrow, stably stratified mixing layer that separates the underlying turbulent cloud from the free atmosphere above [3]. Mixing, or entrainment, of environmental fluid into the underlying stratocumulus determines the evolution of the whole boundary layer and thus it is critical to the subsequent evolution of that same cloud deck [4]. However, because such processes occur on scales and in environments that are difficult to measure or simulate, our understanding of entrainment processes and the dynamics of the entrainment interface remains primitive.
Field measurements by Caughey et al [5] most clearly demonstrate that the stratocumulustopped boundary layer is predominantly forced by the thermal radiation from a layer some tens of meters thick at the top of the cloud. This net emission cools the layer and drives the convective motions. However, pioneering studies by Lilly [2] have already realized that much finer scale processes, associated with small-scale mixing directly at the cloud-top interface, may also play a role in destabilizing the system. The absorption of latent heat (enthalpy of vaporization) as a result of local small-scale mixing between the cool and moist (condensate-laden) cloud below and the dry and relatively warm air above can, under certain circumstances, dominate thermal heating. In such situations, mixtures of the two fluids can exhibit net cooling, thereby becoming heavier than the unmixed environment (i.e. negatively buoyant), a situation known as buoyancy reversal [6]. It has been hypothesized that in such cases the whole stratocumulustopped boundary layer becomes unstable, as mixing leads to more mixing and hence the rapid dissolution of the cloud deck [7]. To the extent that such processes are important, it is thought 3 that they may be decisive in determining the distribution of stratocumulus, and hence the planetary albedo, on much larger scales.
In addition to the frequent reviews of this process from the point of view of the dynamics of the planetary boundary layer as a whole, where typical length scales are of the order of 1 km, Siems et al [6] and co-workers proposed idealized small-domain two-layer configurations to isolate and study specific cloud-top mixing processes. This line of inquiry has recently been continued by Mellado et al [8]- [10]. In the first of this last series of works, the Boussinesq formulation based on a mixture fraction variable usually employed in the literature is rationalized for this particular configuration, the so-called cloud-top mixing layer, in trying to understand clearly the assumptions implicit in such an approach in order to know its limitations. In section 2, we review its basic aspects. Molecular mixing is at the core of the problem because, as explained before, the relative strength of the molecular transfer of heat and mass between the cloud and the cloud-free regions of the system may lead to buoyancy reversing conditions at the cloud boundary, and then the corresponding buoyancy-reversal instability leads to a turbulent state. However, linear stability analysis [9] also suggests that the instability may not be strong enough to break up the cloud, as conjectured by other authors in the past. The reason is as follows. The linear analysis demonstrates that there exist two modes in the problem, one corresponding to the stable interfacial gravity waves of the inversion, defined by a buoyancy difference b 1 > 0, and a second one due to a relatively thin layer with characteristic buoyancy b s = −Db 1 , where the non-dimensional buoyancy-reversal parameter D is defined. When D > 0, the Randall-Deardorff criterion, this second mode is unstable and the buoyancyreversal instability issues. It is shown that the ratio of the time scales associated with those two modes, respectively, is proportional to √ D, and according to measurements, D 1 in typical atmospheric conditions. This means that the inversion returns to the equilibrium position quickly compared with the time that the heavy mixture below needs to move downwards a distance comparable to the perturbation length scale, and that is the reason for the cloud top not being broken. This hypothesis is confirmed by the direct numerical simulations presented by Mellado [10], where the structure and evolution of the system within the self-preserving turbulent regime are explained. We describe its main features in section 3.
The present paper adds to this previous body of work by providing further information about the turbulent convection layer of the system with an emphasis on the probability density functions (p.d.f.s) of the velocity vector and scalar fields. That is done in section 4. It is anticipated that this analysis, which is the first one we are aware of to explore in detail the structure of the turbulence at the cloud top interface, will provide the basis for the development of better models of the interfacial turbulence and for further improving our understanding of entrainment processes at cloud boundaries.

Formulation
The cloud-top mixing layer consists of a layer of warm and unsaturated moist air on top of a second layer, cool and saturated (condensate-laden). A two-fluid formulation based on equilibrium conditions of water leads to a description of the flow in terms of one conserved scalar, the mixture fraction χ. This scalar is equal to the normalized total water mass fraction and enthalpy under the low Mach number conditions typical of the system. It can be shown that the major assumptions for (and therefore limitations of) this approach are (a) liquid-phase continuum, (b) local thermodynamic equilibrium and (c) liquid-phase diffusivity equal to the vapor and thermal diffusivities. Estimates indicate that the conditions required for the validity of these hypotheses are not generally met at the top of stratocumuli, the main reasons being that droplets, with a number density below 1 mm −3 and a typical diameter of 10 µm, are too scarce and too large for that two-fluid approximation. Nevertheless, this physical model has provided great insight during recent years into latent heat effects at the boundary between cloud and cloud-free air. More details can be found in Mellado et al [8] and references therein.
The lower layer (cloud) is chosen to correspond to χ = 0, and the upper layer (cloud-free) corresponds to χ = 1 (see figure 1(a)). The frame of reference is set with the axis Oz along the vertical direction and pointing upwards, perpendicular to the two horizontal layers parallel to the plane x 1 O x 2 . The system is statistically homogeneous inside the horizontal planes, and the data inside these planes are used to construct the different statistics, statistics that depend then on the vertical position z and the time t.
The governing equations are The velocity vector is v = (v 1 , v 2 , w), with w being the vertical component, the kinematic viscosity is ν, κ is the scalar diffusivity, p is a kinematic pressure and k is the unit vector along the vertical Oz. Only the case of Prandtl number ν/κ equal to one is considered in this work.
The equilibrium assumption provides the buoyancy mixing function b e (χ) and closes system (1). The scalar χ gives the local values of enthalpy and total water content, and, once a reference value for the thermodynamic pressure is prescribed, the equilibrium state can be computed, and in particular the density and therefore the buoyancy. Figure 1 shows the buoyancy mixing function as obtained from the equilibrium calculation corresponding to the research flight RF01 of the second Dynamics and Chemistry of Marine Stratocumulus (DYCOMS-II) measurement campaign [4,11], which is the reference configuration to be discussed in this paper. Linear thermodynamic analysis shows that this function can be very well approximated by a piecewise linear profile, as observed in that figure, and we choose to write it as This functional relation is defined in terms of the saturation mixture fraction χ s and the buoyancy-reversal parameter, where b s = b e (χ s ). There is a first branch χ < χ s in which b varies between 0 and the buoyancy b s in saturation conditions χ s and which corresponds to the cloud. There is a second branch χ > χ s where the buoyancy increases up to b 1 , the strength of the stable inversion due to the density difference between the lower and upper horizontal layers. These two branches are finally smoothed over an interval δ s in mixture fraction space due to smoothness requirements imposed by the numerical scheme. When thermodynamic conditions are such that there exists an interval of mixtures over which the buoyancy anomaly b (relative to the lower layer) is negative, as in figure 1, buoyancy reversal is said to occur. Then, a derived quantity of relevance is the cross-over mixture fraction, which partitions the field at each time into a negatively buoyant region {x : χ(x) ∈ (0, χ c )} and a positively buoyant region {x : The equations are discretized on a regular Cartesian mesh of size 2048 × 2048 × 1536 (with stretching in the vertical direction) and numerically solved using finite differences; sixthorder compact Padé schemes for the first-and second-order spatial derivatives and a fourthorder Runge-Kutta scheme for the time advancement. The Poisson equation is solved using Fourier decomposition along the periodic horizontal planes in order to reduce it to a set of one-dimensional second-order equations, which are solved using the above described schemes. No-penetration free-slip boundary conditions are imposed at the top and the bottom. Further details of the numerical algorithm, particularly regarding the treatment of phase changes and the representation of the equation of state, can be found in Mellado et al [8].

Description of the turbulent regime
The consequence of the nonlinear relation between χ and other thermodynamic variables, like the buoyancy in (2), is that the two-layer structure in terms of the conserved quantities implies a three-layer structure in terms of the buoyancy, which is unstable if D > 0, as in figure 1(a). A temporally evolving turbulent flow develops from this instability, and figure 2 depicts that turbulent state.
The vertical structure observed in that figure consists of a relatively thin inversion layer of constant thickness h, dominated by molecular transport, and a turbulent convection layer that grows below, dominated by turbulent transport. The separation of both is well represented by the height of maximum mean stratification (maximum mean buoyancy gradient), a position that corresponds to z = 0 in figure 2 and that we refer to as the inversion base.
Analysis of the problem shows that this turbulent regime is controlled by molecular processes at the inversion base [10]. In particular, the most important quantity is the reference buoyancy flux, written here in terms of the entrainment velocity w e , the rate of displacement of the inversion upwards, which can be expressed as w e = 0.1 f 1 (κ|b s |χ 2 c ) 1/3 in terms of the parameters defining the problem, where the non-dimensional function f 1 (D, χ s ) is of the order of 1. This result stems fundamentally from two observations. Firstly, the saturation surface χ(x, t) = χ s , which is where the condition b = b s is imposed, remains very close to the inversion base at z = 0. Hence, the problem is analogous to free convection below a cold plate, at least for the small values of D corresponding to atmospheric conditions. Secondly, the mean entrainment velocity w e , as calculated from the simulations, is approximately constant. This steadiness can be explained in terms of a dynamical balance between diffusion, on the one hand, and the motion induced by the buoyancy, on the other hand, in a way similar to that employed in classical free convection [1]. The thickness is then given by h = κ/w e .
It is mainly the turbulent zone that we want to further characterize in this paper. The turbulent convection layer thickens continuously in time downwards into the cloud and it is characterized by the outer or convection scales constructed with the reference buoyancy flux B s 7 and the length scale, obtained from the depth-integrated turbulent buoyancy flux B = w b . A self-preserving state is approached where the thickness of the turbulent region grows according to z * ∼ and the buoyancy fluctuation is of the order of b * = B s /w * = (B 2 s /z * ) 1/3 . The non-dimensional functions f 1 (D, χ s ) and f 2 (D, χ s ) are approximately independent of the buoyancy-reversal parameters within the range of interest D 1 and χ s 1 of the parameter space {D, χ s } that defines the problem when Pr = 1.
Both inversion and convection layers are separated by a relatively thin transition or buffer region of approximate thickness 0.1z * . This zone is dynamically very active, with a strong transfer of vertical to horizontal motion as the updrafts approach the inversion from below. This leads to the cellular pattern depicted in figure 3 in terms of the horizontal divergence ∂v 1 /∂ x 1 + ∂v 2 /∂ x 2 and the formation of sheet-like plumes, in agreement with observations and laboratory experiments of related problems [12,13]. Convergent regions are shown in blue and can be associated with downdrafts. The upwelling parcels of fluid impinging upon the interface, in red, are redirected horizontally by pressure forces, and broken down by mixing events, which fuel the narrow zones of convergence.
The turbulent region inside the cloud is characterized in figure 4 in terms of the kinetic energy and the vorticity magnitude, plotting the horizontal and the vertical components separately. The vertical profiles are presented in terms of the convection scales w * and z * and the approximate collapse of the curves at the three different times illustrates the selfpreserving behavior that is achieved after an initial transient. The large anisotropy of the flow is manifest, the intensity of the vertical velocity fluctuation w rms = √ w w being stronger than the horizontal counterpart u rms = v 1 v 1 + v 2 v 2 in the core of the turbulent region due to the turbulent buoyancy flux term B responsible for the generation of turbulent kinetic energy. Close to the inversion base, however, the transfer of energy to the horizontal mode leads to a clear maximum of the horizontal component inside the buffer zone, where the cellular pattern of figure 3 is found. Figure 4(b) shows the vorticity fluctuation, in terms of the root mean square (rms) of the vertical component, ω 3 ω 3 , and the horizontal one, ω 1 ω 1 + ω 2 ω 2 . Note that the symmetries and initial conditions of the problem imply zero mean velocity and zero mean vorticity at all times. The thickness of the turbulent region is of the order of 1.5z * -2.0z * below the inversion base at z = 0, although the velocity fluctuations extend further down due to the irrotational motion induced by the pressure field. It is observed that, inside the bulk region, the fluctuation intensities of the three components of the vorticity vector are comparable among them, since ω 1 ω 1 = ω 2 ω 2 ω 3 ω 3 , but a strong anisotropy develops again inside the buffer zone. The horizontal vorticity fluctuation is significantly stronger there due to the production term ∇ × (bk), a term that has only a horizontal component and that peaks close to the inversion due to the strong local gradient of the mixture fraction (and therefore b, recall figure 1). In fact, the profiles of the terms in the enstrophy transport equation (not shown) indicate that this baroclinic term dominates over the vortex stretching term inside the buffer zone, whereas this relation is reversed at the core of the turbulent region, where vortex stretching controls enstrophy production.
The scalar field inside the convection layer is presented also in terms of the rms in figure 5(a) and its gradient fluctuation intensity, in particular the mean profile ε χ of the scalar dissipation rate 2κ|∇χ | 2 , in figure 5(b). The fact that the saturation surface χ(x, t) = χ s remains close to the inversion base compared to the thickness z * implies that, according to (2), there exists a linear relation between χ and b over most of the turbulent convection layer, and the characteristic value of the mixture fraction inside the turbulent region is related to the convection scale b * by Scaling of the mixture fraction field with χ * yields the values of order 1 exposed in figure 5 within the turbulent region. At the upper end, for values of z/z * close to zero, the scalar fluctuation and its gradient show a conspicuous peak that is a consequence of the small oscillations around the inversion base, where the mean gradient is very strong, as observed in figure 2.
(a) The Reynolds number z * w * /ν achieved in the simulation discussed here is of the order of 5000, which corresponds to a mixing region depth 2z * 5 m for the atmospheric conditions measured in the flight RF01 of the DYCOMS-II campaign. The corresponding Rayleigh number z * 3 |b s |/(νκ) is about 0.5 × 10 9 and the Nusselt number is z * w e /κ = z * / h 24. Reynolds and Nusselt numbers show a power law dependence on the Rayleigh number with exponents 4/9 and 1/3, respectively, which corresponds to the classical limit of no mean wind reported in free convection problems and that is consistent with the open domain configuration investigated in this work (there is no external length scale imposed on the turbulence).
The convection Richardson number Ri * = z * b 1 /w * 2 is larger than 500 and increases with time as t 1/2 , which explains the inability of the turbulent layer to penetrate through the capping inversion. In this respect, however, the internal Richardson number Ri (I ) = hb 1 /w * 2 might be more interesting to discuss (see [14] and references therein). This variable is the ratio between the thickness of the inversion h and the vertical displacement of that inversion due to the kinetic energy of the turbulent motion beneath it; let us call it δ = w * 2 /b 1 . For Ri (I ) = h/δ 1, we have a thick inversion with wavy motion inside it, and that is the situation observed in figure 2, where Ri (I ) 25. However, as the system develops and more latent heat is released into the kinetic energy mode, Ri (I ) decreases as t −1 and there is a critical time beyond which a flapping regime is established, with a typical amplitude δ of the perturbation of the inversion larger than the inversion depth h. The reasonable question is then whether the self-preserving regime presented here, and therefore the scalings, change or not-the answer is not straightforward. One possibility is that, indeed, non-linear interactions among these waves change the mixing processes at the inversion and the rates are no longer controlled by molecular processes. However, on the other hand, the aspect ratio of these interfacial waves is δ/z * = 1/Ri * , that is, increasingly small (already 1/500 for the case in figure 2), meaning that the linear regime might prevail and, locally, mixing occurs as described above in this paper. Larger simulations would be needed to reach a definitive answer. It is noted, however, that the turbulent convection layer has to become relatively thick before this change in regime might happen. It can be shown that δ/ h ∝ (z * / h) 2/3 , and the prefactor 0.1 f 1 (χ c /χ s ) 2/3 D is a small number. For the case RF01 considered here, this prefactor is about 5 × 10 −3 , which means that δ h first when z * reaches values about 3000 times h, about 300 m for the atmospheric case considered here, where h 0.1 m. This is a typical stratocumulus thickness, which implies that the formulation employed in this study is no longer valid because the variation in the background reference pressure needs then to be retained in the equilibrium calculation, and a different analysis is needed anyhow (see e.g. [15]).
It also seems appropriate at this point to comment as well on the possible influence of the real conditions at the stratocumulus top being different from the hypotheses implicit in the mixture-fraction formulation, as presented in the beginning of section 2. The stronger departure probably comes from the gradient transport model used for the liquid phase. The turbulent velocity is large compared to the settling velocity, in particular for the case RF01, about 30 mm s −1 compared to about 3 mm s −1 , which suggests that the suspension is well maintained by the turbulent motion. However, droplet dynamics next to the inversion base, where velocity fluctuations drop to small values, might alter the details of the mixing process and therefore the scaling laws. This issue links directly with the question about the validity to model the disperse and dilute two-phase flow in terms of a two-fluid formulation, and it represents an important question to be addressed in future work. It can be argued, however, that the tendency would be to have a slower mixing rate of cloud droplets with the upper drier region, which implies less evaporation and therefore weaker turbulent motion-the cloud-top therefore remains unbroken. The lack of thermodynamic equilibrium would also tend to bias the results in this direction, since the latent heat is the only driving mechanism in this problem, but, again, the particular corrections that these phenomena might introduce in the scaling laws discussed previously remain to be investigated.

Probability density functions
In this section, we explore the p.d.f.s of the velocity components and the mixture fraction as a function of height z and time t. Throughout, as a matter of notational convenience, we do not distinguish between the random variable itself and the corresponding sample space variable. Values are scaled by their standard deviation, already presented in figures 4 and 5. The full p.d.f. for each variable (namely, v 1 , w and χ) is calculated at three heights: firstly, at z = −0.05z * , which is inside the buffer zone beneath the inversion at a level where the horizontal velocity fluctuations are a maximum; secondly, at z = −0.40z * , which is well within the convection layer at a level where the vertical velocity fluctuations are a maximum; and thirdly, at z = −1.00z * , where the intermittent character of the large scales is clearly observed in figure 2. In addition, a higher moment is presented as a function of height for each variable.

Vertical velocity
The p.d.f.s of w evolve self-similarly over the range of Reynolds numbers sampled by this temporally evolving flow, but are robustly non-Gaussian in a way that depends markedly on the non-dimensional distance from the interface, as measured by z/z * . The collapse of the curves in figure 6(a) at each height for different times indicates that the shape of the distribution remains approximately the same as the flow develops in the self-preserving regime considered in this work, meaning that P w (w; z, t) f (w; z/z * (t)). This suggests that higher moments, and not only the standard deviation, are well characterized by the magnitude of the convection velocity w * -at least within the accuracy allowed by the statistical convergence achieved with the domain size considered in the simulations. However, the degree to which one wishes to extend this proposition (to a statement of complete Reynolds number independence) must be tempered by the realization that the range of Reynolds numbers embraced in the interval of time shown in those figures is only a factor of 2.3 (about one order of magnitude in terms of the Rayleigh number), a moderate value.
The distributions of the vertical velocity are predominantly non-Gaussian, and departures from Gaussianity are most pronounced at the upper and lower ends of the turbulent region, as 12 shown by the profile of the skewness S w = w 3 / w 2 3/2 in figure 6(b) (recall that w = 0). The negative skewness in the lower part lends itself to a geometric interpretation, negative skewness corresponding to less frequent but more intense downdrafts. Hence, the tendency of the skewness to be negative below the buffer zone is consistent with our expectation of a flow accelerated downwards as a result of density anomalies arising from mixtures originating at the cloud top. The dominance of these plummeting structures is also seen in the profile of the turbulent buoyancy flux B = w b (not shown in this paper), which is positive all across the turbulent zone. The positive correlation between w and b implies accelerating plumes whose strong vertical stretching implies small areas of fast downdrafts and, by continuity, a large horizontal divergence leading to relatively broad areas of relatively slower upwelling.
This behavior differs from measurements in Rayleigh-Bénard cells, where the vertical velocity is reported to be near Gaussian except close to the solid boundary [16]. However, the comparisons between both configurations are limited because of the difference in boundary conditions, and because of the asymmetry and the external intermittency in the cloud-top configuration considered here.
In this respect, our interpretation is consistent with a study of convective flows by Moeng and Rotunno [17]. They explored the skewness profiles for Rayleigh-Bénard convection, and for convection driven by cooling at the upper boundary, upon reaching a state of stationarity. They find, similarly to Belmonte et al [16], that w is relatively Gaussian away from the boundaries in the Rayleigh-Bénard system, but that negative skewness predominates for the top-cooled configuration. Unlike in their system, in which the lower solid boundary forces S w to drop to zero there, the magnitude of the negative skewness here increases significantly as we move downwards within the turbulent convection layer. Although there is some scatter of the curves in figure 6(b) for z −1.5z * , which, we believe, simply reflects sampling limitations because in the outer boundary region only a few structures are retained in the simulation (see figure 2), the magnitude is clearly larger than that at z = −0.4z * , which is still deep within the turbulent layer and is where the maximum of the turbulent kinetic energy is found.
Another difference is that our simulations have a free boundary at the top of the turbulent region. In the relatively thin area where the acceleration of fluid away from the interface induces large horizontal velocities, the behavior of the skewness is reversed as a layer of S w > 0 becomes apparent in figure 6(b). The thickness of this layer of positive skewness is commensurate with the thickness of the transition or buffer zone, as identified in plots of the velocity variances. This zone of positive skewness is not found in any of the simulations performed by Moeng and Rotunno [17], irrespective of the velocity boundary conditions they imposed. It indicates that within that transition layer the strongest events are associated with upward motion. The dissipation of updraft energy within the transition zone, along with the fact that the turbulent buoyancy flux is close to zero there because the motion is predominantly in the horizontal direction, likely explains this thin region of positive skewness.

Horizontal velocity
The p.d.f.s of v 1 , one horizontal component of the velocity vector, also show evidence of selfsimilar evolution, although, unlike in terms of the vertical velocity, the turbulent regions are more Gaussian. This is demonstrated with the help of figure 7, wherein the density functions P v 1 (v 1 ; z, t) have been obtained at the same heights and times as those shown in figure 6. The two upper curves, which correspond to the transition layer beneath the inversion and the core of the turbulent convection layer, respectively, show a nearly Gaussian distribution. This behavior of the horizontal velocity has also been found in Rayleigh-Bénard convection [16,18]. However, further down inside the region of external intermittency, two exponential tails develop. These large excursions of horizontal velocity correspond to the strong negative vertical velocities of P w (w) in figure 6(a).
The effect of these long tails is manifest in the higher order moments. In figure 7(b), the flatness F v 1 = v 4 1 / v 2 1 2 is plotted relative to the reference value 3. We focus on the flatness because by symmetry the odd moments of the horizontal velocity must vanish. Below the inversion base and down to z = −0.5z * , the flatness is close to 3, which corresponds to the near Gaussian behavior observed in the two upper sets in figure 7(a). Below, the flatness increases markedly. This figure also shows that there is a very thin region at the inversion base z = 0 where the flatness also increases.
This increase in flatness both at the lower boundary of the turbulent layer and at the top of the transition layer can be explained in terms of external intermittency. Both regions demarcate zones of the flow where the turbulent plumes are being redirected by pressure forces arising from their interaction with the irrotational flow. The strength of the inversion concentrates this effect in a very thin region at the top of the turbulent convection layer, essentially the domes of the upwelling fluid, as illustrated in figure 8, where a section of the saturation surface χ(x, t) = χ s is shown.

Mixture fraction
Unlike the velocity components, the scalar field χ is bounded between 0 and 1. Inside the turbulent region, this scalar is small, smaller than the saturation mixture fraction χ s = 0.09 (see figure 8), and hence the values in the p.d.f.s that we analyze are close to the lower bound χ = 0 corresponding to pure fluid from the cloud layer.
Like in the case of the velocity fields, the density functions show evidence of selfpreserving evolution and vary mainly as a function of the normalized height z/z * (t). The  distributions of scalar fluctuation are plotted in figure 9(a) normalized by the rms values presented before in figure 5(a). While the distributions are, for the most part, independent of time (Reynolds number), there is some suggestion that the tails are tending toward an exponential distribution as time increases. Again, this interpretation should be tempered by the realization that to establish better statistical convergence would require larger simulations.
In terms of spatial variation, there is also a strong dependence of the distributions on the vertical position. This is evident in figure 9(a). The p.d.f. at z = −1.0z * (the lower set of curves) shows a strong peak at χ = 0, which simply quantifies the strong incursions of pure fluid from the lower layer into the turbulent region-the external intermittency. Then, a long tail forms, which seems to be well approximated, to a certain extent, by an exponential tail-increasingly so as the Reynolds number increases. Exponential tails in the p.d.f. of the temperature field have often been reported in the bulk region of Rayleigh-Bénard convection [16,18,19], although here there is a strong asymmetry imposed by the one-sided boundary conditions that define the problem and comparisons are limited, as already mentioned.
As we move upwards across the turbulent boundary layer to examine the p.d.f. at z = −0.4z * , the previous peak at χ = 0 seems to decrease in favor of intermediate values of the scalar, and in the upper curve the mode of the distribution is located at a non-zero value of χ ; there seems to be a tendency to bimodal distributions, which would restrict the number of possible models to be employed (e.g. a beta distribution is often used in the case of bounded scalars, but a well-known limitation of it is that it does not reproduce bimodal behaviors). This broad shape is similar to that observed by Belmonte et al [16] at the edge of, and reported by Verdoold et al [18] to occur inside of, the thermal boundary layer of the Rayleigh-Bénard cell, in spite of the different boundary conditions (the cloud top iso-surface χ(x, t) = χ s where the buoyancy is equal to b s is here fluctuating instead of being a fixed horizontal surface).
The vertical variation in the skewness of the distributions is plotted in figure 9(b) in terms of the self-similar variable z/z * . These curves confirm the strong variation across the turbulent region already identified before in terms of the velocity field, and properties of the p.d.f., in this case the skewness, show a very significant vertical variation due to the external intermittency at the lower and upper boundaries of the turbulent region.

Conclusions
The turbulent convection region that forms beneath the inversion in the cloud-top mixing layer as a consequence of evaporative cooling has been further characterized in terms of the p.d.f.s of the velocity and mixture fraction.
A self-preserving behavior is observed in these distributions, in the sense that the density functions remain constant in time for a fixed self-similar position z/z * when plotted in terms of the corresponding variance, and these variances are again just a function of z/z * when normalized by the convection scales: w * in the case of the velocity and χ * in the case of the mixture fraction. Further corrections in the form of the tails of the distributions could indeed occur as the Rayleigh number increases with time, but the interval of scales achieved in the simulation (< 10 m vertical thickness) and the statistical convergence obtained do not allow any definitive answer.
On the other hand, a very strong variation within this turbulent region is found. Inside the upper transition or buffer zone, where there is a strong transfer of vertical to horizontal motion imposed by the inversion on top, the horizontal component of the velocity is close to Gaussian, the vertical component is positively skewed and the mixture fraction is close to bimodal with a non-zero mode. Within the turbulent core, where the turbulent kinetic energy achieves its maximum, the horizontal component remains nearly Gaussian and the vertical component develops the negative skewness expected from the downward forcing imposed by gravity. The scalar field already displays a peak at χ = 0, corresponding to pure fluid that is being entrained from the lower cloud layer into the turbulent zone, and develops a long exponential tail that corresponds to fluid regions 'dropping' from the cloud top above, the so-called sheet-like plumes typical of these free convection configurations. The lower end of the turbulent layer is dominated by the external intermittency, which leaves a strong imprint in the distributions in the form of very large values of skewness and flatness factors.