Influences of an entrainment–mixing parameterization on numerical simulations of cumulus and stratocumulus clouds

Different entrainment–mixing processes can occur in clouds; however, a homogeneous mixing mechanism is often implicitly assumed in most commonly used microphysics schemes. Here, we first present a new entrainment–mixing parameterization that uses the grid mean relative humidity without requiring the relative humidity of the entrained air. Then, the parameterization is implemented in a microphysics scheme in a large eddy simulation model, and sensitivity experiments are conducted to compare the new parameterization with the default homogeneous entrainment–mixing parameterization. The results indicate that the new entrainment–mixing parameterization has a larger impact on the number concentration, volume mean radius, and cloud optical depth in the stratocumulus case than in the cumulus case. This is because inhomogeneous and homogeneous mixing mechanisms dominate in the stratocumulus and cumulus cases, respectively, which is mainly due to the larger turbulence dissipation rate in the cumulus case. Because stratocumulus clouds break up during the dissipation stage to form cumulus clouds, the effects of this new entrainment–mixing parameterization during the stratocumulus dissipation stage are between those during the stratocumulus mature stage and the cumulus case. A large aerosol concentration can enhance the effects of this new entrainment–mixing parameterization by decreasing the cloud droplet size and evaporation timescale. The results of this new entrainment–mixing parameterization with grid mean relative humidity are validated by the use of a different entrainment–mixing parameterization that uses parameterized entrained air properties. This study sheds new light on the improvement of entrainment–mixing parameterizations in models.

1 Introduction

The process of entrainment and subsequent mixing between clouds and their environment is one of the most uncertain processes in cloud physics, which is thought to be crucial to many outstanding issues, including warm rain initiation and subsequent precipitation characteristics, cloud–climate feedback, and evaluating the indirect effects of aerosol (Paluch and Baumgardner, 1989; Yum, 1998; Ackerman et al., 2004; Kim et al., 2008; Huang et al., 2008; Del Genio and Wu, 2010; Lu et al., 2011, 2014; Kumar et al., 2013; Zheng and Rosenfeld, 2015; Fan et al., 2016; Gao et al., 2020, 2021; Zhu et al., 2021; Xu et al., 2021; Yang et al., 2016, 2021). The most well-studied concepts are homogeneous/inhomogeneous entrainment–mixing mechanisms. During homogeneous mixing, all droplets experience evaporation, and no droplet is evaporated completely. During extremely inhomogeneous mixing, some droplets near the entrained air evaporate completely, while the remaining droplets maintain their original sizes. If the situation is somewhere between these two extreme scenarios, an inhomogeneous mixing process occurs. Some studies suggest that homogeneous mixing is likely to be typical (Jensen et al., 1985; Burnet and Brenguier, 2007; Lehmann et al., 2009), whereas others have claimed that an extremely inhomogeneous scenario is dominant (Pawlowska et al., 2000; Burnet and Brenguier, 2007; Haman et al., 2007; Freud et al., 2008, 2011). Different mechanisms can be undistinguishable when the relative humidity in the entrained air is high (Gerber et al., 2008).

Some sensitivity studies assuming homogeneous or extremely inhomogeneous mixing have found that different mixing mechanisms can significantly influence the microphysics and radiative properties of clouds (Lasher-Trapp et al., 2005; Grabowski, 2006; Chosson et al., 2007; Slawinska et al., 2008). For example, Grabowski (2006) used a cloud-resolving model and found that the amount of solar energy reaching the surface in the pristine case, assuming the homogeneous mixing scenario, is the same as in the polluted case with extremely inhomogeneous mixing. This result was verified by Slawinska et al. (2008) using a large-eddy simulation (LES) model. Although the influence of different mixing mechanisms in simulations is lower when two-moment microphysics schemes are used (Hill et al., 2009; Grabowski and Morrison, 2011; Slawinska et al., 2012; Xu et al., 2020), Hill et al. (2009) also claimed that there are still many uncertainties in the entrainment–mixing process, and the effect of different mixing mechanisms can be more important over the entire cloud life-cycle.

In recent years, methods have been developed to describe general entrainment–mixing processes, with homogeneous and extremely inhomogeneous scenarios as special cases (Andrejczuk et al., 2006, 2009; Lehmann et al., 2009; Lu et al., 2011). Hoffmann et al. (2019) and Hoffmann and Feingold (2019) conducted LESs at the subgrid scale with turbulent mixing using a linear eddy model. Andrejczuk et al. (2009) used the results of direct numerical simulation (DNS) to establish a relationship between instantaneous microphysical properties and Damköhler number (D_a; Burnet and Brenguier, 2007) and developed a parameterization of the entrainment–mixing process. Lu et al. (2013) developed a parameterization of the entrainment–mixing process based on the relationship between the homogeneous mixing degree (ψ) and transition scale number (N_L) in the explicit mixing parcel model (EMPM), as well as aircraft observation data. Gao et al. (2018) investigated how ψ is related to D_a and N_L in a DNS to improve the parameterization of the entrainment–mixing process. Luo et al. (2020) simulated more than 12 000 cases with EMPM by changing a variety of parameters affecting entrainment–mixing processes and developed a parameterization that improved the one proposed by Lu et al. (2013).

