The Flux-Differencing Discontinuous Galerkin Method Applied to an Idealized Fully Compressible Nonhydrostatic Dry Atmosphere

Dynamical cores used to study the circulation of the atmosphere employ various numerical methods ranging from finite-volume, spectral element, global spectral

10.1029/2022MS003527 3 of 24 see, for example, Boyd (2001) for numerous reasons why it can be better to create a well-posed mathematical problem and use an optimally convergent numerical method. A robust numerical method saves human time since it is common in geophysical simulations to include the minimal necessary dissipation for stability; see Winters et al. (2021) for comments with respect to engineering and astrophysical examples. Tuning numerical filters to achieve a desired level of fidelity requires substantial effort, one that must often be repeated upon any change to model configuration. The automation of this effort through a well-chosen numerical method allows model developers to focus on the physics.
The rest of the paper is organized as follows. In Section 2 we give a brief overview of what makes FDDG different from traditional methods. We describe the mimetic discrete properties of the FDDG formulation and our choice of numerical fluxes. This section requires knowledge of spectral collocation methods and is safely skipped for those more interested in the results.
In Section 3 we give three examples of applying the FDDG method to the compressible Euler equations with gravity in both rotating and non-rotating reference frames, which we take as a model for dry atmospheres. First, we examine a large-eddy simulations (LES) of a dry convective boundary layer in a box with rigid walls at the top and bottom and doubly periodic horizontal boundary conditions. Second, we explore the use of FDDG for an idealized dry global Earth configuration (Held & Suarez, 1994). Third, we perform a simulation of an atmosphere in a "small-planet" configuration where the scale separation between convective scales and large scales is reduced (Kuang et al., 2005;Wedi & Smolarkiewicz, 2009). The same equation set and computational kernels are used for all three cases.

