## Abstract

This paper discusses an extension of the approach proposed previously to represent the delay of cloud water evaporation and buoyancy reversal due to the cloud–environment mixing in bulk microphysics large-eddy simulation of clouds. In the original approach, an additional equation for the mean spatial scale of cloudy filaments was introduced to represent the progress toward microscale homogenization of a volume undergoing turbulent cloud–environment mixing, with the evaporation of cloud water allowed only when the filament scale approached the Kolmogorov microscale. Here, it is shown through a posteriori analysis of model simulations that one should also predict the volume fraction of the cloudy air that was diagnosed in the original approach. The resulting model of turbulent mixing and homogenization, referred to as the *λ*–*β* model, is applied in a series of shallow convection simulations using various spatial resolutions and compared to the traditional bulk model. This work represents an intermediate step in the development of a modeling framework to simulate characteristics of microphysical transformations during entrainment and subgrid-scale turbulent mixing.

## 1. Introduction

Recent modeling studies (e.g., Chosson et al. 2004, 2007; Grabowski 2006; Slawinska et al. 2008) demonstrate that assumptions concerning the microphysical evolution of natural clouds, in particular the homogeneity of cloud–environment mixing, significantly affect the albedo of a field of shallow convective clouds, such as subtropical stratocumulus and trade-wind cumulus. It follows that microphysical properties of such clouds have important implications for the role of clouds in the climate system. Because these clouds are strongly diluted by entrainment, the focus should be on modeling dynamical, thermodynamic, and microphysical processes associated with entrainment. This paper extends an approach for modeling subgrid-scale processes associated with entrainment and mixing proposed in Grabowski (2007, hereafter G07). As in G07, we limit the discussion to bulk representation of cloud microphysics and present analysis of a series of simulations of trade-wind shallow nonprecipitating convection observed during the Barbados Oceanographic and Meteorological Experiment (BOMEX; Holland and Rasmusson 1973) recently used in the model intercomparison study described in Siebesma et al. (2003); the same model setup was used in section 4b in G07. The longer-term goal is to combine the approach discussed here with the double-moment bulk microphysics scheme of Morrison and Grabowski (2007, 2008) to locally *predict* the homogeneity of mixing [i.e., the parameter *α* in Eq. (11) in Morrison and Grabowski 2008] using the results of a large set of direct numerical simulations (DNS) with detailed (bin) microphysics discussed in Andrejczuk et al. (2009).

The next section provides a brief summary of the method developed in G07 and presents a modification (or extension) of this approach necessitated by a posteriori analysis of simulations with and without the extension. Selected results from a series of numerical simulations of shallow nonprecipitating convection are then discussed in section 3. A brief summary and outlook in section 4 concludes this note.

## 2. Modeling evaporation of cloud water resulting from entrainment and mixing

The essence of the approach developed in G07 is to supplement the standard thermodynamic grid-averaged equations for the bulk advection–diffusion–condensation problem:

where *θ* is the potential temperature; *q _{υ}* and

*q*are the water vapor and cloud water mixing ratios, respectively;

_{c}*C*is the condensation rate;

*D*terms represent subgrid-scale turbulent transport terms; and all other symbols are exactly as in G07. These thermodynamic equations are supplemented with the evolution equation for the scale (or width) of cloudy filaments,

*λ*. The evolution of

*λ*is supposed to represent the progress of subgrid-scale turbulent mixing toward the microscale homogenization (Broadwell and Breidenthal 1982; Jensen and Baker 1989). When extended into the multidimensional framework and written in the conservative (flux) form analogous to (1), the equation for

*λ*takes the form

where the first term on the right-hand side (rhs) describes the decrease of *λ* as the turbulent mixing progresses [*ε* is the local dissipation rate of the turbulent kinetic energy (TKE) and *γ* ∼ 1 is a nondimensional parameter taken as *γ* = 1.8],^{1 }*S _{λ}* is the source/sink term, and

*D*is the subgrid transport term analogous to the

_{λ}*D*terms in (1). The source/sink term

*S*considers three processes that affect the scale

