Large-eddy simulations with ClimateMachine v0.2.0: a new open-source code for atmospheric simulations on GPUs and CPUs

We introduce ClimateMachine, a new open-source atmosphere modeling framework which uses the Julia language and is designed to be scalable on central processing units (CPUs) and graphics processing units (GPUs). ClimateMachine uses a common framework both for coarser-resolution global simulations and for high-resolution, limited-area large-eddy simulations (LESs). Here, we demonstrate the LES configuration of the atmosphere model in canonical benchmark cases and atmospheric flows using a total energy-conserving nodal discontinuous Galerkin (DG) discretization of the governing equations. Resolution dependence, conservation characteristics, and scaling metrics are examined in comparison with existing LES codes. They demonstrate the utility of ClimateMachine as a modeling tool for limited-area LES flow configurations.

1 Introduction

Hybrid computer architectures and the need to exploit the power of graphics processing units (GPUs) are increasingly driving developments in atmosphere and climate modeling (e.g., Schalkwijk et al., 2012; Palmer, 2014; Schalkwijk et al., 2015; Marras et al., 2015; Abdi et al., 2017b, a; Fuhrer et al., 2018; Schär et al., 2020). The sheer computing power available on modern hardware architectures presents opportunities to accelerate atmosphere and climate modeling. However, exploiting this computing power requires re-coding atmosphere and climate models to an extent not seen in decades, and portable performance and scaling across different platforms remain difficult to achieve (Fuhrer et al., 2014; Balaji, 2021).

In this paper, we introduce ClimateMachine, a new open-source atmosphere model written in the Julia programming language (Bezanson et al., 2017) using the Julia Message Passing Interface (MPI) to provide a computational framework that is portable across CPU and GPU architectures. The use of Julia aims to increase accessibility and utility of ClimateMachine as a simulation tool. Julia, developed with the aim of bridging the gap between productivity (e.g., Python) and performance (e.g., Fortran) languages (Bezanson et al., 2018), is a strong candidate for an application such as atmospheric large-eddy simulations (LESs). Features such as package composability, GPU programming tools, and interactive programming interfaces contribute to its usability in scientific computing. Compilation time, incompatibilities between Julia MPI programs and interactive read–eval–print loops (REPLs), and the discovery of code issues during runtime due to Julia's dynamic typing are general drawbacks for new users of this language. Despite these drawbacks, a recent assessment by Lin and McIntosh-Smith (2021) further supports the suitability of Julia as a high-performance computing language. The atmospheric model presented here is designed to be usable across a range of physical process scales from large-eddy simulations (LESs) with meter-scale resolution to global circulation models (GCMs) with horizontal resolutions of tens of kilometers, as in a few other recent models (Dipankar et al., 2015). We focus on the LES configuration of ClimateMachine in this paper.

Since the pioneering work on turbulence in stratified flows by Lilly (1962), on horizontal subgrid-scale dissipation in GCMs by Smagorinsky (1963), and the subsequent extension of such models to LESs (Lilly, 1966), several models have been developed to improve the ability of LESs to model atmospheric turbulence: from the extensive work by Deardorff in the 1970s and 1980s (Deardorff, 1970, 1974, 1976, 1980), Moeng in the 1980s and beyond (Moeng, 1984; Moeng and Wyngaard, 1988; Sullivan et al., 1994; Moeng et al., 2003), to Stevens, Teixeira, Mellado, and others in the last 2 decades (Stevens et al., 2003, 2005; Savic-Jovcic and Stevens, 2008; Matheou et al., 2011; Pressel et al., 2015; Matheou, 2016; Matheou and Teixeira, 2019; Mellado, 2017; Mellado et al., 2018). LES results in canonical flows are sensitive to the fine details of the equations used to represent the flow dynamics, the viscous dissipation, the thermodynamics, and the numerical methods used to solve them (Ghosal, 1996; Chow and Moin, 2003; Kurowski et al., 2014), especially in the case of cloud simulations (Stevens et al., 2005; Siebesma et al., 2003; Schalkwijk et al., 2012, 2015; Schneider et al., 2019; Pressel et al., 2015, 2017).

One distinguishing aspect of the ClimateMachine LES is that it uses a nodal discontinuous Galerkin (DG) formulation to approximate the Navier–Stokes equations for compressible flow (Giraldo et al., 2002; Hesthaven and Warburton, 2008a; Giraldo and Restelli, 2008; Kopriva, 2009; Kelly and Giraldo, 2012; Giraldo, 2020). The DG method is a spectral-element generalization of finite-volume methods. It lends itself well to modern high-performance computing architectures because its communication overhead is low, enabling scaling on many-core processors including GPUs (Abdi et al., 2017b). Another important consideration within ClimateMachine is the use of total energy of moist air as a prognostic variable, ensuring energetic consistency of the simulations. We demonstrate that the ClimateMachine LES can be successfully used to simulate canonical LES benchmarks, including simulations of flows over mountains and different cloud and boundary layer regimes (e.g., Straka et al., 1993; Schär et al., 2002; Stevens et al., 2005).

In what follows, we describe the conceptual and numerical foundations and governing equations of ClimateMachine and demonstrate the model in a set of standard two- and three-dimensional benchmark simulations. Section 2 begins by highlighting the governing equations. Their numerical approximation through the DG representation is described in Sect. 3. Section 4 presents subgrid-scale models used in the LES to represent under-resolved flow physics, with results from key benchmarks presented in Sect. 5. Conservation properties are examined in Sect. 6, and performance on CPU and GPU hardware is described in Sects. 7 and 8, respectively. Section 9 contains closing remarks. Additional details about the model, boundary conditions, statistical definitions, and computer hardware are summarized in the Appendices.

2 Governing equations

2.1 Working fluid

The working fluid of the atmosphere model is moist, potentially cloudy air, considered to be an ideal mixture of dry air, water vapor, and condensed water (liquid and ice) in clouds. Dry air and water vapor are taken to be ideal gases. The specific volume of the cloud condensate is neglected relative to that of the gas phases (it is a factor of 10³ less than that of the gas phases). All gas phases are assumed to have the same temperature and are advected with the same velocity $\mathit{u}=(u,v,w{)}^T$. Cloud condensate is assumed to sediment relative to gaseous phases slowly enough to be in thermal equilibrium with the surrounding fluid.

The density of the moist air is denoted by ρ. We use the following notation for the mass fractions of the moist air mixture (mass of a constituent divided by the total mass of the working fluid).

• q_d: dry air mass fraction

• q_v: water vapor specific humidity

• q_l: liquid water specific humidity

• q_i: ice specific humidity

• $q_{\mathrm{c}}=q_{\mathrm{l}}+q_{\mathrm{i}}$: condensate specific humidity

• $q_{\mathrm{t}}=q_{\mathrm{v}}+q_{\mathrm{c}}$: total specific humidity

Because this enumerates all constituents of the working fluid, we have $q_{\mathrm{t}}+q_{\mathrm{d}}=\mathrm{1}$. In Earth's atmosphere, the water vapor specific humidity q_v dominates the total specific humidity q_t and is usually 𝒪(10⁻²) or smaller; the condensate specific humidity is typically 𝒪(10⁻⁴). Hence, water is a trace constituent of the atmosphere, and only a small fraction of atmospheric water is in condensed phases. The working fluid pressure is the sum of the partial pressures of dry air and water vapor such that $p=\mathit{\rho }(R_{\mathrm{d}}{q}_{\mathrm{d}}+R_{\mathrm{v}}{q}_{\mathrm{v}})T$, where R_d is the specific gas constant of dry air, and R_v is the specific gas constant of water vapor.

2.2 Mass balance

Moist air mass satisfies the conservation equation:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{u}\right)=\mathit{\rho }{\widehat{\mathcal{S}}}_{q_{\mathrm{t}}}.\end{array}$$

Moist air mass is not exactly conserved where precipitation forms, sublimates, or evaporates, where water diffuses, or where condensate sediments relative to the gas phases (Bott, 2008; Romps, 2008). The right-hand side involves the local source and/or sink of water mass ${\widehat{\mathcal{S}}}_{q_{\mathrm{t}}}$ owing to such nonconservative processes, which we take into account although it is small because water is a trace constituent of the atmosphere.

2.3 Total water balance

Total water satisfies the balance equation:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathit{\rho }{q}_{\mathrm{t}}\right)}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }{q}_{\mathrm{t}}\mathit{u}\right)& =\mathit{\rho }{\mathcal{S}}_{q_{\mathrm{t}}}-\mathrm{\nabla }\cdot \left(\mathit{\rho }{\mathit{d}}_{q_{\mathrm{t}}}\right)\\ & +\mathrm{\nabla }\cdot \left(\mathit{\rho }{q}_{\mathrm{c}}{w}_{\mathrm{c}}\widehat{\mathit{k}}\right)\equiv \mathit{\rho }{\widehat{\mathcal{S}}}_{q_{\mathrm{t}}}.\end{array}\end{array}$$

Here, the source–sink ${\mathcal{S}}_{q_{\mathrm{t}}}$ arises from evaporation or sublimation of precipitation and formation of precipitation. Diffusive fluxes of moisture are captured by ${\mathit{d}}_{q_{\mathrm{t}}}$. The effective sedimentation velocity of cloud condensate w_c is defined such that

$$\begin{array}{}\text{(3)}& q_{\mathrm{c}}{w}_{\mathrm{c}}=q_{\mathrm{l}}{w}_{\mathrm{l}}+q_{\mathrm{i}}{w}_{\mathrm{i}},\end{array}$$

with w_l and w_i defined to be positive downward ($\widehat{\mathit{k}}$ being the upward-pointing unit vector). The right-hand side $\mathit{\rho }{\widehat{\mathcal{S}}}_{q_{\mathrm{t}}}$ of the total water balance equation is the same as the right-hand side of the mass balance in Eq. (1).

2.4 Momentum balance

The coordinate-independent form of the conservation law for momentum is

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathit{\rho }\mathit{u}\right)}{\partial t}}& +\mathrm{\nabla }\cdot \left[\mathit{\rho }\mathit{u}\otimes \mathit{u}+(p-p_{\mathrm{r}}){\mathit{I}}_{\mathrm{3}}\right]=-(\mathit{\rho }-{\mathit{\rho }}_{\mathrm{r}})\mathbf{\nabla }\mathrm{\Phi }\\ & -\mathrm{2}\mathit{\Omega }\times \mathit{\rho }\mathit{u}-\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{\tau }\right)-\mathrm{\nabla }\cdot \left({\mathit{d}}_{q_{\mathrm{t}}}\otimes \mathit{\rho }\mathit{u}\right)\\ & +\mathrm{\nabla }\cdot \left(q_{\mathrm{c}}{w}_{\mathrm{c}}\widehat{\mathit{k}}\otimes \mathit{\rho }\mathit{u}\right)+\mathit{\rho }{\mathit{F}}_{\mathit{u}},\end{array}\end{array}$$

where I₃ is the rank-3 identity matrix, Φ is the effective gravitational potential including centrifugal accelerations, τ is a viscous and/or subgrid-scale (SGS) momentum flux tensor, and F_u (typically with ${\mathit{F}}_u\cdot \mathit{u}<\mathrm{0}$ so that F_u represents a momentum sink) is any other drag force per unit mass that may be applied, for example, at the lower boundary. The term involving the planetary angular velocity Ω accounts for Coriolis forces. To improve numerical stability, we have factored out a reference state with a pressure p_r(z) and density ρ_r(z) that depend only on altitude z and are in hydrostatic balance so that they satisfy

$$\mathbf{\nabla }{p}_{\mathrm{r}}=-{\mathit{\rho }}_{\mathrm{r}}\mathbf{\nabla }\mathrm{\Phi }.$$

The tensor involving the diffusive flux ${\mathit{d}}_{q_{\mathrm{t}}}$ of water on the right-hand side of Eq. (4) represents the momentum flux carried by water that is diffusing; this term is usually very small, but we take it into account.

2.5 Energy balance

The specification of a thermodynamic or energy conservation equation closes the equations of motion for the working fluid. We use the total specific energy, e^tot, as the prognostic variable. Total energy is conserved in reversible moist processes such as phase transitions of water.

Total energy satisfies the conservation law (Romps, 2008; Bott, 2008):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{\partial \left(\mathit{\rho }{e}^{\mathrm{tot}}\right)}{\partial t}}& +\mathrm{\nabla }\cdot \left((\mathit{\rho }{e}^{\mathrm{tot}}+p)\mathit{u}\right)=-\mathrm{\nabla }\cdot \left(\mathit{\rho }{\mathit{F}}_{\mathrm{R}}\right)\\ & -\mathrm{\nabla }\cdot \left[\mathit{\rho }\right(\mathit{J}+\mathit{D}\left)\right]+\mathit{\rho }Q+\mathrm{\nabla }\cdot \left(\mathit{\rho }{W}_{\mathrm{c}}\widehat{\mathit{k}}\right)\\ & -\mathrm{\nabla }\cdot (\mathit{u}\cdot \mathit{\rho }\mathit{\tau }\mathbf{)}\\ & -\sum _{\mathrm{j}\in \mathit{\{}\mathrm{v},\mathrm{l},\mathrm{i}\mathit{\}}}(I_{\mathrm{j}}+\mathrm{\Phi })\mathit{\rho }C(q_{\mathrm{j}}\to q_{\mathrm{p}})-M,\end{array}\end{array}$$

