## Abstract

This study investigates the dynamics of the subsiding shell at the lateral boundary of cumulus clouds, focusing on the role of evaporative cooling. Since the size of this shell is well below what large-eddy simulations can resolve, the authors have performed direct numerical simulations of an idealized subsiding shell. The system develops a self-similar, Reynolds number–independent flow that allows for the determination of explicit scaling laws relating the characteristic length, time, and velocity scales of the shell. It is found that the shell width grows quadratically in time, linearly with the traveled distance. The magnitude of these growth rates shows that evaporative cooling, in its most idealized form, is capable of producing a fast-growing shell with numbers that are consistent with observations of the subsiding shell around real shallow cumulus clouds: for typical thermodynamic conditions in cumulus clouds, a velocity on the order of 1 m s^{−1} and a thickness on the order of 10 m are established in about 2 min. This fits well within the typical cloud lifetime, suggesting that this idealization is an adequate framework for the analysis of relevant aspects in the subsiding shell associated with buoyancy reversal. It also indicates that the scaling laws derived here can be used to estimate the potential strength of a subsiding shell and the mean lateral entrainment associated with it, once an estimate of the local thermodynamical state of the cloud boundary is provided. It is shown that the dominant parameter of this system is the saturation buoyancy, whereas the effect of the saturation mixing fraction is minor.

## 1. Introduction

In parameterizations of (shallow) cumulus convection, the representation of mixing between cloud and environment plays a key role, usually in the form of entrainment and detrainment coefficients (e.g., Tiedtke 1989; Kain and Fritsch 1990; Neggers et al. 2009). Many currently operational parameterizations have a bulk approach, where the effects of the entire cloud field are accounted for by a few effective equations for the transport of heat and moisture. However, as operational models increase their resolution, this bulk approach is being replaced and clouds are represented on a more individual basis. Thus, understanding of the detailed structure of mixing events around single clouds becomes essential.

Lateral entrainment is the overwhelmingly dominant source of entrainment in shallow cumulus convection (Heus et al. 2008). The mixing process at the sides of the cloud is surprisingly complex because of a thin subsiding shell around the cloud. This shell has been found in observations (e.g., Rodts et al. 2003; Siebert et al. 2006; Heus et al. 2009b; Wang and Geerts 2010), and its origin has been studied using large-eddy simulations (LES; Heus and Jonker 2008). They emphasized that the descending shell is a result of the buoyancy reversal because of evaporative cooling; in other words, the shell is a result of lateral mixing between cloudy and environmental air. This view is supported by the observations of Siebert et al. (2006), who found increased turbulent mixing within this shell. Gerber et al. (2008), Jarecka et al. (2009), and Slawinska et al. (2012) showed that the shell has its impact on the microphysical structure of the cloud. Dawe and Austin (2011) argued that the entrainment rates that are typically used as a closure in large-scale models can only be validated through LES if one takes the preconditioning of the air by the shell into account. Similarly, Rio et al. (2010) were able to improve the performance of their boundary layer parameterization by including the transport through the shell, and Cooper et al. (2013) stress the importance of detailed knowledge of the entrained air for in-cloud processes, including the formation of rain.

Although the relevance of lateral mixing seems essential in understanding cumulus dynamics, the detailed structure of the cloud-edge region is not very clear. Heus et al. (2009b) showed that a spatial resolution of 25 m may be just sufficient to obtain a converging mass flux through the shell, but one would need a finer resolution to study the internal dynamics. Heus and Jonker (2008) extend the two-layer model of Asai and Kasahara (1967) to a three-layer model but need crucial quantities like the width of the shell as input parameters. It seems that LES-type resolutions do not suffice to fully resolve the mixing processes at the cloud edge and the full turbulence spectrum is relevant to understand this problem. In this study, therefore, we use direct numerical simulations (DNS) to completely resolve the turbulence structure of the subsiding shell.

With DNS, the full Navier–Stokes equations are solved down to the Kolmogorov scale, so that no subgrid-scale turbulence parameterization is necessary. With the current computational resources, about three orders of magnitude of scale separation can be resolved. Given a Kolmogorov scale typical for the atmospheric boundary layer, this would result in a domain size on the order of meters. While this scale is certainly relevant for cloud-edge mixing, the relevancy can be increased based on Reynolds number similarity. In such a regime, some important characteristics of the flow do not depend on the Reynolds number anymore, so, in other words, the Reynolds number is high enough (Tennekes and Lumley 1972). The existence of such a regime in a complex flow like the one being studied here needs to be shown in an empirical way—the results presented here do show such a regime.

A second aspect of the current limitation in computational resources is that DNS can only resolve the flow directly around the edge of the cloud. In other words, the structures of the cloud and of the far environment are boundary conditions to the setup and assumptions need to be made there. Also, by formulating the problem in a nondimensional form, the results of a single simulation can be rescaled to be valid across a certain range of atmospheric values. This step is precisely the opposite to LES, where the large scales are retained and the effect of the small scales is modeled, and underscores the complementary relationship between the two techniques. The approach used in the current study is comparable to an earlier analysis of local evaporative cooling effects at the stratocumulus cloud top (Mellado 2010). In that research, a Reynolds number–independent regime was found, and scaling laws were presented for the cloud-top entrainment rate because of the turbulence produced by evaporative cooling. The essential difference between that and this work is that, here, the interface between the mixing fluids is aligned with gravity instead of perpendicular to it.

In short, the goal is to study how turbulent mixing induces a subsiding shell, to find the most simplified flow that still bears the characteristics of mixing at the side of a shallow cumulus cloud, and to obtain scaling laws to estimate the relevance of those processes.