Although several entrainment–mixing parametrizations have been proposed, to the best of our knowledge, only one study (Jarecka et al., 2013) has coupled an entrainment–mixing parameterization with cloud microphysics to consider the change in cloud droplet concentration during the entrainment–mixing process. Jarecka et al. (2013) applied an entrainment–mixing parameterization, in terms of the Damköhler number, to a two-moment microphysics scheme and found small impacts of entrainment–mixing parameterization in shallow cumulus clouds. To further explore the influences of entrainment–mixing processes, this study first modifies the entrainment–mixing parameterization in terms of the transition scale number proposed by Luo et al. (2020) to couple it more easily with microphysics schemes. The parameterization is then implemented in the two-moment Thompson aerosol-aware scheme (Thompson and Eidhammer, 2014). Finally, the effects of parameterization on the physical properties of clouds are examined in both cumulus and stratocumulus clouds.

The rest of this paper is organized as follows: Sect. 2 describes the new entrainment–mixing parameterization, simulated cases, and modeling setup. The major results are presented and discussed in Sect. 3. The influences of the new entrainment–mixing parameterization on cloud physics and the underlying mechanisms are examined, and the effects of turbulence dissipation rate (ε) and aerosol concentration are also discussed. Some concluding remarks are presented in Sect. 4.

2 Parameterization, simulated cases, and modeling setup

2.1 The new entrainment–mixing parameterization

According to Morrison and Grabowski (2008), the effect of the entrainment–mixing process on cloud microphysical properties can be expressed as follows:

$$\begin{array}{}\text{(1)}& N_{\mathrm{c}}=N_{\mathrm{c}\mathrm{0}}{\left({\displaystyle \frac{q_{\mathrm{c}}}{q_{\mathrm{c}\mathrm{0}}}}\right)}^{\mathit{\alpha }},\end{array}$$

where N_c and N_c0 are the cloud droplet number concentrations after and before the evaporation process, respectively, and q_c and q_c0 represent the corresponding cloud water mixing ratios. It is noteworthy that when a new saturation is achieved after evaporation, q_c is determined by q_c0, relative humidity (RH), air pressure, and temperature. The parameter α can be pre-set to any value between 0 (homogeneous mixing) and 1 (extremely inhomogeneous mixing) to represent a different degree of subgrid-scale mixing homogeneity. In this study, instead of specifying α as a predetermined constant, here it is determined through the following expressions (Lu et al., 2013; Luo et al., 2020):

$$\begin{array}{}\text{(2a)}& {\displaystyle }\mathit{\alpha }=\mathrm{1}-\mathit{\psi },\text{(2b)}& {\displaystyle }\mathit{\psi }=c\exp \left(a{N}_{\mathrm{L}}^b\right),\end{array}$$

where a, b, and c are the three fitting parameters (Luo et al., 2020). The dimensionless number N_L is a dynamical measure of the degree of subgrid-scale mixing homogeneity (Lu et al., 2011) defined by

$$\begin{array}{}\text{(3a)}& {\displaystyle }{N}_{\mathrm{L}}={\displaystyle \frac{L^{\ast }}{\mathit{\eta }}},\text{(3b)}& {\displaystyle }\mathit{\eta }={\left({\mathit{\nu }}^{\mathrm{3}}/\mathit{\epsilon }\right)}^{\mathrm{1}/\mathrm{4}},\text{(3c)}& {\displaystyle }{L}^{\ast }={\mathit{\epsilon }}^{\mathrm{1}/\mathrm{2}}{\mathit{\tau }}_{\mathrm{evap}}^{\mathrm{3}/\mathrm{2}},\end{array}$$

where L^∗ is the transition length (Lehmann et al., 2009), η is the Kolmogorov microscale, ν is the kinematic viscosity, and ε is calculated from the subgrid turbulent kinetic energy (Deardorff, 1980):

$$\begin{array}{}\text{(4)}& \mathit{\epsilon }=C{E}^{\mathrm{3}/\mathrm{2}}/L,\end{array}$$

where C=0.70 is an empirical constant, E is the subgrid turbulent kinetic energy, and L is the model grid size. The evaporation timescale (τ_evap) is defined as the time taken for droplets to evaporate completely in an unsaturated environment and is calculated as

$$\begin{array}{}\text{(5a)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{evap}}=-{\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{2}A{S}_{\mathrm{e}}}},\text{(5b)}& {\displaystyle }A={\displaystyle \frac{\mathrm{1}}{\left[\left(\frac{L_{\mathrm{h}}}{R_{\mathrm{v}}T}-\mathrm{1}\right)\frac{L_{\mathrm{h}}{\mathit{\rho }}_{\mathrm{L}}}{KT}+\frac{{\mathit{\rho }}_{\mathrm{L}}{R}_{\mathrm{v}}T}{D{e}_{\mathrm{s}}\left(T\right)}\right]}},\end{array}$$

where r is the volume mean radius of cloud droplets, A is a function of pressure and temperature, S_e is the supersaturation (RH − 1) of entrained air, L_h is the latent heat, R_v is the specific gas constant for water vapor, T is air temperature, ρ_L is the density of liquid water, K is the coefficient of thermal conductivity of air, D is the diffusion coefficient of water vapor in the air, and e_s(T) is the saturation vapor pressure over a plane water surface at temperature T.

