Numerics and subgrid-scale modeling in large eddy simulations of stratocumulus clouds.

Stratocumulus clouds are the most common type of boundary layer cloud; their radiative effects strongly modulate climate. Large eddy simulations (LES) of stratocumulus clouds often struggle to maintain fidelity to observations because of the sharp gradients occurring at the entrainment interfacial layer at the cloud top. The challenge posed to LES by stratocumulus clouds is evident in the wide range of solutions found in the LES intercomparison based on the DYCOMS-II field campaign, where simulated liquid water paths for identical initial and boundary conditions varied by a factor of nearly 12. Here we revisit the DYCOMS-II RF01 case and show that the wide range of previous LES results can be realized in a single LES code by varying only the numerical treatment of the equations of motion and the nature of subgrid-scale (SGS) closures. The simulations that maintain the greatest fidelity to DYCOMS-II observations are identified. The results show that using weighted essentially non-oscillatory (WENO) numerics for all resolved advective terms and no explicit SGS closure consistently produces the highest-fidelity simulations. This suggests that the numerical dissipation inherent in WENO schemes functions as a high-quality, implicit SGS closure for this stratocumulus case. Conversely, using oscillatory centered difference numerical schemes for momentum advection, WENO numerics for scalars, and explicitly modeled SGS fluxes consistently produces the lowest-fidelity simulations. We attribute this to the production of anomalously large SGS fluxes near the cloud tops through the interaction of numerical error in the momentum field with the scalar SGS model.

Introduction

The representation of boundary layer clouds and their climate feedbacks persists as a major source of uncertainty in GCM predictions of climate sensitivity [e.g., Cess et al., 1990, 1996; Bony and Dufresne, 2005; Webb et al., 2006; Dufresne and Bony, 2008; Vial et al., 2013; Brient et al., 2016]. Stratocumulus are the most prevalent boundary layer cloud type [Wood, 2012], and they play a central role in controlling the global energy budget because of their large area coverage [Stephens and Greenwald, 1991; Hartmann et al., 1992]. So important are stratocumulus to the global energy budget that a relatively small change in their area coverage has an effect on the global energy budget similar in magnitude to that of the anthropogenic emissions of greenhouse gases [Hartmann and Short, 1980; Randall et al., 1984; Slingo, 1990]. For example, Slingo [1990] estimates that a 15%–20% change in the coverage of boundary layer clouds, of which stratocumulus are the dominant component, could offset the radiative effects of doubling the atmospheric concentration.

The importance of stratocumulus clouds to climate has motivated studies of their dynamics and microphysics, both observational [e.g., Albrecht et al., 1988; Lenschow et al., 1988; Stevens et al., 2003a; Haman et al., 2007; Kalmus et al., 2014] and numerical [e.g., Moeng et al., 1996; Chlond and Wolkau, 2000; Moeng, 2000; Moeng et al., 2005; Stevens et al., 2005; Kirkpatrick et al., 2006; Yamaguchi and Randall, 2008; Kurowski et al., 2009a, 2009b; Ackerman et al., 2009; Mellado et al., 2010; Yamaguchi and Feingold, 2012; Blossey et al., 2013]. From these studies, a clear conceptual picture of stratocumulus dynamics has emerged. In this picture, stratocumulus clouds appear as the saturated portion of a turbulent boundary layer driven from above by longwave radiative cooling at the tops of the clouds. The radiatively driven turbulence entrains mass from the overlying statically stable and relatively dry lower troposphere into the boundary layer. This process of radiative cooling and entrainment establishes a layer of large gradients in thermodynamic quantities such as temperature and humidity called the entrainment interfacial layer (EIL). The persistence of stratocumulus boundary layers requires that the entrainment of dry lower‐tropospheric air across the EIL does not dry the cloud layer so much that cloud‐top radiative cooling weakens and turbulent motions in the cloud layer decouple from the near‐surface layer with its moisture source. Therefore, accurate representation of the dynamical and thermodynamical processes within the EIL is critical for sustaining stratocumulus in numerical models.

In recent decades, increasing high‐performance computing capabilities have made LES an important tool for understanding the dynamics and microphysics of many cloud types [e.g., Matheou et al., 2011; Matheou and Chung, 2014], a trend that is likely to continue [e.g., Schneider et al., 2017]. This is particularly true for stratocumulus clouds. Nonetheless, faithful simulation of stratocumulus remains a challenge for LES codes [e.g., Stevens et al., 2005; Yamaguchi and Feingold, 2012]. Much of the difficulty arises from feedbacks between the strength of the EIL, turbulent mixing, radiative cooling, and the amount of liquid water in the clouds. These feedbacks can amplify numerical and modeling errors, making LES sensitive to details such as thermodynamic formulation [e.g., Xiao et al., 2015], subgrid‐scale (SGS) models [e.g., Stevens et al., 2005; Kirkpatrick et al., 2006; Kurowski et al., 2009a], and grid resolution [Pedersen et al., 2016]. However, a systematic study of the sensitivity of stratocumulus LES to numerical formulation is lacking. The importance of the numerical formulation in LES is well‐understood already in the context of idealized turbulent flows [e.g., Ghosal, 1996; Chow and Moin, 2003], which lack the strong feedbacks associated with stratocumulus‐topped boundary layers. The numerical formulation can be expected to be at least as important for stratocumulus clouds. Indeed, Stevens et al. [2005] cite the numerical formulation and SGS closures as significant contributors to the large spread in an LES intercomparison study of a stratocumulus boundary layer.

Modern LES codes typically incorporate a modular design, in which multiple numerical schemes that may have different numerical properties are available alongside multiple SGS closures [e.g., Heus et al., 2010; Pressel et al., 2015]. In PyCLES [Pressel et al., 2015], the code used here, both weighted essentially non‐oscillatory (WENO) [Liu et al., 1994; Jiang and Shu, 1996; Balsara and Shu, 2000] and central difference [Wicker and Skamarock, 2002] schemes are available for the discretization of momentum and scalar transport. Conceptually, the fundamental difference between centered and WENO schemes is that the numerical error of the former is dispersive, while the numerical error of the latter is largely dissipative. The dispersive numerical error associated with centered schemes produces grid‐scale oscillations in the numerical solution [e.g., Matheou and Dimotakis, 2016], leading to overestimation of small‐scale gradients. On the other hand, the numerical dissipation inherent in WENO schemes tends to smooth small‐scale gradients, though, as their name implies, they are not entirely free from spurious oscillations. In the remainder of this article, we refer to centered schemes as “oscillatory” and WENO schemes as “dissipative” to emphasize this fundamental difference between the two classes of schemes. This terminology is not meant to imply that all nominally dissipative numerical schemes will produce results comparable to those obtained here, given that WENO schemes are designed to remain essentially non‐oscillatory even in the vicinity of discontinuities in the flow field.