A Brief Overview of Flux-Differencing
A detailed explanation of the flux-differencing method used in this work is in Chan (2018) and Waruszewski et al. (2022), but here we focus on a heuristic overview. Furthermore we concentrate on the aspect that makes flux-differencing different from traditional DG methods, the volume term calculations. The core idea behind flux-differencing comes from recognizing that many properties of continuous calculus, such as the product or chain rule, fail to hold upon discretization. Thus algebraically equivalent formulations of an equation set yield different numerical discretizations with different discrete properties, for example, Zang (1991).
We briefly digress to examine why one would expect different numerical properties. The core discrepancy between the continuum calculus and the discrete calculus is the continuum identity for all smooth functions u(x). We choose this example because it is the simplest manifestation of aliasing errors. A grid that represents a polynomial of degree n cannot represent a polynomial of degree n + 1, as would be the case upon multiplication by x.
We illustrate the difference between continuous and discrete calculus using discrete operators. We use three grid points x i at points x 1 = −1, x 2 = 0, x 3 = 1 and let D represent the differentiation matrix corresponding to those points. The differentiation matrix is defined by requiring a discrete power rule, that is, −1 = ∑ ′ ′ ( ′ ) for k = 1, 2 and ∑ i′ D ii′ = 0. The latter requirement states that the row sum of the matrix is zero, or equivalently, that the derivative of a constant is zero. The derivative operator D and position operator X = Diagonal(x) are (2) The discrete analog of Equation 1 is where u is an arbitrary vector with components u = [u 1 , u 2 , u 3 ]. Equation 3 implies − = where is the identity operator, yet which is not the identity matrix. This discrepancy is at the heart of why diverse numerical formulations of equations have different properties.
To make connections with partial differential equation discretizations we use Burgers' equation as a further example, where we recast the right hand side in multiple ways to emphasize potential discretizations. The conservative and advective forms are familiar, and the entropy conservative formulation, in this case a weighted average of the previous two, satisfies a discrete entropy inequality. Satisfying the discrete entropy inequality yields robust numerical solutions and is the most stable of the three potential discretizations. See G. J. Gassner and Winters (2021). We present the differences between the formulations in Equation 5 pictorially in Figure 2. We see that the different formulations, although qualitatively similar, yield quantitatively different tendencies. Upon grid refinement the discrepancies between the discretizations are eliminated, but we seldom have the luxury of resolving all the necessary scales and thus choose to emphasize the differences at a coarse resolution.
All of the tendencies are, in this case, split-forms, an algebraic rewriting of the equations. All split-forms and entropy stable methods involve flux-differencing, but entropy stable methods need not be a split-form. Entropy stability requires satisfying a discrete entropy equality for a given system, which for Burgers' equation is a weighted average but in general can require different kinds of averaging, for example, Winters et al. (2021).
We now rewrite the tendencies of Equation 5 to establish connections with flux-differencing methods. We focus on the "volume" (i.e., within each cell) contribution to the derivatives, since the "interface" contributions (i.e., across neighboring cell interfaces) are handled by standard methods. Letting D ii′ be a collocation differentiation matrix with components i, i′, and u i denote the ith component of the solution vector, the components of the tendency  for Equation 5 are written as We use the primed indices, such as i′, as dummy indices for summation. We rewrite Equations 6-8 in a numerical flux formulation through the Rosetta stone provided by G. J. Gassner et al. (2016a). Let {⋅} be shorthand notation for Since the derivative of a constant is zero, that is, the row sum of a differentiation matrix is zero, We write the conservative form of the equations as Furthermore, we observe to write the advective form as and the entropy stable form as The terms acted upon by the differentiation matrix are now viewed as flux-estimates, that is, All of these numerical fluxes satisfy the symmetry condition and the consistency condition In general, the consistency condition yields the underlying flux of a conservation law.
Interpreting the different split-forms of Equation 5 as different choices in numerical fluxes as in Equation 18 through Equation 20 is more than algebraic sorcery. All of the discretizations are conservative. In general, different flux-differencing methods are conservative if there is an underlying conservation law and the discretized operators satisfy appropriate properties. See G. J. Gassner (2013), G. J. Gassner et al. (2016a), Winters et al. (2021), G. J. Gassner and Winters (2021) for insights on why equations remain conservative for operators that satisfy a summation by parts property and Chan (2018), Waruszewski et al. (2022) for more general operators. Deciding between which form of equations to use is not a matter of guesswork or experimentation; rather, they are derived based on which discrete properties one wants to preserve. Deciding between different forms does not need to be a binary decision. In principle, one can use a "weighted" methodology and dynamically pick between different numerical fluxes to satisfy given criteria, similar to weighted essentially non-oscillatory schemes (Liu et al., 1994). Another option is to switch between different flux formulations at a fixed schedule in time for efficiency.
Transitioning from split forms to flux-differencing allows for extra flexibility in the design of conservative numerical methods. In particular, it is possible to take other averages, including logarithmic or density-weighted averaging, for example, and retain conservation as well as higher-order accuracy, see Fisher and Carpenter (2013). Furthermore, one can refrain from using a higher-order method within a control volume and use lower-order methods instead. This perspective is a valuable consideration to keep in mind when implementing subgrid-scale parameterizations or adaptive shock-capturing schemes.
The generalization to multiple dimensions, systems of equations, vector fields, and curvilinear coordinates is straightforward under a proper interpretation of the averaging operator {⋅} and index juggling. The most general method is outlined by Chan (2018) and Waruszewski et al. (2022), but here we focus on a simplified two-dimensional scenario. For example, let u ij represent the nodal values of a two dimensional vector u whose components are =̂+̂ with similar considerations for a vector field v in an element of size [−1, 1] 2 , then where and represent the components of the differentiation matrix in the x-direction and y-direction, respectively. Similarly [∇ ⋅  ] and [∇ ⋅  ] represent the flux-divergence of the tensor flux in the x-direction and y-direction, respectively. In summary, the averaging prescription changes to "average along the coordinate direction" depending on which derivative is being taken. There are restrictions in the treatment of metric terms for curvilinear coordinates, see Winters et al. (2021), Waruszewski et al. (2022) for details.
The flux-differencing formulation extends to non-conservative terms such as ρ∂ x ϕ as shown by Renac (2019) and Waruszewski et al. (2022). To see this we introduce the jump notation Letting ρ i be a discretized density field and Φ i be a discretized geopotential, we get which is a discretization of ρ∂ x Φ. We warn that the "numerical flux"  = ({ } − ⟦ ⟧)⟦Φ⟧ does not satisfy the symmetry condition Equation 21 and thus we will call it a "non-conservative flux." We include this as a part of our total numerical flux, with the observation that any "jump terms" in the volume contribution to the numerical flux are distinguished as being associated with non-conservative terms.

