the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Improving the representation of shallow cumulus convection with the simplified-higher-order-closure–mass-flux (SHOC+MF v1.0) approach

### Maria J. Chinita

### Mikael Witte

### Marcin J. Kurowski

### Joao Teixeira

### Kay Suselj

### Georgios Matheou

### Peter Bogenschutz

Parameterized boundary layer turbulence and moist convection remain some of the largest sources of uncertainty in general circulation models. High-resolution climate modeling aims to reduce that uncertainty by explicitly attempting to resolve deep moist convective motions. An example of such a model is the Simple Cloud-Resolving E3SM Atmosphere Model (SCREAM) with a target global resolution of 3.25 km, allowing for a more accurate representation of complex mesoscale deep convective dynamics. Yet, small-scale planetary boundary layer turbulence and shallow convection still need to be parameterized, which in SCREAM is accomplished through the turbulent-kinetic-energy-based (TKE-based) simplified higher-order closure (SHOC) – a simplified version of the assumed-double-Gaussian-PDF (probability density function) higher-order-closure method. In this paper, we implement a stochastic-multiplume mass-flux (MF) parameterization of dry and shallow convection in SCREAM to go beyond the limitations of double-Gaussian-PDF closures and couple it to SHOC (SHOC+MF). The new parameterization implemented in a single-column model type version of SCREAM produces results for two shallow cumulus convection cases (marine and continental shallow convection) that agree well with the reference data from large-eddy simulations, thus improving the general representation of the thermodynamic quantities and their turbulent fluxes as well as cloud macrophysics in the model. Furthermore, SHOC+MF parameterization shows weak sensitivity to the vertical grid resolution and model time step.

In general circulation models (GCMs), subgrid physical processes need to be
parameterized due to the typical horizontal resolutions of
GCMs – 𝒪(10^{2}) km. Traditionally, the
turbulent transport in the dry planetary boundary layer (PBL) is represented
by a downgradient eddy-diffusivity approach sometimes combined with a
countergradient flux term to account for the strong nonlocal transport in
the dry convective boundary layer (CBL)
(e.g.,
Deardorff, 1966; Han and Pan, 2011; Teixeira et al., 2004; Holtslag and Moeng,
1991; Stevens, 2000). For shallow cumulus, the transport is often represented
by a separate cumulus parameterization based on the mass-flux approach
(e.g.,
Betts, 1973; Tiedtke, 1989; Yoshimura et al., 2015; Beljaars et al., 2018). Such
parameterizations often require cloud-base closures and trigger functions.
This, combined with the standard GCM modular structure (i.e., GCMs resort to
several independent parameterizations to represent the transport that
happens continuously in the real atmosphere) increases uncertainties and
biases in GCMs
(e.g.,
Teixeira et al., 2008; Sherwood et al., 2014; Schneider et al., 2017).

During the last 2 decades, unified parameterizations have been proposed and implemented in GCMs to reduce some of the issues associated with conventional modular approaches. Unified parameterizations aim to represent the continuous and evolving turbulent transport across the different PBL regimes, e.g., from dry to shallow cumulus convection, in a consistent manner. Two promising approaches emerged to unify boundary layer turbulence and moist convection: eddy-diffusivity mass-flux (EDMF) methods and higher-order closures (HOCs) based on assumed probability density functions (PDFs). Examples of assumed-PDF schemes include the cloud layers unified by binormals (CLUBB; Golaz et al., 2002) and intermediately prognostic higher-order closure (IPHOC; Cheng and Xu, 2006, 2008), where both schemes assume a double-Gaussian PDF to represent subgrid-scale variability of vertical velocity, temperature, and moisture and, therefore, parameterize PBL turbulence and clouds. A key advantage of HOC-PDF schemes is that cloud macrophysical properties and higher-order moments are diagnosed from the joint PDF in a self-consistent manner. A critical downside is that most HOCs are usually computationally expensive as they require at least seven prognostic equations for second- and third-order moments depending on the chosen PDF. To reduce computational costs, the simplified higher-order closure (SHOC; Bogenschutz and Krueger, 2013) was proposed, for which the higher-order moments needed to construct the PDF are diagnosed instead of prognosed.

The EDMF approach is based on the unification of concepts typically used for the parameterization of boundary layer turbulence (eddy diffusivity) and moist convection (mass flux). It was first proposed for dry convective PBLs (Siebesma and Teixeira, 2000; Teixeira and Siebesma, 2000; Siebesma et al., 2007) and later extended to shallow (e.g., Soares et al., 2004; Neggers, 2009; Rio and Hourdin, 2008; Suselj et al., 2013, 2019a; Tan et al., 2018) and deep convection (Suselj et al., 2019a; Cohen et al., 2020), with the latter representing fully unified parameterizations. In a nutshell, the EDMF approach combines the eddy-diffusivity (ED) and mass-flux (MF) parameterizations, where ED represents the local non-convective mixing and MF represents the nonlocal transport via coherent motions such as updrafts. The stochastic-moist-multiplume mass-flux approach (Suselj et al., 2013, 2019a, b) consists of a fully unified EDMF parameterization of PBL turbulence, as well as dry and moist convection (both shallow and deep) without the usage of trigger functions or cloud-base closures. In the EDMF approach's most recent version, the updrafts are coupled to a simple microphysical scheme allowing for precipitating updrafts. A portion of the updrafts' precipitation falls to the surface, and the remaining forms downdrafts that may lead to cold pools. Although the precipitating EDMF version is somewhat complex, especially in comparison with its non-precipitating version, it is still fairly computational efficient, making it a strong parameterization candidate for any GCM. Several EDMF versions have been successfully implemented and evaluated in both climate GCMs (Kurowski et al., 2019b; Witte et al., 2022) and operational numerical weather prediction models (e.g., Kohler et al., 2011; Suselj et al., 2014, 2021; Han et al., 2016; Olson et al., 2019).

Despite the recent advances in unified parameterizations, parameterized convection remains one of the largest sources of uncertainty in GCMs. Thus, high-resolution climate modeling (e.g., cloud-resolving models, CRMs) is emerging as a pathway to reduce that uncertainty by explicitly resolving some deep convection. An example of such a model is the Simple Cloud-Resolving E3SM Atmosphere Model (SCREAM) with a target global resolution of 3.25 km (Caldwell et al., 2021), allowing for a more accurate representation of complex mesoscale deep convective dynamics. Nevertheless, small-scale PBL turbulence and shallow convection still need to be parameterized, which is accomplished in SCREAM using the HOC-PDF scheme SHOC for computational efficiency.

Recent studies (Firl and Randall, 2015; Fitch, 2019) showed that shallow cumulus convection is not properly represented by HOC-PDF schemes due to limitations of the assumed-double-Gaussian PDF in representing high skewness and kurtosis of the distributions. Moreover, higher-order moments and cloud statistics appear to only be properly represented when a larger number of Gaussians is used in the joint PDF, which increases its already expensive computational cost. An alternative solution was recently proposed in Witte et al. (2022), where CLUBB is combined with multiple stochastic MF plumes, leading to a modified CLUBB+MF parameterization where the plumes represent the extreme tail of the joint distribution, which is not represented by CLUBB (see their Fig. 3). Furthermore, their results showed a large improvement of the higher-order moments for two benchmark shallow cumulus convection cases. Thus, the multiple MF plumes offer a physics-based and cost-effective solution by representing the extreme values of the joint distribution not well captured by the assumed PDF.

Here, our main goal is to improve the representation of shallow cumulus convection in SCREAM by merging SHOC with multiple stochastic MF plumes, thereby creating a unified simplified-higher-order-closure–mass-flux (SHOC+MF) parameterization. In our framework, SHOC represents the local mixing and MF the strong nonlocal mixing. The details of the implementation are described in Sect. 2, and the large-eddy simulation (LES) data are described in Sect. 3. In Sect. 4, we discuss its performance in single-column model (SCM) mode for quasi-steady-state trade-wind maritime shallow cumulus convection and a land diurnal cycle of shallow convection. Lastly, sensitivity tests to vertical grid resolution and model time steps are also carried out. Conclusions are presented in Sect. 5.

We combine a stochastic-moist-multiplume MF scheme with SHOC (Bogenschutz and Krueger, 2013) in SCREAM. The coupling of MF and SHOC has the potential to improve the representation of the mean thermodynamic structure, higher-order moments (e.g., the turbulent fluxes), and cloud macrophysics quantities by adding the contribution of nonlocal transport during intense convection.