where the total specific energy e^tot is defined by

$$\begin{array}{}\text{(6)}& e^{\mathrm{tot}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\Vert \mathit{u}{\Vert }^{\mathrm{2}}+\mathrm{\Phi }+I.\end{array}$$

The constituents (dry air and moisture components) here are assumed to be moving with the same velocity u (that is, we neglect, as is common, the diffusive and sedimentation fluxes of water in the kinetic energy). The constituents are also assumed to be in thermal equilibrium at the same temperature T so that the specific internal energy of moist air is the weighted sum of the specific energies of the constituents dry air (I_d), water vapor (I_v), liquid water (I_l), and ice (I_i):

$$\begin{array}{}\text{(7)}& I(T,q)=(\mathrm{1}-q_{\mathrm{t}})I_{\mathrm{d}}\left(T\right)+q_{\mathrm{v}}{I}_{\mathrm{v}}\left(T\right)+q_{\mathrm{l}}{I}_{\mathrm{l}}\left(T\right)+q_{\mathrm{i}}{I}_{\mathrm{i}}\left(T\right),\end{array}$$

with

$$\begin{array}{}\text{(8a)}& {\displaystyle }{I}_{\mathrm{d}}\left(T\right)& {\displaystyle }=c_{\mathrm{vd}}(T-T_{\mathrm{0}}),\text{(8b)}& {\displaystyle }{I}_{\mathrm{v}}\left(T\right)& {\displaystyle }=c_{\mathrm{vv}}(T-T_{\mathrm{0}})+I_{\mathrm{v},\mathrm{0}},\text{(8c)}& {\displaystyle }{I}_{\mathrm{l}}\left(T\right)& {\displaystyle }=c_{\mathrm{vl}}(T-T_{\mathrm{0}}),\text{(8d)}& {\displaystyle }{I}_{\mathrm{i}}\left(T\right)& {\displaystyle }=c_{\mathrm{vi}}(T-T_{\mathrm{0}})-I_{\mathrm{i},\mathrm{0}}.\end{array}$$

Here, c_vk for k∈{d, v, l, i} represents isochoric specific heat capacities for the appropriate species denoted by k; they are taken to be constant. The reference specific internal energy I_v,0 is the difference in specific internal energy between vapor and liquid at the arbitrary reference temperature T₀; I_i,0 is the difference in specific internal energy between ice and liquid at T₀ (Romps, 2008). The reference internal energies are related to specific latent heats of vaporization and fusion, L_v,0 and L_f,0, at the reference temperature T₀ through

$$\begin{array}{}\text{(9)}& {\displaystyle }{I}_{\mathrm{v},\mathrm{0}}& {\displaystyle }=L_{\mathrm{v},\mathrm{0}}-R_{\mathrm{v}}{T}_{\mathrm{0}},\text{(10)}& {\displaystyle }{I}_{\mathrm{i},\mathrm{0}}& {\displaystyle }=L_{\mathrm{f},\mathrm{0}}.\end{array}$$

The values of the thermodynamic constants we use are listed in Table 1.

Table 1 Thermodynamic constants in CLIMAParameters.

Furthermore, the flux F_R is the radiative energy flux per unit mass; J is the conductive energy flux per unit mass, and D is the specific enthalpy flux associated with the diffusive flux of water:

$$\begin{array}{}\text{(11)}& \mathit{D}=(e_{\mathrm{v}}^{\mathrm{tot}}+R_{\mathrm{v}}T){\mathit{d}}_{q_{\mathrm{v}}}+e_{\mathrm{l}}^{\mathrm{tot}}{\mathit{d}}_{q_{\mathrm{l}}}+e_{\mathrm{i}}^{\mathrm{tot}}{\mathit{d}}_{q_{\mathrm{i}}}.\end{array}$$

The flux u⋅ρτ is the energy flux associated with the viscous and/or SGS turbulent momentum flux; Q is any internal energy source (e.g., external diabatic heating). The flux

$$\begin{array}{}\text{(12)}& W_{\mathrm{c}}=q_{\mathrm{l}}{e}_{\mathrm{l}}^{\mathrm{tot}}{w}_{\mathrm{l}}+q_{\mathrm{i}}{e}_{\mathrm{i}}^{\mathrm{tot}}{w}_{\mathrm{i}}\end{array}$$

represents the downward energy flux due to sedimenting condensate.

The terms involving ρC(q_j→q_p) (j∈{ v, l, i }) represent the loss of internal and potential energy of moist air masses owing to precipitation formation; the kinetic energy loss is neglected, consistent with the neglect of the source and/or sink associated with precipitation formation in the momentum balance in Eq. (4). Additional energy sinks involve the energy loss owing to heat transfer from the working fluid to precipitation as it falls through air and possibly melts at the freezing level (Raymond, 2013); the associated energy sources and sinks are generally provided by a microphysics parameterization and are subsumed in the term M.

2.6 Equation of state

Pressure p is calculated from the ideal gas law:

$$\begin{array}{}\text{(13)}& p=\mathit{\rho }{R}_{\mathrm{m}}T,\end{array}$$

where R_m is the gas “constant” of moist air,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{R}_{\mathrm{m}}\left(q\right)& =R_{\mathrm{d}}(\mathrm{1}-q_{\mathrm{t}})+R_{\mathrm{v}}{q}_{\mathrm{v}}\\ & =R_{\mathrm{d}}\left[\mathrm{1}+({\mathit{\epsilon }}_{\mathrm{dv}}-\mathrm{1})q_{\mathrm{t}}-{\mathit{\epsilon }}_{\mathrm{dv}}{q}_{\mathrm{c}}\right],\end{array}\end{array}$$

with the ratio of the gas constants of water vapor and of dry air ${\mathit{\epsilon }}_{\mathrm{dv}}=R_{\mathrm{v}}/R_{\mathrm{d}}$.

2.7 Saturation adjustment

Gibbs' phase rule states that in thermodynamic equilibrium, the temperature T and liquid and ice specific humidities q_l and q_i can be obtained from the three thermodynamic state variables density ρ, total water specific humidity q_t, and internal energy I. Thus, the above equations suffice to completely specify the thermodynamic state of the working fluid given ρ, q_t, and I, with the latter obtained from the total energy via its definition in Eq. (6).

Obtaining the temperature and condensate specific humidities from the state variables ρ, q_t, and I is the problem of finding the root T of

$$\begin{array}{}\text{(15)}& I^{*}(T;\mathit{\rho },q_{\mathrm{t}})-I=\mathrm{0},\end{array}$$

where $I^{*}(T;\mathit{\rho },q_{\mathrm{t}})$ is the internal energy at equilibrium, when the air is either unsaturated and there is no condensate (q_v=q_t), or water vapor is in saturation and the saturation excess q_t−q_v is apportioned, according to temperature T, among the condensed phases q_l and q_i. We solve this nonlinear “saturation adjustment” problem by Newton iterations with analytical gradients (see Tao et al., 1989; Pressel et al., 2015). To obtain the saturation vapor pressure and derived functions needed in this calculation, we assume all isochoric heat capacities to be constant (i.e., we assume the gases to be calorically perfect); with this assumption, the Clausius–Clapeyron can be integrated analytically, resulting in a closed-form expression for the saturation vapor pressure (Romps, 2008).

This procedure allows the use of total moisture q_t as the sole prognostic variable but confines the system to the assumption of equilibrium thermodynamics. Alternatively, using explicit tracers for the condensate specific humidities q_l and q_i allows non-equilibrium thermodynamics to be considered and mixed-phase processes to be explicitly modeled.

3 Discretization of the governing equations

3.1 Space discretization

The governing equations are discretized in space via a nodal DG approximation. To describe the DG procedure, we recast the Eqs. (1)–(5) in divergence form as

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\partial \mathit{Y}}{\partial t}}=-\mathrm{\nabla }\cdot ({\mathit{F}}_{\mathrm{1}}+{\mathit{F}}_{\mathrm{2}})+\mathcal{S}\left(\mathit{Y}\right),\end{array}$$

where $\mathit{Y}=[\mathit{\rho },\mathit{\rho }\mathit{u},\mathit{\rho }{e}^{\mathrm{tot}},\mathit{\rho }{q}_{\mathrm{t}}{]}^T$ is an abstract vector of state variables, F₁ contains the fluxes not involving gradients of state variables and functions thereof, F₂ contains the fluxes involving gradients of state variables (e.g., diffusive fluxes), and 𝓢(Y) contains the sources.

The DG solution of Eq. (16) is approximated on the finite-dimensional counterpart Ω^h of the flow domain Ω, which consists of $N_{{\mathrm{\Omega }}_{\mathrm{e}}}$ non-overlapping hexahedral elements Ω_e such that

$${\mathrm{\Omega }}^h=\bigcup _{e=\mathrm{1}}^{N_{{\mathrm{\Omega }}_{\mathrm{e}}}}{\mathrm{\Omega }}_{\mathrm{e}},$$

where a superscript h indicates the discrete analog of a continuous quantity. By virtue of tensor-product operations allowed on hexahedral elements and the ability to rely on inexact quadrature when elements of order greater than 3 are utilized, high-order Galerkin methods are particularly attractive for operation-intensive solutions (Kelly and Giraldo, 2012). Within each element, the finite-dimensional approximation of Y(x,t) is given by the expansion

$$\begin{array}{}\text{(17)}& {\mathit{Y}}_{\mathrm{e}}^h(\mathit{x},t)=\sum _{\mathit{\alpha }=\mathrm{1}}^{(N+\mathrm{1}{)}^{\mathrm{3}}}{\mathit{\psi }}_{\mathit{\alpha }}^{\mathrm{e}}\left(\mathit{x}\right){\mathit{Y}}_{\mathit{\alpha }}^{\mathrm{e}}\left(t\right),\end{array}$$

where (N+1)³ is the number of collocation points within the three-dimensional element of order N, and ${\mathit{\psi }}_{\mathit{\alpha }}^{\mathrm{e}}$ represents the interpolation polynomials evaluated at local point α inside element e.

From now on, the subscript or superscript e is omitted with the understanding that all operations are executed element-wise unless otherwise stated. Furthermore, the physical elements in the $\mathit{x}=(x,y,z)$ space are mapped to a reference element $\mathit{\xi }=(\mathit{\xi },\mathit{\eta },\mathit{\zeta })$. The three-dimensional basis functions ψ_α result from the one-dimensional functions L_i(ξ), L_j(η), and L_k(ζ) as the tensor product:

$${\mathit{\psi }}_{\mathit{\alpha }}\left(\mathit{\xi }\right)=L_i\left(\mathit{\xi }\right)\otimes L_{\mathrm{j}}\left(\mathit{\eta }\right)\otimes L_k\left(\mathit{\zeta }\right),\forall i,j,k=\mathrm{1},\mathrm{\dots },N+\mathrm{1}.$$

Each function L is a one-dimensional (1D) Lagrange polynomial defined on the 1D reference element $[-\mathrm{1},\mathrm{1}]$. The Lagrange function evaluated at points i along the ξ direction within the element is

$$L_i\left(\mathit{\xi }\right)=\prod _{l=\mathrm{1},l\ne i}^{N+\mathrm{1}}{\displaystyle \frac{\mathit{\xi }-{\mathit{\xi }}_{\mathrm{l}}}{{\mathit{\xi }}_i-{\mathit{\xi }}_{\mathrm{l}}}},$$

where ξ_i represents the N+1 co-located interpolation points along ξ. The polynomials L_j and L_k in the other two directions η and ζ are built in the same way. The N+1 interpolation points may be chosen in a variety of ways (Deville et al., 2002; Karniadakis and Sherwin, 1999); here we choose Legendre–Gauss–Lobatto (LGL) points (Giraldo and Restelli, 2008). The Kronecker δ property of the Lagrange polynomials is such that

$${\mathit{\psi }}_{\mathrm{i}}\left({\mathit{\xi }}_{\mathrm{l}}\right)={\mathit{\delta }}_{\mathrm{il}}$$

in 1D, which, in three dimensions (3D), translates to

$$\begin{array}{}\text{(18)}& {\mathrm{\Psi }}_{\mathit{\alpha }}({\mathit{\xi }}_a,{\mathit{\eta }}_b,{\mathit{\zeta }}_c)={\mathit{\delta }}_{ia}\otimes {\mathit{\delta }}_{jb}\otimes {\mathit{\delta }}_{kc}.\end{array}$$