Numerical Fluxes for the Compressible Euler Equations With Gravity
The prognostic variables of the compressible Euler equations are density, momentum, and, as the thermodynamic variable, total energy. The equations are where Φ is the geopotential,  are momentum sources (e.g., the Coriolis force),  constitutes sources of energy (e.g., radiation), and is the rank-3 identity matrix. Total energy is defined as the sum of kinetic, potential, and internal energy, where c v is the specific heat capacity of dry air at constant volume. Temperature is diagnosed from the prognostic variables and pressure using the ideal gas law, that is, This set of equations includes processes often filtered out in atmospheric general circulation models (AGCMs), such as sound waves. Retaining additional physics is key if the model is used for coarse resolution AGCM simulations, cloud-resolving high-resolution AGCM simulations, and high-resolution LES of boundary layers. The flexibility is especially crucial for simulating other planetary bodies or analogous "small-planet" versions of Earth.
To solve the compressible Euler equations in three-dimensional domains, we use the FDDG formulation of Chan (2018) and Waruszewski et al. (2022) and construct metric terms as outlined by D. A. Kopriva (2006). See the review by G. J. Gassner and Winters (2021) for a general overview of the FDDG method.
The choice of numerical flux is critical in guaranteeing the stability of the simulations. As mentioned, FDDG allows for a selection of numerical fluxes for the interior of the control volume. In addition, there is flexibility in the choice of numerical flux for any interface between elements, that is, a flux along the gravity-aligned direction need not be the same as a flux orthogonal to the direction of gravity.
A kinetic energy preserving (KEP) volume flux greatly enhances the flow's nonlinear stability; see G. J. Gassner et al. (2016a) for an explanation of this property. This is especially important for simulating highly underresolved turbulent flows, as is typical in geophysical fluid dynamics. We find the KEP property to be the key feature that greatly increases the robustness of simulations. Stated succinctly, a numerical flux satisfies the KEP property if the discrete kinetic energy equation mimics the continuous kinetic energy equation. The importance of preserving the discrete algebraic properties of the kinetic energy equation has been commented on before (Zang, 1991).
Numerical fluxes that do not satisfy the KEP property can have terms in the discrete kinetic energy equation that correspond to energy injection due to transport, a manifestation of aliasing errors. It is serendipitous that there are a large class of numerical fluxes that satisfy this property, but it is especially worth noting that traditional DG methods do not have the KEP property when applied to geophysically relevant simulations, leading to stability problems in underresolved flows. In order to control this error, past methods had to use numerical filters, explicit dissipation, or overintegration strategies. None of these corrections are necessary if one just simply uses an FDDG formulation that automatically satisfies the KEP property.
Due to the presence of gravity in geophysical flows, it is natural to require a flux with a kinetic plus potential energy preserving property, which we will denote as KPEP. This choice allows for the discrete equations to be mimetic of the kinetic + potential energy equations (in the absence of sources in the momentum equation) ( We build a discretization that satisfies this property in two steps. The first step is to choose a KEP flux for the conservative terms. The second step is to modify the non-conservative gravity term so that the discrete equations are mimetic. Specifically, we choose the Kennedy-Gruber flux (Kennedy & Gruber, 2008) with the following modification to the gravity source term for all simulations, where the fluxes  ,  , and  are the associated numerical fluxes for density, momentum, and total energy. The modification to the gravity source term was motivated by combining the entropy stable scheme of Waruszewski et al. (2022) with the Kennedy-Gruber flux. Ultimately the justification is that the flux satisfies the KPEP property, as shown in A3 in Appendix A. The non-conservative flux is equivalent to rewriting the gravity source term as For stability the KEP property is the key, but without the KPEP property there is a spurious source of gravitational energy. The KEP flux works against this spurious source to stabilize the simulations in the absence of the KPEP property. Thus we use a KPEP flux to prevent any potential issues.
Less innovation is required for the interface fluxes since simulations are robust with respect to many choices, for example, a central flux with a Rusanov penalty suffices for all simulations; however, we use the fluxes from Equations 40-42 as the "central" part of the interface flux and add a penalty term either as a Roe flux or a Rusanov flux (Hesthaven & Warburton, 2007;Roe, 1981). The specific choice is given in the respective sections. Further numerical details and formulas are given in Appendix A.

Numerical Experiments
Simulations are performed in a Julia-based open-source codebase that can exploit heterogeneous and distributed CPU/GPU architectures (Besard, Churavy, Edelman, & Sutter, 2019;Besard, Foket, & De Sutter, 2019;Bezanson et al., 2017). Although the DG method is well suited for parallel-computing architectures (D. S. Abdi et al., 2019;Sridhar et al., 2021), the scale of our problem allowed us to perform all simulations on a single Nvidia Titan V GPU. All plots in this text were generated using the Julia package Makie.jl (Danisch & Krumbiegel, 2021).
In the following subsections, we illustrate the ability of FDDG methods to simulate 1. Convection in a dry boundary layer. 2. The dry atmospheric circulation in an Earth-like domain. 3. The dry atmospheric circulation on a small Earth.
The domain for the first simulation is a horizontally periodic Cartesian box, for the second simulation an Earthlike thin spherical shell, and for the third simulation a spherical shell with a 20 times decreased planetary radius and increased rotation rate.
Nonetheless, the same computational kernels are used for all simulations. The connectivity between the elements and metric terms is the only change to transform from one domain to another. Explicit Runge-Kutta timestepping is used for the convective boundary layer and the small-Earth simulations. In addition, for computational efficiency, we implicitly timestep vertical acoustic and gravity wave modes in the Earth-like domain.

Dry Convection in the Atmospheric Boundary Layer
We start by simulating a dry atmospheric boundary layer. The following simulation is similar in spirit to Margolin et al. (1999), but with additional simplifications. All parameters for the simulation and their physical meaning are summarized in Table 1. We use a cubic domain of volume L 3 with periodic boundary conditions in the horizontal direction and no-flux, no-penetration boundary conditions in the vertical direction. The geopotential is Φ = gz where z is the vertical coordinate and x, y are the horizontal coordinates.
For our initial condition, we choose a linear potential temperature profile which, when combined with the ideal gas law and hydrostatic balance, implies that pressure is We also apply a random perturbation to the initially zero velocity field to induce a rapid transition to turbulence, where  is a random normal variable at each grid point. Thus the initial condition for total energy is where c v is the specific heat capacity of dry air at constant volume.
We apply a radiative forcing to drive convective instability. The resulting equations are We use a Kennedy-Gruber flux for the volume terms and a Kennedy-Gruber flux with a Roe flux penalty term for the interface numerical fluxes (Kennedy & Gruber, 2008); see Appendix A for details. For time-stepping, the fourth-order low storage 14-stage Runge-Kutta method of Niegemann et al. (2012) is employed since it offers the fastest time-to-solution. The sound waves are resolved in the simulation. We emphasize that we have not included any viscosity or diffusivity and solely rely on the numerical dissipation of the FDDG method for stability. Although we use implicit LES (ILES) for modeling convection in the interior of the domain, it is doubtful that ILES can capture wall-effects, and a more realistic simulation would need to be supplemented with wall-models.
The domain is partitioned into 24 3 elements, each of which has three-dimensional fourth-order polynomials, leading to a total of 120 3 degrees of freedom. The smallest grid spacing between interpolation points within an element is 21 m, leading to a timestep size of Δt = 0.11 s to ensure compliance with the acoustic CFL limit.
The radiative heating is strongest near the surface, leading to air parcels to become buoyant and rise. As the plumes rise, they laterally entrain air from the surrounding environment; we expect the fluid to develop a well-mixed region of potential temperature near the surface. As the plumes move through the well-mixed layer, they eventually reach a stably-stratified region and overshoot their level of neutral buoyancy. The plumes drum on the stratified layer above, developing a layer of downward potential temperature fluxes and high potential temperature variance. This process erodes the stratification, leading to diffusive growth of the well-mixed region over time.
We estimate the growth of the well-mixed region from classic energetic arguments as done by, for example, Stull (1988). First, we observe that the flux of potential temperature is approximately . We define the boundary layer height to be the height of maximum horizontally-averaged stratification. The boundary layer height at a given moment in time, t, is given by the empirical scaling law where the entrainment layer modifies the constant of proportionality. Without accounting for the entrainment layer, one derives ( ) = √ 2  Δ as in Stull (1988). Accounting for the entrainment layer seems to only modify the constant "2", for example, (Van Roekel et al., 2018), as opposed to modifying the scaling law.
Specifically, we compare the boundary layer height given by ( ) = √  Δ , with C = 3 as in A. N. Souza et al. (2020), to that of the simulation in Figure 3. We see that the simulation agrees well with the empirical scaling law. This agreement suggests that the implicit dissipation mechanisms of the FDDG method enable subgrid-scale modeling, similar to other methods such as a Smagorinsky closure or a non-oscillatory scheme (Margolin et al., 1999;Van Roekel et al., 2018). An instantaneous snapshot of the simulation after 5 hr is typified by Figure 4. The three-dimensional figure shows the mixed layer potential temperature as transparent, thereby emphasizing potential temperature anomalies. The visualization reveals the three-dimensional convective structure and small scorching plumes emanating from the surface. The top of the domain is chosen to be the height at which the horizontally averaged potential temperature flux is most negative.
To the right of the three-dimensional figure are horizontal averages of potential temperature (top), vertical potential temperature flux (middle), and potential temperature variance (bottom). The horizontal average of potential temperature displays a well-mixed layer in the bottom kilometer of the domain, capped by an entrainment layer of enhanced stratification before easing into the background stratification. The vertical advective flux exhibits the expected linear structure in the mixed layer and is negative in the entrainment region. The negative flux arises from an anti-correlation between the vertical velocity and potential temperature, associated with plumes overshooting their region of neutral buoyancy. On average, this entrainment produces a negative flux whose maximum is approximately 17% of the input heat flux Q θ . The negative flux minima is consistent with those commonly found in the literature, for example (Margolin et al., 1999;Siebesma et al., 2007;Van Roekel et al., 2018), where the most negative flux is between 10% and 20% of the heat input. The oscillations above the entrainment layer are due to gravity waves reflecting from the top of the domain. Furthermore, the plot shows that the temperature variance is largest in the entrainment layer.

Atmospheric Dynamical Core: The Held-Suarez Test
We next consider the GCM benchmark test proposed by Held and Suarez (1994), HS94 hereafter. The formulation of the problem allows for flexibility in hydrostatic versus non-hydrostatic dynamics, dissipation mechanisms, prognostic variables, and boundary conditions. We choose to use an equation set that retains fully compressible dynamics and is formulated in terms of density, total energy, and Cartesian momentum as the prognostic variables, yielding the equations where Φ = 2 −1 planet − −1 is the geopotential, = Ω̂ is the planetary angular velocity, and ̂ is the direction of the planetary axis of rotation. We do not make the traditional approximation, which assumes a thin atmospheric shell in which the distance from any point in the atmosphere to the center of the planet is taken to be equal to the planetary radius, leading to the Coriolis force having only horizontal components.
The HS94 forcing is applied to momentum and energy as follows where T equilibrium is the radiative equilibrium temperature depending on latitude (φ) and pressure σ = p/p 0 , and the parameters k v , k T are the inverse timescales for momentum damping and temperature relaxation, respectively, with = Δ and = + ( − )Δ cos 4 ( ), with Δ = max{0, ( − )∕(1 − )} . The temperature and pressure are diagnosed from total energy and the ideal gas law, The forcing terms differ only in quantitatively irrelevant aspects from the original formulation in HS94. In particular, we choose a constant pressure p 0 in the definition of σ instead of the instantaneous surface pressure.
The parameter values are summarized in Table 2.
The domain is a piecewise polynomial approximation to a thin spherical shell of radius r planet and height z top . The thin spherical domain is partitioned into curved elements and uses an isoparametric representation of the domain and the cubed sphere mapping by Ronchi et al. (1996). In essence, this choice represents the domain as a piecewise polynomial function where the order of the polynomial corresponds to the order of the discretization . The metric terms are treated as in D. A. Kopriva (2006) and satisfy the discrete property that the divergence of a constant vector field is zero, that is, the metric terms are free-stream preserving. The use of an isoparametric representation of the sphere with free-stream preserving metrics has a few subtleties. Since the vertical and horizontal directions are no longer discretely orthogonal, one must distinguish covariant and contravariant vertical directions.
We use no-flux boundary conditions for density and total energy. We use free-slip boundary conditions for the horizontal momenta and no-penetration boundary conditions for the vertical momentum. Our initial condition is a fluid that starts from rest in an isothermal atmosphere. We take the global temperature to be T I = 285 K, leading to ) and ( ) = 1 ( ).
We use implicit time-stepping in order to numerically filter vertically propagating sound waves and gravity waves. Specifically, we use the second-order Runge-Kutta IMEX scheme of Giraldo et al. (2013), but modify Equation 3.9 of their work by choosing a 32 = 1/2 for an enhanced stability region (see ARKB(2, 3, 2) in Giraldo (2020)). We use the Jacobian of both the surface and volume flux in the vertical for the implicit time-stepping component; see A2 in Appendix A for details. We linearize about the previous timestep, update the Jacobian for every column every 20 min of simulated time, and factorize it using a banded LU decomposition (Golub & Loan, 2013).
Horizontal acoustic modes then limit the timestep. The largest Mach number, the ratio of the advective speed and the soundspeed, for the flow is roughly 0.25 in this simulation.
Aside from the inherent numerical dissipation resulting from the interface flux terms and implicit time-stepping, we use no additional forms of damping such as those in Jablonowski and Williamson (2011). In particular, we do not use any form of viscosity/hyperviscosity for small-scale damping. Furthermore, we do not include any divergence damping or filters. The method remains conservative up to rounding errors from finite-precision arithmetic. For the Held-Suarez benchmark, only density is conserved since it has no sources.
We run the HS94 test case with 6 × 10 2 elements in the horizontal on an equiangular cubed sphere, eight evenly spaced elements in the vertical, polynomial order four within each element, totaling at 6 × 50 2 degrees of freedom in the horizontal and 40 degrees of freedom in the vertical. The minimum grid spacing is 120 km in the horizontal and 650 m in the vertical. We choose a timestep of 55 s to keep within the horizontal acoustic CFL limit. We discard the first 200 days of the simulation as spinup and average over the last 1,000 days, as in HS94. We gather statistics by interpolating the cubed sphere grid to spherical coordinates and converting the Cartesian momentum to spherical velocities. As usual, we denote the zonal velocity component by u, the meridional velocity by v, and the vertical velocity by w. We gather statistics in height coordinates and for plotting we use the zonal and temporal average of pressure at the equator as the height.
In Figure 5 we show the long-time average of the zonal-mean zonal wind 〈u〉, temperature 〈T〉, temperature variance 〈T′T′〉, eddy momentum flux 〈u′v′〉, eddy heat flux 〈v′T′〉, and horizontal eddy kinetic energy 0.5〈u′u′ + v′v′〉. The choice of fields is to directly compare with Figure 1 of Wan et al. (2008). The results here are in agreement with those reported in the literature (Chen et al., 1997;Held & Suarez, 1994;Ringler et al., 2000;Ullrich & Jablonowski, 2012). For example, the peak in westerly winds, temperature variance, and eddy kinetic energy are all within 10% of published results. Perhaps the largest difference is in the meridional heat transport. In our simulations, the 〈v′T′〉 = −9 K m s −1 contour remains disconnected above and below the "stretched height" = 400 hPa line. This difference could be due to the use of height coordinates for averaging rather than pressure coordinates, since a zonal average over a surface of constant height is different than that of constant pressure. For a fully compressible code it is more natural to use density-weighted averages (Favre averages), thus we also present those statistics in Figure 6. The color scale is the same as that of Figure 5, allowing for a direct comparison. We see that the density weighted statistics for mean quantities and eddy-statistics associated only with momentum are relatively unchanged with respect to the unweighted versions; however the eddy fluxes that include temperature appear to be noisier. For example, the density weighted eddy-heat flux exhibits oscillations near the equator, perhaps due to the need for longer averaging over the same time interval.

