## Abstract

The role of nonlocal transport on the development and maintenance of marine stratocumulus (Sc) clouds in coarse-resolution models is investigated, with a special emphasis on the downdraft contribution. A new parameterization of cloud-top-triggered downdrafts is proposed and validated against large-eddy simulation (LES) for two Sc cases. The applied nonlocal mass-flux scheme is part of the stochastic multiplume eddy-diffusivity/mass-flux (EDMF) framework decomposing the turbulent transport into local and nonlocal contributions. The complementary local turbulent transport is represented with the Mellor–Yamada–Nakanishi–Niino (MYNN) scheme. This EDMF version has been implemented in the Weather Research and Forecasting (WRF) single-column model (SCM) and tested for three model versions: without mass flux, with updrafts only, and with both updrafts and downdrafts. In the LES, the downdraft and updraft contributions to the total heat and moisture transport are comparable and significant. The WRF SCM results show a good agreement between the parameterized downdraft turbulent transport and LES. While including updrafts greatly improves the modeling of Sc clouds over the simulation without mass flux, the addition of downdrafts is less significant, although it helps improve the moisture profile in the planetary boundary layer.

## 1. Introduction

Stratocumulus (Sc) clouds are one of the most common cloud types on Earth (Hahn and Warren 2007). They form under strong temperature inversions and are prevalent off the western coast of continents, on the descending side of the Hadley cell. Their impact on Earth’s energy budget is significant as they strongly reflect incoming solar radiation, with a much weaker effect on outgoing longwave radiation (Wood 2012). Accurate modeling of Sc clouds has high importance for several reasons: (i) they are one of the key sources of uncertainty in climate predictions (Bony and Dufresne 2005; Zelinka et al. 2017), (ii) they affect solar power integration into the electric grid (Yang and Kleissl 2016; Zhong et al. 2017; Wu et al. 2018), and (iii) they impact aviation by hindering the takeoff and landing of flights (Reynolds et al. 2012).

Physical processes governing the evolution of the stratocumulus-topped boundary layer (STBL)—such as cloud-top radiative cooling, entrainment, evaporative cooling, surface fluxes, wind shear, and precipitation—widely range on spatial and temporal scales, and modeling Sc clouds is quite challenging as a result (e.g., Lilly 1968; Stevens 2002; Wood 2012). Efforts through both observational campaigns (e.g., Stevens et al. 2003; Malinowski et al. 2013; Crosbie et al. 2016) and high-resolution numerical modeling (e.g., Stevens et al. 2005; Kurowski et al. 2009; Yamaguchi and Randall 2012; Chung et al. 2012; Blossey et al. 2013; de Lozar and Mellado 2015; Pedersen et al. 2016; Mellado et al. 2018; Matheou and Teixeira 2019) have significantly advanced our understanding of the physics of Sc clouds. These physical insights are important for numerical weather prediction (NWP) and general circulation models (GCMs) where grid resolution is coarse.

The picture emerging from those studies is that cloud-top radiative cooling is a critical source of STBL turbulence (Matheou and Teixeira 2019), contributing to cloud-top entrainment (Mellado 2017). The combined effect of both evaporative and radiative cooling—the former typically enhanced by wind shear (Mellado et al. 2014)—destabilizes the top of cloud layer through buoyancy reversal that leads to the formation of negatively buoyant weak downdrafts. This process is often considered responsible for the generation of cloud holes in largely unbroken Sc clouds (Gerber et al. 2005; Kurowski et al. 2009). Many small-scale phenomena (e.g., entrainment, shear, evaporative cooling, cloud microphysics) are at play in the origin of downdrafts and can strongly influence vertical mixing (Mellado 2017). Exactly how these processes interact with each other remains a research challenge.

Turbulent transport in the STBL is the main driver to the formation, maintenance, and dissipation of Sc clouds. In coarse-resolution models, turbulent transport is typically parameterized using simplified one-dimensional planetary boundary layer (PBL) schemes. Global NWP models (e.g., Teixeira 1999) and climate models tend to underestimate Sc clouds (Teixeira et al. 2011; Lin et al. 2014), although there is an improvement in the representation of the radiative properties by a newer generation of climate models (Engström et al. 2014). In terms of mesoscale models, Ghonima et al. (2017) compared three different PBL schemes in the Weather Research and Forecasting (WRF) Model and found that they all underestimate entrainment, producing too moist and cold STBLs. Huang et al. (2013) compared five different WRF PBL parameterizations and highlighted the difficulties of simulating the STBL. Recent studies supported the importance of downdrafts in transporting turbulent heat and moisture flux in the PBL (Chinita et al. 2018; Davini et al. 2017; Brient et al. 2019) through analyzing LES of STBL. Brient et al. (2019) concluded that for a more accurate parameterization of turbulence within STBL, downdrafts should be explicitly included in climate models. Downdrafts were recently implemented by Han and Bretherton (2019) in a turbulent kinetic energy (TKE)-based moist eddy-diffusivity/mass-flux (EDMF) parameterization within the GFS model, and they found more accurate liquid water and wind speed profiles for marine STBLs.

This study introduces parameterized downdrafts into NWP and aims at investigating their impact on the evolution of the STBL. To test whether convective downdrafts are necessary to properly represent Sc clouds, we implement a new downdraft parameterization in WRF based on the EDMF approach that uses Mellor–Yamada–Nakanishi–Niino (MYNN) as the ED component. This differs from Han and Bretherton (2019) where different ED and MF models were used and additional features were implemented to advance the vertical turbulence mixing parameterization for not only STBL but also other conditions. We place a special emphasis on evaluating the role of nonlocal transport in STBL with gradual changes to the model in order to separate the effects coming from convective downdrafts. The new parameterization is evaluated in two typical STBL cases, frequently used in modeling studies.

Section 2 describes the EDMF and MYNN schemes as well as the updraft and downdraft implementation in WRF. The numerical design of the LES setup, WRF single-column model (SCM), and updraft and downdraft properties are described in section 3. WRF SCM results for both STBL cases are shown in section 4. Finally, conclusions are presented in section 5.

## 2. PBL scheme with downdrafts