## 2.1 Host model description

SCREAM emerged as one of the next-generation development efforts of the Energy Exascale Earth System Model (E3SM) project led by the U.S. Department of Energy (DOE) to help guide future energy-sector decisions in light of the current long-term trends due to global warming (Golaz et al., 2019; Caldwell et al., 2021). SCREAM is presently still in development but in its final form aims to represent the next generation of global convection-permitting models (GCPMs) by running faster than previous GCPMs due to its performance portability from CPU to GPU machines. To achieve this, the GCM E3SMv1 model serves as a template and is being rewritten from Fortran to C$++$. SCREAM is based on a nonhydrostatic spectral element dycore and parameterizes turbulence, shallow moist convection, microphysics, radiation, and aerosols (see Caldwell et al., 2021, for a detailed model description). Its target global resolution is 3.25 km.

Here, we use the SCREAM version dyamond2-try1 still written in Fortran and released in October 2020 (https://github.com/E3SM-Project/scream/releases/tag/dyamond2-try1, last access: 20 March 2023), with two modifications: (1) multiplume MF scheme implemented in SHOC's code base and (2) correction of a bug in the SCM spectral element dynamical core that was producing a strong unphysical temperature cold bias. This bug has been fixed in the current development code base (https://github.com/E3SM-Project/E3SM/pull/4027, last access: 20 March 2023).

Our initial assessment of the MF implementation is performed using SCREAM in a SCM framework (Bogenschutz et al., 2020). We are currently migrating the MF component module of the SHOC+MF parameterization to SCREAMv0 (version used in Caldwell et al., 2021), and preliminary results show no significant differences relative to the results presented here.

## 2.2 EDMF parameterization

In weather and climate models, the prognostic equation of the thermodynamic variables depends on the vertical divergence of the turbulent flux in addition to the advective tendencies and diabatic processes:

where $\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}$ represents the prognostic horizontally averaged
thermodynamic variable, here taken as the liquid water potential temperature
and total water mixing ratio, $\mathit{\varphi}=\left\{{\mathit{\theta}}_{l},{q}_{t}\right\}$;
*w* is the vertical velocity; and the primes denote fluctuations with
respect to the mean $\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}$. In the convective boundary layer, the
turbulent flux corresponds to a combination of small-scale and large-scale
coherent turbulent structures and can be decomposed as

where the subscripts “e” and “u” denote the environment and the strong updrafts,
respectively. In the EDMF approach, the following approximations are usually
made: (1) the first term is parameterized with the ED approach, (2) the second
term is neglected because the environmental and grid-mean values are
approximately equal (i.e., ${w}_{\mathrm{e}}\approx \stackrel{\mathrm{\u203e}}{w}$) following the assumption
of small updraft area (i.e., *a*_{u}≪1), and (3) the third term vanishes
because the updrafts are assumed horizontally homogeneous and their internal
covariances are zero. The fourth term is commonly known as the mass-flux
contribution since ${M}_{u}\equiv {a}_{u}\left({w}_{u}-\stackrel{\mathrm{\u203e}}{w}\right)$. Thus,
Eq. (2) can be simplified to

which encapsulates the eddy-diffusivity mass-flux (EDMF) approach
(e.g.,
Siebesma et al., 2007; Suselj et al., 2013). Here, the eddy-diffusivity
coefficient, *K*_{ϕ}, is defined according to SHOC's formulation
(Bogenschutz and Krueger, 2013),
and the MF contribution follows the stochastic-moist-multiplume MF scheme
introduced in Suselj et al. (2019a). Thus, the updraft horizontal grid area is
partitioned into multiple updrafts, and Eq. (3) is rewritten as

where ${\sum}_{n=\mathrm{1}}^{N}{M}_{n}\left({\mathit{\varphi}}_{n}-\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}\right)={\sum}_{n=\mathrm{1}}^{N}{a}_{n}{w}_{n}({\mathit{\varphi}}_{n}-\stackrel{\mathrm{\u203e}}{\mathit{\varphi}})$, *N* is the
user-selected total number of updrafts (here, *N*=40 updrafts), *a*_{n} is
the area fraction of the *n*th updraft, and *w*_{n} and *ϕ*_{n} are the
vertical velocity and thermodynamic property of the *n*th updraft. The updraft
properties are defined according to the updraft model described below
(Sect. 2.3).

It is becoming more common to include downdrafts in EDMF parameterizations (e.g., Wu et al., 2020; Han and Bretherton, 2019) mostly due to their relevance to turbulent transport in stratocumulus-topped boundary layers (Chinita et al., 2018; Brient et al., 2019). Despite this, Wu et al. (2020) showed that the inclusion of updrafts is sufficient to represent the vertical thermodynamic structure and turbulent fluxes of non-precipitating stratocumulus, which is in agreement with the findings reported in Matheou and Teixeira (2019), where the authors showed using LES results that the surface buoyancy and wind shear are as important for turbulence production as cloud-top radiative cooling. Combined with a need for computational efficiency, these recent findings led us to neglect downdrafts in our current MF implementation.

## 2.3 Updraft model

The updraft model closes the multiplume EDMF parameterization and defines
the vertical evolution of an updraft's vertical velocity and thermodynamic
properties. Here, we follow the updraft model described in
Suselj et al. (2019a). Accordingly, at the
surface, we release *N* independent, steady-state buoyancy-driven updrafts with
surface vertical velocities sampled from the right tail of an assumed-Gaussian PDF, with values ranging between *w*_{min} and *w*_{max}, here
defined as $\mathrm{1.5}{\mathit{\sigma}}_{w}<\phantom{\rule{0.125em}{0ex}}{w}_{n}<\mathrm{3}{\mathit{\sigma}}_{w}$, where *σ*_{w} is
the vertical velocity standard deviation (note that the interval [1.5*σ*_{w}, 3*σ*_{w}] corresponds to a total updraft surface fraction area
equal to 6.65 %, in agreement with the sensitivity analysis to the surface
updraft area presented in Suselj et al., 2019a). The tail of the velocity PDF, i.e., the interval [1.5*σ*_{w}, 3*σ*_{w}], is discretized into *N* equidistant bins, and the mean
vertical velocity value of each bin is associated with a corresponding
updraft (*N* is the total number of updrafts). The surface thermodynamic
properties of each updraft are computed by integrating the joint-normal
PDF (*θ*_{lu}, *q*_{tu}, *w*_{u}) over the updraft's velocity bin
(see Suselj et al., 2019b, for details on the
joint-normal PDF characterization). Here, we use *N*=40 updrafts. The
number of updrafts was chosen based on a sensitivity analysis of SHOC+MF
to its value in which SHOC+MF showed weak sensitivity to *N*>30 updrafts
(not shown). Note that a small number of updrafts can produce noisier
results due to the lateral entrainment's stochasticity
(Suselj et al., 2019b).

The vertical evolution of the *n*th updraft depends on surface properties and
lateral entrainment as follows:

where $\mathit{\varphi}=\left({\mathit{\theta}}_{l},{q}_{t}\right)$ and *ε*_{n} is
the entrainment rate of the *n*th updraft. Thus, Eq. (5) represents the
dilution of *ϕ*_{n} by lateral entrainment of environmental air
$\stackrel{\mathrm{\u203e}}{\mathit{\varphi}}$; the environmental air properties are assumed equal to the
grid-mean values following Kurowski et al. (2019a). The vertical velocity of the *n*th updraft is determined by