Small-Planet Held-Suarez
In addition to the typical HS94 configuration, we simulate a small planet with a large-scale climatology similar to that of HS94 by rescaling the equations in a manner similar to a DARE/hypohydrostatic rescaling of the equations as done by Kuang et al. (2005) and Pauluis et al. (2006), respectively. This rescaling is an exact similarity transformation of the hydrostatic primitive equations using the traditional approximation and thus only affects the balance between the non-hydrostatic and hydrostatic components of the flow. This test is similar to that proposed by Wedi and Smolarkiewicz (2009) with minor modifications.
We decrease the planetary radius by a factor of  = 20 compared to Earth, increase the rotation rate by a factor of  , and decrease the mass of the planet by a factor of  2 . Furthermore, we increase all relaxation timescales in the problem by a factor of  . The atmospheric height and temperature equilibrium remain the same. The parameter values are tabulated in Table 2. We will justify these choices shortly.
Changing the planetary radius, increasing the rotation rate, and keeping the same temperature equilibrium results in a planetary model with a similar thermal wind. This result is a natural consequence of the rescaling being an exact similarity transformation for the hydrostatic primitive equations. Indeed the thermal wind, u thermal scales like thermal ∼ Δ ΩΔ (61) where ΔT/ΔH is the latitudinal gradient of temperature. Observe that ΔH ∝ r planet and recall that the equilibrium temperature distribution is unchanged from the original configuration. Thus both ΔT and ΩΔH remain constant, and the resulting thermal wind is approximately the same across the two simulations. Consequently, the Rossby number Ro ≡ u thermal /(2Ωr planet ) remains the same.
Changing the planetary mass is necessary to retain an Earth-like hydrostatically balanced state. The gradient of the geopotential scales like ∇Φ ∼ −2 planet and thus the planetary mass must scale by a factor of  −2 to maintain the same force. We could have achieved a similar result by simply taking the geopotential to be Φ = gr, but we saw no need to use this linearization.
We keep the same number of grid points, 6 × 50 2 × 40 degrees of freedom, leading to a minimum grid spacing of 6 km in the horizontal and 650 m in the vertical. For the small planet, we use explicit time-stepping-the same low storage Runge-Kutta method of Niegemann et al. (2012)-which affords timesteps of size dt = 6.5 s, which corresponds to an acoustic Courant number of 3.6 in the vertical and 0.38 in the horizontal. Small timesteps are less of a limitation because planetary-scale dynamics are  = 20 times faster than Earth's. Thus we only need to simulate 60 Earth days, which corresponds to 1,200 small-planet days. The initial condition uses the same formula as before, Equation 61. We discard the first 20% of the simulation and average over the rest. Figure 7 shows that statistics are relatively unchanged with respect to those in Figure 5, except for the zonal velocity, which has a vigorous easterly flow along the equator. We attribute the change in the zonal mean climatology of the zonal velocity to the increased vertical velocity, which in turn affects the non-traditional terms in the Coriolis force; these terms are not negligible in the small planet. See Marshall et al. (1997) for an explanation of the underlying physics in the ocean context. An enhanced easterly flow in the small planet configuration has been observed before. For example, see Figure 18 of Wedi and Smolarkiewicz (2009). We confirm this statement by neglecting the non-traditional components of the planetary angular velocity, and comparing the zonal mean velocity statistics of the three simulations in Figure 8. We do not modify the metric terms thus the approximation is inconsistent, nonetheless it serves to illustrate the point. We see that the zonal mean velocity statistic of the original HS94 setup corresponds to that of the small planet with the "traditional" planetary angular velocity but not that of the small planet with the full-planetary angular velocity. This effect is a consequence of the decreased aspect ratio of the vertical versus horizontal domain in the small planet, which in turn increases the magnitude of the vertical velocity by a factor  . Stated differently, even though the full Coriolis  force is present in the Earth-like domain, the vertical velocity component is too weak to make a substantial difference, as expected for this test case.
We reemphasize no further code tuning is required to retain stability. Upon modification of the domain and appropriate parameters, the only necessary change was a reduction of timestep to stay within the acoustic CFL of the small planet. The ability to easily change planetary parameters allows for a systematic investigation of scaling laws of planetary systems with respect to rotation rates, planetary radii, and atmospheric heights.