As SGS models are designed to represent the dissipative nature of unresolved turbulent processes, it is reasonable to conjecture that dissipative and oscillatory errors may interact with SGS models in different ways. This sort of interaction is alluded to in the work of Stevens et al. [2005] and Savic‐Jovcic and Stevens [2010], where configuring an LES code to rely on numerical dissipation in the advection scheme as the sole source of SGS scalar transport allowed it to produce high‐fidelity simulations. LES configured in this way were shown by Stevens et al. [2005] to significantly outperform LES with oscillatory treatment of momentum advection and dissipative treatment of scalar advection. Moreover, others have shown that realistic LES of neutral, convective, and stratocumulus boundary layers can be performed using dissipative numerical schemes for both momentum and scalar transport, with or without explicit SGS closures [e.g., Margolin et al., 1999; Brown and MacVean, 2000; Stevens et al., 2002; de Szoeke and Bretherton, 2004; Stevens et al., 2005; Kurowski et al., 2009a; Ackerman et al., 2009; Pedersen et al., 2016]. We perform a set of LES experiments with different combinations of momentum and scalar advection schemes, with and without explicit SGS closures, to understand how numerical error in the resolved transport scheme and its interaction with SGS closures determines the fidelity of LES of stratocumulus.

We specifically focus on the nonprecipitating stratocumulus case described in the Stevens et al. [2005] LES intercomparison study, which sought to verify LES against observations. The intercomparison exhibits a wide spread of LES solutions, with liquid water path and cloud fraction varying by factors of roughly 12 and 5, respectively. This spread in results suggests strong feedbacks between entrainment, liquid water path, and radiation. However, some simplification is provided by the fact the DYCOMS‐II RF01 case is nonprecipitating. This eliminates the need to consider sensitivities to microphysical model formulation. In particular, any potentially confounding interactions between numerical error and the SGS closure with highly uncertain microphysical models are excluded.

In section 2, we describe the formulation of the LES code, the DYCOMS‐II RF01 case, and the numerical and SGS configurations used in each of the numerical simulations. In section 3, we show the results of the numerical simulations and assess their fidelity to observations. In section 4, we discuss the results, giving particular attention to understanding how the interaction of numerical error with the SGS closure can degrade simulation accuracy. In section 5, we conclude by providing guidance for best practices when configuring LES for the simulation of stratocumulus. We also discuss the potential extension of these findings to LES of other cloud types.

Methods and Simulations

LES Description

The simulations described in this study are performed using PyCLES, a recently developed Python‐based LES code described in detail in Pressel et al. [2015]. PyCLES solves an energetically consistent form of the anelastic equations of motion using total water specific humidity q and moist specific entropy s as prognostic thermodynamic variables [Pauluis, 2008]. PyCLES has multiple options for discretizing the equations of motion, including standard central difference schemes and more sophisticated WENO schemes.

The transport schemes used in this study are implemented almost exactly as described in Pressel et al. [2015], with one exception. The staggered arrangement of variables on an Arakawa C‐grid [Arakawa and Lamb, 1977] requires an additional interpolation of the flux velocity in the momentum transport schemes. In Pressel et al. [2015], an interpolation scheme with an order of accuracy commensurate with the order of the other components of the advection scheme was used. In the current work, however, second‐order interpolation of the flux velocity is used in conjunction with all advection schemes. This makes the centered schemes used here identical to those described by Wicker and Skamarock [2002], which are widely used in atmospheric simulations [e.g., Skamarock and Klemp, 2008; Heus et al., 2010]. This change to the WENO schemes described in Pressel et al. [2015] does not reduce their order of accuracy, as in either case the WENO schemes, like the centered schemes [Skamarock and Klemp, 2008], are second‐order accurate for nonlinear problems on staggered grids. Indeed, sensitivity tests (not shown) indicate that the results are not sensitive to this detail of the numerical schemes. While all of the nominally higher‐order schemes exhibit second‐order convergence, the nominal orders of the schemes are indicative of their relative accuracy at a given mesh resolution. In other words, the nominally higher‐order schemes have lower total error all other things being equal.

All simulations employ the third‐order strong stability preserving time stepping scheme of Shu and Osher [1988]. The combination of time stepping and scalar transport schemes used here is not positivity preserving, thus clipping is employed to guarantee positivity of q in thermodynamic computations. The results reported here use a nonconservative approach to clipping, which leads to a positive moisture bias in the simulations where clipping is active. Clipping is rarely if ever active in simulations using WENO schemes but is more active when using centered schemes. Sensitivity tests with conservative clipping (not shown) indicate that the lack of conservation does not impact our conclusions. The time step is adjusted dynamically to maintain a Courant‐Friedrichs‐Lewy (CFL) number of approximately 0.3. The anelastic continuity equation is enforced using a predictor‐corrector method as in Pressel et al. [2015], which is second‐order accurate in space. The use of a second‐order accurate predictor‐corrector method is consistent with the order of accuracy of the momentum and scalar transport schemes, which only guarantee greater than second‐order accuracy for simple linear advection [e.g., Skamarock and Klemp, 2008].

In all simulations using explicit closure of SGS fluxes, a Smagorinsky‐Lilly model [Smagorinsky, 1958, 1963] with stability correction [Lilly, 1962] is used as described in Pressel et al. [2015]. The effect of the stability correction is to shrink the eddy viscosity toward zero as the stratification becomes increasingly stable. The Smagorinsky constant c is specified to be 0.17 in all simulations. The eddy viscosity diagnosed by the Smagorinsky‐Lilly model is used to compute the eddy diffusivity by assuming a turbulent Prandtl number . Simulations with a TKE‐based closure led to similar results.

Case Description

The initial condition, large‐scale forcing, radiative forcing, and surface boundary conditions of the DYCOMS‐II RF01 case are as described in Stevens et al. [2005]. Because PyCLES uses q and s as prognostic thermodynamic variables, special care must be taken to ensure that the initial condition prescribed through this set of variables is consistent with an initial condition given in terms of q and the liquid water potential temperature θ as in Stevens et al. [2005]. To do this, we assume that θ is given by as in Tripoli and Cotton [1981], where θ is the dry potential temperature, L is the latent heat of vaporization, T is the temperature, c is the specific heat of dry air at constant pressure, and q is the specific humidity of liquid water. The thermodynamic constants used to compute θ from equation (1) are listed in Table 1 and are identical to those given in Stevens et al. [2005]. A saturation adjustment procedure is performed to determine initial values of T and q. Once T and q are known, the specific entropy s is computed using where the dry air specific entropy s and water vapor specific entropy s are given by and

Table 1

Thermodynamic Constants Used in Computation of θl, Taken From Stevens et al. [2005].

Parameter	Value
Gas constant, dry air	R_d=287.0 J kg^−1 K^−1
Specific heat dry air at constant pressure	c_pd=1015.0 J kg^−1 K^−1
Latent heat of vaporization	L_v=2.47×10^6 J kg^−1

The partial pressures of dry air and water vapor, p and p, are determined as in Pressel et al. [2015]. The thermodynamic constants used to compute s are identical to those given in Pressel et al. [2015] and are reproduced here in Table 2. To perturb the initial conditions, fluctuations in θ with maximum amplitude of 0.1 K are drawn from a uniform distribution and added at every grid point at or below 200 m altitude. To ensure thermodynamic consistency, the fluctuations are applied prior to the saturation adjustment from which T and q are determined.

Table 2

Thermodynamic Constants Used in Computation of s, Taken From Pressel et al. [2015]