In coarse-resolution atmospheric models, the PBL scheme determines turbulent flux profiles within the PBL as well as the overlying air, providing tendencies of temperature, moisture, and horizontal momentum due to mixing and turbulent transport for the entire atmospheric column. This section first gives an overview of the EDMF framework, then the details of ED and MF models are presented (sections 2b and 2c). The properties of updrafts and downdrafts are diagnosed using LES and presented in section 3 in order to quantify the validity of the parameterized mass-flux model.

### a. The eddy-diffusivity/mass-flux (EDMF) approach

Siebesma and Teixeira (2000), Teixeira and Siebesma (2000), and Siebesma et al. (2007) introduced the eddy diffusivity/mass-flux (EDMF) approach for parameterizing turbulence in the dry convective boundary layer, and additional improvements have been made by Witek et al. (2011). The idea behind EDMF is to parameterize the turbulent fluxes as a sum of local (diffusive) transport through ED and nonlocal (convective) transport through the mass-flux contribution. The EDMF approach has been extended to represent moist convection since then (e.g., Soares et al. 2004; Neggers et al. 2009; Neggers 2009; Angevine et al. 2010, 2018; Sušelj et al. 2013; Suselj et al. 2019a,b). For the moist extention, the updrafts start out as dry and condense when the saturation conditions are met. In other words, dry updrafts can continuously change into moist updrafts, without assuming any coupling between those two layers.

### b. ED scheme: The Mellor–Yamada–Nakanishi–Niino (MYNN)

The ED component we use is the level 2.5 Mellor–Yamada–Nakanishi–Niino (MYNN) model, which is a modified Mellor–Yamada turbulence closure scheme originally developed by Mellor and Yamada (1982), with significant improvements made over the years (Nakanishi and Niino 2006, 2009). In MYNN, vertical turbulent fluxes are modeled according to *K*-theory:

where eddy diffusivity *K* is parameterized as a function of the TKE (*q*), master length scale *L*, and stability correction functions *S*_{h,m}, which differ for heat and momentum:

The prognostic thermodynamic equations in MYNN use moist conserved variables: liquid water potential temperature *θ*_{l} and total water mixing ratio *q*_{t}. The prognostic dynamic variables are the horizontal components of wind *u* and *υ*. An additional prognostic equation of the MYNN Level 2.5 model solves the (doubled) subgrid TKE: $q2=2\xd7TKE=u\u20322\xaf+\upsilon \u20322\xaf+w\u20322\xaf$, and is formulated as

Equation (3) describes the tendency of TKE, due to turbulent and pressure transport, shear production, buoyant production, and turbulent dissipation. *L* is the master length scale as in Eq. (2), and *S*_{q} = 3*S*_{m} is the stability correction function for TKE [see Nakanishi and Niino (2009) for detailed formulations]; *L* is designed such that the smallest length scale out of three different formulations dominates at a given level. The first formulation is the surface length scale *L*_{sfc}, which is the Prandtl mixing length corrected for stability. It is small near the surface, but increases rapidly with height. The second one, the turbulent length scale for a well-mixed layer *L*_{turb}, is formulated as a function of the vertically integrated TKE, independent of height. Finally, the buoyancy length scale *L*_{buoy} is computed as a function of local stratification (i.e., ∂*θ*_{υ}/∂*z*), and it decreases with increasing stratification. The buoyancy length scale is only active in stable conditions. The stability functions for heat and moisture *S*_{h,m} contain empirical constants, which generally decrease with increasing stability, as they are inversely related to the Richardson number [Eqs. (27) and (28) in Nakanishi and Niino 2009]. Finally, the dissipation rate is parameterized as *ε* = *q*^{3}/*B*_{1}*L*, where *B*_{1} is a closure constant (*B*_{1} = 24 in the MYNN scheme).

### c. Adding mass flux to MYNN

The MYNN Level 2.5 ED model determines turbulent mixing at each vertical level based on the gradients in scalars between immediately adjacent vertical levels [Eq. (1)]. When deep mixing due to larger eddies becomes important, the MYNN scheme has been shown to produce erroneous thermodynamic profiles (Huang et al. 2013). Nonlocal models, such as the YSU and ACM2 schemes, account for this deep mixing by using a countergradient term (Hong et al. 2006) or a transilient mass-flux matrix (Pleim 2007). Another common approach is the EDMF framework, which decomposes the subgrid vertical mixing into local mixing through ED and nonlocal [mass flux (MF)] transport through convective plumes. Traditionally, PBL schemes, such as MYNN, model the turbulence within the PBL through only turbulent diffusion. In the EDMF framework, ED is used to model the turbulent transport in the nonconvective environment, with an additional contribution from the convective mass flux.

#### 1) Mass-flux model overview

To represent nonlocal transport, we use the stochastic multiplume EDMF model. The idea behind this model is that the horizontal subgrid domain is composed of an ensemble of convective plumes and the remaining nonconvective environment. The multiplume approach is designed to account for the nonlinear interactions between the plumes and the environment through the stochastic lateral entrainment. Following the same notation as Suselj et al. (2019a,b), the grid-mean value of any variable *φ* can be written as

where *N*/*M* is the total number of updrafts/downdrafts. The subscripts *u*_{n}, *d*_{m}, and *e* denote mean values from the *n*th updraft, *m*th downdraft, and the environment, while $aun$, $adm$, and *a*_{e} are the corresponding areas. In WRF, assuming the fractional area of updraft and downdraft are small, we approximate $\phi \xaf\u2248\phi e$, and the turbulent flux can be written as [see Eq. (7) in Suselj et al. (2019b)]:

where the vertical transport of nonconvective environment $w\u2032\phi \u2032\xaf|e$ is modeled using Eq. (1).

#### 2) Surface-driven updrafts

A version of EDMF including surface-driven updrafts (Olson et al. 2019) has been implemented as an add-on option in MYNN since WRF v3.8 and is used for NOAA’s operational Rapid Refresh (RAP; Benjamin et al. 2016) and High Resolution Rapid Refresh (HRRR) forecast systems. The original version of this dynamic multiplume mass-flux scheme in WRF v3.8 (bl_mynn_edmf = 1) followed Sušelj et al. (2013), but the version in the current WRF v4.0 contains considerable changes from the original form. We do not base our EDMF implementation (bl_mynn_edmf = 3) on what is currently available in WRF, but instead follow Sušelj et al. (2013) and Suselj et al. (2019a,b). The numerical implementation is documented in Suselj et al. (2019b) (their appendix B).