Conclusions
We have presented the application of a FDDG method to an idealized dry atmosphere for local LES and global circulation modeling. We have shown that the statistics generated from using a fully-compressible code, with density, total energy, and Cartesian momentum as prognostic variables, are similar to other models in local and global settings. Furthermore, we did not require stabilization mechanisms outside those naturally afforded by the FDDG numerical method and time-stepping.
The main limitations of the numerical method are not associated with the spatial discretization per se but rather the need to develop efficient time-stepping strategies for modern computer architectures that can overcome limitations induced by acoustic waves, especially in the presence of topography. Different architectures may necessitate different algorithms to achieve an optimal time-to-solution. There are many approaches for obtaining a better time-to-solution that are worth exploring, for example, fully implicit time stepping (Nguyen et al., 2020) and multi-rate methods (Knoth & Wensch, 2014). Furthermore, switching between different flux-differencing methods in the vertical versus horizontal may yield larger timesteps due to better linearization properties (G. Gassner et al., 2020;Ranocha & Gassner, 2021). An alternative option is to use lower order methods, such as staggered grid finite volume or lower polynomial orders, for the implicit vertical discretization, which may yield a faster time-to-solution.
The present study is limited to an idealized dry atmosphere, and moisture, topography, and radiation are necessary for realistic simulations. Positivity-preserving methods such as those outlined in Light and Durran (2016) need to be used, and topographic effects can also be handled Baldauf (2021). Previous studies of DG methods have involved designing numerical fluxes that preserve desired discrete properties. It would also be interesting to compare candidate methods for geophysical flows.
It is possible to bridge the gap between existing parameterizations and novel numerics by leveraging the sub-cell finite-volume interpretation of the FDDG method. This interpretation is similar to using a "physics grid" as in Herrington et al. (2019) but simpler in its implementation. Another option is to develop new parameterizations that leverage the subgrid-scale shape functions of the spectral element method, akin to using a higher-order moment closure.
Flux-Differencing DG methods are an interesting alternative discretization for Earth system modeling. They enable LES modeling with its natural subgrid-scale dissipation mechanisms, allow for flexible representation of the domain, and can exploit parallel hardware architectures. Developing efficient implicit timesteping methods in order to overcome the limitations due to gravity and sound waves are a remaining challenge, but we hope that the extra robustness and higher-order accuracy provided by FDDG methods will eventually allow for an overall simpler and more accurate method.

