An Assessment of Nonhydrostatic and Hydrostatic Dynamical Cores at Seasonal Time Scales in the Energy Exascale Earth System Model (E3SM)

In global atmospheric modeling, the differences between nonhydrostatic (NH) and hydrostatic (H) dynamical cores are negligible in dry simulations when grid spacing is larger than 10 km. However, recent studies suggest that those differences can be significant at far coarser resolution when moisture is included. To better understand how NH and H differences manifest in global fields, we perform and analyze an ensemble of 28 and 13 km seasonal simulations with the NH and H dynamical cores in the Energy Exascale Earth System Model global atmosphere model, where the differences between H and NH configurations are minimized. A set of idealized rising bubble experiments is also conducted to further investigate the differences. Although NH and H differences are not significant in global statistics and zonal averages, significant differences in precipitation amount and patterns are observed in parts of the tropics. The most prominent differences emerge near India and the Western Pacific in the boreal summer, and the central‐southern Indian Ocean and Pacific in the boreal winter. Tropical differences influence surrounding regions through modification of the regional circulation and can propagate to the extratropics, leading to significant temperature and geopotential differences over the middle to high latitudes. While the dry bubble experiments show negligible deviation between H and NH dynamics until grid spacing is below 6.25 km, precipitation amount and vertical velocity are different in the moist case even at 25 km resolution.

Many studies have sought to identify the grid spacing at which NH effects emerge, but the conclusions have varied depending on the model they used and the specific phenomena under consideration. Most previous studies have taken place in the context of idealized simulations such as the Held-Suarez dry test (Wedi & Smolarkiewicz, 2009) or for regional short-term numerical weather forecasts (Dudhia, 1993). Traditionally, this grid spacing has been estimated as 10 km in the absence of latent heat release (Kato & Saito, 1995). Orlanski (1981) showed that the hydrostatic approximation is invalid for grid resolutions finer than about 8 km in the simulation of moist convection in a mesoscale model; however, Dudhia (1993) compared approximately 1-day simulations using both NH and H versions of the Penn State-NCAR Mesoscale Model version 5 and found insignificant NH effects at a horizontal grid spacing of 6.67 km. In conducting numerical weather prediction experiments for a realistic case over California, Janjic et al. (2001) found that NH effects had only a weak influence at 8 km resolution. However, they observed a visible effect on orographic precipitation and the NH deviation of pressure made a significant small-scale contribution to the pressure gradient force. Jeevanjee (2017) examined the convective vertical velocity in the simple doubly periodic radiative-convective equilibrium simulations using NH and H dynamical cores. Differences between NH and H solvers were virtually undetectable at resolutions of 2 km or coarser, while the hydrostatic solver was found to overestimate convective vertical velocity by a factor of 2-3 when the resolutions are finer than 2 km.
Although the aforementioned work showed that H and NH differences were small on the timescales of numerical weather prediction, few studies have examined NH effects in the context of moist seasonal to decadal simulations. In fact, a number of recent studies have suggested that seasonal or climatological differences between NH and H dynamical cores could be significant even at resolutions much coarser than 10 km. Gao et al. (2017) simulated two summer seasons in the years 2006 and 2007 over the continental U.S. and the surrounding land and ocean areas in the Weather Research and Forecasting model (WRF). They found that NH simulations reduce the large wet bias over the western U.S. at 12 and 36 km resolutions when a convection scheme is used. Maurya et al. (2020) compared the performance of NH and H dynamical cores in a regional climate model over the India at 12 and 36 km resolutions in three monsoon seasons. The results revealed strong regionality with an inhomogeneous pattern emerging over the domain, as well as inconsistency among the three monsoon seasons, demonstrating the complexity of real-world simulations compared to the idealized simulations. Yang et al. (2017) compared the NH and H solvers of the WRF model in 3-month regional climate simulations using an aquaplanet framework with and without mountains. In the simulations without mountains, the differences in the time-and domain-averaged precipitation between NH and H dynamical cores was negligible at 36 and 12 km while more precipitation was found in the H than in the NH simulations over the intertropical convergence zone (ITCZ) at 12 km resolution when the cumulus parameterization (CP) is activated. They found that simulated differences between NH and H dynamics exhibit a strong sensitivity to whether a CP is active or not. In the simulations with an idealized mountain, the differences between NH and H are significant when it came to precipitation and circulations from 36 to 4 km over the tropics. Notably, these differences are more sensitive to the mountain half width than to the model grid spacing. A potential caveat is that some differences related to damping, diffusion, and filtering associated with the dynamical cores may have contributed to the differences between NH and H simulations, although attempts were made to minimize these differences. One advantage of our study is that the NH and H models used in the comparison are essentially identical except for the addition of the NH term. This study is able to more definitively attribute the simulation differences to NH terms because the hydrostatic version of the model explored in this study is implemented by simply excluding the NH term.
Although several NH limited-area simulations were discussed above, studies employing non-idealized NH models globally and on seasonal to decadal timescales are rare (Satoh et al., 2008). Given that NH effects vary across climate regimes, regional simulations are insufficient to fully characterize the differences between NH and H dynamical cores. Consequently, several unresolved questions remain when it comes to differences between NH and H dynamical cores. For example, are there particular regions where differences between NH and H dynamical cores are more significant? How do the results from theoretical studies (Jeevanjee, 2017;Yang et al., 2017) manifest in real-world simulations? At what grid spacing should hydrostatic models no longer be employed for seasonal to decadal simulations? Given the rapid progress being made toward producing ensembles of model simulations with grid spacing near 10 km, there is significant value in investigating the differences between these dynamical core formulations in this regime.