This allows us to reduce the operation count as follows. We construct the space and time derivatives as

$$\begin{array}{}\text{(19)}& {\displaystyle \frac{\partial {\mathit{Y}}^h(\mathit{x},t)}{\partial \mathit{x}}}& {\displaystyle }=\sum _{\mathit{\alpha }=\mathrm{1}}^{(N+\mathrm{1}{)}^{\mathrm{3}}}{\displaystyle \frac{\partial {\mathit{\psi }}_{\mathit{\alpha }}\left(\mathit{x}\right)}{\partial \mathit{x}}}{\mathit{Y}}_{\mathit{\alpha }}\left(t\right),\text{(20)}& {\displaystyle \frac{\partial {\mathit{Y}}^h(\mathit{x},t)}{\partial t}}& {\displaystyle }=\sum _{\mathit{\alpha }=\mathrm{1}}^{(N+\mathrm{1}{)}^{\mathrm{3}}}{\mathit{\psi }}_{\mathit{\alpha }}\left(\mathit{x}\right){\displaystyle \frac{\partial {\mathit{Y}}_{\mathit{\alpha }}\left(t\right)}{\partial t}}.\end{array}$$

By virtue of the 3D Kronecker δ property, the spatial derivatives of the basis functions appearing here are given by

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial {\mathit{\psi }}_{\mathit{\alpha }}}{\partial \mathit{\xi }}}({\mathit{\xi }}_a,{\mathit{\eta }}_b,{\mathit{\zeta }}_c)={\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{i}}\left({\mathit{\xi }}_a\right)}{\partial \mathit{\xi }}}\otimes {\mathit{\delta }}_{jb}\otimes {\mathit{\delta }}_{ck},\text{(22)}& {\displaystyle \frac{\partial {\mathit{\psi }}_{\mathit{\alpha }}}{\partial \mathit{\eta }}}({\mathit{\xi }}_a,{\mathit{\eta }}_b,{\mathit{\zeta }}_c)={\mathit{\delta }}_{ia}\otimes {\displaystyle \frac{\partial {\mathit{\psi }}_{\mathrm{j}}\left({\mathit{\eta }}_b\right)}{\partial \mathit{\eta }}}\otimes {\mathit{\delta }}_{ck},\text{(23)}& {\displaystyle \frac{\partial {\mathit{\psi }}_{\mathit{\alpha }}}{\partial \mathit{\zeta }}}({\mathit{\xi }}_a,{\mathit{\eta }}_b,{\mathit{\zeta }}_c)={\mathit{\delta }}_{ia}\otimes {\mathit{\delta }}_{jb}\otimes {\displaystyle \frac{\partial {\mathit{\psi }}_k\left({\mathit{\zeta }}_c\right)}{\partial \mathit{\zeta }}}.\end{array}$$

Using this property reduces the operation count significantly since we only require 3N operations instead of the N³ operations otherwise needed to compute the derivatives at a given node (Abdi et al., 2017b).

The operators defined on the reference elements are mapped onto the physical space by means of the transformation

$$\begin{array}{}\text{(24)}& \mathrm{\nabla }\mathit{\psi }={\mathit{J}}^{-\mathrm{1}}\widehat{\mathrm{\nabla }}\mathit{\psi },\end{array}$$

where $\mathrm{\nabla }=({\partial }_x,{\partial }_y,{\partial }_z{)}^T$ and $\widehat{\mathrm{\nabla }}=({\partial }_{\mathit{\xi }},{\partial }_{\mathit{\eta }},{\partial }_{\mathit{\zeta }}{)}^T$, and

$${\mathit{J}}^{-\mathrm{1}}=\left[\begin{array}{ccc}{\mathit{\xi }}_x& {\mathit{\xi }}_y& {\mathit{\xi }}_z\\ {\mathit{\eta }}_x& {\mathit{\eta }}_y& {\mathit{\eta }}_z\\ {\mathit{\zeta }}_x& {\mathit{\zeta }}_y& {\mathit{\zeta }}_z\end{array}\right]$$

is the inverse Jacobian of the transformation from physical space to the reference element.

The DG approximation of the differential Eq. (16) is constructed by multiplying, within each element, the equation by the test function ψ_α and then integrating over the element volume Ω_e such that

$$\begin{array}{}\text{(25)}& \begin{array}{rl}\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathit{\psi }}_{\mathit{\alpha }}& \left({\partial }_{\mathrm{t}}\mathit{Y}+\mathrm{\nabla }\cdot {\mathit{F}}_{\mathrm{1}}\left(\mathit{Y}\right)+\mathrm{\nabla }\cdot {\mathit{F}}_{\mathrm{2}}(\mathit{Y},\mathrm{\nabla }\mathit{Y})\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\\ & =\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathit{\psi }}_{\mathit{\alpha }}\mathcal{S}\left(\mathit{Y}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}},\end{array}\end{array}$$

where ψ_α within each element belongs to the function space of square integrable piecewise polynomials of order N (i.e., ψ∈L²). By definition, these functions are discontinuous across element boundaries; differentiability is not globally required but only within each element (Hesthaven and Warburton, 2008b). Integrating the divergence term by parts yields

$$\begin{array}{}\text{(26)}& \begin{array}{rl}\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathrm{\Psi }}_{\mathit{\alpha }}{\partial }_{\mathrm{t}}\mathit{Y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}& +\underset{{\mathrm{\Gamma }}_{\mathrm{e}}}{\oint }{\mathrm{\Psi }}_{\mathit{\alpha }}\mathit{n}\cdot {\mathit{F}}_{\mathrm{1}}^{*}\left(\mathit{Y}\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}\\ & -\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathrm{\nabla }{\mathrm{\Psi }}_{\mathit{\alpha }}\cdot {\mathit{F}}_{\mathrm{1}}\left(\mathit{Y}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\\ & -\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathrm{\Psi }}_{\mathit{\alpha }}\mathrm{\nabla }\cdot {\mathit{F}}_{\mathrm{2}}(\mathit{Y},\mathbf{\nabla }\mathit{Y})\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}\\ & =\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathrm{\Psi }}_{\mathit{\alpha }}\mathcal{S}\left(\mathit{Y}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}},\end{array}\end{array}$$

where Ω_e and Γ_e are the volume and boundary of each element, respectively, n is the outward-facing unit vector orthogonal to each element face, and ${\mathit{F}}_{\mathrm{1}}^{*}$ is a numerical flux. The imposition of the numerical fluxes across element boundaries is the numerical mechanism that promotes continuity of the discontinuous solution across the elements. The numerical fluxes are calculated as the approximate solution to a Riemann problem across two neighboring elements. ClimateMachine currently implements the Rusanov (1961), Roe (1981), and Harten–Lax–van Leer contact (HLLC) (Toro et al., 1994; Harten, 1983) numerical fluxes. The Rusanov flux, for instance, is constructed as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\mathit{n}\cdot {\mathit{F}}_{\mathrm{1}}^{*}\left(\mathit{Y}\right)& ={\displaystyle \frac{\mathit{n}}{\mathrm{2}}}\cdot \left[{\mathit{F}}_{\mathrm{1}}\left({\mathit{Y}}^{-}\right)+{\mathit{F}}_{\mathrm{1}}\left({\mathit{Y}}^{+}\right)\right]\\ & +\mathit{n}{\mathit{\lambda }}_{{\mathrm{\Gamma }}_{\mathrm{e}}}\left({\mathit{Y}}^{-}-{\mathit{Y}}^{+}\right),\end{array}\end{array}$$

where Y⁻ is the state at the internal interface of element e, Y⁺ is the state at the external interface of e, and ${\mathit{\lambda }}_{{\mathrm{\Gamma }}_{\mathrm{e}}}({\mathit{Y}}^{-},{\mathit{Y}}^{+})$ is an estimate of the maximum flow speed (e.g., the maximum eigenvalue of the Jacobian of the flux F₁ with respect to the state variables, which is the speed of sound).

Because the second-order derivatives in ∇⋅F₂ cannot be directly built with the weak variational formulation if a discontinuous function space is used (Bassi and Rebay, 1997), an auxiliary variable $\stackrel{\mathrm{̃}}{\mathit{Y}}$ is introduced such that

$$\begin{array}{}\text{(28)}& {\displaystyle }\mathrm{\nabla }\mathit{Y}& {\displaystyle }=\stackrel{\mathrm{̃}}{\mathit{Y}}\text{(29)}& {\displaystyle }\mathrm{\nabla }\cdot \left(\mathit{\mu }\mathrm{\nabla }\mathit{Y}\right)& {\displaystyle }=\mathrm{\nabla }\cdot \left(\mathit{\mu }\stackrel{\mathrm{̃}}{\mathit{Y}}\right),\end{array}$$

which can then be discretized via DG as

$$\begin{array}{}\text{(30)}& \begin{array}{rl}\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathrm{\Psi }}_{\mathit{\alpha }}\mathrm{\nabla }\cdot \mathrm{\nabla }\mathit{Y}\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}& \approx \underset{{\mathrm{\Gamma }}_{\mathrm{e}}}{\oint }{\mathrm{\Psi }}_{\mathit{\alpha }}\mathit{n}\cdot \left(\mathit{\mu }{\stackrel{\mathrm{̃}}{\mathit{Y}}}^{*}-\mathit{\mu }\stackrel{\mathrm{̃}}{\mathit{Y}}\right)\mathrm{d}{\mathrm{\Gamma }}_{\mathrm{e}}\\ & +\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathrm{\Psi }}_{\mathit{\alpha }}\mathrm{\nabla }\cdot \left(\mathit{\mu }\stackrel{\mathrm{̃}}{\mathit{Y}}\right)\mathrm{d}{\mathrm{\Omega }}_{\mathrm{e}}.\end{array}\end{array}$$

Here, ${\stackrel{\mathrm{̃}}{\mathit{Y}}}^{*}$ is approximated via centered flux like in Bassi and Rebay (1997). We also refer to Abdi et al. (2017b) for more details.

For algorithmic efficiency, inexact quadrature is used to calculate the integrals above. By virtue of inexact integration and of Eqs. (17), (19), and (24), the variational DG equations yield the semi-discrete matrix problem:

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\mathrm{d}{\mathit{Y}}_{\mathrm{i}}^{\mathrm{e}}}{\mathrm{d}t}}=-\left({\mathrm{\nabla }}_{j,i}^T\right){\mathit{F}}_{\mathrm{j}}^{\mathrm{e}}+{\mathcal{S}}_{\mathrm{i}}^{\mathrm{e}}+{\displaystyle \frac{{w_{\mathrm{i}}}^{\mathrm{s}}|\mathit{J}{|}_{\mathrm{i}}^{\mathrm{s}}}{{w_{\mathrm{i}}}^{\mathrm{e}}|\mathit{J}{|}_{\mathrm{i}}^{\mathrm{e}}}}{\mathit{n}}_{\mathrm{i}}^{\mathrm{s}}{\left({\mathit{F}}^{\mathrm{e}}-{\mathit{F}}^{*}\right)}_{\mathrm{i}},\end{array}$$

where $w_{\mathrm{i}}^{\mathrm{s}}$ represents interpolation weights. The algebraic details to obtain this expression can be found in Giraldo and Restelli (2008), where the s superscript indicates a value that is defined on the element boundary surface. The system (Eq. 31) is integrated on each element with respect to time.

In order to achieve good parallel scaling it is necessary to overlap communication and computation to the fullest extent possible. With DG (and all element-based Galerkin methods) this can be naturally achieved by splitting Eq. (31) into terms that arise from the approximation of volume integrals and surface integrals. All volume contributions can be calculated independently of element-to-element communication regardless of the order of the spatial approximation, as can surface integrals that are not on elements which share boundaries across MPI ranks. Thus, in the code, we start with Message Passing Interface (MPI) communication, do all volume calculations and surface calculations for elements not on boundaries shared across ranks, and then apply surface calculations for elements on the rank boundaries after communication operations have been completed. This approach makes DG naturally effective with respect to parallel computing as previously shown by, e.g., Müller et al. (2018) on CPUs and Abdi et al. (2017b) on GPUs. At high order, element-based Galerkin methods such as DG require fewer neighboring degrees of freedom than high-order finite-difference and finite-volume discretizations.

3.2 Time discretization

ClimateMachine provides a suite of time integrators consisting of explicit Runge–Kunge methods as well as low-storage (Carpenter and Kennedy, 1994; Niegemann et al., 2012), strong stability-preserving (Shu and Osher, 1988), and additive Runge–Kutta (ARK) implicit–explicit (IMEX) methods (Giraldo et al., 2013; Kennedy and Carpenter, 2019).

The benchmarks presented in this paper with isotropic grid spacing are run using the fourth-order 14-stage method of Niegemann et al. (2012), which has a large explicit time-stepping stability region. One of the benchmarks, however, uses a highly anisotropic grid, which benefits from the use of a 1D IMEX approximation; there we use a variant of the horizontally explicit, vertically implicit (HEVI) schemes by Bao et al. (2015). While 3D IMEX is also an option, its performance in terms of time to solution is ultimately limited by the availability of scalable 3D implicit solver algorithms.

