Adapting the MPAS Dynamical Core for Applications Extending Into the Thermosphere

To extend the nonhydrostatic global Model for Prediction Across Scales (MPAS) for deep‐atmosphere (geospace) applications, we have modified the model equations and numerics to include variable atmospheric composition and (potentially large) molecular viscosity and thermal conductivity. The split‐explicit numerical integration techniques in MPAS remain stable in idealized test cases for atmospheric domains extending into the upper thermosphere and continue to provide an efficient numerical framework for nonhydrostatic simulations. Variations in the atmospheric constituents influence the dynamical equations by altering the heat capacity and ideal gas constants. These feedbacks require little alteration of the dynamical equations although our testing reveals that the amplitude of disturbances may be sensitive to even small variations in the thermodynamic coefficients. Although the potential temperature is no longer formally conserved for adiabatic flow, it remains effective as a prognostic thermodynamic variable in the model equations. Molecular viscosity and thermal conductivity are dominant influences in the upper thermosphere and are represented implicitly in the model numerics. Because of the large magnitude of these terms, their treatment, though stable, may significantly under represent the true magnitude of their damping effects. Further consideration of these deep‐atmosphere extensions to MPAS will be explored in more realistic simulations of thermospheric dynamics.

In recent modifications to the MPAS dynamical core, Skamarock et al. (2020) have removed the shallow-atmosphere approximation and adapted the model numerics for the traditional deep-atmosphere equations (White et al., 2005) in which the depth of the atmosphere is not negligible in comparison to the radius of the earth. For this purpose, the MPAS equations and numerics were modified to use the actual geocentric distance r instead of the Earth radius e r in the governing equations and the grid mesh configuration, to allow gravity to vary with height, and to include the full Coriolis force terms containing vertical velocity components in the momentum equations. With these augmentations to the dynamical core, the MPAS solver remains valid for applications extending into the lower thermosphere (100-200 km).
Further extension of MPAS for use throughout the thermosphere raises additional challenges. The model numerics must accommodate extreme variations in the atmospheric parameters as the model domain may span 25-30 scale heights, with the density decreasing by 12 orders of magnitude, and thermospheric temperatures that often reach 1,000-2,000 K. The numerical approach for integrating the prognostic equations in the current MPAS utilizes a split-explicit technique for integrating the prognostic equations in which the terms responsible for the propagation of acoustic and gravity-wave modes are advanced on smaller time steps (limited by the horizontal speed of sound), while the remaining terms and parameterized physics are stepped forward with a larger time increment appropriate for the physical modes of interest (Klemp et al., 2007;Skamarock et al., 2012). Thus, a basic question arises as to whether the split-explicit numerics used in the model integration remain viable when applied throughout the upper atmosphere.
In addition to issues related to the basic model numerics, a number of additional modifications to the MPAS dynamical equations and numerics are required to render the model potentially suitable for geospace applications. In the thermosphere, the composition of the atmosphere can no longer be considered independent of altitude; atmospheric constituents vary with scale heights inversely proportional to their molecular weight, and atomic ions become increasingly dominant. Changes in the atmospheric composition alter the thermodynamic properties such as heat capacity and the ideal gas constant, which also require modification of the thermodynamic equation. Although the kinematic viscosity and thermal diffusivity are negligible in the troposphere and stratosphere, these effects become dominant in the upper thermosphere.
To evaluate the viability of the MPAS model numerics with the inclusion of these deep-atmosphere modifications, herein we conduct our initial assessment in a 2-D vertical slab configuration on a flat earth before implementing these extensions in the full 3-D global MPAS code. Our emphasis is on assessing how this increased complexity can be accommodated while maintaining the accuracy and stability of the overall MPAS numerical framework. For this reason, we do not consider variations in the grid configuration with altitude, which have already been implemented in the full 3-D MPAS code (Skamarock et al., 2020). Confining our focus to the numerics for the basic dynamical system, we also defer any treatment of the added complexity of upper atmosphere chemistry and physics such as solar and Joule heating, ion drag, photoionization/ photodisassociation of molecules, etc.
We begin in Section 2 with a description of modifications of the dynamical equations to include variable atmospheric composition, molecular viscosity and thermal conductivity, and to improve the treatment of pressure gradients in deep-atmosphere applications. In Section 3, we explain how our split-explicit numerical integration technique is adapted to work effectively in deep-atmosphere environments, and we describe our initial method for including implicit treatment of the molecular viscosity and thermal conductivity terms. Testing of the model numerics is presented in Section 4 for idealized simulations in a 2-D slab configuration for an elevated oscillating heat source and for mountain wave simulations at several different horizontal scales. Summary remarks are included in Section 5. numerical error, negative pressures can also arise if, for example, the hydrostatic equation is integrated vertically from the surface to balance the initial thermodynamic fields. For these reasons, we have altered the representation of pressure gradients in the dynamical equations, writing them in terms of   ln( / ) p p 0 , which tends to vary more linearly with height and avoids the possibility of producing negative pressures. With this modification the vertical pressure gradient and buoyancy terms in the vertical momentum equation are represented simply as where d M is the molecular weight of the dry air mixture, and  * 8.314 R J/(Kmol) is the universal gas constant (product of the Boltzmann constant and the Avogadro constant). Here, the subscript j refers to the th j atmospheric constituent having a mixing ratio  j , molecular weight j M , and j f molecular degrees of freedom. For the purpose of including the coupling of the dynamics with composition, we assume the atmosphere is comprised of nitrogen ( 2 N ), oxygen ( 2 O ), and free oxygen (O), which are the three major neutral constituents of the thermosphere (Dickinson et al., 1984). Since the mixing ratio for a constituent is defined as the density of the constituent divided by the density of dry air,   