Unfortunately, S_e in Eq. (5a) is generally unavailable in atmospheric models, including LES models. Thus, the entrainment–mixing parameterization developed by Luo et al. (2020) based on the properties of entrained air cannot be used directly. To solve this problem, we modify the entrainment–mixing parameterization of Luo et al. (2020) by replacing S_e with the domain mean RH in the EMPM, after entrainment but before evaporation, based on 12 218 cases:

$$\begin{array}{}\text{(6)}& \mathit{\psi }=\mathrm{107.19}\exp \left(-\mathrm{1.99}{N}_{\mathrm{L}}^{-\mathrm{0.29}}\right).\end{array}$$

Figure 1 shows the fitting results of the modified new entrainment–mixing parameterization. Compared to the parametrization proposed by Luo et al. (2020), the modified parameterization has similar ψ−N_L distributions but with a larger N_L for the same ψ because the EMPM domain mean RH is larger than the entrained air RH. With this modification, N_L, ψ, and thus the effect of the entrainment–mixing processes on droplet concentration can be directly calculated using the LES grid mean RH. It is important to note that the parameterization does not mean that the entrained air RH is equal to that of the LES grid mean RH. It is also worth noting that a wide range of ε, S_e, and fraction of entrained air (f) are taken into account when establishing the parameterization with the EMPM. The details of the EMPM simulations and related calculations are provided by Luo et al. (2020).

Figure 1 Parameterization of cloud entrainment–mixing mechanisms by relating homogeneous mixing degree (ψ) to transition scale number (N_L) from EMPM. The contours represent the joint probability distribution function (PDF) of ψ vs. N_L. The magenta dots and error bars are mean values and standard deviations of ψ in each N_L bin, respectively. The mean values are fitted using a weighted least squares method with the number of data points in each N_L bin as with the weight. The fitting equation, coefficient of determination (R²), and p values are also given. N_L is calculated with the domain-averaged relative humidity after entrainment but before evaporation in the EMPM.

2.2 LES model, simulation cases, and modeling setup

The LES model is built by applying the large-scale forcing module presented in Endo et al. (2015) to the Weather Research and Forecasting (WRF) model tailored for solar irradiance forecasting (WRF-Solar; Hacker et al., 2016; Haupt et al., 2016). The large-scale forcing data (VARANAL) used in this process are derived from the constrained variational analysis (CVA) approach developed by Zhang et al. (2001) and provided by the U.S. Department of Energy's Atmospheric Radiation Measurement Program (http://www.arm.gov, last access: 20 March 2022). The modified entrainment–mixing parameterization is implemented in the two-moment Thompson aerosol-aware scheme (Thompson and Eidhammer, 2014).

To investigate the behaviors of the new entrainment–mixing parameterization in different cloud types, cumulus and stratocumulus cases are simulated. For both the cumulus and stratocumulus cases, the horizontal resolution of the model is 100 m × 100 m with a domain area of 14.4 km × 14.4 km. The vertical direction is divided into 225 layers with a resolution of 30 m.

Table 1 Summary of names and corresponding descriptions of the four experiments for each case of cumulus and stratocumulus. The meaning of each symbol for each experiment can be found in the text.

For each cloud case, ψ is first set to 1 for the default experiment because most LES models assume a homogeneous entrainment–mixing mechanism. The simulation with the new entrainment–mixing parameterization (Eqs. 1–6) is hereafter referred to as new. First, N_L is diagnosed for each grid, and ψ is then calculated using Eq. (6). Finally, the variation in N_c during entrainment–mixing is obtained using Eqs. (1) and (2a). To examine the influence of the aerosol number concentration on the entrainment–mixing process, we conduct the numerical experiments default_10 and new_10 by multiplying the initial aerosol number concentrations for the default and new models, respectively, by a factor of 10. Thus, four sets of numerical experiments are conducted for both the cumulus and stratocumulus cases; the names of the experiments and corresponding descriptions are summarized in Table 1.

3 Results

3.1 Cumulus case

For the cumulus case, the simulation starts at 12:00 UTC on 11 June 2016 and ends at 03:00 UTC on 12 June 2016 with an output interval of 10 min and spin times of 3 h. To demonstrate the utility of the model, Fig. 2 compares the temporal evolution of the observed and simulated cloud fraction (Fig. 2a) and solar irradiance (Fig. 2b) from the default experiment. Grid points with q_c larger than 0.01 g kg⁻¹ are defined as “cloudy areas”. Also shown for comparison are observational data with a 1 h temporal resolution, which are provided by the LES Atmospheric Radiation Measurement Symbiotic Simulation and Observation (LASSO) campaign (Gustafson et al., 2020). The observations show that the cloud forms at 12:00 UTC on 11 June and dissipates completely by 01:00 UTC on 12 June with a maximum cloud fraction of 0.47 at 16:00 UTC on 11 June. Considering the difference between the solar irradiances obtained from point measurements and the value representing the simulation domain, the observed solar irradiance at the Southern Great Plains (SGP) Central Facility are compared with the results of the central grid point in simulation (Fig. 2b). Evidently, although the results of the simulation do not fluctuate as much as the observations, the model captures the general behaviors of both cloud fraction and solar irradiance. The general agreement between the simulations and observations lends credence to using the model in further study.

Figure 2 Time series of (a) domain-averaged cloud fraction and (b) total downward irradiance at the central point from the observation and the default experiment in the cumulus case.