4 Subgrid-scale models

The governing equations are resolved with the discretizations presented in Sect. 3. This leaves unresolved but dynamically significant scales on the computational grid that must be modeled with three-dimensional SGS models. We note that, while filtering over the grid volume is typically a formal requirement for LES variables, no explicit filter kernel is applied to the equations here. The resolution of the prognostic variables on a discretized grid constitutes an implicit filtering operation, with an equivalent length scale given by the grid resolution. In general, SGS fluxes are modeled as diffusive fluxes, which capture down-gradient transport of conservable scalar quantities assuming that mixing lengths are small compared with the scales over which the gradients of the scalars vary. We address the physical form of the diffusive flux components in Eqs. (1)–(5), following which we describe standard models of subgrid-scale turbulence available for use in ClimateMachine.

The diffusive turbulent stress tensor τ is represented in terms of the symmetric rate of strain tensor S such that

$$\begin{array}{}\text{(32)}& \mathit{S}\left(\mathit{u}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{\nabla }\mathit{u}+{\left(\mathbf{\nabla }\mathit{u}\right)}^T\right),\end{array}$$

with

$$\begin{array}{}\text{(33)}& \mathit{\tau }=-\left(\mathrm{2}{\mathit{\nu }}_{\mathrm{t}}\mathit{S}\right).\end{array}$$

Here, ν_t is a turbulent viscosity tensor whose components are typically orders of magnitude larger than the molecular viscosity and are a function of the velocity gradient tensor.

The diffusive flux ${\mathit{d}}_{q_{\mathrm{t}}}$ of total water specific humidity in Eq. (2) is modeled as

$$\begin{array}{}\text{(34)}& {\mathit{d}}_{q_{\mathrm{t}}}=-\left({\mathcal{D}}_{\mathrm{t}}\mathbf{\nabla }{q}_{\mathrm{t}}\right),\end{array}$$

where 𝓓_t is a turbulent diffusivity vector. The turbulent diffusivity 𝓓_t is related to the turbulent viscosity tensor ν_t via the turbulent Prandtl number such that

$$\begin{array}{}\text{(35)}& {\mathcal{D}}_{\mathrm{t}}={\displaystyle \frac{\mathrm{diag}\left({\mathit{\nu }}_{\mathrm{t}}\right)}{{\mathrm{Pr}}_{\mathrm{t}}}},\end{array}$$

where Pr_t takes a typical value of ${\mathrm{Pr}}_{\mathrm{t}}=\mathrm{1}/\mathrm{3}$.

The unresolved flux of total enthalpy h^tot results in a diffusive subgrid flux term of the form

$$\begin{array}{}\text{(36)}& \mathit{J}+\mathit{D}=-\left({\mathcal{D}}_{\mathbf{t}}\mathbf{\nabla }{h}^{\mathrm{tot}}\right),\end{array}$$

where J is the thermal diffusion flux analogous to the molecular conductive heat flux, and D is the energy flux carried by water vapor, defined in Eq. (11). For energetic consistency, we use the same turbulent diffusivity 𝓓_t for moist enthalpy and water.

4.1 Smagorinsky–Lilly model

The turbulent eddy viscosity ν_t in the model by Lilly (1962) and Smagorinsky (1963) (SL henceforth) is defined by means of the magnitude of the rate of the strain tensor S, whose components are $S_{ij}=\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$, according to

$$\begin{array}{}\text{(37)}& {\mathit{\nu }}_{\mathrm{t}}=(C_{\mathrm{s}}\mathrm{\Delta }{)}^{\mathrm{2}}\sqrt{\mathrm{2}{S}_{ij}{S}_{ij}}\end{array}$$

for $i,j=\mathrm{1},\mathrm{2},\mathrm{3}$; C_s is a constant Smagorinsky coefficient usually within the range $\mathrm{0.12}<C_{\mathrm{s}}<\mathrm{0.21}$, and Δ is the LES filter width. Inside each hexahedral element of order N and side lengths $L_{x,y,z}$ along the x,y, and z directions, the effective grid resolution is $\mathrm{\Delta }(x,y,z)=L_{x,y,z}/N$, which is the average distance between two consecutive nodal points. We use an isotropic eddy viscosity tensor ν_t in LESs, with its components defined by Eq. (37). Similar eddy viscosity models can be found in other recent models of compressible atmospheric flows (e.g., Jähn et al., 2015; Shi et al., 2018).

4.2 Vreman eddy viscosity model

The SGS model developed by Vreman (2004) is of interest because of its robustness across flow regimes and because it has low dissipation near wall boundaries and in transitional flows. Its computational complexity is similar to the classical SL model. While the Vreman model is extensively used in engineering LESs, it is uncommon in atmospheric flows, for which a constant coefficient SL or the one-equation turbulent kinetic energy (TKE) model by Deardorff (1970, 1980) is the most common choice (see, e.g., Stevens et al., 2005).

The turbulent eddy viscosity of this model depends on first-order derivatives of velocities and is given by

$$\begin{array}{}\text{(38a)}& {\mathit{\nu }}_{\mathrm{t}}=\mathrm{2.5}{C}_{\mathrm{s}}^{\mathrm{2}}\sqrt{{\displaystyle \frac{B_{\mathit{\beta }}}{u_{i,j}{u}_{i,j}}}},\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }{B}_{\mathit{\beta }}& {\displaystyle }={\mathit{\beta }}_{\mathrm{11}}{\mathit{\beta }}_{\mathrm{22}}+{\mathit{\beta }}_{\mathrm{11}}{\mathit{\beta }}_{\mathrm{33}}+{\mathit{\beta }}_{\mathrm{22}}{\mathit{\beta }}_{\mathrm{33}}\\ \text{(38b)}& {\displaystyle }& {\displaystyle }-({\mathit{\beta }}_{\mathrm{13}}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{12}}^{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{23}}^{\mathrm{2}}),\text{(38c)}& {\displaystyle }{\mathit{\beta }}_{ij}& {\displaystyle }={\mathrm{\Delta }}_m^{\mathrm{2}}{u}_{i,m}{u}_{j,m},\text{(38d)}& {\displaystyle }{u}_{i,j}& {\displaystyle }={\displaystyle \frac{\partial u_i}{\partial x_j}}.\end{array}$$

Here, summation over repeated indices $i,j\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3}\mathit{\}}$ is implied, C_s is the constant Smagorinsky coefficient, and u_i represents the components of the resolved-scale velocity vector so that u_i,j is the velocity gradient tensor. The mixing lengths Δ_m can be determined as the grid spacing in the direction implied by subscript m; in this paper ${\mathrm{\Delta }}_m\equiv L_{x,y,z}/N$. Near-surface corrections to the characteristic length (e.g., Mason and Callen, 1986) have not been applied to either subgrid-scale model in the present work.

Richardson correction in stable regions of the atmosphere. To account for the atmospheric stability and, effectively, reduce turbulence generation to zero in stably stratified atmospheres, we multiply the eddy viscosity by the correction factor

$$\begin{array}{}\text{(39)}& f_b=\left\{\begin{array}{ll}\mathrm{1}& Ri\le \mathrm{0},\\ max(\mathrm{0},\mathrm{1}-Ri/P{r}_{\mathrm{t}}{)}^{\mathrm{1}/\mathrm{4}}& Ri>\mathrm{0},\end{array}\right.\end{array}$$

where

$$Ri={\displaystyle \frac{\left(\frac{g}{{\mathit{\theta }}_{\mathrm{v}}}\frac{\partial {\mathit{\theta }}_{\mathrm{v}}}{\partial z}\right)}{|\mathit{S}{|}^{\mathrm{2}}}}.$$

Here, $P{r}_{\mathrm{t}}=\mathrm{1}/\mathrm{3}$ is a constant turbulent Prandtl number, and θ_v is the virtual potential temperature for a given specific humidity q_t and specific liquid water content q_l (see, e.g., Deardorff, 1980). This temperature is computed consistently with the moist phase partitioning (saturation adjustment) in the system and is therefore a reasonable measure of the effect of buoyancy on flow stratification. This is applied to both turbulence closures discussed in this paper.

4.3 Numerical stability

When high-order Galerkin methods are used to solve nonlinear advection-dominated problems, spurious Gibbs oscillations may affect the solution and need to be addressed. ClimateMachine provides a set of spectral filters, cut-off filters, and artificial diffusion methods to remove these oscillations. While filters may be effective, we found that stabilizing the LES solution by means of the SGS eddy viscosity alone is effective and robust; this is in agreement with results shown by Marras et al. (2015), Marras and Giraldo (2015), and Reddy et al. (2022) in the case of continuous Galerkin methods. Other methods of stabilizing numerical solutions are discussed by Light and Durran (2016) and Yu et al. (2015). This approach stems from the idea that under-resolved scales are responsible for the numerical oscillations observed in solutions through the introduction of unphysical artifacts in properties of the solution variables. Detailed analyses of the interactions of subgrid-scale models and filtering techniques with DG numerics in the context of atmospheric flows will be presented in a forthcoming paper.

5 Numerical experiments and discussion

We first demonstrate the convergence of the numerical solution to the Euler equations with an isentropic vortex advection problem. Following this, we demonstrate results from the ClimateMachine using standard benchmark problems including (1) dry rising thermal bubble in a neutrally stratified atmosphere, (2) dry density current, (3) hydrostatic and nonhydrostatic mountain-triggered linear gravity waves, (4) the Barbados Oceanographic and Meteorological Experiment (BOMEX), and (5) decaying Taylor–Green vortex in a triply periodic domain.

5.1 Advection of an isentropic vortex

To demonstrate convergence of the numerical solutions, we consider the two-dimensional dry isentropic vortex advection problem with fourth-order polynomials, consistent with the benchmark cases shown in the sections that follow. We consider the pure advection of a vortex in a domain with edge lengths L=0.1 m, with a vortex radius R=0.005 m. The initial velocity profile is given by

$$\begin{array}{}\text{(40)}& \mathit{u}=(u_{\mathrm{0}}\cos \mathit{\alpha }+{\mathit{\delta }}_u,u_{\mathrm{0}}\sin \mathit{\alpha }+{\mathit{\delta }}_{\mathrm{v}}),\end{array}$$

with a prescribed translational velocity u₀=150 m s⁻¹, translation angle $\mathit{\alpha }=\mathit{\pi }/\mathrm{4}$, and a vortex speed v_s=50 m s⁻¹ so that the velocity perturbations are given by

$$\begin{array}{}\text{(41)}& {\mathit{\delta }}_u=-v_{\mathrm{s}}x{\displaystyle \frac{y^{\prime }}{R}}\exp ({\displaystyle \frac{r}{\mathrm{2}R}}{)}^{\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(42)}& {\mathit{\delta }}_{\mathrm{v}}=v_{\mathrm{s}}x{\displaystyle \frac{x^{\prime }}{R}}\exp ({\displaystyle \frac{r}{\mathrm{2}R}}{)}^{\mathrm{2}}.\end{array}$$

The perturbations are computed with an offset coordinate given by

$$\begin{array}{}\text{(43)}& {\displaystyle }{x}^{\prime }=x-\mathrm{2}L\mathrm{floor}\left({\displaystyle \frac{x+L}{\mathrm{2}L}}\right),\text{(44)}& {\displaystyle }{y}^{\prime }=y-\mathrm{2}L\mathrm{floor}\left({\displaystyle \frac{y+L}{\mathrm{2}L}}\right),\end{array}$$

where x and y represent Cartesian coordinates, to ensure that the resulting functions are periodic. The initial temperature profile is given by

$$\begin{array}{}\text{(45)}& T=T_{\mathrm{\infty }}(\mathrm{1}-{\mathit{\kappa }}_{\mathrm{d}}{\displaystyle \frac{v_{\mathrm{s}}^{\mathrm{2}}}{\mathrm{2}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{\infty }}}{p_{\mathrm{\infty }}}}\exp (-{\displaystyle \frac{r}{R}}{)}^{\mathrm{2}}),\end{array}$$

with a reference temperature of T_∞=300 K, reference pressure p_∞=10⁵ Pa, and the corresponding density ρ_∞ computed from the ideal gas law. The simulation time is given by $t_{\mathrm{end}}=L/\left(\mathrm{10}{u}_{\mathrm{0}}\right)$. We consider fourth-order polynomials, with increasing refinement given by the addition of elements in the domain discretization, to generate the Euclidean distance between the initial state and the time-stepped solution state at $t_{\mathrm{end}}={\mathrm{10}}^{-\mathrm{4}}$ s, as shown in Fig. 1. The results demonstrate consistent convergence across three numerical flux methods for element numbers ranging from 5–320. The convergence rate approaches 5 for the fourth-order polynomials considered in this study.