where *a*_{w}=1 and *b*_{w}=1.5 are constants and *B*_{n} is the updraft's
buoyancy given by ${B}_{n}=g\left({\mathit{\theta}}_{\mathit{\upsilon},n}/{\stackrel{\mathrm{\u203e}}{\mathit{\theta}}}_{\mathit{\upsilon}}-\mathrm{1}\right)$, where
*θ*_{υ} is the virtual potential temperature. The boundary condition
values needed to integrate Eqs. (5) and (6), i.e., the surface
thermodynamic properties (*w*_{n}|_{s}, ${\mathit{\theta}}_{\mathit{\upsilon},n}{|}_{\mathrm{s}}$, ${q}_{t,n}{|}_{\mathrm{s}}$), are computed as in
Suselj et al. (2019a), and their standard
deviation values (*σ*_{w}, ${\mathit{\sigma}}_{{\mathit{\theta}}_{\mathit{\upsilon}}}$, ${\mathit{\sigma}}_{{q}_{t}}$) follow
Suselj et al. (2019b). Note that
${\mathit{\theta}}_{l,n}{|}_{\mathrm{s}}$ is defined with respect to ${\mathit{\theta}}_{\mathit{\upsilon},n}{|}_{\mathrm{s}}$ as ${\mathit{\theta}}_{l,n}{|}_{\mathrm{s}}={\mathit{\theta}}_{\mathit{\upsilon},n}{|}_{\mathrm{s}}/(\mathrm{1}+\mathrm{0.61}{q}_{t,n}{|}_{\mathrm{s}}$) assuming ${\mathit{\theta}}_{l,n}{|}_{\mathrm{s}}\equiv {\mathit{\theta}}_{n}{|}_{\mathrm{s}}$ (subscript
“s” denotes surface). The numerical discretization of Eqs. (5) and (6) follows
that described in Suselj et al. (2014).

Lastly, the lateral entrainment of the *n*th updraft is defined as a stochastic
process (Romps and Kuang, 2010;
Suselj et al., 2019b):

where *ε*_{0} is the fraction of entrained mass flux during each
entrainment event, here set to *ε*_{0}=0.2; 𝒫_{n} is a random
number drawn from the Poisson distribution that represents the number of
entrainment events for a given average event frequency equal to
*L*_{ε}, and Δ*z* is the thickness of the respective layer.
Note that we evaluate 𝒫_{n} and *ε*_{n} for each updraft
independently. Following Suselj et
al. (2019b), the entrainment length scale is defined as a function of the
depth of the CBL including the cloud layer when present, *h*_{CBL}:

where *a*=1.25 m${}^{\mathrm{1}/\mathrm{2}}$ is a constant and *h*_{CBL} is defined as the
model level where the vertical heat flux vanishes ($\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\theta}}_{l}^{\prime}}\approx \mathrm{0}$). In agreement with previous studies
(e.g., Böing et al., 2012; Takahashi et al., 2021), the entrainment length scale *L*_{ε} is larger for deeper clouds (i.e., higher *h*_{CBL}) as these tend to be
wider and thus better protected from the environment, leading to smaller
entrainment rates. Note that diagnosing *L*_{ε} as the
square root of *h*_{CBL} allows for continuous adjustment of *ε*_{n} as a function of the CBL state; i.e., the entrainment rate is reduced
for deeper CBLs, allowing the updrafts to reach higher vertical levels and
vice versa for shallower CBLs, which is particularly important to represent
the strong diurnal cycle over land while remaining insensitive to small
oscillations of *h*_{CBL}.

Condensation in each updraft takes place if the updraft water vapor reaches saturation. The MF contribution to the total cloud fraction corresponds to the sum of the area fraction of the updrafts that condense, and the MF contribution to the total cloud water is defined as the area average of the cloud water of all moist updrafts.

## 2.4 SHOC

In SCREAM, boundary layer turbulence and moist shallow convection are parameterized by SHOC (Bogenschutz and Krueger, 2013). SHOC is considered a simplified assumed-PDF-based scheme because the second-order moments needed to construct the PDF are diagnosed instead of prognosed to increase computational efficiency. Accordingly, the turbulent fluxes $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\theta}}_{l}^{\prime}}$ and $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{q}_{t}^{\prime}}$ are estimated following an eddy-diffusivity approach:

where $\mathit{\varphi}=\left\{{\mathit{\theta}}_{l},{q}_{t}\right\}$ and *K*_{ϕ} represents
the eddy-diffusivity coefficient for heat. It is important to note that SHOC
has been updated since
Bogenschutz and Krueger (2013)
to improve numerical stability and overall performance to better represent
the various regimes present in a GCM. For instance, the formulation of the
turbulence length scale has been revised and now follows a continuous
formulation instead of two separate definitions for the sub-cloud and cloud
layers as documented in
Bogenschutz and Krueger (2013).
Nevertheless, the SHOC version used in SCREAM exhibits similar scientific
performance to the original formulation. For completeness, the turbulence
length scale is defined as

where *k* is the von Karman constant, *e* is the turbulent kinetic energy, and
*l*_{c}=0.5 is a tunable length scale factor. The constant *δ* is
defined as *δ*=1 if the Brunt–Väisälä frequency *N*^{2}>0 or *δ*=0 if *N*^{2}≤0, where ${N}^{\mathrm{2}}=\frac{g}{\stackrel{\mathrm{\u203e}}{{\mathit{\theta}}_{\mathit{\upsilon}}}}\frac{\partial \stackrel{\mathrm{\u203e}}{{\mathit{\theta}}_{\mathit{\upsilon}}}}{\partial z}$. The asymptotic value of the length scale
*L*_{∞} is defined following Blackadar (1962) as
${L}_{\mathrm{\infty}}=\mathrm{0.1}{\int}_{\mathrm{0}}^{\mathrm{\infty}}\sqrt{e}z\mathrm{d}z/{\int}_{\mathrm{0}}^{\mathrm{\infty}}\sqrt{e}\mathrm{d}z$. Lastly, the eddy turnover timescale *τ* is defined as $\mathit{\tau}=h/{w}_{\ast}$, where *h* is
the boundary layer depth calculated according to
Holtslag and Boville (1993), and the
convective velocity scale ${w}_{\ast}^{\mathrm{3}}=\mathrm{2.5}\frac{g}{\stackrel{\mathrm{\u203e}}{{\mathit{\theta}}_{\mathit{\upsilon}}}}{\int}_{\mathrm{0}}^{h}\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\theta}}_{\mathit{\upsilon}}^{\prime}}\mathrm{d}z$. If the boundary
layer is stable (i.e., ${w}_{\ast}^{\mathrm{3}}<\mathrm{0}$), then *τ*=100 s.

## 2.5 Coupling of SHOC and multiplume mass-flux parameterizations

We implement the stochastic-multiplume MF scheme in SCREAM by coupling it to SHOC. Thus, the multiplume MF contribution (second term of the right-hand side of Eq. 11) is added to SHOC's numerical solver for the mean thermodynamic variables, $\mathit{\varphi}=({\mathit{\theta}}_{l},{q}_{t}$), according to the following one-dimensional prognostic equation:

where *K*_{ϕ} is the eddy-diffusivity coefficient, *a*_{n} is the area
fraction of the *n*th updraft, *w*_{n} and *ϕ*_{n} are the vertical
velocity and the *ϕ* value in the *n*th plume, and the overbar denotes a
grid-mean value. The SHOC term (first term of the right-hand side of
Eq. 11) represents the time tendency of *ϕ* due to the
downgradient diffusion of the mean field, and the MF term takes into account
the nonlocal transport due to strong convection as discussed in Sect. 2.2.
The prognostic Eq. (11) is discretized according to the semi-implicit
forward-in-time centered-in-space scheme and solved using
the Richtmyer and Morton (1967) method (see
Kurowski et al., 2019b, for the discretized
form of Eq. 11). Note that the surface boundary conditions of Eq. (11) (i.e., the surface fluxes of the thermodynamic variables, ${\stackrel{\mathrm{\u203e}}{{w}^{\prime}{\mathit{\varphi}}^{\prime}}}_{\mathrm{s}}$) are either calculated by the surface layer parameterization or are
prescribed and not modified by the MF component. Results of the coupling of
MF and SHOC are denoted as “SHOC+MF” in the next sections.

The shortwave and longwave radiation and the large-scale cloud microphysics parameterizations are switched off for these experiments – basically because these processes are believed to be of secondary importance for these shallow convection case studies and as such are also off for the LES experiments. However, cloud fraction and water are calculated at every model level for diagnostic purposes. In SCREAM, the cloud macrophysical properties cloud fraction and liquid water mixing ratio are estimated using the SHOC PDF. Here, the combination of the grid-mean cloud properties is done as a simple weighted sum of the SHOC and MF contributions:

where CF_{SHOC} and ${q}_{{l}_{\mathrm{SHOC}}}$ are diagnosed from SHOC's assumed-double-Gaussian distribution, and *a*_{n} and ${q}_{{l}_{n}}$ are the fractional
area and condensate loading of the *n*th updraft. Note that, in practice,
although the algorithm also imposes CF≤1, this value is not reached in
these simulations because of the low values typical of these shallow
convection case studies. These low cloud values, the overall cloud vertical
structure of these shallow convection regimes, and the fact that the
radiation parameterization is off are the key reasons for not using more
complex cloud overlap algorithms in these estimates.