Appendix A: Discontinuous Galerkin Details
In this appendix, we collect choices of numerical fluxes and linear models. To highlight our choices, we use the compressible Euler equations with gravity, where γ = 7/5 and Φ is the geopotential. The source terms that do not involve gradients are collocated with grid-points and require no further description.
To describe the numerical fluxes we use the same notation as G. J. Gassner et al. (2016a). Thus for a scalar field ψ with + as the "exterior" value and − as the "interior" value (G. J. Gassner et al., 2016a;Hesthaven & Warburton, 2007), we take the averaging operator {⋅} and jump operator ⟦⋅⟧ to mean The averaging and jump operators are applied componentwise for vector and tensor Cartesian fields. We point out that our definition of jump, ⟦⋅⟧, has a factor of two that is different from most other conventions.
The flux-differencing and metric term implementations are done in skew-symmetric form as outlined by Chan (2018) and Waruszewski et al. (2022). The metric terms are constructed to be free-stream preserving (D. A. Kopriva, 2006;D. Kopriva, 2009).

A1. Numerical Fluxes
We decompose the numerical flux normal to an interface between elements into two components by using the flux above as the "central" component and a penalty term, which adds dissipation in a manner similar to upwinding. For the "central part" of the interface terms we use the Kennedy-Gruber flux (Kennedy & Gruber, 2008), where is the identity matrix. The nonconservative term is inconsequential since the geopotential Φ is continuous along an interface, hence ⟦Φ⟧ = 0 on an interface. We choose different penalty terms for the vertical versus horizontal directions when evolving the compressible Euler-Equations on the sphere. Distinguishing between vertical and horizontal fluxes is natural given the anisotropy of the Earth-like computational domain: a spherical shell with radius  ( 10 4 ) kilometers and height (10) kilometers. This domain typically leads to pancake-like grid elements whose breadth is roughly 100 times its height.
In the direction associated with vertical face normals we use a Rusanov penalty whose wavespeed is based on a reference pressure and reference density. The reference density and pressure are updated every 20 min of simulated time with the instantaneous values. Specifically we add the following numerical fluxes, where { } ∞ = max{ + , − } . In the directions orthogonal to the vertical direction we use Roe fluxes.
where the averaging, {⋅} ρ is for all fields ψ except for ρ in which case The variable = ( ) is the normal vector to a point on an element face (unit vectors of the contravariant basis) and = ⋅̂ is the velocity component normal to a face.
On the boundaries of the sphere we set the density and energy fluxes to zero and for momentum we use the exterior + state and interior − state as where is the wall-normal unit vector. We then use central fluxes to compute the flux on the boundary. Equation A21 amounts to using the reflection principle on the wall-normal velocity, while also implementing no-flux boundary conditions for the tangential velocities. See Hesthaven and Warburton (2007) for further clarification on the reflection principle.