Figure 1 Euclidean distance between initial and final solution vectors of the prognostic variables generated using Rusanov (octagons), Roe (squares), and Harten–Lax–van Leer contact (HLLC) numerical fluxes. Final solutions are evaluated at $t_{\mathrm{end}}={\mathrm{10}}^{-\mathrm{4}}$ s based on a translational velocity of 150 m s⁻¹. For each numerical flux method, we consider $N_{\mathrm{e}}=[\mathrm{5},\mathrm{10},\mathrm{20},\mathrm{40},\mathrm{80},\mathrm{160},\mathrm{320}]$ elements with fourth-order polynomials. The dashed line represents $N_{\mathrm{e}}^{-\mathrm{5}}$ convergence.

Tests presented in Sect. 5.2–5.5 are executed in a 3D domain even when the problem is two-dimensional, in which case we have effectively zero tendencies in the third dimension; this setup is identified as 2.5D in what follows.

5.2 2.5D rising thermal bubble in a neutrally stratified atmosphere

A neutrally stratified atmosphere with uniform background potential temperature θ₀=300 K is perturbed by a circular bubble of warmer air. The hydrostatic background pressure decreases with z as

$$\begin{array}{}\text{(46)}& p=p_{\mathrm{0}}{\left(\mathrm{1}-{\displaystyle \frac{g}{c_{\mathrm{pd}}{\mathit{\theta }}_{\mathrm{0}}}}z\right)}^{c_{\mathrm{pd}}/R_{\mathrm{d}}}\end{array}$$

in a domain $\mathrm{\Omega }=[\mathrm{0},\mathrm{10}\mathrm{000}]\times [-\mathrm{\infty },\mathrm{\infty }]\times [\mathrm{0},\mathrm{10}\mathrm{000}]$ m³. The perturbation is as defined in Ahmad and Lindeman (2007),

$$\begin{array}{}\text{(47)}& \mathrm{\Delta }\mathit{\theta }={\mathit{\theta }}_{\mathrm{c}}\left[\mathrm{1.0}-{\displaystyle \frac{r}{r_{\mathrm{0}}}}\right]\text{if}r\le r_{\mathrm{0}}=\mathrm{2000}\mathrm{m},\end{array}$$

where $r=\sqrt{(x-x_{\mathrm{c}}{)}^{\mathrm{2}}+(z-z_{\mathrm{c}}{)}^{\mathrm{2}}}$, $(x_{\mathrm{c}},z_{\mathrm{c}})=(\mathrm{5000},\mathrm{2000})$ m, and θ_c=2 K. The initial velocity field is zero everywhere. Periodic boundary conditions are used along y, and solid walls with impenetrable, free-slip boundary conditions are used in the x and z directions. Detailed information on boundary conditions for all test cases is provided in Appendix A. Five runs are performed at effective uniform resolutions $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{250}$, 125, 62.5, and 31.25, and 15.625 m, with polynomial order N=4. Potential temperature θ and the two velocity components u and w are plotted at t=1000 s in Fig. 2 for the grid resolution of 15.625 m, which represents a reference solution for comparison with solutions at coarser resolutions. The value of the maximum potential temperature perturbation Δθ_max and of the horizontal and vertical velocity components agree with the 125 m resolution results shown by Ahmad and Lindeman (2007). The grid dependence of the solution is shown in Fig. 3, where potential temperature is plotted for $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{31.25}$, 62.5, 125, and 250 m. While the solution is visibly more dissipative at coarser resolutions, the bubble's leading edge position (and hence propagation speed) is not sensitive to the grid resolution.

The SL and Vreman closures are used to model diffusive fluxes in this problem. The solutions show no discernible differences, and only the SL solution is shown. A visual comparison of the two becomes more meaningful when shear triggers mixing, which is shown for the density current test in Sect. 5.3. Although DG inherits an implicit numerical diffusion as an effect of the numerical flux calculation across elements, under-resolved advection-dominated problems still require a dissipation or filtering mechanism to preserve the solution's stability. In the case of Marras et al. (2015), a dynamically adaptive SGS model was used, whereas a Boyd–Vandeven filter (Boyd, 1996; Vandeven, 1991) was used by Giraldo and Restelli (2008).

Figure 2 2.5D rising thermal bubble with effective resolution $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{15.625}$ m and N=4. Panels depict (a) potential temperature θ, (b) horizontal velocity u, and (c) vertical velocity w at t=1000 s in $\mathrm{\Omega }=\mathrm{10}\times \mathrm{10}$ km².

Figure 3 2.5D rising thermal bubble solution with N=4 at decreasing effective resolution. Grid convergence of potential temperature θ at four different resolutions to be compared against the 15.625 m resolution results shown in Fig. 2. (a) $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{31.25}$ m, (b) 62.5 m, (c) 125 m, and (d) 250 m. Although the solution is visibly more dissipated at coarser resolution, the bubble's leading edge (and hence propagation speed) is not affected.

5.3 2.5D density current in a neutrally stratified atmosphere

The density current problem by Straka et al. (1993) is used to test the LES framework in a flow with Kelvin–Helmholtz instabilities. As for the rising thermal bubble, the background initial state is in hydrostatic equilibrium at uniform potential temperature θ₀=300 K. A perturbation of θ centered on $(x_{\mathrm{c}},z_{\mathrm{c}})=(\mathrm{0},\mathrm{3000})$ m and with radii $(r_x,r_z)=(\mathrm{4000},\mathrm{2000})$ m is given by the function

$$\begin{array}{}\text{(48)}& \mathrm{\Delta }\mathit{\theta }={\displaystyle \frac{{\mathit{\theta }}_{\mathrm{c}}}{\mathrm{2}}}\left[\mathrm{1}+\cos \left({\mathit{\pi }}_{\mathrm{c}}r\right)\right]\text{if}r\le \mathrm{1},\end{array}$$

where ${\mathit{\theta }}_{\mathrm{c}}=-\mathrm{15}$ K and $r=\sqrt{(x-x_{\mathrm{c}})/r_x^{\mathrm{2}}+(z-z_{\mathrm{c}})/r_z^{\mathrm{2}}}$ in the domain $\mathrm{\Omega }=[\mathrm{0},\mathrm{25}\mathrm{600}]\times [-\mathrm{\infty },\mathrm{\infty }]\times [\mathrm{0},\mathrm{6400}]$ m³. Periodic boundary conditions are used along y; impenetrable free-slip conditions are imposed in x and z. The flow is initially stationary.

To reach solution grid convergence, this test is classically executed with a constant kinematic viscosity ν=75 m² s⁻¹. Increasingly finer structures are resolved when the resolution increases (Marras et al., 2012, 2015). A measure of solution fidelity is the front position, which we compare against other models in Appendix D for different resolutions. The ClimateMachine results show quantitative agreement with respect to the frontal location from a range of models with varying spatial discretizations at resolutions ranging from 12.5 to 100 m. This demonstrates the scope for capturing small-scale flow features with the numerics described in Sect. 3.

The structure of potential temperature at the final time t=900 s is shown in Fig. 4 for the Vreman and SL solutions. When Rusanov is the chosen numerical flux (Fig. 4a and b), the solutions are very similar although Vreman is visibly less dissipative for a prescribed value of the Smagorinsky coefficient C_s=0.18, with finer scales of motion apparent in the contours of potential temperature. Further quantitative analysis of the dissipative properties of the SGS models is reported in Sect. 5.6. Despite Vreman being less dissipative than SL, the Roe (Fig. 4c) and HLLC (Fig. 4d) fluxes contribute to additional numerical diffusion when compared to Rusanov fluxes. Detailed analysis of the interaction between numerical fluxes and subgrid-scale models will be presented in future articles.

Figure 4 2.5D density current. Potential temperature θ (K) at t=900 s computed with an effective DG resolution of $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{12.5}$ m with domain extents shown in meters. (a) Rusanov numerical flux with SL SGS. (b) Rusanov numerical flux with Vreman SGS. (c) Roe numerical flux with Vreman SGS. (d) HLLC numerical flux with Vreman SGS. The color scale ranges from θ=285 to 300 K. A shared color bar for all plots is shown in panel (a) for clarity.

5.4 Passive transport over warped grids

To verify the correct behavior of the DG implementation in the presence of topographic features, the simple passive advection test described by Schär et al. (2002) is used. The conservation law for the diffusive transport of a passive tracer χ is

$$\begin{array}{}\text{(49)}& {\displaystyle \frac{\partial \left(\mathit{\rho }\mathit{\chi }\right)}{\partial t}}+\mathrm{\nabla }\cdot \left(\mathit{\rho }\mathit{u}\mathit{\chi }\right)=-\mathrm{\nabla }\cdot \left(\mathit{\rho }{\mathit{d}}_{\mathit{\chi }}\right),\end{array}$$

which is approximated via DG in the same way as Eq. (16). For a scalar tracer variable χ, we model diffusive fluxes d_χ such that

$$\begin{array}{}\text{(50)}& {\mathit{d}}_{\mathit{\chi }}=-\left({\mathit{\delta }}_{\mathit{\chi }}{\mathcal{D}}_{\mathrm{t}}\mathbf{\nabla }\mathit{\chi }\right),\end{array}$$

where δ_χ relates the ratio of turbulent diffusivity of the tracer to that of the energy and moisture variables. For this test case, however, tracer diffusivity is set to zero to assess the stability and transport properties when using a warped grid. The volume grid in ClimateMachine is built by stacking elements above the surface and warping them around the terrain profile. To reduce the element distortion across the domain, a linear grid damping function is used such that a topography-conforming surface of nodal points near the domain's bottom surface decays to a horizontal plane at higher altitudes (Gal-Chen and Somerville, 1975).

The initial scalar field χ is described by an elliptical perturbation centered on $(x_{\mathrm{c}},z_{\mathrm{c}})=(\mathrm{25},\mathrm{9})$ km with radii $(r_x,r_z)=(\mathrm{25},\mathrm{3})$ km such that