The surface-driven updrafts are represented by an ensemble of steady-state plumes with different initial conditions and stochastic entrainment rate profiles. The thermodynamic and dynamic properties of the *n*th updraft $\phi un={\theta l,un,qt,un,uun,\upsilon un}$ follow:

where $\epsilon un$ is the entrainment rate. Note that an additional source term, due to microphysical processes in Suselj et al. (2019b), is not included here as it has no effect in nonprecipitating STBL. The number of updrafts is fixed to ten (*n* = 1, …, *N*; *N* = 10). The steady-state equation of the updraft velocity is

where *a*_{w} = 1, *b*_{w} = 1.5 are model constants (de Roode et al. 2012; Sušelj et al. 2013; Suselj et al. 2019b). Variable $Bun=g\u2061(\theta \upsilon ,un/\theta \xaf\upsilon \u22121)$ is the updraft buoyancy, and $\theta \upsilon =\theta \u2061(1+0.61q\upsilon \u2212ql)$ is the virtual potential temperature. $Pwud$ represents the dynamical pressure effects as updrafts approach the inversion and is parameterized as

where *z*_{00} denotes the distance from *z*_{i} when $pwud$ starts to be in effect. For this work, we use *z*_{00} = 100 m. Assuming a normal distribution of the vertical velocity near the surface, the updrafts are thought to represent the positive tail of the distribution, between one and three standard deviations, divided into *N* bins. This results in a total updraft area of approximately 15% near the surface. The thermodynamic surface conditions for the updrafts are identical to Suselj et al. (2019a) (their appendix A); $\epsilon un$ is the stochastic entrainment rate, computed as