Of course shallow cumulus clouds are a complex phenomenon, and this paper will ignore many processes, ranging from microphysics and cloud shape to the in-cloud buoyancy, turbulence, and vertical velocity. Still, given that the shell is governed by the local reversal of buoyancy, it is interesting to see up to what level a highly simplified model can produce a shell with numbers that compare with observations, at least in order of magnitude. Ignoring these processes and allowing the shell to live in an idealized environment is likely to create an overestimation of the strength of the subsiding shell. The aim of this study is more to investigate whether a setup that is exclusively based on evaporative cooling is capable of explaining the subsiding shell as it was reported by Heus and Jonker (2008) and other studies.

The outline of the paper is as follows. Section 2 formulates the setup and discusses its strengths and limitations. Section 3 discusses the results in nondimensional form to derive general expressions for the characteristic buoyancy, velocity, and length scales explicitly as a function of the parameters of the problem. In the discussion in section 4, these results are applied to typical atmospheric conditions and compared to LES and observations.

## 2. Formulation

We are interested in obtaining the scaling laws describing, to leading order, the evolution of the subsiding shell under buoyancy-reversing conditions. For that purpose, we model the lateral mixing between the cloud and its environment as a two-layer system, with the interface aligned with gravity in the vertical direction Oz (see Fig. 1). For our current experiments, no additional external forcings are applied and only evaporative cooling is considered. This evaporative cooling induces a vertical descending motion and the cloud boundary starts to broaden in time because of turbulence. We assume that the dominant mixing occurs locally, far enough from the cloud base and the cloud top to neglect their influence. According to Heus et al. (2008), this is typically true for roughly 50% of the cloud. This means that the system is statistically identical as a function of height and, once the background hydrostatic equilibrium is subtracted in the equations and if the domain size is large compared to the autocorrelation length, the top and bottom boundaries can be set as periodic. For similar reasons, periodicity is also assumed in the spanwise horizontal direction Ox along the cloud boundary. Hence, the statistics can be calculated by averaging within the planes xOz, and they depend on the time *t* and on the crosswise coordinate *y* varying perpendicularly to the interface. The underlying assumption in this problem definition is that the size of the mixing layer is much smaller than the size of the cloud. While this assumption clearly breaks for long simulation times, it seems to be a reasonable approach on the time scales of one mixing event—the scaling laws derived later in the analysis confirm this assumption.

Assuming thermodynamical equilibrium, the state of the fluid particle is set by the enthalpy *h* and the total water content *q _{t}* together with the pressure

*p*. The thermodynamic pressure can be considered as constant over domains on the order of a few tens of meters. Such an approach has often been used in the investigation of evaporative cooling effects at the stratocumulus top (Albrecht et al. 1985; Bretherton 1987; Shy and Breidenthal 1990; Grabowski and Clark 1993; Mellado 2010). This methodology sacrifices the details of droplet dynamics in favor of a reduced set of equations that allow the simulations to reach higher Reynolds numbers. This provides a broad spectrum of well-resolved turbulence scales interacting with evaporative cooling effects. Assuming identical diffusivity for heat and moisture,

*q*and

_{t}*h*can be normalized and described in terms of a single scalar, the mixture (or mixing) fraction

The subscript 0 corresponds to reference minimum in-cloud values and the subscript 1 to the reference maximum environmental values within this idealized two-layer configuration, so that *χ* = 0 holds for pure cloudy air, and *χ* = 1 for pure environmental air.

In the Boussinesq limit, the governing equations are

The velocity vector is **v** = (*u*, *υ*, *w*), *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 direction *z*.

The buoyancy *b* = −*g*(*ρ* − *ρ*_{0})/*ρ*_{0} is provided by *b*(**x**, *t*) = *b ^{e}*[

*χ*(

**x**,

*t*)], where the mixing function

*b*(

^{e}*χ*) is obtained from the thermodynamical equilibrium calculation given

*χ*locally at (

**x**,

*t*) (i.e., the enthalpy and total water content) and a global reference thermodynamic pressure level. However, this function is very well approximated by

This is essentially the piecewise linear approximation of the buoyancy function *b ^{e}*(

*χ*) =

*b*min{

_{s}*χ*/

*χ*, (1 −

_{s}*χ*)/(1 −

*χ*)}, varying from

_{s}*b*= 0 at

*χ*= 0 inside the cloud, through

*b*at the saturation mixture fraction

_{s}*χ*=

*χ*, and back to the environmental value

_{s}*b*= 0 at

*χ*= 1 (see Fig. 1). To overcome numerical side effects associated with the high-order schemes, the function (2) slightly smooths the discontinuity of derivatives at

*χ*over a distance

_{s}*δ*in the mixing-fraction space small enough for not altering the results discussed in this paper. Further details on the bulk model discussed in this section can be found in Mellado et al. (2010).

_{s}One assumption that becomes obvious from Eq. (2) is that the inner part of the cloud is taken to be neutrally buoyant with respect to the environment. Furthermore, the vertical velocity is initially zero both at *χ* = 0 and *χ* = 1. These assumptions need to be put into context. As an example, one could look at conserved variable mixing diagrams (Paluch diagrams), for instance such as published by Neggers et al. (2002), republished in Fig. 2. Looking at the distribution of cloudy parcels in Fig. 2 (upper left part of the distribution, above the dashed saturation line), a significant amount of cloudy air is neither very buoyant nor has a strong vertical velocity. Presumably, this is the air outside of the cloud core, the air that is the most prone to being detrained. Therefore, in designing the simplest realistic setup, the vertical velocity and buoyancy are taken to be equal between *χ* = 0 and *χ* = 1. However, especially the assumption of zero velocity is a simplification that warrants some further research, should one be interested in using the setup beyond the current idealizations, and in a context that tries to capture all the complexity of the shallow cumulus cloud edge.