Model Setup
A series of global simulations are conducted using a newly developed global cloud permitting model called the Simple Cloud-Resolving E3SM Atmosphere Model (SCREAM), as part of the Energy Exascale Earth System Model (E3SM) project. Prescribed SST and sea ice extent are used for all simulations. The atmospheric dynamics is solved using the High Order Method Modeling Environment (HOMME; J. Dennis et al., 2005; J. M. Dennis et al., 2012;Evans et al., 2013), which provides both the hydrostatic (H) spectral element dynamical core used by E3SMv1 Golaz et al., 2019;Rasch et al., 2019) and the Community Earth System Model (Lauritzen et al., 2018;S. Zhang et al., 2020;Small et al., 2014), as well as the new HOMME NH (HOMME-NH; Bertagna et al., 2020;Taylor et al., 2020) spectral element dynamical core developed for SCREAM. In this section, we focus specifically on the differences between the H and NH model. More details about SCREAM are provided in Caldwell et al. (2021).
As mentioned in Yang et al. (2017), some differences related to model configuration (e.g., damping, diffusion, and filtering associated with the dynamical cores) may hinder the comparison of H and NH dynamical cores. Small deviations could be further amplified by feedbacks associated with moist processes. Thus, an effective and targeted study of the differences between H and NH solvers should aim to keep all other possible cofactors fixed. SCREAM, with the HOMME dynamical core in its hydrostatic and NH configuration, has such a capability.
The HOMME-NH dynamical core employs the equations described in Taylor et al. (2020). It includes prognostic equations for the three components of the velocity field, the mass-coordinate pseudo-density, the geopotential height, and the virtual potential temperature. Specifically, the prognostic equations for the vertical velocity, w, and the geopotential height, ϕ, are where u = (u, v) is the horizontal velocity, s is a terrain following coordinate = ∕ , where D/Dt is the standard material derivative, and g is the gravitational acceleration. The deviation from the hydrostatic pressure, μ, is defined as where p is the NH pressure and p h is the hydrostatic pressure. In the hydrostatic version of HOMME (denoted HOMME-H), there is no prognostic equation for w and the geopotential is diagnosed from the hydrostatic equation.
For the time-stepping algorithm, the H simulations utilize the third order KGU53 explicit method described in Guerra and Ullrich (2016). NH simulations use a Horizontally Explicit Vertically Implicit approach (HEVI; Satoh, 2002), discretized with an IMplicit-EXplicit (IMEX) Runge Kutta method (Ascher et al., 1997), which combines an explicit KGU53 method for all horizontal derivatives and an implicit method for terms responsible for vertically propagating acoustic waves, similar to the IMEX methods constructed in Guba et al. (2020).
Two horizontal dynamics grids are used in this study, both based on the unstructured cubed-sphere finite element grid with four Gauss-Lobatto-Legendre nodes per element's edge. Our coarser resolution simulation uses 120 finite elements along the edge of each cube face, corresponding to a horizontal grid spacing of 0.25° (∼28 km). This resolution is abbreviated as ne120. Our finer resolution simulation uses 256 finite elements along the edge of each cube face, corresponding to a horizontal grid spacing of 0.12° (∼13 km). The corresponding fine resolution is abbreviated as ne256. At ne120 resolution, the dynamics timestep and the physics timestep are 75 and 900 s, respectively. Simulations at ne256 resolution use a 37.5 s dynamics timestep and a 600 s physics timestep. Native output is saved every 3 hr and later remapped to a uniform 1° finite volume grid using the TempestRemap software suite (Ullrich, Devendran, & Johansen, 2016;Ullrich & Taylor, 2015) before calculating derived variables or performing analyses. Most simulations use the low-resolution (LR) configuration, which is typical for simulations conducted with unresolved convection; the model is run with 72 vertical levels with a top at 60 km, and the Zhang-McFarlane deep convection scheme (G. J. Zhang & McFarlane, 1995) is applied. A short test using the high-resolution configuration (where deep convection is turned off) is also conducted and discussed in Supporting Information S1. These two configurations are described in Caldwell et al. (2019Caldwell et al. ( , 2021, respectively. SCREAM at ne120 and ne256 resolutions using the LR configuration has not been tuned to minimize climate biases. Therefore, this study is better suited for exploring the NH and H differences rather than assessing the ability of SCREAM to match observations. Four sets of ne120 (28 km) simulations are designed as shown in Table 1. Similar to Gao et al. (2017), we conduct continuous simulations (i.e., Summer-NH and Summer-H) from May to August 31 of 2006 and analyzed only the summer period (1 June to 31 August) with May as a spin-up period. In addition, we performed 10 ensemble simulations for each ne120 configuration to validate our results since the simulations could be sensitive to small uncertainties in the initial conditions. The 10 ensemble members are generated by initializing the model on successive days (1 May to 10 May). Past studies Lebassi-Habtezion & Diffenbaugh, 2013) Table 1). The simulations at ne256 resolution only contain three ensemble members considering the high computational cost. The abbreviation NH-H, is used throughout this manuscript to represent the difference between simulations with the NH dynamical core and simulations with its H counterpart.
Atmospheric initial conditions and SSTs are derived from the fifth generation of atmospheric reanalysis, ERA5 (Hersbach et al., 2019), produced by the European Centre for Medium-range Weather Forecasts. The land initial file is generated from a 6-month spinup simulation prior to the initial date.

Data and Statistical Tests
We rely on the hourly ERA5 reanalysis at 0.25° resolution to evaluate our model performance. For precipitation, we use daily Global Precipitation Climatology Project (GPCP; Adler et al., 2018;Huffman et al., 1997) remapped to 1° horizontal resolution and the 3-hr Tropical Rainfall Measurement Mission Multi-satellite Precipitation Analysis (TRMM; Huffman et al., 2007) at 0.25° resolution.
We evaluate significance throughout this manuscript at the 99% confidence level using Student's t-test. The significance of the seasonal mean is calculated by examining the 3-hr time series at each 1° grid point. Additionally, two statistical methods are tested for the purpose of validating our findings. First, the Student's t-test is conducted on daily time series, which results in a smaller sample size than the one using 3-hr data. Second, we use the effective sample size to replace the traditional sample size given the fact that precipitation time series are usually highly autocorrelated. The effective sample size, n′, can be estimated using the approximation (Wilks, 2011) where n is the sample size and ρ 1 is the lag-1 autocorrelation coefficient for the time series. These two additional statistical methods are found to have no significant impact on the detection of the hot spots in NH-H precipitation.

Realistic Seasonal Simulations
We begin our discussion with the differences between NH and H dynamical cores in the real-world simulations on the seasonal timescale. Following the three scientific questions raised in Section 1, we develop hypotheses based on past studies. The hypotheses are tested in the SCREAM simulations and answers to the scientific questions are provided. The following Sections 3.1 and 3.2 provide analyses during the boreal summer and winter seasons, respectively. The majority of the simulations are at ne120 (28 km) resolution with some sensitivity tests at ne256 (13 km) resolution in Section 3.1.8.

Summer Simulations
This section is organized as follows: simulated precipitation is evaluated and analyzed in Sections 3.1.1-3.1.3; NH and H's differences in other fields and regional circulations are analyzed in Sections 3.1.5 and 3.1.6; Section 3.1.7 presents changes in the tropical cyclones (TCs); and Section 3.1.8 describes the sensitivity to horizontal resolutions.

Model Performance
To first understand model skill, The precipitation bias in the simulation compared to the GPCP is shown in Figures 1e and 1f. The model precipitation bias is larger than the difference between NH and H simulations ( Figure 1d) in general. Both NH and H simulations have wet biases over the ITCZ, in particular near the Maritime Continent and western Pacific, which remains a persistent problem in high resolution climate models including HighResMIP Caldwell et al., 2019Caldwell et al., , 2021McClean et al., 2011;Roberts et al., 2019). NH does not exhibit better performance than its H counterpart from the global perspective at 28 km. Specifically, both NH and H models show a global pattern correlation between simulated precipitation and GPCP of approximately 0.87. NH reduces the wet bias