where *ε*_{0} = 0.2 is the fractional mass of air entrained in a single entrainment event; *P*(*λ*) is a random number drawn from the Poisson distribution with parameter *λ* = (Δ*z*/*L*_{ε}), which represents the number of entrainment events a single updraft experiences over height Δ*z*. *L*_{0} = 100 m denotes the distance a plume needs to travel to entrain once. The exponential term in Eq. (10) represents the dynamic effect near strong temperature inversion, as the updrafts cannot penetrate above that layer and are assumed to entrain more and disintegrate; *c*_{ent} = 0.5 is a model constant controlling how fast entrainment length decreases with height. For STBL, we use the cloud-top height *z*_{i} (also known as the inversion height) to denote where this dynamic effect is at its strongest. *z*_{i} is defined as the last point near the PBL height where *q*_{l} > 10^{−6} kg kg^{−1}, and cloud fraction is greater than 50%. This definition of *z*_{i} is identical to that in Olson et al. (2019), where they included an option for top-down buoyancy production in ED when Sc clouds were present. In the MYNN parameterization, there are three options to represent subgrid cloudiness, which are controlled by bl_mynn_cloudpdf parameter. In this work, bl_mynn_cloudpdf = 1, for which a statistical partial condensation cloud scheme based on joint-Gaussian probability distribution function of *θ*_{l} and *q*_{t} is used (Kuwano-Yoshida et al. 2010). By default, the Gaussian PDFs are applied to the whole grid box (i.e., including nonconvective environment and convective updrafts and downdrafts). We thus assume that Gaussian distributions of the thermodynamic variables (cf. Fig. 1) yield reasonably accurate cloud cover and liquid water values for STBL. Cloud fraction would ideally be computed from Eq. (4), and we use this approximation for simplicity. Note that for STBL, saturation conditions are usually met for most of the PDFs area.

While Sušelj et al. (2013) did not include either the dynamical pressure effect [i.e., $Pwud$ term in Eq. (7)] or modification of entrainment length (*L*_{ε}) by proximity of inversion for the STBL simulation, we find that those modifications yield results that are more consistent with the plume statistics in LES, as discussed further in section 3. The entrainment rate is the same for all variables ($\theta l,un$, $qt,un$, $uun$, and $\upsilon un$). Although Suselj et al. (2019b) used $\u2061(1/3)\epsilon un$ for $uun$ and $\upsilon un$, we find that equal entrainment rate results in more consistent *u* and *υ* profiles.

Since each updraft is characterized by different surface conditions and entrainment rates, the thermodynamic properties and termination heights also differ. Each plume is integrated independently in the vertical until the vertical velocity becomes negative. As mentioned above, we assume Gaussian distribution of the subgrid-scale *θ*_{l} and *q*_{t}, which determines the grid-mean liquid water and cloud cover. We then compute the condensation in the convective plumes assuming an uniform distribution of *θ*_{l} and *q*_{t} within each of the plumes. Therefore, the condensation within a plume occurs when *q*_{t} exceeds the saturated mixing ratio (*q*_{s}). The precise calculation of condensation within each plume is necessary because of its impact on plume buoyancy. Finally, there exist dry and partly moist plumes among the *N* updrafts, and the fate of each plume is determined by its initial conditions, dynamical pressure effect, and lateral entrainment. Since each individual updraft is integrated independently, whenever the vertical velocity becomes negative and the updraft terminates, the total updraft area is reduced. This is common in regions with strong lateral entrainment rates.

#### 3) Cloud-top-triggered downdrafts

Several important physical processes are at play near the STBL top. Radiative and evaporative cooling produces cooled downdrafts and drives buoyant production of turbulence in the PBL. Entrainment from the free troposphere can impact downdrafts near the cloud top: warm air from the free troposphere counteracts the radiative cooling and buoyant production of turbulence. When the PBL is less turbulent, the entrainment rate decreases, indicating a negative feedback loop (Wood 2012). Surface-driven updrafts may also affect the downdrafts. As updrafts approach the inversion, they begin to diverge and can help initiate or enhance downdrafts (Kurowski et al. 2009; Davini et al. 2017). This enhances the downdraft vertical velocity and, in turn, the turbulence in the PBL. In the proposed 1D parameterization of downdrafts, those dependencies are important for the formulation of the downdraft initial conditions. Our downdraft parameterization in MYNN can be activated by specifying bl_mynn_edmf_dd = 1 in the namelist. The numerical implementation follows Suselj et al. (2019b) (see their appendix C).

Similar to the surface-driven updrafts, downdrafts are also represented by an ensemble of steady-state plumes with stochastic lateral entrainment. The thermodynamic and dynamic properties of the *m*th downdraft $\phi dm={\theta l,dm,qt,dm,udm,\upsilon dm}$ follow:

$\epsilon dm=(\epsilon 0/\Delta z)P\u2061(\Delta z/L\epsilon )$ is the entrainment rate similar to Eq. (9), where *L*_{ε} = *L*_{0}, and the values of *L*_{0} and *ε*_{0} are the same as for the updrafts. The entrainment rate is same for *θ*_{l,dm} and *q*_{t,dm}, however, it is increased 1.4 times for *u*_{dm} and *υ*_{dm}. We find that increasing entrainment rate for momentum results in more consistent *u* and *υ* profiles. The additional source term due to microphysical processes described in Suselj et al. (2019b) is neglected here. The number of downdrafts is fixed to ten (*m* = 1, …, *M*; *M* = 10). The steady-state equation of the downdraft velocity is identical to Suselj et al. (2019b):

where $pwdd$ represents the dynamical pressure effects as downdrafts approach the surface and is parameterized as

where *z*_{00} = 100 m. This is equivalent to the dynamical pressure effect in updraft, except we replace *z*_{i} with 0.

We assume that downdrafts are initiated randomly in the upper half of the cloud layer. We avoid starting all downdrafts at *z*_{i} to avoid numerical instabilities in this region during the model spinup time (not shown; see next section for details). Similarly to the updraft parameterization, we assume that the downdrafts represent the negative tail of the vertical velocity distribution that is assumed to be normal (between negative one and three standard deviations), resulting in a total downdraft area of approximately 15% slightly below cloud top. The formulation of cloud-top conditions for downdrafts is similar to the formulation for surface-driven updrafts (Suselj et al. 2019a). The difference lies in the parameterization of the variances of vertical velocity *σ*_{w}, total water mixing ratio $\sigma qt$, and virtual potential temperature $\sigma \theta \upsilon $. The strength of downdraft vertical velocity is proportional to *σ*_{w}:

where *c*_{1} = 0.3 is a model constant, and *w*_{*,dd} is the the convective vertical velocity scale, which takes into account both the intensity of surface-driven updrafts and cloud-top radiative cooling and is similar to the entrainment parameterization in Ghonima et al. (2017):

where $w*\u2261(g/\theta \upsilon )w\u2032\theta \upsilon \u2032\xaf|sztop$ is the Deardorff convective velocity scale, $u*$ is the surface friction velocity, and $wrad\u2261(g/\theta \upsilon )w\u2032\theta \upsilon \u2032\xaf|radztop$ is a velocity scale based on the net radiative flux divergence at the cloud top where $w\u2032\theta \upsilon \u2032\xaf|rad=Frad/\rho cp$ (Lock and Macvean 1999). In WRF, *F*_{rad} is defined as the radiative flux divergence between cloud top and cloud base.

The framework for parameterizing $\sigma qt$ and $\sigma \theta \upsilon $ is similar to that described by Köhler (2006). The downdraft initial total mixing ratio deficit is proportional to $\sigma qt$:

where *c*_{2} = 30 is a model constant, and $q*\u2261w\u2032qt\u2032\xafent/wrad$ is the moisture scale due to mixing with entrained air. The entrainment fluxes $w\u2032\phi \u2032\xafent$ are modeled according to the flux-jump relation $w\u2032\phi \u2032\xafent=we\Delta \phi zinv$ (Lilly 1968), where $\Delta \phi zinv=\phi zinv+1\u2212\phi zinv$ represents the jump value of the scalar *φ* across the inversion; *w*_{e} is the entrainment velocity and is parameterized following Ghonima et al. (2017):

In WRF, the jump in moisture, Δ*q*_{t}, is defined as the difference in *q*_{t} at 700 hPa and the surface.

The downdraft initial virtual potential temperature is proportional to $\sigma \theta \upsilon $:

where *c*_{3} = 1 is a model constant, and $\theta \upsilon ,*\u2261w\u2032\theta \upsilon \u2032\xaf\u2009ent/w*,rad$ is the buoyancy scale due to mixing with entrained air and radiative cooling. The jump in heat, Δ*θ*_{υ}, is similar to (Wood and Bretherton 2006):

where *θ*_{υ,700} is *θ*_{υ} at *p* = 700 hPa, *θ*_{υ,0} is *θ*_{υ} at the surface, Γ_{FT} is the free tropospheric adiabat, and *z*_{700} is the height of the *p* = 700-hPa surface. Since difference in *θ*_{υ} at 700 hPa and the surface is a combination of temperature increase across the capping inversion and the accumulated static stability between this inversion and the 700-hPa reference level, we subtract Γ_{FT}(*z*_{700} − *z*_{inv}) to focus on temperature jump across the inversion. We find this definition of inversion jumps to be more systematic and consistent than attempting to diagnose the exact point where the temperature inversion begins and ends.

Similar to the updrafts, equations for each downdraft are independently integrated in the vertical until the downdraft velocity vanishes. Condensation occurs within a downdraft if its total water mixing ratio exceeds the saturated water mixing ratio. Similarly to updrafts, there can exist dry and partly moist plumes among the *M* downdrafts, and the fate of each plume is determined by its initial conditions, dynamical pressure effect near the surface, and lateral entrainment. We assume that both updrafts and downdrafts interact with the mean field only and not with other updrafts or downdrafts. This assumption has been tested in Kurowski et al. (2019) for updrafts in shallow convection, and further investigation for downdrafts should be explored. Since each individual downdraft is integrated independently, whenever vertical velocity becomes zero/positive and the downdraft terminates, the total downdraft area is reduced.

## 3. Design of numerical experiments

### a. LES setup

Large eddy simulations are performed using the UCLA-LES model (Stevens 2010) and treated as “ground truth.” Two idealized nondrizzling marine Sc cases are chosen as baseline simulations: the DYCOMS-II RF01 (Stevens et al. 2005) and CGILS S12 Control (Blossey et al. 2013) (hereafter DYCOMS and CGILS). The experiments are set up following the respective intercomparison studies. Interactive radiation is treated differently in the two cases. Specifically, a simplified model of radiative forcing matching the *δ*-four stream transfer code (Stevens et al. 2005) is used in DYCOMS. As for CGILS, a full radiative transfer code is used, which utilizes Monte Carlo sampling of the spectral integration (Pincus and Stevens 2009). The DYCOMS case is run for 4 h, and the CGILS case is run for 24 h. While we focus our analysis of the updraft and downdraft properties on nocturnal quasi-steady conditions (first 4 h), the 24-h simulation of CGILS provides reference to the generalization of the parameterization during the day. In both experiments, a nonuniform vertically stretched grid is used with 5-m resolution around the inversion, and a several times coarser resolution in the horizontal. This LES setup is identical to that in Ghonima et al. (2017). A summary of the model setups is provided in Table 1.

#### Determining plume properties

Simulation outputs are stored at 1 min intervals from hours 3 to 4 in order to diagnose updraft and downdraft properties. The statistics are averaged over 1 h. We use the joint normal probability density function (PDF) between vertical velocity *w*, total water mixing ratio (*q*_{t} = *q*_{υ} + *q*_{l}), virtual potential temperature [*θ*_{υ} = *θ*(1 + 0.61*q*_{υ} − *q*_{l})], and liquid water potential temperature [*θ*_{l} = *θ* − (*L*_{υ}*q*_{l})(*c*_{p}*π*)^{−1}] to define LES updrafts and downdrafts. *L*_{υ} is the latent heat of vaporization, *c*_{pd} is the specific heat of dry air at constant pressure, *π* is the Exner function, and subscripts are *υ* for vapor, *l* for liquid. We define the normalized variable to be $\phi \u2032=\phi \u2212\phi \xaf/\sigma \phi $, where $\phi \xaf$ is the slab mean and *σ*_{φ} is the standard deviation of *φ*. By carefully investigating the joint PDFs, we define updrafts to be the LES grid points that conform to the following conditions: *w*′ > 1, $qt\u2032>0$, and either $\theta l\u2032<0$ or $\theta \upsilon \u2032>0$. We define downdrafts to be *w*′ < 1, $qt\u2032<\u22121$, and $\theta l\u2032>0$. Specifically, this definition of downdrafts captures the negative tail in the joint normal PDF. Figure 1 shows the joint normal PDF for DYCOMS at a normalized height close to the cloud top (*z*/*z*_{i} = 0.97). A strong negative tail is observed in Fig. 1a, where *w*′ < 0 and $qt\u2032<\u22121$. We also confirm that grid points satisfying these criteria correspond well with negatively buoyant ($\theta \upsilon \u2032<0$) parcels that are warmer in terms of the liquid water potential temperature ($\theta l\u2032>0$). While the definitions of updraft and downdraft used here are not as rigorous as in Chinita et al. (2018), Davini et al. (2017), and Brient et al. (2019), we find that the overall properties are consistent with their study.