## 2.6 Single-column model simulations

The ability of the unified SHOC+MF parameterization to represent shallow
cumulus boundary layers has been investigated by simulating benchmark cases
including the shallow cumulus Barbados Oceanographic and Meteorological
Experiment (BOMEX; Siebesma et al., 2003), quasi-steady-state warm maritime
shallow convection over the Atlantic Ocean in June 1969, and the Atmospheric
Radiation Measurement (ARM) shallow cumulus case (Brown et al., 2002), diurnal cycle of warm shallow convection over land at the Southern
Great Plains site of the ARM program on 21 June 1997. The two cases were
simulated using SCREAM in a SCM framework in which we used the intensive-observation-period (IOP) forcing files available in the E3SM SCM library
(Bogenschutz et al., 2020) with prescribed
horizontal large-scale forcing and surface turbulent fluxes. It is worth
noting that we modified the ARM case forcing file to run the model with a
30 min time step (i.e., the ARM forcing file available in the E3SM SCM
library contains values at every 20 min). Also, the SCM reads the wind
information from the forcing file at every host model time step; however, for
the ARM case, the large-scale advective tendencies of the winds were not
available when the case was set up (Brown et al., 2002), and, consequently, the
time-varying ** u** profiles in the forcing file were set equal to the initial
profiles which are constant with height. Resetting the

**profile to the initial vertically constant profile at every host model time step interferes with the development of the turbulent-kinetic-energy (TKE) field through the shear production term. To circumvent this issue, we replaced the**

*u***profiles in the forcing file with the**

*u***profiles from our LES reference data; the meridional wind component**

*u***is zero in the ARM case. Note that this is specific to the SCM used here and to the ARM case as the large-scale advective tendencies of the winds were not available when the ARM case was set up (Brown et al., 2002).**

*v*We kept the default host model setup but deactivated the deep-convection,
large-scale microphysics and radiation schemes to allow for a more
straightforward comparison with our LES reference data. The dynamic and
physics time steps are equal to 30 and 5 min, respectively. For
our initial implementation and performance evaluation, we used a 72-layer
vertical grid (L72) with 21 levels resolving the bottom 3 km. In Sect. 4.3, we conduct tests to quantify sensitivity to the vertical grid
resolution and to the time step using BOMEX. Thus, for the vertical grid
resolution, we assess the sensitivity of the results using L72 and a
relatively finer 128-layer vertical grid (L128) with twice as many
grid cells resolving the bottom 3 km (40 grid cells). For the time step
sensitivity, we compare the results using the L128 and dynamics and physics
time steps equal to 30 and 5 min (300 s), respectively,
with dynamics and physics time steps both equal to 75 s, which
resembles the configuration used in Caldwell et al. (2021) for the first
global results of SCREAM in convection-permitting mode (Δ*x*=3.25 km).

We evaluate our SHOC+MF parameterization by comparing it to LES output of the same benchmark cases. These LES reference data are acquired with the LES model of Matheou and Chung (2014). Table 1 summarizes the LES runs and their configurations. The computational domain is doubly periodic in the horizontal directions, and all grids are uniform and isotropic ($\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{\Delta}z$). The simulations have different domain sizes in the vertical adjusted to their respective boundary layer depths. A Rayleigh damping layer is imposed near the domain top to mitigate gravity wave reflection, and the surface turbulent fluxes are prescribed as in the SCREAM SCM. The momentum and scalar advection terms are discretized according to the sixth-order fully conservative centered scheme of Morinishi et al. (1998) adapted for the anelastic approximation (Matheou et al., 2016). The subgrid-scale (SGS) turbulence is represented by the buoyancy-adjusted stretched-vortex SGS model (Chung and Matheou, 2014). Precipitation is neglected in the LES model according to the case descriptions (Brown et al., 2002; Siebesma et al., 2003), and all water condensate is assumed suspended using an all-or-nothing saturation adjustment scheme based on the local grid-mean state. The simulations are carried out in the frame of reference of the domain-mean horizontal wind to reduce numerical errors (Lamaakel and Matheou, 2021). The LES model has been successfully used in previous studies spanning several meteorological conditions (Chung et al., 2012; Matheou and Chung, 2014; Matheou, 2018; Matheou and Teixeira, 2019; Couvreux et al., 2020; Chinita et al., 2022a, b).

We compare the results of SHOC and SHOC+MF against the LES reference data
for the benchmark cases listed in Sect. 2.6. Note that in SHOC+MF, we
reduced SHOC's turbulence mixing length scale relative to SHOC alone to
prevent excessive mixing. This was done by increasing the tunable length
scale factor *l*_{c} from 0.5 to 1 in Eq. (10). Thus, *l*_{c}=0.5 in the
SHOC experiments, and *l*_{c}=1 in the SHOC+MF experiments. All other
model and parameterization configurations were kept the same for all
simulations shown here. Lastly, for the simulations of SHOC alone, we used
the same tunable constants from the global high-resolution simulation
presented in
Caldwell
et al. (2021).

## 4.1 Trade-wind maritime shallow cumulus

Figure 1 shows results for BOMEX averaged over simulated hours 4 to 6. The thermodynamics profiles are generally similar but with some noticeable differences: SHOC is colder above cloud base and warmer near the cloud top relative to the LES (Fig. 1a), and it is moister above cloud base and drier near the cloud top (Fig. 1b). This is because SHOC mixes excessively up to ∼ 1 km and does not reproduce a shallow cumulus layer (Fig. 1d and e). Consequently, moisture does not reach the levels where the cloud top should be located and instead it gets trapped between 0.5 and 1 km. In contrast, the turbulent transport of SHOC+MF is very similar to the LES, leading to comparable thermodynamic profiles, except near the cloud top (∼ 1.5 km), where SHOC+MF mixes slightly less, leading to a drier (warmer) layer relative to the LES. The partitioning of turbulent transport between local and nonlocal mixing in SHOC+MF is similar to previous EDMF studies (e.g., Suselj et al., 2013; Kurowski et al., 2019b); i.e., the local mixing dominates in the subcloud layer, and MF takes over in the cloud layer.

Because of the excessive humidity between 0.5 and 1 km in SHOC, the cloud fraction and liquid water mixing ratio are overestimated relative to the LES (Fig. 1c and f). On the other hand, SHOC+MF captures the profiles of both cloud fraction and cloud liquid water content fairly well due to the adequate vertical distribution of the thermodynamic quantities. Note that the cloud fraction and liquid water content of SHOC+MF shown in Fig. 1c and f are calculated as the sum of SHOC and MF contributions.

A key aspect in simulating shallow cumulus with an MF-type parameterization
like SHOC+MF is the accurate representation of the moist updraft
properties (i.e., updraft area, vertical velocity, and the excess of moist
conserved variables). Figure 2 shows the moist updraft properties of the
mass flux of SHOC+MF and the respective LES values based on cloud
(${q}_{l}>\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ kg kg^{−1}) and cloud core (${q}_{l}>\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ kg kg^{−1}, *w*>0, *θ*_{υ}>〈*θ*_{υ}〉, where the angle brackets denote the instantaneous horizontal average of
the LES domain) samplings (Siebesma and Cuijpers,
1995). Since the SHOC+MF turbulent transport is controlled mostly by MF in
the cloud layer (dashed profiles in Fig. 1d and e), the moist updraft
properties should lie close to the LES cloud and cloud core values
(Couvreux
et al., 2010; Suselj et al., 2013; Kurowski et al., 2019b). The SHOC+MF
updraft area agrees reasonably well with the LES values, especially when
considering the relatively coarse vertical grid used here, resulting in just
two grid levels to resolve the sharp increase near the cloud base. The
vertical velocity and the excess of updraft moist conserved variables
relative to the grid-mean values of SHOC+MF are close to the LES profiles,
except in the middle of the cloud layer (∼ 1 km), where the *θ*_{lu} and *q*_{tu} excesses are underestimated. This is due to the slight
overestimation of the SHOC+MF grid-mean *q*_{t} (Fig. 1b) relative to
the LES. Nevertheless, the SHOC+MF moist updraft properties agree well
with the LES, which confirms the suitable behavior of our MF scheme. Note
that the updraft (second Gaussian) moist properties of the SHOC's PDF
are not shown because they are quite small (e.g., maximum *w*_{u}≈ 0.3 m s^{−1}) and vanish around 800 m in agreement with Fig. 1d–e, where the
MF contribution makes up the total turbulent fluxes.