Regions of Significant NH-H Precipitation Difference
Figure 1d displays the precipitation difference between Summer-NH and Summer-H simulations. To simplify the results, only the regions with significant differences at the 99% confidence level (Section 2.2) are shown, while the insignificant regions are masked as white. Additionally, regionally averaged precipitation for selected regions in the H simulation and associated NH-H difference is quantified in Table 2.
Past studies Janjic et al., 2001;Lebassi-Habtezion & Diffenbaugh, 2013) argue that the NH and H models largely differ because of topographically induced differences. Therefore, one might initially expect that hotspots in the precipitation difference between NH and H simulation are located over topographic features such as the mountainous Western U.S. However, this is clearly not the case in our simulations-none of our hotspots can be found over the mountain ranges (e.g., the Rocky Mountains). Although drying is statistically significant over the Western-Southern U.S. and Alpine Europe in the NH simulation, the magnitude of the rainfall differences is small (e.g., 0.12 mm/day over the Western U.S. and 0.14 mm/day over the Alpine Europe; Table 2). Figure 2a shows zonally averaged precipitation in Summer-H, Summer-NH, GPCP, and ERA5. The meridional distribution of precipitation is well reproduced by the models relative to observations. The simulations are wetter than GPCP by 0-2 mm/day and agree well with ERA5 except for a wet bias near 40°S by 1 mm/day. NH and H do not show significant differences in their global zonal means because the complicated zonal structures in the NH-H precipitation leads to cancellation of terms ( Figure 1d).
In examining Figure 1d, absolute precipitation differences are small nearly everywhere except in the vicinity of the northern hemisphere summer ITCZ. And within that band, the most prominent differences emerge near India and the Western Pacific. An analysis region which includes these two hotspots is defined as (60°E-165°W and 5°-40°N; shown as the black box in Figure 1d). These two identified hotspots appear even when using six-member subsets of the complete ensemble (not shown), which increases the confidence in our results. The hotspot Alpine Europe (0°-30°E and 40°-50°N) 2.22 −0.14 −6.31 regions are associated with pronounced differences in not only the seasonal averaged precipitation but also the standard deviation of the precipitation (not shown), although this is perhaps unsurprising as precipitation standard deviation generally scales with its mean. The precipitation in these regions is significant: precipitation over the northern Indian Ocean is relevant to prediction and analysis of the Indian Monsoon from synoptic to seasonal timescales (Ashok et al., 2001;Schott & McCreary, 2001), and the western Pacific is also an influential region in the development of extreme rainfall events (e.g., typhoons and subsequent flooding), which can have significant impacts, particularly on southeastern Asian countries (Li & Zhou, 2012;Troselj et al., 2017). However, these two hotspots need to be analyzed separately as, curiously, they exhibit opposing patterns of enhanced and suppressed precipitation.

Table 2 Regionally Averaged Precipitation (mm/day) in the Summer-H Simulation and NH-H Difference in the
The northern Indian Ocean-India region is characterized by a meridional dipole in the NH-H precipitation difference anomaly. Specifically, NH produces 20.5% more precipitation over the Central India and the surrounding oceans and is drier than the H simulation in the South India-northern Indian Ocean by 4.8% (Table 2). For the western North Pacific region, precipitation maxima in the H and NH simulations occur at the same location ∼ (135°E and 10°N) while the NH model shows a higher rainfall rate than H by 11.6%. There is a transition region (145°E-170°W and 20°-35°N) from the tropics to subtropics where the NH simulations display less precipitation than H by 16.12%. While the precipitation differences in these two hotspots is significant, precipitation amounts in the analysis region in both H and NH simulations are nearly identical (7.71 mm/day and 7.74 mm/ day, respectively). This is likely because, on climatological time scales, precipitation amounts are constrained by available energy in the system (Held & Soden, 2006). Consequently, it makes sense that NH and H differences in precipitation primarily manifest as shifts in precipitation patterns rather than changes in precipitation amounts. In this case, those shifts are manifest via a northward shift in the ITCZ in the Indian Ocean-India region and a southward shift in the western Pacific.
Figures 2b and 2c show the zonally averaged precipitation averaged over the western Pacific (120°E-170°W) and Indian Ocean (60°-100°E), respectively. The averaging zones are selected to focus on the regions where NH-H precipitation is significant ( Figure 1d). The yellow shading represents the standard deviation of the 10 ensembles in the NH simulations. There is significantly higher rainfall rate in the NH simulations over 10°-20°N and less rainfall over 20°-35°N in the western Pacific with anomalous magnitudes larger than the standard deviation of the ensembles (Figure 2b). In the Indian Ocean (Figure 2c), the NH precipitation maximum location is shifted northward by about 3°, even though the magnitudes of the maxima in NH and H simulations are similar (∼16 mm/day).
To summarize Figures 1 and 2, although the differences on the global scale are not significant at 28 km horizontal resolution, there are significant precipitation differences on the regional scale over the moist tropics. The precipitation difference is presented by the spatial patterns rather than the averaged rainfall rates. Unlike past studies Janjic et al., 2001;Lebassi-Habtezion & Diffenbaugh, 2013) showing the NH effects are more influential over the topography, the tropical oceans where heavy rainfall occurs are more sensitive to the selection of dynamical cores in our study.

Precipitation Differences From Convective Versus Large-Scale Processes
There are two components to total precipitation (PRECT) in the model: large-scale (resolved) precipitation (PRECL) and convective precipitation (PRECC), which is unresolved precipitation from the subgrid convective parameterization. Table 3 displays the percentage of the total precipitation that is convective precipitation and large-scale precipitation in the NH and H simulations averaged over the hotspots identified in Figure 1d. The precipitation is averaged over the regions before calculation of the percentage. Over the four selected regions, total precipitation is dominated by the convective precipitation with percentages ranging from 55% to 78% while the large-scale precipitation only contributes to the total precipitation by 22%-45%. Comparing the contributions of convective/large-scale precipitation to the total precipitation between NH and H simulations, the differences are minor (<4%).  The total precipitation difference between NH and H models is written as where subscript H and NH denoting the simulation ensemble. Using Δ to denote the difference between NH and H simulations (NH minus H), this equation can also be written as To further understand the precipitation differences between NH and H simulations, Table 4 displays regionally averaged NH-H total precipitation (ΔPRECT), convective precipitation (ΔPRECC), and large-scale precipitation (ΔPRECL), and their percentages. The convective and large-scale precipitation partitions contribute roughly equally to the total precipitation difference between NH and H simulations, with contributions ranging from 45% to 55%. The magnitude difference between ΔPRECC and ΔPRECL is less than 0.16 mm/day (<10%).

Precipitation Differences From Moisture Versus Uplift
Precipitation variability is also influenced by both moisture and circulation (Chou et al., 2009;Dong et al., 2018;Held & Soden, 2006;Huang & Xie, 2015;Huang et al., 2013;McClenny et al., 2020;Seager et al., 2010). Consequently, we now turn our attention to addressing whether the NH-H precipitation differences are caused by low-level moisture difference or changes in the circulation. Using a simplified decomposition (Dong et al., 2018;Huang & Xie, 2015;Huang et al., 2013), the precipitation difference between NH and H simulations is written as where P is precipitation, ω is 500 hPa pressure vertical pressure velocity (negative ω means upward flow), q is surface specific humidity, and Δ denotes the difference between NH and H. The two terms on the right-hand side can be interpreted as follows: q ⋅ Δω is a dynamics component due to differences in vertical motions and Δq ⋅ ω is a component due to the differences in surface moisture. Note that the units on the two sides of Equation 7 differ by a constant coefficient, which is omitted for simplicity as in previous studies (Huang & Xie, 2015;Huang et al., 2013).    term (q ⋅ Δω) resembles the spatial pattern of NH-H precipitation ( Figure 1d) with a correlation ∼−0.77. The regional pattern correlation over the analysis region (the black box in Figure 1d) is approximately −0.9. By contrast, the moisture term Δq ⋅ ω is not significantly correlated to NH-H precipitation (correlation magnitude <0.13). Moreover, the global mean of Δq ⋅ ω is one order smaller than q ⋅ Δω in absolute value. Comparison of the two terms indicates that the precipitation differences between NH and H simulations are related to differences in vertical velocity fields rather than changes in the surface moisture.

