Hydrostatic and Non‐Hydrostatic Baroclinic Instability in the Dynamical Core of the DOE Global Climate Model

The dynamical core of the Department of Energy global climate model is used to understand the role of non‐hydrostatic dynamics in the simulation of dry and moist baroclinic waves. To reduce computational cost, the Diabatic Acceleration and REscaling approach is adopted. Scale analysis and numerical simulations suggest that the model solution is not distorted by the change of spherical metric terms due to the change of Earth radius, neither are the inter‐scale interactions strongly altered as indicated by the analysis of spectral flux. Compared with hydrostatic simulations, the onset of baroclinic instability is delayed under non‐hydrostatic dynamics as the associated weaker vertical motions tend to increase the critical wavelength and narrow the range of waves that can be baroclinically unstable. During the development of baroclinic waves, non‐hydrostatic dynamics tends to induce vertical motions in the upper troposphere, accelerating the eastward propagation of upper‐level ridges through their impact on local vorticity tendency. These processes reduce the westward tilt of the vertical ridge axes and suppress the conversion of mean flow available potential energy to eddy kinetic energy, leading to weaker baroclinic eddies than in hydrostatic simulations. The contrasts between hydrostatic and non‐hydrostatic settings hold in supplemental experiments that vary the background flow Rossby number, the amount of water vapor content and vertical resolution. We also find that mesoscale and smaller‐scale activities can be considerably under‐represented when the vertical resolution is limited to that typically used in global climate models.

goes to infinity, its vertical scale is characterized with ∼ (Billant & Chomaz, 2001;Lindborg, 2005;Waite & Bartello, 2006). In the literature, the latter scale is sometimes termed as the buoyancy scale and indicates the thickness of the vertical shear layers in stratified turbulence. At even smaller scales, the outer scale of isotropic three-dimensional turbulence is described by the Ozmidov length scale = ( 3 ) 1 2 , where is the dissipation rate of turbulent kinetic energy (KE). In the free troposphere is typically ∼10 m (Luce et al., 2019;Waite, 2011), much larger than the Kolmogorov microscale where energy is dissipated by viscosity. Although physical processes at very small scales are usually parameterized in atmospheric models, inadequate vertical resolutions that under-resolve the relatively large scales tend to produce spurious flow structures and grid-scale noises as reported in earlier investigations (e.g., Lindzen & Fox-Rabinovitz, 1989;Pecnick & Keyser, 1989;Skamarock et al., 2019;Snyder et al., 1993;Waite, 2016). These studies also suggested that the vertical resolution should be further constrained in the presence of fronts and gravity wave critical layers. Moreover, the vertical resolution needs to be consistent with horizontal resolution for faithful simulations. For instance, the use of grids that resolve in the vertical but not in the horizontal can filter out small-scale near-isotropic motions, leading to incorrect representations of Kelvin-Helmholtz instabilities, internal gravity waves, saturating vortex instabili ties, and the subsequent transition to isotropic three-dimensional turbulence (Waite, 2011). As demonstrated by Waite (2011), in the simulation of stratified turbulence an anisotropic grid with Δ Δ (Δ and Δ are, respectively, the horizontal and vertical grid spacings) is inadequate to properly represent non-local interactions, and hence cannot capture the associated energy transfer from the energy-containing scales to the buoyancy scale.
In the atmosphere, water vapor plays an important role in dynamics and thermodynamics through its greenhouse effect and its phase change induced latent heating or cooling. Previous studies have shown that moist processes can lead to additional growth at the synoptic scale (e.g., Snyder & Lindzen, 1991;Whitaker & Davis, 1994), accelerate frontogenesis and organize convection at the mesoscale (e.g., Houze, 2004;Thorpe & Emanuel, 1985), and enhance extreme precipitation at the convective scale (e.g., Loriaux et al., 2013). According to Hamilton et al. (2008) and Waite and Snyder (2013), the energy cascade is also altered when moisture is included in numerical simulations, where the KE spectrum is featured with a larger amplitude than in dry runs as moist processes energize the mesoscale. A question raised by these effects is whether and how the non-hydrostatic dynamics is influenced by moist physics.
As new computing paradigms appear, high-resolution global climate simulations are getting affordable, and it is essential to improve our understanding of non-hydrostatic dynamics. In this study, we examine in detail the non-hydrostatic effect in baroclinic waves, a prototype of midlatitude disturbances. During the life cycle of baroclinic eddies, baroclinic instability drives frontogenesis (e.g., Nadiga, 2014) and acts together with convection to stabilize the atmosphere, while inertia gravity waves are spontaneously generated (e.g., Waite & Snyder, 2013), especially at the surface front (e.g., Snyder et al., 1993) and near the upper-level jet (e.g., Plougonven & Snyder, 2005). Specifically, we examine the resolution at which the non-hydrostatic effect becomes important in dry and moist baroclinic wave simulations, and how this resolution is impacted by the background flow properties, by the amount of moisture in the atmosphere and by the number of vertical levels. In Section 2, the numerical models employed are briefly introduced, followed with a description of simulations in Section 3, where the Earth radius is reduced to explore the effect of non-hydrostatic dynamics at reduced computational cost. When moisture is included in the small-earth experiments, diabatic processes also need to be accelerated as proposed by Kuang et al. (2005). We document relevant technical details in Appendix A. Relying on scale analysis, we focus on non-hydrostatic effects under small-earth setups in Section 4, and report model results in Section 5. Conclusions are given in Section 6. SUN ET AL.

Model Description
The simulations are performed with the dynamical core of the Department of Energy′s (DOE) Energy Exscale Earth System Model (E3SM; Golaz et al., 2019), including the hydrostatic High-Order Methods Modeling Environment (HOMME; Dennis et al., 2012) and its non-hydrostatic enhancement HOMME-NH . In the horizontal directions, HOMME and HOMME-NH are based on continuous Galerkin spectral elements on a cubed-sphere grid and use degree three polynomials to achieve a fourth-order accuracy. In the vertical, both dynamical cores are built upon the terrain-following hybrid sigma-pressure coordinate and use second-order finite differences. Since it is found that collocated vertical grids may introduce spurious waves associated with the top boundary conditions (not shown), the Lorenz staggering is applied. For timestepping, a multi-stage second-order accurate Runge-Kutta method is implemented in HOMME and improved as a multi-stage implicit-explicit Runge-Kutta algorithm in HOMME-NH. As thoroughly evaluated in Dennis et al. (2012), Taylor et al. (2020), and Guba et al. (2020), both HOMME and HOMME-NH can produce well-converged solutions when running the model without forcing or explicit dissipation.
Both HOMME′s hydrostatic and non-hydrostatic formulations use a Hamiltonian structure preserving discretization of the adiabatic component of the dynamics. As such, the adiabatic core of the model conserves total energy in exact timestepping . In the full dynamical core there are three sources of dissipation. The largest source of dissipation is the explicitly added ∇ 4 hyperviscosity (Guba et al., 2014). In addition, when running with E3SM′s full suite of physical parameterizations, the dynamical core adds a sponge layer at the model top. This sponge layer is not needed and turned off for the simulations presented here. The remaining two sources of dissipation in the model come from the Runge-Kutta timestepping, and from the monotone limiter used when remapping the solution on the vertically Lagrangian layers back to the reference layers.
Compared with HOMME, additional prognostic equations of vertical velocity and geopotential (Taylor et al., 2020, their Equations 10 and 11) are introduced in HOMME-NH to account for non-hydrostatic dynamics. Consequently, the budget equation for KE differs substantially in HOMME and HOMME-NH. Taking the dry case as an example and neglecting dissipation and the forcing terms for clarity, the KE equation in HOMME can be written as where is time, is the hybrid sigma-pressure vertical coordinate such that ⋅ = and ∇ is the two-dimensional gradient operator along the surfaces, is the so-called pseudodensity with as the hydrostatic pressure, = ( ) denotes horizontal velocity, is the specific heat of air at constant pressure, is potential temperature, and Π is the Exner pressure. In HOMME-NH, the corresponding KE equation is where is gravitational acceleration, and = with p as the total pressure. Because = 1 2 2 in Equation 1 for HOMME while = 1 2 ( 2 + 2 ) in Equation 2 for HOMME-NH, as well as the fact that is absent on the right-hand side of Equation 1, HOMME only accounts for the contribution of horizontal motions to and omits the contribution from the vertical component. As pointed out by Pauluis et al. (2006), although Equation 1 is valid for large-scale motions with small vertical velocity, it is problematic for convective motions where the vertical velocity can be large.