The initial condition is defined by a mean mixing fraction profile

and a mean vertical velocity following

Equation (3) also defines *δ*_{0}, which will be used to normalize some of the variables through this paper. A broadband perturbation is added to the velocity profile to accelerate the transition into the fully developed turbulent regime. This perturbation is defined by a spectrum proportional to *f*^{4} exp[−2(*f*/*f*_{0})^{2}], *f* being the spatial frequency, with a random phase. In physical space, it is restricted to a thickness 2*δ*_{0} around the maximum mean velocity using a Gaussian shape function.

The reason for such an apparently complex initial condition instead of, for instance, simply setting *w*_{0} = 0 and white random noise is to be able to ascertain the duration of the initial transient before the flow enters into the fully developed turbulent regime. We can do that by comparing simulations with different initial conditions and we have done it by varying {*δ*_{0}, *w*_{0}, *f*_{0}} in the setup described above. After this transient, for long enough time, the details of the initial perturbation are sufficiently forgotten and the parameters {*δ*_{0}, *w*_{0}, *f*_{0}} drop out of the analysis—in particular, the scaling laws described in the following section are characterized solely by {*ν*, *κ*, *b _{s}*,

*χ*}. Hence, the final set of nondimensional parameters is {

_{s}*ν*/

*κ*,

*χ*}. We do not explore Prandtl number effects and Pr =

_{s}*ν*/

*κ*is set to unity. Having used

*b*to nondimensionalize the problem, the rest of the thermodynamical properties are fully covered by varying

_{s}*χ*between 0 and 1. In fact, because of the symmetry of the problem, only the interval between 0 and 0.5 needs to be considered.

_{s}Looking again at Fig. 2, and specifically at the relative distances along the mixing line between the environmental profile, the saturation line, and the zero-buoyancy line, the value *χ _{s}* = 0.3 is taken as reference case. Configurations with

*χ*varying between 0.1 and 0.5 were also investigated, as summarized in Table 1. Since the effect of

_{s}*χ*is relatively small in all of the results discussed in this paper, the figures presented in the following sections only include the cases

_{s}*χ*= 0.1 and

_{s}*χ*= 0.5 in addition to the reference configuration

_{s}*χ*= 0.3, for the sake of clarity of those figures. Table 1, however, includes all the results.

_{s}### a. Numerical algorithm

The numerical algorithm consists of high-order, spectral-like compact finite-differences for the spatial derivatives, and a low-storage fourth-order Runge–Kutta scheme to advance in time. The lateral boundary of the computational domain in Oy is placed far enough from the turbulent region that develops in the center on the domain to avoid any significant interaction. No-penetration free-slip boundary conditions are used in this direction Oy, along with zero flux conditions for the scalar. Further details can be found in Mellado and Ansorge (2012).

### b. Setup of the simulations

The grid is 3584 × 1280 × 3584 for the reference case *χ _{s}* = 0.3. Grid stretching is used along the crosswise direction Oy so that the domain size is 3.5

*L*

_{0}× 1.6

*L*

_{0}× 3.5

*L*

_{0}, where

*L*

_{0}measures the width of the computational domain, and is therefore commensurate with the thickness of the subsiding layer at the final time of the simulation. Smaller domains with a grid size 2048 × 1024 × 2048 are used for the remaining cases studying the influence of the parameter

*χ*(see Table 1). This size was enough to quantify the effect of this parameter, and allowed us to concentrate computational resources in the reference case

_{s}*χ*= 0.3 to achieve better statistical convergence and higher Reynolds numbers.

_{s}The reference Reynolds number is set to 4 × 10^{4}, which guarantees a minimum resolution at the end of the simulation on the order of Δ*x*/*η* ≃ 1.5, Δ*x* being the grid spacing and *η* = (*ν*^{3}/*ɛ*)^{1/4} being the Kolmogorov scale, where *ɛ* is the mean viscous dissipation rate of turbulent kinetic energy. For reference, for a typical value *b _{s}* = −0.03 m s

^{−2}(evaporative cooling of about 1 K), the domain size is

*L*

_{0}≃ 2.3 m, having used a kinematic viscosity

*ν*≃ 1.5 × 10

^{−5}m

^{2}s

^{−1}.

The initial width is such that the cloud interface is resolved by 12 grid points, so that *δ*_{0}/*L*_{0} ≃ 0.0117. The initial velocity is . The broadband perturbation is defined by *f*_{0} = 1/(2*πδ*_{0}), which is within the envelope of linearly unstable modes for the Bickley jet profile (4) at a Reynolds number (Drazin and Reid 2004). The initial fluctuation intensity is , where *k* is the turbulent kinetic energy. The simulations are stopped when the pressure root-mean-square (rms) at the boundaries is about 2%–3% of the maximum within the turbulent region, to assure that the finite size of the computational domain along Oy has a negligible impact on the flow. The smoothing parameter *δ _{s}* in the buoyancy mixing function (2) is

*δ*≃ 0.1/16.

_{s}## 3. Results

Figure 3 depicts the structure of the flow. The left panel shows *q _{l}*

_{,0}max(0, 1 −

*χ*/

*χ*), which is a piecewise linear approximation to the equilibrium liquid water content similar to that in Eq. (2) for the buoyancy. The middle and right panels contain the buoyancy field

_{s}*b*and the scalar dissipation rate

*ε*

_{χ}= 2

*κ*|

**∇**

*χ*|

^{2}, respectively. The figure corresponds to the reference case

*χ*= 0.3 but the other configurations in Table 1 present a similar structure. The snapshot is taken at

_{s}*t*

_{2}≃ 33

*t*

_{0}, where the time scale