Differences in Other Fields
Besides the precipitation differences, there are a few significant differences in other atmospheric variables that should be examined. These differences are of interest for two reasons. First, the analyses of atmospheric variables can provide insights into the physical mechanisms of the precipitation differences and the interactions between atmospheric processes and precipitation. Second, NH effects may emerge over certain regions without significant precipitation, which is not investigated and discussed in previous results. Figure 4 shows differences between NH and H simulations in 850 hPa temperature and specific humidity, vertically integrated precipitable water (TMQ), sea level pressure (PSL), and 500 hPa geopotential height and vertical pressure velocity. In the Indian and western Pacific region, low-level temperature is slightly reduced (0-0.4 K) in the NH simulations ( Figure 4a) which we associate enhanced precipitation. Low-level specific humidity does not differ significantly over the India and western Pacific, consistent with the results in Figure 3. The geopotential height difference over the northern Indian Ocean is about 10 times weaker than the northern Pacific.
Differences in the temperature and geopotential height emerge over middle to high latitudes related to planetary-scale atmospheric Rossby waves (Figures 4a, 4d and 4e). The observed differences in midlatitudinal Rossby wave activity are likely driven by the west Pacific circulation differences, and corresponding propagation of this signal to the extratropics (Yadav, 2017;Yun et al., 2011). Notably, the resulting difference pattern appears similar to the negative phase of the Arctic Oscillation (L'Heureux et al., 2010;Minami & Takaya, 2020). A strong geopotential height dipole exists over the Northern Pacific (130°E-170°W and 20°-60°N) extending vertically from the surface to 500 hPa (Figures 4d and 4e). The anomalous high pressure over 20°-40°N extends the Pacific High westward and is associated with decreased upward vertical motions ( Figure 4f) and less precipitation (Figure 1d) in the NH simulation.
Although there are circulation differences over middle-high latitudes, precipitation differences in the extratropics are smaller as less water vapor is available. With that said, the atmosphere is drier in the NH simulation over the continental U.S. as seen in the 850 hPa specific humidity field (Figure 4b), TMQ (Figure 4c), and 2 m specific humidity (not shown), consistent with Gao et al. (2017). The dry anomaly reduces the wet bias over the continental U.S., which in turn leads to a regional improvement in the model performance by the NH dynamical core. The dryness is related to the regional circulation modifications, especially in the Great Plains low-level jet, which will be discussed in the next section.

Modifications to the Meridional Circulation
Figures 5a and 5b show 850 hPa wind and geopotential height in the H simulation and NH-H difference. Figures 5c-5f display the vertically averaged mean meridional streamfunction (MMS) over the Pacific and Indian Ocean. Given that NH-H precipitation is structured more in the meridional direction (Figure 1d), the MMS index is calculated to quantify the meridional circulations. It is defined as vertically integrated meridional mass flux between a given pressure and the top of the atmosphere (Davis & Birner, 2016;Holton, 1994): where g = 9.81 m/s is the acceleration due to gravity, a = 6,371 km is Earth's mean radius, p is the pressure, ϕ is latitude, [v] is the zonal-mean meridional wind, and Ψ(p, ϕ) is the MMS at the given pressure and latitude. Positive MMS represents a clockwise meridional circulation. Note that the large-scale meridional circulations (such as the Hadley cell) are well captured by the model using [v] averaged over all longitudes in terms of the locations of the cells and the magnitudes of the MMS (not shown). However, [v] is averaged over only certain longitudes in Figures 5c-5f since the NH effect is more on the regional scale than the global scale. The averaging longitudes are 120°E-165°W over the Pacific in Figures 5c and 5d and 60°-120°E over the Indian Ocean in Figures 5e and 5f.
We will start by discussing the western North Pacific region. The moist monsoon flows from the India and easterlies from the Pacific convergent over the western North Pacific region (Figure 5a) producing high rainfall rates in both NH and H simulations (Figures 5a-5c). 850 hPa geopotential height using the NH model shows an anomalous high-pressure system associated with strong anticyclonic circulation over the western Pacific between 10° and 40°N compared to the H counterpart (Figure 5b). The anticyclone is associated with diminished upward vertical velocity and less rainfall. Furthermore, the anticyclonic circulation intensifies the easterly from the Pacific to the Maritime Continent, which leads to more convergence near (140°W and 10°N) and, therefore, increased rainfall. The NH-H MMS clearly show two anomalous circulations with opposite directions over 20°S-10°N and 10°-30°N (Figure 5d). The upward branches of the two meridional circulations are overlapped around 10°-20°N, consistent with the negative anomaly in 500 hPa pressure vertical velocity ω (Figure 4f). In the northern Indian Ocean-India region, the model simulation captures the Indian Monsoon circulation accurately (Figure 5a) including the southeasterly winds from the southern Indian Ocean and southwesterlies after crossing the equator to the northern hemisphere, which supply large amounts of moisture supporting the Indian Monsoon rainfall. The southwest monsoon flows from the Arabian Sea to the central India are intensified in the NH simulations (Figure 5b). The stronger convergences of winds over the India are associated with the increased precipitation in the NH simulations ( Figure 1d) and the poleward shift of the Indian precipitation maximum location (Figure 2c). The NH-H MMS cross section shows an anticlockwise circulation around 10°S-25°N with an upward branch near 20°N, consistent with the negative anomaly in pressure vertical velocity ω (Figure 4f).
From the Gulf of Mexico to the central U.S., the Great Plains low-level jet (LLJ) is simulated well (Figure 5a). The warm season precipitation in the Southeast, South-central, and Midwest U.S. is known to be affected by the LLJ, which supplies moisture originating in the Gulf of Mexico (Pu et al., 2016;Schubert et al., 1998;Weaver & Nigam, 2008). The LLJ extends from about 20°N over the Gulf of Mexico to the northern Great Plains and Midwest, with its core located over the southern Great Plains. The dry anomaly over the U.S. observed in the NH simulation (Figures 4b and 4c) is associated with the modification of low-level circulation. The southerly inflow over the Gulf of Mexico is reduced while there is enhanced southerly winds in the northern portion of the Great Plains LLJ (Figure 5b). Corresponding to these wind anomalies, anomalous negative vorticity is centered over the Gulf States (30°-35°N and 90°-100°W) suppressing precipitation. Meanwhile, the equatorward flows associated with the North Pacific High is enhanced over the California coast in the NH simulation related to the intensified Pacific High.

Tropical Cyclones
As features strongly associated with convection, we expect TCs to be particularly sensitive to the choice of dynamical core (Roberts et al., 2020). To this end, we track TCs using TempestExtremes (Ullrich & Zarzycki, 2017;Ullrich et al., 2021) and assess differences in the composite horizontal structure of these systems between the NH and H model. Notably, Qi et al. (2018) previously investigated NH and H differences in TC structure in WRF at 1.67 km grid spacing. Figure 6 shows composite precipitation and 500 hPa vertical pressure velocity from western Pacific TCs. We focus on this region since it accounts for the vast majority of storms (76 out of 155 in the H ensemble, and 74 out of 139 in the NH ensemble) and is identified as a major hotspot in Figure 1d. The left column shows the composite H fields while the right column shows the difference between NH and H fields. Stippling indicates regions which are significantly different at the 90% level. As in Villarini et al. (2014) Figure 1, we find that the most highly precipitating region occurs to the south-southwest of the eye.
As in the regional climatology, we observe that regions with the strongest vertical velocity are the most sensitive to the choice of NH or H dynamical core. Precipitation is enhanced by around 10% in the southwest and southeast of the TC core, where updrafts appear strongest in the model. Similarly, increased subsidence to the north of the core leads to a reduction of precipitation in this region. As a result, the NH dynamical core produces a greater asymmetry in its vertical structure. Note that while the central sea-level pressure is slightly weaker in the NH dycore (∼40 Pa at the center of the eye), the difference is not significant at the 90% level. This suggests that H and NH differences occur primarily above the surface level. Our results in this section are somewhat difficult to compare to Qi et al. (2018), whose simulations are performed at much finer resolution. However, consistent with our results, they also observe that the NH model produces generally stronger updrafts in the eyewall, with compensating subsidence that suppresses precipitation in the outer bands of the storm. Table 1) and precipitation difference between NH and H simulations at 13 km (ne256) resolution. Our goal in this analysis is to determine if the pattern and magnitude of NH-H precipitation change is preserved as the model is pushed to finer resolution. The precipitation shown is a three-member ensemble mean. Comparing the H precipitation at 28 (Figure 1a) and 13 km (Figure 7a) resolutions, the global precipitation patterns are similar while the rainfall over the western Pacific is reduced at 13 km by up to 3-4 mm/day implying a smaller regional wet bias. The maxima of NH-H precipitation in 13 km simulations are located over the Bay of Bengal and western Pacific (Figure 7b). NH simulation shows lower rainfall rates than H simulation by approximately 4 mm/day over the Bay of Bengal, decreasing the regional wet bias in the simulation compared to the H simulations. The NH-H precipitation patterns are not identical in 28 and 13 km simulations over some regions. For example, the dry anomaly over the northern Pacific (140°-180°E and 20°-40°N) in the 28 km simulation is not present in the 13 km simulation. The northern Pacific dry anomaly does emerge in one member of the 13 km ensemble, but it does not appear to be significant in the three-ensemble mean. Considering the H precipitation also differs at 28 and 13 km resolutions, the different patterns in NH-H precipitation between the two resolutions are likely related to the differences in precipitation climatology, especially near the western Pacific ITCZ. In spite of the differences between the two resolutions, the hotspots of NH-H precipitation still appear near the Indian region and western Pacific indicating the two regions are sensitive to the NH effects in the boreal summer at resolutions from 13 to 28 km.