Figure 3 The temporal evolutions of main cloud microphysical and optical properties in all simulation experiments for the cumulus case, including (a) cloud water mixing ratio (q_c) (g kg⁻¹), (b) cloud droplet number concentration (N_c) (cm⁻³), (c) cloud droplet volume mean radius (r_v) (µm), (d) cloud water path (CWP) (g m⁻²), and (e) cloud optical depth (τ). In (b), N_c values in the experiments default and new are normalized by the maximum cloud droplet concentration (N_cmax) from default; N_c values in the experiments default_10 and new_10 are normalized by N_cmax from default_10. The four experiments are detailed in Table 1.

Figure 3 shows the evolution of the microphysical and optical properties of clouds in the cloudy areas of all simulation experiments, including q_c, N_c, droplet volume mean radius (r_v), cloud water path (CWP), and cloud optical depth (τ). To visually and simultaneously compare the change in cloud droplet concentration under different aerosol concentrations, the maximum cloud droplet concentration (N_cmax) from default is used to normalize N_c in default and new, while N_cmax from new_10 is used to normalize N_c in default_10 and new_10. The CWP is calculated as follows:

$$\begin{array}{}\text{(7)}& \mathrm{CWP}=\underset{\mathrm{0}}{\overset{H}{\int }}{\mathit{\rho }}_{\mathrm{a}}{q}_{\mathrm{c}}\left(z\right)\mathrm{d}z,\end{array}$$

where ρ_a is the air density, q_c(z) is the cloud water mixing ratio at each height (z), and H is the cloud thickness. The optical depth τ is estimated with

$$\begin{array}{}\text{(8)}& \mathit{\tau }={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}}}\underset{\mathrm{0}}{\overset{H}{\int }}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{q}_{\mathrm{c}}\left(z\right)}{r_{\mathrm{e}}\left(z\right)}}\mathrm{d}z,\end{array}$$

where ρ_w is the water density, and r_e(z) is the effective radius of the cloud droplets at each height (z). The time-averaged values of these physical properties of the clouds are listed in Table 2 for convenience.

Table 2 Summary of the case mean values of the key quantities in all the simulations of the cumulus case, containing cloud water mixing ratio (q_c), cloud droplet number concentration (N_c), cloud droplet volume mean radius (r_v), cloud water path (CWP), and cloud optical depth (τ). The experiments are detailed in Table 1.

For the low aerosol number concentration, the simulations with the new entrainment–mixing parameterization have smaller N_c (35.53 cm⁻³) and larger r_v (13.29 µm) than the default homogeneous simulation (35.78 cm⁻³ for N_c and 13.27 µm for r_v in default). However, comparing new to default, the relative changes in N_c, r_v, and τ are very small. When the aerosol concentration increases 10-fold (default_10 and new_10), q_c, CWP, and τ increase according to the aerosol indirect effect (Peng et al., 2002; Wang et al., 2019; Li et al., 2011; Wang et al., 2011). Meanwhile, r_v decreases significantly owing to the larger cloud number concentration. The effects of the new entrainment–mixing parameterization also increase; for example, the change in N_c increases from −0.70 % (new compared to default) to −2.74 % (new_10 compared to default_10), r_v increases from +0.15 % to +0.57 %, and τ increases from −0.38 % to −0.58 %; the reasons for these changes are discussed later. These small changes are similar to those identified in previous cumulus studies (Jarecka et al., 2013; Hoffmann et al., 2019).

3.2 Stratocumulus case

The stratocumulus case is simulated from 09:00 UTC on 19 April 2009 to 03:00 UTC on 20 April 2009; the first 3 h are set to be spin-up times. We examine the stratocumulus region of the cloud base at ∼ 2.1 km and the cloud top at ∼ 2.3 km (cloud thickness of ∼ 200 m). Figure 4 shows the time series of the domain-averaged cloud fraction and total downward irradiance at the central point in the observation and the default experiment from 12:00 to 24:00 UTC. Similar to the cumulus case, the simulations compare favorably with the observations, which further reinforces the utility of the LES model. The observed data show that the cloud fraction increases with time and peaks at 16:00 UTC. The simulated cloud fraction has a value of 1 before 16:00 UTC, fluctuates from 16:00 to 21:00 UTC, and decreases sharply after 21:00 UTC. This period can be divided into three stages, namely the mature stage, pre-dissipation stage, and dissipation stage.

Figure 4 Time series of (a) domain-averaged cloud fraction and (b) total downward irradiance at the central point from the observation and the default experiment in the stratocumulus case.

Figure 5 The temporal evolutions of main cloud microphysical and optical properties in all simulation experiments for the stratocumulus case, including (a) cloud water mixing ratio (q_c) (g kg⁻¹), (b) cloud droplet number concentration (N_c) (cm⁻³), (c) cloud droplet volume mean radius (r_v) (µm), (d) cloud water path (CWP) (g m⁻²), and (e) cloud optical depth (τ). In (b), N_c values in the experiments default and new are normalized by the maximum cloud droplet number concentration (N_cmax) from default; N_c values in the experiments default_10 and new_10 are normalized by N_cmax from default_10. The four experiments are detailed in Table 1.

As with the cumulus case, the temporal evolutions of the physical properties (q_c, N_c, r_v, CWP, and τ) of the clouds are shown in Fig. 5. In contrast to the oscillating changes exhibited by the physical quantities in the cumulus case (Fig. 3), the physical properties in the stratocumulus case exhibit a mostly smooth temporal evolution. Furthermore, default and new exhibit clear distinctions during the early periods, but these differences decrease during the dissipation stage. This is also the case with default_10 and new_10.