measures the duration of the spinup phase. In the following sections, we are only interested in the fully developed turbulent flow as depicted in Fig. 3, when *t* ≫ *t*_{0} and the precise initial conditions do not matter anymore.

The liquid water content in Fig. 3 exposes very clearly the cloud boundary, with large-scale structures and the typical filamentation of the interface itself. Comparing this field with the other two, it is also observed that the cloud comprises only half of the turbulent flow, approximately, and there is strong turbulent motion in regions with zero liquid water content. The cloud boundary coincides with the saturation surface *χ*(**x**, *t*) = *χ _{s}* in our bulk formulation and, consistently, the buoyancy is most negative for mixtures around that cloud boundary, as observed in the middle panel of Fig. 3.

The system behaves basically as a temporally evolving jet that accelerates downward and broadens under the effect of the negative buoyancy caused by the evaporative cooling. It is interesting to observe that, though the system is buoyancy driven, the shear instability that issues from the resulting velocity field confers many features typical of free shear flows, like jets or wakes (Pope 2000). We observe all the defining characteristics of turbulent mixing: turbulent entrainment, stirring, and molecular mixing (Dimotakis 2005). Large-scale structures exhibit occasionally a very conspicuous meandering of the cloud interface. This sinusoidal morphology is typical in other jetlike flows and entails a strong intermittency, the local signal alternating between turbulent (rotational) and nonturbulent (irrotational) states (Pope 2000).

Stirring is crucial in this system because it increases the extension of *χ*(**x**, *t*) = *χ _{s}* and this surface plays a determining role. This can be seen mathematically through the buoyancy transport equation (Mellado et al. 2009a)

with *S* being the source term

The factor *d*^{2}*b ^{e}*/

*dχ*

^{2}in the expression above implies that

*S*is only nonzero in a thin region around the saturation interface

*χ*(

**x**,

*t*) =

*χ*, given the piecewise linear approximation (2) for the buoyancy function

_{s}*b*(

^{e}*χ*). Physically, evaporation and condensation are still occurring inside the cloud (where

*χ*<

*χ*), but, under the assumption of thermodynamical equilibrium, this process is infinitely fast and molecular mixing fully determines the evaporation/condensation rate. Thus,

_{s}*κ*∇

^{2}

*b*completely describes the process. The dependence of

*S*on

*ε*

_{χ}= 2

*κ*|

**∇**

*χ*|

^{2}, depicted in the right panel of Fig. 3, emphasizes even further the relevance of the small scales for evaporative-cooling processes, in this case the small-scale molecular mixing. The typical lamellar structure of this field, with a sharp definition of the turbulent/nonturbulent interface and a local variation of several orders of magnitude over relatively short distances, is clearly observed in the figure.

### a. Time evolution of the characteristic scales

The previous description of the flow can be quantified in terms of a relatively small set of characteristic scales, namely, the mean minimum buoyancy *b _{c}*, the mean minimum velocity

*w*, and the width of the shell

_{c}*δ*. If the flow reaches a self-similar regime, this limited number of quantities should suffice to set several other quantities one might be interested in.

The magnitude *b _{c}* ≥ 0 of the minimum mean buoyancy provides the characteristic buoyancy. It is defined as

and it gives a measure of the strength of the evaporative cooling. Angle brackets indicate a mean value, calculated by averaging within the statistically homogeneous, vertical planes xOz. The second equality in Eq. (8) defines *y _{b}*(

*t*), the lateral position of that minimum. As observed in Fig. 4, all cases tend toward an approximately constant behavior beyond

*t*

_{1}≃ 15

*t*

_{0}–20

*t*

_{0}. Results from additional simulations with different initial conditions (not shown) approximately coincide with these curves beyond that same interval of time, which is further indicative of a fully developed regime being established beyond

*t*

_{1}(Tennekes and Lumley 1972; Monin and Yaglom 2007). The corresponding values

*b*/|

_{c}*b*| calculated as a mean over that final interval of time vary within the range 0.65–0.71 (see also Table 1). The result that

_{s}*b*is a fraction of the magnitude of the saturation buoyancy |

_{c}*b*| is to be expected because the buoyancy field

_{s}*b*(

**x**,

*t*) varies between 0 and

*b*, by definition, and turbulence brings parcels of fluid with a buoyancy magnitude less than |

_{s}*b*| to any given lateral position

_{s}*y*. The constant behavior can be interpreted in terms of the evolution equation for

*b*,

_{c}obtained by taking the time derivative in the definition (8), using the mean of the transport equation (6), and noting that ∂〈*b*〉/∂*y*[*y _{b}*(

*t*),

*t*] = 0. The prime indicates a turbulent fluctuation. The equation above shows that the steady behavior of

*b*implies a balance between the turbulent transport and the mean source term 〈

_{c}*S*〉, provided that the Reynolds number is large enough for the mean molecular term to be negligible.

A characteristic thickness can be defined as