Winter Simulations
So as to understand if those regions experiencing the most significant differences between NH and H dycores are sensitive to season, we now analyze H and NH ten-member boreal winter simulation ensembles and compare results with the boreal summer. Figure 8 shows precipitation in Winter-H, Winter-NH, and GPCP, and their differences averaged over 1 December to 28 February 28. The simulated precipitation is averaged over the 10 ensembles and the GPCP precipitation is averaged over years 2001-2010. The precipitation patterns in our simulations resemble the GPCP with a global pattern correlation of 0.88 (Figures 8a-8c). The global pattern correlations of GPCP and model rainfall using NH and H dynamical cores are not significantly different (<0.0002). Wet biases are found over the central-southern Indian Ocean (45°-80°E and 30°S-10°N) and South Pacific Convergence Zone (SPCZ), and the northern Pacific and Atlantic over 25°-40°N by 0-6 mm/day (Figures 8e and 8f). The precipitation over the Amazon rain forest is underestimated by 0-3 mm/day, similar to other climate models (Yin et al., 2013). Figure 8d with regionally averaged precipitation quantified in Table 5. Consistent with the boreal summer season (Figure 1d), the hotspots of NH-H precipitation appear unrelated to topography. The NH-H precipitation maxima are located in the southern Indian Ocean and central-southern Pacific. Simulations using the NH dynamical core show significantly heavier precipitation over the southern Indian Ocean (45°-90°E and 10°-30°S) and the central Pacific (165°E-150°W and 0°-10°N) than the H counterpart by 14%-16%.

NH-H precipitation in the boreal winter season is shown in
In contrast, NH model experiences less precipitation over the central Indian Ocean (40°-100°E and 0°-10°S), the Maritime Continent (120°-150°E and 15°S-10°N), and the southern Pacific (120°-180°W and 10°-25°S) by 7%-11%. In fact, the pattern of drier/wetter regions along the ITCZ is flipped compared to the summer simulation in the western Pacific and Indian Ocean regions, again suggesting that the ITCZ is more heavily tilted poleward in the NH model in this region as one goes from west to east.
In the remainder of this section we will thus focus on the analysis region (30°E-120°W and 30°S-10°N; shown as the black box in Figure 8d). The boreal winter simulations support the conclusion of the summer simulations (Section 3.1.2) that the precipitation differences between NH and H dynamical cores in the seasonal simulations occur primarily on the regional scale, whereas differences on the global scale are not significant at ∼28 km horizontal resolution. Clearly, the dynamical core largely yields different spatial patterns over the hotspot regions and the inclusion of a land surface or topography confounds global-scale statistics compared with the idealized simulations (Yang et al., 2017). Figure 9 shows differences between NH and H simulations in 850 hPa temperature and specific humidity, TMQ, PSL, and 500 hPa geopotential height and vertical pressure velocity using the same color scales as Figure 4. Not surprisingly, the patterns of TMQ ( Figure 4c) and 500 hPa vertical velocity (Figure 4f) are highly correlated to the NH-H precipitation over tropics. Similar to the boreal summer simulation (Section 3.1.5), the temperature and geopotential height anomalies are modest at low latitudes and larger at high latitudes, especially over 40°-90°N (Figures 9a, 9d and 9e). The NH simulation shows a lower 850 hPa temperature (and 2 m temperature; not shown) over the Russia, Alaska, and eastern Canada by up to 2.5 K, associated with a negative geopotential height anomaly. Similar patterns are also present in the 200 hPa geopotential height (not shown) showing a modification of wintertime Rossby wave activity. The magnitudes of anomalies are larger in the winter season than summer. Those geopotential height differences in turn do produce changes in the wind fields, especially in meridional winds (not shown).  Figures 10e and 10f). Over the southern Indian Ocean, the Mascarene High is centered at (85°E and 28°S) in the H simulation (Figure 10a). The westerly originating from the retreating Indian Monsoon and the southeasterly from the Mascarene High strongly converge over the central Indian Ocean from 10° to 20°S (Figure 10a), associated with the regional precipitation maximum (Figure 8d). The positive geopotential height anomaly over the subtropical Indian Ocean from 30° to 40°S (Figure 10b) in the NH simulation indicating that the Mascarene High extends southwestward compared to the H simulations. The meridional dipole of anomalous geopotential height produces a convergence zone centered at (80°E and 20°S), associated with the higher rainfall rate in the NH simulation. The meridional dipole pattern is also presented in the MMS cross section (Figure 10f), where two anomalous circulations emerge from the equator to 30°S with    Figure 9f) and increased precipitation (Figure 8d).
The western Pacific is characterized by the West Pacific Warm Pool and associated low-level wind convergence (Figure 10a). The equatorial easterly/southeasterly trade winds over the central Pacific (120°-170°W and 20°S-10°N) is enhanced in the NH simulation, associated with the increased precipitation (Figure 8d). The anomalous MMS displays upward vertical speeds over the central Pacific (Figure 10d). On the other side, an anomalous downward vertical motion exist near 20°S over the Pacific corresponding to the less rainfall in the NH simulation. Notably, there are middle-latitude waves over 40°-60°S accompanied by local changes in circulation (Figure 10b), although they have little impact on the precipitation due to the lack of moisture supply.