Parameter	Value
Gas constant, dry air	R_d=287.1 J kg^−1 K^−1
Gas constant, water vapor	R_v=461.5 J kg^−1 K^−1
Specific heat, dry air	c_pd=1004 J kg^−1 K^−1
Specific heat, water vapor	c_pv=1859 J kg^−1 K^−1
Standard temperature	T∼=298.15 K
Standard pressure	p∼=10^5 Pa
Standard entropy, dry air	s∼_d=6699.1 J kg^−1 K^−1
Standard entropy, water vapor	s∼_v=10450.8 J kg^−1 K^−1
Latent heat of vaporization	L_v=2.501×10^6 J kg^−1

The DYCOMS‐II RF01 case specification includes large‐scale subsidence forcing and an empirical representation of radiative transfer [Stevens et al., 2005]. The subsidence forcing is applied directly to the prognostic thermodynamic variables and to the horizontal momentum components in advective form using first‐order upwind finite differencing. The radiative heating rate determined from the empirical radiation scheme is applied as an external moist entropy source term, as described in Pressel et al. [2015].

Experiments

Our goal is to understand how the numerical errors produced by various forms of scalar and momentum advection interact with SGS closures to control the fidelity of LES simulations of stratocumulus. To this end, we construct a set of 17 LES performed using various discretizations for scalar and momentum advection, with and without (when possible) SGS closures. The LES can be logically grouped into the four sets described below, which we refer to as experiments. The experiments are summarized in Table 3, and the results from each set of experiments are plotted with the line labels given in the table.

Table 3

Experiments Included in This Study

Experiment	Case Name	Momentum Advection	Scalar Advection	Momentum SGS	Scalar SGS
Mixed‐SGS	25MS	2nd Central	5th WENO	Yes	Yes
	45MS	4th Central	5th WENO	Yes	Yes
	65MS	6th Central	5th WENO	Yes	Yes
Paired‐SGS	22MS	2nd Central	2nd Central	Yes	Yes
	44MS	4th Central	4th Central	Yes	Yes
	66MS	6th Central	6th Central	Yes	Yes
	55MS	5th WENO	5th WENO	Yes	Yes
	77MS	7th WENO	7th WENO	Yes	Yes
	99MS	9th WENO	9th WENO	Yes	Yes
Mixed‐NSGS	25MN	2nd Central	5th WENO	Yes	No
	45MN	4th Central	5th WENO	Yes	No
	65MN	6th Central	5th WENO	Yes	No
Paired‐NSGS	55NN	5th WENO	5th WENO	No	No
	77NN	7th WENO	7th WENO	No	No
	99NN	9th WENO	9th WENO	No	No

All simulations are performed on a domain extending 3360 m in the horizontal directions and 1500 m in the vertical direction, with 35 m horizontal and 5 m vertical resolution. This domain size and resolution is consistent with the guidelines for the Stevens et al. [2005] LES intercomparison.

Mixed‐SGS

The first experiment (Mixed‐SGS) uses the common LES configuration that pairs an oscillatory advection scheme for momentum with dissipative advection for scalars. Explicit closures are used for SGS transport of both momentum and scalars. The combination of oscillatory numerics for momentum with a dissipative scheme for scalars has become common in LES of atmospheric flows, because in some cases it is possible to conserve kinetic energy using oscillatory numerics [e.g., Harlow and Welch, 1965; Morinishi et al., 1998], while dissipative numerics for scalars avoid spurious scalar oscillations [e.g., Matheou and Chung, 2014]. Here the use of an explicit SGS model for momentum is required to prevent kinetic energy from accumulating near the grid scale, which otherwise occurs because the momentum scheme lacks sufficient numerical dissipation. Following common practice, we also apply an explicit closure for SGS scalar transport, although it may not be needed to maintain a stable solution because of the numerical dissipation inherent in dissipative schemes.

Paired‐SGS

The second experiment (Paired‐SGS) pairs similar momentum and scalar numerics, meaning that the scalar and momentum scheme are either both dissipative or both oscillatory. Explicit SGS closures are used for both momentum and scalars. Simulations using oscillatory numerics for both momentum and scalar transport along with explicit SGS closures have shown significant success for numerous boundary layer and cloud types [e.g., Matheou et al., 2011; Matheou and Chung, 2014]. So too have simulations of neutral and dry convective boundary layers using dissipative scalar and momentum numerics [Brown and MacVean, 2000]. Similarly, simulations performed with the Distributed Hydrodynamic Aerosol and Radiative Modeling Application (DHARMA) [Stevens and Bretherton, 1996; Stevens et al., 2003a], using dissipative scalar and momentum numerics and SGS closures, were among the best performing simulations in the Stevens et al. [2005] intercomparison study. More recently, Yamaguchi and Feingold [2012] describe high‐fidelity simulations of DYCOMS‐II RF01 using this approach. Their simulations differ from most of those presented in Stevens et al. [2005] and those we consider here by including cloud droplet sedimentation and drizzle processes. Our own tests (not shown) indicate that including these microphysical processes increases the fidelity of LES of this case, consistent with the findings of an LES intercomparison study of the related DYCOMS‐II RF02 case [Ackerman et al., 2009].

Mixed‐NSGS

The third experiment (Mixed‐NSGS) pairs various oscillatory momentum discretizations with a dissipative scalar advection scheme. In contrast to the Mixed‐SGS experiment, the SGS closure is disabled for scalars (but not momentum). A similar approach was taken by Stevens et al. [2005] and Savic‐Jovcic and Stevens [2010] and was found to successfully limit spurious entrainment across the EIL and increase simulation fidelity to observations. The numerical dissipation inherent in the scalar advection scheme serves as an implicit SGS model for scalars. While no scalar SGS model is used for the interior of the flow, an eddy diffusivity closure is used to represent the turbulent fluxes in the surface layer, where all scales of turbulent motion are too small to be resolved explicitly [Stevens et al., 2005; Savic‐Jovcic and Stevens, 2010]. Here, the eddy diffusivity is modeled using a Smagorinsky‐Lilly closure as in the Mixed‐SGS and Paired‐SGS experiments, but only within the surface layer. The eddy diffusivity decreases linearly with height to zero at the top of the surface layer, here defined as the first model level whose height exceeds the vertical grid spacing divided by the von Kármán constant.

Paired‐NSGS

The fourth experiment (Paired‐NSGS) pairs dissipative momentum and scalar discretizations. No explicit SGS model is used for either scalars or momentum. This approach has been used in the simulation of dry and neutrally stratified atmospheric boundary layers [Margolin et al., 1999; Brown and MacVean, 2000] and also for stratocumulus boundary layers [Kurowski et al., 2009a; Pedersen et al., 2016]. This method has come to be known as implicit large eddy simulation, referring to the fact that SGS closures are implicitly provided by the dissipation inherent in the numerical schemes. Implicit LES has been shown to be successful in a broad class of shear flows [e.g., Oran and Boris, 1993; Fureby and Grinstein, 1999; Margolin et al., 2002].

As in the Mixed‐NSGS experiment, an eddy diffusivity closure is used to model turbulent fluxes in the surface layer.

Results

In the results presented below, all time series, for example, of cloud cover and liquid water path (LWP), are sampled every 60 s of simulated time. Profiles of the mean thermodynamic variables and higher‐order vertical velocity statistics are also collected every 60 s, then are averaged over the last two hours of the simulation to produce the plots shown in what follows.

Cloud Cover