where the integral extends across the complete domain. Figure 5 plots the temporal evolution of this thickness normalized with *δ*_{0}. Several features of this evolution are readily noticed. First, like for the minimum buoyancy, there is a transient during the first 10*t*_{0} followed by a well-defined behavior beyond *t*_{1} ≃ 15*t*_{0}–20*t*_{0} until the final time of the simulation *t*_{2}. Second, the growth rate *dδ*/*dt* increases during this second interval of time. Third, the shell has broadened by more than a factor of 2 between *t*_{1} and *t*_{2} for the reference case *χ _{s}* = 0.3, attaining values of

*δ*much larger than

*δ*

_{0}.

We will see below that the length *δ* provides the thickness of the shell not only in terms of the buoyancy but also in terms of other statistical properties like the mean mixture fraction 〈*χ*〉 and the mean vertical velocity 〈*w*〉. It can be argued that the mean mixture fraction (i.e., the normalized total enthalpy and total water content) is then fully characterized because it simply varies monotonously between 0 and 1. However, we still need a measure of the velocity scale. For this, we use the magnitude of the minimum vertical velocity,

As in the case of the buoyancy, the second equality defines the lateral position *y _{w}*(

*t*) of the minimum and the minus sign is introduced in the definition above to obtain

*w*> 0. As seen in Fig. 5, this velocity measure also shows an initial transient and a subsequent well-defined behavior, in this case, an apparently linear evolution in time, and it becomes large compared to the velocity scale associated with the initial condition [see Eq. (4)]. The evolution in time can also be analyzed in terms of the corresponding transport equation

_{c}The mean buoyancy is the source term in this equation, though it is noted that the magnitude of 〈*b*〉[*y _{w}*(

*t*),

*t*] can be less than

*b*, defined before in Eq. (8), if

_{c}*y*does not coincide with

_{w}*y*, which is actually the case if

_{b}*χ*≠ 0.5, as shown later. The turbulent transport partly compensates this forcing but, unlike in the case of

_{s}*b*, Fig. 5 shows that this balance is not complete and the shell continuously accelerates. This is a direct consequence from the setup of the experiment: the producing mechanism (evaporative cooling) is included, but the infinite cloud assumption removes all restraining mechanisms that could bound the width and depth of the shell. From Figs. 4 and 5 it can be observed that

_{c}*w*grows linearly, that

_{c}*δ*grows quadratically, and that the buoyancy minimum becomes constant in time. These are important findings and their significance is discussed in section 3c.

We conclude this section by noting that the Reynolds numbers

attained in the simulations, on the order of 10^{4}, are large enough for Reynolds number similarity to be observed. One manifestation thereof is the inviscid scaling of the viscous dissipation rate , where *s _{ij}* is the rate-of-strain tensor. This means that

*ɛ*scales proportionally to for large-enough Reynolds numbers (Tennekes and Lumley 1972). Figure 6 depicts this property of the flow in terms of the temporal evolution of the integral of the self-similar profile as a function of Re

_{c}. A plateau is clearly observed beyond 2–4 × 10

^{3}, approximately. The separation of scales quantified in terms of the ratio between the thickness and the minimum Kolmogorov length scale

*η*= (

*ν*

^{3}/

*ɛ*)

^{1/4}at the end of the simulation in case

*χ*= 0.3 is about 350 (see Table 1). This relatively large-scale separation suggests again that the current DNS studies on moderate Reynolds numbers are applicable to the high-Reynolds-number atmosphere.

_{s}### b. Self-similar profiles

Figure 7 plots the profiles for the mean vertical velocity and mean buoyancy at two times *t*_{1} and *t*_{2} normalized with the characteristic scales discussed in the previous section and as a function of the variables

respectively. Despite the fact that the shell has thickened by more than a factor of 2 within that time interval, it is clear that the normalized profiles collapse on top of each other. This result indicates that the flow develops self-similarity and the mean profiles can be written as

where *f _{w}* and

*f*are solely a function of the so-called self-similar variables in Eq. (14) and independent of time (Tennekes and Lumley 1972; Monin and Yaglom 2007; Pope 2000).

_{b}It is interesting to notice that there exists a small shift between the two profiles. The shell does not only broaden in time but also migrates sideways into the cloud for the asymmetric cases *χ _{s}* < 0.5. The lateral positions of the velocity and buoyancy minima

*y*and

_{b}*y*, defined by Eqs. (8) and (11), respectively, drift toward negative values, as shown in Fig. 8. The buoyancy minimum shifts farther into the cloud than the velocity minimum.

_{w}One possible physical interpretation of this shifting of buoyancy minimum relative to the velocity minimum is that the thermodynamics impose the asymmetry in the buoyancy field, which is the driving mechanism of the system and pulls down the shell off center. The mean velocity follows but, at the same time, the turbulence homogenizes the momentum excess, tending to shift the maximum velocity toward the exterior of the shell. These two competing processes were already identified in section 3a, Eq. (12), being quantified by the mean buoyancy term and the divergence of the turbulent transport, respectively.

Consistent with the self-similar behavior identified before, the displacement is completely determined by the local *δ*, so that the ratios *y _{w}*/

*δ*and

*y*/

_{b}*δ*become constant in time after the initial transient. In particular, the distance

*y*−

_{w}*y*becomes a constant fraction of the shell width and relates both self-similar variables by

_{b}It is at its maximum for *χ _{s}* = 0.1, the most asymmetric case, with a value on the order of (

*y*−

_{w}*y*)/

_{b}*δ*≃ 0.1 (see Table 1).

Second-order moments display also a self-similar behavior, but later in time. This is shown in Fig. 9, where the nonzero components of the Reynolds-stress tensor are plotted in self-similar variables. The variance , that is, the turbulent fluctuations along the horizontal direction inside the shell, is the slowest in achieving self-similarity and at the early time *t*_{1} = 15*t*_{0} (not shown) is still about 20% smaller than its final value. However, beyond *t* = 20*t*_{0}–25*t*_{0}, all of the profiles approximately collapse on top of each other. It is also interesting to note that the numerical values obtained here in the subsiding shell agree with results from other free shear flows, like jet flows (Pope 2000). Also, Fig. 9 shows that the asymmetry of the system is more clearly exposed in these statistics, with a slightly stronger variance toward the cloud side, the side where the buoyancy flux peaks.

### c. Scaling laws

The observed self-similarity can be used to obtain the scaling laws for *δ*(*t*), *b _{c}*(

*t*), and