Table 3 Summary of the case mean values of the key quantities in all the simulations of the stratocumulus case, including cloud water mixing ratio (q_c), cloud droplet number concentration (N_c), cloud droplet volume mean radius (r_v), cloud water path (CWP), and cloud optical depth (τ). The numbers in and out of the parentheses are the results at the mature and dissipation stages, respectively. The experiments are detailed in Table 1.

To compare the different behaviors of the simulation experiments at different stages, the results at the mature and dissipation stages are analyzed in detail. The mean values of the main microphysical and optical properties of the clouds are summarized in Table 3. As expected, the cloud microphysical and optical properties at the mature stage are all larger than those at the dissipation stage. The effects of the new entrainment–mixing parametrization are also more significant at the mature stage. Compared to default, the new model results in a 7.36 % smaller N_c, 3.20 % larger r_v, and 5.98 % smaller τ during the mature stage. During the dissipation stage, the changes in N_c, r_v, and τ are −4.76 %, +2.12  %, and −2.56 %, respectively. The largest influence of the new entrainment–mixing parametrization occurs during the mature stage when the aerosol concentration is 10 times greater. The differences in N_c, r_v, and τ between new_10 and default_10 are −9.69 %, +3.88 %, and −5.85 %, respectively, averaged over the mature stage. These differences are much larger than those reported by Hill et al. (2009), who found that assuming extremely inhomogeneous mixing has a negligible effect on stratocumulus simulations. Our results also prove the speculation of Hill et al. (2009) that the mixing process might play an important role when the stratocumulus is thin (∼ 200 m in this study). Furthermore, implementing the new entrainment–mixing parameterization has similar effects on cloud properties to those described by Hoffmann and Feingold (2019), who used the linear eddy model to represent subgrid-scale turbulent mixing. Note that stratocumulus clouds occur in most regions around the world and are important contributors to the surface radiation budget (Wood, 2012; Zheng et al., 2016; Wang et al., 2021; Wang and Feingold, 2009). Stratocumulus clouds dominate in some regions and occur over 60 % of the time as vast long-lived sheets, such as the semipermanent subtropical marine stratocumulus sheets (Wood, 2012). In these regions, a nearly 6 % decrease in τ caused by the new entrainment–mixing parameterization is expected to have significant effects on the simulation of regional radiative properties and climate change.

The averaged influences of the new entrainment–mixing parametrization over all the simulation periods are also examined (Table 4). Quantitatively, the effect of the new entrainment–mixing parameterization is much greater on stratocumulus clouds than on cumulus clouds. Compared to default, new has an average change of −6.20 % in N_c, +2.01 % in r_v, and −3.23 % in τ. When the aerosol concentration increases 10-fold, the differences in N_c, r_v, and τ between default_10 and new_10 are −9.00 %, +3.16 %, and −4.14 %, respectively. These differences are larger than the largest changes in the cumulus case.

Table 4 Cloud water mixing ratio (q_c), cloud droplet number concentration (N_c), cloud droplet volume mean radius (r_v), cloud water path (CWP), and cloud optical depth (τ) in all simulations for the entire lifetime of the stratocumulus case. The experiments are detailed in Table 1.

3.3 Mechanisms of the effects of the new entrainment–mixing parameterization

The different effects of the new entrainment–mixing parameterization on different types of clouds and even on different stages of stratocumulus clouds are likely be related to variations in the dominant mixing mechanism. To confirm this, we calculate the average ψ at all grid points experiencing evaporation, the proportion of inhomogeneous mixing grid points to all grid points experiencing evaporation, and the average ψ at the inhomogeneous mixing grid points in new and new_10 (Table 5) for the cumulus case, as well as the mature and dissipation stages in the stratocumulus case.

Table 5 Homogeneous mixing degree (ψ) at all grid points experiencing evaporation, the proportion of inhomogeneous mixing grid points to all grid points experiencing evaporation, and ψ at the inhomogeneous mixing grid points in the experiments new and new_10 (Table 1) for the cumulus (Cu) and stratocumulus (St) cases. The numbers in and out of the parentheses are the results at the mature and dissipation stages in the stratocumulus (St) case, respectively. The experiments are detailed in Table 1.

For the cumulus case, simulations exhibit large ψ and a small proportion of inhomogeneous mixing, indicating that homogeneous mixing is the dominant entrainment–mixing mechanism (Luo et al., 2020; Lu et al., 2013). Correspondingly, the influences of the new entrainment–mixing parameterization on the cloud physical properties are not significant, as shown in Fig. 3 and Table 2. The new_10 model exhibits a smaller average ψ and a larger proportion of inhomogeneous mixing than new, which results in larger changes in cloud physics, as mentioned in Sect. 3.1.

Figure 6 Probability distribution functions (PDFs) of (a) turbulence dissipation rate (ε), (b) cloud droplet volume mean radius (r_v), (c) evaporation timescale (τ_evap), and (d) transition scale number (N_L) of cloud grids experiencing the entrainment–mixing process in the simulations with the new entrainment–mixing parameterization for the cumulus case (Cu, the solid lines) and the stratocumulus case (St, the dashed lines), respectively. The experiments are detailed in Table 1.