Vertical profiles of cloud liquid water q are shown in Figure 1 for each of the four sets of experiments, along with the observed q at four levels within the cloudy layer [Stevens et al., 2005]. The results show a striking sensitivity to the numerical treatment of the equations of motion and to the use of SGS models. The maximum q varies by more than a factor of 3 across all simulations, a range similar to that seen in the Stevens et al. [2005] LES intercomparison. The observed maximum in q occurs just below 800 m, although with such sparse observations it is difficult to know if this actually samples the level of maximum q. Qualitatively, all simulations produce similarly shaped profiles of q, with maximum values occurring near cloud top, and they generally agree on the heights of cloud base and top. Nonetheless, the simulated q profiles show substantial quantitative differences.

Figure 1

Cloud liquid water q for each of the four experiments (see Table 3). Observations from Stevens et al. [2003b, 2005] are indicated with black circles.

The quantitative differences in q among the various cases are also clearly reflected in the time series of cloud fraction and LWP shown in Figures 2 and 3. In these figures, the Stevens et al. [2005] intercomparison mean, interquartile range, and absolute range are shown by the solid black line, dark gray shading, and light gray shading. Additionally, the mean values of LWP and cloud fraction over the last 2 h of each simulation are reported in Table 4. Across all simulations, liquid water path varies by a factor of more than 6, and cloud fraction varies by a factor of almost 2, essentially spanning the range of intercomparison results. While dependable observations of liquid water path are not available for this case, cloud fraction was observed to be near one during the field campaign [Stevens et al., 2005].

Figure 2

Time series of cloud fraction for each of the experiments. The absolute and interquartile ranges of the Stevens et al. [2005] LES intercomparison are shown in light and dark shading, respectively. The solid back line shows the Stevens et al. [2005] intercomparison mean.

Figure 3

Same as in Figure 2 but for liquid water path.

Table 4

Liquid Water Path (LWP) and Cloud Fraction (CF) for Each of the Cases

Experiment	Case Name	CF	LWP (μm)
Mixed‐SGS	25MS	0.54	9.1
	45MS	0.52	9.3
	65MS	0.51	9.9
Paired‐SGS	22MS	0.75	12.4
	44MS	0.97	30.1
	66MS	0.98	33.4
	55MS	0.98	39.3
	77MS	0.99	42.0
	99MS	0.99	45.1
Mixed‐NSGS	25MN	0.80	19.2
	45MN	0.87	24.0
	65MN	0.88	23.3
Paired‐NSGS	55NN	1.0	56.3
	77NN	1.0	55.8
	99NN	1.0	53.6

The simulations of the Mixed‐SGS experiment produce the sparsest cloud fields by far (Figure 2), and their mean profiles of q deviate substantially from observations at all levels (Figure 1). Relative to the results of the Stevens et al. [2005] intercomparison, the LWP predicted by simulations in the Mixed‐SGS experiment systematically fall on the low side of the predicted range (Figure 3). The simulations here show little dependence on the accuracy of the momentum discretization, suggesting that either the numerical errors in the momentum field are not of particular importance in determining the simulated q profile, or that other sources of error dominate.

The simulations in the Paired‐SGS experiment feature much denser cloud fields (Figure 2). Maxima in the q profiles are more comparable to observed values, although q values in the lower portion of the cloud layer underpredict the observations (Figure 1). Both the LWP (Figure 3) and cloud fraction predicted here (with the exception of the 22MS case) fall within their interquartile ranges from the Stevens et al. [2005] intercomparison. The 22MS, 44MS, and 66MS cases in this experiment differ from the 25MS, 45MS, and 65MS cases in the Mixed‐SGS experiment only by the replacement of dissipative scalar numerics with oscillatory scalar numerics (refer to Table 3 for a listing of case names and configurations). Yet this change alone leads to substantial increases in the fidelity of the paired 44MS and 66MS cases over the mixed 45MS and 65MS cases. The 22MS case, which uses second‐order centered numerics, does not show as marked an increase in fidelity as the other two paired oscillatory cases, likely because of the low numerical accuracy of the second‐order schemes. Ignoring the 22MS case, there is only a slight sensitivity of LWP and cloud fraction to the use of oscillatory or dissipative numerics in the Paired‐SGS experiment, with dissipative cases producing somewhat denser cloud fields.

The simulations of the Mixed‐NSGS experiment feature slightly higher cloud fractions (Figure 2) than the simulations of the Mixed‐SGS experiment. However, in either experiment they are not as high as the cloud fractions simulated in the Paired‐SGS experiment. The 25MN, 45MN, and 65MN cases differ from the 25MS, 45MS, and 65MS of the Mixed‐SGS experiment only by elimination of the closure for SGS scalar fluxes. Yet, merely turning off the explict SGS scalar flux closure leads to a near doubling of LWP (Figure 3) and a substantial increase in cloud fraction.

The simulations in the Paired‐NSGS experiment produce the densest cloud fields among the four experimental cases, with LWP a factor of almost 6 larger than the highest value obtained in the Mixed‐SGS experiment (Figure 3). In all cases, the cloud cover averaged over the last 2 h of the simulations is 100 percent, consistent with observations (Figure 2). Not surprisingly, then, the profiles of q (Figure 1) show the highest fidelity to the observed values, particularly in the interior of the cloudy layer. They offer a substantial improvement in the lower parts of the cloud over the 55MS, 77MS, and 99MS cases of the Paired‐SGS experiment.

In summary, we find:

The simulations in the Paired‐NSGS experiment produce the densest cloud fields, followed by the Paired‐SGS, Mixed‐NSGS, and Mixed‐SGS cases.

Liquid water path is increased by using paired scalar and momentum numerics and is further increased by turning off explicit SGS closures when using dissipative numerics.

For the Paired‐SGS and Mixed‐NSGS experiments using centered schemes, the density of the cloud fields is increased by using higher‐accuracy numerics, but similar sensitivity is not seen using dissipative numerical schemes.

Conserved Variables

Mean vertical profiles of q and θ are shown in Figures 4 and 5, along with the observed values from Stevens et al. [2005]. Here the observed profiles of q and θ are suggestive of a largely well‐mixed boundary layer extending from the surface layer to the EIL around 850 m [Stevens et al., 2003a, 2005]. All cases in the Paired‐SGS experiment (with the exception of 22MS) and the Paired‐NSGS experiment maintain a well‐mixed boundary layer up to the inversion. Qualitatively it is difficult to differentiate between any of the cases in these two experiments. Both the Paired‐SGS and Paired‐NSGS experiments maintain high fidelity to the observations. On the other hand, the results in the Mixed‐SGS and Mixed‐NSGS experiments differ significantly from the observed profiles and are not well mixed throughout the boundary layer. In these cases, the departure from a well‐mixed state is consistent with decoupling of the cloudy layer from the surface moisture source, driven by excessive entrainment across the EIL. This is consistent with the reduction of LWP and cloud fraction seen in the Mixed‐SGS and Mixed‐NSGS experiments. It is further supported by the moistening of the subcloud layer relative to observations shown by the q profiles (Figure 4) and the stabilization of the cloud layer evident in the θ profiles (Figure 5).

Figure 4

Same as in Figure 1 but for total water specific humidity q.

Figure 5

Same as in Figure 1 but for liquid water potential temperature θ.