*w*(

_{c}*t*). By substituting the definitions in Eqs. (14) and (15) into the mean transport equation for the vertical velocity, it can be shown that a self-similar regime can exist if the ratios

are constant. (See the appendix for a full derivation.) This is confirmed in Fig. 10. This means that the local scales *b _{c}* and

*δ*determine the dynamics within the self-similar regime and the details of the initial conditions are approximately forgotten, which is one of the main features of this type of turbulent regime (Tennekes and Lumley 1972; Monin and Yaglom 2007). Indeed, the first two terms imply that the acceleration is set by the buoyancy

*b*, and that the resulting velocity scale is proportional to , which resembles the free-fall law for a given reduced gravity

_{c}*b*—though we still need to find an explicit expression for

_{c}*b*(

_{c}*t*). The third term tells us that the growth rate

*dδ*/

*dt*is proportional to the maximum mean velocity difference between the shell and the environment. This is a typical relation in free shear flows, since that velocity scale

*w*is also the velocity scale of the large-scale engulfment processes responsible for the mean rate of turbulent entrainment into the shell (Pope 2000).

_{c}It is also shown in the appendix that the self-similar parameters satisfy

for any value of *χ _{s}*. This constraint implies in turn that

Consequently, the condition in Eq. (17) that the group remain constant within the self-similar regime implies that the buoyancy scale *b _{c}* itself remains constant. This is in agreement with what we saw earlier in Fig. 4. From Eq. (17), we obtain similarly that

*w*grows linearly in time and

_{c}*δ*grows quadratically in time. Again, these behaviors are in agreement with Fig. 5.

We can now derive explicit expressions for these scaling laws. Equation (18) shows that only two self-similar parameters are independent. It proves convenient to formulate the analysis in terms of

both a function of *χ _{s}* solely. Then, Eq. (18) provides the relation

and we obtain

explicitly as a function of the thermodynamic parameters *b _{s}* and

*χ*. If an initial condition

_{s}*δ*

_{1}at a time

*t*

_{1}inside of the self-similar regime is available, we can integrate the last equation to obtain explicitly the thickness at a later time according to

This quadratic growth law is included in Fig. 5 and the good agreement is evident.

Values of the constants *c*_{1}(*χ _{s}*) and

*c*

_{2}(

*χ*) are summarized in Table 1 for the different configurations. The dependence on

_{s}*χ*is negligibly small for all of the cases

_{s}*χ*≥ 0.1 studied here and typical values are

_{s}*c*

_{1}≃ 0.25 and

*c*

_{2}≃ 0.1. The differences observed in Fig. 5 are due to differences in the early nonlinear transient. Note, however, that

*χ*does not drop out of the problem completely and that it affects the system in terms of a migration of the shell into the cloud and asymmetries in some self-similar profiles, for example, through the difference in location between the buoyancy and the vertical velocity minima.

_{s}We conclude this section by explaining the particular choice of *c*_{1} and *c*_{2}. The first parameter measures how turbulence reduces the dynamical effect of *b _{s}*. If the flow would be laminar, then

*b*≃ |

_{c}*b*| and

_{s}*dw*/

_{c}*dt*≃ |

*b*|, except for molecular processes, which would be small for moderate velocities [see Eq. (12)]. Turbulence modifies these values in two ways. First, it reduces the characteristic buoyancy

_{s}*b*to a fraction of |

_{c}*b*| of about 0.7 simply because of the turbulent fluctuations at a given lateral position

_{s}*y*. Second, the lateral turbulent transport of vertical momentum toward the sides of the shell reduces the maximum acceleration even more, to a value of

*c*

_{1}≃ 0.25.

The second parameter *c*_{2} is related to the mean vertical distance *h* traveled by an observer moving with the characteristic velocity *w _{c}*, that is,

*w*=

_{c}*dh*/

*dt*. Indeed, from Eq. (20) we obtain

as the distance traveled by the shell during the interval of time required to increase the width from *δ*_{1} to *δ*_{2} > *δ*_{1}. Then, Eqs. (21) and (23) become simply the well-known expressions corresponding to free-fall motion with a reduced gravity *c*_{1}|*b _{s}*| and a distance

*δ*/

*c*

_{2}.

Equation (24) also shows that the parameter *c*_{2} is the spreading rate (*δ*_{2} − *δ*_{1})/*h*. Values on the order of 0.1 like the ones obtained here in the subsiding shell are common in other shear-driven configurations (Pope 2000). Similarly, *dδ*/*dt* can be interpreted as a mean entrainment velocity in the process of lateral mixing between the shell, and therefore the cloud, and the environment. Then, according to Eq. (20), this second parameter *c*_{2} quantifies this entrainment velocity as a fraction of the mean vertical velocity measured inside the shell.

## 4. Implications under atmospheric conditions

### a. Strength of the subsiding shell

The nondimensional formulation is useful for the derivation of the scaling laws, but dimensional variables corresponding to a relevant reference case become appropriate for further discussion of the results. Using a typical value *b _{s}* = −0.03 m s

^{−2}for the saturation buoyancy, which corresponds to an evaporative cooling of about 1 K, the simulations presented here cover the first 23–31 s of the mixing event, depending on the particular case in Table 1. In that time, a shell thickness between

*δ*≃ 0.36 and 0.54 m is established. The corresponding velocity scales vary between

*w*≃ 0.23 and 0.28 m s

_{c}^{−1}, respectively. The higher values correspond to case

*χ*= 0.3, which has been intentionally studied in a larger domain, on the order of 8 m in the vertical direction, because this case is more representative of typical conditions inside the shell, as argued in section 2.

_{s}Rewriting Eq. (23) as