$$\begin{array}{}\text{(51)}& \mathit{\chi }=\left\{\begin{array}{ll}{\mathit{\chi }}_{\mathrm{0}}{\cos }^{\mathrm{2}}\left(\frac{\mathit{\pi }r}{\mathrm{2}}\right),& \text{for }r\le \mathrm{1},\\ \mathrm{0}\text{otherwise },\end{array}\right.\end{array}$$

where χ₀=1 and $r=\sqrt{(x-x_{\mathrm{c}})/r_x^{\mathrm{2}}+(z-z_{\mathrm{c}})/r_z^{\mathrm{2}}}$ in the domain $\mathrm{\Omega }=[\mathrm{0},\mathrm{150}\mathrm{000}]\times [-\mathrm{\infty },\mathrm{\infty }]\times [\mathrm{0},\mathrm{30}\mathrm{000}]$ m³. An effective uniform grid resolution $\mathrm{\Delta }x=\mathrm{\Delta }z=\mathrm{500}$ m is used for this test. The initial velocity profile is given by

$$\begin{array}{}\text{(52)}& u\left(z\right)=u_{\mathrm{0}}\left\{\begin{array}{ll}\mathrm{1},& \text{for }z\ge z_{\mathrm{2}}\\ \sin \left(\frac{\mathit{\pi }}{\mathrm{2}}\frac{z-z_{\mathrm{1}}}{z_{\mathrm{2}}-z_{\mathrm{1}}}\right),& \text{for }{z}_{\mathrm{1}}\le z<z_{\mathrm{2}},\\ \mathrm{0},& \text{for }z<z_{\mathrm{1}},\end{array}\right.\end{array}$$

where u₀=10 m s⁻¹, z₁=4 km, and z₂=5 km.

The topography is defined by the function

$$\begin{array}{}\text{(53)}& \begin{array}{rl}& z_{\mathrm{sfc}}\left(x\right)=\\ & \left\{\begin{array}{ll}{h}_{\mathrm{0}}{\cos }^{\mathrm{2}}\left(\frac{\mathit{\pi }(x-x_{\mathrm{0}})}{\mathrm{2}a}\right){\cos }^{\mathrm{2}}\left(\frac{\mathit{\pi }(x-x_{\mathrm{0}})}{\mathit{\lambda }}\right)& \text{for }|x-x_{\mathrm{0}}|\le a\\ \mathrm{0},& \text{for }|x-x_{\mathrm{0}}|>a,\end{array}\right.\end{array}\end{array}$$

where h₀=3 km, a=25 km, λ=8 km, and x₀=75 km.

The contours of χ in Fig. 5 show minimal distortion in spite of the warped elements directly above the topographical feature, indicating that the DG transformation metrics from physical to logical space in the presence of topography do not adversely affect the solution. Deviation from the initial profile of tracer magnitudes is between −4 % and +2 % when the tracer is above the topographical feature, with maximum deviation amplitudes of −5 % and +3 % at the end of the test, showing a favorable comparison with the hybrid and SLEVE coordinate results presented in Schär et al. (2002).

Figure 5 Solution of the passive transport of scalar χ at three different time instances, with contours of tracer quantity χ (from 0, outermost contour, to 1, innermost contour) overlaid on a representation of nodal points on the underlying mesh. Tracer advection is driven by a prescribed velocity profile from left to right. As it crosses the deformed grid above the mountain ridge, only a minimal distortion of tracer contours is observed, which is completely recovered back to a smooth solution downwind of the ridge.

5.5 Mountain-triggered gravity waves

To assess the correct implementation of a Rayleigh sponge layer to attenuate fast, upward-propagating gravity waves before they reach the top of the domain, two steady-state mountain-triggered gravity wave problems suggested by Smith (1980) are solved. The sponge layer is described in Appendix A2. These problems consist of a flow that moves eastward with uniform horizontal velocity $\mathit{u}=(u,\mathrm{0},\mathrm{0})$ m s⁻¹ in a doubly periodic domain. The flow impinges against a mountain of height h_m and base length a centered at x_c as

$$\begin{array}{}\text{(54)}& z_{\mathrm{sfc}}\left(x\right)={\displaystyle \frac{h_{\mathrm{m}}{a}^{\mathrm{2}}}{(x-x_{\mathrm{c}}{)}^{\mathrm{2}}+a^{\mathrm{2}}}}.\end{array}$$

The background state is in hydrostatic balance with Brunt–Väisälä frequency 𝒩 such that

$$\mathit{\theta }={\mathit{\theta }}_{\mathrm{sfc}}\exp \left({\displaystyle \frac{{\mathcal{N}}^{\mathrm{2}}}{g}}z\right)$$

for a given surface potential temperature θ_sfc=T_sfc. The hydrostatically balanced pressure is

$$\begin{array}{}\text{(55)}& p=p_{\mathrm{sfc}}{\left[\mathrm{1}+{\displaystyle \frac{g^{\mathrm{2}}}{c_{\mathrm{pd}}{\mathit{\theta }}_{\mathrm{sfc}}{\mathcal{N}}^{\mathrm{2}}}}\left(\exp \left({\displaystyle \frac{-z{\mathcal{N}}^{\mathrm{2}}}{g}}\right)-\mathrm{1}\right)\right]}^{c_{\mathrm{pd}}/R_{\mathrm{d}}},\end{array}$$

which yields, by means of the ideal gas law, the background density:

$$\begin{array}{}\text{(56)}& \mathit{\rho }={\displaystyle \frac{p_{\mathrm{sfc}}}{\mathit{\theta }{R}_{\mathrm{d}}{\left(\frac{p}{p_{\mathrm{sfc}}}\right)}^{c_{\mathrm{vd}}/c_{\mathrm{pd}}}}}.\end{array}$$

These tests are affected by spurious oscillations that appear approximately 5000 s into the simulation. In the absence of shear, because of the free-slip bottom and top boundaries, the SGS models are unable to introduce sufficient diffusion to remove the Gibbs modes so that an exponential filter (Hesthaven and Warburton, 2008b) of order 64 was applied on the velocity field to remove spurious modes. The filter assumes the form of

$$\begin{array}{}\text{(57)}& \mathit{\sigma }\left(\mathit{\eta }\right)=e^{-\mathit{\alpha }{\mathit{\eta }}^s},\end{array}$$

where s is the filter order, η is a function of the polynomial order, and $\mathit{\alpha }=-\log \left({\mathit{\epsilon }}_{\mathrm{M}}\right)$ is a parameter that controls the smallest value of the filter function for machine precision ε_M. In double precision, α≈36. The filter in this form is applied to perturbations of the prognostic variables from the balanced background state.

5.5.1 Linear hydrostatic

The linear hydrostatic case proposed by Smith (1979) consists of a neutrally stratified isothermal atmosphere with $\mathit{\theta }={\mathit{\theta }}_{\mathrm{sfc}}=\mathrm{250}$ K. The background atmosphere is isothermal with temperature T₀, resulting in a Brunt–Väisälä frequency of

$$\mathcal{N}={\displaystyle \frac{g}{\sqrt{c_{\mathrm{pd}}{T}_{\mathrm{0}}}}}.$$

The flow moves in a periodic channel along the x direction with velocity $\mathit{u}=(\mathrm{20}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},\mathrm{0},\mathrm{0})$ over a mountain with h_m=1 m and a=10 000 m. A Rayleigh absorbing layer is added at z_s=25 km with relaxation coefficient α=0.5 s⁻¹, power γ=2, and domain top z_top=30 km (see Appendix A2 for details). The domain extends from 0 to 240 km in the horizontal direction. The steady-state solution at t=15 000 s is shown in Fig. 6a. It is consistent with the DG results shown by Giraldo and Restelli (2008).

5.5.2 Linear nonhydrostatic

The linear nonhydrostatic mountain waves are forced by a flow of uniform horizontal velocity $\mathit{u}=(\mathrm{10}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}},\mathrm{0},\mathrm{0})$ over a mountain with h_m=1 m and a=1000 m. The domain extends from 0 to 144 km in the horizontal direction and is 30 km high.

The steady-state solution at t=18 000 s is shown in Fig. 6b. It is consistent with results shown by Giraldo and Restelli (2008).

Figure 6 Vertical velocity w for two linear mountain wave tests. Velocity contours are shown in the range from $-\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ (blue) to $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ (red). (a) Hydrostatic test at t=15 000 s, with an effective horizontal resolution Δx=500 m and effective vertical resolution Δz=240 m with N=4. (b) Nonhydrostatic test at t=18 000 s, with an effective horizontal resolution Δx=200 m and effective vertical resolution Δz=160 m with N=4.

5.6 Decaying Taylor–Green vortex

The decaying Taylor–Green vortex (TGV) is a classical test to estimate the dissipative properties of turbulence models in the absence of solid boundaries. The gravity-free flow is initialized in a triply periodic cube of dimensions $[-\mathit{\pi },\mathit{\pi }{]}^{\mathrm{3}}$. The solenoidal initial velocity field ${\mathit{u}}_{\mathbf{0}}=(u_{\mathrm{0}},v_{\mathrm{0}},w_{\mathrm{0}})$ is defined as

$$\begin{array}{}\text{(58)}& {\displaystyle }{u}_{\mathrm{0}}& {\displaystyle }=U_{\mathrm{0}}\sin \left(kx\right)\cos \left(ky\right)\cos \left(kz\right),\text{(59)}& {\displaystyle }{v}_{\mathrm{0}}& {\displaystyle }=U_{\mathrm{0}}\cos \left(kx\right)\sin \left(ky\right)\cos \left(kz\right),\text{(60)}& {\displaystyle }{w}_{\mathrm{0}}& {\displaystyle }=\mathrm{0},\end{array}$$

with initial pressure

$$\begin{array}{}\text{(61)}& p_{\mathrm{0}}=p_{\mathrm{\infty }}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}{U}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{16}}}(\mathrm{2}+\cos (kz\left)\right)\left(\cos \right(\mathrm{2}kx)+\cos (\mathrm{2}ky\left)\right),\end{array}$$

where k is the wavenumber, U₀=100 m s⁻¹, ρ₀=1.178 kg m⁻³, and p_∞=101 325 Pa. Fourth-order polynomials are used for all simulations considered in this section.

We first consider the volume-averaged kinetic energy, which provides insight into the dissipation characteristics of the flow with respect to nondimensionalized time $t^{*}=k{U}_{\mathrm{0}}t$. In integral form, the kinetic energy can be written as

$$\begin{array}{}\text{(62)}& E_k={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\langle \mathrm{|}\mathit{u}{\mathrm{|}}^{\mathrm{2}}\rangle ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathrm{\Omega }}^h}}\underset{{\mathrm{\Omega }}^h}{\int }\mathit{u}\cdot \mathit{u}\mathrm{d}{\mathrm{\Omega }}^h,\end{array}$$

where 〈⋅〉 denotes a volumetric average over the volume Ω^h. If the flow is inviscid, the kinetic energy should be conserved. This is only valid if the numerics or SGS models do not introduce numerical dissipation or if all flow scales are well resolved. As such, the time series of kinetic energy is a metric that shows the point along the simulation at which the solution becomes under-resolved. The kinetic energy dissipation rate is the second quantity of interest, which allows us to quantify the rate of decay of kinetic energy over time. This is defined as

$$\begin{array}{}\text{(63)}& \mathit{ϵ}=-{\displaystyle \frac{\mathrm{d}{E}_k}{\mathrm{d}t}}.\end{array}$$

A third quantity of interest for this analysis is enstrophy, which is defined as the square of the vorticity norm:

$$\begin{array}{}\text{(64)}& \langle {\mathit{\omega }}^{\mathrm{2}}\rangle =\langle \mathrm{|}\mathrm{|}\mathbf{\nabla }\times \mathit{u}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}\rangle .\end{array}$$

The enstrophy of a fully resolved flow should go to infinity if the flow is inviscid. Therefore, enstrophy can be used as a criterion to estimate the effect of numerical dissipation.

By means of a three-dimensional fast Fourier transform (FFT) of the velocity field, the kinetic energy spectrum is calculated as

$$\begin{array}{}\text{(65)}& E\left(k\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{K}{\int }}A(k_x,\mathit{\upsilon },\mathit{\zeta })k^{\mathrm{2}}\sin \left(\mathit{\zeta }\right)\mathrm{d}k\mathrm{d}\mathit{\upsilon }\mathrm{d}\mathit{\zeta },\end{array}$$

where $K=\mathrm{2}\mathit{\pi }/L$, L is the characteristic length, A is a three-dimensional array of Fourier mode amplitudes, $\mathit{\upsilon }=k_z/k$, $\mathit{\zeta }={\tan }^{-\mathrm{1}}(k_y/k_x)$, and $k=\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}+k_z^{\mathrm{2}}}$. The TGV flow is simulated using both the SL and Vreman models on grids with 32³, 64³, 128³, and 192³ points. Figure 7 shows a 3D visualization of the flow at two different nondimensional times using zero Q-criterion isosurfaces, which identify balance between rotation and shear in the flow (Hunt et al., 1988). As the flow evolves, the flow generates smaller and smaller-scale vortices. Eventually, the flow becomes under-resolved, making it impossible to conserve kinetic energy. As the flow continues to evolve, an instability occurs, which causes the disintegration of the vortex sheet. After this point, the TGV's dynamics are controlled by the interaction of small-scale vortical structures formed by vortex stretching.

Figure 7 Taylor–Green vortex. Isosurfaces of zero Q-criterion (these identify surfaces where the vorticity norm is identical to the strain rate magnitude) on a 192³ grid, corresponding to 48 elements of order 4 at (a) $t^{*}=\mathrm{4}$ and (b) $t^{*}=\mathrm{60}$. Plots are shown for the Vreman solution. The isosurfaces of the Q-criterion are colored by dimensionless kinetic energy.

Results for the coarse-resolution simulations are presented in Fig. 8, and those for the fine-resolution simulations are shown in Fig. 9. Figure 8a shows that the kinetic energy changes with time for the 32³ resolution simulations are distinctly different from their higher-resolution counterparts in Fig. 9a. The severe under-resolution of the flow seems to generate much larger amounts of dissipation early on. This is further demonstrated when comparing Figs. 8b and 9b. Beyond 64³ resolution, the peaks in dissipation are larger as the resolution increases, since smaller vortices can be resolved before the instability eventually happens. The 32³ simulations (Fig. 8b) have larger peaks than even the 192³ (Fig. 9b) simulations, demonstrating that the under-resolution of the vortex structures leads to different, more dissipative early-flow behavior.

Figure 9b shows that the 128³ and 192³ simulations using both SGS closures are characterized by a peak in the kinetic energy dissipation at $t^{*}=\mathrm{9}$. Brachet et al. (1983) and Brachet (1991) demonstrate this result for their direct numerical simulations (DNSs) of the Taylor–Green vortex for a Reynolds number $R_{\mathrm{e}}=U_{\mathrm{0}}/k\mathit{\nu }\ge \mathrm{3000}$. Examining Fig. 8b, the 64³ simulations suggest a dissipation peak time of $t^{*}=\mathrm{8}$, and the 32³ simulations suggest a dissipation peak time of $t^{*}=\mathrm{6}$. The underprediction of the time at which the peak occurs is due to an inability to resolve the vortices that appear early on in the flow's evolution at extremely coarse resolutions. This leads to the early appearance of the instability which causes the dissipation peak. Furthermore, we see that the kinetic energy decay occurs sooner for the lower-resolution simulations as a result of increased dissipation from the SGS models.

Figure 8 Evolution of volumetrically averaged (a) kinetic energy represented on a semi-logarithmic scale, (b) kinetic energy dissipation, and (c) enstrophy, each computed on 32³ and 64³ grids with the Vreman and SL SGS models.