## 4.2 Continental shallow cumulus

The ARM shallow cumulus case represents a diurnal cycle of warm convection over land at the Southern Great Plains site of the ARM program on 21 June 1997 (Brown et al., 2002). The case starts with a morning transition at 11:30 UTC (5:30 LST – local standard time) from a stable boundary layer (negative surface heat flux until 1.5 simulated hours) to a fully developed CBL with a top close to 2.5 km at around 21:00 UTC (15:00 LST). The case represents a typical buoyancy-driven shallow cumulus case where convection is primarily forced by the surface sensible and latent heat fluxes. The nonstationary conditions of the ARM shallow cumulus case make it more challenging to properly simulate than the quasi-steady-state BOMEX case.

Figure 3 shows that the diurnal evolution of hourly mean thermodynamic
profiles is well represented by SHOC+MF, whereas SHOC produces a warm
(dry) bias in the subcloud layer and a cold (moist) bias near the cloud top.
To illustrate the magnitude and temporal evolution of these biases, panels a
and b of Fig. 4 show the temperature and moisture differences relative to
the LES fields, and by the end of the simulation, the temperature (moisture)
bias exceeds 0.5 K (1 g kg^{−1}) in the subcloud layer and 1.5 K (4 g kg^{−1}) near the cloud top. On the other hand, SHOC+MF is able to
reproduce the diurnal evolution of the PBL and cloud layer remarkably well
for both thermodynamic quantities. The largest bias is located near the
surface for the total water mixing ratio with maximum deviations from LES
around +0.5 g kg^{−1} (Fig. 4d).

Figure 5 shows the diurnal evolution of hourly mean profiles of the
turbulent fluxes of liquid water potential temperature (top row) and total
water mixing ratio (bottom row). In agreement with the vertical distribution
of the thermodynamic quantities, SHOC+MF represents both
turbulent transports fairly well, whereas SHOC produces excessive turbulent mixing and
consequently places the cloud top near 3 km at 20:30 UTC – about 1 km deeper
than the LES. Figure 6 shows the differences between the temperature and
moisture turbulent fluxes of SHOC and SHOC+MF relative to the LES. Both
temperature and moisture panels (a–b) confirm SHOC's excessive mixing and
cloud layer deepening and also reveal oscillations, particularly on the
$\stackrel{\mathrm{\u203e}}{{w}^{\prime}{q}_{t}^{\prime}}$ field after hour 5. These oscillations may be due to
the eddy turnover timescale *τ* used in the calculation of SHOC's
turbulence mixing length scale (Eq. 10) in this SCREAM version since
these are not present when a constant timescale is used (e.g., *τ*= 400 s; not shown). However, SHOC+MF is able to reproduce the
turbulent fluxes without these oscillations (Fig. 6c–d) while using the
dynamic timescale, matching the LES reasonably well, except for the last 4 h of simulation where $\stackrel{\mathrm{\u203e}}{{w}^{\prime}{q}_{t}^{\prime}}$ decreases slower (faster)
than the LES in the upper half (lower half) of the boundary layer (Fig. 5h).

Note that the turbulent transport partition between local and nonlocal mixing in SHOC+MF is similar to BOMEX when the cloud layer forms (from simulated hour 5 to hour 12); i.e., the transport is mostly controlled by the local mixing in the subcloud layer, whereas MF dominates in the cloud layer (e.g., dashed red profiles in Fig. 5c and g). Before cloud formation, the local mixing contribution to the turbulent transport is larger and the MF contribution is only significant near the surface. A similar behavior was observed for a dry convection case (case 1 of Siebesma et al., 2007 – not shown), where SHOC properly represented the turbulent transport including the PBL growth but developed a warm bias near the surface – this is a typical pattern of ED-type schemes without MF (e.g., Teixeira and Cheinet, 2004; Siebesma et al., 2007; Witek et al., 2011). The inclusion of MF partially reduced this warm bias by slightly enhancing the turbulent mixing near the surface whilst adjusting its contribution to negligible values away from the surface.

Figure 7 shows the temporal evolution of the cloud fraction. The LES cloud
fraction field is smoother than SHOC and SHOC+MF due to the LES higher
temporal resolution (Δ*t*=1 min vs. 30 min). Nevertheless,
SHOC+MF cloud cover (Fig. 7c) roughly follows the LES, except between
hours 10 and 12 because of its excessive turbulent flux in the cloud layer
(Fig. 5h), resulting in significantly higher cloud fraction values. The
cloud fraction values of SHOC surpass those from LES because of the
excessive moisture content in the cloud layer (Fig. 4b), and its onset
happens about 1 h earlier than in the LES. Conversely, SHOC+MF captures
the cloud layer evolution reasonably well due to a better representation of
the heat and moisture turbulent transports (Figs. 5 and 6c–d).

## 4.3 Sensitivity to vertical grid resolution and time step

The sensitivity of SHOC and SHOC+MF results to the vertical
grid resolution and time step is explored here using the BOMEX case. Results
are similar for the ARM case and are thus omitted. We compare the results of
BOMEX discussed in Sect. 4.1 using the default vertical grid (72 vertical
levels; L72) and a 128-layer vertical grid (L128). This vertical resolution
increase translates to about twice as many grid cells within the CBL
(including the cloud layer). All other aspects of the model configuration
were held unchanged. Lastly, we use the L128 grid to explore the sensitivity
of the model to the time step comparing the default 30 min time step to a
75 s time step. The vertical grid and the time step used in SCREAM's
global simulations will likely be close to L128 and Δ*t*=75 s
(Caldwell et al., 2021), thus the importance of exploring the sensitivity of
SHOC+MF to both configurations.

Figure 8 shows the sensitivity of the temporally averaged vertical profiles of the thermodynamic and cloud macrophysics variables, as well as turbulent fluxes with respect to the vertical grids. The results of SHOC+MF demonstrate low sensitivity to the grid resolution but still a slight improvement when using L128; e.g., unsurprisingly, the sharp increase in cloud fraction near the cloud base is better resolved with L128 (Fig. 8c although not clearly visible). On the other hand, the results of SHOC show a strong sensitivity to the grid resolution – specifically, both heat and moisture turbulent fluxes roughly double in magnitude in the cloud layer (Fig. 8d–e). On a positive note, this increase in turbulent transport warms up the cloud layer relative to the results using L72, which improves the cloud fraction.

The sensitivity with respect to the time step is shown in Fig. 9. The results of SHOC+MF also show low sensitivity to the time step, while SHOC seems to be slightly sensitive to it. Overall, the results of SHOC+MF do not depend on the vertical grid resolution or on the time step. Thus, further tuning does not seem to be necessary for shallow convection when using a different vertical grid or time step.

This study documents the implementation of the stochastic-multiplume mass-flux (MF) parameterization (Suselj et al., 2013, 2019a, b) in the Simple Cloud-Resolving E3SM Atmosphere Model (SCREAM) by coupling it to the simplified-higher-order-closure (SHOC) turbulence and cloud macrophysics scheme. The MF contribution to the total turbulent transport is added to SHOC's numerical solver for the moist conserved thermodynamic variables.

SHOC is a unified assumed-PDF-based scheme that represents both boundary layer turbulence and cloud macrophysics, and while it satisfactorily represents dry convection and stratocumulus layers, it struggles to adequately represent shallow cumulus convection (Firl and Randall, 2015; Fitch, 2019). Following a recent study that showed promising results in solving this issue by combining an MF parameterization with the assumed-PDF scheme CLUBB (Witte et al., 2022), we coupled MF to SHOC to improve the representation of shallow cumulus convection in SCREAM.

Our new scheme (SHOC+MF) was evaluated in a single-column model (SCM) simulation
framework against LES reference data for two shallow cumulus convection
cases: BOMEX – quasi-steady-state warm maritime shallow convection – and
ARM – diurnal cycle of warm shallow convection over land. We also compared
the SHOC+MF results with standard SHOC. In general, SHOC+MF represents
the mean and flux profiles of moist conserved thermodynamic variables well
(liquid water potential temperature, *θ*_{l}, and total water mixing
ratio, *q*_{t}), as well as the cloud macrophysics properties (cloud fraction
and cloud water mixing ratio, *q*_{c}) for shallow cumulus boundary layers.
This represents an improvement versus SHOC alone, since for BOMEX, SHOC does
not reproduce a shallow cumulus layer but rather simulates a structure
similar to a stratocumulus boundary layer, and for ARM, SHOC mixes
excessively up to 3 km, producing a cloud layer too deep, and also
overestimates the cloud macrophysical properties.