For the stratocumulus case, Table 5 shows the average ψ at all grid points experiencing evaporation, the proportion of inhomogeneous mixing grid points to all grid points experiencing evaporation, and the average ψ at the inhomogeneous mixing grid points during the two stages. The mature stage always has a smaller ψ but a larger proportion of inhomogeneous mixing than the dissipation stage. The inhomogeneous mixing process dominates the mature stage in new because more than 60 % of the grid points experience inhomogeneous mixing. The inhomogeneous mixing process is more dominant in new_10 because less than 3 % of the cloudy grid points experience a homogeneous mixing process during the mature stage, which explains why new_10 has the largest influence when implementing the new entrainment–mixing parametrization. Meanwhile, the average ψ in both stages is smaller than that in the cumulus case for the same simulation configuration. Thus, the effects of the new entrainment–mixing parameterization are more significant for stratocumulus than for cumulus clouds, especially at the mature stage. It is noted that the average ψ and the proportion of inhomogeneous mixing at the dissipation stage of new in the stratocumulus case are very close to the results of new_10 in the cumulus case. This is because the cloud fraction decreases sharply during the dissipation stage; the stratocumulus clouds break up and produce cumulus clouds with small cloud droplet radii.

3.4 The effects of dissipation rate and aerosol concentration on the entrainment–mixing process

Previous studies have shown the notable effects of the dissipation rate and aerosol concentration on the entrainment–mixing process. For example, Luo et al. (2020) changed ε from 10⁻⁵ to 10⁻² m² s⁻³ and noted huge differences in the corresponding ψ. Small et al. (2013) compared aircraft observations with different background concentrations and found that higher-pollution flights tended to slightly more inhomogeneous mixing; Jarecka et al. (2013) also showed various homogeneities of subgrid mixing when aerosol concentration increases 10-fold. To explain the different behaviors of different simulations with the new entrainment–mixing parameterization, the influences of ε and aerosol concentration are examined. Figure 6 shows the probability distribution functions (PDFs) of ε, r_v, τ_evap, and N_L for cloud grids experiencing entrainment–mixing processes in new and new_10 for the cumulus and stratocumulus cases, respectively. The PDFs from the mature and dissipation stages of the stratocumulus case are shown in Fig. 7.

Figure 7 Probability distribution functions (PDFs) of (a) turbulence dissipation rate (ε), (b) cloud droplet volume mean radius (r_v), (c) evaporation timescale (τ_evap), and (d) transition scale number (N_L) of cloud grids experiencing the entrainment–mixing process at the mature stage from 12:00 to 16:00 UTC (the red lines) and the dissipation stage from 21:00 to 24:00 UTC (the green lines) in new for the stratocumulus case. The experiment is detailed in Table 1.

3.4.1 Dissipation rate

According to Eq. (3), N_L is a function of ${\mathit{\epsilon }}^{\mathrm{3}/\mathrm{4}}$; hence, the PDF of ε directly affects N_L and further results in different ψ. For the cumulus case, the mean ε of 0.0043 m² s⁻³ in new is similar to those obtained for cumulus clouds in previous studies (e.g. Lu et al., 2016; Hoffmann et al., 2019). As shown in Fig. 1, cloud grids experience a homogeneous mixing process if N_L is larger than ∼ 10⁵, the limited distribution of N_L values less than 10⁵ in new results in a very small number of cloud grid points undergoing an inhomogeneous mixing process. Even at the cloud grid points that undergo inhomogeneous mixing, the average ψ is large (98.62 %) because most of the N_L values are larger than 10³. Therefore, the cloud properties in new are close to those in default.

For the stratocumulus case, the mean value of ε (2.9 × 10⁻⁴ m² s⁻³ in new) is an order of magnitude less than those in the cumulus case. Therefore, compared with the cumulus case, N_L is reduced in the stratocumulus case, while the peak value of new almost reaches the criterion of inhomogeneous mixing (∼ 10⁵). For the two stages of stratocumulus clouds, ε is an order of magnitude smaller, but r_v was larger (Fig. 7) during the mature stage than during the dissipation stage. According to Eq. (5a), droplets with smaller r_v are more prone to complete evaporation and have a smaller τ_evap. The combination of smaller ε and larger r_v results in a smaller N_L (Eq. 3). This is the reason for the new entrainment–mixing parametrization having more significant effects during the mature stage than during the dissipation stage. In addition, the similarity of the ε and r_v values during the dissipation stage of the stratocumulus case in new, compared to the cumulus case in new_10 (Fig. 6a and b), explains the similar average ψ values of these scenarios and the proportion of inhomogeneous mixing (Table 5).

Therefore, the distribution of ε has a vital impact on the influence of the new entrainment–mixing parameterization. Smaller values of ε result in the new entrainment–mixing parameterization having a more significant influence. Moreover, the r_v in the stratocumulus case is smaller than that in the cumulus case, which is also conducive to a more inhomogeneous mixing process. These are the reasons why the implementation of the new entrainment–mixing parameterization has a larger influence in the stratocumulus case than in the cumulus case when compared to a homogeneous mixing mechanism.

3.4.2 Aerosol concentration