Figure 9 Evolution of volumetrically averaged quantities in the numerical solution to the Taylor–Green vortex problem computed on 128³ and 192³ grids with the Vreman and SL SGS models.

Figure 9c also shows that the flow's enstrophy behaves as expected, with peak values at $t^{*}=\mathrm{9}$ coinciding with peaks in kinetic energy. The higher-resolution simulations are able to reach a higher enstrophy than the lower-resolution simulations as they are naturally able to resolve more vortical motion. On the other hand, the choice of SL or Vreman models seems to have very little impact on the ability to resolve more small-scale eddies for low resolutions. However, the SL SGS scheme leads to higher enstrophy than the Vreman SGS scheme, increasingly so as resolution increases.

Figure 10 shows the kinetic energy spectra obtained for the higher-resolution simulations used for this test. All simulations present peaks in their respective spectra at k=2 to 4, which persist even as the flow evolves over time. These peaks are explained by Drikakis et al. (2007) as being imprinted on the spectra by the initial velocity field. Furthermore, as the flow evolves, all spectra show a slope close to the theoretical $k^{-\mathrm{5}/\mathrm{3}}$ of homogeneous turbulence and eventually decay towards a k⁻³ slope at higher wavenumbers. This behavior is consistent with the DNS results of Brachet et al. (1983), who showed, using DNS, that this transition occurs around k=60.

Figure 10 Kinetic energy spectra obtained using (a) 128³ points and SL, (b) 128³ points and Vreman, (c) 192³ points and SL, and (d) 192³ points and Vreman.

5.7 Barbados Oceanographic and Meteorological Experiment (BOMEX)

BOMEX features a shallow-cumulus-topped boundary layer as described in Holland and Rasmusson (1973). The setup of this test follows Siebesma et al. (2003). The initial profiles are characterized by a well-mixed sub-cloud layer below 500 m, a cumulus layer between 500 and 1500 m, an inversion layer up to 2000 m, and a free troposphere above. Large-scale forcing includes prescribed large-scale subsidence, horizontal advective drying, radiative cooling, and pressure gradient effects. Sensible and latent heat fluxes at the surface are prescribed to $\mathrm{SHF}=\mathit{n}\cdot \left(\mathit{\rho }{\mathit{J}}_{\mathrm{sfc}}\right)=\mathrm{9.5}$ W m⁻² and $\mathrm{LHF}=\mathit{n}\cdot \left(\mathit{\rho }{\mathit{D}}_{\mathrm{sfc}}\right)=\mathrm{147.2}$ W m⁻². Additional detail on the application of boundary conditions is presented in Appendix A. The domain, $\mathrm{\Omega }=\mathrm{6400}\times \mathrm{6400}\times \mathrm{3000}$ m³, is doubly periodic in the x and y directions. A Rayleigh sponge layer (see Appendix A2 for details) is applied along the z direction to damp upward-propagating gravity waves. On the bottom surface, a momentum drag forcing is applied (see Appendix A1). The effective horizontal and vertical resolutions are $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{100}$ m and Δz=40 m, respectively. The simulation time is 6 h.

Figure 11 shows the vertical profiles of the domain mean thermodynamic and turbulence properties over the last hour of the simulations. Figure 12 shows the time series of liquid water paths (LWPs), cloud cover, and turbulence kinetic energy. The time-averaged results of the vertical profiles of θ_l,q_l, cloud fraction, and $q_{\mathrm{v}}=q_{\mathrm{t}}-q_{\mathrm{l}}$ during the last hour are in good agreement with the same quantities presented by Siebesma et al. (2003). The SGS model does not have much effect on the simulation characteristics, except that the Vreman SGS model produces a stronger peak in the variance of vertical velocity, $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$, near the cloud top. Although the difference is mild, it is possibly due to the low dissipation nature of the Vreman's model. Excluding the first hour of flow spin-up, the results compare well with PyCLES (Pressel et al., 2015) and fall within the ensemble range shown in Fig. 2 of Siebesma et al. (2003). PyCLES solves the anelastic conservation equations with entropy as the prognostic thermodynamic variable using weighted essentially non-oscillatory (WENO) numerics. Differences in the initial conditions (as seen in the initial LWP in Fig. 12) are likely contributing to the initial mismatch in horizontally averaged properties of moisture and turbulence. Details on the computation of horizontally averaged profiles can be found in Appendix B.

A large-domain simulation of BOMEX with effective horizontal resolution $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{50}$ m and vertical resolution Δz=20 m in a $\mathrm{30}\times \mathrm{30}\times \mathrm{3.6}$ km³ domain was executed using 16 GPUs on the Google Cloud Platform. The simulation was executed using 1D IMEX time integration (1D implicit in the vertical direction and 2D explicit in the horizontal direction) at maximum horizontal advective Courant number C=0.9. A visualization of instantaneous shallow cumulus structures is shown in Fig. 13.

Figure 11 BOMEX. Profile of the mean state of liquid potential temperature, total specific humidity, cloud fraction, liquid water specific humidity, and variance of the vertical velocity fluctuations averaged along the last hour of the simulation. The solutions with PyCLES and ClimateMachine were calculated with effective grid resolution $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{100}$ m and Δz=40 m. Results calculated with Vreman (blue lines) and SL SGS (orange lines) models in ClimateMachine are compared against results of PyCLES (grey lines) in its paired SGS modality using the SL model.

Figure 12 BOMEX. From left to right, time series of horizontally averaged LWP, cloud cover, and turbulence kinetic energy diagnosed from ClimateMachine using Vreman and SL. These results are consistent with ensemble results presented in Fig. 2 of the intercomparison study by Siebesma et al. (2003). Line colors are as in Fig. 11.

Figure 13 BOMEX. Instantaneous visualization of the shallow cumulus structures on a $\mathrm{30}\mathrm{km}\times \mathrm{30}\mathrm{km}\times \mathrm{3.6}$ km domain. The white shading of the volume rendering of the cloud corresponds to a maximum value of $q_{\mathrm{l}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{4}}$ kg kg⁻¹; a maximum $q_{\mathrm{t}}=\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{2}}$ kg kg⁻¹ is shown on the bottom surface in light yellow shading. The simulation was executed using 16 GPUs on the Google Cloud Platform.

6 Verification of conservation property of prognostic variables

We respectively define the time-dependent normalized total mass and energy changes as

$$\begin{array}{}\text{(66)}& \mathrm{\Delta }M\left(t\right)={\displaystyle \frac{{\int }_{\mathrm{\Omega }}\left[\mathit{\rho }\left(t\right)-\mathit{\rho }\left(t_{\mathrm{0}}\right)\right]\mathrm{d}\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}\mathit{\rho }\left(t_{\mathrm{0}}\right)\mathrm{d}\mathrm{\Omega }}}\end{array}$$

and

$$\begin{array}{}\text{(67)}& \mathrm{\Delta }\mathit{\rho }{e}^{\mathrm{tot}}\left(t\right)={\displaystyle \frac{{\int }_{\mathrm{\Omega }}\left[\mathit{\rho }{e}^{\mathrm{tot}}\left(t\right)-\mathit{\rho }{e}^{\mathrm{tot}}\left(t_{\mathrm{0}}\right)\right]\mathrm{d}\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}\mathit{\rho }{e}^{\mathrm{tot}}\left(t_{\mathrm{0}}\right)\mathrm{d}\mathrm{\Omega }}},\end{array}$$

where t₀ indicates the initial time and Ω is the full domain. Figure 14 shows ΔM(t) and Δρe^tot(t) for a 1 h simulation of a moist rising thermal bubble. The SL eddy viscosity model was used to represent under-resolved diffusive fluxes. This simulation was run using free-slip boundary conditions with adiabatic walls. We note that the loss of energy and mass in the system is contained to 𝒪(10⁻¹⁵): that is, numerical roundoff error. This result highlights a key benefit of the general formulation of the prognostic conservation equations in flux form, which guarantees conservation properties up to source or sink contributions.

Figure 14 Time evolution of the mass and total energy loss (relative change when compared against initial conditions) for a moist thermal bubble simulation. Blue line: relative change in energy; orange line: relative change in mass.

7 CPU strong scaling

Demonstration of favorable scaling capabilities across multiple hardware types is critical to the utility of ClimateMachine as a competitive tool for large-eddy simulations. Toward this, we first examine strong scaling on CPU architectures. The rising thermal bubble problem described in Sect. 5.2 is used as the test problem, with its domain extents modified to form an 8.192 km³ cube, with an effective nodal resolution of 32 m to ensure that CPU memory on a single rank is maximally loaded. For tests with multiple MPI ranks, each rank resides on a unique node to ensure communication overhead is appropriately represented; in practice, one would expect to use multiple ranks per node. Scaling across multiple threads is not assessed in the present work. Figure 15 shows the speed-up in time to solution for 10 time integration steps of the test problem with $N_{\mathrm{ranks}}\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{4},\mathrm{8},\mathrm{16},\mathrm{32}\mathit{\}}$ for both dry and moist simulations. We exclude checkpoint, diagnostic, and periodic runtime output steps from time-to-solution measurements. In both dry and moist simulations, we see a speed-up of approximately 19.7 when using 32 ranks compared with the corresponding single-rank simulation.

A single-rank GPU run of the test problem on a 6.144 km³ domain with 32 m effective resolution (restricted by GPU memory capacity) has a wall-clock time for 10 integration steps of 314 s. The wall-clock time for a 32-rank CPU run was 449 s for a problem 2.37 times larger.

This provides an estimate for a comparison between CPU and GPU hardware performance. However, the balance between memory bandwidth limits and computing operation limits guides the maximum scaling possible on the GPU hardware relative to its CPU counterpart, so this cannot be interpreted as a direct comparison across hardware types. Based on the present results, we conclude that it is more feasible to pursue strong scaling improvements on CPU hardware than on GPU hardware. Further optimization and exploration of scaling in ClimateMachine are ongoing work. Additional details on the hardware used for scaling tests can be found in Appendix C.

Figure 15 Speed-up of the time integration (solver) step relative to the time to solution for a single-rank simulation of the rising thermal bubble problem in an 8.192 km³ domain with a uniform effective resolution of 32 m (with fourth-order polynomials). Blue circles: dry simulation; orange squares: moist simulation.

8 GPU weak scaling

To test the multi-GPU scalability of ClimateMachine, we first execute a BOMEX setup that is sufficiently large to saturate one GPU. The single-GPU execution represents the baseline from which we calculate the average time per time step denoted by t₁. We then expand the domain size to match an increase in the number of GPUs and measure the average time per time step. Our scaling is then obtained as the ratio $t_{\mathrm{1}}/t_n$, with t_n being the average time per time step obtained with n GPUs. The results are obtained using up to 16 NVIDIA Tesla V100 GPUs running Julia version 1.4.2, CUDA 10.0, and CUDA-aware OpenMPI 4.0.3. Figure 16 shows excellent weak scaling for up to 16 GPUs on Google Cloud Platform resources. Over 95 % weak scaling was achieved with 1D IMEX time integration and over 98 % for the simulation with explicit time stepping. This is an encouraging result and supports the ability to prototype smaller problem setups and deploy larger simulations in the ClimateMachine limited-area configuration at identical resolutions without significantly compromising the time to solution.

Figure 16 BOMEX weak scaling using 1D IMEX versus fully explicit time integration.

9 Conclusions

This paper introduced and assessed the LES configuration of ClimateMachine, a new Julia-language simulation framework designed for parallel CPU and GPU architectures. Notable features of this LES framework are the following:

• conservative flux form model equations for mass, momentum, total energy, and total moisture to ensure global conservation of dynamical variables of interest (up to nonconservative source or sink processes);

• discontinuous Galerkin discretization with element-wise evaluation of the approximations to volume and interface integrals, resulting in reduced time to solution due to MPI operations;

• application of model equations to the solution of benchmark problems in typical LES codes, including atmospheric flows in the shallow cumulus regime (BOMEX); and

• demonstration of strong scaling on CPUs with up to 32 MPI ranks (speed-up of 19.7 in time to solution) and weak scaling with up to 16 GPUs (95 %–98 %) in both dry and moist simulation configurations.

Appendix A: Boundary conditions

A1 Solid walls and wall fluxes

A1.1 Momentum

Rigid surfaces are considered impenetrable such that the wall-normal component of the velocity vanishes at rigid boundaries by imposing $\mathit{u}\cdot \mathit{n}=\mathrm{0}$. The viscous sublayer is not explicitly resolved, and a momentum sink is applied to model the effect of wall-shear stresses. While the wall-normal advective momentum flux vanishes, the wall-normal viscous or SGS momentum flux, also known as the bulk surface stress (units of Pa),

$$\begin{array}{}\text{(A1)}& \mathit{n}\cdot \left(\mathit{\rho }\mathit{\tau }\right)=-\mathit{n}\cdot \left[\mathrm{2}\mathit{\rho }{\mathit{\nu }}_{\mathrm{t}}\mathit{S}\left({\mathit{u}}_{\mathrm{p}}\right)\right],\end{array}$$