We performed a sensitivity analysis to the vertical grid resolution, as well as dynamic and physics time steps for SHOC and SHOC+MF. While SHOC seems to be sensitive to both grid resolution and time step, SHOC+MF showed weak sensitivity to both. Thus, SHOC+MF appears to be robust to changes in the vertical resolution and time step, suggesting there is no need for additional parameter optimization.

In summary, the results of SHOC+MF in SCREAM demonstrate good performance by improving the representation of shallow cumulus convection. Furthermore, the SHOC+MF configuration introduced here is robust enough to properly represent two different shallow cumulus convection cases (i.e., quasi-stationary and non-stationary) regardless of the vertical grid resolution and time step used. Based on these encouraging results, we are currently expanding the evaluation of SHOC+MF to both stratocumulus and deep-convection regimes, as well as to global simulations.

In this study, we used the E3SM model (https://doi.org/10.11578/E3SM/dc.20210927.1, E3SM Project, DOE, 2021), specifically the E3SM SCREAM version dyamond2-try1 released in October 2020 (https://github.com/E3SM-Project/scream/releases/tag/dyamond2-try1, last access: 20 March 2023). The modified code shoc.F90 and shoc_intr.F90, the mass_flux.F90, and the scripts and respective IOP files used to generate the present SCM simulations are archived at https://doi.org/10.5281/zenodo.7011628 (Chinita, 2022a). The SCM and LES output data are archived at https://doi.org/10.5281/zenodo.7011652 (Chinita, 2022b).

MJC: conceptualization; methodology; formal analysis; visualization. MW: conceptualization; methodology. MJK: conceptualization; methodology. JT: conceptualization; methodology; supervision. KS: conceptualization; methodology. GM: provided the large-eddy simulation output data. PB: provided crucial SCREAM code assistance. All authors contributed to interpreting the results, as well as writing and reviewing the manuscript.

The contact author has declared that none of the authors has any competing interests.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). We gratefully acknowledge the support of the U.S. Department of Energy, Office of Biological and Environmental Research, Earth System Modeling (DE-SC0019242). This research used resources from the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under contract no. DE-AC02-05CH11231.

This research has been supported by the Office of Science (grant no. DE-SC0019242).

This paper was edited by Simon Unterstrasser and reviewed by three anonymous referees.

Beljaars, A., Balsamo, G., Bechtold, P., Bozzo, A., Forbes, R., Hogan, R. J., Köhler, M., Morcrette, J. J., Tompkins, A. M., Viterbo, P., and Wedi, N.: The Numerics of Physical Parametrization in the ECMWF Model, Front. Earth Sci., 6, 1–18, https://doi.org/10.3389/feart.2018.00137, 2018.

Betts, A. K.: Non-precipitating cumulus convection and its parameterization, Q. J. Roy. Meteor. Soc., 99, 178–196, https://doi.org/10.1002/qj.49709941915, 1973.

Blackadar, A. K.: The vertical distribution of wind and turbulent exchange in a neutral atmosphere, J. Geophys. Res., 67, 3095–3102, 1962.

Bogenschutz, P. A. and Krueger, S. K.: A simplified PDF parameterization of subgrid-scale clouds and turbulence for cloud-resolving models, J. Adv. Model. Earth Sy., 5, 195–211, https://doi.org/10.1002/jame.20018, 2013.

Bogenschutz, P. A., Tang, S., Caldwell, P. M., Xie, S., Lin, W., and Chen, Y.-S.: The E3SM version 1 single-column model, Geosci. Model Dev., 13, 4443–4458, https://doi.org/10.5194/gmd-13-4443-2020, 2020.

Böing, S. J., Jonker, H. J. J., Siebesma, A. P., and Grabowski, W. W.: Influence of the subcloud layer on the development of a deep convective ensemble, J. Atmos. Sci., 69, 2682–2698, https://doi.org/10.1175/JAS-D-11-0317.1, 2012.

Brient, F., Couvreux, F., Villefranque, N., Rio, C., and Honnert, R.: Object-Oriented Identification of Coherent Structures in Large Eddy Simulations: Importance of Downdrafts in Stratocumulus, Geophys. Res. Lett., 46, 2854–2864, https://doi.org/10.1029/2018GL081499, 2019.

Brown, A. R., Cederwall, R. T., Chlond, A., Duynkerke, P. G., Golaz, J. C., Khairoutdinov, M., Lewellen, D. C., Lock, A. P., MacVean, M. K., Moeng, C. H., Neggers, R. A. J., Siebesma, A. P., and Stevens, B.: Large-eddy simulation of the diurnal cycle of shallow cumulus convection over land, Q. J. Roy. Meteor. Soc., 128, 1075–1093, https://doi.org/10.1256/003590002320373210, 2002.

Caldwell, P. M., Terai, C. R., Hillman, B., Keen, N. D., Bogenschutz, P., Lin, W., Beydoun, H., Taylor, M., Bertagna, L., Bradley, A. M., Clevenger, T. C., Donahue, A. S., Eldred, C., Foucar, J., Golaz, J. C., Guba, O., Jacob, R., Johnson, J., Krishna, J., Liu, W., Pressel, K., Salinger, A. G., Singh, B., Steyer, A., Ullrich, P., Wu, D., Yuan, X., Shpund, J., Ma, H. Y., and Zender, C. S.: Convection-Permitting Simulations With the E3SM Global Atmosphere Model, J. Adv. Model. Earth Sy., 13, e2021MS002544, https://doi.org/10.1029/2021MS002544, 2021.

Cheng, A. and Xu, K. M.: Simulation of shallow cumuli and their transition to deep convective clouds by cloud-resolving models with different third-order turbulence closures, Q. J. Roy. Meteor. Soc., 132, 359–382, https://doi.org/10.1256/qj.05.29, 2006.

Cheng, A. and Xu, K. M.: Simulation of boundary-layer cumulus and stratocumulus clouds using a cloud-resolving model with low- and third-order turbulence closures, J. Meteorol. Soc. Japan, 86A, 67–86, https://doi.org/10.2151/jmsj.86A.67, 2008.

Chinita, M.: Code for SHOC+MF_v1.0, Zenodo [code], https://doi.org/10.5281/zenodo.7011628, 2022a.

Chinita, M.: Data for SHOC+MF_v1.0, Zenodo [data set], https://doi.org/10.5281/zenodo.7011652, 2022b.

Chinita, M. J., Matheou, G., and Teixeira, J.: A joint probability density-based decomposition of turbulence in the atmospheric boundary layer, Mon. Weather Rev., 146, 503–523, https://doi.org/10.1175/MWR-D-17-0166.1, 2018.

Chinita, M. J., Matheou, G., and Miranda, P. M. A.: Large-eddy simulation of very stable boundary layers. Part I: Modeling methodology, Q. J. Roy. Meteor. Soc., 148, 1805–1823, https://doi.org/10.1002/qj.4279, 2022a.

Chinita, M. J., Matheou, G., and Miranda, P. M. A.: Large-eddy simulation of very stable boundary layers. Part II: Length scales and anisotropy in stratified atmospheric turbulence, Q. J. Roy. Meteor. Soc., 148, 1824–1839, https://doi.org/10.1002/qj.4280, 2022b.

Chung, D. and Matheou, G.: Large-eddy simulation of stratified turbulence. Part I: A vortex-based subgrid-scale model, J. Atmos. Sci., 71, 1863–1879, https://doi.org/10.1175/JAS-D-13-0126.1, 2014.

Chung, D., Matheou, G., and Teixeira, J.: Steady-state large-eddy simulations to study the stratocumulus to shallow cumulus cloud transition, J. Atmos. Sci., 69, 3264–3276, https://doi.org/10.1175/JAS-D-11-0256.1, 2012.