A2. Jacobian for Implicit Timestepping
To calculate the Jacobian of the compressible Euler equations with gravity it suffices to focus on the numerical flux, First we make the observation that variables such as u, e, and p are nonlinear functions diagnosed from the prognostic variables ρ, ρu, and ρe, = , = , and = ( − 1) Thus the linearization of Equations A22-A24 will involve linearizations of u, e, and p. Furthermore, we can make use of the identities { + } = { } + { } since we are using simple averages for the numerical flux. For example, the linearization of the mass conservation flux with respect to reference states ρ r and (ρu) r is calculated by including infinitesimal perturbations ρ and ρu, for example, where in the last line we made use of We condense Equation A29 by defining the reference velocity u r and linearized velocity u L as so that Similarly we define linearized and reference values as In total, the Jacobian of Equations A22-A24 with respect to a reference state ρ r , (ρu) r , ρe r , yields the linearized numerical fluxes We see by inspection that the above system is indeed linear with respect to ρ, ρu, and ρe.
For the interface term component of the numerical flux, we use linearized versions of the interface flux used in the full equations plus a reference state based Rusanov flux for the penalty term. Each column has its own reference state and the resulting linear systems are factored and solved directly. The reference state itself is constructed from instantaneous values of density, horizontal-momentum, and total-energy. Projecting out the vertical momentum from the reference state makes the method slightly more robust.