The mean downdraft and updraft properties are shown in Fig. 2 for DYCOMS and Fig. 3 for CGILS. Updraft and downdraft areas are comparable in the middle of the PBL (Figs. 2a and 3a), with updrafts decreasing near cloud top and downdrafts decreasing before reaching the surface. Figures 2b,c and 3b,c show partial contributions to the total heat and moisture fluxes from the environment, updrafts, and downdrafts. Similar results are found in both STBL cases: cloud-top entrainment heat flux is largely from updrafts; the peak in downdraft heat and moisture transport is slightly below the peak in updrafts (≈100 m lower); heat and moisture transport from downdrafts is stronger than updrafts in cloudy region; environmental mean of *w*, *θ*_{l}, *θ*_{υ}, *q*_{t}, and *q*_{l} is very close to the grid mean. Both cases have similar updraft and downdraft properties: downdrafts terminate before reaching the surface (Figs. 2a and 3a); updraft and downdraft vertical velocity are approximately a mirror image of each other (Figs. 2d and 3d); downdrafts become negatively buoyant ($\theta \upsilon \u2032<0$) slightly below cloud top (Figs. 2f and 3f); updrafts correspond to thicker cloud regions and downdrafts are collocated with cloud holes (Figs. 2h and 3h). Since the peak in downdraft heat and moisture transport is slightly below the peak in updraft, the choice of starting downdrafts randomly between cloud top and halfway through cloud base is consistent with the findings in LES.

The properties shown in these two STBL cases compare well to the case in Brient et al. (2019), where the First ISCCP Regional Experiment (FIRE) study was simulated for 24 h to study the diurnal cycle of coherent updraft and downdraft properties. Specifically, the nighttime results of Brient et al. (2019) show that the areas of updrafts and downdrafts are comparable in the middle of the PBL (around 12%) and the downdraft area decreases quickly to zero below 100 m, which corresponds well with our findings for DYCOMS. CGILS results show a slightly smaller downdraft area in the middle of the PBL (around 9%). The turbulent heat flux in Brient et al. (2019) shows that the transport of heat by updrafts is the strongest at cloud top, the peak of the downdraft heat transport is slightly below that for the updrafts (≈50 m lower), and the heat transport by updrafts in cloudy region is nearly zero when downdrafts dominate. This corresponds well with DYCOMS, while updrafts in CGILS have a slightly positive heat transport in the cloudy region. As for the turbulent moisture flux, Brient et al. (2019) shows that updrafts dominate from the surface up to slightly above cloud base, while downdrafts dominate in the cloud layer. Moisture flux is similar in DYCOMS and CGILS, but our results show a positive peak of updraft moisture flux near cloud top, making the updraft contribution to the moisture flux a dominating term around cloud top. Chinita et al. (2018) shows large differences in the contribution of updrafts and downdrafts to total flux for DYCOMS in the cloud layer. In general, they find that updrafts account for most of the organized motions near the surface, while downdrafts are more important near the boundary layer top. While the overall properties are similar, updraft and downdraft areas in Chinita et al. (2018) are 5%–10% larger.

### b. WRF single-column model