Rising Bubble Experiments
Given that the real-world simulations in Section 3 suggest that the most significant NH-H precipitation differences in E3SM are connected to regions of convective activity, in this section we employ idealized rising thermal bubble experiments to understand the essential distinction between NH and H simulations. Specifically, this section will focus on the questions: What are the differences between the NH and H rising bubble experiments, especially in the vertical velocity and precipitation? How does this connect to the real-world simulations? How does the HOMME dynamical core behavior compared to other theoretical studies (Herrington & Reed, 2018;Jeevanjee & Romps, 2016;Kato & Saito, 1995;Orlanski, 1981)? In Sections 4.1 and 4.2, we discuss experiments without and with moisture, respectively.

Experiments Without Moisture
We begin with a set of idealized rising thermal bubbles using the doubly periodic HOMME-NH and H dynamical cores without moist processes (i.e., specific humidity is zero), analogous to the dry framework of Herrington and Reed (2018). This allows us to identify a resolution threshold where the NH and H models diverge in the simplest framework. Although there is no precipitation in the dry experiments, vertical velocity can be used in the comparison.
The initial state consists of a uniform potential temperature of θ = 300 K. A standard surface pressure of p 0 = 1000 hPa is assumed. A warm bubble perturbation is applied to the θ field and is given by where Δθ = 2 K and where x r is the horizontal bubble radius and z r is the vertical bubble radius. The warm bubble is initially located at (x c , z c ) = (0 m, 2,000 m) with z r = 500 m. The domain width is fixed as 1,200 km horizontally and 20 km vertically. The average horizontal grid spacing Δx is computed as where L = 1,200 km is the domain width and n x is the number of grid points in the x direction. In the experiments, the horizontal grid resolution is controlled by changing n x since the domain width L is fixed. Table 6 provides a list of horizontal grid spacing (Δx), dynamics timestep (Δt), and horizontal bubble radius (x r ) used in the dry experiments without moist processes. To compare NH and H dynamical cores and their differences at various horizontal grid spacings, we conduct the simulations at resolutions ranging from ne4 (Δx = 100 km) to ne512 (Δx = 0.781 km). The horizontal bubble radii and dynamics timesteps are scaled to the resolution. Another series of simulations using a constant dynamics timestep (i.e., Δt = 0.781 s) among all resolutions are performed to test the sensitivity to the dynamics timesteps (not shown). These results are not sensitive to the dynamics timestep-differences between the simulations using scaled and uniform dynamics timesteps are smaller than 0.1%. Figure 11a shows the 900 hPa vertical pressure velocity ω, at 3 hr into the simulations as a function of horizontal grid spacing, while Figure 11b shows the time series of 900 hPa vertical pressure velocity ω simulated in the NH (solid) and H (dashed) simulations. The ω shown is averaged underneath the bubble over (−x r and x r ). The magnitude of the H vertical velocity is inversely proportional to the grid spacing (Figure 11a), consistent with past studies (e.g., Herrington & Reed, 2018;Jeevanjee & Romps, 2016), while the NH vertical velocity plateaus at finer resolutions. Differences in ω between NH and H simulations are negligible (<10%) from 100 to 12.5 km resolutions (Figure 11a). However, ω in the H simulation is significantly greater than its NH counterpart at resolutions finer than 6.25 km, indicating an overestimate of the upward vertical motion in the H dynamical core (Figure 11a). A strong upward vertical velocity develops abruptly in the H dynamical core at the beginning of the simulation (0-10 min), reflecting the model's attempt to rapidly recover hydrostatic balance (Figure 11b). The threshold where the NH and H difference emerges (Δx ∼ 6.25 km) is close to the value observed in past studies (Orlanski, 1981;Kato & Saito, 1995), although the conclusions about the specific threshold in the literature vary as discussed in Section 1. For example, Jeevanjee (2017) shows the differences between hydrostatic and NH solvers are virtually undetectable at resolutions of 2 km or coarser (which is about two times finer than our result) and indicates that this number is almost certainly dependent on the model and the phenomena under consideration. In  our dry experiments at fine resolutions (Δx < 1.56 km), the H dynamical core overestimates the ω by a factor of 2-3, consistent with past studies (Jeevanjee, 2017;Weisman et al., 1997). While this "dry" threshold of 6.25 km cannot be used to explain the NH and H differences observed in Section 3, which are observed at much coarser resolution, it is nonetheless illuminating to investigate the drivers of this difference further.
In the HOMME dynamical core, vertical velocity w advection (−u ⋅∇w) occurs in the NH model while it is absent in the H model. One hypothesis is that the differences in the simulated vertical velocities are caused by the presence of w advection. To test the hypothesis, we conduct another series of dry experiments similar to those in Table 6 but with the w advection term turned off by setting the term as zero in the source code. The differences between simulations with and without w advection are negligible (<0.25%) at all horizontal resolutions from 100 to 0.78 km, suggesting w advection is not responsible for the different vertical motions between NH and H simulations. We also test this hypothesis in the real-world simulations (as in Section 3) by comparing the results with and without w advection term. The realistic simulations also confirm that the w advection is not the reason for the different behaviors between NH and H simulations.
Instead, the overestimation of vertical velocity by the hydrostatic solver is better explained using the theory of effective buoyancy illustrated in Figure 6 of Jeevanjee (2017). In the NH model, the total pressure (p in Equation 3) has two components: the hydrostatic pressure p h and a residual NH pressure component ′ , that is, where the prime is used to help distinguish the total pressure in the NH model (p) and the NH pressure perturbation ( ′ ) . The effective buoyancy, β, is defined as in both NH (β nh ) and H (β h ) systems, and represents the net vertical acceleration due to density anomalies. By setting the wind field u = 0, the effective buoyancy neglects the "dynamic" pressure gradients arising from advection of momentum (Jeevanjee & Romps, 2016;Markowski & Richardson, 2011). The effective buoyancy pressure, p β , is defined as which is the NH pressure field associated with zeroing out the wind field. In the case of a positively buoyant bubble, a negative horizontal gradient of hydrostatic pressure (∇ h p h < 0) is produced forcing the mass convergence into the bubble. In the hydrostatic systems where ′ = 0 , mass continuity is enforced by demanding that vertical divergences cancel all horizontal convergences. The vertical acceleration β h must compensate for the convergence from the horizontal gradient of ∇ h p h . In contrast, in the NH systems, the NH p β generates both horizontal and vertical divergence. Therefore, the divergence from β nh needs only to compensate for the net convergence in the horizontal direction, which exhibits a partial cancellation between ∇ h p h and ∇ h p β , resulting in β h > β nh .
Furthermore, ω using NH dynamical core is found to converge to a constant value (∼−0.45 Pa/s) at the resolution of 1.56 km or finer while the H analog does not converge (Figure 11a). This grid spacing of convergence is consistent with the theoretical threshold from the framework of convecting parcels in Jeevanjee (2017). Namely, we would expect the NH ω to converge to a constant realistic value when Δ where H = 1,000 m is the bubble height and n gray = 2x r /Δx = 1. Equation 15 yields 1,000 m using parameters in our case, close to the convergence threshold of 1.56 km observed in Figure 11a.