Thermodynamic Equation
The potential temperature coupled with the dry-air density is utilized as a prognostic variable in MPAS as it is a conserved quantity for adiabatic flow in a shallow atmosphere with constant composition. This KLEMP AND SKAMAROCK 10.1029/2021MS002499 3 of 19 representation also facilitates the formation of the vertically implicit numerics that is an essential aspect of the split-explicit numerical scheme (Klemp et al., 2007). However, in the deep atmosphere, as discussed above, the atmospheric composition is not constant, which alters the conservative properties of this formulation. This is illustrated by writing the internal energy equation (5) in terms of the potential temperature where Q is the specific heating and   / d p R c . In the absence of diabatic heating (  0 Q ), potential temperature is not conserved since  will vary due to variations in the atmospheric composition. This correction to the potential temperature equation has also been proposed for the WACCM-X 2.0 model for extension of the finite volume (FV) dycore to the deep atmosphere (Liu et al., 2018). Technically, the thermodynamic coefficients in  are not precisely constant even in the shallow atmosphere equations in the presence of moisture. However, these variations in  are typically less than a percent and are usually ignored (Emanuel, 1994).
Since the atmospheric composition varies dramatically across the deep atmosphere, we retain this term to further assess its significance.

Deep Atmosphere Equation Set
The prognostic equations in MPAS for a shallow atmosphere are written in flux form using a height-based terrain-following vertical coordinate    ( , ) H x z as presented by Skamarock et al. (2012). As mentioned above, we do not consider here vertical variations of the grid-cell area that would arise in a deep atmosphere equations discretization on the globe. These metric adjustments have been derived and implemented in the full global MPAS code, as described by Skamarock et al. (2020). Flux variables are thus defined here using a transformation-adjusted dry density      d d z , such that     d d represents the mass in a grid cell. In terms of the flux variables the nonhydrostatic equations can then be expressed as Here, the notation for differentiation follows Dutton (1986) as summarized by Skamarock et al. (2012). The horizontal momentum Equation 8 is written in vector invariant form (Skamarock et al., 2012;Thuburn et al., 2009)  Equations 8-13 are the same as the deep atmosphere equations in Skamarock et al. (2020) except for the modified form of the pressure gradient terms in 8 and 9, the more general form of the potential temperature Equation 10 that permits  to vary in response to changes in the atmospheric composition, and the addition of prognostic equations for the individual atmospheric constituents 13. To be consistent with the flux-form transport of other prognostic variables, the  / d dt term from 6 is also written in flux form in 10, which is the same form as suggested by Liu et al. (2018) for their extensions of the FV dycore. Since the influences of molecular viscosity  and thermal conductivity c k may become dominant in the deep atmosphere, diffusive terms for  and c k are also included in the momentum 8 and 9 and thermodynamic 10 equations, respectively.
The above equations include a number of simplifying approximations. Here, the viscosity terms in 8 and 9 have been approximated using the full velocity gradient (as in its incompressible form) instead of the symmetric stress tensor (that is appropriate for compressible flow). Because the influences on viscosity become large in the upper thermosphere, the viscosity terms must be solved implicitly (discussed further below) and thus using the full stress tensor would add significant complexity to the numerics. However, this simplification has been shown to have a noticeable effect in global climate simulations (Becker, 2001;Becker & Burkhardt, 2007). Although we have added dissipation terms in Equations 8-10 and 13, we have not included dissipative heating terms in the thermodynamic Equation 10 that arise from those terms. In addition, the atmospheric composition is represented with the simplifying assumptions that all of the constituents are in mechanical and thermal equilibrium, and thus require computing budgets only for the mass of each constituent [as in 13] without the need for constituent budgets for momentum and entropy. The simplifications mentioned above are consistent with representations in existing deep-atmosphere models [for example, the hydrostatic WACCM-X model (Liu et al., 2010) and its predecessors (Dickinson et al., 1984) and the non hydrostatic GITM model (Ridley et al., 2006),] but will require further consideration in our future work to enhance the accuracy of MPAS for deep-atmosphere applications.
The equations for the composition of the atmospheric constituents 13 follow those derived and applied in thermospheric models by Dickinson and Ridley (1972), Dickinson et al. (1984) and R. G. Roble et al. (1987). The first three terms on the RHS of 13 represent the advective transport, the molecular diffusive transport, and the vertical eddy mixing, respectively. In the absence of turbulent mixing, the molecular constituents tend toward an equilibrium state (called aerostatic balance) in which the partial log pressure of each constituent  i varies vertically in proportion to its molecular weight i M , such that  The molecular diffusive transport term in 13 governs the rate at which the atmosphere relaxes toward this aerostatic balance. In this term,  is a diagonal matrix that represents the relative diffusive rates of the constituents (Dickinson et al., 1984) and  j L , given by reflects the departures from diffusive equilibrium, such that   0 j L corresponds to aerostatic balance. Here,  thermopause (∼85 km), the eddy coefficient  K is expected to decrease with height such that it has little influence above the middle mesosphere.
As mentioned above, we consider here only the three major atmospheric constituents that comprise the thermosphere, 2 N , 2 O , and O. Since the sum of the mixing ratios  j must equal unity, we solve the prognostic To illustrate the effect of the eddy diffusivity on the atmospheric constituents, we have integrated 13 for an atmosphere at rest, using a horizontally homogeneous vertical temperature profile adapted from a representative sounding from the WACCM-X 2.0 Model (Liu et al., 2018) as shown in Figure 1. The model top is set at 500 km with   2 z km. The three atmospheric constituents  Figure 2. Notice that in the initial aerostatically balanced state the mixing ratio of molecular oxygen decreases rapidly with height in the lower atmosphere while the nitrogen correspondingly increases.
This occurs because the density scale height for 2 O is smaller than that for 2 N due to its greater molecular weight. This behavior is obviously unrealistic and is removed by the vertical eddy mixing, which promotes well mixed constituents below about 100 km. The resulting vertical distribution of constituents in this simple idealization is qualitatively similar to the global mean mixing ratios obtained in the NCAR TCGM model for the thermosphere (Dickinson et al., 1984; Figure 4). Variations in the atmospheric composition may contain considerable additional complexity through the source S and removal R terms in 13 that govern the disassociation and recombination of the various atmospheric constituents. For our purposes in evaluating the numerics of the dynamical equations, we ignore these effects and include only the transport and diffusive processes for these constituents. Although changes in composition may have significant influences in regulating the behavior of the deep atmosphere, it is worth noting that the only feedbacks of atmospheric composition to the dynamical Equations 8-11 occur through variations in the thermodynamic coefficients d R and p c as defined in 3 and 4.