DYCOMS and CGILS case are simulated using the Weather Research and Forecasting (WRF) v4.0 single-column model (SCM) and compared against LES. Initial conditions and forcing are identical to that in LES (i.e., fixed surface fluxes for DYCOMS and CGILS, large-scale subsidence as in Table 1) and was used previously in Ghonima et al. (2017). The SCM vertical domain includes 116 levels to resolve the lowest 12 km of the troposphere, which comes out to be Δ*z* ≈ 20 m in the first 1 km. A simulation time step of 40 s is used. In section 4c, we show that results are insensitive to time step between 5 and 60 s. Three different versions of one PBL scheme are used to determine the importance of the introduced changes: 1) the original Mellor–Yamada–Nakanishi–Niino scheme (MYNN; hereafter ED) (Nakanishi and Niino 2006, 2009), 2) MYNN with updrafts (EDMF_{U}), and 3) MYNN with updrafts and downdrafts (EDMF_{UD}). For EDMF_{U} and EDMF_{UD}, the MYNN scheme is used as a parameterization of local transport in the nonconvective environment. The radiation scheme is RRTMG (Iacono et al. 2008). No microphysics or cumulus schemes are used since both cases represent nonprecipitating STBL.

## 4. Results

### a. DYCOMS-II RF01

Figure 4 shows the mean fields of *θ*_{l}, *q*_{t}, *q*_{l}, cloud fraction, *u*, *υ*, heat flux ($\rho cpw\u2032\theta l\u2032\xaf$), and moisture flux ($\rho L\upsilon w\u2032qt\u2032\xaf$). Figure 5 shows the time series of liquid water path (LWP), boundary layer averaged heat (*θ*_{l}), and moisture (*q*_{t}) for the three tested PBL schemes and LES. ED has a cold and moist bias in the PBL (Figs. 5b,c), resulting in an overestimation of LWP for the entire simulation. The underestimation of entrainment flux is likely the cause of this behavior as ED fails to model heat and moisture transport between the free-troposphere and the PBL (Figs. 4g,h). Moreover, ED does not have a transition in horizontal wind between the PBL and the free troposphere, indicating that ED does not capture the momentum transport properly (Figs. 4e,f). EDMF_{U} has a weaker cold and moist bias, and the bias in LWP is minimal during hours 3 to 4. However, inversion base height is slightly lower than ED. This is a result of updrafts overshooting into the free troposphere in the early time of the simulation, mixing out the initial inversion base height. EDMF_{UD} has a much smaller bias in boundary layer averaged heat and moisture and has a more well-mixed profile in *q*_{t} than EDMF_{U}. Inversion base height is also slightly lower in EDMF_{UD}. Both EDMF_{U} and EDMF_{UD} capture the entrainment heat and moisture flux well. Among the three tested PBL schemes, EDMF_{U} has the best match in horizontal wind in the PBL, and EDMF_{UD} overestimates *u* but underestimates *υ* in the PBL.

Figures 6 and 7 show the vertical flux contribution from the individual components: environment (ED), updraft, and downdraft. Figure 6 is for EDMF_{U}, which includes only ED and updraft. Note that LES transport in 6A and D includes LES environmental and downdraft transport because in the case of updrafts only, the remaining area is considered to be the environment and should therefore be modeled by ED. Updraft contribution to the heat flux matches the profile in LES well, however it is overestimated in most of PBL and the cloud-top entrainment heat flux is too strong. It is important to note that cloud-top entrainment is not fully understood even in LES. We find here that even though entrainment heat flux appears to be strong, boundary layer averaged temperature in EDMF_{U} is still too cold compared to LES (Fig. 5b). However, EDMF_{U} produces a warmer boundary layer compared to ED, which strongly underestimates entrainment heat flux. Updraft contribution to the moisture flux is overestimated throughout the PBL, but ED component is underestimated and the total moisture flux matches LES well. The initial updraft starting *θ*_{l} and *q*_{t} are stronger than LES (not shown) and eventually leads to overestimation of moisture flux. This indicates that the formulation of updraft surface condition in STBL may be different from shallow convection since we retain the same updraft starting condition used in Suselj et al. (2019a). In shallow convection, surface fluxes are the main driver for updraft surface conditions. Whether other physical processes are at play in the parameterization of updraft surface conditions in STBL should be investigated in the future. We find that in the current configuration, ED compensates for the overestimation of updraft moisture flux, resulting in a good match with LES in the total moisture flux.

Based on 800 additional simulations, exploring the parameter space, with different lateral entrainment rates and dynamical effects (varying *L*_{0} and *c*_{ent} in Eq. (10) from and 10 to 100 m 0.5 to 5 m^{−1}, as well as varying *z*_{00} in Eq. (8) from 50 to 200 m; not shown), we observe that the most important impact of the updraft is the transport near cloud top because ED models an insufficient heat and moisture transport in this location, causing a cold and moist bias. Additionally, ED does not accurately represent a well-mixed layer, while EDMF_{U} has a better well-mixed profile in both *θ*_{l} and *q*_{t}. The final configuration was chosen to have the best match in the mean field of *θ*_{l}, *q*_{t}, and total heat and moisture transport with LES.

For EDMF_{UD}, Fig. 7 shows partial contributions to the total transport from ED, updrafts, and downdrafts. Comparing Figs. 6 and 7, we argue that the downdraft transport is implicitly included in the ED contribution in EDMF_{U} (Figs. 6a,d) as the sum of heat and moisture transport for EDMF_{U} versus EDMF_{UD} is similar. Averaged plume properties from EDMF_{UD} are shown in Fig. 8. For downdraft contribution to total fluxes, EDMF_{UD} underestimates the strength in heat and moisture flux. More spefically, downdraft heat transport decreases too quickly before reaching the surface (Fig. 7c). For moisture tansport, downdraft *q*_{t} also decrease quickly, and the starting downdraft *q*_{t} is underestimated (Fig. 8c). Updraft contribution to heat transport (Fig. 7b) is similar to that in EDMF_{U}, and they both slightly overestimate compared to LES in terms. This can be seen in the overestimation of updraft area and vertical velocity (Figs. 8a,b), and is a result of the positive bias in updraft starting surface conditions, espicially updraft starting vertical velocity. For updraft moisture transport, updrafts in EDMF_{UD} do not overestimate as strongly as EDMF_{U}. This is likely due to downdrafts transporting dry and warm air in the PBL and causing updrafts to mix differently. On top of that, the mean fields of *θ*_{l} and *q*_{t} are different in EDMF_{U} and EDMF_{UD}. Note that since the definition of updrafts and downdrafts in LES is somewhat arbitrary, the total transport should be the main indicator of success for a parameterization. Nevertheless, the definition of updrafts and downdrafts as in section 3 is a reference point for bench-marking updraft and downdraft parameterizations. Overall, general agreement of plume properties are found between the SCMs and LES. For DYCOMS, downdraft transport decreases too quickly for both heat and moisture. We find that modeling downdraft transport in the upper part of the boundary layer correctly is more important than retaining downdraft throughout the PBL. The mean fields respond more to changes in turbulent transport in the upper part of the PBL. Indeed, the *q*_{t} profile is most well mixed in EDMF_{UD}, signaling the importance of downdraft moisture transport. This is consistent with the hypothesis in Sušelj et al. (2013), suggesting that the inclusion of downdrafts could increase vertical mixing in the upper part of the boundary layer. In STBL, mixing from the surface provides moisture and entrainment from the free troposphere dries the boundary layer. However, in the heat profile, both the surface and entrainment from the free troposphere heats the boundary layer. We find here that downdrafts help provide stronger moisture mixing near cloud top and keep the bias in total moisture low. In addition, EDMF_{UD} has the least bias in boundary layer averged *θ*_{l}, as downdrafts also contribute to transporting warm air in the PBL.