a characteristic time *t* − *t*_{1} necessary for the shell to broaden a multiple *δ*/*δ*_{1} of the initial thickness *δ*_{1} can be obtained explicitly in terms of the thermodynamic parameters of the cloud-environment system *b _{s}* and

*χ*. For instance, for a strength of evaporative cooling in the range |

_{s}*b*| = 0.02–0.04 m s

_{s}^{−2}, we find that growing from

*δ*= 1 to 10 m requires about 100–140 s, reaching a velocity in the range

*w*= 1.0–1.4 m s

_{c}^{−1}. Those estimates agree well with the typical time scales of a shallow cumulus cloud, where an in-cloud updraft has a lifetime on the order of 400 s (Heus et al. 2009a). In other words, a subsiding shell that is based on the processes as described in the current study will grow fast enough to reach the observed size and strength, but further growth will be inhibited by the larger-scale motions of the cloud. The mean vertical distance traveled during this 120 s is about 90 m, according to Eq. (24). This is again consistent with previous work, where it is observed that parcels are transported down over 200 m because of the subsiding shell (Heus et al. 2008). Thus, the current framework seems to be capturing the turbulent mixing in the subsiding shell reasonably well.

There is one number where our results do not match with the results of Heus and Jonker (2008); that study reported a vertical velocity minimum around 0.1 m s^{−1}, apparently much smaller than our estimate of over 1 m s^{−1}. There are a few reasons for the discrepancy, beyond the fact that our results are likely a drastic idealization of the atmospheric situation. First, the sampling technique used in Heus and Jonker (2008) tends to smooth the shell out, as it does not align the data at the edge of the cloud. This is evidenced by Jonker et al. (2008), who focused their sampling at the cloud edge, obtaining values of around 0.5 m s^{−1}. Second, as suggested in Heus et al. (2009b), and also by our shell size *δ* < 25 m, the minimum velocity is not yet fully resolved in their simulations, although they claim that the mass flux through the shell is. It is important to realize that we are talking about the velocity minimum of a fully developed shell, not the velocity averaged over the shell width and over all life stages of the shell. Finally, the study of Heus and Jonker (2008) investigated clouds and their subsiding shell in a sheared environment. While their sensitivity studies showed limited change in the vertical velocity minimum under changes in shear, two effects where notable: 1) evaporative cooling tends to happen predominantly at the downshear side of the cloud, and 2) in the case of a strong shear, the shell of slanted clouds will not align, thus smoothing out the vertical velocity profile. Obviously, these effects have not been incorporated in our simulations, this being a study on the evolution of the flow given a localized event of evaporative cooling. The impact of shear, together with the irregular interface of the cloud, on the evolution of the shell is a topic for future studies.

### b. Turbulent entrainment and the cloud boundary

Last, we return to the discussion of entrainment. We have already used Fig. 3 to introduce the intermittency—the alternation between turbulent and nonturbulent regions—and the corresponding turbulent–nonturbulent interface. We have also indicated that the cloud boundary is different from the turbulent–nonturbulent interface and the former is typically inside the turbulent region. Hence, it seems appropriate to distinguish between turbulent entrainment and cloud entrainment. This distinction is further quantified in Fig. 11 and now discussed.

Turbulent entrainment is defined as the process by which the fluid particle acquires vorticity fluctuations as it moves from the nonturbulent (irrotational) region to the turbulent (rotational) region (Pope 2000). The topic has attracted considerable attention in the last decade, in terms of the local description of the turbulent–nonturbulent interface (Westerweel et al. 2005; Holzner and Lüthi 2011), bulk zonal statistics (Mellado et al. 2009b), or subgrid-scale modeling (da Silva 2009). One major statistic is the intermittency factor *γ*(*y*, *t*), which measures the relative frequency of occurrence of a turbulent signal at a given lateral position *y* and time *t*. It is simply defined by the vorticity rms being larger that a prescribed threshold, and the curve in Fig. 11 corresponds to a threshold 10^{3} times smaller than the maximum value inside the domain. The values *γ* = 0 and 1 indicate fully nonturbulent and fully turbulent regions, respectively. The sensitivity of *γ* to variations in the threshold between 10^{2} and 10^{4} is relatively small for the resolutions employed in these simulations. The mean position of the turbulence interface is therefore at a distance 0.5*δ*–1.0*δ* from the center of the shell. It is also observed in Fig. 11 that the profile of the viscous dissipation rate follows relatively closely the average position of the turbulence boundary and hence it can be used as a good indicator of that mean position. Occasionally, the scalar is also used as a surrogate to study some bulk properties of turbulent entrainment. In our case, an intermittency factor based on a threshold of the magnitude of the scalar gradient (Fig. 3c) is indistinguishable from that in Fig. 11 based on the enstrophy.

The average position of the cloud boundary is provided in Fig. 11 in terms of the cloud fraction. That figure confirms that the cloud boundary lies well inside the turbulent region for the case of the subsiding shell considered here. This is of interest because small-scale motions behave differently whether being at the turbulent–nonturbulent interface or inside the core of the turbulent region. It is also inferred from that figure that the dependence of this mean position on the saturation mixing fraction is relatively small: about 30% of *δ* when *χ _{s}* varies by a factor of 5. The figure also shows that the cloud boundary extends over an interval of order

*δ*in the center of the shell and the flow is strongly inhomogeneous there, as inferred from Figs. 7 and 9. Hence, the subsiding-shell framework presented in this paper might be an interesting configuration to extend previous work on cloud microphysics and its interaction with small-scale turbulent motions based on homogeneous turbulence (Lanotte et al. 2009; Andrejczuk et al. 2009; Kumar et al. 2013) and complement observations (Siebert et al. 2006; Lehmann et al. 2009) and laboratory studies (Malinowski et al. 2008; Korczyk et al. 2012). As a first contribution, we have characterized the system in terms of integral scales and bulk properties. However, further discussion of this issue, even within the context of the simplified two-fluid formulation employed here, is beyond the scope of this paper.

## 5. Conclusions

We have presented a framework to study the turbulent mixing locally at the lateral boundary of shallow cumulus clouds. We have used a bulk formulation to describe the cloud physics together with direct numerical simulations to investigate the turbulent flow caused by evaporative cooling. In this formulation, the state of air is set by the mixture fraction *χ*, which represents both the normalized total water content and the normalized enthalpy. An idealization of the shell is made by assuming both in-cloud and environmental air to be equally buoyant.

The system develops a subsiding shell that constantly widens and accelerates because of the evaporative cooling. A self-similar regime is established after an initial transient. The minimum mean buoyancy becomes a constant fraction of the saturation buoyancy *b _{s}*, varying between 0.65

*b*and 0.70

_{s}*b*as the saturation mixture fraction

_{s}*χ*varies between 0.5 and 0.1. The shell thickness

_{s}*δ*grows quadratically in time according to

*δ*≃ (

*c*

_{1}

*c*

_{2}

*b*

_{2}/2)

*t*

^{2}and the magnitude of the minimum mean velocity

*w*increases linearly in time according to . The coefficients

_{c}*c*

_{1}(

*χ*) and

_{s}*c*

_{2}(

*χ*) are approximately independent of

_{s}*χ*, equal to 0.25 and 0.1, respectively.

_{s}When put in an atmospheric perspective, these scaling laws imply that the shell around typical shallow cumulus clouds develops within 2 min to a width on the order of 10 m and a velocity of 1 m s^{−1}, which is mostly in agreement with observations. This agreement gives us confidence in the described framework, and in the use of the scaling laws summarized in the previous paragraph to estimate the potential strength of a subsiding shell.

This study shows that evaporative cooling is capable of creating a shell of sufficient size in a sufficiently small time. However, especially in the light of the quadratic growth, the 2-min time scale and the 10-m width are somewhat arbitrary. No mechanism is built in the current model setup to limit the size of the shell. Such a limiting mechanism has to be provided by larger scales. What this mechanism should be depends on the circumstances. For small or short-lived clouds, the assumption of an infinite cloud size will quickly become a limiting factor, or the cloud will simply reach the end of its life time before further development can occur. More generally speaking, the time scales of the larger scale in-cloud turbulence will at some point limit the entrainment event. Finally, stratification of the environment will always limit the vertical extent of the shell, and through that the width and vertical velocity of the shell. In typical shallow cumulus conditions, this would be a vertical extent on the order of 300 m and a maximum width of 30 m. After reaching this level of neutral buoyancy, the shell may propagate a bit farther based on its inertia, to extend away from the cloud and form the gravity wave that was suggested by Verzijlbergh et al. (2009). Because of the inviscid nature of the scaling laws describing the dynamics of the shell, LES could be used to study this interaction between the shell and the cloud.

Other mechanisms beyond our current focus are likely to play some role in the physics of the subsiding shell. We would like to stress that the exact numbers being retrieved should not be taken as a precise prediction of the characteristics of the shell in all its atmospheric complexity. A less idealized study of the subsiding shell could involve, for instance, a cloud with positive vertical velocity, positive buoyancy, and perhaps a finite size. The microphysical structure of the mixing event might be of interest as well; specifically, a finite condensation/evaporation rate may influence the turbulent mixing here. These kinds of effects should all be the subject of future research. However, given the similarities with observations and the parameter studies that we have conducted, it seems likely that the key elements of the subsiding shell have been captured in the current study.

## Acknowledgments

Support from the Max Planck Society through its Max Planck Research Groups program is gratefully acknowledged. Dick Abma was funded by the *ZEIT-Stiftung Ebelin und Gerd Bucerius* and the Max Planck Foundation. Thijs Heus was funded by the Deutscher Wetterdienst (DWD) through the Hans Ertel Centre for Weather Research. Computational time was provided by the Jülich Supercomputer Centre and the German Climate Computing Center.

### APPENDIX

#### Self-Similarity Analysis

We summarize here the intermediate steps of the self-similar analysis used in the previous sections to obtain the different scaling laws. [Details of this approach can be found, for instance, in Tennekes and Lumley (1972) or Pope (2000).]

We look for solutions in the self-similar forms of Eqs. (14) and (15). From these definitions we have

where

Substituting into the mean-momentum transport equation

and dividing by *b _{c}* yields

having defined

We have observed in section 3b that the ratio *y _{w}*/

*δ*approaches a constant after the initial transient. Hence, if the Reynolds number Re

_{c}=

*w*/

_{c}δ*ν*is large enough, the self-similar behavior documented in that same section requires that the groups within curly brackets are constant for Eq. (A4) to depend only on the self-similar variable

*ξ*[or

*ζ*, since both are related through the constant (

*y*−

_{w}*y*)/

_{b}*δ*]. These three ratios are those in Eq. (17).

Also, the integral form of the mean-momentum equation

or, equivalently,

provides a constraint that reduces to two the number of independent self-similar parameters identified in Eq. (A4). For that purpose, we need to eliminate . This can be achieved by using the equation above to write

Since is constant in the self-similar regime, we obtain

## REFERENCES

*Hydrodynamic Stability*. Cambridge University Press, 628 pp.

*Statistical Fluid Mechanics: Mechanics of Turbulence*. Dover Publications, 896 pp.

*Turbulent Flows*. Cambridge University Press, 806 pp.

*A First Course in Turbulence*. MIT Press, 300 pp.

*Atmos. Chem. Phys.,*

**9,**1289–1302.