Numerical Integration of the MPAS Deep-Atmosphere Equations
The split-explicit numerical integration of the dynamical equations in the shallow-atmosphere MPAS is a fundamental component of the model numerics, and is thus a critical aspect to evaluate for deep-atmosphere configurations. In the split explicit formulation, the pressure gradient and buoyancy terms in 8 and 9 are stepped forward on the smaller acoustic time steps, along with the flux divergence terms in the potential temperature and density Equations 10 and 11. In evaluating these terms, the horizontal gradients are computed using explicit (forward-backward) numerics. However, using explicit numerics in the vertical direction would KLEMP AND SKAMAROCK 10.1029/2021MS002499 6 of 19  typically impose severe restrictions on these acoustic time steps. Therefore, numerical efficiency is achieved by writing the vertical pressure gradient, buoyancy, and vertical flux divergence terms implicitly and solving within each vertical column by inverting a simple tridiagonal matrix. (See Klemp et al., 2007 for further details.) Thus the restriction on the acoustic time step should be limited only by the horizontal sound wave propagation on the horizontal grid. On these acoustic time steps, the prognostic variables are cast as perturbations from the most recent full time step. In this manner, essentially all of the slow-mode contributions to these terms are evaluated on the large time steps using a third order Runge-Kutta integration scheme. With these numerics, the model equations can be integrated on large time steps comparable to hydrostatic thermospheric models such as the TGCM (Dickinson et al., 1981) and WACCM-X (Liu et al., 2018), which use time steps on the order of 5 min, in comparison to the nonhydrostatic GITM that is limited to a 2 s time step due to its fully explicit time integration (Ridley et al., 2006).
Although we have written the pressure gradient terms in 8 and 9 in terms of   0 ln( / ) p p as discussed in Section 2.1 to improve their numerical accuracy, this representation is applied only to the portion of the terms advanced on the large time steps. Thus, the pressure gradient terms in the small time step equations do not need to be altered from their form in the shallow-atmosphere implementation; they are the same as in Equations 13-16 presented in Klemp et al. (2007) except that reference profile variables in the vertical momentum Equation 14 are replaced by the full variables at time t. Since the small time step variables are expressed as perturbations from the most recent large time step, the specific numerics used in representing the small time step terms have little impact on the overall accuracy of the numerical integration. Additional terms that have been included for deep atmosphere applications are accommodated in the large time step calculations.
The density-coupled mixing ratios for the atmospheric constituents 13 are advanced in flux form over the Runge-Kutta time steps in a similar manner to the moist constituents 12. After completing the Runge-Kutta steps, d R and p c are updated to the new time level using 3 and 4, respectively. The variability in   / p R c directly influences the potential temperature Θ m through the second term on the rhs of 10. This is evaluated following the update of the thermodynamic coefficients using the expression and applied as a forcing term in the next time step for Θ m . Here, the overbar indicates an average between times t and   t t, and a V is the momentum that has been time-averaged over the acoustic steps, which is also used to advect the density weighted scalars to maintain consistency with the numerics of the continuity equation. Although the term in 15 is quite small, the integration is sensitive to the specific numerical representation of this term. For example, in the flux divergence term in 15, representing  at either time t or   t t (instead of averaging in time) can occasionally promote an numerical instability in the numerical integration.
The viscosity and diffusivity terms in 8-10 require some special treatment. Because the molecular viscosity and thermal conductivity are relatively constant throughout the atmosphere, the kinematic viscosity k c increase exponentially with height in inverse proportion to the density, such that these terms become dominant in their influences on v and  m the upper thermosphere. Because of their large magnitude, these terms must be computed implicitly to maintain stability. At present, this is accommodated with incremental updates through directional splitting. At the beginning of each large time step the diffusion terms are evaluated implicitly, first in the vertical and then in the horizontal. The tendencies from these updates are then combined with the tendencies from the remaining terms in the equations in completing the time step.
This approach is easy to implement and produces the desired effect of strong damping in the upper thermosphere. However, it raises the question as to whether these terms are being treated with sufficient accuracy. Setting aside the issue of evaluating the terms with incremental updates, consider the influence of a basic one-dimensional damping term ), it is much less than the correct analytic value given by  . The deficiencies of Crank-Nicolson for large values of  are discussed by Osterby (2003), who also proposes several alternatives besides forward centering to reduce Crank-Nicolson oscillations. If more accurate treatment of the molecular viscosity and thermal conductivity is needed in thermospheric applications, further measures, such as temporal substepping, could be considered.
Because atmospheric quantities such as pressure and density may vary by many orders of magnitude across a deep-atmosphere domain, we expect that numerically integrating the equation set 8-13 will require double precision calculations. In practice, we have found that double precision is needed when the model top is significantly greater than about 200 km to avoid the appearance of small scale noise in the model fields.
In the current implementation of MPAS for shallow atmospheres, the thermodynamic variables in the governing equations are cast in terms of perturbations ( d  ,  p , Θ m  ) from a representative reference sounding ( d , p, Θ m ) to reduce truncation and roundoff errors in the numerics (Klemp et al., 2007). However, for deep-atmosphere applications, the perturbations variables might typically become large in comparison to their reference values, thus negating much of the motivation for their implementation (and further obviated by the use of double precision). Thus, double precision and the discrete Equations 8-13 employing full variables are used in obtaining the results reported in this paper.