Simulations
The initial condition of our baroclinic wave simulations is derived in Ullrich et al. (2014) and implemented as in Ullrich et al. (2016, henceforth U16). Specifically, the base state virtual temperature, pressure and zonal wind fields are specified per Equations 5-13 of U16, while the reference meridional wind is set as zero. After 4 of 26 initializing the model with zero surface geopotential and constant surface pressure, the reference state is under geostrophic and hydrostatic balances but baroclinically unstable (Figure 1a). To trigger baroclinic instability, a Gaussian-shaped perturbation centered at (40°N, 20°E) is applied to the zonal wind field using Equations 14-16 of U16. When moisture is considered, specific humidity is prescribed as Equations 17 and 18 of U16, yielding a maximum relative humidity of ∼85% in the lower troposphere of the midlatitudes (Figure 1b). To keep the reference state in geostrophic balance in the presence of moisture, the base state temperature is adjusted and colder than the temperature one would obtain in a dry atmosphere. However, given that the pressure gradient term in the momentum equations is determined by virtual temperature and moist pressure, which are identical to the temperature and pressure in the dry case, the dry and moist configurations would produce identical solutions in the absence of diabatic processes.
As an adjunct to investigating non-hydrostatic effects with increased effective horizontal resolution, we consider the Diabatic Acceleration and REscaling (DARE) approach (Kuang et al., 2005). To this end, we conduct baroclinic wave simulations on a small and rapidly rotating Earth with enhanced diabatic forcing. After reducing the Earth radius as = , where X is the scaling or reduction factor, the following modifications are made: (a) the dynamics timestep is divided by X; (b) the Earth rotation rate Ω is multiplied by X; and (c) the diffusion coefficient is divided by 2 2 −1 , where 2 = 2 corresponds to second-order hyperviscosity. For the moist experiments involved in this study, the diabatic forcing is limited to the Kessler microphysics (Kessler, 1969) and scaled by a factor of X as described in Appendix A to correctly represent its interaction with the scaled large-scale flow. In DARE, the scale of convective motions remains unaffected, while the horizontal spatial scale of large-scale circulations is reduced by a factor of X due to the reduction of Rossby radius. Consequently, the timescale and horizontal length scale are reduced by a factor of X, while the vertical velocity is enhanced X times since the vertical length scale remains unchanged. With → and → , the vertical momentum equation is scaled as where ρ is air density, z is model height and B is buoyancy, using overbar to indicate hydrostatic balance and prime to denote deviations. Equation 3 exactly resembles the vertical momentum equation in hypohydrostatic   (Garner et al., 2007, their Equation 9; Pauluis et al., 2006, their Equation 8), and these two approaches are in fact mathematically equivalent for inviscid adiabatic flows in a dry atmosphere (Garner et al., 2007;Pauluis et al., 2006;Semane & Bechtold, 2015;Wedi & Smolarkiewicz, 2009). However, it should be noted that, when the same number of grid cells are used in the model for different values of X, the simulations relying on DARE share the same grid and are computationally characterized with the same grid spacing as under real-earth setups (Δ ). The differences of model results from simulations with different X values stem from non-hydrostatic dynamics due to the change of effective aspect ratio Boos et al., 2016;Garner et al., 2007;Hsieh et al., 2020;Peters & Bretherton, 2006).
Our simulations are integrated for 100 days and conducted at the ne30 resolution, where "ne" denotes the number of elements and ne30 corresponds to a resolution of approximately one degree or ∼111 km horizontal grid spacing in real-earth setups with = 6, 376 km . For small-earth simulations, we consider values of the scaling factor X of 10, 50, 100, 200, 300, 400, and 600; the resultant effective grid spacings are about 11.1, 2.2, 1.1, 0.6, 0.4, 0.3, and 0.2 km, respectively. To examine the role of vertical resolution, the number of vertical levels is varied (with 30, 60, 120, 180, and 300 levels), with a model top at ∼35 km. The maximum vertical grid spacing is up to ∼3.5 km in the stratosphere when 30 vertical levels are used and reduced to ∼0.4 km in the simulations with 300 vertical levels. In the upper troposphere and lower stratosphere, the vertical levels are rather uniformly distributed, characterized with a typical vertical grid spacing of ∼1.1 km in the 30-level simulations and ∼0.11 km in the 300-level runs. The convention used to name the experiments consists of an identifier for the presence of moisture (D for dry and M for moist), followed by an indicator of dynamics (H for hydrostatic and NH for non-hydrostatic), the scaling factor X, and the number of vertical levels. For instance, a dry non-hydrostatic experiment with a scaled planet radius 100 times smaller than the Earth radius using 30 vertical levels is labeled as D_NH100X30Z, while its moist variant would read as M_NH100X30Z.
To explore the sensitivity to flow properties, we conducted another set of experiments where the Rossby number of the background flow is increased (Figure 1c). This is accomplished by increasing the surface temperature contrast between the equator ( ) and the poles ( ) in Equations 5-13 of U16. In this way, the strength of the background zonal flow and its Rossby number are increased through thermal wind adjustment, while the reference state is still in geostrophic and hydrostatic balances. In the model, the default values of and are 310 and 240 K, respectively. They are modified as 315, 320, and 325 K for and 235, 230, and 225 K for , so that Δ = − is changed from the default 70 K to 80, 90, and 100 K accordingly. These experiments are named with one more identifier for Δ . For instance, D_NH100X30Z_80K and M_NH100X30Z_80K denote the same simulations as D_NH100X30Z and M_NH100X30Z, respectively, but with Δ = 80 K . As in Waite and Snyder (2013), we also examined the sensitivity to atmospheric moisture content. Specifically, in our moist runs the maximum specific humidity q 0 in Equation 18 of U16 is changed from the default 18 to 10 g kg −1 , which yields a maximum low-level relative humidity of ∼44% in the midlatitudes (Figure 1d). These modifications are reflected when labeling relevant simulations, for example, M_NH100X30Z_Q10 is the same as M_NH100X30Z, but with 0 = 10 gkg −1 instead of the default 0 = 18 gkg −1 .

Scale Analysis of Small-Earth Simulations
As both HOMME and HOMME-NH use the shallow-atmosphere approximation, our scale analysis is based on the shallow-atmosphere momentum equations that may be written as where = cos and =