_{λ}*λ*: (i) initial formation of cloudy volumes due to grid-scale condensation, (ii) disappearance of cloudy volumes due to complete evaporation of cloud water, and (iii) homogenization of a cloudy volume. The source/sink term resets the current value

*λ*to either the scale comparable to the size of the grid box—say, Λ ≡ (Δ

*x*Δ

*y*Δ

*z*)

^{1/3}(where Δ

*x*, Δ

*y*, and Δ

*z*are the model grid lengths in the

*x*,

*y*, and

*z*directions, respectively)—or to zero. Because uniform cloudy grid boxes are characterized by

*λ*≡ Λ and cloud-free grid boxes have

*λ*≡ 0, the former is applied in cases (i) and (iii) and the latter in case (ii). Microscale homogenization of the cloudy grid box is assumed once the scale predicted by (2) falls bellow the threshold value

*λ*

_{0}, taken as 1 cm.

The overall motivation behind the approach is to represent the chain of events characterizing turbulent mixing—from the initial engulfment of the ambient fluid by an entraining eddy to the small-scale homogenization—and to include a corresponding delay in the saturation adjustment until the grid box can be assumed to be homogenized. This is schematically illustrated in Fig. 1. In practice, when *λ* = Λ or *λ* ≤ *λ*_{0}, the evaporation rate is exactly the same as in the traditional bulk model; that is, it is given by the saturation adjustment. However, when Λ > *λ* > *λ*_{0} (i.e., the turbulent mixing has not reached scales characterizing the small-scale homogenization), the adiabatic condensation rate (which depends on the resolved vertical velocity; see appendix in G07) is only used over the cloudy part of the grid box: *C* = *βC ^{a}*, where

*β*is the fraction of the grid box with cloudy air and

*C*is the adiabatic condensation rate.

^{a}It was suggested in G07 that *β* can be diagnosed locally based on the mean relative humidity of a gridbox RH and on the environmental relative humidity RH^{e} at this level as

where in the simulations and in the analysis here RH^{e} is taken from the initial RH profile, and additional limiting is used to avoid unphysical values of *β*. Arguably, (3) may provide a reasonable approach for the case of a convective cloud. For stratocumulus, on the other hand, entrainment and cloud–environment mixing takes place primarily at the cloud top where vertical gradients of environmental profiles are large because of the presence of boundary layer inversion. It follows that (3) is most likely of limited use when modeling stratocumulus. Because the accuracy of (3) is uncertain, we decided to make *β* an additional model variable and to locally predict its value together with *λ*. The advection–diffusion equation for *β* is

where *S _{β}* is the source/sink term and

*D*is the subgrid transport term. Similarly to the source term for

_{β}*λ*in (2),

*S*resets

_{β}*β*to 1 if either grid-scale condensation or microscale homogenization takes place. Otherwise,

*β*is advected and diffused in a manner similar to

*λ*. Additional justification for the approach in which

*β*is predicted by (4) rather than diagnosed from (3) comes from a posteriori comparison between the values of

*β*obtained from the diagnostic formulation (3) and the prognostic formulation (4) discussed in the next section.

## 3. Application of the λ–β model to BOMEX shallow convection

Equations (3) and (4) were included in the anelastic semi-Lagrangian/Eulerian cloud model EULAG, documented in Smolarkiewicz and Margolin (1997; model dynamics), Grabowski and Smolarkiewicz (1996; model thermodynamics), and Margolin et al. (1999; subgrid-scale model) applied in G07. The Eulerian version of the model is used to simulate the quasi-steady-state trade-wind shallow nonprecipitating convection observed during BOMEX (Holland and Rasmusson 1973; Siebesma et al. 2003), as in G07. In this case, the 1.5-km-deep trade-wind convection layer overlays the 0.5-km-deep mixed layer near the ocean surface and is covered by the 500-m-deep trade-wind inversion layer. The cloud cover is about 10% and quasi-steady conditions are maintained by prescribed large-scale subsidence, large-scale moisture advection, surface heat fluxes, and radiative cooling.