Elevated Heating/Cooling Function
For an initial test of the model numerics, we simulate the 2-D atmospheric response in a non-rotating framework (  0 f ) to an oscillating heating/cooling source term located in the middle thermosphere. The model domain has a depth  500 t z km and is initialized using a horizontally homogeneous sounding with the temperature profile shown in Figure 1 and where the heating rate is  0.015 h C K 1 s and the radial nondimensional distance q L from the center of the heating function at ( , c c x z ) is given by The heating function is centered at  300 c z km, with c x located in the middle of the model domain. For this heating function, the thermal forcing oscillates with a period of 2 h.
Because the influences of the molecular viscosity  and thermal conductivity c k become large in the thermosphere, we tested the model numerics with several configurations for the model dissipation. We first simulated the case with    0 c k to observe the inviscid response. Since this case will produce significant reflection of gravity waves from the rigid lid upper boundary, we also ran the test case with    0 c k but included a gravity wave absorbing layer near the model top. The absorbing layer is specified as proposed by (Klemp et al., 2008) . For the final test configuration, we ran with no absorbing layer but specified  and c k at their appropriate atmospheric values as described by Liu et al. (2010)  km due to the large decreases in density, confirming the necessity to use implicit numerics. The same considerations apply to the mixing due to the thermal conductivity.
For each of these cases, we set the horizontal and vertical mesh spacings to   20 x km and   2 z km, respectively, and integrated the split-explicit model numerics with a large time step   80 t s and a smaller, acoustic time step     / 6 t out to 24 h. These integrations exhibit good stability with no additional explicit filters besides those specifically included as modifications in this test case. For this temperature sounding (Figure 1), the speed of sound increases by a factor of 2.8 from the surface to the top of the model domain (347986 m 1 s ), such that the horizontal acoustic Courant number at upper levels is about 0.65. In treating the atmospheric composition, we avoid the complication of spinning up the atmosphere from an aerostatic balance by initializing the constituents with vertical profiles that approximate those shown by the solid lines in Figure 2.
The time series for the maximum vertical velocity for these three test simulations are displayed in Figure 3. In each of these cases, the flow achieves a nearly steady 2 h periodicity, consistent with the thermal forcing period, by the end of the 24 h integration. As expected, the strongest response occurs in the inviscid case. The perturbations amplify with increasing height above the forcing level, and are then reflected from the rigid lid upper boundary, as is evident in the plots of the vertical velocity and perturbation temperature at  24 t h in Figures 4a and 4b. Including a gravity wave absorbing layer above  400 d z km, the maximum vertical velocity (Figure 3) closely follows the inviscid profile over the first 12 h and then rapidly levels off to a steady periodic behavior as the wave energy reaches the absorbing layer. With this absorbing layer, the vertical velocity and perturbation temperature fields in Figures 4c and 4d are very similar to those in the inviscid case in the upward propagating branch of the disturbance below 400 km.
Including the molecular viscosity and thermal conductivity terms, the perturbations produced by the thermal forcing are strongly damped. The periodic oscillations become essentially steady by 6 h, with significant amplitude only in the immediate vicinity of the forcing (Figures 4e and 4f). Notice that the contour intervals for the vertical velocity and temperature in Figures 4e and 4f are an order of magnitude smaller than in Figures 4a-4d.