of 26
In these equations, φ is latitude, λ is longitude and = 2Ω sin is the Coriolis parameter, using F rx , F ry , and F rz to denote the frictional forces in the x, y, and z directions, respectively. To gauge the scale of each term more accurately, the vertical momentum equation is written in its perturbed form (Equation 6), where 0 = 0( ) is the reference air density that is in hydrostatic balance with the reference pressure 0 = 0( ) , while the prime indicates deviations from the reference state. Compared with Holton (2004, his Equations 2.19-2.21; henceforth H04), the metric terms not involving tanφ and the Coriolis terms that contain a factor of cosφ are dropped to make the conservation properties of energy, angular momentum, and potential vorticity analogous to the non-approximated deep-atmosphere equations (see White et al., 2005 for a thorough discussion). In HOMME and HOMME-NH, this treatment is needed to preserve the Hamiltonian structure.
In small-earth simulations with a large reduction factor, the shallow-atmosphere assumption breaks down, as ≫ is not valid anymore. This effect, however, has limited impact on the simulation of synoptic systems. Taking X = 1,000 for an instance and using the characteristic scales in H04 (p. 39), except for friction, the scale of the rest of the terms in Equations 4 and 5 increases X = 1,000 times as documented in Table 1. Thus, as in the real-earth setups, the dominant balance in the horizontal is still between the Coriolis terms (2Ω sin and − 2Ω sin ) and the pressure gradient force, and the retained metric terms ( tan and − 2 tan ) remain at least one order of magnitude smaller than other terms, except for the trivial F rx and F ry . In the vertical, is scaled up by a factor of X 2 = 10 6 , making it comparable with the terms involved in hydrostatic balance, including − 1 0 ′ and − ′ 0 . These results imply that, as we shall see in Section 5, the large-scale baroclinic waves in small-earth experiments would not be distorted by the change of spherical metric terms even after dramatically reducing the Earth radius 1,000 times.
For processes at smaller scales, where the vertical acceleration is no longer negligible, their flow properties can differ substantially from the large-scale environment. When the large-scale reference state is in hydrostatic balance and characterized with zero vertical velocity, the perturbation form of the continuity equation may be written as where = ( ) and the term 1 is neglected as it is typically one order of magnitude smaller than the retained terms. Since and in the large-scale environment tend to be of equal magnitude but with opposite sign, where is the vertical velocity scale, imposing an upper bound for ′ .
Alternatively, if we take the large-scale hydrostatic reference sate as in Lin (2007, p. 14), such that ) is the large-scale geostrophic wind, the overbar indicates the basic state, and the general quantity refers to ρ, p, θ, and T, the continuity equation becomes ( is the speed of sound) and (J is diabatic heating rate) are neglected due to their small magnitude. We therefore get the constraint between ′ and ′ as which, in combination with Equation 8, suggest the following scenarios for shallow and deep convection at small scales (e.g., = 10 4 m , the typical size of cumulus convection): In Equation 12, the upper bound of ′ in Equation 8 is taken and the large-scale = 10 ms −1 is used as in H04 (p. 39).
Using the same large-scale reference state (Equation 9) as in Lin (2007), we next write Equations 4-6 in their perturbed form as where the f-plane approximation is applied, not only for simplicity but also because it is a fairly good approximation for small-scale processes. Based on Equation 12 and the relationship ′ ∼ 2 = 100 m 2 s −2 ( ′ indicates horizontal pressure fluctuations) (Markowski & Richardson, 2010, their Equation 1.18), the characteristic magnitude of each term in Equations 13-15 is presented in Table 2. As in the real-earth setups, the metric terms ( tan and − 2 tan ) in Equations 13 and 14 are still at least two orders of magnitude smaller than other terms for both shallow and deep convection in small-earth experiments (X = 1,000), if the very small F rx and F ry , which can be occasionally large near the edge of clouds or in rising thermals (Markowski & Richardson, 2010), are not considered. In the vertical (Equation 15), shallow convection is switched from hydrostatic dynamics

SC
Real Earth 10 −2 10 −2 10 −2 10 −2 10 −5 10 −3 10 −2 10 −10 Small Earth (X = 1,000) 10 10 10 10 10 −2 1 10 10 −10 DC Real Earth 10 −2 10 −2 10 −2 10 −2 10 −5 10 −3 10 −2 10 −12 Small Earth (X = 1,000) 10 10 10 10 10 −2 1 10 10 −12 Real Earth 10 −3 10 −3 10 −3 10 −2 10 −2 10 −13 Small Earth (X = 1,000) 10 3 10 3 10 3 10 −2 10 −2 10 −10 DC Real Earth 10 −2 10 −2 10 −2 10 −2 10 −2 10 −15 Small Earth (X = 1,000) 10 4 10 4 10 4 10 −2 10 −2 10 −12   Table 2. These results demonstrate that, in a model with shallow-atmosphere approximation, the radius-reduction induced change of small-scale processes is not due to the change of spherical metric terms even after reducing the Earth radius 1,000 times. Moving from the real-earth setup to small-earth configurations, shallow convection may change more than deep convection. Admittedly, the inter-scale interactions may be affected when X is very large due to the change of scale separation between the large-scale and convective motions. The spectral budget analysis reported later in Section 5, however, indicates that this is not the case for the idealized simulations presented here. For small-earth runs with the same scaling factor X, non-hydrostatic dynamics is the fundamental contributor to the differences between the hydrostatic HOMME and non-hydrostatic HOMME-NH simulations. Figure 2 shows the evolution of eddy kinetic energy (EKE) for the latitude band from 30° to 60°N based on the simulations using 120 vertical levels. This quantity is calculated as a mass-weighted vertical average using Equation 1 of Waite and Snyder (2013), but only retains the horizontal component as in Simmons and Hoskins (1978). In the dry and hydrostatic runs (Figure 2a), baroclinic instability starts developing around day 8, saturates around day 25, and decays thereafter as no external forcings exist to maintain the eddies (see Figure S1 and relevant discussions in Supporting Information S1). Compared with other idealized studies (e.g., Gutowski et al., 1992;Hsieh et al., 2020;Waite & Snyder, 2013), the simulated baroclinic waves are featured with a relatively large amplitude. The fact that the baroclinic eddies are indistinguishable in the dry hydrostatic experiments with varied reduction factor X confirms that the simulations are not distorted by the change of spherical metric terms under shallow-atmosphere approximation, in agreement with the scale analysis in Section 4. This result also follows the nature of DARE, which only modifies the effective aspect ratio to accentuate non-hydrostatic effects, rather than changing the computational grid spacing (Section 3). When non-hydrostatic dynamics is not considered, properly implemented dry hydrostatic DARE experiments would produce very similar, if not exactly the same, results under different values of X as shown in Figure 2a. In the dry non-hydrostatic experiments (Figure 2b), the onset of the eddies is delayed. This delay is up to several days and more prominent as the reduction factor X increases. In the experiment with X = 600, the eddies are also weaker under non-hydrostatic dynamics compared to the hydrostatic case (cf., the magenta line in Figures 2a  and 2b). When moisture is included (Figures 2c and 2d), both the hydrostatic and non-hydrostatic setups produce stronger eddies than in dry simulations. This is consistent with previous studies (e.g., Booth et al., 2013;Gutowski et al., 1992;Waite & Snyder, 2013), and the maximum EKE in some moist experiments nearly doubles that in the corresponding dry experiment (cf., D_H010X120Z in Figure 2b and M_H010X120Z in Figure 2d). Because of the physical constraints on moisture variables (Equations A14-A16, Appendix A), EKE from moist hydrostatic simulations varies with the reduction factor X as shown in Figure 2c, in contrast to Figure 2a where the DARE scaling produces indistinguishable EKE for dry hydrostatic runs. As in the dry non-hydrostatic experiments, the onset of baroclinic instability is also delayed in the moist non-hydrostatic simulations, with the weakest storm produced in the M_NH600X120Z experiment (the magenta line in Figure 2d).