Downdraft model coefficients and final lateral entrainment configuration are chosen to have the best match against LES in the mean field of *θ*_{l}, *q*_{t}, *u*, and *υ*. EDMF_{U} and EDMF_{UD} have the same updraft lateral entrainment configuration.

Comparing EDMF_{U} with SCM results from Sušelj et al. (2013), a resemblance of the updraft transport of heat and moisture is found. The formulations of updrafts are identical except for the added entrainment and dynamical pressure effect near cloud top in EDMF_{U}. It is no surprise that some differences are seen, given the different assumptions made in ED. Specifically, the vertical transport in the middle of the boundary layer is different in the two models. While EDMF_{U} shows positive transport from updraft in the cloudy region for heat, the updraft model in Sušelj et al. (2013) shows a negative heat transport. For moisture, EDMF_{U} produces stronger transport. This is likely due to the added entrainment dynamic effect in our updraft model, different subgrid cloud assumption, and different ED model for the nonconvective environment. In the end, the total heat and moisture transport is similar between the two models as ED compensates for the difference, and they both match LES well.

Comparing EDMF_{UD} with SCM results from Han and Bretherton (2019), we found contrary conclusions for the effect of the downdraft parameterization. While Han and Bretherton (2019) found a slight overprediction for *θ*_{l} and overmixing for *q*_{t} in their DYCOMS experiment, we found slight underprediction for *θ*_{l} and undermixing for *q*_{t}.

### b. CGILS S12 control

Figure 9 shows the mean fields of *θ*_{l}, *q*_{t}, *q*_{l}, cloud fraction, *u*, *υ*, heat flux ($\rho cpw\u2032\theta l\u2032\xaf$), and moisture flux ($\rho L\upsilon w\u2032qt\u2032\xaf$) during h 3 to 4, and the 24-h time series of liquid water path (LWP), boundary layer averaged heat (*θ*_{l}), and moisture (*q*_{t}) for the three tested PBL schemes are shown in Fig. 10. ED shows a strong cold and moist bias throughout the entire simulation. For EDMF_{U}, boundary averaged heat and moisture both follow LES closely up to h 10, then the moisture does not increase as much as in LES. Around h 15, EDMF_{U} begins to cool when compared to LES. This is likely a result of different radiation treatment used in LES and WRF. For EDMF_{UD}, similar trend is observed. Boundary layer averaged heat is warmer and moisture is direr than EDMF_{U}. Both EDMF_{U} and EDMF_{UD} match LWP in LES well. EDMF_{UD} produces a slightly thinner cloud in the first half of the simulation, while EDMF_{U} produces a slightly thicker cloud in the second half of the simulation.

During hours 3 to 4, EDMF_{U} and EDMF_{UD} show small bias in heat and moisture profile, whereas ED is too cold and too moist. This causes the overestimation of LWP in ED. The cloud-top height in EDMF_{UD} is one grid point above ED, likely due to the stronger entrainment flux near cloud top from mass flux. EDMF_{UD} overestimates *u* and underestimates *υ* in the PBL. ED shows similar results as DYCOMS, where the horizontal wind does not have a strong transition between the PBL and the free troposphere. EDMF_{U} shows a very good match in total heat and moisture transport, while EDMF_{UD} has a slightly stronger moisture transport near cloud top. Similar to DYCOMS, ED does not capture cloud-top entrainment flux. Figures 11 and 12 show the vertical flux contribution from each individual component: environment (ED), updrafts, and downdrafts. In both EDMF_{U} and EDMF_{UD}, updraft heat and moisture transport are overestimated. However, in the presence of downdrafts, updraft moisture transport decreases more strongly in-cloud. Downdrafts in EDMF_{UD} partially compensate for these changes, resulting in a similar total transport. Averaged plume properties from EDMF_{UD} are shown in Fig. 13. In CGILS, a good agreement of plume properties is found between the SCMs and LES. Again, we find that simulation results are more sensitive to the modeling of downdraft transport in the upper part of the PBL. In the end, we select model parameters that result in realistic mean profiles of *θ*_{l}, *q*_{t}, *u*, and *υ* for both DYCOMS and CGILS. While the parameterized downdrafts terminate too quickly in DYCOMS, we find that they mostly reach the surface in CGILS.

In the present study, we develop our updraft and downdraft parameterization using their nocturnal properties. The 24-h simulation of CGILS suggests that updrafts and downdrafts may play different roles during the daytime. This is also observed in the study done by Brient et al. (2019). Parameterization of updrafts and downdrafts during the day should be investigated in the future.

### c. Simulation time step and run time

To test the numerical stability of the scheme and the convergence of the results, we also run simulations with different time steps: 5, 10, 20, 30, and 40, and 60 s as shown in Fig. 14. Note that the figures shown in this study use a time step of 40 s. The obtained results confirm that both EDMF_{U} and EDMF_{UD} are not sensitive to the imposed time step changes, proving high robustness of the scheme. The LWP, and the boundary layer-averaged heat and moisture amounts all converge to the same values at the end of the simulation. Additionally, we record simulation run times normalized by the ED-simulation run time for different time steps (Table 2). On average, including the updrafts slows the simulation down by approximately 5%, while including both updrafts and downdrafts slows it down by 7%. This indicates that EDMF is a numerically inexpensive scheme.