Just above the inversion, the profiles of q for the paired 22MS, 44MS, and 66MS cases show a marked drying that does not appear in other experiments (Figure 4). This anomalous drying is related to the interaction of the oscillatory scalar numerics with the steep gradients at the EIL, and has been observed by others [e.g., Stevens et al., 2005; Matheou and Chung, 2014]. A similar, albeit less marked, effect is also seen in the profiles of θ for these three cases (Figure 5). The propensity of oscillatory scalar numerics to produce warmer, drier cloud layers associated with decoupling and reduced cloud amount is noted in the Stevens et al. [2005] LES intercomparison. It is a primary motivation for using dissipative scalar numerics, which can reduce or eliminate this drying.

In summary, the conserved‐variable profiles show:

The simulations with the least dense cloud fields (all simulations in the Mixed‐SGS experiments and the 22MS case in the Paired‐SGS experiment) show clear evidence of decoupling from the surface moisture and enthalpy source.

The simulations with the most dense clouds field feature a well‐mixed boundary layer up to the EIL.

Oscillatory scalar advection schemes lead to anomalous drying and warming above the EIL.

Vertical Velocity Statistics

Vertical Velocity Variance

Profiles of resolved vertical velocity variance are shown in Figure 6 along with the observed values from Stevens et al. [2005]. Only the resolved component of vertical velocity variance and other higher‐order statistics is shown here, as SGS statistical quantities are unavailable with implicit LES. However, the SGS variance is typically small relative to the resolved variance. The observed profiles of have a single maximum near 600 m, roughly at the cloud base, and are consistent with a well‐mixed stratocumulus boundary layer with strong coupling between the cloud layer and the surface.

Figure 6

Same as in Figure 1 but for resolved vertical velocity variance . The observations shown here combine the in situ and radar‐derived observations plotted separately in Stevens et al. [2003b, 2005].

All simulations in the Mixed‐SGS and Mixed‐NSGS experiments significantly underpredict relative to the observations and also fail to replicate the shape of the observed profile. All simulated profiles of have two local maxima: the first near 200 m and well below the cloud base, and the second near the EIL. The presence of two maxima in the variance profile is consistent with a transition to a decoupled and more cumulus‐like boundary layer and also agrees with the large reduction in cloud fraction seen in these experiments.

All simulations in the Paired‐SGS experiment again underpredict relative to the observations, although to a lesser extent than seen in the Mixed‐SGS experiment (with the exception of the 22MS case). Again excepting the 22MS case, all simulations have a single maximum in the resolved vertical velocity variance profile near the EIL. The simulations using oscillatory schemes in this case, namely 22MS, 44MS, and 66MS, predict lower than those using dissipative schemes. The 22MS case appears to have transitioned to a cumulus‐like state. There is clear evidence that increasing the nominal numerical accuracy of the advection schemes improves the fidelity of the simulated profiles to that observed for both oscillatory and dissipative schemes in the Paired‐SGS experiment and, to some extent, in the Mixed‐NSGS experiment; a similar trend is not seen in the Mixed‐SGS experiment.

All simulations in the Paired‐SGS experiment using oscillatory schemes produce significant nonzero variance above the cloud layer, a result that is only seen in this experiment. This occurs because the oscillatory scalar numerics generate a nonphysical local minimum in q and maximum in θ just above the EIL, which reduces the stratification above the cloud layer and allows production of variance there.

The profiles in the Mixed‐NSGS experiment, where dissipative scalar numerics are used without an SGS closure, are quite similar to variance profiles associated with the 22MS, 44MS, and 66MS cases in the Mixed‐SGS experiment. Again, the presence of two maxima in the profile indicates a tendency toward decoupling of the cloud layer from the surface.

Results from the Paired‐NSGS experiment most closely follow the observed profile, correctly locating the maximum variance near the cloud base as well as producing variance magnitudes in good agreement with observations. Like the observations, the simulated profile of indicates a well‐mixed stratocumulus boundary layer, with active turbulent mixing coupling the surface and cloud layers.

All cases underpredict in the subcloud layer relative to the observations, which appears systematically in LES of this case [Stevens et al., 2005; Kirkpatrick et al., 2006; Matheou and Chung, 2014]. While it is tempting to ascribe this underprediction to a systematic bias in LES or lack of resolution, it is difficult to determine if the bias reflects an issue in the simulations or is a result of observational limitations.

In summary, the profiles show:

All simulations in the Paired‐NSGS experiments achieve good agreement with observed profiles of .

All simulations in the Mixed‐SGS, Paired‐SGS, and Mixed‐NSGS cases underpredict at all levels below the EIL.

Simulations with the sparsest cloud fields most significantly underpredict variance relative to the observations and tend to produce profiles with a local minimum in variance near cloud base, suggesting a transition toward a cumulus‐like state.

Vertical Velocity Skewness S

Simulated profiles of the resolved vertical velocity skewness, are shown in Figure 7 along with the observations from Stevens et al. [2003b].

Figure 7

Same as in Figure 6 but for the resolved vertical velocity skewness S.

The observations suggest the S profile has an S shape with a maximum in the interior of the subcloud layer, a secondary maximum near cloud top, and minimum near the cloud base. The minimum in S at the cloud base corresponds roughly to the location of the maximum in observed . Positive values of S near the surface are indicative of convective motion driven by the surface fluxes of sensible and latent heat; negative values in and just below the cloud layer are characteristic of negatively buoyant plumes driven by cloud‐top radiative cooling [e.g., Hogan et al., 2009].

In the Mixed‐SGS experiment, all simulations fail to produce the observed negative values of Sw in the cloud layer. Consistent with the profiles of (Figure 6) and the reduced cloud fraction (Figure 2), the positive values of S seen here suggest a transition to a more cumulus‐like boundary layer driven by surface fluxes rather than cloud‐top radiative cooling.

The 22MS simulation in the Paired‐SGS experiment has a similar structure to the simulations in the Mixed‐SGS experiment, underpredicting the minimum values of S near the cloud base. However, the underprediction is slightly less severe. The 44MS and 66MS simulations maintain greater fidelity to observations than the simulations of the Mixed‐SGS experiment or the 22MS simulation. There is clear evidence that the higher‐accuracy oscillatory numerics provide better predictions of S, with the 66MS case producing nearly zero skewness at cloud base, consistent with the observations. The simulations using dissipative numerics (55MS, 77MS, and 99MS cases) more closely follow the observations, particularly close to the surface, and all outperform the oscillatory schemes in predicting the cloud base minimum in skewness.

The S profiles from the Mixed‐NSGS experiment maintain greater fidelity to observations than those of the 22MS, 44MS, and 66MS cases in the Paired‐SGS experiment; however, they still fail to produce the negative skewness seen in the observations.

The S profiles from the Paired‐NSGS experiment are similar to those of the simulations using dissipative numerics in the Paired‐SGS experiment, but with better agreement to the observations in the interior of the subcloud layer. Overall S profiles from the Paired‐NSGS case maintain the greatest fidelity to the observations.

In summary, the S profiles show:

Simulations using dissipative numerics both in the Paired‐SGS experiment and in the Paired‐NSGS experiment reproduce the observed profiles reasonably well.

Simulations using oscillatory numerics generally fail to produce the negative values of S seen in the cloud layer, consistent with the reduced density of the cloud fields in these cases and the commensurate reduction in turbulence driven by cloud‐top radiative cooling.

Discussion

Qualitative Interpretation