A3. Kinetic + Potential Energy Preservation
As mentioned in Section 2.1, our choice of numerical flux is mimetic of The calculation is as follows. We start by noticing where the flux  includes the non-conservative gravity term. We now replace the above equation with the discrete analog, using a natural extension of the notation in Section 2, The ∇ i′ operator is acting in the same capacity as the 2D ii′ operator from Section 2, but without explicitly writing out all the summations and different component directions. Hence, it is no longer a differential operator but a matrix that implicitly sums over the i′ index. Choosing the numerical flux from Equations 40-42 yields Note the following identities to write Equation A44 as where the first term is the divergence of a conservative flux and the second term is the non-conservative pressure work. Since the left term has been expressed as a flux, the equation is now mimetic of the continuous kinetic + potential energy equations, With regards to checking that one has a conservative flux, it is not necessary to express the terms as averages, rather, one can simply see if a resulting expression is invariant with respect to interchanging i and i′.
As our final comment, we state that kinetic energy preservation amounts to noticing that is expressible as a conservative flux, where ma is the advection part of the momentum flux. This is always possible, as noted by (G. J. Gassner et al., 2016a), if we choose the advection part of the momentum flux as as we have seen from the prior calculations. The analogous statement for the potential energy term is where Φ is the non-conservative flux associated with gravity. For example, a central flux for  , implies the use of the usual gravity source term, for the non-conservative flux, so that in the kinetic + potential energy equations.

Data Availability Statement
The software to plot all the figures is found in the GitHub repository (https://github.com/sandreza/DryAtmos-phereFluxDifferencingVisualization) archived at Zenodo (A. Souza, 2022a). The data files are found via figshare at (A. Souza, 2022b) along with the software used to produce the data.