In the three-dimensional simulations presented here, the model setup is as described in Siebesma et al. (2003) but different domain sizes and model grid lengths are applied (i.e., increasing the model resolution but keeping the same number of grid points in the horizontal, 128 × 128; and adjusting the number of grid points in the vertical to maintain the 3-km vertical extent of the domain). Three different model grid lengths were considered: 100 m × 40 m (horizontal × vertical) i.e., as in Siebesma et al. 2003 and in G07), 50 m × 40 m, and 25 m × 25 m. Note that decreasing the horizontal grid length with the same number of grid points results in progressively decreasing horizontal extent of the domain and thus poorer cloud statistics. However, this does not seem to affect the results discussed below. Results from both the traditional bulk model and the *λ*–*β* model will be presented. As in Siebesma et al. (2003) and G07, the model is run for 6 h and snapshots of model results for hours 2 to 6, archived every 4 min, are used in the analysis.

Figure 2 shows the results from the simulation applying the *λ*–*β* model with a grid length of 25 m in the horizontal and vertical directions. The figure compares values of *β* predicted by (4) to corresponding values obtained using the diagnostic relationship (3) at the 1250-m level (i.e., about half of the cloud field depth). As the figure illustrates, there is a significant difference between the predicted and diagnosed values of *β*, most likely reflecting the fact that assumptions leading to (3) are seldom justified. More importantly, however, the diagnostic approach leads to a significant overestimation of *β* relative to that obtained from the prognostic approach. This is because the air entrained into a cloud is typically more humid than the environmental profile; that is, the environmental air is premoistened before it is entrained into a cloud. This is consistent with the fact that reducing RH^{e} in (3) indeed results in an increase of the diagnosed *β*. One way to show that the entrained air is premoistened is to consider the plot of the ratio RH/RH^{e} versus *β*. Equation (3) implies that RH/RH^{e} → 1 (i.e., RH approaches RH^{e}) when *β* → 0. To investigate this, we plot RH/RH^{e} versus *β* in Fig. 3. The figure shows that the ratio RH/RH^{e} is typically larger than 1 when *β* is close to 0. This implies that the air entrained into the cloud is, on average, more humid than environmental air at this level. Similar analyses can be performed separately for the temperature and water vapor mixing ratio. They show that the water vapor is typically higher near the cloud than the environmental value at a given level, whereas the temperature close to the cloud is usually lower than the environmental temperature. The latter is consistent with the fact that premoistening is associated with the evaporative cooling.

Another way to demonstrate that the air in the vicinity of clouds is typically more humid than in the far environment is to plot RH/RH^{e} as a function of the distance to the nearest cloud. Such plots are shown in Fig. 4 for both the traditional bulk and the *λ*–*β* models and the simulations applying 25-m grid lengths. As the figure shows, RH in the vicinity of a cloud is indeed higher than further away. In general, this is reminiscent of both numerical (Jonker et al. 2008; Heus and Jonker 2008) and observational (Heus et al. 2009) studies demonstrating the existence of a thin shell of subsiding air near edges of shallow convective clouds, with thermodynamic properties different than the far environment. Figure 4 also shows that the cloud edge (i.e., the zero distance in Fig. 4) corresponds to only one value of RH for the traditional bulk model (i.e., the saturated one), whereas significant scatter exists for the *λ*–*β* model. The latter is consistent with the expectation that a grid box at the cloud edge can be partially cloudy with RH below saturation. In summary, the above results provide a posteriori support for the approach where *β* is predicted using (4) rather then diagnosed from (3).

Figure 5 presents contoured frequency by altitude diagrams (CFADs) of the vertical velocity for simulations with the traditional bulk model and the *λ*–*β* model using the 25-m grid lengths. Only grid points with a cloud water mixing ratio larger than 10^{−2} g kg^{−1} are included in the analysis. Although both CFADs are similar, there are some differences. First, clouds in the *λ*–*β* model seem slightly deeper than in the bulk approach. This is consistent with the lower-resolution results discussed in section 4b of G07 and can be argued to result from delayed evaporation of cloud water due to entrainment and mixing, resulting in more cloud buoyancy. Second, in the *λ*–*β* model, points with vertical velocities in the range of 0 to 1 m s^{−1} contribute more to the distribution than points in the range of 1 to 2 m s^{−1} across most of the depths of the cloud field. This is consistent with the expectation that cloud evaporation (and thus buoyancy reversal and subsequent downward acceleration) is delayed when the *λ*–*β* approach is used.