is not necessarily negligible. Here, $\mathit{S}\left({\mathit{u}}_{\mathrm{p}}\right)=(\mathrm{\nabla }{\mathit{u}}_{\mathrm{p}}+\mathrm{\nabla }{\mathit{u}}_{\mathrm{p}}^T)/\mathrm{2}$ is the strain rate tensor of the near-surface wall-parallel velocity, u_p. We note that, throughout Appendix A, n refers to the inward-pointing normal vector at domain boundaries, distinct from the prior definition of the element-interface normal vector in Sect. 3.

In the case of free-slip conditions at a solid surface (indicated by subscript “sfc”), there is no viscous or SGS momentum transfer between the atmosphere and the surface, such that

$$\mathit{n}\cdot \left(\mathit{\rho }\mathit{\tau }\right){|}_{\mathrm{sfc}}=\mathrm{0}.$$

Because the momentum flux tensor depends linearly on velocity derivatives, this amounts to homogeneous Neumann boundary conditions on velocity components parallel to the surface. On the other hand, viscous drag is imposed by the classical aerodynamic drag law,

$$\begin{array}{}\text{(A2)}& \mathit{n}\cdot \left(\mathit{\rho }\mathit{\tau }\right)=-\mathit{\rho }{C}_{\mathrm{d},\mathrm{int}}\Vert {\mathit{u}}_{\mathrm{p},\mathrm{int}}\Vert {\mathit{u}}_{\mathrm{p},\mathrm{int}},\end{array}$$

where the quantities with subscript int are evaluated at an interior point x_int adjacent to the surface. The drag coefficient,

$$C_{\mathrm{d},\mathrm{int}}=C_{\mathrm{d}}(\mathit{Y};{\mathit{x}}_{\mathrm{int}}),$$

can depend parametrically on state variables Y(x,t) and on the position x_int of the interior point relative to the surface. In the present implementation, the plane of interior points relevant to boundary flux evaluation is interpreted as the first layer of interior nodes in the surface-adjacent elements. The drag law boundary condition amounts to inhomogeneous Neumann boundary conditions on velocity components parallel to the surface.

A1.2 Specific humidity

As for momentum, the advective specific humidity fluxes normal to a rigid surface vanish, but the diffusive or SGS specific humidity fluxes normal to the surface may not vanish. Normal components of SGS fluxes of condensate, q_l, are generally set to zero at boundaries,

$$\begin{array}{}\text{(A3)}& \mathit{n}\cdot \left(\mathit{\rho }{\mathit{d}}_{q_{\mathrm{l}}}\right){|}_{\mathrm{sfc}}=\mathrm{0},\end{array}$$

implying homogeneous Neumann boundary conditions ($\mathit{n}\cdot \mathbf{\nabla }{q}_{\mathrm{l}}=\mathrm{0}$) on the condensate specific humidities. The total SGS specific humidity flux then reduces to the vapor flux at the surface. SGS turbulent deposition of condensate (fog) on the surface can in principle occur; representing this would require nonzero condensate fluxes at the surface. With the assumption of zero condensate boundary fluxes, we have

$$\begin{array}{}\text{(A4)}& \mathit{n}\cdot \left(\mathit{\rho }{\mathit{d}}_{q_{\mathrm{t}}}\right){|}_{\mathrm{sfc}}=\mathit{n}\cdot \left(\mathit{\rho }{\mathit{d}}_{q_{\mathrm{v}}}\right){|}_{\mathrm{sfc}},\end{array}$$

where evaporation (measured in kg m⁻² s⁻¹), or condensation if negative, is given by

$$\begin{array}{}\text{(A5)}& E=\mathit{n}\cdot \left(\mathit{\rho }{\mathit{d}}_{q_{\mathrm{v}}}\right){|}_{\mathrm{sfc}}=-\mathit{n}\cdot \left(\mathit{\rho }{\mathcal{D}}_{\mathrm{t}}\mathbf{\nabla }{q}_{\mathrm{t}}\right){|}_{\mathrm{sfc}}.\end{array}$$

Evaporation can be zero (water impermeable) or it can be given as a function of Y at the surface according to

$$\begin{array}{}\text{(A6)}& E=\mathit{n}\cdot \left(\mathit{\rho }{\mathit{d}}_{q_{\mathrm{v}}}\right){|}_{\mathrm{sfc}}=E(\mathit{Y};{\mathit{x}}_{\mathrm{sfc}},t),\end{array}$$

which, numerically, translates into an inhomogeneous Neumann boundary condition on the vapor specific humidity.

A1.3 Energy

As for momentum and humidity, the advective energy fluxes normal to a rigid surface vanish, but the diffusive or SGS flux of total enthalpy, h^tot, normal to the surface (units of W m⁻²),

$$\begin{array}{}\text{(A7)}& \mathit{n}\cdot \mathit{\rho }(\mathit{J}+\mathit{D})=-\mathit{n}\cdot \left(\mathit{\rho }{\mathcal{D}}_{\mathbf{t}}\mathbf{\nabla }{h}^{\mathrm{tot}}\right),\end{array}$$

may not vanish. Because the kinetic energy contribution to the total enthalpy flux near a surface is usually 3–4 orders of magnitude smaller than the thermal and potential energy components, it is generally neglected so that the total enthalpy flux J+D reduces to a flux of moist static energy $\mathrm{MSE}=h+\mathrm{\Phi }$. The surface can be insulating, in which case the SGS transfer of total enthalpy between the atmosphere and the surface is zero:

$$\begin{array}{}\text{(A8)}& \mathit{n}\cdot \mathit{\rho }(\mathit{J}+\mathit{D}{)}_{\mathrm{sfc}}=\mathrm{0},\end{array}$$

which, from a numerical point of view, translates to a homogeneous Neumann condition on the total enthalpy ($\mathit{n}\cdot \mathbf{\nabla }{h}^{\mathrm{tot}}=\mathrm{0}$) or, by neglecting kinetic energy, on MSE such that ($\mathit{n}\cdot \mathbf{\nabla }\mathrm{MSE}=\mathrm{0}$). With known MSE, the moist static energy flux (MSEF) at the surface is given by

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}\mathit{n}& \cdot \mathit{\rho }(\mathit{J}+\mathit{D}{)}_{\mathrm{sfc}}=\mathrm{MSEF}(\mathit{Y};{\mathit{x}}_{\mathrm{sfc}},t)\\ & =-\mathit{\rho }{C}_{\mathrm{h},\mathrm{int}}\Vert {\mathit{u}}_{\mathrm{p},\mathrm{int}}\Vert ({\text{MSE}}_{\mathrm{int}}-{\text{MSE}}_{\mathrm{sfc}}).\end{array}\end{array}$$

with C_h,int the thermal exchange coefficient.

The value of $\mathit{\rho }\mathit{n}\cdot (\mathit{J}+\mathit{D})$ can also be assigned by the summation of given LHF and SHF, as done in the case of BOMEX described in Sect. 5.

A2 Non-reflecting top boundary

To prevent the reflection of fast, upward-propagating gravity waves at the top boundary, a Rayleigh damping sponge is added to the right-hand side of the momentum equation (see Sect. 5). The damping in the momentum equations takes the form

$$\begin{array}{}\text{(A10)}& {\displaystyle \frac{\partial \mathit{\rho }\mathit{u}}{\partial t}}=\mathrm{\dots }-{\mathit{\tau }}_{\mathrm{s}}\mathit{\rho }(\mathit{u}-{\mathit{u}}_{\mathrm{relax}}),\end{array}$$

where u_relax is a specified background velocity to which the flow is relaxed within the absorbing layer, with a characteristic relaxation coefficient τ_s. Of the many alternative options known for τ_s (e.g., Durran and Klemp, 1983), the default in ClimateMachine is

$$\begin{array}{}\text{(A11)}& {\mathit{\tau }}_s=\mathit{\alpha }{\sin }^{\mathit{\gamma }}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{z-z_{\mathrm{s}}}{z_{\mathrm{top}}-z_{\mathrm{s}}}}\right)\text{ for }z>z_{\mathrm{s}},\end{array}$$

where the absorbing sponge layer starts at z=z_s, γ is a positive even power typically set to 2, and α>0 is a relaxation coefficient typically of order 𝒪(1 s⁻¹).

A3 Numerical implementation

For the boundary velocity corresponding to the impenetrable wall condition, we use the following reflecting condition:

$$\begin{array}{}\text{(A12)}& {\displaystyle }{\mathit{u}}_{\mathrm{bc}}& {\displaystyle }={\mathit{u}}^{-}-(\mathit{n}\cdot {\mathit{u}}^{-}\mathit{n}),\text{(A13)}& {\displaystyle }{\mathit{u}}^{+}& {\displaystyle }=\mathrm{2}{\mathit{u}}_{\mathrm{bc}}-{\mathit{u}}^{-}.\end{array}$$

The no-slip condition follows from Eq. (A12) by setting all components of u_bc=0. Boundary conditions on a scalar χ are similarly specified as follows:

$$\begin{array}{}\text{(A14)}& {\mathit{\chi }}^{+}=\mathrm{2}{\mathit{\chi }}_{\mathrm{bc}}-{\mathit{\chi }}^{-}.\end{array}$$

Nonzero mass flux boundary conditions can be imposed at penetrable or free surfaces by applying the transmissive boundary condition:

$$\begin{array}{}\text{(A15)}& {\mathit{u}}^{+}={\mathit{u}}^{-}.\end{array}$$

Diffusive fluxes are applied by a direct specification of the wall-normal fluxes, and over-specified boundary conditions are avoided by using only the interior (⁻) gradients and first-order fluxes.

Appendix B: Statistics

Since the flow is compressible, we use density-weighted Favre averages following Canuto (1997) when computing horizontally averaged statistics. For a scalar ϕ, the density-weighted average $\overline{\mathit{\varphi }}$ at a given height level z is defined by

$$\begin{array}{}\text{(B1)}& \overline{\mathit{\varphi }}={\displaystyle \frac{\langle \mathit{\rho }\mathit{\varphi }\rangle }{\langle \mathit{\rho }\rangle }},\end{array}$$

where 〈⋅〉 denotes a horizontal mean. All calculations of horizontal statistics are done on the DG nodal mesh to avoid introducing interpolation errors (Yamaguchi and Feingold, 2012), with metric terms accounted for in the descriptions of diagnostic variables. The density-weighted vertical eddy flux for a variable ϕ is defined by

$$\begin{array}{}\text{(B2)}& \stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\varphi }}^{\prime }}={\displaystyle \frac{\langle \mathit{\rho }{w}^{\prime }{\mathit{\varphi }}^{\prime }\rangle }{\langle \mathit{\rho }\rangle }},\end{array}$$

where

$$\begin{array}{}\text{(B3)}& {\mathit{\varphi }}^{\prime }=\mathit{\varphi }-\overline{\mathit{\varphi }}\end{array}$$

denotes the deviation from the density-weighted mean. The variance can then be defined as

$$\begin{array}{}\text{(B4)}& \stackrel{\mathrm{‾}}{{{\mathit{\varphi }}^{\prime }}^{\mathrm{2}}}={\displaystyle \frac{\langle \mathit{\rho }{{\mathit{\varphi }}^{\prime }}^{\mathrm{2}}\rangle }{\langle \mathit{\rho }\rangle }}.\end{array}$$

Appendix C: Hardware

This section summarizes the hardware characteristics for the primary computing resources used in tests throughout this paper. This is particularly relevant to the data presented in Sect. 7. Compute nodes for the CPU tests were 14-core Intel Xeon (2.4 GHz), with a maximum memory capacity of 1.5 TB. GPU nodes on this cluster were of 14-core Intel Broadwell (2.4 GHz) type with 28 cores per node and 256 GB of memory per node. Compute nodes on the Google Cloud Platform leverage Tesla V100 GPUs available for general-purpose use.

Appendix D: Supplementary results

This section provides additional information on the comparison of the density current benchmark in ClimateMachine with existing literature references in Table D1.

(Giraldo and Restelli, 2008)(Giraldo and Restelli, 2008)(Marras et al., 2015)(Marras et al., 2015)(Marras et al., 2015)(Marras et al., 2015)(Marras et al., 2015)(Marras et al., 2015)(Marras et al., 2015)(Marras et al., 2013)(Marras et al., 2013)(Marras et al., 2013)(Marras et al., 2013)(Ahmad and Lindeman, 2007)(Straka et al., 1993)

Table D1 Summary of frontal locations for the density current test case from existing literature. Results tabulated are of the front location at t=900 s. The results are reported for the following models: ClimateMachine with SL and Vreman, finite-element model (FEM) variational multiscale stabilization (VMS), f-wave, filtered spectral elements (SEs), filtered discontinuous Galerkin (DG), and the piecewise parabolic method (PPM). NUMA is the Nonhydrostatic Unified Model of the Atmosphere.