The aerosol concentration affects the entrainment–mixing process by decreasing the cloud droplet radius. As r_v decreases, the distributions of τ_evap in new_10 move to a smaller overall value, while the mean value is an order of magnitude smaller than that in new, which causes a much smaller N_L because N_L is proportional to ${\mathit{\tau }}_{\mathrm{evap}}^{\mathrm{3}/\mathrm{2}}$ (Eqs. 3a and 3b). The larger percentage of smaller N_L values indicates that in new_10, more grid points undergo an inhomogeneous mixing process, and the proportion of such grid points is much larger than in the new model (Table 5). Therefore, compared to new, new_10 exhibits a smaller ψ, and the effects of the new entrainment–mixing parameterization on cloud properties are more obvious for both the cumulus and stratocumulus cases.

3.5 Verification by the simulations with a different parameterization using entrained air relative humidity

In the above simulations, the new entrainment–mixing parameterization is based on the grid mean RH. This section serves to verify these simulations using the entrainment–mixing parameterization proposed by Luo et al. (2020),

$$\begin{array}{}\text{(9)}& \mathit{\psi }=\mathrm{107.96}\exp \left(-\mathrm{0.95}{N}_{\mathrm{L}}^{-\mathrm{0.35}}\right),\end{array}$$

which was developed using the entrained air RH in the EMPM. This parameterization needs the entrained air RH within each grid in WRF, which is estimated following Grabowski (2007) and Jarecka et al. (2009). Briefly, assuming that RH mixes linearly when the dry air entrains into the cloud, then entrained air RH can be simply calculated by

$$\begin{array}{}\text{(10)}& {\mathrm{RH}}_{\mathrm{entrained}}={\displaystyle \frac{{\mathrm{RH}}_{\mathrm{grid}}-(\mathrm{1}-f){\mathrm{RH}}_{\mathrm{cloud}}}{f}},\end{array}$$

where the subscripts entrained, cloud, and grid indicate the RH of the entrained, cloudy, and grid point air, respectively. In Eq. (10), although the cloudy air RH is approximately 100 % and grid mean RH is predicted in the model, the entrained air fraction f needs to be further parameterized. To obtain f at 100 m, a parameterization of f is developed based on the simulations for both the cumulus and stratocumulus cases with a higher resolution of 10 m; the other configurations are the same as those in the experiment default. The 10 m resolution simulation results are then averaged to the resolution of 100 m. Following Xu and Randall (1996), “1−f” can be fitted by the function

$$\begin{array}{}\text{(11)}& \mathrm{1}-f={\mathrm{RH}}_{\mathrm{grid}}^{\mathit{\gamma }}\left[\mathrm{1}-\exp \left(-\mathit{\beta }{q}_{\mathrm{c}}\right)\right],\end{array}$$

where γ and β are empirical parameters. Figure 8 shows that the parameterization with γ=8.72 and β=1.47 × 10⁴ can reproduce well the simulated values of “1−f”, with the correlation coefficient (R) of 0.89 and significant level p value < 0.01. Considering that local shear ($\mathrm{d}w/\mathrm{d}z$) and buoyancy (B) may drive turbulence generation and entrainment for a microscale process, the two quantities are also used to fit “1−f” except for RH_grid and q_c. However, the addition of $\mathrm{d}w/\mathrm{d}z$ and B to Eq. (11) does not increase R. Therefore, using RH_grid and q_c to parametrize “1−f” is good and reasonable for a microscale process.

Figure 8 The fitted 1−f as a function of the calculated 1−f. The fitted 1−f is obtained by the fitting functions with grid mean relative humidity (RH_grid) and cloud water mixing ratio (q_c). The black line denotes the 1 : 1 line.

Figure 9 The temporal evolutions of main cloud microphysical and optical properties in all simulation experiments for the cumulus case, including (a) cloud water mixing ratio (q_c) (g kg⁻¹), (b) cloud droplet number concentration (N_c) (cm⁻³), (c) cloud droplet volume mean radius (r_v) (µm), (d) cloud water path (CWP) (g m⁻²), and (e) cloud optical depth (τ). In (b), N_c values in the experiments default and new_f are normalized by the maximum cloud droplet concentration (N_cmax) from default; N_c values in the experiments default_10 and new_f_10 are normalized by N_cmax from default_10; new_f and new_f_10 are the experiments using entrained air relative humidity.

Figure 10 The temporal evolutions of main cloud microphysical and optical properties in all simulation experiments for the stratocumulus case, including (a) cloud water mixing ratio (q_c) (g kg⁻¹), (b) cloud droplet number concentration (N_c) (cm⁻³), (c) cloud droplet volume mean radius (r_v) (µm), (d) cloud water path (CWP) (g m⁻²), and (e) cloud optical depth (τ). In (b), N_c values in the experiments default and new are normalized by the maximum cloud droplet number concentration (N_cmax) from default; N_c values in the experiments default_10 and new_10 are normalized by N_cmax from default_10; new_f and new_f_10 are the experiments using entrained air relative humidity.

Table 6 Cloud water mixing ratio (q_c), cloud droplet number concentration (N_c), cloud droplet volume mean radius (r_v), cloud water path (CWP), and cloud optical depth (τ) in new_f and new_f_10 for the cumulus (Cu) and stratocumulus (Sc) cases. The results of new and new_10 in Tables 2 and 4 are shown in the parentheses.