Mountain Waves
For a mountain wave test case, we are faced with the complication that the wave amplitude increases with increasing height in proportion to the inverse square root of density. For the temperature sounding in Figure 1 this would result in an amplification factor of ∼ 5 5 x 10 between the surface and a height of 500 km. This amplification is largely offset in the mid to upper thermosphere by the molecular viscosity and thermal conductivity, which produce strong dissipation at these altitudes. To test the viability of the numerics for KLEMP AND SKAMAROCK 10.1029/2021MS002499 9 of 19  domains extending through the thermosphere, we first consider deep atmosphere mountain wave solutions in the absence of molecular viscosity and thermal conductivity. For this purpose, we specify an admittedly artificial vertical eddy mixing that is large enough to keep the waves below overturning amplitudes at high altitudes. This eddy mixing combines a small constant coefficient with a Richardson number i R dependent local eddy mixing that loosely follows the scheme proposed by Hong et al. (2006) for free atmosphere diffusion: Here, we set the length scale    z and define  0.2 c , which is smaller than in Hong's formulation in order activate the local eddy mixing over a broader range of Richardson numbers. The constant background eddy mixing is set at    2 / b K z t to control small scale noise in the simulation while allowing the local eddy mixing to stabilize the wave structure in regions of small Richardson number.
The mountain wave simulations are conducted for a horizontally homogeneous atmosphere characterized by the temperature sounding displayed in Figure 1 and a uniform wind km with a constant vertical grid size of   2 z km and a gravity wave absorbing layer applied above 400 km as described in the previous section. This vertical resolution is appropriate given that waves generated in this atmosphere will have a vertical wavelength (  2 / U N) that ranges between about 14 and 40 km over the depth of the domain.
The mountain terrain is defined as the frequently used bell-shaped mountain: where  100 m h m and a is the mountain half width. The horizontal resolution is set to   / 5 x a and the domain width (  300 x) is sufficiently large that the periodic lateral boundary conditions do not adversely affect the results. The equations are integrated using a time step in seconds of four times the horizontal grid size in kilometers and six acoustic sub-steps per time step. As in the previous test case, we initialize the atmospheric constituents with vertical profiles similar to those shown in Figure 2, which are close to diffusive equilibrium.
Introducing the terrain into the horizontally homogeneous atmosphere at the initial time produces unbalanced disturbances that propagate away from the mountain. In mountain wave simulations, these disturbances typically diminish rapidly during the transient response and do not adversely affect the evolving wave structure. However, for the very deep atmospheric domain in our simulations these acoustic disturbances can amplify rapidly as they propagate upward into the thermosphere, producing spurious results and often numerical instability. These disturbances are particularly troublesome in the absence of molecular viscosity and thermal conductivity. To damp out initially induced perturbations, we include Rayleigh damping terms in the prognostic equations at early times to constrain deviations from their initial state. These terms have the form: for a prognostic variable  having a value  0 in the initially undisturbed state. This procedure appears to work well for these idealized mountain wave simulations, but leaves open the question of how best to deal with these initial disturbances (if necessary) in real-atmosphere simulations.
In conducting these mountain wave test simulations, we consider different horizontal scales of motion by specifying several different mountain half-widths a. Setting  100 a km, the wave response should be hydrostatic (  / 20 aN U ) and rotational effects should be weak ( of the waves should be nearly vertical, and this behavior is confirmed in Figure 5a, which displays the vertical velocity at an integration time of about 3 days. By this time, the wave structure is essentially steady with the wave amplitude concentrated directly above the terrain. To control small-scale noise in this simulation, we have specified the dimensionless constant eddy mixing to be   0.03. (The triangular shape at bottom of Figure 5 and subsequent figures are included only to indicate the location of the mountain, which is actually only 100 m in height.) At larger horizontal scales rotational influences become stronger and the waves propagate downstream as they progress upward. This is apparent in the vertical velocity field in Figure 5b for the simulation with  500 a km at about 16 days, in which the wave train extends noticeably downstream at upper levels. This wave structure is not completely steady at this time as it is slowly progressing farther downstream with time.
For this simulation we are able to reduce the constant vertical mixing to   0.005, which permits the Richardson number dependent portion of the eddy mixing to play a more significant role in influencing the wave structure. The spatial distribution of the mixing coefficient m K , displayed in Figure 6a, reveals that this dissipation is active throughout the wave train, with the maximum dimensionless value reaching To avoid an artificial spinup of the atmospheric constituents from an aerostatic balance, for these simulations we initialized these constituents with vertical profiles close to those depicted by the solid lines in Figure 2.
Since the molecular diffusive transport and vertical eddy mixing terms are nearly in balance, perturbations of the constituents are modified primarily by advective transport (as confirmed by comparison with a simulation in which the diffusive transport and eddy mixing terms are turned off). Figure 6b illustrates  bations in the ideal gas constant  R that result from the variations in these constituents. Notice that the perturbations in  R are confined primarily to the layer between about 150 and 350 km, which is consistent with the region in which most of the variations in the mean profile for R occur (see Figure 1). Although the maximum deviations from the initial profile in Figure 1 are less than one percent (also true for p c , not shown), even these small variations in the thermodynamic properties can have a noticeable effect on the results. For example, holding R and p c fixed at their initial values, the simulated wave amplitude for this  500 a km case increases by about 30%. Since the molecular diffusive transport and vertical eddy mixing terms in the atmospheric constituent Equation 13 produce a small amount of vertical redistribution of constituents during the integration, the perturbations in R shown in Figure 6b are displayed as the deviation from the undisturbed vertical profile far from the mountain.
At smaller horizontal scales, nonhydrostatic effects begin to influence the mountain wave structure. This is evident in our test simulation for  25 a km shown in Figure 7a after 1 day, conducted with   0.027 and with no molecular viscosity or thermal conductivity. Due to the nonhydrostatic influences, the wave energy propagates downstream as it progresses upward, contributing to the significant downstream extent of the waves. Here the conditions that contribute most to the nonhydrostatic response are located at the levels of smaller atmospheric stability in the upper thermosphere (due to the high temperatures). Thus, the downstream refraction of the waves are strongest at levels above 200 km.
As mentioned above, the influences of molecular viscosity and thermal conductivity become large in the upper thermosphere and absorb most of the gravity wave energy propagating up from below. This is illustrated in the physically more realistic simulation for  25 a km displayed in Figure 7b with the molecular viscosity and thermal conductivity specified according to 21 and 22. In this case the wave amplitude is confined almost entirely below about 200 km, removing any need for an artificial damping layer to absorb wave energy approaching the upper lid of the domain. Below about 150 km, the viscosity and conductivity have little effect on the mountain waves, and the vertical velocity fields in Figures 7a and 7b below this level are essentially the same. (Note that the contour interval in Figure 7b is only one fifth the interval shown in Figure 7a.) For the mountain wave test cases discussed above, the mean wind speed was set at a uniform   1 50 ms U , which enabled testing of the model numerics for mountain waves propagating into the upper thermosphere in the absence of the large dissipation produced by molecular viscosity and thermal conductivity. However, the actual horizontal winds in the upper thermosphere can often reach several hundred meters per second, particularly during geomagnetically disturbed periods (Akmaev, 2011).
To test the model numerics in the presence of stronger thermospheric winds, we reran these mountain wave cases specifying a mean wind  / ( ) m s K for the  500 a km mountain wave simulation shown in Figure 5b.
where N is the Brunt-Väisälä frequency and k and  are the horizontal and vertical wave numbers, respectively (Klemp & Lilly, 1980). For a given k, as U increases, the vertical wave number  decreases as the denominator in 27 increases. In addition, the nonhydrostatic influences increase [through the numerator in 27] while the rotational influences decrease [through the denominator in 27]. Since the vertical wave length increases significantly as U approaches 300 1 ms in the upper thermosphere, vertical diffusion is much less effective in dissipating vertically propagating wave energy. Thus, it is not practical to simulate this case using only the eddy mixing 23 to control the wave amplitude. For this reason, we only simulate these cases with the inclusion of molecular viscosity and thermal conductivity. Although the mountain waves are almost completely damped out below the middle thermosphere in the cases with a constant   1 50 ms U , with the stronger horizontal winds shown in Figure 8a the waves can penetrate much farther into the thermosphere due to their greater vertical wave length. With these stronger winds, we find that numerical stability can be maintained by reducing the large time step from    4 t x (km) to    2.5 response to become evanescent above the middle mesosphere as aN U / .  0 8 (Klemp & Lilly, 1980). Thus, most of the wave energy is trapped in the lower mesosphere, which displaces the wavefield downstream from the mountain.
Increasing the mountain half width to  100 a km, the stationary waves penetrate further into the upper mesosphere, as shown in Figure 8c. The nonhydrostatic influences become significant as U increases with height and as waves are turned toward the horizontal, they are absorbed by the molecular viscosity and thermal conductivity. This behavior is in contrast to that for this case with   1 50 ms U where wave structure was strongly hydrostatic (Figure 5a).
With a mountain half width further increased to  500 a km, the waves regain a strong hydrostatic character, propagating to the top of the model domain where they are attenuated by the absorbing layer, which is included above  400 z km. For this case, the constant portion of the vertical eddy mixing coefficient has been reduced to   0.025. Rotational influences are greatly reduced from those apparent in the   1 50 ms U case ( Figure 5b); the propagation of energy is nearly vertical and the waves remain directly above the mountain terrain. Note that in Figure 8, the scaling for the horizontal x axes are expanded by a factor of 2 from those in Figures 5-7 for the   1 50 ms U cases to better display the wave structure.
When the waves have the capability to propagate to high altitudes, the kinematic viscosity and thermal diffusivity can influence the structure of the waves in addition to dissipating them. In a detailed analysis of the linear dispersion characteristics in the presence of kinematic viscosity and thermal diffusivity, Vadas and Fritts (2005) confirmed that the waves with the largest vertical wave lengths dissipate at the highest altitudes. They found that as these diffusive influences increase with height, the intrinsic frequency decreases and causes the vertical wave length to decrease, and with sufficiently large diffusion the wave may be reflected. For this case with  500 a km, the large vertical wave length resulting from the strong horizontal winds allows the waves to propagate to the top of the simulated domain. However, it is difficult KLEMP AND SKAMAROCK 10.1029/2021MS002499 15 of 19 to distinguish the relative influences of the increasing winds with height and the increasing dissipation on the wave structure.

Summary
Adapting MPAS for applications that extend well into the thermosphere requires some significant modifications to the model equations and the model numerics. As mentioned in the Introduction, some of these extensions have been implemented in MPAS by Skamarock et al. (2020) in removing the shallow atmosphere approximations from the model equations applied to the sphere. With these modifications the MPAS equations use the actual geocentric distance from the center of the Earth to allow the mesh configuration and gravity to vary with height and include the full Coriolis terms in the momentum equations. To adapt MPAS for model domains that span a large number of vertical scale heights, we have further modified the model numerics and generalized the model equations to include variable composition, and potentially large molecular viscosity and thermal conductivity.
Allowing the atmospheric constituents to vary can have potentially significant influences on the model dynamics. These feedbacks occur through changes in the ideal gas constant d R and the heat capacity p c , which depend solely on the relative concentrations of these constituents. With variable thermodynamic coefficients, potential temperature is no longer strictly conserved for adiabatic flow. Nevertheless, in our numerical testing with MPAS it appears that potential temperature can still be effectively used as a prognostic variable with an additional term that accounts for the variability in   / d p R c .
The height-based hybrid terrain-following coordinate in MPAS seems well suited for deep atmosphere domains, and avoids the complications that arise in a pressure-based vertical coordinate when gravity becomes a function of height. Also, the split-explicit numerical integration used in MPAS appears to remain viable in the limited deep atmosphere testing we have conducted across a range of horizontal scales and that includes the influences of variable atmospheric composition. As expected, the presence of molecular viscosity and thermal conductivity in the mid and upper thermosphere strongly damps disturbances propagating up from the lower atmosphere. However, it is interesting to see in our mountain wave simulations that waves having large horizontal wave lengths can penetrate throughout the thermosphere in the presence of strong thermospheric winds that promote large vertical wave lengths that are not as strongly damped by the molecular viscosity and thermal conductivity.
There are a number of issues that may require further consideration. In writing the governing Equations 8-13 we have incorporated simplifying approximations as discussed above that will require further scrutiny in our continuing efforts to improve the accuracy of the modeling system. We have also mentioned that the numerical integration appears to be sensitive to the specific numerics used to represent the variations in the thermodynamic coefficients in the potential temperature equation. While the numerics described by 15 appear to maintain numerical stability for this feedback, further testing may be needed to determine whether this treatment is sufficiently robust. Another issue is the rapid growth of acoustic disturbances that arise in our deep-atmosphere mountain wave test cases when the terrain is introduced at the initial time. While including Rayleigh damping at early times in the simulation appears to remove these disturbances, it remains unclear if these initial imbalances will be an even larger problem in real data simulations with much larger amplitude terrain. Also, as discussed above, the large values of molecular diffusivity and thermal conductivity in the upper thermosphere make it difficult to accurately represent these damping influences. Although implicit representation of these terms maintains numerical stability, the magnitude of damping in the current implementation may significantly under represent their true influences. Further consideration of these issues will be explored in more realistic simulations of thermospheric dynamics, and with the inclusion of appropriate thermospheric physics.
KLEMP AND SKAMAROCK 10.1029/2021MS002499 18 of 19 Figure A1. Time series of the maximum vertical velocity for the (Klemp, 2011) resting-atmosphere test case using the original horizontal pressure gradient formulation A1 from (Klemp, 2011) (black line), and using the revised formulation A2 currently implemented in MPAS (red line), using either full thermodynamic variables x p or perturbation variables p x  .