## 5. Summary and conclusions

In this study, we investigated the role of nonlocal transport on the development and maintenance of the STBL in coarse-resolution atmospheric models. A special emphasis has been put on the evaluation of downdraft contribution, recently suggested as an important missing element of convection/turbulence parameterizations (Chinita et al. 2018; Davini et al. 2017; Brient et al. 2019). A new parameterization of cloud-top-triggered downdrafts has been proposed along with a complementary parameterization of surface-driven updrafts. The parameterization was validated against large-eddy simulations of two marine stratocumulus cases: DYCOMS and CGILS. The applied nonlocal mass-flux scheme is part of the stochastic multiplume EDMF approach decomposing the turbulent transport into the local and nonlocal contributions. The local transport in the boundary layer is represented by the MYNN scheme. The EDMF scheme has been implemented and tested in the WRF single-column modeling framework.

In the new parameterization, the thermodynamic and dynamic properties of downdrafts are controlled by stochastic lateral entrainment affecting their dilution along the vertical development. The number of downdrafts is fixed to 10, and all downdrafts are assumed to start randomly in the upper half of cloud layer, with the total starting area of approximately 15%, similarly to updrafts. The strength of the downdraft vertical velocity is formulated as a combined effect of the intensity of the surface-driven updrafts and cloud-top radiative cooling. The starting downdraft thermodynamic properties are proportional to the entrainment flux at the STBL top, which is determined by the jump values of heat or moisture across the inversion.

To evaluate the importance of the updraft and downdraft contributions, we run three different SCM simulations for the tested STBLs: without mass flux (ED), with updrafts only (EDMF_{U}), and with both updrafts and downdrafts (EDMF_{UD}). When there is no mass-flux (neither updraft nor downdraft), ED underestimates the cloud-top entrainment flux, yielding a cold and moist bias that leads to a strong overestimation of LWP. The inclusion of updrafts increases the cloud-top entrainment flux and keeps the mean STBL profiles more well-mixed and reduces the temperature and moisture biases. We find that including downdrafts increases vertical mixing in the upper part of the boundary layer especially for *q*_{t}, and it results in a warmer and drier STBL than for EDMF_{U}. Overall, the proposed parameterization reproduces the LES profiles because of the addition of downdraft heat and moisture transport in the WRF SCM. However, we find that differences in EDMF_{U} and EDMF_{UD} are not significant.

Based on the results from the two STBL cases, we conclude that, for the tested version of the WRF Model, it is necessary to include updrafts as part of the nonlocal mass-flux as ED does not represent correctly the cloud-top entrainment flux. An addition of downdrafts shows some improvements in these two cases. However, further investigations are needed to determine whether downdrafts play a greater role in different meteorological conditions. We hypothesize that ED would have a better match with LES when there is less cloud-top entrainment (e.g., when the PBL is less turbulent), and that the inclusion of downdrafts would be necessary when surface fluxes are small. A recent study by Matheou and Teixeira (2019) compared various LESs of STBL for different physical processes included and concluded that surface fluxes, surface shear, and cloud-top radiative cooling all contribute substantially to the turbulence in STBL. Whether the EDMF parameterization responds similarly in such conditions will be investigated in the future.

## Acknowledgments

Parts of this research were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, and were supported by the U.S. Department of Energy, Office of Biological and Environmental Research, Earth System Modeling; and the NASA MAP Program. E. Wu acknowledges the support provided by the JPL Education Office. We thank Mónica Zamora Zapata and Thijs Heus for constructive comments. We also thank Minghua Ong for proofreading the manuscript. The WRF Model is freely available at https://github.com/wrf-model/WRF, the modifications made in this paper can be found at https://github.com/elynnwu/EDMF_JPL.

## REFERENCES

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Adv. Model. Earth Syst.*

*Geophys. Res. Lett.*

*Geophys. Res. Lett.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*J. Climate*

*J. Atmos. Sci.*

*J. Adv. Model. Earth Syst.*

*Wea. Forecasting*

*Mon. Wea. Rev.*

*Mon. Wea. Rev.*

*J. Geophys. Res.*

*Proc. Workshop on Parametrization of Clouds in Large-Scale Models*, Reading, United Kingdom, ECMWF, 93–101

*Quart. J. Roy. Meteor. Soc.*

*Geophys. Res. Lett.*

*Quart. J. Roy. Meteor. Soc.*

*Quart. J. Roy. Meteor. Soc.*

*J. Climate*

*Quart. J. Roy. Meteor. Soc.*

*Atmos. Chem. Phys.*

*Mon. Wea. Rev.*

*Annu. Rev. Fluid Mech.*

*J. Atmos. Sci.*

*J. Adv. Model. Earth Syst.*

*Rev. Geophys. Space Phys.*

*Bound.-Layer Meteor.*

*J. Meteor. Soc. Japan*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Adv. Model. Earth Syst.*

*J. Adv. Model. Earth Syst.*

*J. Appl. Meteor. Climatol.*

*Bull. Amer. Meteor. Soc.*

*14th Symp. on Boundary Layer and Turbulence*, Aspen, CO, Amer. Meteor. Soc., 4.16, https://ams.confex.com/ams/AugAspen/techprogram/paper_14840.htm

*J. Atmos. Sci.*

*Quart. J. Roy. Meteor. Soc.*

*Quart. J. Roy. Meteor. Soc.*

*Bull. Amer. Meteor. Soc.*

*Mon. Wea. Rev.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*Quart. J. Roy. Meteor. Soc.*

*14th Symp. on Boundary Layer and Turbulence*, Aspen, CO, Amer. Meteor. Soc., P4.12, https://ams.confex.com/ams/AugAspen/techprogram/paper_14915.htm

*J. Climate*

*J. Atmos. Sci.*

*Mon. Wea. Rev.*

*J. Climate*

*Sol. Energy*

*J. Atmos. Sci.*

*Sol. Energy*

*Nat. Climate Change*

*Sol. Energy*