Equations (9)–(11) are applied in the simulations for both the cumulus and stratocumulus cases with different aerosol background (hereafter new_f and new_f_10). The same as for Figs. 3 and 5, the temporal evolutions of the cloud physical properties (q_c, N_c, r_v, CWP, and τ) in default, default_10, new_f, and new_f_10 are shown in Figs. 9 and 10. The results are similar to Figs. 3 and 5. The mean values of these properties of new_f and new_f_10 for the cumulus and stratocumulus cases are also shown in Table 6, and the results of new and new_10 are also shown in the parentheses for the convenience of comparison. The results of new_f and new are very similar, with the maximum difference being no more than 1 %, and so are the results of 10-fold aerosol background. Such a close agreement suggests that the results of the new entrainment–mixing parametrization with grid mean RH are reliable.

It is worth noting that instead of parameterizing f, Jarecka et al. (2009, 2013) added an equation to predict f for each grid. In principle, this is a good choice if this method is available in models. Our method shown here is an alternative way to represent the entrainment–mixing process when the prognostic f is not available.

4 Concluding remarks

The entrainment–mixing process near cloud edges has important effects on cloud microphysics, but the most commonly used microphysics schemes simply assume one extreme mechanism, that is, homogeneous entrainment–mixing. This study first improves the entrainment–mixing parameterization proposed by Luo et al. (2020), which connects the homogeneous mixing degree and transition scale number to estimate the homogeneity of the subgrid mixing process and its impact on the droplet number concentration. The improved parameterization uses grid mean relative humidity and can be implemented directly into microphysics schemes; there is no need to know the relative humidity of the entrained air. Second, the modified entrainment–mixing parameterization is implemented in the two-moment Thompson aerosol-aware scheme of the LES version of WRF-Solar to examine its effects on the microphysical and optical properties of cumulus and stratocumulus clouds. Third, several sensitivity experiments are conducted to investigate the effects of the new entrainment–mixing parameterization under different conditions of turbulence dissipation rate and aerosol number concentration. The results of implementing the new entrainment–mixing parameterization are finally verified by the results using entrained air properties.

Unlike the commonly assumed homogeneous mixing scenario, the new entrainment–mixing parameterization produces a smaller cloud droplet number concentration and larger cloud droplet radius, with the degree of difference depending on cloud types and stages. Sensitivity tests show that in the cumulus case, the largest average influence of the new entrainment–mixing parameterization occurs under a high aerosol background but results in only a 2.74 % decrease in cloud droplet number concentration and a 0.57 % increase in cloud droplet volume mean radius. The changes become even smaller with a low aerosol background because of the larger cloud droplet radius. In contrast, the new entrainment–mixing parameterization has a larger influence on the microphysical and optical properties of stratocumulus clouds, especially under a high aerosol background and during the mature stage, with a cloud fraction equal to 1. The largest changes resulting from the new entrainment–mixing parameterization are −9.69 %, +3.88 %, and −5.85 % for cloud number concentration, cloud droplet volume mean radius, and cloud optical depth, respectively. The new entrainment–mixing parameterization has less of an influence on the dissipation stage than on the mature stage of the stratocumulus case but affects this case more than the cumulus case.

The varying effects of the new entrainment–mixing parameterization are caused by variations in the dominant entrainment–mixing mechanism between different cloud types and stages. Compared to the cumulus case, the stratocumulus case has a much smaller homogeneous mixing degree and a larger proportion of inhomogeneous mixing grid points, especially during the mature stage, which indicates that the inhomogeneous mixing mechanism dominates in the stratocumulus case, while the homogeneous mixing mechanism dominates in the cumulus case. As mentioned above, the changes in physical properties of stratocumulus clouds in the dissipation stage are between those in the mature stage and those of the cumulus case; this is because stratocumulus clouds dissipate sharply to form small cumulus clouds, and the degree of homogeneous mixing during the dissipation stage is therefore between that which occurs during the mature stage and the cumulus case.

Sensitivity studies show that turbulence dissipation rate and aerosol concentration can have notable effects on the subgrid homogeneity of the mixing process. A larger dissipation rate can accelerate the mixing process, which results in a larger transition scale number and homogeneous mixing degree and, therefore, a mostly homogenous mixing mechanism. This is why the cumulus case exhibits smaller changes than the stratocumulus case after the new entrainment–mixing parameterization is implemented. Larger aerosol number concentrations cause a smaller cloud droplet radius. Smaller droplets evaporate more easily, which leads to a smaller transition scale number and a smaller homogeneous mixing degree. Thus, the entrainment–mixing mechanism tends to be inhomogeneous. Therefore, a larger aerosol number concentration increases the influence of the new entrainment–mixing parameterization in both the cumulus and stratocumulus cases.

The influences of implementing the new entrainment–mixing parameterization with grid mean relative humidity have been verified by simulations with entrained air properties. The entrained air properties are obtained and calculated from simulations with a finer resolution (10 m). Sensitivity tests show similar cloud microphysical and optical properties in the two different methods, which suggests that the new entrainment–mixing parameterization with grid mean relative humidity is convincing.

Note that the new entrainment–mixing parameterization could be more important in the models if the relative humidity near the cloud is more accurately simulated because numerical diffusion may spuriously humidify the entrained air (Hoffmann and Feingold, 2019). The artificially increased relative humidity limits the influences of the new entrainment–mixing parameterization because homogeneous and inhomogeneous entrainment–mixing processes are close to each other under conditions of high relative humidity.