The numerical results reported here show the sensitivity of LES of stratocumulus boundary layers to the discretizations of the equations of motion and to the inclusion or exclusion of closures representing SGS transport of scalars and momentum. The results span a range similar to that seen in the Stevens et al. [2005] intercomparison, with liquid water paths varying by more than a factor 6 across all simulations and cloud fraction varying by almost a factor 2.

Comparison of the simulations to observations makes it clear that LES results from the Paired‐NSGS experiment (where dissipative numerics for scalar and momentum transport are used without an explicit SGS closure) possess the highest fidelity. This suggests that the dissipation inherent in the numerical schemes used in this experiment serves as an adequate SGS model for both scalar and momentum transport. Similar findings have been reported previously for dry and neutral boundary layers [e.g., Brown and MacVean, 2000] and for stratocumulus layers [Kurowski et al., 2009a; Pedersen et al., 2016]. However, it is remarkable that the implicit SGS closure works so well for stratocumulus, given that it knows nothing about the saturated versus unsaturated state of air in a grid box, nor has it any explicit corrections for buoyancy effects.

The fact that dissipative numerics act as an adequate SGS model illuminates why the simulations in the Mixed‐SGS experiment (where dissipative numerics are used for scalar transport, oscillatory numerics are used for momentum transport, and an SGS closure is used for both momentum and scalars) perform so poorly. The oscillatory numerics used for the momentum field generate grid‐scale oscillations, causing gradients estimated from the grid‐scale momentum field to be erroneously large. The effect this has on the modeled SGS scalar fluxes is evident upon considering the form of typical SGS closures. For example, eddy diffusivity closures model the SGS scalar flux as where D is an eddy diffusivity that must be determined from resolved fields. Commonly, the eddy diffusivity is not computed directly but rather is obtained from the eddy viscosity νt by assuming a constant turbulent Prandtl number (typically ) and using the identity . In the Smagorinsky‐Lilly model used here, the eddy viscosity νt is modeled as where c is the Smagorinsky constant, is the geometric mean of the grid spacings in each direction, f is a buoyancy factor which tends to zero with increasing stratification [e.g., Pressel et al., 2015], and is the magnitude of the strain rate tensor S given by

This suggests the following pathway though which errors in grid‐scale gradients of the momentum fields act on the SGS scalar flux closure:

Oscillatory momentum numerics lead to grid‐scale oscillations and thus overestimation of grid‐scale gradients.

Large gradients in the momentum field lead to large values of .

Large values of lead to large values of ν, which then propagate to D through the Prandtl number relation.

Large values of D lead to anomalously large values of the SGS scalar flux and thus lead to excessive turbulent mixing.

While large values of eddy viscosity are needed in the Mixed‐SGS experiment to provide sufficient dissipation to prevent a pile up of kinetic energy at the grid scale, the resulting eddy diffusivity values are not appropriate for representing SGS fluxes of a scalar field that is already being acted upon by the dissipation inherent in the scalar transport scheme. Effectively two scalar SGS transport schemes are being used, one that is explicit in the SGS closure and a second that is implicit in the dissipative scalar transport scheme. The net effect of the two sources of SGS transport is to introduce too much mixing at the EIL, eventually leading to decoupling of the cloud layer from the surface.

The results from the Paired‐SGS and the Mixed‐NSGS experiments provide additional support for this mechanism. In the Paired‐SGS experiment, oscillatory transport schemes are used for momentum and scalar fields, so both fields are likely affected by grid‐scale numerical oscillations. In the Paired‐SGS experiment, as in the Mixed‐SGS experiment, the numerical oscillations in the momentum field lead to large values of D through its dependence on the magnitude of the strain rate tensor. But in the Paired‐SGS experiment, the larger values of D play an important role in dissipating numerical oscillations imparted by the oscillatory scalar transport schemes; the scalars are not subject to implicit numerical dissipation. Hence, there is effectively only one scalar SGS scheme, and the fidelity of the Paired‐SGS experiment is increased over that of the Mixed‐SGS experiment. Similarly, the Mixed‐NSGS experiment shows increased fidelity to observations over the Mixed‐SGS experiment. Here, the numerical dissipation inherent in the scalar discretization serves as an implicit SGS closure, and thus SGS fluxes are not explicitly closed in terms of grid‐scale gradients in the velocity fields. This elimination of explicit dependence of the SGS scalar fluxes on gradients in the velocity field removes the possibility that spuriously large gradients arising from oscillatory treatments of momentum will produce excessively large values of D that lead to excess entrainment and a degradation of the fidelity of the LES.

A further conclusion that can be drawn from the Paired‐SGS experiment (in particular the 55MS, 77MS, and 99MS cases) is that there is compensation between the effects of the SGS closure and dissipative numerics when dissipative numerics are used for both scalar and momentum transport. This is evident from the relatively high fidelity of these simulations, despite having both numerical dissipation and an explicit SGS closure. The likely mechanism here is that numerical dissipation provided by the transport scheme reduces grid‐scale gradients in the momentum field, yielding reduced eddy diffusivity values that are more appropriate for a scalar field that has already been affected by numerical dissipation. This effect was first noted in the atmospheric context by Brown and MacVean [2000] in their studies of neutral and convective boundary layers but clearly holds here in simulations of stratocumulus.

Here we have considered the interaction of numerical error with SGS closures only in terms of a Smagorinsky‐Lilly closure. Similar conclusions can be drawn from simulations using a TKE‐based closure (not shown). However, this is not surprising given that both closures are based on gradients of the resolved momentum field, either directly (the Smagorinsky‐Lilly closure) or indirectly (through the shear production term of the SGS TKE equation). We have also investigated the use of a Smagorinsky‐Lilly closure that has anisotropic eddy diffusivities (corresponding to the anisotropic discretization) but still guarantees a symmetric SGS stress tensor, similar to that used in the System for Atmospheric Modeling [Khairoutdinov and Randall, 2003] code as described in the Stevens et al. [2005] intercomparison. LES using this modified closure (not shown) only show modest differences from simulations using the standard isotropic closure.

While less widely used than Smagorinksy‐Lilly and TKE type closures, several other forms of SGS models exist [e.g., Kirkpatrick et al., 2003, 2006; Chung and Matheou, 2014; Matheou and Chung, 2014], and it is useful to consider our results in the context of these alternatives. Kirkpatrick et al. [2006] shows that high‐fidelity simulations of the DYCOMS‐II RF01 case are achievable using a dynamic SGS model along with dissipative scalar and momentum numerics (in this case, the Modified Utopia scheme of Stevens and Bretherton [1996]). In a dynamic modeling approach, SGS model coefficients are locally estimated from the small scales of the resolved fields based on scale similarity arguments, and the effects of numerical error are implicitly included in the estimated coefficients. As their simulations use both SGS closures and dissipative numerics, their configuration is most analogous to our 55MS, 77MS, and 99MS cases. However, the simulations using a dynamic SGS model show greater fidelity to observations than those three cases, and in fact are most similar to our Mixed‐NSGS and Paired‐NSGS experiments. This suggests that their dynamic procedure increases the ability of the SGS scheme to compensate for the numerical dissipation implicit in the transport schemes.