Experiments With a Simple Condensation Scheme
As our dry bubble experiments suggest that the H and NH models should diverge at grid spacing of 6.25 km, it's clear that dry processes alone are insufficient to explain the differences observed in our realistic simulations, especially over the moist tropics. Consequently, we now turn our attention to rising bubble experiments with moisture included. To simplify the moist processes, the dynamical cores are coupled to the large-scale condensation scheme implemented in the "simple physics" package of Reed and Jablonowski (2012), which assumes all condensates are immediately removed from the column with no cloud phase. An implementation of this scheme is available as part of the DCMIP2016 test case document (Ullrich, Jablonowski, et al., 2016). Simulations were also attempted with Kessler-type microphysics, but with the complexity of this scheme it was difficult to determine an idealized setup that exhibited consistent results.
In the moist experiments, the potential temperature is initialized using the same setting as in the dry experiments (Equation 9). Relative humidity is constant, except for a relative humidity perturbation collocated with the potential temperature perturbation of the form where Δh is the relative humidity perturbation maximum and R is a normalized distance from the center of the bubble, as defined in Equation 10. Two combinations (abbreviated as RH1 and RH2) of background relative humidity h 0 (i.e., the relative humidity outside the bubble where R > 1) and relative humidity perturbation maximum Δh are employed in the moist experiments as shown in Table 7. They are designed to mimic the relative humidity over the tropics in realistic simulations and to test the sensitivity of the results to the moisture initialization. The saturation specific humidity q s is given by Equation 13 of Reed and Jablonowski (2012). The specific humidity q in the initialization is calculated as where q s is the saturation specific humidity, RH′ is the relative humidity perturbation given in Equation (16), and h 0 is the background relative humidity.   & Reed, 2018;Kato, 1997), we investigate multiple bubble radii ranging from 0.5Δx-3Δx at each resolution to test the sensitivity of the results to the bubble widths. Figure 12 shows the time series of 900 hPa pressure vertical velocity ω, precipitation, and 250 hPa temperature averaged over the spatial footprint of the bubble (−x r and x r ) using various bubble radii and two humidity initializations (Table 7) at horizontal grid spacings of 12.5 and 25 km. The time series of precipitation and 900 hPa vertical velocity are highly correlated. As the bubble rises vertically from the initialized location, the precipitation rate increases associated with upward vertical velocity. The rising speeds of the bubbles vary depending on the bubble sizes and humidity initializations. Following the peaking of precipitation, the bubble develops near the top of the atmosphere, leading to subsidence on either side of the center, while the upper-level temperature increases.

Differences in Vertical Velocity and Precipitation
In contrast to the dry experiments (Section 4.1) where no notable difference is found in vertical velocity between NH and H dynamical cores at 25 km resolution (Figure 11), the ω and precipitation differences are detectable even at 25 km resolution in the moist experiment. Over the period prior to peak, the upward vertical velocity using the H dynamical core is larger than its NH counterpart by up to 10% at 25 km resolution and 65% at 12.5 km resolution (Figures 12a and 12d). The H simulation is also associated with higher upper-level temperature than the NH simulation, consistent with more latent heat release. The exaggerated upward motion in the H Figure 12. Time series of (a) 900 hPa pressure vertical velocity ω (Pa/s), (b) precipitation (mm/day), and (c) 250 hPa temperature (K) averaged underneath the bubble over (−x r and x r ) using the NH (solid) and H (dashed) dynamical cores at 25 km resolution. Panels (d-f) are the same as panels (a-b), but from the simulations at 12.5 km resolution. RH1 and RH2 are defined in Table 7. simulation is consistent with the theory of effective buoyancy described in Jeevanjee (2017). The comparison of moist and dry experiments indicates that the difference between NH and H dynamical cores is amplified by the condensational heating in the moist physics. The conclusion is similar to Yang et al. (2017), which compared results between moist and pseudodry simulations and found much smaller NH-H difference in the pseudodry simulation. Our real-world simulations (Section 3) at 28 km resolution also show larger NH-H difference over the tropics in the presence of adequate moisture suggesting that physics-dynamics coupling particularly related to moist physics is a key reason for the significant NH-H difference at coarse grid spacing (e.g., 25-28 km).
To quantify the precipitation difference, Table 9 presents three precipitation-related indices among various combinations of bubble radius and humidity initialization in 25 km simulations. In the first and second indices, the ratio of accumulated precipitation in H simulation compared to that in NH simulation is calculated over two accumulation time periods: one for the entire time series shown in Figure 12b (Period 1; 0-270 min) and the other for the period preceding H's precipitation maximum (Period 2). Given that the bubbles rise at different speeds depending on the bubble sizes and humidity initializations, the absolute length of Period 2 is subject to the variations. The third index measures the ratio of precipitation rate maximum in the H simulation compared to its NH counterpart. Table 10 is the same as Table 9, but from simulations at 12.5 km resolution (Figure 12e).
NH and H simulations generate approximately same amount of precipitation (with difference smaller than 2%) over the entire time series (Period 1). This appears consistent with our observations from the realistic simulations, where total global precipitation was essentially the same in both dynamical cores. However, the precipitation maximum in the H simulation is higher than NH by 3%-4% at 25 km and 11%-18% at 12.5 km. In Period 2, accumulated precipitation in the H simulation exceeds that in NH by 8%-21% at 25 km and 36%-106% at 12.5 km. Differences in precipitation between NH and H become more pronounced when resolution becomes finer. The accumulated precipitation ratio in our simulation at 25 km resolution agrees roughly with Kato (1997), who found ratios ranging from 1.06 to 1.25 at 20 km resolution depending on the presence of hydrostatic water loading but applied only one bubble configuration (vertical bubble radius z r = 1.5 km and horizontal bubble radius x r = 8Δx). Figure 13 shows time series of μ − 1 at 900 hPa averaged underneath the bubbles in the moist experiments with simple physics at 25 and 12.5 km resolutions. Only NH simulations are compared as μ is a constant (i.e., 1) in all H simulations. According to the definition of μ in Equation (3), the deviation of μ from 1 is caused by the inequality of NH and hydrostatic pressures. The evolution of μ − 1 is strongly correlated with the precipitation and vertical velocity ( Figure 12). The magnitude of μ − 1 ranges from −0.0004 to 0.001 at 25 km resolution and becomes approximately two times larger at 12.5 km resolution. In contrast, μ in the dry experiment is closer to 1 than in the moist experiment, indicating that enhanced nonhydrostaticity is caused in the presence of moisture. Specifically, the deviation of μ from 1 (i.e., μ − 1) is on the scale of 10 −5 at 25 km resolution in the dry experiment. As the grid spacings become finer, μ − 1 has a wider range indicating that the importance of NH pressure significantly increases as the grid spacing decreases.