To document the impact of the *λ*–*β* approach on the vertical velocity field inside clouds, Fig. 6 shows a scatterplot of the vertical velocity versus *λ* at an altitude of 1250 m for the simulation using 25-m grid lengths (i.e., as in Figs. 2 –5). It needs to be kept in mind that the cloud water in a grid box with *λ*_{0} < *λ* < Λ = 25 m would immediately evaporate in the traditional bulk model. As Fig. 5 shows, there are many grid boxes with intermediate values of *λ* and they are characterized by small positive and negative values of the vertical velocity, with the mean around zero. A clear tendency toward positive bias is apparent for *λ* values approaching Λ. The values at *λ* = Λ correspond to the homogeneous cloudy grid boxes and thus are characterized by both positive and negative values, with a bias toward positive values at this particular height.

Figure 7 shows profiles of the horizontally averaged (i.e., over the entire horizontal domain) cloud water mixing ratios (also averaged in time between hours 2 and 6) for all simulations (i.e., 100 m × 40 m, 50 m × 40 m, and 25 m × 25 m for *λ*–*β* and traditional models). Spatial resolution impacts the depth of the cloud layer, and this impact is significantly stronger for the *λ*–*β* model. Arguably, the latter is consistent with the delayed cloud water evaporation, which allows clouds to longer maintain positive buoyancy despite the dilution by environmental air. Note that the cloud water profiles near the cloud base and the height of the cloud base are similar in all simulations. The cloud water increases with height for the *λ*–*β* model with the lowest resolution, but for higher-resolution simulations this is reversed, in agreement with the bulk model results. A local maximum in the cloud water profiles at the bottom of the trade-wind inversion (i.e., just above 1.5-km height) is significantly stronger for the *λ*–*β* model than for the bulk model. This is again consistent with more gradual evaporation of cloud water once the vertical development of a cloud is arrested by the inversion. The differences between the traditional and *λ*–*β* models decrease as the model resolution increases, as one might expect.

The difference in the cloud field depth for low- and high-resolution *λ*–*β* simulations can be understood using the following argument, supported by selected snapshots of cloud field shown in Fig. 8. Entrainment and cloud dilution in numerical simulation is a combination of resolved and parameterized dynamics. When model spatial resolution increases, entrainment becomes better resolved. As Fig. 8 shows, the cloud in the low-resolution simulation is wide and features just a single updraft, and its cloud–environment interface is relatively smooth. In contrast, clouds in the high-resolution simulation typically contain multiple turrets and show complicated cloud–environment structures responsible for the entrainment. Arguably, clouds in high-resolution simulation are more realistic, and it is reassuring to see that the differences between the traditional and *λ*–*β* models diminish as the model spatial resolution increases.

## 4. Summary and outlook

This paper presents results using an approach developed in G07 to represent the delay of cloud water evaporation and buoyancy reversal associated with the cloud–environment mixing. This is accomplished by including an equation for the scale (width) of cloudy filaments, *λ*. A conceptual picture assumes that turbulent mixing is initiated by an engulfment of environmental fluid by an entraining eddy. It progresses through the gradual filamentation of the initial coarse mixture of cloudy and cloud-free air (represented by the decrease of *λ*) and reaches the final homogenization once *λ* approaches spatial scales close to the Kolmogorov microscale. Using this approach, evaporation of cloud water due to saturation adjustment is delayed until the microscale homogenization stage. This is schematically illustrated in Fig. 1. Here, we extend this approach and include an additional prognostic model variable, the fraction of the grid box covered by the cloudy air, *β*. The parameter *β* was diagnosed in G07, but here we show that the diagnostic approach results in significant inaccuracies when modeling shallow convective clouds. These inaccuracies exist because the air entrained into the cloud has thermodynamic properties significantly different from the environmental air at the entrainment level, in contrast to the assumptions underlying the diagnostic approach (3). The diagnostic approach is even less likely to work in the case of stratocumulus, where steep gradients of environmental profiles in the vicinity of boundary layer inversion make the validity of (3) questionable. We refer to the extended approach as the *λ*–*β* model.

