## Abstract

Direct numerical simulations resolving meter and submeter scales in the cloud-top region of stratocumulus are used to investigate the interactions between a mean vertical wind shear and in-cloud turbulence driven by evaporative and radiative cooling. There are three major results. First, a critical velocity jump exists, above which shear significantly broadens the entrainment interfacial layer (EIL), enhances cloud-top cooling, and increases the mean entrainment velocity; shear effects are negligible when the velocity jump is below . Second, a depletion velocity jump exists, above which shear-enhanced mixing reduces cloud-top radiative cooling, thereby weakening the large convective motions; shear effects remain localized within the EIL when the velocity jump is below . The critical velocity jump and depletion velocity jump are provided as a function of in-cloud and free-tropospheric conditions, and one finds and for typical subtropical conditions. Third, the individual contributions to the mean entrainment velocity from mixing, radiative cooling, and evaporative cooling strongly depend on the choice of the reference height where the entrainment velocity is calculated. This result implies that the individual contributions to the mean entrainment velocity should be estimated at a comparable height while deriving entrainment-rate parameterizations. A strong shear alters substantially the magnitude and the height where these individual contributions reach their maxima, which further demonstrates the importance of shear on the dynamics of stratocumulus clouds.

## 1. Introduction

Wind shear in the cloud-top region can significantly alter the temporal evolution of the stratocumulus-topped boundary layer (STBL), as shown by a number of observational studies (Brost et al. 1982; Caughey et al. 1982; Driedonks and Duynkerke 1989; Faloona et al. 2005; de Roode and Wang 2007; Katzwinkel et al. 2012; Malinowski et al. 2013; Jen-La Plante et al. 2016) and some numerical experiments (Wang et al. 2008, 2012; Kopec et al. 2016). However, the aspect of meter- and submeter-scale mixing processes has obtained less attention, even though former studies have shown that these small-scale processes are crucial for the dynamics of the cloud in general and for shear effects in particular (Katzwinkel et al. 2012; Malinowski et al. 2013; Mellado 2017). Here, direct numerical simulations (DNS) are employed to explicitly resolve these small-scale processes. For a cloud top solely driven by evaporative cooling and shear, shear effects are found to enhance mixing mainly within a shallow layer (Mellado et al. 2014). The thickness of this layer is typically a few tens of meters or less, confirming the importance of small-scale processes. However, radiative cooling has been neglected in the former study, which motivates us to investigate how a vertical wind shear alters the dynamics of a radiatively and evaporatively driven stratocumulus cloud top.

The first goal is to identify when shear effects become relevant. Shear can enhance the entrainment of tropospheric air (cf. studies cited above), can thicken the entrainment interfacial layer (EIL; Wang et al. 2008; Katzwinkel et al. 2012; Jen-La Plante et al. 2016), and can change the budget of turbulent kinetic energy (TKE; Caughey et al. 1982; Kopec et al. 2016; Jen-La Plante et al. 2016). Nonetheless, it remains unclear at which minimal shear strength significant changes in these quantities occur. We answer this question by deriving a critical velocity jump, below which shear effects are negligible, and a depletion velocity jump, below which shear effects remain localized within the cloud-top region and in-cloud turbulence remains unaffected.

The second goal is to quantify shear effects on the mean entrainment velocity. The magnitude of radiative and evaporative cooling drastically varies with height within a few meters, which renders shear broadening of the EIL, despite being small (~10 m), a crucial process for understanding shear effects. Especially with respect to the mean entrainment velocity *w*_{e}—here defined as the time rate of change of a reference height marking the inversion atop the cloud (Lilly 1968)—resolving these small-scale processes is critical, as different definitions of the reference height differ only by a few meters. Previous measurements and numerical studies (e.g., Stevens et al. 2003; Faloona et al. 2005; Gerber et al. 2016) indicate that these small height differences might be crucial for entrainment velocity parameterizations. This motivates us to investigate how *w*_{e} depends on the choice of the reference height.

A related observation is that most local analyses of cloud-top entrainment assume a quasi-steady state (i.e., a state in which the in-cloud and free-tropospheric conditions change slowly, compared to the cloud-top processes). Nonetheless, it is known that this is not always the case (e.g., during transients), and, at least for a dry atmospheric boundary layer, unsteady effects are reported to affect the entrainment velocity substantially (Sullivan et al. 1998). For such an unsteady state, the shape of the mean profiles changes significantly in time; therefore, different reference heights can evolve differently in time, which implies that the magnitude of *w*_{e} depends on the choice of the reference height. Still, an explicit quantification of unsteady effects is missing to the best of our knowledge. Here, we provide a first attempt to quantify unsteady effects by analyzing the corresponding term in the entrainment-rate equation.

The paper is structured as follows. Section 2 introduces the cloud-top mixing layer (CTML), defines the simulation setup, and reviews some of the fundamental concepts and quantities needed for the analysis. Section 3 investigates how shear affects the vertical structure of the cloud top, while section 4 investigates when shear effects start to become significant. Section 5 discusses shear effects on the entrainment velocity. Results are summarized and discussed in section 6.

## 2. The cloud-top mixing layer

The CTML mimics the upper part of the STBL and consists of a region of warm and dry air, representing the free troposphere, and a region of moist, relatively cold air, representing the cloud below (cf. Fig. 1). The formulation of the CTML is identical to the one used by de Lozar and Mellado (2015), where a CTML solely driven by evaporative and radiative cooling has been investigated by means of DNS. Here, we extend this work by imposing a vertical wind shear. For conciseness, the detailed formulation is presented in appendix A, and this section only includes the description of the parameters and variables needed for the discussion of the results.

### a. Description of the simulations

#### 1) Simulation parameters

