Galilean invariance of shallow cumulus convection large-eddy simulations
aa r X i v : . [ phy s i c s . a o - ph ] A ug Galilean invariance of shallow cumulus convection large-eddysimulations
Oumaima Lamaakel and Georgios MatheouDepartment of Mechanical Engineering, University of Connecticut, Storrs, ConnecticutAugust 20, 2020
Abstract
In large-eddy simulations (LES) a computational-domain translation velocity can be usedto improve performance by allowing longer time-step intervals. The continuous equations areGalilean invariant, however, standard finite-difference-based discretizations are not discretelyinvariant with the error being proportional to the product of the local translation velocity andthe truncation error. Even though such numerical errors are expected to be small, it is shownthat in LES of buoyant convection the turbulent large-scale flow organization can modulateand amplify the error. Galilean invariance of global flow statistics is observed in well-resolveddirect numerical simulations (DNS). In LES of single-phase convection under an inversion, flowstatistics are nearly Galilean invariant and do not depend on the order of accuracy of the finitedifference approximation. In contrast, in LES of cloudy convection, flow statistics show strongdependence on the frame of reference and the order of approximation. The error with respectto the frame of reference becomes negligible as the order of accuracy is increased from second tosixth in the present LES. Schemes with low resolving power can produce large dispersion errorsin the surface-fixed frame that can be amplified by large-scale flow anisotropies, such as strongupdrafts rising in a non-turbulent free troposphere in cumulus-cloud layers. Interestingly, in thepresent large-eddy simulations, a second-order discretization in the proper Galilean frame canyield comparable accuracy as a high-order scheme in the surface-fixed frame.
The effects of shallow clouds are one of the largest sources of uncertainty in climate projections(Bony and Dufresne, 2005; Rieck et al., 2012; Klein et al., 2017; Vial et al., 2017; Zelinka et al.,2017). Shallow clouds form in the atmospheric boundary layer, the lower part of the atmospherethat is in contact with the surface, and can reach heights up to 4 km. Large-eddy simulation (LES)is currently the best-available cloud modeling technique, since the range of flow scales is too largefor direct numerical simulation (DNS) methods.LES resolves all dynamically important flow scales and models the smaller, more “generic”in nature (Pope, 2004). In LES of atmospheric boundary layers grid spacings typically range in5–50 m enabling simulations to explicitly resolve individual cloud shapes and the structure of up-drafts, downdrafts, and the entrainment process. LES applications typically include atmosphericphysics investigations, development and evaluation of weather and climate model parameteriza-tions using LES as a reference model, and, more recently, a promising weather forecast model (e.g.,1chalkwijk et al., 2015). The first LES were simulations of atmospheric boundary layers (Deardorff,1972; Sommeria, 1976) and the ubiquitous Smagorinsky subgrid scale (SGS) model was first for-mulated in the context of a groundbreaking atmospheric general circulation model (Smagorinsky,1963).LES simulations of boundary layer clouds are challenging because of the multiscale spatialorganization of convection, which requires large computational domains, and the long time inte-grations needed to capture the evolution of convection and diurnal cycle effects. Often, integra-tions with explicit time marching schemes are used, e.g., Runge–Kutta methods (Heus et al., 2010;Fuka and Brechler, 2011; Maronga et al., 2015; van Heerwaarden et al., 2017; Lac et al., 2018) orAdams–Bashforsh (Khairoutdinov and Randall, 2003; Basu and Port´e-Agel, 2006; Huang and Bou-Zeid,2013), which adhere to an advection-dominated time-stability constraint. The time step interval,∆ t , depends on the Courant–Friedrichs–Lewy (CFL) number,CFL = ∆ t (cid:18) | u | ∆ x + | v | ∆ y + | w | ∆ z (cid:19) , (1)where [ u, v, w ] is the velocity vector, and ∆ x , ∆ y , and ∆ z are the spatial grid spacings. Theinterval of the next time step ∆ t n +1 can be estimated based on the previous ∆ t n and CFL n , anda scheme-dependent maximum CFL number, CFL max ,∆ t n +1 = min ijk CFL max
CLF n ∆ t n , (2)where the minimum is taken over all grid points.For a given computational domain size and simulation time length, the computation expensedecreases as the grid spacing increases, i.e., ∆ t is proportional to the spatial grid size ∆ x in(2). Coarser grids are preferred for computational efficiency. However, because a larger fractionof spatial flow scales remains unresolved when coarse grids are used, LES models require skillfulturbulence parameterizations to maintain the fidelity of the simulation as the grid becomes coarser(Matheou and Chung, 2014; Matheou and Teixeira, 2019).Another approach to gain computational advantage exploits the dependence of ∆ t on the ab-solute value of the velocity field components in (2): ∆ t increases for smaller absolute velocitycomponents. Also, the equations of motion are Galilean invariant and do not depend on the abso-lute value of the velocity. That is, the equations of motion are invariant under the transformation τ = t (3) y = x − u t (4) v ( y , τ ) = u ( x − u t, t ) + u , (5)where u is a constant velocity vector. In general, the boundary conditions are not Galileaninvariant. However, there are particular cases, such as in LES with doubly periodic boundaryconditions in the horizontal directions and spatially homogeneous surface conditions where a movingreference frame can be used.LES of atmospheric boundary layers in a moving reference frame is equivalent to the compu-tational domain translating with a constant horizontal velocity with components [ u , v ]. Onlyhorizontal translation vectors u can be employed because the vertical velocity must be zero at the2urface (no penetration condition). When a Galilean transformation is applied, the surface bound-ary condition is modified using (4). The Galilean frame is chosen such that the domain average ofthe horizontal wind speed ( u + v ) / is minimized across all grid points, i.e., the computationaldomain translating with the domain-mean horizontal wind. Different domain translation velocities u result in changes in local velocity vector u ( t, x, y, z ) and, following the term of Wyant et al.(2018), differences in horizontal-direction cross-grid flow. As shown in the present results, the com-putational cost savings can be significant when a suitable frame of reference is chosen. Differencesin computer program execution speed by about a factor of two between LES in the surface-fixedand Galilean frames can be achieved.In contrast to the continuous equations of motion, not all discrete numerical approximationsare Galilean invariant. The breakdown of Galilean invariance can originate from non-linear discretedifference operators where the error is introduced as artificial numerical dissipation proportionalto a non-dimensional absolute velocity, such as a CFL number (e.g., Lomax et al., 2003, p. 195).Presently, more subtle issues related to the Galilean invariance of atmospheric boundary layer LESare investigated. Even though finite difference methods can preserve the Galilean invariance of thespatial derivatives, when finite differences are applied in discretizations of the advection term , thesemi-discrete form of the momentum includes an additional term, the third term in the left handside of d v i d τ + N ( v ) i − u ,i (cid:20)(cid:18) ∂v∂y (cid:19) i − D ( v ) i (cid:21) = RHS( v ) i , (6)where N is the discrete advection term operator, D is the discrete first derivative approximations,and the right hand side (RHS) includes all remaining diffusive terms, resolved and subgrid scale(Bernardini et al., 2013, eq. 4). As discussed in Bernardini et al. (2013), the nature of the error isdispersive and the error is locally proportional to the product of the translation velocity components, u ,i , and the truncation error of the first derivative. Special symmetry preserving schemes havebeen developed to discretely satisfy Galilean invariance (e.g., Bihlo and Nave, 2014).Previous investigations of Galilean effects using finite difference schemes focused on simplified-flow models, such as the one-dimensional Burgers equation, and on small-scale flow features in DNS.For instance, Bernardini et al. (2013) discuss errors on the high wavenumber content of the flowfield. In atmospheric boundary layer LES, Matheou et al. (2011) and Wyant et al. (2018) discussGalilean invariance effects on the mean flow and domain-averaged statistics (e.g., turbulent fluxesand cloud properties). The mechanisms leading to the relatively large differences in the simulationsof Matheou et al. (2011) and Wyant et al. (2018) are not clear because of the use of non-invariantdissipative numerical discretizations. Even though the implications of using monotone schemes onflow statistics can be significant, and a violation of the Galilean invariance property of the equationsof motion, the use of such schemes is often a necessary tradeoff because of the preservation of thephysical bounds of scalar variables, e.g., temperature and humidity. In astrophysical simulations,Springel (2010) shows how the impact of numerical dissipation can be modulated using movinggrids and, at the same time, illustrates notable changes in the numerical solution with respect tothe frame of reference.LES of shallow convection includes the added challenges of active scalar turbulent transport ina multi-phase flow. For instance, Grabowski and Smolarkiewicz (1990) discuss numerical artifactsthat can develop in an advection–condensation problem. Usually, in LES, the mean state of the the term “advection” is used to refer to the non-linear term, u · ∇ u , as customary in the atmospheric modelingliterature. The term “convection” refers to buoyant convection in the fluid. a )the interaction of dispersion errors with the large-scale flow anisotropy, and ( b ) multi-phase floweffects that can potentially weaken some of the flow field smoothness assumptions of the numericalapproximation. Focusing on shallow convection, two types of simulations are carried out: ( a )DNS, where the flow is fully resolved and no SGS model is used, i.e., the physical viscosity of thefluid provides all dissipation; and ( b ) LES, where the flow is not fully resolved and relatively largegradients are present at the grid scale. No explicit filtering is performed in either LES and DNS.Three convection cases are explored: shallow cumulus convection based on the conditions observedduring the Cumulus over the Ocean (RICO) campaign (Rauber et al., 2007; vanZanten et al., 2011),buoyant bubble simulations (a simple model for cumulus-topped updrafts) and dry (cloud free)convection.The focus of the analysis is on differences with respect to the frame of reference of domain-averaged boundary layer statistics (e.g., cloud cover vs. time, turbulent fluxes vs. height) becausethese statistics are used to tune and evaluate weather and climate model convection parameteriza-tions (e.g., Siebesma and Holtslag, 1996; Siebesma et al., 2007; Teixeira et al., 2008; Neggers et al.,2009; Witek et al., 2011).The outline of the study is as follows. The model formulation, the simulation cases, and themain flow statistics are presented in Section 2. The results are discussed in Section 3. Section 4provides support for the large-scale flow anisotropy hypothesis and discusses the mechanism of erroramplification. Finally, conclusions are summarized in Section 5. A unified numerical model is used to perform both DNS and LES. When LES is carried out, aturbulence SGS model is used to account for the effects of the unresolved motions on the resolved-scale variables. The contribution of the resolved-scale viscous terms is neglected in LES, i.e., aninfinite Reynolds number flow is considered. When DNS is carried out, the SGS model terms arenot computed and all dissipation is provided by the viscous terms.The LES model of Matheou and Chung (2014) is used with the addition of viscous terms. Theconservation equations for mass, momentum, liquid water potential temperature θ l , and total water4 t mixing ratio on an f -plane, are, respectively, ∂ ¯ ρ ˜ u i ∂x i = 0 , (7) ∂ ¯ ρ ˜ u i ∂t + ∂ (¯ ρ ˜ u i ˜ u j ) ∂x j = − θ ¯ ρ ∂ ¯ π ∂x i + δ i g ¯ ρ ˜ θ v − h ˜ θ v i x θ − ǫ ijk ¯ ρ f j (˜ u k − u g,k ) − ∂τ ij ∂x j + ∂d ij ∂x j , (8) ∂ ¯ ρ ˜ θ l ∂t + ∂ ¯ ρ ˜ θ l ˜ u j ∂x j = − ∂σ θ,j ∂x j + ∂∂x j (cid:18) ρ D θ ∂θ l ∂x j (cid:19) + ˜ S θ , (9) ∂ ¯ ρ ˜ q t ∂t + ∂ ¯ ρ ˜ q t ˜ u j ∂x j = − ∂σ q,j ∂x j + ∂∂x j (cid:18) ρ D q ∂q t ∂x j (cid:19) + ˜ S q . (10)The Cartesian coordinates are ( { zonal , meridional , vertical } = { x , x , x } = { x, y, z } ) and thecomponents of the velocity vector and geostrophic wind, are u i and u g,i , respectively. The Coriolisparameter is f = [0 , , f ], θ is the constant basic-state potential temperature, ρ ( z ) is the density, π is the dynamic part of the Exner function that satisfies the anelastic constraint (7), and θ v isthe virtual potential temperature. The angled brackets h•i x denote an instantaneous horizontalaverage.When LES is performed, the prognostic variables u i , θ l , and q t are defined as Favre-filteredvariables ˜ φ ≡ ρφ/ ¯ ρ , where ρ is the density and the overbar denotes a spatially filtered variable.When the flow is fully resolved (i.e., in DNS), u i , θ l , and q t correspond to the local values (withoutany filtering or averaging), thus tildes and overbars are not needed.The viscous stress tensor is d ij = 2 µ D ij (11)where µ is the dynamic viscosity coefficient, which is assumed constant presently, and D ij is thedeviatoric rate of strain tensor, D ij = 12 (cid:18) ∂u i ∂x j + ∂u j ∂x i (cid:19) − δ ij ∂u k ∂x k . (12)The Fickian diffusion coefficients D θ and D q are related to the momentum coefficient through thePrandtl and Schmidt numbers, P r ≡ ν D θ = 0 . , (13) Sc ≡ ν D q = 1 , (14)where ν = µ/ρ is the kinematic viscosity.The subgrid-scale (SGS) stress tensor and scalar flux are modeled using an eddy-diffusivityassumption τ ij = − ρ ν t ˜ D ij , (15)and σ j,φ = − ¯ ρ ν t Pr t ∂ ˜ φ∂x j . (16)5he eddy diffusivity for all scalar variables is related to the SGS momentum diffusivity, ν t , throughthe constant model turbulent Prandtl and Schmidt numbers, Pr t = 0 .
33, Sc t = 0 . ν t = ∆ | ˜ D | f m (Ri) , (17)where ∆ = C s ∆ x is the characteristic SGS length scale, | ˜ D | = (2 ˜ D ij ˜ D ij ) / is the resolved-scale deformation, and f m a stability correction function (Lilly, 1962; Matheou, 2016). The value C s = 0 . θ l and q t through the source terms S θ and S q , which havea case-dependent form.Condensation is modeled based on the mean thermodynamic state in each grid cell. For all butone pair of simulations an “all or nothing” scheme is used, i.e., no partially saturated air in eachgrid cell is assumed. Two runs use a modified saturation scheme that allows for the formation ofliquid when the mean state is not saturated.The liquid water mixing ratio q l in the “all or nothing” scheme is q l = max(0 , q t − q s ) , (18)where q s ( p, T ) is the saturation mixing ratio. An ad hoc subgrid condensation scheme is constructedby using a smoother transition between the unsaturated and saturated regimes, q l = q t − q s < − . − q t − q s ) + 0 . q t − q s ) + 0 . − . ≤ q t − q s ≤ . − q t − q s if q t − q s > . − . (19)In the modified saturation scheme (19) a second degree polynomial is used to transition betweenthe two branches of (18).Liquid water is assumed suspended (i.e., no drizzle or precipitation is present) in all simula-tions, even though for the shallow cumulus case precipitation should develop as the boundary layerdeepens (vanZanten et al., 2011).At the surface, turbulent fluxes are estimated using either Monin–Obukhov similarity theory(MOST) or bulk aerodynamic formulae (vanZanten et al., 2011, eq. 1–4). In both approaches, thewind vector at the first model half level is needed. In the surface-fixed (non-moving) frame the local u ( t, x, y, ∆ z/
2) is used. In the Galilean frame, the translation u is added to form v ( t, x, y, ∆ z/
2) = u ( t, x, y, ∆ z/
2) + u , which is used to estimate the turbulent fluxes. All simulations are preformedin a doubly-periodic domain in the horizontal directions. A Rayleigh damping layer is used at thetop of the domain to limit gravity wave reflection.Spatial derivatives are approximated with centered finite differences. The family of fully con-servative schemes of Morinishi et al. (1998) adapted for the anelastic approximation is used to ap-proximate the momentum and scalar advection terms. The globally conserved quantities are P ρu i P ρφ , where φ is a passive scalar, and the sum taken over all grid points, see Morinishi et al.(1998) for further details. The second-, fourth-, and sixth-order approximations are used. Theproperties of the advection schemes are discussed in Matheou (2016) and Matheou and Dimotakis(2016). The derivatives in the SGS model and viscous terms are estimated using second-ordercentered differences.A third-order Runge–Kutta method is used for time integration (Spalart et al., 1991). The LESmodel was successfully used in several previous studies spanning a diverse set of meteorological con-ditions (Matheou et al., 2011; Inoue et al., 2014; Matheou and Chung, 2014; Matheou and Bowman,2016; Matheou, 2016; Thorpe et al., 2016; Matheou, 2018; Matheou and Teixeira, 2019; Jongaramrungruang et al.,2019; Couvreux et al., 2020). To create similar wind profiles and a uniform baseline for comparison, all simulations use thegeostrophic wind forcing the Cumulus over the Ocean (RICO) case (vanZanten et al., 2011). Thecomponents of the geostrophic wind are u g ( z ) = − . . × − z m s − and v g = − . − .The latitude is 18 ◦ N. Different initial temperature and humidity profiles and surface fluxes areused to initiate various types of convection. All cases are run in pairs: with and without a Galileantranslation velocity u . The translation velocity for the dry convection and shallow cumulus casesis ( − , −
4) m s − . Because the buoyant bubble case does not include surface shear, the mean windis somewhat different and the translation velocity is ( − , − .
8) m s − . The translation velocity isset based on the domain-averaged mean wind.For all cases the grid spacing is typical of similar studies in the literature (e.g., Margolin et al.,1999; Sullivan and Patton, 2011; vanZanten et al., 2011; Seifert and Heus, 2013). Grid spacing isuniform and isotropic ∆ x = ∆ y = ∆ z . The time step is adjusted to maintain CFL ≈ .
2. Table 1summarizes the LES runs, including the number of grid points used, grid spacing, and u . Resultsfrom additional sensitivity runs are discussed in the appendices. Table 2 compares the execution times and total number of time steps between the simulations inthe fixed and Galilean frames. The execution time in Table 2 is the wall-clock time length fromthe beginning to the end of the computer program, and includes all computation, I/O, calculationof flow statistics, and other synchronization and setup tasks. The total number of steps is aperformance metric that only depends on u , ∆ t , and ∆ x . All Galilean frame LES execute abouttwice as fast and require about half the number of time steps to complete compared to the fixedframe runs. The Galilean frame DNS completed 2.7 times faster than the fixed frame run. Ingeneral, the computational savings when using the Galilean frame are significant in all simulations. The temperature and humidity initial profiles of the RICO case are used. An initial sphericalpositively buoyant region with radius r = 200 m and center at z = r is created by increasing thevalues of θ l and q t by 10 % with respect to the standard (horizontally uniform) initial condition.7able 1: Summary of the cases simulated. The first and second columns correspond to the shortenedform of the simulation case and the convection type, respectively. Fully resolved simulations,without any SGS model, are denoted as DNS, whereas simulations of the full dynamics using aSGS model are labeled as LES. The grid spacing is denoted by ∆ x , N x = N y and N z are numberof horizontal and vertical grid points, respectively, u the the Galilean translation velocity, and“Advection” corresponds to the order of the advection scheme. For all runs ∆ x = ∆ y = ∆ z .The star ( ∗ ) denotes that a 10-member ensemble were carried out. The dagger ( † ) denotes thatadditional simulations with variable CFL numbers were carried out.Run Description Model ∆ x (m) N x N z u (m s − ) AdvectionBD2f † Buoyant bubble DNS 5 512 600 (0 ,
0) secondBD2g Buoyant bubble DNS 5 512 600 ( − , − .
8) secondBD4f Buoyant bubble DNS 5 512 600 (0 ,
0) fourthBD4g Buoyant bubble DNS 5 512 600 ( − , − .
8) fourthBHf Buoyant bubble LES 10 256 300 (0 ,
0) fourthBHg Buoyant bubble LES 10 256 300 ( − , − .
8) fourthBLf Buoyant bubble LES 20 128 150 (0 ,
0) fourthBLg Buoyant bubble LES 20 128 150 ( − , − .
8) fourthD2f Dry convection LES 40 512 100 (0 ,
0) secondD2g Dry convection LES 40 512 100 ( − , −
4) secondD4f ∗ Dry convection LES 40 512 100 (0 ,
0) fourthD4g Dry convection LES 40 512 100 ( − , −
4) fourthD6f Dry convection LES 40 512 100 (0 ,
0) sixthD6g Dry convection LES 40 512 100 ( − , −
4) sixthC2f Shallow Cu LES 40 1024 100 (0 ,
0) secondC2g Shallow Cu LES 40 1024 100 ( − , −
4) secondC4f ∗ †
Shallow Cu LES 40 1024 100 (0 ,
0) fourthC4g Shallow Cu LES 40 1024 100 ( − , −
4) fourthC6f Shallow Cu LES 40 1024 100 (0 ,
0) sixthC6g Shallow Cu LES 40 1024 100 ( − , −
4) sixthCSf Shallow Cu (mod. saturation) LES 40 1024 100 (0 ,
0) fourthCSg Shallow Cu (mod. saturation) LES 40 1024 100 ( − , −
4) fourthThe initial condition is given by φ ( x, y, z ) = [0 .
05 erf (0 . r − r )) + 1 . φ i ( z ) (20)where r = (cid:0) x + y + ( z − r ) (cid:1) / is the distance from the center of the sphere in meters, φ denoteseither θ l or q t , and φ i ( z ) the corresponding initial profile of the RICO case. The large-scale forcingof the RICO case is not included in the buoyant bubble simulations, i.e., S θ = 0 and S q = 0.Sensible and latent heat surface fluxes are set to zero. LES and DNS types of simulations arecarried out.In the DNS run, viscosity is set to µ = 2 kg m − s − . The Reynolds number, defined as8able 2: Comparison of execution time and number of time steps n of simulations in the fixedand Galilean frames. The number of CPU cores used is N r , t fixed and t Galilean are the wall clocktimes for runs in the fixed and Galilean frames, respectively, and Speedup is defined as the ratio t fixed /t Galilean .Run N r t fixed (h) t Galilean (h) Time Speedup n fixed n Galilean n fixed /n Galilean
BD4f/g 64 31.5 11.5 2.7 14863 5408 2.7BHf/g 64 0.79 0.35 2.3 2671 1137 2.3BLf/g 64 0.06 0.03 2 1327 554 2.4C2f/g 64 37.9 20.2 1.9 33291 17662 1.9C4f/g 64 44.8 21.5 2.1 31418 17678 1.7C6f/g 64 60.7 34.5 1.8 30537 17128 1.8CSf/g 144 5.3 2.8 1.9 29842 15793 1.9D4f/g 64 1.6 0.86 1.9 7191 3647 2 Re ≡ r w max ρ/µ , where w max is the maximum vertical velocity in the updraft, is Re ≈
200 at t = 0 .
25 h, when the first instance of cloud forms. The SGS model terms are set to zero in the DNSruns. The DNS resolution is chosen (by trial of different ∆ x since the flow is not fully turbulentand Kolmogorov scaling does not apply) such that the fourth-order scheme fully resolves the flowand the second-order marginally resolves the flow. The grid spacing is ∆ x = 5 m.In the LES runs the viscosity is set to zero and grid resolutions ∆ x = 10 m and ∆ x = 20 m areused.Because surface shear cannot be resolved, a slip (no penetration, no stress) surface condition isused. The DNS simulations are run for 1 h and the LES simulations for 0 . w ( t, x, y, z = 0) = 0 and ∂ { u, v, q t , θ l } /∂z = 0 at z = 0 is applied),they can be used to verify that the breakdown of Galilean invariance is not a boundary conditionartifact. The cloud-free (i.e., “dry”) convection case of Matheou et al. (2011) is modified by the addition ofgeostrophic wind forcing. The resulting case is somewhat unphysical because the wind profile doesnot correspond to the boundary layer thermodynamic profiles. The initial potential temperaturelapse rate is 2 K km − , with θ ( z = 0) = 297 K. The initial total water mixing ratio lapse rate is − .
37 g kg − km − up to z = 1350 m and − .
94 g kg − km − higher up with q t ( z = 0) = 5 g kg − .The temperature and humidity surface fluxes are 0 .
06 K m s − and 2 . × − kg m (kg s) − , re-spectively. The surface shear stresses are computed in each grid cell using the Monin–Obukhovsimilarity theory. The simulations are run for 4 h. The shallow cumulus convection simulations follow the setup of the RICO case but do not includethe process of precipitation. The RICO conditions are chosen because convection is more vigorous9ompared to other cases of non-precipitating shallow convection, e.g., Siebesma et al. (2003), there-fore, it is expected to be a more stringent case. The initial θ and q t profiles have a mixed layer depthof 740 m and linearly vary above the mixed layer. Large-scale subsidence, moisture and humidityadvection, and a uniform clear sky radiative cooling are included in the simulations. The surfacefluxes are parameterized using bulk transfer coefficients and a constant sea surface temperature298 . Key parameterization-relevant boundary layer statistics are considered, because often LES of at-mospheric boundary layers is viewed as a reference model for Reynolds Averaged Navier–Stokes(RANS) turbulence closures. In cloud-free convection the depth of the boundary layer, z i ( t ), is de-fined as the height of the minimum of the buoyancy flux. In cloudy cases, the cloud-top height, z c ( t ),and cloud-base height, z b ( t ), are used as reference depths. Cloud cover, cc , is defined as the fractionof model columns with at least one model level with liquid water mixing ratio q l > − kg kg − .Cloud cover is sensitive to small fluctuations of the thermodynamic variables because regions withsmall amounts of water content are counted as cloudy columns.Turbulent fluxes are estimated in the LES as the sum of horizontally-averaged fluctuations andthe subgrid-scale stress. Fluctuations are denoted by primes, e.g., u ′ ( t, x, y, z ) = u ( t, x, y, z ) −h u ( t, x, y, z ) i x . Additionally, to increase the statistical sample, the turbulent flux is averaged overa short time interval, T = 0 . h ww i ( t, z ) = 1 T Z tt − T (cid:0) h ˜ w ′ ˜ w ′ i x + h τ i x (cid:1) d t (21)The turbulent kinetic energy does not include the subgrid-scale contribution, because in the Smagorin-sky model only the deviatoric stress is included in τ . The vertically-integrated turbulent kineticenergy VTKE( t ) = Z L z ρ ( z ) (cid:0) h ˜ u ′ ˜ u ′ i x + h ˜ v ′ ˜ v ′ i x + h ˜ w ′ ˜ w ′ i x (cid:1) d z, (22)provides a bulk measure of TKE in the boundary layer, and liquid water pathLWP( t ) = Z L z ρ ( z ) h q l i x d z, (23)a bulk measure of cloud liquid water content. The integrals at taken from the surface to the top ofthe computational domain. Turbulence is only present in the boundary layer. The buoyant bubble case simulations are a simplified model of convection. The simple configurationallows for fully-resolved simulations using constant viscosity/diffusivity coefficients, thus creatingan effective DNS. Figure 1 shows the evolution of the flow in the DNS case with the fourth-orderscheme, Case BD4g. The flow is initially driven by potential energy. The rising bubble in a flow10igure 1: Direct numerical simulation of a buoyant bubble (Case BD4g). Color contours show theevolution of total water mixing ratio. The colorbar units are kg kg − . Black contour correspondsto the saturation mixing ratio, denoting the cloud boundary. Second orderPSfrag replacements BD2fBD2gBD4fBD4g t (h) t (h) t (h) t (h) V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) t (h) t (h) t (h) t (h) V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 2: Comparison of the time evolution of vertically integrated turbulent kinetic energy(VTKE), liquid water path (LWP), cloud cover cc , and cloud base z b and height z c for the buoyantbubble DNS cases. Tow row panels correspond to the second-order scheme (Cases BD2f/g) andbottom row to the fourth-order scheme (Cases BD4f/g).11 Sfrag replacements BHfBHg t (h) t (h) t (h) t (h) V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 3: Time evolution of the vertically integrated turbulent kinetic energy (VTKE), liquidwater path (LWP) cloud cover, and cloud base and top heights for the high resolution (∆ x = 10 m)buoyant bubble LES in the fixed (BHf) and Galilean (BHg) frames. PSfrag replacements BLfBLg t (h) t (h) t (h) t (h) V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 4: Time evolution of the vertically integrated turbulent kinetic energy (VTKE), liquid waterpath (LWP) cloud cover, and cloud base and top heights for the low resolution (∆ x = 20 m) buoyantbubble LES in the fixed (BLf) and Galilean (BLg) frames.with mean shear creates a fairly complex flow that is not fully captured in the vertical planes shownin Fig. 1. At about t = 0 . z ≈ . t = 0 . t = 2280 s in Fig. 1.The preceding panel, t = 1320 s, shows the dissipation phase of the first cloud. The buoyantbubble DNS cases were set up such that the fourth-order scheme fully resolves the flow whereassecond-order scheme creates sufficiently large errors to excite the second term in (6). The timetraces of vertically integrated turbulent kinetic energy (VTKE), liquid water path (LWP), cloudcover, and cloud base z b and cloud top z c height of Cases BD2f/g and BD4f/g are shown Fig. 2.The results confirm that the model behaves as (6) predicts: the fourth-order results are Galileaninvariant whereas the second-order pair of solutions decorrelates. In Cases BD2g/f the truncationerror [angled brackets in (6)] is essentially the same. However, the coefficient u is larger in thefixed frame making the overall error large. In Cases BD4g/f the truncation error is sufficiently smalland does not significantly increase when multiplied by u . Case BD2g agrees with the fourth-orderresults., thus the error is in the fixed frame BD2f run as predicted by (6).The buoyant bubble LES results with ∆ x = 10 and 20 m are shown in Figs. 3 and 4, respectively.Results for both grid resolutions are not Galilean invariant. The only common pattern in all panelsof Figs. 3 and 4 is that differences appear after about t > .
25 h, when the flow develops rich12
Sfrag replacements D2fD2gD4fD4gD6fD6g t (h) t (h) V T K E ( k g s − ) z i ( k m )
00 11 22 33 44 0300600900120000.30.60.91.2
Figure 5: Time evolution of the boundary layer height z i and vertically integrated turbulent kineticenergy (VTKE) for the dry convective boundary layer cases.three-dimensional structure. There is no significant trend of the differences with respect to gridresolution. Overall, VTKE and LWP differences are larger for ∆ x = 10 m compared to the LESpair with ∆ x = 20 m. The reverse is observed for cloud cover and cloud-top height where overalldifferences are larger in the course grid runs. Figure 5 shows time traces of boundary layer height z i and VTKE. Figure 6 shows profiles averagedbetween t = 3 . u and v profiles. Thetemperature structure and entrainment rate (see time evolution of z i in Fig. 5) are nearly identical.The differences of VTKE in Fig. 5 are comparable to the random statistical variability (see AppendixB). Even though the results of the dry convection case are not identical between the two frames,as for instance the buoyant bubble DNS, any differences are relatively small. Figure 7 shows time traces of VTKE, LWP, cc , z b and z c for the shallow cumulus cases. ThreeLES pairs are shown corresponding to second-, fourth-, and sixth-order accurate schemes. Thetime traces in Fig. 7 correspond to one-hour moving averages starting at t = 1 . u with differences depending on the advection scheme order.The differences are larger when the second-order scheme is used and very small or negligible forthe sixth-order scheme. Moreover, the sensitivity depends on the flow statistic: LWP and theboundary layer depth (defined here as z c ) are less sensitive to the frame of reference and onlysimulations using the second-order scheme show significant differences. VTKE and cloud cover arethe most sensitive quantities. Cloud cover is very sensitive to horizontal fluctuations of small cloud13 PSfrag replacements z ( k m ) z ( k m ) DfDg u (m s − ) v (m s − ) θ (K)TKE (m s − ) h ww i (m s − ) h wθ v i ( × − m K s − )000 000 0.50.50.5 0.50.50.5 111 111 1.51.51.5 1.51.51.5 222 222 − − − − − − . − − . − − Figure 6: Dry convective boundary layer profiles of zonal wind u , meridional wind v , potential tem-perature θ , vertical velocity variance h ww i , resolved-scale turbulent kinetic energy, and buoyancyflux h wθ v i averaged in t = 3 . u .Profiles of the shallow cumulus runs averaged in t = 17 . u , v , θ l , and q t ) are negligible but turbulent fluxesand q l differ with respect to the frame of reference. Particularly TKE and h ww i exhibit significantdifferences in the two frames. With the exception of q l and h ww i near the cloud top, all Galileanframe profiles are in good agreement. In Figs. 7 and 8 the surface-fixed frame results differ andappear to converge to the Galilean frame results, which do not exhibit significant sensitivity to theadvection scheme.Larger differences are observed in TKE profiles of fixed frame LES, particularly in the cloudlayer. The amount of VTKE error for the second-order scheme is somewhat surprising, even afteraccounting of a larger error expectation in (6). The vertical velocity variance h ww i of fixed frameruns is also different in the cloud layer and the lower half of the mixed layer. These observationssuggest that discrepancies with respect to u are mostly found in the cloud layer and they are causedby the fluctuating character of the flow, because the mean profiles of the conserved thermodynamicvariables are similar in the two frames. Presently, LES results are compared only with respectto the frame of reference, therefore, errors related to the SGS model might still be present in the14 Sfrag replacements C2fC2gC4fC4gC6fC6g t (h) t (h) t (h) t (h) V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m )
00 00 66 66 1212 1212 1818 181818 0200400600800 0510152000.0750.1500.2150.3 00.751.502.153
Figure 7: Time evolution of the vertically integrated turbulent kinetic energy (VTKE), liquid waterpath (LWP) cloud cover, and cloud base and top heights for the shallow cumulus cases.Galilean frame results.An additional pair of LES using a modified condensation scheme (19) was carried out to assess ifthe observed differences are because of condensation/evaporation effects. Figure 9 shows time tracesfor the pair of LES with the modified condensation scheme (Cases CS4f/g). The time averagingprocedure used in traces of Fig. 7 is also followed in the traces shown in Fig. 9. The differences inVKTE with respect to the frame of reference are similar to the corresponding runs using the “all ornothing” condensation scheme, C4f/g. LWP and cloud cover in runs CSf/g increase with respectto runs C4f/g since additional partial condensation occurs in the modified scheme.Radial spectra computed on horizontal planes at the end of the simulations ( t = 18 h) are shown15 PSfrag replacements C2fC2gC4fC4gC6fC6g z ( k m ) z ( k m ) z ( k m ) u (m s − ) v (m s − ) θ l (K)TKE (m s − ) h ww i (m s − ) h wθ l i (m K s − ) q l (g kg − ) q t (g kg − ) h wq t i (m s − )000 000 000 111 111 111 222 222 222 333 333 3334 − − − − − − − − − − . − . − .
02 0 0.02 − . Figure 8: Shallow cumulus cases profiles of zonal wind u , meridional wind v , liquid water potentialtemperature θ l , total water mixing ratio q t , liquid water mixing ratio q l , vertical velocity variance h ww i , resolved-scale turbulent kinetic energy, temperature flux h wθ l i , and total water flux h wq t i .The profiles are time averaged in t = 17 . u and q t are shown at mid-height in the subcloud layer( z = 250 m) and at the middle of the cloud layer ( z = 1500 m). The frame of reference affects thesmall scales, corroborating the analysis and findings of Bernardini et al. (2013) in LES modeling.Spectra of the fixed frame LES exhibit a “pile up” of energy before the spectral roll-off at high16 Sfrag replacements t (h) t (h) t (h) t (h)CSfCSg V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 9: Time evolution of the vertically integrated turbulent kinetic energy (VTKE), liquid waterpath (LWP) cloud cover, and cloud base and top heights for the shallow cumulus simulations withthe modified saturation scheme in the fixed (CSf) and Galilean (CSg) frames.wavernumbers. Overall the results of Fig. 10 follow the observations of turbulent fluxes. The energypump is larger in the could layer and the effect is smaller when the sixth-order scheme is used. Theenergy pump is not only smaller when the high-order scheme is used but also confined to higherwavenumbers. Only the q t spectra at z = 250 m exhibit power-law scaling of about a decade witha somewhat shallower slope than − /
3. As will be discussed in the next section, the flow is notuniformly turbulent in the cloud layer, thus classical turbulent scaling is not expected.
The present results show that Galilean invariance in LES can be flow dependent. Well-resolvedDNS confirms that the error can be negligible, or at least controlled. Further, LES results with amodified condensation scheme suggest that the error is not a result of variations of the buoyancyforcing due to condensation and evaporation. The differences with respect to the frame of referencein “dry” (i.e., cloud-free) convective boundary layers are small and comparable to the range ofrandom statistical variability of the present results. In cloudy convection, the error depends on theorder of accuracy (resolving power) of the scheme and it is not uniformly distributed in the flowdomain. Most of the error is in the cloud layer.Earlier studies of Galilean invariance examined non-turbulent flows or continuously turbulentflows using DNS (e.g., Bernardini et al., 2013; Bihlo and Nave, 2014). Presently, we explore LES-SGS modeling of turbulent flows and, additionally, the cumulus convection cases have intermittentlyturbulent regions, i.e., in the cloud layer. The aforementioned points lead to the main question ofthe present study: can we formulate a mechanistic physical-space view of the observed Galileaninvariance error?We consider the skewness of the vertical velocity for a dry (Case D4g) and cumulus convection(Case C4g) in Fig. 11 as a simple measure of the large-scale flow anisotropy. In Fig. 11, height isnormalized with z i for both dry and cumulus convection, which scales z with the depth of the mixedlayer. As a consequence, in the dry convection case, turbulence is confined in the layer < . z/z i ,whereas in the shallow cumulus case, the boundary layer grows in time to reach about 2 z/z i atthe end of the run. In the cumulus case, z i scales the w skewness well for z/z i <
1. As expected(Heus and Jonker, 2008; Chinita et al., 2018), the vertical velocity distribution is positively skewedin the cloud layer because of the strong updrafts in the cloud cores. In the subcloud layer and in17
Sfrag replacements k r ∆ x/ πk r ∆ x/ π k r ∆ x/ πk r ∆ x/ π E uu E uu E qq E qq C2fC2gC6fC6g 10 − − − − − − − − − − − −
11 11 10 − − − − − − − − − − − − − − − − − − − − − − − − − z = 250 m z = 250 m z = 1500 m z = 1500 m − / − / − / − / Figure 10: Radial spectra of zonal wind (top row) and total water mixing ratio at mid-height in thesubcloud layer (left column) and about mid-height in the cloud layer. The second and sixth-ordershallow cumulus cases are shown using instantaneous horizontal planes at the end, t = 18 h, of thesimulation.dry convection cases, the positive bias of the w distribution is significantly less, implying a flowwith comparatively more symmetric structure.Figures 12 and 13 show conceptual mechanisms of how the flow structure can modulate initialdispersion errors to inhibit or allow their growth. In the dry convective case, Fig. 13, the flow inthe boundary layer is continuously turbulent (e.g., Chinita et al., 2018; Haghshenas and Mellado,2019). Thus, the SGS model is expected to rapidly dissipate any dispersive oscillations. In the18 Sfrag replacements D4g C4g4 h3 h2 h1 h 16 h12 h8 h4 h zz − i zz − i w skewness w skewness − − − . − . Figure 11: Vertical velocity skewness profiles at different times for the dry convection (D4g) andshallow cumulus (C4g) cases. The vertical axis is scaled by the height of the minimum buoyancyflux z i .cumulus cloud layer, there is no dissipation mechanism in the free troposphere and dispersiveoscillations can be long lived, c.f., Fig. 12 of Matheou and Dimotakis (2016). Also, as shown inFig. 11 for the present flow, the contrast between updrafts and downdrafts is significantly larger inthe cloud layer compared to the mixed layer. In the surface-fixed frame (Fig. 12) cloudy updraftsleave a trail of dispersive oscillations in the free troposphere. In the Galilean frame the developingturbulent cloud does not translate significantly on the grid and can engulf the spurious oscillations,which will then be dissipated by the action of the SGS. Figure 14 shows horizontal slices of q t atabout the middle of the cloud layer ( z = 1500 m) at t = 18 h for cases C2f and C2g. Dispersiveoscillations downwind of the clouds are present in the C2f case, whereas dispersive oscillationsare comparatively far fewer in C2g. The LES fields (Fig. 14) support the conceptual mechanismdepicted in Fig. 12: dispersive oscillations are most prominent in run C2f compared to C2g. The Galilean invariance properties of shallow convection LES are explored. In the past, a computa-tional domain translation velocity u was used to improve computational performance by allowinglarger time steps. Simulations carried out in a translating with the domain-mean wind frame ofreference completed in about half the time compared to simulations in the surface-fixed frame.Even though the equations of motion are Galilean invariant, i.e., they do not depend on u ,LES results have been observed to depend on u (Matheou et al., 2011; Wyant et al., 2018). Cen-tered fully-conservative finite difference schemes are used in the present simulations. For DNScases the analysis and prediction of Galilean invariance errors of Bernardini et al. (2013) are con-firmed. In LES, velocity and scalar spectra also corroborate the conclusions of Bernardini et al.19 alilean frame InitialLater updraft downdraft
Fixed frame
InitialLater dispersionerrorsinversion
Figure 12: Dispersion error growth in dry convection LES in the fixed and Galilean frames. Circlesrepresent dispersion errors. Larger error are depicted with larger circles.
Fixed frame
InitialLater
Galilean frame
InitialLater dispersionerrors
Figure 13: Dispersion error growth in cumulus convection LES in the fixed and Galilean frames.Circles represent dispersion errors. Larger error are depicted with larger circles.(2013). However, in LES the Galilean invariance error was found to strongly depend on the flowconfiguration. The error in the dry convection case is small. In the cumulus convection case, theerror depends on the resolving power of the scheme and it is mostly present in the cloud layercompared to the subcloud layer. The error significantly decreases as the order or accuracy of theadvection scheme is increased from second to sixth. The second-order scheme results in significantdiscrepancies between the Galilean and fixed-frame LES, with differences between the two frames20 ixed frame (C2f) Galilean frame (C2g) q t (g kg –1 ) Figure 14: Total water mixing ratio q t contours on horizontal planes at z = 1 . t = 18 h for the shallow cumulus case and second-order advectiondiscretization. The left panel corresponds to the fixed frame simulation (C2f) and the right to theGalilean frame (C2g). Black contour corresponds to the saturation mixing ratio, denoting the cloudboundary. Black arrows show some instances of dispersive oscillations downwind of the clouds.becoming negligible when the sixth-order accurate advection scheme was used.The error mostly affects second-order statistics in cumulus convection, including liquid waterprofiles and liquid water path, and, to a lesser extent, boundary-layer growth rates. In the presentsimulations, the most sensitive quantity is the turbulent kinetic energy (TKE). The vertical integralof TKE was found to differ by as much as 20% between surface-fixed frame and Galilean frameLES. The present results suggest that biases in finite difference dispersion errors can be amplified bylarge-scale flow asymmetries, such as strong updrafts rising in the non-turbulent free troposphere incumulus-cloud layers. The strong dissipative action of the SGS model in the continuously-turbulentmixed layer in dry convection can control the error growth.One of the most interesting findings is that using a second-order discretization in the properGalilean frame can result in comparable accuracy as a high-order scheme in the surface-fixed frame.The aforementioned observation is a likely explanation of a long standing conundrum in LES of con-vection in the atmospheric boundary layer: atmospheric LES results have been sufficiently accurateeven when low-order schemes are used (e.g., Siebesma et al., 2003), whereas in general computa-tional fluid dynamics modeling many high-order schemes have been developed to improve accuracy(e.g., Lele, 1992; Kravchenko and Moin, 1997; Laizet and Lamballais, 2009; Pirozzoli, 2011). Un-fortunately, the frame of reference of the LES is not a quantity that is often reported in theatmospheric boundary layer literature. In summary, it appears that a technique primarily used forcomputational performance gain can have significant accuracy advantages as well.Presently, we use a translating frame where the domain-volume-mean wind is nearly zero.However, we do not expect that the global error is necessarily minimum in this frame (i.e., this21hoice may not be globally optimal). Furthermore, the choice of the of the SGS model can affectthe error in the LES as has been shown in several past studies (Ghosal, 1996; Vreman et al.,1996; Kravchenko and Moin, 1997; Fedioun et al., 2001; Chow and Moin, 2003; Geurts, 2009). Thesimple Smagorinsky–Lilly model is presently used because of its prevalence in LES of atmosphericboundary layers (e.g., Stevens et al., 2001; Siebesma et al., 2003; vanZanten et al., 2011). Acknowledgments
This work was funded by the National Science Foundation via Grant NSF-AGS-1916619. Theresearch presented in this paper was supported by the systems, services, and capabilities providedby the University of Connecticut High Performance Computing (HPC) facility.22
Sfrag replacements t (h) t (h) t (h) t (h)CSfCSg V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 15: Time traces of vertically integrated turbulent kinetic energy, liquid water path, cloudcover cc , cloud base z b and cloud top height z c for four buoyant bubble DNS with the second-orderscheme (Case BD2f) with CFL = 0 .
2, 0.4, 0.8, and 1.2. The results do not depend on CFL, thusindividual lines are not labeled. In each panel, four coinciding lines are plotted.
PSfrag replacements t (h) t (h) t (h) t (h)CSfCSg V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 16: Time traces of vertically integrated turbulent kinetic energy, liquid water path, cloudcover cc , cloud base z b and cloud top height z c for four shallow cumulus LES with CFL = 0 . t = 1 . A Dependence on numerical time step interval
DNS and LES results were not found to depend on the CFL number, or equivalently, on the timestep length ∆ t . Sensitivity to CFL for the buoyant bubble DNS with the second-order scheme isshown in Fig. 15. The time traces of VTKE, LWP, cc , z b and z c for four BF2 runs with CFL = 0 . t = 1 . cc , z b and z c of four simulations with CFL = 0 .
2, 0.4, 0.8, and 1.2are plotted. The moving average only removes the short-time variability of the turbulent flow anddoes not effect the nature of the comparison. No sensitivity to CFL is observed in Fig. 16 since alltraces are within the statistical variability of the present simulations.23
Statistical variability
Two ten-member ensembles are carried out to estimate the statistical variability of the LES resultsin dry convection and shallow cumulus cases. Because of the finite computational domain, acomplete sample of the flow states is not accomplished and instantaneous horizontal averages arenot fully converged statistically. Ensemble simulations were initialized by applying different randomtemperature and humidity perturbations in the LES near the surface.Figure 17 shows the band of VTKE variability of the dry convection Case D4f ensemble. Asthe boundary layer deepens, the convection cells become larger and fewer in the fixed domain size.Thus, the sampling of the flow declines with time and the band of VTKE variability widens. After t = 3 h, VTKE is uncertain by about 50 kg m − or 4%.Figure 18 shows time traces of VTKE, LWP, cc , z b and z c for a ten-member Case C4f ensem-ble. Because the LES computational domain is somewhat large, about 16 times the depth of theboundary layer, the statistical variability of the quantities in Fig. 18 is relatively small. PSfrag replacements t (h) V T K E ( k g s − ) z i (km) 0 1 2 3 40300600900120000.30.60.91.2 Figure 17: Spread of vertically integrated turbulent kinetic energy vs time for the ten-member-ensemble dry convective boundary layer (case Df).
PSfrag replacements t (h) t (h) t (h) t (h)CSfCSg V T K E ( k g s − ) L W P ( g m − ) cc ( f r a c t i o n ) z b , z c ( k m ) Figure 18: Spread of vertically integrated turbulent kinetic energy, cloud cover cc , cloud base z b and cloud top height z c vs time for the ten-member-ensemble shallow cumulus case C4f. For eachsimulation, a one-hour moving average starting at t = 1 . eferences Basu, S. and Port´e-Agel, F. (2006). Large-eddy simulation of stably stratified atmospheric boundarylayer turbulence: a scale-dependent dynamic modeling approach.
J. Atmos. Sci. , 63:2074–2091.Bernardini, M., Pirozzoli, S., Quadrio, M., and Orlandi, P. (2013). Turbulent channel flow simula-tions in convecting reference frames.
J. Comput. Phys. , 232(1):1–6.Bihlo, A. and Nave, J.-C. (2014). Convecting reference frames and invariant numerical models.
J.Comput. Phys. , 272:656–663.Bony, S. and Dufresne, J.-L. (2005). Marine boundary layer clouds at the heart of tropical cloudfeedback uncertainties in climate models.
Geophys. Res. Lett. , 32:L20806.Chinita, M. J., Matheou, G., and Teixeira, J. (2018). A joint probability density-based decompo-sition of turbulence in the atmospheric boundary layer.
Mon. Weather Rev. , 146(2):503–523.Chow, F. K. and Moin, P. (2003). A further study of numerical errors in large-eddy simulations.
J. Comput. Phys. , 184(2):366–380.Couvreux, F., Bazile, E., Rodier, Q., Maronga, B., Matheou, G., Chinita, M. J., Edwards, J.,Stratum, B. J. H. V., Heerwaarden, C. C. V., Huang, J., Moene, A. F., Cheng, A., Fuka, V.,Basu, S., Bou-Zeid, E., Canut, G., and Vignon, E. (2020). The GABLS4 experiment: intercom-parison of large-eddy simulation models of the antarctic boundary layer challenged by very stablestratification.
Boundary-Layer Meteorol.
Deardorff, J. W. (1972). Numerical investigation of neutral and unstable planetary boundary layers.
J. Atmos. Sci. , 29(1):91–115.Fedioun, I., Lardjane, N., and G¨okalp, I. (2001). Revisiting numerical errors in direct and largeeddy simulations of turbulence: physical and spectral spaces analysis.
Journal of ComputationalPhysics , 174(2):816–851.Fuka, V. and Brechler, J. (2011). Large eddy simulation of the stable boundary layer. In
Finitevolumes for complex applications VI problems & perspectives , pages 485–493. Springer.Geurts, B. J. (2009). Analysis of errors occurring in large eddy simulation.
Philosophical Trans-actions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences ,367(1899):2873–2883.Ghosal, S. (1996). An analysis of numerical errors in large-eddy simulations of turbulence.
J.Comput. Phys. , 125:187–206.Grabowski, W. W. and Smolarkiewicz, P. K. (1990). Monotone finite-difference approximations tothe advection–condensation problem.
Mon. Weather Rev. , 118(10):2082–2097.Haghshenas, A. and Mellado, J. P. (2019). Characterization of wind-shear effects on entrainmentin a convective boundary layer.
J. Fluid Mech. , 858:145–183.25eus, T., Heerwaarden, C. C. V., Jonker, H. J. J., Siebesma, A. P., Axelsen, S., Dries, K. V. D.,Geoffroy, O., Moene, A. F., Pino, D., Roode, S. R. D., and de Arellano, J. V.-G. (2010). Formu-lation of the Dutch atmospheric large-eddy simulation (DALES) and overview of its applications.
Geosci. Model Dev. , 3:415–444.Heus, T. and Jonker, H. J. J. (2008). Subsiding shells around shallow cumulus clouds.
J. Atmos.Sci. , 65:1003–1018.Huang, J. and Bou-Zeid, E. (2013). Turbulence and vertical fluxes in the stable atmosphericboundary-layer. Part I: A large-eddy simulation study.
J. Atmos. Sci. , 70:1513–1527.Inoue, M., Matheou, G., and Teixeira, J. (2014). LES of a spatially developing atmospheric bound-ary layer: Application of a fringe method for the stratocumulus to shallow cumulus cloud tran-sition.
Mon. Weather Rev. , 142(9):3418–3424.Jongaramrungruang, S., Frankenberg, C., Matheou, G., Thorpe, A. K., Thompson, D. R., Kuai,L., and Duren, R. M. (2019). Towards accurate methane point-source quantification from high-resolution 2-D plume imagery.
Atmospheric Measurement Techniques , 12(12).Khairoutdinov, M. F. and Randall, D. A. (2003). Cloud resolving modeling of the arm summer 1997iop: Model formulation, results, uncertainties, and sensitivities.
J. Atmos. Sci. , 60(4):607–625.Klein, S. A., Hall, A., Norris, J. R., and Pincus, R. (2017). Low-cloud feedbacks from cloud-controlling factors: a review. In
Shallow Clouds, Water Vapor, Circulation, and Climate Sensi-tivity , pages 135–157. Springer.Kravchenko, A. G. and Moin, P. (1997). On the effect of numerical errors in large eddy simulationsof turbulent flows.
J. Comput. Phys. , 131(2):310–322.Lac, C., Chaboureau, J.-P., Masson, V., Pinty, J.-P., Tulet, P., Escobar, J., Leriche, M., Barthe,C., Aouizerats, B., Augros, C., et al. (2018). Overview of the meso-nh model version 5.4 and itsapplications.
Geosci. Model Dev. , 11(5):1929.Laizet, S. and Lamballais, E. (2009). High-order compact schemes for incompressible flows: Asimple and efficient method with quasi-spectral accuracy.
J. Comput. Phys. , 228(16):5989–6015.Lele, S. K. (1992). Compact finite difference schemes with spectral-like resolution.
J. Comput.Phys. , 103:16–42.Lilly, D. K. (1962). On the numerical simulation of buoyant convection.
Tellus , 14(2):148–172.Lilly, D. K. (1966).
On the application of the eddy viscosity concept in the inertial sub-range ofturbulence , volume 123. National Center for Atmospheric Research.Lilly, D. K. (1967). The representation of small-scale turbulence in numerical simulation experi-ments. In
Proc. IBM Sci. Computing Symp. Environmental Sci , pages 195–210.Lomax, H., Pulliam, T. H., and Zingg, D. W. (2003).
Fundamentals of Computational FluidDynamics . Scientific Computation. Springer. 26argolin, L. G., Smolarkiewicz, P. K., and Sorbjan, Z. (1999). Large-eddy simulations of convectiveboundary layers using nonoscillatory differencing.
Physica D , 133(1):390–397.Maronga, B., Gryschka, M., Heinze, R., Hoffmann, F., Kanani-S¨uhring, F., Keck, M., Ketelsen, K.,Letzel, M. O., S¨uhring, M., and Raasch, S. (2015). The parallelized large-eddy simulation model(PALM) version 4.0 for atmospheric and oceanic flows: model formulation, recent developments,and future perspectives.
Geosci. Model Dev. , 8:2515–2551.Matheou, G. (2016). Numerical discretization and subgrid-scale model effects on large-eddy simu-lations of a stable boundary layer.
Q. J. R. Meteorol. Soc. , 142:3050–3062.Matheou, G. (2018). Turbulence structure in a stratocumulus cloud.
Atmosphere , 9(10):392.Matheou, G. and Bowman, K. W. (2016). A recycling method for the large-eddy simulation ofplumes in the atmospheric boundary layer.
Environmental Fluid Mechanics , 16:69–85.Matheou, G. and Chung, D. (2014). Large-eddy simulation of stratified turbulence. Part II: Ap-plication of the stretched-vortex model to the atmospheric boundary layer.
J. Atmos. Sci. ,71(12):4439–4460.Matheou, G., Chung, D., Nuijens, L., Stevens, B., and Teixeira, J. (2011). On the fidelity of large-eddy simulation of shallow precipitating cumulus convection.
Mon. Weather Rev. , 139:2918–2939.Matheou, G. and Dimotakis, P. E. (2016). Scalar excursions in large-eddy simulations.
J. Comput.Phys. , 327(2):97–120.Matheou, G. and Teixeira, J. (2019). Sensitivity to physical and numerical aspects of large-eddysimulation of stratocumulus.
Mon. Weather Rev. , 147:2621–2639.Morinishi, Y., Lund, T. S., Vasilyev, O. V., and Moin, P. (1998). Fully conservative higher orderfinite difference schemes for incompressible flow.
J. Comput. Phys. , 143(1):90–124.Neggers, R. A. J., Koehler, M., and Beljaars, A. C. M. (2009). A dual mass flux framework forboundary layer convection. Part I: Transport.
J. Atmos. Sci. , 66:1465–1487.Oberlack, M. (1997). Invariant modeling in large-eddy simulation of turbulence. In
Annual ResearchBriefs , pages 3–22. Stanford University.Pirozzoli, S. (2011). Numerical methods for high-speed flows.
Annu. Rev. Fluid Mech. , 43:163–194.Pope, S. B. (2004). Ten questions concerning the large-eddy simulation of turbulent flows.
NewJournal of Physics , 6:35.Rauber, R. M., Stevens, B., Ochs III, H. T., Knight, C., Albrecht, B. A., Blyth, A. M., Fairall,C. W., Jensen, J. B., Lasher-Trapp, S. G., Mayol-Bracero, O. L., Vali, G., Anderson, J. R.,Baker, B. A., Bandy, A. R., Burnet, E., Brenguier, J. L., Brewer, W. A., Brown, P. R. A.,Chuang, P., Cotton, W. R., Girolamo, L. D., Geerts, B., Gerber, H., Goke, S., Gomes, L.,Heikes, B. G., Hudson, J. G., Kollias, P., Lawson, R. P., Krueger, S. K., Lenschow, D. H.,Nuijens, L., O’Sullivan, D. W., Rilling, R. A., Rogers, D. C., Siebesma, A. P., Snodgrass, E.,Stith, J. L., Thornton, D. C., Tucker, S., Twohy, C. H., and Zuidema, P. (2007). Rain in shallowcumulus over the ocean: The RICO campaign.
Bull. Amer. Meteor. Soc. , 88:1912–1928.27ieck, M., Nuijens, L., and Stevens, B. (2012). Marine boundary layer cloud feedbacks in a constantrelative humidity atmosphere.
J. Atmos. Sci. , 69(8):2538–2550.Schalkwijk, J., Jonker, H. J. J., Siebesma, A. P., and Meijgaard, E. V. (2015). Weather forecastingusing GPU-based large-eddy simulations.
Bull. Amer. Meteor. Soc. , 96(5):715–723.Seifert, A. and Heus, T. (2013). Large-eddy simulation of organized precipitating trade windcumulus clouds.
Atmos. Chem. Phys. , 13:5631–5645.Siebesma, A. P., 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.(2003). A large eddy simulation intercomparison study of shallow cumulus convection.
J. Atmos.Sci. , 60:1201–1219.Siebesma, A. P. and Holtslag, A. A. M. (1996). Model impacts of entrainment and detrainmentrates in shallow cumulus convection.
J. Atmos. Sci. , 53:2354–2364.Siebesma, A. P., Soares, P. M. M., and Teixeira, J. (2007). A combined eddy-diffusivity mass-fluxapproach for the convective boundary layer.
J. Atmos. Sci. , 64:1230–1248.Smagorinsky, J. (1963). General circulation experiments with the primitive equations. I. The basicexperiment.
Mon. Weather Rev. , 91:99–164.Sommeria, G. (1976). Three-dimensional simulation of turbulent processes in an undisturbed tradewind boundary-layer.
J. Atmos. Sci. , 33:216–241.Spalart, P. R., Moser, R. D., and Rogers, M. M. (1991). Spectral methods for the Navier–Stokesequations with one infinite and two periodic directions.
J. Comput. Phys. , 96(2):297–324.Springel, V. (2010). E pur si muove: Galilean-invariant cosmological hydrodynamical simulationson a moving mesh.
Mon. Not. R. Astron. Soc. , 401(2):791–851.Stevens, B., Ackerman, A. S., Albrecht, B. A., Brown, A. R., Chlond, A., Cuxart, J., Duynkerke,P. G., Lewellen, D. C., Macvean, M. K., Neggers, R. A. J., Sanchez, E., Siebesma, A. P., andStevens, D. E. (2001). Simulations of trade wind cumuli under a strong inversion.
J. Atmos.Sci. , 58:1870–1891.Stevens, B., Moeng, C.-H., Ackerman, A. S., Bretherton, C. S., Chlond, A., De Roode, S., Edwards,J., Golaz, J.-C., Jiang, H. L., Khairoutdinov, M., Kirkpatrick, M. P., Lewellen, D. C., Lock, A.,Muller, F., Stevens, D. E., Whelan, E., and Zhu, P. (2005). Evaluation of large-eddy simulationsvia observations of nocturnal marine stratocumulus.
Mon. Weather Rev. , 133:1443–1462.Sullivan, P. G. and Patton, E. G. (2011). The effect of mesh resolution on convective boundary layerstatistics and structures generated by large-eddy simulation.
J. Atmos. Sci. , 68(10):2395–2415.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. (2008). Parameterization of the atmospheric boundary layer: A view from just above theinversion.
Bull. Amer. Meteor. Soc. , 89:453–458.28horpe, A. K., Frankenberg, C., Green, R. O., Thompson, D. R., Aubrey, A. D., P.Mouroulis,Eastwood, M. L., and Matheou, G. (2016). The Airborne Methane Plume Spectrometer (AMPS):Quantitative imaging of methane plumes in real time. In , pages1–14.van Heerwaarden, C., Van Stratum, B. J., Heus, T., Gibbs, J. A., Fedorovich, E., and Mellado, J.-P.(2017). MicroHH 1.0: a computational fluid dynamics code for direct numerical simulation andlarge-eddy simulation of atmospheric boundary layer flows.
Geosci. Model Dev. , 10:3145–3165.vanZanten, M. C., Stevens, B., Nuijens, L., Siebesma, A. P., Ackerman, A., Bogenschutz, P.,Burnet, F., Cheng, A., Couvreux, F., Jiang, H., Khairoutdinov, M., Lewellen, D. S., Noda, A.,Mechem, D., Shipway, B., Slawinska, J., and Wang, S. (2011). Controls on precipitation andcloudiness in simulations of shallow fair-weather cumulus.
J. Adv. Model. Earth Syst. , 3:Art.M06001.Vial, J., Bony, S., Stevens, B., and Vogel, R. (2017). Mechanisms and model diversity of trade-windshallow cumulus cloud feedbacks: A review. In
Shallow Clouds, Water Vapor, Circulation, andClimate Sensitivity , pages 159–181. Springer.Vreman, B., Geurts, B., and Kuerten, H. (1996). Comparison of numerical schemes in large-eddysimulation of the temporal mixing layer.
Int. J. Numer. Methods Fluids , 22(4):297–311.Witek, M. L., Teixeira, J., and Matheou, G. (2011). An eddy diffusivity–mass flux approach tothe vertical transport of turbulent kinetic energy in convective boundary layers.
J. Atmos. Sci. ,68(10):2385–2394.Wyant, M. C., Bretherton, C. S., and Blossey, P. N. (2018). The sensitivity of numerical simulationsof cloud-topped boundary layers to cross-grid flow.
J. Adv. Model. Earth Syst. , 10(2):466–480.Zelinka, M. D., Randall, D. A., Webb, M. J., and Klein, S. A. (2017). Clearing clouds of uncertainty.