The simulations described by Matheou and Chung [2014] use a stretched‐vortex SGS model [Chung and Matheou, 2014] along with oscillatory numerics to simulate the DYCOMS‐II RF01 case. They perform a convergence study on isotropic grids with resolutions ranging from 10 to 2.5 m. The difference in resolution between their simulations and ours makes direct comparisons difficult; however, several points are worth noting. First, there is clear evidence that numerical oscillations produced by the oscillatory scalar advection schemes used in their simulations generate spurious mixing just above cloud top, which is manifested most clearly by a local minimum of q. Second, it is clear that with increasing resolution, their simulations are converging toward the observed values, and thus also toward our Mixed‐NSGS and Paired‐NSGS experiments results. However, their study shows that even at 2.5 m resolution, their simulations have not converged in any meaningful way, underscoring the difficulty of simulating stratocumulus clouds. That being said, we have not considered grid convergence as part of the present study, and it is likely that the conclusions regarding the relative fidelity of the various approaches are dependent on grid resolution as well as the cell aspect ratio [e.g., Pedersen et al., 2016].

A Simplified Quantitative Model

To gain deeper insight into the interaction of various forms of numerical error arising from the discretization of the advection terms with SGS closures, we consider a simplified model in which these interactions can be characterized quantitatively (see Appendix A for details). The model is based on simplified governing equations (A1), which preserve aspects of the anelastic equations used in the LES and consist of governing equations for a constant density fluid with velocity field u and a generic scalar . The scalar can be viewed as an analog of the specific entropy s or total water specific humidity q in the LES. However, unlike s and q, which are coupled to the velocity field through buoyancy, in the simplified model is only allowed to interact with the velocity field through the advection.

In Appendix A, equivalent equations [e.g., LeVeque, 2002] are used to represent the various advection discretizations and SGS closures described in section 2.3. From the master equivalent equations (A8) for the discretized momentum field and discretized scalar field at mesh spacing h, the kinetic energy dissipation relation can be derived (Appendix A). This dissipation relation states that in discrete solutions to the simplified model, dissipation of kinetic energy arises from two sources: dissipation of kinetic energy by the SGS closure (the first term on the RHS), and dissipation of kinetic energy by the numerics (the second term on the RHS). Here p is the order of accuracy of the momentum advection scheme and the exponent r > 0 models the spurious dispersive oscillations that can arise near sharp velocity gradients, due to centered finite difference discretizations of derivatives and due to the possible vanishing of the SGS terms in regions of stably stratified flow as expected at and above the inversion that marks the EIL. The determination of a numerical value of r is discussed in Appendix A. Similarly, a relation for scalar variance dissipation can be derived from the master equivalent equation for the scalar:

The scalar variance dissipation likewise arises from two sources: dissipation by the SGS closure (the first term on the RHS), and dissipation by the numerics (the second term on the RHS). Here, s is the order of accuracy of the scalar advection scheme. The constants depend on the combination of numerical schemes and SGS closures being used. For example, the constants are nonzero when an SGS closure is being used, and the constants are nonzero when dissipative numerics are being used. The constants and the RHS terms of the dissipation relations are described in greater detail in Appendix A. Note that the integrand of the second term on the RHS of equation (10) takes the form of the magnitude of the scalar gradient raised to the power of the order of the scalar advection scheme plus one, which suggests that this term may exhibit anomalous scaling [e.g., Shraiman and Siggia, 2000] with increasing accuracy of the scalar advection scheme. Given the numerical schemes used in this study and the focus of our analysis on the area surrounding the inversion, anomalous scaling is likely not of consequence. However, all of the arguments made here can be generalized to include a correction for anomalous scaling.

The dissipation relations (9) and (10) play a key role in what follows. The numerical properties of each experiment described in section 2.3 can be mimicked in the dissipation relations by selecting appropriate constants . From the resulting scalar variance dissipation relation for each of the experiments, we can estimate the magnitude of the dissipation and its convergence properties with mesh refinement. The dissipation relations provide an analytical basis for understanding the interaction of numerical error and the SGS closure, on which the fidelity of our various LES experiments depended. Our goal here is to use the simplified model to rank the various LES configurations according to how much dissipation and hence spurious turbulent mixing they are likely to produce at the EIL.

In the Paired‐NSGS experiment, because no explicit SGS closure is included and because both the scalar and momentum advection schemes are dissipative. This means that the kinetic energy and scalar variance dissipation relations can be written as and

Following the Onsager interpretation of the Kolmogorov (K41) hypothesis [Kolmogorov, 1991] for homogeneous, isotropic turbulence [Robert, 2003], we assume that the velocity field is Hölder continuous with an exponent α such that for each velocity component u of u and for any from which we deduce that the underlying velocity gradient must satisfy

The K41 hypothesis predicts at scales in the inertial subrange of homogeneous isotropic turbulence [Kolmogorov, 1991]; however, the arguments made here hold for a range of α, which is important near the EIL, where turbulence is neither homogeneous nor isotropic. Therefore, since h is not a function of x, the dissipation relations (11) and (12) can be simplified to and

Thus, as long as (which is the case for the higher‐order WENO schemes in PyCLES) [see Parés Pulido, 2016], the exponents θ1 and θ2 are greater than zero, so the dissipation of both kinetic energy and scalar variance decreases with decreasing mesh spacing.

In the Mixed‐SGS experiment, oscillatory numerics and an SGS closure are used for the momentum equation, so and . For the scalar equations, dissipative numerics and an SGS closure are used, so . First, we consider the kinetic energy dissipation by substituting C 1 and C 2 into the dissipation relation (9) and integrating in time. Making use of the finiteness of kinetic energy, we obtain where C is a constant. If we assume the discrete velocity gradients are of the same scale across the domain D and throughout the time interval , i.e., this expression can be substituted into equation (17), giving

Because we have assumed that r > 0, comparison of equation (19) with the discrete velocity gradients predicted by the Onsager interpretation of the Kolmogorov (K41) hypothesis given in equation (14) makes it clear that using oscillatory momentum numerics produces discrete velocity gradients that are not consistent with the Onsager interpretation. The discrete velocity gradients produced by oscillatory schemes coupled with Smagorinsky‐Lilly type closures for the eddy viscosity are likely to be large, and they are increasing with mesh refinement relative to those predicted by the Onsager interpretation. This strongly suggests the potential for the oscillatory schemes to produce dispersive oscillations, particularly in the part of the domain where the eddy viscosity vanishes.

For the scalar variance dissipation, substituting equation (19) into the dissipation relation (10) and using equation (17) and the proportionality of scalar and velocity gradients, as deduced in Appendix A, then

In this expression, the dissipation from the SGS scheme is a constant C 3 and is independent of the mesh spacing h. However, the numerical dissipation due to the WENO scheme scales with , where

Given that , the two contributions to the dissipation of are in balance and independent of mesh spacing as long as

This suggests that if a formally high‐order WENO scheme with is used for the scalar transport, the dominant contribution to the scalar variance dissipation stems from the eddy diffusivity. On the other hand, a low‐order WENO scheme with s < 2 would lead to even greater numerical diffusion of the scalar.

In the Mixed‐NSGS experiment, the momentum equation is discretized with oscillatory schemes and the eddy viscosity term is retained. Thus the equivalent equation for momentum is identical to that for the Mixed‐SGS scheme, meaning that discrete velocity gradients behave as in equation (19); the velocity field exhibits spurious dispersive oscillations induced by the oscillatory numerics. The dissipation relation for the scalar variance is deduced by substituting equations (19) and (20) into equation (12) and using , leading to