As far as the representation of entrainment and mixing is concerned, it would be desirable to compare predictions of the *λ*–*β* model with very-high-resolution numerical simulations and/or with cloud observations. For the modeling, the rate of decrease of the filament scale can be deduced from DNS of the cloud–clear air interfacial mixing of the type discussed in Andrejczuk et al. (2004, 2006, 2009). Although results of DNS simulations qualitatively agree with the Broadwell and Breidenthal model, the dynamic range of these simulations (i.e., the ratio between the scale at which energy is introduced to the system and the scale at which it is dissipated) is insufficient for a quantitative comparison between predictions of (2) and the DNS. Applying the linear eddy model (LEM; Kerstein 1988, 1991) as in Krueger (1993) might be more practical. Although the parameter *γ* in (2) was selected based on LEM, a detailed comparison between the prediction of the filament scale evolution given by (2) and the results of LEM, involving idealized mixing events, would still be useful. For the observations, one may attempt to compare statistical characteristics of *λ* and *β* predicted by the model with very-high-spatial-resolution aircraft observations, for instance such as reported in Gerber et al. (2008). If successful, such comparisons would provide strong support for the *λ*–*β* model.

Application of the *λ*–*β* model to the same BOMEX case as in G07 shows that predicting *β* in addition to *λ*, rather than diagnosing it, results in small changes of mean cloud properties, such as profiles of the cloud water content and cloud fraction, profiles of the water vapor and cloud water fluxes, liquid water path, cloud cover, etc. In other words, bulk properties of simulated clouds seem unaffected whether *β* is predicted or diagnosed. However, an important motivation to accurately predict *β*—not obvious in the current study—is that *β* will play an important role when the *λ*–*β* model is combined with the double-moment bulk microphysics scheme of Morrison and Grabowski (2007, 2008). The overall strategy for predicting spectral changes of cloud droplets is to locally diagnose the homogeneity of mixing and to decrease only the number of droplets for the extremely inhomogeneous mixing, only their sizes in the case of homogeneous mixing, or both the number and the sizes for the intermediate mixing. We have recently developed an approach to relate the homogeneity of mixing to the local intensity of turbulence, the humidity of the entrained air, and the size of cloud droplets (Andrejczuk et al. 2009). Combining all these developments is the subject of ongoing research.

## Acknowledgments

This work was partially supported by NOAA Grants NA05OAR4310107 and NA08OAR-4310543 (DJ and WWG), NCAR’s Geophysical Turbulence Program (DJ), and by Polish MNiSW Grant ACCENT/295/2006 (DJ and HP). Critical comments on the manuscript by Steve Krueger and an anonymous reviewer are acknowledged.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Wojciech W. Grabowski, NCAR/MMM, P.O. Box 3000, Boulder, CO 80307–3000. Email: grabow@ucar.edu

^{1}

This value was suggested by Prof. Steven Krueger based on the theoretical argument and his simulations applying the linear eddy model (Kerstein 1988, 1991; Krueger 1993), a 1D analog of turbulent stirring and molecular diffusion. As mentioned in G07, LES results did show some sensitivity to the value of this parameter in low-resolution simulations reported in G07. One might argue that this sensitivity should diminish as model resolution is increased because of a better representation of entrainment dynamics with higher spatial resolution. Also note that this parameter was marked *α* in G07. We change the notation to avoid conflict with parameter *α* that describes the homogeneity of mixing in Morrison and Grabowski [2008, Eq. (11)].

* The National Center for Atmospheric Research is sponsored by the National Science Foundation.