Nonhydrostatic Pressure
Although the magnitudes in the precipitation amounts show some parallels, there are differences between the idealized rising bubble experiments and real-world simulations. These differences appear rooted in the idealized bubble experiment being more illustrative of short term effects. The thermal bubble rises vertically without horizontal propagation, whereas in the realistic global simulations, enhanced or suppressed convection in one region could influence the surrounding regions by modifying the atmospheric circulations and lead to the spatial patterns over hotspots.
The rising bubble experiments do suggest that the inclusion of moist processes can lead to significant divergence between the NH and H models even at resolutions of 25 km, and that the 10%-20% deviation in total precipitation   Table 10 Same as Table 9, But From the Simulations at 12.5 km Resolution amount observed in the realistic experiments is indeed reasonable in that case. With that said, these bubble experiments do not provide a completely satisfactory explanation for our observations in the realistic experiments: Namely, in the realistic experiments we observe that the NH model actually tends to develop stronger convective plumes in some regions, contrary to the theory in the bubble experiments and past studies (Section 4.1; Yang et al., 2017). This suggests that further work is needed to understand the differences between H and NH dynamics on seasonal to climatological time scales.
An understanding of the influence of NH terms on the climatology is complicated by sensitivities in the model to small differences in forcing. The degree to which these terms influence the results depends strongly on the configuration of the solver, the coupling to physics, implementation of topography and parameterizations, and the problem and region being examined. For example, Dueben et al. (2020) compared NH and H forecast simulations at 1.45 km resolution and showed the difference between H and NH simulations is smaller when compared with other changes in the model configuration (e.g., topography and timesteps). Our results do not contradict this observation and simply suggest that model sensitivities from moist physics will amplify the effect of the forcing terms in the NH formulation, and can result in significant differences in the climatological variables that are far larger than dry NH-H studies would suggest. Notably, the specific differences in spatial patterns that emerge in our realistic simulations are indeed sensitive to the physics timesteps, although the regions of greatest sensitivity are largely the same. In a study of TC forecasts, Chen et al. (2019) found the initial condition is more important than the choice of the NH/H dynamical core. Hence, the model configurations must be taken into account and the NH effect should be further investigated in future studies with different modeling frameworks.

Conclusions
This study explores the differences between NH and hydrostatic (H) models using seasonal SCREAM simulations with the HOMME dynamical core. We employ ensembles of realistic global simulations at seasonal timescales in both boreal summer and winter seasons at ne120 (28 km) and ne256 (13 km) resolutions. Significant differences in precipitation amount (on the order of 10%-20%) and associated precipitation patterns are observed in some regions, and for features with strong convections such as TCs. These differences appear robust across ensemble members and so cannot be attributed to internal variability. Results from idealized dry and moist rising bubble experiments are used to validate our observations in the realistic simulations; while the dry bubble test cases agree well with theory and show minimal deviation between H and NH dynamics until grid spacing is below 6.25 km, precipitation amount and vertical velocity are different in the moist case even with a grid spacing of 25 km.
Notably, the differences between NH and H dynamics are not significant when purely examining global statistics at 28 km resolution. The precipitation difference is 10%-20% over the hotspots in the tropical oceans, 10%-20% for TCs, and 5%-10% over the topographic features such as the mountainous Western U.S. (Tables 2  and 5). Notably, in highly convective regions, it is the NH model that is wetter on average, which is contrary to Yang et al. (2017) and our expectations given that NH models tend to produce weaker vertical motion (Jeevanjee, 2017). Nonetheless, the differences in precipitation totals are generally in agreement with Yang et al. (2017) (i.e., up to 7% at 36 km and 15%-34% at 12 km resolutions over some areas in the aquaplanet simulations with idealized mountains). The most prominent differences emerge near India and the Western Pacific in the boreal summer (Figure 1), and central-southern Indian Ocean and Pacific the in the boreal winter ( Figure 8). Although past studies Janjic et al., 2001;Lebassi-Habtezion & Diffenbaugh, 2013) showed the NH effects are more influential over the topography, the tropical oceans near the ITCZ where heavy rainfall exists are more likely influenced by the NH dynamics and their interactions with latent heating.
The differences between NH and H dynamical cores over a specific area influence the surrounding regions by modifying surrounding atmospheric circulation. These differences are particularly pronounced when moist processes are present. Differences in the temperature and geopotential height emerge over middle to high latitudes related to planetary-scale atmospheric Rossby waves (Figures 4 and 9). The observed differences in midlatitudinal Rossby wave activity are likely driven by the tropical circulation differences, and corresponding propagation of this signal to the extratropics, via a Gill-type response (Gill, 1980). The patterns of precipitation differences are related to the modifications to the local circulations, especially in the meridional direction ( Figures 5 and 10).
In the idealized rising thermal bubbles experiments without moist processes, differences in vertical velocity between NH and H simulations are negligible (<10%) from 100 to 12.5 km resolutions ( Figure 11). However, a significant overestimate of the upward vertical motion by the H dynamical core is displayed at resolutions finer than 6.25 km, close to the past theoretical studies (Kato & Saito, 1995;Orlanski, 1981), which could be explained using the theory of effective buoyancy illustrated in Jeevanjee (2017). Vertical velocity advection is found to be unimportant in both idealized experiments and realistic simulations.
In the moist rising bubble experiments, the dynamical cores are coupled to the large-scale condensation scheme implemented in the "simple physics" package of Reed and Jablonowski (2012). In contrast to the dry experiments where no notable difference is found in vertical velocity between NH and H dynamical cores at 25 km resolution (Figure 11), the vertical velocity difference is detectable (up to ∼10%) even at 25 km resolution in the moist experiments, indicating that the difference is amplified by the latent heating in the moist physics. The difference between NH and H dynamical cores in accumulated precipitation is found to be 8%-21% at 25 km resolution and 36%-106% at 12.5 km. Comparing μ (Equation 3) at the same resolution between dry and moist experiments, the deviation of μ from 1 in the moist experiment is 10 times stronger than the dry experiment, illustrating that the stronger deviations of μ are caused in the presence of simple physics.
In summary, our study demonstrates that the difference between NH and hydrostatic dynamical cores, when moist processes are included, is significant at resolutions as coarse as 28 km. This grid spacing is much coarser than the 5-10 km threshold obtained from our dry rising bubble experiments and past theoretical idealized analysis (Kato & Saito, 1995;Orlanski, 1981). For the cases where tropical precipitation near the ITCZ is of interest, a NH dynamical core may be needed in climate simulations even at coarse grid spacing.