Results
Since the overall behavior of the baroclinic waves has limited sensitivity to the number of vertical levels as discussed later and because non-hydrostatic effects are most obvious when the earth-radius reduction factor is large (Figure 2), we analyze the non-hydrostatic effect in Section 5.1 using the 120-level X = 600 runs. These Evolution of eddy kinetic energy (EKE) (m 2 s −2 ) for the latitude band from 30° to 60°N. They are mass-weighted vertical averages from the simulations as indicated in the legend. As Diabatic Acceleration and REscaling is exact for the dry hydrostatic simulations, their EKE overlaps in panel (a). This is not the case for the moist hydrostatic simulations in panel (c) due to the required physical constraints on moisture variables as detailed in Appendix A. 9 of 26 simulations, including D_H600X120Z, D_NH600X120Z, M_H600X120Z, and M_NH600X120Z, have an effective grid spacing ∼0.2 km that is close to the resolution at which the cloud area in squall line simulations tends to converge (Lebo & Morrison, 2015). As suggested by the probability density functions in Figure 3, in these simulations vertical motions are more intense under hydrostatic approximation than using the non-hydrostatic dynamics (Figure 3a), and the moist hydrostatic runs are consistently featured with enhanced precipitation intensity, rainwater and cloud content (Figures 3b-3d). These results agree well with previous studies focusing on the non-hydrostatic effect (Kato & Saito, 1995;Orlanski, 1981;Sun et al., 2023;Weisman et al., 1997;Yang et al., 2017), even though our simulations are limited by under-resolved convection due to the coarse ne30 mesh and constrained by the Kessler warm rain microphysics. The sensitivity to Rossby number of the background flow, the moisture content in the atmosphere and vertical resolution is reported in Section 5.2.