Once the system has become sufficiently independent of the initial conditions, flow properties only depend on the height *z*, the convective length scale , characterizing the large-scale turbulent motions in the cloud [cf. section 2b(1)], and six nondimensional parameters . The reference Reynolds number, , and the reference Richardson number, , are based on two radiative reference scales, namely, the extinction length and the reference buoyancy flux . is the reference longwave radiative cooling at the cloud top [cf. Eq. (A6)], is the density of cloudy air, is the specific heat capacity of cloudy air, is the temperature of cloudy air, *ν* is the kinematic viscosity, and the subscript 0 indicates reference values. Based on the former two parameters, we can define a reference velocity and a reference buoyancy scale as

respectively. The buoyancy reversal parameter compares the buoyancy at saturation conditions with the buoyancy jump across the inversion , where the superscript *d* indicates dry conditions and the superscript *c* indicates cloudy conditions. Buoyancy reversal instabilities are associated with (Randall 1980; Deardorff 1980a). The parameter indicates the mixture fraction at saturation conditions, while *β* describes how enthalpy changes translate into buoyancy changes. The last two parameters are explained in more detail in appendix A.

To characterize shear effects, we introduce a shear number *S* as

where specifies a constant vertical velocity jump across the cloud top. The vectors and represent the mean velocity in the dry free troposphere and in the cloud, respectively. Since we can always choose a reference frame that moves with the mean velocity and that is aligned with , the parameter *S* is sufficient to characterize shear effects in the CTML.

#### 2) Simulation setup

In this work, we fix all parameters according to the first research flight (RF01) of the DYCOMS II campaign (Stevens et al. 2003, 2005; see Table 1), and we vary the shear number by varying the initial velocity jump (see Table 2). Since we consider RF01 of DYCOMS II as a reference, our simulation setup resembles subtropical clouds, which are characterized by relatively large jumps in total-water specific humidity and temperature across the inversion. We match all parameters of the RF01 of the DYCOMS II campaign, except the Reynolds number; therefore, we need to study the dependence of our results on the Reynolds number. This dependence is discussed in appendix B. For the Reynolds numbers reached in our simulations, the properties relevant for the discussion in this paper show only a weak dependence on the Reynolds number. This tendency toward Reynolds number similarity, which is a general characteristic of turbulent flows (Dimotakis 2005; Mellado et al. 2018), partly justifies the extrapolation of our results to atmospheric conditions.

The grid spacing is uniform and isotropic in the region of the computational domain where the turbulent flow develops. The ratio between the grid spacing and the Kolmogorov length *η* is approximately 1.5, which is sufficient for the statistical properties of interest to depend less than 5% on the grid spacing, which is comparable to or less than the statistical uncertainty of the properties considered in this work (Mellado 2010; Mellado et al. 2014). For the conditions of RF01 of DYCOMS II, the corresponding grid spacings vary between 16 and 32 cm, depending on (see Table 2). For the compact schemes used in these simulations, about four points per wavelength provide 99% accuracy in the transfer function of the derivative operator, which implies that we reach submeter-scale resolution in these studies. [For comparison, second-order central schemes need about eight points per wavelength to reach 90% accuracy, which is the motivation to employ compact schemes despite being computationally more demanding (Lele 1992).] The size of the computational domain in the horizontal direction is , except for the low Reynolds number cases with and , where we doubled the domain size to improve statistical convergence. In the vertical direction, we stretch the grid spacing to separate the boundaries of the computational domain while reducing the computational costs (Mellado 2010; Mellado et al. 2014). The resulting vertical domain size is approximately 600 m for the cases and approximately 300–400 m for all other cases. Further simulation details are given in de Lozar and Mellado (2015), and details about the numerical algorithm can be found in Mellado and Ansorge (2012).

### b. Description of the vertical structure

#### 1) In-cloud convective scalings

The prevalence of free convection in the cloud suggests introducing a convective length scale to characterize the depth of the convective region and the size of the large-scale motions in the cloud. According to Deardorff (1980b) and following Mellado et al. (2014), we define as

where denotes the Heaviside function, is the turbulent buoyancy flux, is the maximum of *B* within the cloud, angle brackets indicate a horizontal average, and an apostrophe indicates the turbulent fluctuation field. The Heaviside function ensures that only the positive part of the turbulent buoyancy flux profile, which generates turbulence, is retained. We will show later (Fig. 9) that the height of zero buoyancy flux is located near the height of zero mean buoyancy , and thus, is mainly associated with a negatively buoyant region.

A measure of the intensity of the in-cloud turbulence is provided by the convective velocity scale (Deardorff 1970):

and, indeed, the maximum of the TKE within the cloud follows the scaling law for (not shown). Note that our definition of is smaller by a factor of , compared to previous work that considers the whole STBL (Deardorff 1980b; Wood 2012). The reason is that the buoyancy flux in the CTML setup does not have the linear vertical variation characteristic of the subcloud layer in the STBL, which justifies the factor of 2.5.