Cohen, Y., Lopez-Gomez, I., Jaruga, A., He, J., Kaul, C. M., and Schneider, T.: Unified Entrainment and Detrainment Closures for Extended Eddy-Diffusivity Mass-Flux Schemes, J. Adv. Model. Earth Sy., 12, e2020MS002162, https://doi.org/10.1029/2020MS002162, 2020.

Couvreux, F., Hourdin, F., and Rio, C.: Resolved versus parametrized boundary-layer plumes. Part I: A parametrization-oriented conditional sampling in large-eddy simulations, Bound.-Lay. Meteorol., 134, 441–458, https://doi.org/10.1007/s10546-009-9456-5, 2010.

Couvreux, F., Bazile, E., Rodier, Q., Maronga, B., Matheou, G., Chinita, M. J., Edwards, J., van Stratum, B. J. H., van Heerwaarden, C. C., Huang, J., Moene, A. F., Cheng, A., Fuka, V., Basu, S., Bou-Zeid, E., Canut, G., and Vignon, E.: Intercomparison of Large-Eddy Simulations of the Antarctic Boundary Layer for Very Stable Stratification, Bound.-Lay. Meteorol., 176, 369–400, https://doi.org/10.1007/s10546-020-00539-4, 2020.

Deardorff, J. W.: The Counter-Gradient Heat Flux in the Lower Atmosphere and in the Laboratory, J. Atmos. Sci., 23, 503–506, https://doi.org/10.1175/1520-0469(1966)023<0503:tcghfi>2.0.co;2, 1966.

E3SM Project, DOE: Energy Exascale Earth System Model v2.0. (Computer Software), DOE [software], https://doi.org/10.11578/E3SM/dc.20210927.1, 2021.

Firl, G. J. and Randall, D. A.: Fitting and analyzing les using multiple trivariate Gaussians, J. Atmos. Sci., 72, 1094–1116, https://doi.org/10.1175/JAS-D-14-0192.1, 2015.

Fitch, A. C.: An improved double-Gaussian closure for the subgrid vertical velocity probability distribution function, J. Atmos. Sci., 76, 285–304, https://doi.org/10.1175/JAS-D-18-0149.1, 2019.

Golaz, J. C., Larson, V. E., and Cotton, W. R.: A PDF-based model for boundary layer clouds. Part I: Method and model description, J. Atmos. Sci., 59, 3540–3551, https://doi.org/10.1175/1520-0469(2002)059<3540:APBMFB>2.0.CO;2, 2002.

Golaz, J. C., Caldwell, P. M., Van Roekel, L. P., Petersen, M. R., Tang, Q., Wolfe, J. D., Abeshu, G., Anantharaj, V., Asay-Davis, X. S., Bader, D. C., Baldwin, S. A., Bisht, G., Bogenschutz, P. A., Branstetter, M., Brunke, M. A., Brus, S. R., Burrows, S. M., Cameron-Smith, P. J., Donahue, A. S., Deakin, M., Easter, R. C., Evans, K. J., Feng, Y., Flanner, M., Foucar, J. G., Fyke, J. G., Griffin, B. M., Hannay, C., Harrop, B. E., Hoffman, M. J., Hunke, E. C., Jacob, R. L., Jacobsen, D. W., Jeffery, N., Jones, P. W., Keen, N. D., Klein, S. A., Larson, V. E., Leung, L. R., Li, H. Y., Lin, W., Lipscomb, W. H., Ma, P. L., Mahajan, S., Maltrud, M. E., Mametjanov, A., McClean, J. L., McCoy, R. B., Neale, R. B., Price, S. F., Qian, Y., Rasch, P. J., Reeves Eyre, J. E. J., Riley, W. J., Ringler, T. D., Roberts, A. F., Roesler, E. L., Salinger, A. G., Shaheen, Z., Shi, X., Singh, B., Tang, J., Taylor, M. A., Thornton, P. E., Turner, A. K., Veneziani, M., Wan, H., Wang, H., Wang, S., Williams, D. N., Wolfram, P. J., Worley, P. H., Xie, S., Yang, Y., Yoon, J. H., Zelinka, M. D., Zender, C. S., Zeng, X., Zhang, C., Zhang, K., Zhang, Y., Zheng, X., Zhou, T., and Zhu, Q.: The DOE E3SM Coupled Model Version 1: Overview and Evaluation at Standard Resolution, J. Adv. Model. Earth Sy., 11, 2089–2129, https://doi.org/10.1029/2018MS001603, 2019.

Han, J. and Bretherton, C. S.: TKE-based moist eddy-diffusivity mass-flux (EDMF) parameterization for vertical turbulent mixing, Weather Forecast., 34, 869–886, https://doi.org/10.1175/WAF-D-18-0146.1, 2019.

Han, J. and Pan, H. L.: Revision of convection and vertical diffusion schemes in the NCEP Global Forecast System, Weather Forecast., 26, 520–533, https://doi.org/10.1175/WAF-D-10-05038.1, 2011.

Han, J., Witek, M. L., Teixeira, J., Sun, R., Pan, H. L., Fletcher, J. K., and Bretherton, C. S.: Implementation in the NCEP GFS of a hybrid eddy-diffusivity mass-flux (EDMF) boundary layer parameterization with dissipative heating and modified stable boundary layer mixing, Weather Forecast., 31, 341–352, https://doi.org/10.1175/WAF-D-15-0053.1, 2016.

Holtslag, A. A. M. and Boville, B. A.: Local versus nonlocal boundary-layer diffusion in a global climate model, J. Climate, 6, 1825–1842, https://doi.org/10.1175/1520-0442(1993)006<1825:LVNBLD>2.0.CO;2, 1993.

Holtslag, A. A. M. and Moeng, C.-H.: Eddy diffusivity and countergradient transport in the convective atmospheric boundary layer, J. Atmos. Sci., 48, 1690–1698, https://doi.org/10.1175/1520-0469(1991)048<1690:EDACTI>2.0.CO;2, 1991.

Kohler, M., Ahlgrimm, M., and Beljaars, A.: Unified treatment of dry convective and stratocumulus-topped boundary layers in the ECMWF model, Q. J. Roy. Meteor. Soc., 137, 43–57, https://doi.org/10.1002/qj.713, 2011.

Kurowski, M. J., Suselj, K., and Grabowski, W. W.: Is Shallow Convection Sensitive to Environmental Heterogeneities?, Geophys. Res. Lett., 46, 1785–1793, https://doi.org/10.1029/2018GL080847, 2019a.

Kurowski, M. J., Thrastarson, H. T., Suselj, K., and Teixeira, J.: Towards unifying the planetary boundary layer and shallow convection in CAM5 with the eddy-diffusivity/mass-flux approach, Atmosphere, 10, 484, https://doi.org/10.3390/atmos10090484, 2019b.

Lamaakel, O. and Matheou, G.: Galilean invariance of shallow cumulus convection large-eddy simulations, J. Comput. Phys., 427, 110012, https://doi.org/10.1016/j.jcp.2020.110012, 2021.

Matheou, G.: Turbulence Structure in a Stratocumulus Cloud, Atmosphere, 9, 392, https://doi.org/10.3390/atmos9100392, 2018.

Matheou, G. and Chung, D.: Large-eddy simulation of stratified turbulence. Part II: Application of the stretched-vortex model to the atmospheric boundary layer, J. Atmos. Sci., 71, 4439–4460, https://doi.org/10.1175/JAS-D-13-0306.1, 2014.

Matheou, G. and Teixeira, J.: Sensitivity to physical and numerical aspects of large-eddy simulation of stratocumulus, Mon. Weather Rev., 147, 2621–2639, https://doi.org/10.1175/MWR-D-18-0294.1, 2019.

Matheou, G., Chung, D., and Teixeira, J.: On the synergy between numerics and subgrid scale modeling in LES of stratified flows: Grid convergence of a stratocumulus-topped boundary layer., in: Eighth Int. Symp. on Stratified Flows, Vol. 1, San Diego, CA, NASA Jet Propulsion Laboratory, http://hdl.handle.net/2014/46171 (last access: 20 March 2023), 2016.

Morinishi, Y., Lund, T. S., Vasilyev, O. V, and Moin, P.: Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow, J. Comput. Phys., 143, 90–124, https://doi.org/10.1006/jcph.1998.5962, 1998.

Neggers, R. A. J.: A dual mass flux framework for boundary layer convection. Part II: Clouds, J. Atmos. Sci., 66, 1489–1506, https://doi.org/10.1175/2008JAS2636.1, 2009.