Non-Hydrostatic Effects in Baroclinic Waves
The delayed onset of baroclinic eddies is further illustrated in Figure 4, showing the fluctuations of surface pressure from its zonal mean. Besides in the Northern Hemisphere on day 10 ( Figure 4, left), the delay in the Southern Hemisphere (Figure 4, right) is also evident: on day 25, baroclinic waves are prominent in the hydrostatic D_H600X120Z and M_H600X120Z experiments along the vicinity of 45°S but still absent in the non-hydrostatic D_NH600X120Z and M_NH600X120Z simulations. In a supplemental experiment, we removed the initial perturbations that were used to trigger baroclinic waves, but only relied on numerical truncation errors to excite the eddies. The development of baroclinic waves, including in the Southern Hemisphere, was delayed compared to the simulations reported in Figure 4 (not shown), implying that the Southern Hemisphere eddies are more likely associated with the propagation of disturbances from the Northern Hemisphere to the Southern Hemisphere, rather than stemming from numerical truncation errors.
The contrast in Figure 4 may be related to the weaker vertical motions in the non-hydrostatic simulations than in the hydrostatic runs, as indicated by the root-mean-square (RMS) vertical pressure velocity shown in Figure 5. Following H04, we qualitatively interpret the delayed onset of baroclinic instability using an idealized geopotential distribution on a midlatitude β-plane written as which consists of the sum of a zonally averaged part ϕ 0 − f 0 Uy and a zonally varying part 0 sin cos that depends sinusoidally on x and y. In Equation 16, ϕ 0 is the latitude at which the Coriolis parameter f 0 is evaluated, k x and k y are, respectively, the zonal and meridional wavenumbers, and the parameters ϕ 0 , U and A only depend on pressure. A distinct advantage of this simplification is that the geostrophic vorticity can then be defined as = − ( 2 + 2 ) sin cos . Proceed to obtain the omega equation as in H04 (his Section 6.4.1), it can be shown that the required vertical pressure velocity ω to maintain a hydrostatic temperature field, where temperature and thickness are proportional, follows the relationship Recall that is proportional to the square of the horizontal wavenumber, Equation 17 suggests that the weaker vertical velocity in the non-hydrostatic simulations corresponds to a smaller horizontal wavenumber or a longer wavelength than in the hydrostatic runs characterized with more intense vertical motions. In the case of this wavenumber corresponds to the critical wavelength, below which the flow is always stable, a narrower range of wavelengths are unstable under non-hydrostatic dynamics, leading to delayed development of baroclinic waves than using the hydrostatic approximation. This analysis is also consistent with the stronger vertical ageostrophic flow in the hydrostatic experiments as indicated by the horizontal divergence field (not shown). As a secondary circulation that arises in response to the disruption of thermal wind balance by the primary geostrophic motion, the ageostrophic flow acts to restore this base state balance. Intuitively, the diagnosed vertical ageostrophic circulation is expected to be strong under hydrostatic approximation where the hydrostatic condition is instantaneously and rigidly enforced. Figure 6 shows the KE spectra computed as in Jakob et al. (1993, their Equation B.15). They are calculated at each 6-hourly time interval using a spherical harmonic transform of the instantaneous horizontal velocity field and averaged for days 5-15, 20-30, and 90-100. Due to the limitation of the ne30 mesh, the spherical wavenumber k stops at 120, corresponding to a resolvable resolution that is three times of the computational grid spacing. Rather than resembling the −5/3 mesoscale scaling (Nastrom et al., 1984), the slope of the total KE spectra ( Figure 6, solid lines) follows the synoptic k −3 relationship at relatively small wavenumbers even though the Earth radius has been reduced by a factor of X = 600. This is because the Earth rotation rate is also increased X = 600 times, such that the Rossby number remains unchanged compared with real-earth setups. Consequently, we still term k ≤ 8 as the planetary scale, 9 ≤ k ≤ 40 as the synoptic scale, and k > 40 as the mesoscale and smaller scales in the small-earth simulations. Under real-earth setups, this is equivalent to defining the planetary scale as for wavelengths greater than ∼5,000 km, the synoptic scale for the wavelengths between ∼4,500 and ∼1,000 km, with the rest shorter wavelengths designated as the mesoscale and smaller scales. The −3 slope is more evident in the lower troposphere than at upper levels and more evident later into the simulation (cf., column 1 and 2 of Figure 6), probably related to the limited development of baroclinic waves during days 5-10. During days 90-100, air motion in the model is largely random and turbulent; the spectra are to some extent noisy and show a sign of aliasing for wavenumbers greater than ∼100 (Figure 6, column 3). The lack of the transition from the large-scale −3 spectrum to a mesoscale −5/3 spectrum is likely due to the dominance of rotational KE over divergent KE (Hamilton et al., 2008;Lindborg, 2007;Waite & Snyder, 2009. Specifically, the total KE spectra are dominated by the balanced rotational contribution (Figure 6, dotted lines) up to k ∼ 100, whereafter divergent inertia-gravity waves ( Figure 6, dashed lines) contribute roughly equally as the rotational component to the total KE spectra. Even though the divergent KE spectra are characterized with a shallow slope at large wavenumbers, their amplitude is too small to affect the slope of the total KE spectra. Alternatively, the lack of mesoscale shallowing may be related to the scale at which the transition toward the −5/3 slope occurs. In our simulations the resolvable scale at k = 120 corresponds to a wavelength of ∼333 km under real-earth setups, while the transition scale can be at a smaller scale of ∼160 km as found in Waite (2016). Thus, the shallowing toward a −5/3 spectrum can be absent in Figure 6 because the transition scale has not yet been reached. The steeper slope than −3 at wavenumbers greater than ∼70 is also featured in the real-earth KE spectra ( Figure S2 in Supporting Information S1) and may be related to the parameterization of gravity wave drag and numerical filters (Skamarock et al., 2019). As in Hamilton et al. (2008) and Waite and Snyder (2013), the spectral energy also increases in the presence of moist processes, especially at relatively large wavenumbers. Compared with the hydrostatic runs, more energy is distributed at relatively large wavenumbers in the non-hydrostatic experiments during the last 10 days of model integration when mesoscale and smaller-scale motions are prevalent over the globe (Figure 6, column 3).
To deepen our understanding of the physics underlying the energy cascade in Figure 6, the conversion of available potential energy (APE) to KE, as well as the spectral fluxes of APE and KE, are analyzed in Figure 7. They are calculated as in Augier and Lindborg (2013, their Equations 17, 19 and 23; henceforth AL13) and denoted as , Π , and Π , respectively. To facilitate the interpretation, these quantities are vertically integrated (mass-weighted) and summed up over total wavenumbers (i.e., the spherical wavenumber k) (AL13, their Equation 27). Taking ( ) as an example, it therefore represents the cumulative conversion of APE to KE from all wavenumbers greater than or equal to k. As illustrated in Malardel and Wedi (2016, their Figure 3; henceforth MW16), it follows that: (a) Π ( ) > 0 represents a downscale energy transfer, with the energy in wavenumbers [0, ) being deposited to wavenumbers [ max] ( max is the largest wavenumber considered in the calculation), and vice versa; (b) Π ( ) > 0 ( Π ( ) < 0 ) corresponds to an energy loss (gain) at wavenumber , while a plateau with Π ( ) = 0 indicates a neutral transfer without energy gain or loss at wavenumber k. Given Π ( ) denotes the energy transfer accumulated over all the wavenumbers that are greater than or equal to k, and because the non-cumulative spectral flux is undefined at = 0 (AL13, their Equation 12), Π (1) is theoretically zero as spectral transfer only redistributes energy among wavenumbers. Due to discretization, however, the practically calculated Π (1) can deviate slightly from zero as seen in Figure 7. At small wavenumbers (k ≤ 8), the shape of the curves in Figure 7 (solid green lines) differs from similar analyses in previous studies (AL13, their Figure 1; MW16, their Figure 4). In those studies, the spectral budget analysis is based on comprehensive simulations, and the conversion between APE and KE at the planetary scale is considered resulting from the Hadley and Ferrel circulations (Lambert, 1984). Both circulations, however, are absent in our idealized experiments. In Figure 7, the conversion of APE to KE at the planetary scale (k ≤ 8), as well as at the synoptic scale (9 ≤ k ≤ 40), is thus due to the development of baroclinic instability. Along with the development of the baroclinic eddies, APE is transferred downscale from wavenumbers k ≤ 3 toward subsequent larger wavenumbers (Π in Figure 7, solid blue lines), and the APE spectral flux dominates over the KE spectral flux that mainly transfers energy upscale as indicated by its overall negative values (Π in Figure 7, solid red lines). This upscale KE transfer agrees with the fact that Π is dominated by its rotational component (Π in Figure 7, dotted red lines) because interactions of rotational modes act to transfer energy upscale and feed the large-scale zonal flow (AL13). However, it should be noted that Π is slightly positive for wavenumbers between ∼20 and 30 and ∼110 (Figure 7, the inset plots), transferring KE down the cascade and depositing it at the higher end of this wavenumber range. According to AL13 and MW16, this positive Π is related to the mesoscale −5/3 slope in the KE spectrum, but it is fairly weak in our simulations (up to ∼0.3 W m −2 ), consistent with the absence of the mesoscale shallowing in Figure 6 and Figure S2 in Supporting Information S1. In the moist simulations, the total spectral flux of APE and KE (Π + Π in Figure 7, solid orange lines) extends further into the small scales than in their dry counterparts, consistent with the greater spectral energy at relatively large wavenumbers in the moist simulations ( Figure 6). Compared with Figure S3 in Supporting Information S1, the , Π , and Π curves from the X = 600 and real-earth hydrostatic runs are almost identical. These results demonstrate that the inter-scale interactions, in terms of energy conversion and the transfer of energy flux across wavenumbers, are not distorted in the small-earth experiments, as otherwise the exact dry hydrostatic DARE scaling as indicated by the resemblance between Figure 7a and Figure S3a in Supporting Information S1 would not be possible. The very high similarity between Figure 7c and Figure S3c in Supporting Information S1 further illustrates that the impact from the moisture-induced non-exact DARE scaling is also rather limited in our idealized simulations, although this may not be the case under more complex or comprehensive setups.
Compared with the hydrostatic runs, the conversion of APE to KE for wavenumbers smaller than ∼15 is strongly suppressed in the non-hydrostatic X = 600 simulations (cf., Figures 7a and 7c with Figures 7b and 7d, solid green lines), leading to much weaker baroclinic eddies as shown in Figure 2. The accompanied downscale APE spectral flux Π is also much reduced, even though more APE is retained in the atmosphere due to decreased APE conversion in the non-hydrostatic simulations. The positive APE spectral flux Π , however, spreads more toward the large wavenumbers (up to k ≈ 90) under non-hydrostatic dynamics, and this effect is more prominent than the contrast between the dry and moist simulations. These results point to the importance of non-hydrostatic effects on APE conversion. As demonstrated in previous studies, either using a two-layer model (e.g., H04, his Equation 8.38) or for an Eady atmosphere (e.g., Lackmann, 2011, his Equation 7.86), the conversion of mean flow APE to eddy KE must be accompanied by a net poleward heat transport that requires the ridge and trough axes to tilt westward with height. In the model, the phase line of the baroclinic eddies does slope westward more in the hydrostatic simulations than in the non-hydrostatic runs, as illustrated by the vertical cross sections of zonal geopotential height anomaly along the 45°N latitude on day 20 (Figure 8, contours). Similar contrast between the hydrostatic and non-hydrostatic simulations is also evident most of the time from days 10-30 (Movie S1, contours).
We may understand the change of vertical phase line tilt shown in Figure 8 and Movie S1 from the perspective of vorticity tendency, which is closely related to the vertical acceleration in the non-hydrostatic runs. In these simulations, the vertical momentum equation can be written as Figure 7. Nonlinear spectral fluxes (W m −2 ) of available potential energy (Π , solid blue) and kinetic energy (KE) (Π , solid red), as well as the conversion of available potential energy to KE ( , solid green) from the simulations as indicated in the panel titles. Panels (a and c) are from hydrostatic simulations, while panels (b and d) are from non-hydrostatic runs. The total spectral flux Π + Π is shown as orange lines. Π (dottedred) is the rotational component of Π . They are vertically integrated (mass-weighted) and accumulated over spherical wavenumbers. The calculation is conducted at every 6-hourly time interval and averaged for the period during days 20-30. The inset is to illustrate the small magnitude of Π and Π for wavenumbers between ∼20-30 and ∼110.
where the vertical dissipation term is omitted due to its small magnitude. If we define the non-hydrostatic parameter = 1 = − 1 as in Janjic et al. (2001, their Equation 2.8), it stands for the Lagrangian divided by the gravitational constant. Following the air motion, positive ε values thus correspond to vertical acceleration and negative ε values vertical deceleration. Since = as well, we take percentage as the unit of ε to illustrate how much the non-hydrostatic pressure deviates from the hydrostatic value. According to Figure 8 and Figure S4 in Supporting Information S1 (shaded) and as illustrated spatially in Figures S5 and S6 in Supporting Information S1, in the non-hydrostatic runs non-zero ε values are mainly found in the upper troposphere, and generally indicate Lagrangian deceleration in ridges where = 1 < 0 and Lagrangian acceleration in troughs where = 1 > 0 . In typical baroclinic waves, because the ageostrophic vertical velocity is in phase with the temperature perturbation and lags the geopotential field (see Figure 8.4 of H04 for a schematic description), downward (upward) motion is prevalent on the leading (trailing) edge of ridges, while the opposite occurs in troughs. As displayed in Figure S7 in Supporting Information S1 and Movie S2, this idealized structure is not only in the Eady waves (e.g., Figure 8.10 of H04) but also in our simulations. Since both the downward and upward motions ahead of and behind the ridges would be damped under the influence of negative , the non-hydrostatic Lagrangian deceleration within ridges would induce incremental upward motions on the leading edges and incremental downward motions on the trailing edges, as depicted by the solid vertical green arrows in Figure 9. In troughs We next demonstrate how these non-hydrostatic vertical motions would influence the evolution of baroclinic waves in the context of vorticity budget. Following H04 and under β-plane approximation, the quasi-geostrophic vorticity equation can be written as where = ⋅ ∇ × refers to the geostrophic vorticity ζ g , with the subscript g dropped for clarity, and f 0 is the Coriolis parameter at a reference latitude. Using subscript t to indicate time tendency, the net vorticity tendency = is thus determined by the sum of geostrophic advection of absolute vorticity ( ) = − ⋅ ∇( + ) and the concentration or dilution of vorticity through divergence or convergence ( ) = 0 . Specifically for ( ) , vertical stretching ( > 0 ) would contribute a positive vorticity tendency and this occurs upstream of a vertical velocity anomaly (e.g., behind the open vertical black arrows in Figure 9), while the opposite occurs with vertical shrinking. As noted in Figure 9, given ridges correspond to minimum vorticity whereas troughs are characterized with maximum vorticity, ( ) tends to be negative ahead of ridges but positive ahead of troughs. On the other hand, the vorticity tendency ( ) associated with the ageostrophic circulation (open vertical black arrows in Figure 9) tends to oppose ( ) at upper levels but strengthen ( ) in the lower troposphere. The latter process is essential to maintain the westward tilt of the phase line as otherwise the upper-level ridge-trough patterns would move eastward more rapidly than the lower-level patterns because the zonal mean wind increases with height (the right panel of Figure 9). When the non-hydrostatic vertical motion is included and through vortex stretching and shrinking, vorticity at the very top of the troposphere tends to be more negative ahead of the ridges but more positive behind the ridges, as indicated by the sign of ( ) in Figure 9. Thus, non-hydrostatic dynamics contributes positively to the local net vorticity tendency and acts to accelerate the eastward propagation of the upper-level ridges as denoted by the solid eastward green arrows in Figure 9, reducing the westward tilt of the ridge lines. Relying on the same reasoning but for troughs, it follows that the eastward propagation of upper-level troughs is decelerated as indicated by the solid westward green arrows in Figure 9. It would be ideal  (negative) anomalies. The solid horizontal green arrows indicate acceleration and deceleration of upper-level ridges and troughs due to non-hydrostatic effects. The blue "−" denotes regions with negative , while the red "+" indicates regions with positive . Using subscript t to indicate time tendency, the vorticity tendencies due to advection, divergence or convergence, and non-hydrostatic dynamics are marked as (ζ t ) A , (ζ t ) D , and (ζ t ) NH , respectively. The vertical profile of the zonal mean wind is shown in the right panel. As indicated by the solid vertical green arrows, non-hydrostatic dynamics tends to induce incremental upward (downward) motions on the leading (trailing) edge of both ridges and troughs in the upper troposphere. Through vortex stretching and shrinking, these non-hydrostatic vertical motions impact the local vorticity tendency and act to accelerate the eastward propagation of upper-level ridges and decelerate the eastward propagation of upp-level troughs as denoted by the solid horizontal green arrows. These processes reduce the westward tilt of the vertical ridge axes and suppress the conversion of mean flow available potential energy to eddy kinetic energy, leading to weaker baroclinic eddies than in hydrostatic simulations.
16 of 26 to quantitatively illustrate these processes using vorticity budget analysis, but the delayed onset of baroclinic instability in the non-hydrostatic experiments results in phase differences from the hydrostatic simulations and makes such a comparison difficult. However, the schematic illustration in Figure 9 is quite well reproduced in our non-hydrostatic simulations as shown in Figure 8 and Movie S1, where in the non-hydrostatic runs the upper-level ridges are more to the east than in the lower troposphere as the model integrates forward.
In contrast, without considering the mechanical work that is needed to constrain vertical air motions, in hydrostatic simulations the mass-conservation based Lagrangian acceleration is very strong as reflected by their much larger magnitude of ε compared to the non-hydrostatic experiments (cf., Figures 8a and 8c with Figures 8b  and 8d). The relationship between zonal geopotential height anomalies and ε that is responsible for the reduced amplitude of baroclinic waves in non-hydrostatic runs is also lacking under hydrostatic conditions. Figure 10 shows the evolution of RMS ε in the X = 600 non-hydrostatic runs for the latitude band from 30° to 60°N. Because of spatial averaging the magnitude of ε is reduced, although its instantaneous value can be up to ∼7%-8%. A distinct feature in Figure 10 is that large RMS ε values are mainly found in the middle and upper troposphere, where the contrast between hydrostatic and non-hydrostatic simulations is also most prominent (Figure 8, Movie S1). This consistency corroborates our interpretation that links the non-hydrostatic parameter to the amplitude change of baroclinic waves. However, because the non-hydrostatic parameter is essentially the Lagrangian variation of vertical velocity, which includes contributions from both local vertical acceleration and horizontal advections, ε is expected to be strongly system dependent. Even in the simulations presented here, the magnitude of ε varies with the development of the eddies and becomes quite small after the decay of the baroclinic waves.
In Figure 11, we evaluate the non-hydrostatic effect from the perspective of KE dissipation. As described in Taylor et al. (2020), the dynamical core has detailed diagnostics for the kinetic, potential and internal energy, as well as the changes in these terms due to the adiabatic transfers between the different types of energy. The difference between the observed change in KE from each timestep and the expected change due to the adiabatic transfer terms represents the total dissipation introduced by the non-adiabatic components in the model. These sources come from dissipation in the timestepping method, the monotonicity filter in the vertical remapping and the explicitly added hyperviscosity. In practice, hyperviscosity is always the dominant source of KE dissipation. The hyperviscosity in HOMME and HOMME-NH is applied with a time-split approximation that isolates the application of hyperviscosity to its own timestep. We thus measure the dissipation from hyperviscosity by computing the change in KE from the hyperviscosity timestep. In all the cases presented here, including both dry and moist experiments with various reduction factors X, it is responsible for more than 90% (up to ∼95%) of the total dissipation in KE. Thus in the following, we examine the globally-integrated hyperviscosity-induced horizontal KE dissipation D Kν . This quantity is accumulated over time and plotted as a function of the Earth radius reduction factor X. As the accumulation is for the entire 100-day model integration, the cumulative D Kν is closely related to the simulated synoptic baroclinic eddies, as well as the mesoscale and smaller-scale motions that are secondary flow features driven by synoptic-scale shear and strain (e.g., Nadiga, 2014). Compared with EKE that focuses  Figure 2b showing a delay of eddy development at X = 200 although the impact is still quite limited. This is in contrast with the moist simulations that show a slight difference of cumulative D Kν at a coarser effective resolution with X = 100 (Δx eff ≈ 1.1 km, Figure 11b). Although the non-hydrostatic effect becomes influential in moist runs at a coarser effective resolution than in the dry simulations, the difference between hydrostatic and non-hydrostatic experiments with the same X is more prominent in the dry simulations than in their moist counterparts. Besides demonstrating the role of moisture in non-hydrostatic effects, Figure 11 also suggests that under the idealized setups here the effective horizontal grid spacing, at which the non-hydrostatic effect becomes influential, is ∼0.6-1.1 km, much finer than the expected grid spacing on the order of the scale height of the atmosphere (∼10 km).

Sensitivity to Rossby Number, Moisture Content, and Vertical Resolution
We examine the sensitivity to Rossby number (Figure 12), water vapor content ( Figure 13) and vertical resolution (Figures 14 and 15) using the evolution of the mass-weighted vertically averaged EKE for the latitude band from 30° to 60°N. As shown in Figure 12, the baroclinic waves are strengthened in the simulations with increased Δ , where the background Rossby number is increased along with the enhanced zonal mean flow through thermal wind adjustment. The model behavior, however, largely follows the results in Section 5.1 (e.g., Figure 2), including the delayed onset of baroclinic instability at finer effective resolutions in the non-hydrostatic experiments. In the X = 600 non-hydrostatic simulations, the baroclinic waves are also weaker than in the hydrostatic experiments or non-hydrostatic small X runs (coarse effective resolution) where the non-hydrostatic effect is not important. When the moisture content is reduced, Figure 13 is characterized with a similar EKE evolution as in Figure 2, except that the simulated eddies are weaker than in the simulations with more abundant water vapor (cf., Figures 13a and 13c with Figures 13b and 13d).
In the dry hydrostatic simulations that share the same reduction factor but use different numbers of vertical levels, the EKE of the baroclinic eddies is nearly indistinguishable until day 35, except for the simulations with 30 vertical levels (Figure 14, left). In the dry non-hydrostatic simulations, the sensitivity to vertical resolution is still low, although the EKE deviates earlier on day 25 among the runs using different numbers of vertical levels (Figure 14, right). When moisture is included, the sensitivity to vertical resolution is enhanced compared with the dry runs (cf., Figures 14 and 15), probably because mesoscale and smaller-scale activities are more active due to latent heat release. However, in all simulations the deviation induced by the change of vertical resolution mainly occurs within the mature stage and after the decay of the baroclinic waves, during which mesoscale and smaller-scale motions are ubiquitous. More importantly, the delayed onset of baroclinic instability with increased 18 of 26 reduction factor X and the weaker baroclinic waves in the X = 600 runs are always the case no matter how many vertical levels are used in the model. These results, together with Figures 12 and 13, illuminate the robustness of the findings in Section 5.1.
Although the impact of non-hydrostatic dynamics does not strongly depend on vertical grid spacing, the KE spectra increase in magnitude with increasing vertical resolution at large wavenumbers and roughly converge when 120 vertical levels are used ( Figure S8 in Supporting Information S1). The shallowing of the spectra in Figure 13. As in Figure 2, but for the moist simulations with reduced water vapor content as indicated in the panel titles and the legend. Panels (a) and (c) are, respectively, repeated from Figures 2c and 2d for a convenient comparison to the simulations with reduced moisture content in panels (b) and (d). Figure 2, but for the simulations with enhanced temperature difference between the equator and the poles as indicated in the panel titles and the legend. Panels (a, c, and e) are from dry and hydrostatic simulations, while panels (b, d, and f) are from dry and non-hydrostatic runs. In these experiments, the strength of the background zonal flow and its Rossby number are increased. As Diabatic Acceleration and REscaling is exact for the dry hydrostatic simulations, their eddy kinetic energy (EKE) overlaps in panels (a, c, e). Figure 2, but for the dry simulations with various number of vertical levels as indicated in the panel titles and the legend. Panels (a, c, e, g, i, k, m, o) are from dry and hydrostatic simulations, while panels (b, d, f, h, j, l, n, p) are from dry and non-hydrostatic runs. As Diabatic Acceleration and REscaling is exact for the dry hydrostatic simulations, the left panels are the same, although in each panel the eddy kinetic energy (EKE) slightly varies with vertical resolution. SUN ET AL.

Figure 14. As in
10.1029/2022MS003021 20 of 26 Figure 15. As in Figure 14, but for the moist simulations with various number of vertical levels as indicated in the panel titles and the legend. Panels (a, c, e, g, i, k, m, o) are from dry and hydrostatic simulations, while panels (b, d, f, h, j, l, n, p) are from dry and non-hydrostatic runs. Figure S8 in Supporting Information S1 is consistent with Waite (2016, his Figure 9) and Skamarock et al. (2019, their Figure 2), showing that hyperviscosity tends to damp the large-wavenumber KE spectra when the vertical grid spacing is coarse. The convergence of KE spectra for the simulations using 120 or more vertical levels also agrees well with these two studies. With 120 vertical levels, the vertical grid spacing Δz in HOMME and HOMME-NH ranges from ∼230 to 310 m in the upper troposphere and lower stratosphere (5-17 km in height). This range of vertical grid spacing is very close to the required Δz ∼ 200 m to produce converged spectra as proposed by Waite (2016) and Skamarock et al. (2019). Overall, the sensitivity of KE spectra to vertical resolution is more pronounced in the upper troposphere than at low levels (cf., rows 1 and 2, Figure S8 in Supporting Information S1), and the non-hydrostatic dynamics appears to be more influential than moisture as indicated by Figures S8 in Supporting Information S1 (row 1). Figure 16 shows the cumulative dissipation as a function of the number of vertical levels. Specifically, the simulations sharing the same X are analyzed as a group (X-group), and in each X-group the cumulative D Kν is normalized by their asymptotic value that is generally achieved at the finest vertical resolution with 300 vertical levels. The fitting lines follow the relationship ( ) = 1 − 2 − 3 , where D KνN is the normalized cumulative D Kν , the coefficients a 1 , a 2 , and a 3 are obtained using non-linear least squares, and Z N stands for the number of vertical levels. While it is not clear on why the normalized cumulative dissipation is exponentially related to the number of vertical levels, Figure 16 demonstrates that the typically used vertical resolution in global climate models (30-72 levels) underestimates or poorly captures mesoscale and smaller-scale motions, here diagnosed through hyperviscosity-induced dissipation (up to ∼25%). This sensitivity is more obvious under non-hydrostatic dynamics than using the hydrostatic formulation (cf., Figures 16a and 16c with Figures 16b and 16d), and more obvious in the moist experiments than in their dry counterparts (cf., Figures 16a and 16b with Figures 16c and 16d).

Conclusions and Discussion
As a prototype of midlatitude disturbances, baroclinic waves play an important role in weather and climate. In global climate models, these midlatitude eddies are usually simulated with hydrostatic approximation at coarse resolutions, while in the foreseeable future it becomes affordable to conduct high-resolution climate simulations where non-hydrostatic effects become important. In this study, HOMME and HOMME-NH, the dynamical core of the DOE E3SM global climate model, are used to examine the role of non-hydrostatic dynamics in the simulation of dry and moist baroclinic waves and how they are influenced by the background flow properties, by the amount of water vapor in the atmosphere, and by the number of vertical levels. We conduct the study using the DARE approach proposed by Kuang et al. (2005) and detail the needed scaling of the Kessler microphysics in Appendix A. Building upon small-earth experiments to change the effective horizontal resolution, such that the vertical to horizontal aspect ratio, the appeal of this approach to studying the effect of non-hydrostatic dynamics lies in keeping the computational cost minimal over that required for real-earth high-resolution setups.
According to the detailed scale analysis in Section 4, both the large-scale and small-scale circulations in the small-earth runs are not significantly distorted by the change of spherical metric terms for the range of scaling considered. However, shallow convection may be influenced more than deep convection as the former can be potentially transitioned from the hydrostatic regime to the non-hydrostatic regime (Tables 1 and 2). This conclusion is supported by the fact that the simulated baroclinic waves are identical in the dry hydrostatic runs using different Earth radii. Further, spectral flux analysis (cf., Figure 7, Figure S3 in Supporting Information S1) suggests that the inter-scale interactions are not strongly altered either under our highly idealized small-earth setups, where the impact of non-exact moist DARE scaling is also rather limited. However, it should be noted that this may not be the case for more complex or comprehensive configurations.
Compared with the hydrostatic simulations, the onset of baroclinic instability is delayed in the non-hydrostatic experiments, as the associated weaker vertical motions tend to increase the critical wavelength and hence reduce the range of wavelengths that can be baroclinically unstable (Equations 16 and 17; Figure 5). In the model, the delay is up to several days and more prominent as the Earth radius decreases, or equivalently at finer effective horizontal resolutions. This conclusion is consistent with the results from Tokioka (1970, his Figure 2) and Stone (1971, his Figure 1) that are expanded upon Eady′s model.
At very high effective horizontal resolutions, that is, in the experiments with dramatically reduced Earth radius, the baroclinic eddies under non-hydrostatic dynamics are also weaker than in the hydrostatic simulations as found in Stone (1971, his Figure 1). Spectral budget analysis in Figure 7 shows that this results from the suppressed conversion of mean flow APE to eddy KE under non-hydrostatic dynamics. As suggested by the model output ( Figure 8, Figure S7 in Supporting Information S1, Movies S1 and S2) and summarized in Figure 9, non-hydrostatic dynamics tends to induce incremental upward (downward) motions on the leading (trailing) edge of both ridges and troughs. Through vortex stretching and shrinking, these non-hydrostatic vertical motions impact the local vorticity tendency in the upper troposphere and act to accelerate the eastward propagation of upper-level ridges and decelerate the eastward propagation of upper-level troughs. These processes reduce the westward tilt of the vertical ridge axes, suppressing the conversion of mean flow APE to eddy KE that is indispensable for the development of baroclinic waves. Interestingly, further analysis suggests that the mechanism proposed in Figure 9 primarily functions when the baroclinic eddies are close to and within the saturation stage, and shows much weaker impacts during the amplification period (not shown). This is probably why changes in mode structure from the hydrostatic to the non-hydrostatic case are not reported in Tokioka (1970) and Stone (1971) as both studies focus on the onset and amplification processes.
As both the onset and growth rate of baroclinic waves are related to the magnitude of static stability, the delayed onset and the reduced amplitude of baroclinic waves under non-hydrostatic dynamics than in hydrostatic simulations might be also understood from the perspective of effective static stability (Thuburn, 2022, personal communication). According to the linear mode analysis in Davies et al. (2003, their Equation 5.25), when going from hydrostatic to non-hydrostatic dynamics, the static stability for the internal modes switches from N 2 to N 2 − δ v σ 2 , where δ v is equal to unity under non-hydrostatic conditions and σ is the intrinsic frequency. Since σ is imaginary, the non-hydrostatic case effectively experiences a larger base state static stability compared to the hydrostatic case. Although a similar demonstration for baroclinic instability is not straightforward because the corresponding eigenvalue equations are related to the Richardson number (e.g., Tokioka, 1970, his Equation 2-20″;Stone, 1971, his Equations 2.25 and 3.2) rather than N 2 , further exploration in this direction can be interesting.
The contrasts between hydrostatic and non-hydrostatic simulations hold when moisture is present, and more energy is distributed at relatively large wavenumbers as found in previous studies (e.g., Hamilton et al., 2008;Waite & Snyder, 2013). Using the cumulative hyperviscosity-induced KE dissipation as a measure of flow features, however, our results indicate that the non-hydrostatic effect becomes influential in moist runs at a coarser effective horizontal resolution than in the dry simulations, but the difference between hydrostatic and non-hydrostatic simulations at the same effective horizontal resolution is less obvious when moist processes are accounted for. Also, under the idealized setups here the effective horizontal grid spacing at which the non-hydrostatic effect becomes influential is, respectively, ∼0.6 and 1.1 km for the dry and moist cases, much finer than the expected horizontal grid spacing of ∼10 km that is on the order of the scale height of the atmosphere. It should be noted that as the non-hydrostatic parameter ε, which is essentially the Lagrangian variation of vertical velocity, is system dependent, more work is needed to better understand the exact resolution that is needed to fully simulate non-hydrostatic features. For instance, non-hydrostatic effects become influential in convective self-aggregation at a much coarser effective horizontal grid spacing, ∼7 km as demonstrated in Sun et al. (2023).
We also examined the model behavior under different background flow properties, specifically varying its Rossby number, and conducted simulations with reduced moisture content. The model results largely follow the aforementioned conclusions. When the vertical resolution is changed, although the KE spectra vary accordingly as found in previous studies (e.g., Skamarock et al., 2019;Waite, 2016), the conclusions regarding the role of non-hydrostatic dynamics in the simulated baroclinic waves remain unaltered, showing the robustness of our 23 of 26 findings. The mesoscale and smaller-scale activities, however, can be strongly under-represented when the vertical resolution is insufficient as typically used in global climate models (30-72 levels). This under-representation, as measured by the hyperviscosity-induced cumulative dissipation, is as large as ∼25% for the simulations with 30 vertical levels compared with the same simulation but with 300 vertical levels ( Figure 16). Nevertheless, our simulations are limited to simplified physical parameterizations and constrained by the coarse ne30 mesh, leading to a well-resolved wavenumber of k = 120 that corresponds to a wavelength of ∼333 km under real-earth setups. More work using sophisticated simulations that better resolve the mesoscale and smaller-scales is still needed to verify our findings.

Appendix A: Scaling the Kessler Microphysics for DARE
The cloud microphysics in our moist baroclinic wave simulations is represented with a Kessler-type parameterization (Kessler, 1969), which contains three water species in terms of mixing ratio, including water vapor q v , cloud water q c , and rainwater q r . According to Klemp and Wilhelmson (1978, henceforth KW78), the continuous form of the governing equations that are implemented in numerical models, for example, the E3SM, can be written as where L v is the latent heat of vapourization, q vs is saturation water vapor mixing ratio, ρ is the density of moist air, V r is the terminal velocity of rainwater, with A r , C r , and E r as the rates of autoconversion, collection and evaporation of rain, respectively. At the end of each physics timestep, the Kessler scheme is applied to adjust θ, q v , q c , and q r , and diagnose either accumulated precipitation as in Klemp et al. (2015, their Appendix C) or precipitation rate P as in E3SM using where ρ d is the density of dry air, while ρ w is the density of liquid water. In Equation A5, all variables refer to their values at the lowest model level.
When solving Equations A1-A4 numerically, however, one has to follow a two-step approximation as proposed by Soong and Ogura (1973, their Appendix) to determine . This procedure introduces certain approximations, for example, the use of a linearized Teten′s formula (Equation 3.9 of KW78). As a result, the discrete form of Equations A1-A4 becomes where Δ = * − 1+ * ( = 4,093 (Π * −36) 2 ) and the asterisk denotes the values that an adiabatically lifted air parcel would have at its exact saturation (see Figure A1 of Soong & Ogura, 1973), S stands for the discrete form of the sedimentation term 1 ( ) , and Δ = − − Δ 1+ 2 0.875 Δ (k 2 = 2.2 s −1 ) indicates the change of q r due to autoconversion and collection. It should be noted that in the derivation of Δq rAC , a simple linearization using Taylor′s expansion is applied. Thus, the scaling of cloud microphysics for Diabatic Acceleration and REscaling (DARE) simulations should be based on Equations A5-A9, rather than Equations A1-A5. Specifically, the scaling involves which are further constrained by the following relationships to ensure that the condensation of q v , autoconversion and collection of q c to q r , and evaporation of q c and q r do not exceed the amount of available q v , q c , and q r in the atmosphere. Per Equation A5, the diagnosed precipitation rate is also scaled X times in DARE simulations. However, it should be noted that because of the physical constraints, as reflected by the min and max operators in Equations A14-A16, an exact moist DARE scaling is not necessarily guaranteed.

Data Availability Statement
The E3SM model is available at https://e3sm.org/model/.