Because of a continuous cloud-top cooling, the turbulent buoyancy flux increases with time, and hence, increases with time. We can use this relationship between time and to express the evolution of the system in terms of the nondimensional variable . Introducing this variable has the advantage that represents the scale separation between the integral scale of the in-cloud turbulence and the scale at which the radiative forcing is introduced. In addition, represents the intensity of the in-cloud turbulence, according to Eq. (4). We reach in our simulations, which is within the range of *z*_{*}/*λ* ≃ 3–80 reported in Deardorff (1981). To improve statistical convergence, the results are averaged in time over a period of , which corresponds to approximately five to seven eddy turnover times .

#### 2) The entrainment interfacial layer

In general terms, the EIL refers to the layer where the entrainment of dry and warm tropospheric air takes place and thus represents a transition layer between the cloud and the free troposphere. Therefore, the EIL is characterized by strong vertical variations of temperature, specific humidity, mean vertical velocity, and buoyancy. We hence define the EIL thickness as

where is the height of zero mean buoyancy, and denotes the height where the mean buoyancy has increased by 90% of . According to this definition, the EIL is stably stratified and contains the cloud boundary and the turbulent–nonturbulent interface (see Fig. 2).

For weak-enough shear, the EIL dynamics is determined by the in-cloud turbulent motions penetrating the free troposphere. This process is characterized by the penetration depth, here defined as

where denotes the height of minimum mean buoyancy flux. Previous work on shear-free conditions has shown that the thickness associated with the evaporative cooling caused by diffusion is comparable to for low to moderate Reynolds numbers (Mellado 2010; de Lozar and Mellado 2015). For the sheared configurations considered in this study, this is also the case, and we find *h*_{diff}/*h*_{EIL} = 0.5–1.0, which indicates that needs to be retained in the scaling law of . Figure 3 demonstrates that is scaled by .

For strong-enough shear, however, locally generated turbulence enhances mixing, which can thicken the EIL substantially. This shear effect can be characterized by the reference length (Mellado et al. 2014):

which is henceforth referred to as the critical-shear-layer thickness, where the subscript *S* indicates “shear.” The factor 1/3 in is well established (±15%) from laboratory and numerical experiments of stably stratified shear layers (Smyth and Moum 2000; Brucker and Sarkar 2007) and from observations and numerical experiments of cloud-free sheared convective boundary layers (Mahrt and Lenschow 1976; Fedorovich and Conzemius 2008), and it can be associated with a critical value of the shear Richardson number . A strong shear can broaden the EIL sufficiently for the EIL thickness to be well approximated by the critical-shear-layer thickness . This is shown in Fig. 3 for the case , where we still need to add the diffusion correction to account for the low-Reynolds-number effect, as we did in *δ*.

The physical interpretation of the critical-shear-layer thickness is rationalized as follows. Given a stably stratified shear layer characterized by a buoyancy jump and a velocity jump , if the initial shear layer is thin enough, Kelvin–Helmholtz instabilities will cause an overturning of the stably stratified fluid and a thickening of the shear layer. As the shear layer thickens, overturning the fluid becomes more difficult because the vertical displacement increases, whereas the available kinetic energy, proportional to , remains constant. Once the shear layer has grown to its critical thickness , the available kinetic energy is insufficient to overturn the fluid and turbulence decays.

We emphasize that is an average quantity, and there is a range of smaller motions in the EIL that locally can have a gradient Richardson number below 1/3 (Kurowski et al. 2009; Malinowski et al. 2013). In particular, in-cloud turbulent motions penetrate into the stably stratified EIL, which locally thins the inversion and creates further shear instabilities that increase the shear-layer thickness toward its critical value . In the absence of this thinning process, shear-generated turbulence would decay once the EIL reaches the critical thickness . This dynamical equilibrium between penetrating thermals and shear instabilities permanently maintains mixing within the cloud-top region, despite the mean Richardson number being near the critical value for stability (Mahrt and Lenschow 1976; Fedorovich and Conzemius 2008; Howland et al. 2018).

In this study, we consider shear-free conditions as a reference case and increase the cloud-top velocity jump to study weak shear conditions with and and strong shear conditions with and . Note that even for strong shear conditions, the penetration depth of in-cloud turbulence *δ* is still comparable to (see Table 2, Fig. 3). The limit of very strong shear where is much larger than the penetration depth of in-cloud motions, considered in Mellado et al. (2014) for the solely evaporatively driven case, is not investigated here.

Finally, we note that our definition of the EIL [Eq. (5)] follows the definition of Caughey et al. (1982) as the layer containing most of the temperature jump (buoyancy, in our case). This definition approximately coincides with the region of negative buoyancy flux, which is often referred to as the entrainment zone in cloud-free boundary layers (e.g., Fedorovich and Conzemius 2008). This similarity proved convenient to use results from the study of shear effects in cloud-free conditions—in particular, the relevance of the critical shear layer and its thickness . However, this definition of EIL differs from the one proposed by Malinowski et al. (2013), especially in the definition of the lower boundary of the EIL, which has to be taken into account when comparing results. In Malinowski et al. (2013), the lower boundary of the EIL is defined as the height where the square of the horizontal wind shear reaches 90% of its maximum value. In our definition of the EIL, this height roughly coincides with the height of minimum turbulent buoyancy flux , which is approximately located in the center of the EIL (cf. Fig. 4).

## 3. Vertical extent of wind shear effects

This section studies the vertical extent of wind shear effects. The first part of the analysis is based on the TKE evolution equation, which is given by

where is the turbulent kinetic energy, is the vertical turbulent flux, is the shear-production rate, is the turbulent buoyancy flux, is the viscous dissipation term, and is the fluctuating part of the viscous stress tensor.

Based on the TKE budget in the EIL, we can distinguish between two limiting regimes (cf. Fig. 4a). For , the shear production term is zero, and the turbulent transport term −∂_{z}*T* is the only source of TKE in the EIL, where part of the TKE is dissipated and part is used to entrain warm and dry air from above, as indicated by the negative turbulent buoyancy flux. For , the shear production term is the dominating source of TKE in the EIL, where part of the TKE is dissipated, part is redistributed by the transport term, and part is used to enhance the entrainment of tropospheric air. These observations are in general agreement with previous work (e.g., Kopec et al. 2016; Jen-La Plante et al. 2016). In addition, it is noteworthy that for , the input of TKE is strong enough to generate an area within the EIL where the transport term turns negative, indicating an export of TKE. By increasing the shear strength, we hence change the ratio between the turbulent transport term and the shear-production term, and the system transitions from a transport-dominated regime within the EIL to a shear-dominated regime.

In contrast, within the cloud region below the EIL, the TKE and the terms in its evolution equation remain approximately the same for different shear numbers (cf. Fig. 4b). This independence of in-cloud properties to the strength of cloud-top wind shear is further demonstrated in Fig. 5, where the temporal evolution of the maximum TKE within the cloud is approximately independent of *S*, whereas the maximum within the EIL increases by about 60% when increasing the shear strength from to . Hence, wind shear effects remain localized within the EIL.

Although shear effects remain localized within the EIL for all investigated cloud-top velocity jumps, shear can affect large-scale properties if is large enough for the critical shear layer thickness to become comparable to the cloud thickness *H*. In this condition, shear-enhanced mixing depletes the cloud sufficiently to reduce the net radiative cooling, thus weakening turbulence and favoring a decoupling of the STBL (Wang et al. 2008, 2012; Kopec et al. 2016). In what follows, we estimate the depletion velocity jump at which this shear effect occurs.

The radiative flux difference between cloud top and cloud base is , where

This expression follows from Eq. (A6) and the assumption of an adiabatic liquid lapse rate in the cloud. If the radiative extinction length is much smaller than the cloud thickness *H*, then , and the radiative flux difference between the cloud interior and the free troposphere is . How much can shear reduce the liquid-water specific humidity in the cloud-top region for this radiative flux difference to remain approximately ?

Let us consider that shear-enhanced mixing causes the evaporation of a cloud layer of thickness *d* at the top, reducing the cloud thickness from the initial value *H* to (cf. Fig. 6). This process reduces the maximum liquid-water specific humidity at the cloud top from the initial value to , if an adiabatic liquid lapse rate is assumed. This reduction increases the radiative extinction length from to because the radiative extinction length is inversely proportional to the maximum liquid-water specific humidity (Larson et al. 2007; Wood 2012; Mellado 2017). For this new, partially depleted cloud, the radiative flux difference between cloud top and cloud base is , where

Since , the depleted cloud has a smaller radiative flux difference between the cloud interior and the free troposphere. However, as long as , the radiative difference remains approximately , and the radiative forcing of the boundary layer turbulence remains approximately the same.

To relate *d* to the velocity jump , we assume that *d* is proportional to the thickness of the critical shear layer . As a first approximation, we assume that the critical shear layer is centered at the height of minimum turbulent buoyancy flux *z*_{i,f} and mixes and evaporates a cloud layer of a thickness . Substituting this relation into Eq. (10) and solving for yields

as an estimate for the critical-shear-layer thickness that is necessary to obtain a reduction of the radiative flux difference between the cloud top and the cloud base. The corresponding depletion velocity jump follows from Eq. (7) and is given by

Therefore, only a strong wind shear with can deplete the cloud sufficiently to change the net radiative cooling. Equation (12) shows that increases with increasing inversion strength (i.e., when the buoyancy jump increases). This dependence seems reasonable since increasing the inversion strength hinders entrainment. Equation (12) further shows that increases with the cloud thickness. This result also seems reasonable; increasing the cloud thickness implies that cloud-top depletion needs to extend over a thicker layer to change the net radiative cooling of the cloud. Last, we note that the proportionality used in the derivation of Eq. (12) may depend on the thermodynamic conditions. It is, for example, expected that moistening the free troposphere increases by weakening evaporation. Hence, further assessment of the dependence of on the thermodynamic conditions is necessary. However, as discussed in the following paragraph, Eq. (12) provides a leading-order estimate for the effect of shear broadening on the radiative cooling forcing that is consistent with observations.

For the case RF01 of the DYCOMS II field campaign, we find , where we used , , , and an arbitrary threshold of (i.e., a 5% reduction of the net radiative flux difference across the cloud-top region, compared to the no-shear case). This result is consistent with field measurements reporting a compact cloud layer and a strong radiative forcing in that case, since the velocity jump across the cloud top is only . We can generalize this result and consider an interval of cloud thickness between 100 and 200 m, which yields an interval of depletion velocity jump (Δ*u*)_{dep} ≃ 3–10 m s^{−1}. This range seems consistent with measurements from different field campaigns, which report velocity jumps up to 4–10 m s^{−1} for compact clouds, but not much higher (Brost et al. 1982; Nicholls and Leighton 1986; Faloona et al. 2005; de Roode and Wang 2007; Katzwinkel et al. 2012; Malinowski et al. 2013). However, the range of is partly in disagreement with previous numerical experiments (Wang et al. 2008, 2012; Kopec et al. 2016), where cloud-top depletion and decoupling of initially well-mixed STBLs is observed for smaller velocity jumps. This disagreement might be partly caused by the excessive mixing associated with the subgrid models used in those large-eddy simulations. Other numerical artifacts, like numerical diffusion and effects of the gridbox aspect ratio, could further contribute to this disagreement (Stevens et al. 1999; Pedersen et al. 2016, 2018).

## 4. Strong and weak shear regimes

The previous section (e.g., Fig. 4) indicates that shear effects become relevant for shear numbers larger than *S* ≃ 5–10, which indicates a transition between the weak shear regime and the strong shear regime described in section 2b(2). This section rationalizes this behavior in terms of two length scales and provides an analytic expression for the critical velocity jump , beyond which shear effects in the EIL are significant.

### a. The penetration depth

Considerations of kinetic and potential energy within the EIL allow us to derive a scaling law for the penetration depth *δ*, defined by Eq. (6), as follows. To a first approximation, air parcels with a kinetic energy can penetrate a distance *δ* into the EIL until all their kinetic energy is converted into a potential energy . As a first approximation, the kinetic energy of an air parcel within the EIL is estimated by the sum of a kinetic energy associated with in-cloud convective motions—characterized by —and a kinetic energy associated with the shear production in the EIL—characterized by . Hence, we write

where the ratio represents the gradient of the horizontal mean velocity within the EIL. Likewise, the potential energy at a height *δ* within the EIL is given by

where the ratio characterizes the gradient of the mean buoyancy profile within the EIL. The height that parcels can penetrate into the EIL is hence implicitly given by the energy balance , which allows us to write

A similar equation is analyzed in Haman (2009) for shear-free conditions. Figure 7 supports this linear relationship, and a linear regression to the data provides the parameters and .

In summary, Eq. (15) yields the following expression for the penetration depth *δ*:

Two simplified expressions of Eq. (16) are obtained by introducing different scalings for . For a weak-enough shear, in-cloud turbulence dominates mixing in the EIL, and we can estimate for large-enough Reynolds numbers, as discussed in section 2b(2). In this case, Eq. (16) simplifies to

where the subscript *C* refers to “convection.”

For a strong-enough shear, we can alternatively apply the scaling , as explained in section 2b(2). In this case, Eq. (16) simplifies to

where the subscript *S* refers to “shear.” This scaling is consistent with the definition of the penetration depth by Eq. (A1) in Mellado et al. (2014), which is obtained in the limit of a very strong shear (i.e., ). We find for , which is reasonable since is equally well scaled by and for [cf. section 2b(2)]. We hence conclude that convection and shear are comparably important for .

### b. The critical velocity jump

Shear effects are negligible for , and in-cloud turbulent convection penetrating into the EIL dominates the EIL dynamics (cf. section 3). In this regime, the EIL thickness is well characterized by the penetration depth *δ*, namely, when we extrapolate the result in Fig. 3 to high Reynolds numbers. As shear intensifies, shear and convection become similarly important for , which corresponds to the condition , according to Fig. 3. Therefore, it seems reasonable to define a critical penetration depth by the condition

with . Substituting by Eq. (7) and by Eq. (17) into Eq. (19) yields the following expression for the critical velocity jump:

Shear effects become relevant for shear velocities larger than this critical value, which allows us to distinguish a continuous transition across three regimes as we increase the velocity jump. For , convection dominates and shear effects are negligible; for , convection and shear are equally important; and for , shear dominates. According to Figs. 3 and 7, respectively, we set , , and , which finally results in

Typical values of are given in Fig. 8, where the range of *w*_{*} ≃ 0.2–0.9 m s^{−1} is chosen according to Wood (2012). Figure 8 shows that only the case exceeds the critical velocity jump, which is consistent with the observation that the shear Richardson number Ri_{S} is significantly smaller than unity for this case (cf. Table 2).

It is remarkable that, according to Eq. (21), is independent of the inversion strength as measured by the buoyancy jump . The reason is that and *δ* are both inversely proportional to , which shows that effects of cancel each other to leading order. Still, the convective velocity scale can depend on the thermodynamic parameters , where the dependence on *D* introduces implicitly a dependence on ; these potential effects need to be further investigated. The dependence of Eq. (21) on the Reynolds number also needs to be further assessed. However, the coefficients and change by less than 30% and 15%, respectively, when increasing the Reynolds number up to a factor of 3. This implies that changes by 10% or less, which provides certain support to extrapolate our results to atmospheric conditions.

## 5. Wind shear effects on the entrainment velocity

Following Lilly (1968), we define the mean entrainment velocity as

where is a reference height marking the cloud-top region, and is a mean vertical velocity at . Mixed layer models of the STBL use the same definition of mean entrainment velocity (e.g., Stevens 2002; de Roode et al. 2014), and hence our results can be easily interpreted in such a mixed layer framework. Henceforth, a subscript indicates that the corresponding quantity is evaluated at . Common choices of reference heights are the height of maximum gradient of the mean buoyancy , which sits on top of the height of minimum turbulent buoyancy flux , which sits on top of the height of zero mean buoyancy , which sits on top of the height of zero turbulent buoyancy flux . Figure 9 shows that for RF01 of DYCOMS II, those different heights are separated by only a few meters. Alternative reference heights have also been proposed (e.g., by Malinowski et al. 2013). Hence, in this section, we study how the entrainment velocity depends on an arbitrary definition of a reference height and how wind shear affects this dependence.

To analyze the dependence of on wind shear and on the choice of the reference height , we use the entrainment-rate equation (Mellado 2017; Mellado et al. 2018). The starting point for deriving the entrainment-rate equation is the buoyancy evolution equation:

where and denote the radiative and evaporative source terms, respectively (cf. appendix A). Integrating Eq. (23) from an arbitrary height upward yields the entrainment-rate equation

where , and

The radiative cooling contribution is set by the difference between the net radiative flux above the cloud , and the net radiative flux evaluated at , . Equivalently, the evaporative cooling contribution is set by the difference between the integrated evaporative cooling of the cloud and the integrated evaporative cooling at , , which is given by

Finally, quantifies deformation effects of the mean buoyancy profile (i.e., how temporal changes in the shape of the mean buoyancy profile affect ).

Figure 10 shows that the radiative and evaporative cooling rates strongly vary with height. Radiative cooling peaks at ≃5 K h^{−1} slightly below —outside the EIL—and decays within a few tens of meters inside the cloud. This behavior is consistent with the vertical profiles of radiative cooling reported in previous work (Stevens et al. 2003; Yamaguchi and Randall 2012; Gerber et al. 2014). In contrast, evaporative cooling peaks near —within the EIL—and decays more rapidly. This behavior qualitatively agrees with previous work (e.g., Yamaguchi and Randall 2012; Wood 2012; Gerber et al. 2016). However, the assumption of phase equilibrium (saturation adjustment) often used in numerical simulations might overemphasize this narrow shape of the evaporative cooling profiles, but whether this effect is significant remains to be investigated. Either way, a quantitative comparison with previous work is difficult because, to the best of our knowledge, previous LES studies and measurement campaigns have not provided the profile of evaporative cooling for DYCOMS II–like conditions. Nonetheless, we can gain confidence in our results by noting that the total evaporative cooling is related analytically to the mean entrainment velocity according to (de Lozar and Mellado 2015; Mellado 2017). This relationship implies that evaporative cooling rates are directly proportional to entrainment rates (for a given set of thermodynamic parameters), and correctly simulating the mean entrainment velocity therefore ensures realistic values of total evaporative cooling . For the investigated case RF01 of the DYCOMS II campaign, we obtain a quasi-steady entrainment velocity of for the no-shear case, which is commensurate with the range of *w*_{e} ≃ 3.9–4.7 mm s^{−1} reported in Stevens et al. (2003) and Faloona et al. (2005).

As a consequence of the rapid variation with height of cloud-top properties, the different contributions to vary strongly with height in the EIL, and hence, Eq. (24) depends strongly on the choice of the reference height . The contributions and represent the cumulative effect of evaporative and radiative cooling and increase monotonically from zero, when is far above the cloud, to the maximum values and , when is deep within the cloud. In contrast, the mixing contribution is evaluated locally at *z*_{i} and has a maximum value near . In the remainder of this section, we investigate how shear affects these vertical dependences, focusing on a region of about 10 m around the EIL. To this aim, we compare the different contributions to the mean entrainment velocity for and in Fig. 11; those two cases are chosen since shear effects on the mean entrainment velocity are observed to be negligible for . The contour plots in Fig. 11 show the variation with height in the vertical axis and with the integral scale of the in-cloud turbulence in the horizontal axis. Entrainment velocities are normalized by the reference value , where we have used the corrected buoyancy flux to account for the fact that due to condensational warming, only a fraction *β* of the enthalpy changes induced by radiative cooling translates into buoyancy changes (cf. appendix A).

### a. Radiative cooling effects

The direct contribution of radiative cooling is determined by the cumulative radiative cooling above (or, equivalently, the flux difference between and the free troposphere) and thus by the amount of liquid water above . For a fixed stratification and a fixed velocity jump , cloudy air penetrates deeper into the EIL as in-cloud turbulence intensifies, as discussed in section 4, which explains the increase of with increasing observed in Fig. 11a. In addition, shear broadening brings more cloudy air into the EIL, implying that shear enhances . For a strong shear with , this enhancement can be large in relative terms, as inferred from Fig. 11b; we obtain, for example, an enhancement of 100%–140% with respect to at the height of zero mean buoyancy . The enhancement of radiative cooling in the EIL is accompanied by a weakening of radiative cooling within the cloud (cf. Fig. 10), since the cloud liquid-water specific humidity and thus the net radiative flux difference across the cloud-top region are held constant in our experiments. In-cloud turbulence intensity, however, stays approximately constant, since the shear-induced weakening of radiative cooling is compensated by a shear-induced enhancement of evaporative cooling (cf. section 5c).

The magnitude of and the effect of shear on it depend on the choice of the reference height . The radiative contribution contributes 1%–3% to the sum (depending on *S* and ) at the height of minimum buoyancy flux , whereas this contribution increases to 8%–15% for the height of zero mean buoyancy . This increase is reasonable since more liquid water accumulates above a reference height as the reference height moves downward and (cf. Fig. 9). Nonetheless, since most of the liquid water and thus most of the radiative flux difference across the cloud-top region is located below , the contribution of to the mean entrainment velocity is small when is evaluated above . Hence, at least for DYCOMS II–like conditions, the parameterization of shear effects on is not a priority, as long as is less than the depletion velocity jump necessary to thin the cloud enough to change the net radiative flux difference across the cloud-top region (cf. section 3).

### b. The mixing contribution

The mixing contribution to consists of a turbulent part and a molecular part. As observed in Figs. 11c and 11d, the turbulent part of the mixing contribution and the effect of shear on it depends significantly on the choice of the reference height, since the turbulent buoyancy flux varies substantially in the EIL (cf. Fig. 4). By definition, the magnitude of reaches a maximum near , and contributes 40%–60% to the sum at this height. For reference heights below , the magnitude of decreases, and so does its relative contribution to the sum . In contrast, the relative contribution of increases above (e.g., 60%–80% for ) since and decay more rapidly than in this region. Regarding shear effects, we observe that a strong shear with significantly increases the magnitude of the negative buoyancy flux and therefore also ; for example, shear increases by 35%–45% with respect to at and by approximately 200% at . All these observations indicate that tends to dominate the sum for reference heights located around and above it, and finding a shear-dependent parameterization of is key to understanding shear effects on the evolution of stratocumulus clouds.

The molecular part of the mixing contribution can become larger than the turbulent part at the top of the EIL, since the turbulent fluctuations decay faster than the mean buoyancy gradient in that stably stratified region. In numerical simulations, however, the molecular contribution is artificially exaggerated, and we need to understand this contribution to interpret the numerical results, even if is irrelevant under atmospheric conditions. For the Reynolds numbers achieved in our simulations, the molecular part is comparable to the turbulent part at , where the turbulent part is maximum, for (cf. Fig. 11e). However, a strong shear substantially decreases the magnitude of in most of the EIL (cf. Fig. 11f), whereas the magnitude of is substantially increased. This result indicates that despite the moderate Reynolds numbers achieved in our simulations, the strong effect of shear on , and thereby on , is appropriately represented.

### c. Evaporative cooling effects

The magnitude of depends strongly on the choice of the reference height , since evaporative cooling varies strongly within a few meters in the cloud-top region (cf. Figs. 10, 11g). The evaporative contribution reaches its maximum for reference heights located near , as more mixing of environmental and cloudy air accumulates above such reference heights; contributes, for example, 80%–90% to the sum at , while contributing only 20%–40% at . This implies that is the dominant contribution to the entrainment velocity for reference heights located near , in contrast to the mixing contribution, which dominates for reference heights located above .

The importance of evaporative cooling is confirmed by Fig. 12, which shows that the total evaporative cooling is comparable to the reference buoyancy flux associated with radiative cooling. As already mentioned, we note that —as well as —are cumulative effects of evaporative cooling, as it appears in Eq. (24), and not the local evaporative cooling rate at a given reference height . Locally inside of the cloud, radiative cooling can be larger than evaporative cooling, as shown in Fig. 10. The reason is that evaporative cooling is concentrated within the EIL; for example, for the subtropical conditions considered in this study, the evaporative cooling above contributes 60% of . Anyhow, total evaporative cooling is expected to decrease if the free troposphere is moistened and if the inversion’s strength is weakened, as, for example, observed during the Physics of Stratocumulus Top (POST) campaign (Malinowski et al. 2013). The relative importance of evaporative cooling and radiative cooling is therefore expected to depend on the thermodynamic conditions and needs further assessment.

A strong shear with enhances by increasing the mixing of tropospheric and cloudy air, as shown in Fig. 11h; for example, there is an increase by 30%–50% with respect to at . In terms of the total evaporative cooling, increases by approximately 50%, as indicated in Fig. 12. [The maximum rate of evaporative cooling decreases with increasing shear, but this decrease is overcompensated by a broadening of the evaporative cooling profile (cf. Fig. 10).] Shear, therefore, enhances the total amount of evaporative cooling , which is in contrast to the net radiative flux difference (and hence ) being constant. The shear enhancement of evaporative cooling within the cloud below the EIL is about 15% and compensates for the decrease of radiative cooling noted in section 5a. This explains why in-cloud turbulence remains approximately independent of the wind shear. All these observations stress the importance of evaporative cooling and shear effects on it.

Despite this substantial shear enhancement of evaporative cooling, however, we do not observe a cloud-top entrainment instability, understood as a runaway instability that leads to a rapid desiccation of the cloud. This concept is based on the positive feedback that exists between evaporative cooling and entrainment, since evaporative cooling enhances in-cloud turbulence, which in turn enhances entrainment and hence evaporation (Deardorff 1980b; Randall 1980). However, this feedback seems to be small because the ratio tends toward a constant value for , as shown in Fig. 12.

### d. Deformation effects

The deformation term arises from temporal changes in the shape of the mean buoyancy profile and allows us to distinguish between a quasi-steady and an unsteady state. For a quasi-steady state, is negligibly small, compared to the sum , which implies that the shape of the mean buoyancy profile changes slowly, compared to the deepening of the STBL. Hence, for a quasi-steady state, the magnitude of depends only weakly on , even though the different contributions to in Eq. (24) strongly vary with . In contrast, in an unsteady state, the deformation term is not negligible, and the shape of the mean buoyancy profile changes on time scales comparable to the turnover time of the large-scale turbulent motions in the STBL. Therefore, different choices of reference height can yield substantially different entrainment velocities.

In our simulations, the deformation term is comparable to the sum at , signaling an unsteady state. Deformation effects are indicated in Fig. 11a by the misalignment between contour lines of and the time evolution of the different reference heights (dashed lines), which is more pronounced for reference heights located near . This misalignment indicates that in an unsteady state, the partitioning between the individual contributions of changes significantly as function of and not only as function of height.

Shear weakens the EIL stratification and hence makes the cloud top more susceptible to deformations, since these deformations are created by in-cloud turbulent motions penetrating the stably stratified EIL (cf. section 4). Consistently, a strong shear with is observed to increase , which explains why in Fig. 9 and decrease more rapidly for , compared to . Shear effects on the deformation term are therefore important in unsteady regimes.

## 6. Summary and conclusions

Interactions of a mean vertical wind shear and turbulent convection driven by radiative and evaporative cooling of the stratocumulus cloud top have been studied by means of direct numerical simulations. We have focused on DYCOMS II–like conditions (i.e., subtropical stratocumulus with strong variations of specific humidity and static energy across the cloud top).

Shear effects are only found to be significant if the cloud-top velocity jump exceeds the critical velocity jump , where is the convective velocity scale characterizing in-cloud turbulence. For typical values of in the range of 0.2–0.9 m s^{−1} (Wood 2012), one finds (Δ*u*)_{crit} ≃ 1–4 m s^{−1}. For , shear enhances the entrainment of tropospheric air substantially. However, for , shear effects are negligible, and in-cloud turbulence penetrating the stably stratified entrainment interfacial layer (EIL) dominates the EIL dynamics. This threshold suggests that cloud-top shear associated with large-scale convective motions of the atmospheric boundary layer is unable to enhance cloud-top cooling significantly, since such a shear is typically characterized by velocity jumps commensurate with the convective velocity scale and .

Even for a strong wind shear with , shear effects are found to remain localized within the EIL (i.e., shear does not affect the in-cloud turbulence intensity). Shear only depletes the cloud, reduces the net radiative cooling, and weakens in-cloud turbulence if wind shear thickens the EIL substantially, compared to the cloud thickness. An analytic expression for the corresponding depletion velocity jump is provided, and (Δ*u*)_{dep} ≃ 3–10 m s^{−1} is found for a cloud thickness in the interval 100–200 m. The range of is consistent with measurement campaigns but is partly in disagreement with previous LES studies, where a depletion of the cloud is observed for shear velocities with . This difference is hypothesized to result from the spurious mixing associated with subgrid models and numerical artifacts.

Although cloud-layer properties (e.g., *w*_{*}) remain similar for all shear velocities investigated, a strong shear with enhances the entrainment velocity significantly. This enhancement has been studied by means of an integral analysis of the buoyancy equation, which provides an analytic decomposition of the entrainment velocity into four contributions. The radiative and evaporative cooling contributions and appear as the cumulative radiative and evaporative cooling above the reference height and not as the local cooling rates at this height. In contrast, the mixing contribution , which consists of the sum of a turbulent buoyancy flux and a molecular flux, is locally evaluated at the height . The deformation contribution describes changes in the shape of the mean buoyancy profile, which are generated by changes in the in-cloud turbulent convection penetrating the EIL.

The turbulent buoyancy flux contribution maximizes near the height of minimum turbulent buoyancy flux , which renders a significant contribution to at this height. As the reference height moves downward toward the height of zero mean buoyancy , the direct contributions from radiation and evaporation and monotonically increase, since more liquid water accumulates above and more mixing of tropospheric and cloudy air occurs above . We find that can significantly exceed and near , even though total evaporative cooling remains comparable to total radiative cooling (and hence to the reference buoyancy flux ) for all cloud-top velocity jumps investigated, namely, *E*_{0}/*B*_{1} ≃ 1.0–1.5. Therefore, at least for DYCOMS II–like conditions, the entrainment velocity might be well approximated by near but not so near , where the sum needs to be considered. This demonstrates that the partitioning of strongly depends on the choice of the reference height , even though the magnitude of does not in a quasi-steady state, and even though different definitions of the reference height only differ by a few meters. Entrainment-rate parameterizations should consistently reflect this dependence when estimating the different contributions to as needed in mixed-layer models.

A strong shear with enhances , , and significantly, indicating that shear effects should be considered in entrainment-rate parameterizations. For example, a strong shear with enhances the sum by 60%–100%, compared to the shear-free case at . This enhancement remains finite, however, and we do not observe a cloud-top entrainment instability whereby enhanced evaporative cooling leads to more in-cloud turbulence, which in turn promotes more evaporative cooling.

In unsteady cases, the deformation contribution is nonzero, and the magnitude of (and not only the various contributions to it) depends on the choice of the reference height. We find that decreases as the reference height moves upward in the direction of , and increases slightly for strong shear conditions (i.e., ). In general, deformation effects are argued to be important when the height of the atmospheric boundary layer varies significantly on time scales comparable to the large-eddy turnover time. Hence, deformations are expected to matter during cloud formation processes and during transients, such as during the transition from stratocumulus to shallow cumulus. Within those regimes, the deepening of the atmospheric boundary layer might not be well described by one single reference height.

We finally note that performing a conditional analysis of the cloud-top region might provide further insights into the entrainment processes, complementing the conventional analysis based on mean quantities described in this paper. For example, it has been observed that that entrainment operates differently in updraft and downdraft regions (Gerber et al. 2005), and incorporating that knowledge into mean quantities might help to further understand and parameterize the entrainment rates. We note, however, that both a conventional analysis based on mean quantities and a conditional analysis based on the distance to the cloud boundary are complementary since the former can be mathematically related to the latter, and results should therefore be complementary.

## Acknowledgments

We thank H. Gerber, S. Malinowski, B. Stevens, and C. Bretherton for motivating and instructive discussions on the topic. The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC). Funding was provided by the Max Planck Society through its Max Planck Research Groups program. Primary data and scripts used in the analysis and other supporting information that may be useful in reproducing the author’s work are archived by the Max Planck Institute for Meteorology and can be obtained by contacting publications@mpimet.mpg.de.

### APPENDIX A

#### Linearized Formulation

The formulation of the CTML is based on the conservation of the specific humidity and specific enthalpy *h*. These variables, as well as the liquid-water specific humidity , can be expressed in terms of three nondimensional variables , , and :

where the superscripts *c* and *d* refer to cloudy and dry air, respectively. The mixing fraction *χ* defines the hypothetical process of adiabatically mixing two air parcels in the mass fraction (Albrecht et al. 1985; Nicholls and Leighton 1986; Bretherton 1987). The variable describes diabatic deviations introduced by radiative effects (Moeng et al. 1995; Shao et al. 1997; Vanzanten 2002; Yamaguchi and Randall 2012), is the specific heat capacity of cloudy air, and is the temperature of cloudy air. The evolution equations for and are (Mellado et al. 2010; de Lozar and Mellado 2015; Mellado 2017)

where is the thermal diffusivity, and microphysical effects are neglected. Here, **u** is a velocity vector, and is the one-dimensional longwave radiative forcing based on Larson et al. (2007), with **k** being a unit vector pointing in the vertical direction. The net longwave radiative flux can be well approximated by

where is the net radiative flux cooling the cloud-top region, and is the extinction length.

Moreover, we apply the Boussinesq approximation to the Navier–Stokes equation

where *ν* refers to kinematic viscosity and *b* to buoyancy , with *ρ* being density. We assume that the Prandtl number is equal to 1 (i.e., ). To complete this set of equations, we still need an expression for the normalized liquid water and the buoyancy *b*. We can write analytic expressions for these variables when assuming phase equilibrium (infinitely fast thermodynamics) and linearizing the caloric and thermal equation of state. Under these assumptions, and *b* can be diagnosed from the prognostic variables according to

as discussed in Bretherton (1987), Pauluis and Schumacher (2010), and de Lozar and Mellado (2015). Here, and *ξ* are given by and , respectively, and the condition defines the saturation surface, which can be used to define the cloud boundary. The parameters and characterize how variations in enthalpy translate into changes of buoyancy and liquid water, respectively. The function *f* tends to a piecewise linear function in the limit , but has a finite second-order derivative of order , which is convenient for the numerical calculations. For our simulations, we apply since the obtained results become independent of