As before, we assume that if equation (23) holds, the numerical dissipation of the scalar is constant and independent of mesh spacing. However, if either s < 2 (e.g., lower‐order numerical accuracy) or (large dispersive oscillations because of central differencing), the dissipation is even greater.

The Paired‐SGS experiment can be divided into two groups: oscillatory numerics, denoted Paired‐SGS(O), and dissipative numerics, denoted Paired‐SGS(D).

Paired‐SGS(O)

In this set of experiments, and because central difference schemes are used to represent the advective terms, and the Smagorinsky‐Lilly closure is used to represent subgrid‐scale fluxes. As the momentum equation is discretized in the same manner as in the Mixed‐SGS and Mixed‐NSGS schemes, the discrete velocity gradients behave as in equation (19). Substituting equation (19) into (12) and using gives that is, the scalar variance dissipation is constant and independent of mesh spacing.

Paired‐SGS(D)

In this set of experiments, because WENO schemes are used to approximate the momentum and scalar transport equations, and eddy viscosity and diffusivity terms are retained in both. As before, we assume the discrete velocity gradients follow equation (14). Substituting equation (14) into (9) gives the energy dissipation relation where and . Here, for and r > 0, and as long as . For , the dissipation is dominated by the eddy viscosity term, which is proportional to C 1. Hence, the momentum discretization in the Paired‐SGS(D) experiment using dissipative numerics will be considerably more dissipative than in the Paired‐NSGS experiment.

The scalar variance dissipation relation, which can be derived by substituting equations (14) and (20) into equation (10) and using , is

Because , r > 0, and , both terms on the right hand side of the scalar variance dissipation relation decrease with mesh refinement; it is possible to choose a mesh resolution for which the numerical dissipation of the scalar transport is small.

Interpretation of LES Results With Simple Model

We have argued that LES of stratocumulus clouds are degraded by spurious numerical entrainment at the cloud top, which can be driven by diffusive scalar fluxes arising from the interaction of SGS closures with numerical errors. The simple model suggests an ordering of the fidelity of the various LES experiments. We use the notation if experiment A diffuses scalars more than experiment B, meaning that B is superior to A. We propose the following ordered relation:

This ordering can be justified on the basis of the dissipation relations:

In the Mixed‐SGS experiment, the numerical dissipation (21) has contributions from both the eddy diffusivity and numerical dissipation terms. As long as and , both contributions are approximately of the same order and do not decrease with mesh refinement, leading to poor numerical performance.

In the Mixed‐NSGS experiment, the scalar variance dissipation (24) also does not decrease with mesh refinement as long as . However, the Mixed‐SGS experiment is even more diffusive because . However, if r is large (but r < 1), implying that the central difference approximation of the momentum equation is more oscillatory, the WENO dissipation term proportional to C 4 is much greater than the eddy diffusivity term. In this case, one would expect only minor differences between the Mixed‐SGS and Mixed‐NSGS experiments.

In the Paired‐SGS experiment with oscillatory numerics, the scalar variance dissipation (25) differs from its counterpart (24) for the Mixed‐SGS experiment through the relative size of the constants C 3 and C 4, provided and . However, if r > 0, the numerical dissipation arising from the dissipative numerics will dominate the mesh‐independent dissipation due to the eddy diffusivity. In this case, the Paired‐SGS(O) experiment will outperform the mixed‐NSGS experiment.

In the Paired‐SGS experiment with dissipative numerics, the scalar variance dissipation (27) decreases under mesh refinement as long as , and . Hence, it is generally smaller than in the Paired‐SGS(O) experiment. If r is not very small, the eddy diffusivity term dominates the numerical dissipation.

In the Paired‐NSGS experiment, the scalar variance dissipation (12) is smallest because, unlike in the scalar variance dissipation (27) for the Paired‐SGS(D) experiment, there is no additional contribution from an eddy diffusivity term. The latter dominates in Paired‐SGS(D) if r is not small.

In summary, the scalar variance dissipation relations from the simplified model provide a basis for ranking (28) the experiments according to how much scalar mixing they produce. This ranking is consistent with the numerical results in section 3 and provides quantitative validation of the qualitative explanation given in section 4.1.

Summary and Conclusion

We have shown that LES of stratocumulus are highly sensitive to numerical error in the discretization of transport terms and to the interaction of that numerical error with SGS closures. By comparing numerical results to observations, we have shown that the configuration that performs best for this stratocumulus case (Paired‐NSGS experiment) uses dissipative scalar and momentum numerics and no explicit SGS closures. Within the Paired‐NSGS experiment, the simulations show little sensitivity to the accuracy of the numerical schemes. However, as all schemes used are based on high‐order WENO reconstructions, the numerical accuracy is already high even for the lowest‐accuracy scheme considered here (WENO5).

Using oscillatory momentum transport, schemes with dissipative scalar transport schemes and standard turbulence closures (Mixed‐SGS experiment) leads to significant reductions in the fidelity of LES. We argue that the grid‐scale numerical oscillations generated by the oscillatory momentum numerics lead to excessively large values of the eddy diffusivity, which augments the dissipation already inherent in the scalar transport scheme. This redundant additional SGS scalar diffusion drives excessive cloud top entrainment, leading to thinning of the cloud, reduced buoyancy production, and eventual decoupling of the cloud layer from the surface.

A simplified analytical model quantitatively characterizes the interaction of numerical error with the SGS closure. The analysis of the simplified model allowed us to order the various numerical experiments according to how much scalar variance dissipation they produce. The ordering (28) obtained from the simplified model and the ordering obtained empirically from the LES experiments are identical, supporting our argument that the interaction of numerical error with the SGS closure controls the fidelity of LES.

Despite the modular design of modern atmospheric LES codes, which makes it easy to mix and match various aspects of their numerical formulations and SGS closures, the results here make it clear that special care must be taken when doing so. SGS closures should not be selected independently of the choices made for the scalar and momentum transport numerics. The large sensitivities seen here also underscore the importance of detailed reporting of the numerical schemes and SGS closures used in LES studies to ensure the reproducibility of published results.

The conclusions from the numerical experiments and the simplified model strongly support the use of implicit LES (as in our Paired‐NSGS experiment) for the simulation of stratocumulus clouds, corroborating and extending the results of Kurowski et al. [2009a] and Pedersen et al. [2016]. Implicit LES avoids spuriously large SGS mixing resulting from the interaction of numerical error in the vicinity of the inversion with the SGS closure. This is essential for the ability of implicit LES with WENO schemes to yield high‐fidelity simulations of stratocumulus clouds, with a sharp transition in mixing strength across the inversion. In fact, implicit LES with WENO schemes achieves high fidelity already at the relatively modest resolution we considered, at a fraction of the computational expense of approaches that require much higher resolution to achieve comparable results.

Generalizing our results to simulations of other cloud types and to LES using different numerical schemes and SGS closures is difficult because the optimal numerical and SGS modeling approach may be flow and grid‐dependent [e.g., Pedersen et al., 2016]. Moreover, we have considered only a few of the numerical schemes and SGS closures used in LES. Nonetheless, a key result of this work, that care must be taken to consider the interaction of numerical error with the SGS closures used in atmospheric LES, especially where sharp transitions in mixing strength occur, still holds true and should serve as a motivation for further studies of this type.