Olson, J. B., Kenyon, J. S., Angevine, W. A., Brown, J. M., Pagowski, M., and Sušelj, K.: A Description of the MYNN-EDMF Scheme and the Coupling to Other Components in WRF–ARW, NOAA Technical Memorandum OAR GSD-61, https://doi.org/10.25923/n9wm-be49, 2019.

Richtmyer, R. and Morton, K.: Difference methods for initial value problems, 2nd edn., Wiley-Interscience, New York, 405 pp., 1967.

Rio, C. and Hourdin, F.: A thermal plume model for the convective boundary layer: Representation of cumulus clouds, J. Atmos. Sci., 65, 407–425, https://doi.org/10.1175/2007JAS2256.1, 2008.

Romps, D. M. and Kuang, Z.: Nature versus nurture in shallow convection, J. Atmos. Sci., 67, 1655–1666, https://doi.org/10.1175/2009JAS3307.1, 2010.

Schneider, T., Teixeira, J., Bretherton, C. S., Brient, F., Pressel, K. G., Schär, C., and Siebesma, A. P.: Climate goals and computing the future of clouds, Nat. Clim. Change, 7, 3–5, https://doi.org/10.1038/nclimate3190, 2017.

Sherwood, S. C., Bony, S., and Dufresne, J. L.: Spread in model climate sensitivity traced to atmospheric convective mixing, Nature, 505, 37–42, https://doi.org/10.1038/nature12829, 2014.

Siebesma, A., Bretherton, C. S., Brown, A., Chlond, A., Cuxart, J., Duynkerke, P. G., Jiang, H., Khairoutdinov, M., Lewellen, D., Moeng, C. H., Sanchez, E., Stevens, B., and Stevens, D. E.: A large eddy simulation intercomparison study of shallow cumulus convection, J. Atmos. Sci., 60, 1201–1219, https://doi.org/10.1175/1520-0469(2003)60<1201:ALESIS>2.0.CO;2, 2003.

Siebesma, A. P. and Cuijpers, J. W. M.: Evaluation of parametric assumptions for shallow cumulus convection, J. Atmos. Sci., 52, 650–666, https://doi.org/10.1175/1520-0469(1995)052<0650:EOPAFS>2.0.CO;2, 1995.

Siebesma, A. P., Soares, P. M. M., and Teixeira, J.: A combined eddy-diffusivity mass-flux approach for the convective boundary layer, J. Atmos. Sci., 64, 1230–1248, https://doi.org/10.1175/JAS3888.1, 2007.

Siebesma, P. and Teixeira, J.: An advection-diffusion scheme for the convective boundary layer: Description and 1D-results, in: Proc. 14th Symp. on Boundary Layers and Turbulence, Aspen, CO, Amer. Meteor. Soc., 133–136, 2000.

Soares, P. M. M., Miranda, P. M. A., Siebesma, A. P., and Teixeira, J.: An eddy-diffusivity/mass-flux parametrization for dry and shallow cumulus convection, Q. J. Roy. Meteor. Soc., 130, 3365–3383, https://doi.org/10.1256/qj.03.223, 2004.

Stevens, B.: Quasi-steady analysis of a PBL model with an Eddy-diffusivity profile and nonlocal fluxes, Mon. Weather Rev., 128, 824–836, https://doi.org/10.1175/1520-0493(2000)128<0824:QSAOAP>2.0.CO;2, 2000.

Suselj, K., Teixeira, J., and Chung, D.: A unified model for moist convective boundary layers based on a stochastic eddy-diffusivity/mass-flux parameterization, J. Atmos. Sci., 70, 1929–1953, https://doi.org/10.1175/JAS-D-12-0106.1, 2013.

Suselj, K., Hogan, T. F., and Teixeira, J.: Implementation of a stochastic eddy-diffusivity/mass-flux parameterization into the Navy Global environmental model, Weather Forecast., 29, 1374–1390, https://doi.org/10.1175/WAF-D-14-00043.1, 2014.

Suselj, K., Kurowski, M. J., and Teixeira, J.: A unified eddy-diffusivity/mass-flux approach for modeling atmospheric convection, J. Atmos. Sci., 76, 2505–2537, https://doi.org/10.1175/JAS-D-18-0239.1, 2019a.

Suselj, K., Kurowski, M. J., and Teixeira, J.: On the factors controlling the development of shallow convection in eddy-diffusivity/mass-flux models, J. Atmos. Sci., 76, 433–456, https://doi.org/10.1175/JAS-D-18-0121.1, 2019b.

Suselj, K., Teixeira, J., Kurowski, M. J., and Molod, A.: Improving the representation of subtropical boundary layer clouds in the NASA GEOS model with the eddy-diffusivity/mass-flux parameterization, Mon. Weather Rev., 149, 793–809, https://doi.org/10.1175/MWR-D-20-0183.1, 2021.

Takahashi, H., Luo, Z. J., and Stephens, G.: Revisiting the Entrainment Relationship of Convective Plumes: A Perspective From Global Observations, Geophys. Res. Lett., 48, e2020GL092349, https://doi.org/10.1029/2020GL092349, 2021.

Tan, Z., Kaul, C. M., Pressel, K. G., Cohen, Y., Schneider, T., and Teixeira, J.: An Extended Eddy-Diffusivity Mass-Flux Scheme for Unified Representation of Subgrid-Scale Turbulence and Convection, J. Adv. Model. Earth Sy., 10, 770–800, https://doi.org/10.1002/2017MS001162, 2018.

Teixeira, J. and Cheinet, S.: A simple mixing length formulation for the eddy-diffusivity parameterization of dry convection, Bound.-Lay. Meteorol., 110, 435–453, https://doi.org/10.1023/B:BOUN.0000007230.96303.0d, 2004.

Teixeira, J. and Siebesma, P.: A mass flux/K-diffusion approach to the parameterization of the convective boundary layer: Global model results, in: Proc. 14th Symp. on Boundary Layers and Turbulence, Aspen, CO, Amer. Meteor. Soc., 231–234, 2000.

Teixeira, J., Ferreira, J. P., Miranda, P. M. A., Haack, T., Doyle, J., Siebsema, A. P., and Salgado, R.: A new mixing-length formulation for the parameterization of dry convection: Implementation and evaluation in mesoscale model, Mon. Weather Rev., 132, 2698–2707, https://doi.org/10.1175/MWR2808.1, 2004.

Teixeira, J., Stevens, B., Bretherton, C. S., Cederwall, R., Doyle, J. D., Golaz, J. C., Holtslag, A. A. M., Klein, S. A., Lundquist, J. K., Randall, D. A., Siebesma, A. P., and Soares, P. M. M.: Parameterization of the atmospheric boundary layer: A view from just above the inversion, B. Am. Meteorol. Soc., 89, 453–458, https://doi.org/10.1175/BAMS-89-4-453, 2008.

Tiedtke, M.: A Comprehensive Mass Flux Scheme for Cumulus Parameterization in Large-Scale Models, Mon. Weather Rev., 117, 1779–1800, 1989.

Witek, M. L., Teixeira, J., and Matheou, G.: An Integrated TKE-based eddy diffusivity/mass flux boundary layer closure for the dry convective boundary layer, J. Atmos. Sci., 68, 1526–1540, https://doi.org/10.1175/2011JAS3548.1, 2011.

Witte, M. K., Herrington, A., Teixeira, J., Kurowski, M., Chinita, M. J., Storer, R. L., Suselj, K., Matheou, G., and Bacmeister, J.: Augmenting the double-Gaussian representation of atmospheric turbulence and convection via a coupled stochastic multi-plume mass flux scheme, Mon. Weather Rev., 150, 2339–2355, https://doi.org/10.1175/MWR-D-21-0215.1, 2022.

Wu, E., Yang, H., Kleissl, J., Suselj, K., Kurowski, M. J., and Teixeira, J.: On the parameterization of convective downdrafts for marine stratocumulus clouds, Mon. Weather Rev., 148, 1931–1950, https://doi.org/10.1175/MWR-D-19-0292.1, 2020.

Yoshimura, H., Mizuta, R., and Murakami, H.: A spectral cumulus parameterization scheme interpolating between two convective updrafts with semi-lagrangian calculation of transport by compensatory subsidence, Mon. Weather Rev., 143, 597–621, https://doi.org/10.1175/MWR-D-14-00068.1, 2015.