## Abstract

The formulation of a fully compressible nonhydrostatic atmospheric model called the Model for Prediction Across Scales–Atmosphere (MPAS-A) is described. The solver is discretized using centroidal Voronoi meshes and a C-grid staggering of the prognostic variables, and it incorporates a split-explicit time-integration technique used in many existing nonhydrostatic meso- and cloud-scale models. MPAS can be applied to the globe, over limited areas of the globe, and on Cartesian planes. The Voronoi meshes are unstructured grids that permit variable horizontal resolution. These meshes allow for applications beyond uniform-resolution NWP and climate prediction, in particular allowing embedded high-resolution regions to be used for regional NWP and regional climate applications. The rationales for aspects of this formulation are discussed, and results from tests for nonhydrostatic flows on Cartesian planes and for large-scale flow on the sphere are presented. The results indicate that the solver is as accurate as existing nonhydrostatic solvers for nonhydrostatic-scale flows, and has accuracy comparable to existing global models using icosahedral (hexagonal) meshes for large-scale flows in idealized tests. Preliminary full-physics forecast results indicate that the solver formulation is robust and that the variable-resolution-mesh solutions are well resolved and exhibit no obvious problems in the mesh-transition zones.

## 1. Introduction

Computational capabilities that are expected to be generally available in the near future will enable global atmospheric model applications permitting explicitly resolved nonhydrostatic motions that will require the solution of the nonhydrostatic equations. The applications will likely include both variable-resolution mesh modeling with limited areas allowing nonhydrostatic motions, such as the exploratory efforts in the global nonhydrostatic variable-resolution modeling of Yeh et al. (2002), and global high-resolution uniform-mesh modeling where nonhydrostatic scales are permitted everywhere on the globe, as in Satoh et al. (2008). Development efforts for atmospheric models capable of these applications are in their early stages and face significant challenges in discretizing the nonhydrostatic fluid-flow equations on the sphere (Williamson 2007), in allowing selective refinement, and in achieving efficient scaling on massively parallel computer architectures.

We have constructed a global fully compressible nonhydrostatic model using finite-volume numerics discretized on centroidal Voronoi (nominally hexagonal) meshes using C-grid staggering of the prognostic variables based on the work of Thuburn et al. (2009) and Ringler et al. (2010). This model is called the Model for Prediction Across Scales (MPAS), and its development is a collaborative project being led by the National Center for Atmospheric Research (NCAR) and Los Alamos National Laboratory (LANL). NCAR is responsible for developing the MPAS atmospheric component, LANL is responsible for the ocean component, and the shared software infrastructure is being developed jointly.

Quasi-uniform centroidal Voronoi meshes are similar to icosahedral (hexagonal) meshes, such as those used in the nonhydrostatic icosahedral atmospheric model (NICAM; Satoh et al. 2008), and they provide nearly uniform resolution over the globe. In contrast, global atmospheric models have typically employed latitude–longitude grids for their discretization (Williamson 2007). Latitude–longitude mesh implementations of finite-difference, finite-volume, and spherical transform methods do not scale well on the latest generations of supercomputers that rely on large numbers of distributed-memory processing elements; the domain decomposition needed for efficient distributed-memory parallelism conflicts with the need to use global transforms of some sort. Finite-difference and finite-volume methods used on latitude–longitude meshes need polar filters and other extensions to the numerical schemes because of the convergence of grid lines at the poles. Spherical transform methods, while not requiring polar filtering, also do not scale optimally with increasing resolution, and the semi-Lagrangian transport with which they are often paired needs special attention in the polar regions in addition to polar filtering. As with the icosahedral meshes, Voronoi meshes do not require polar filtering, and parallelization based on standard horizontal (2D) domain decomposition is possible and appears to allow good weak scaling performance on massively parallel architectures in our early tests.

The centroidal Voronoi meshes also allow for local refinement, and the variable-resolution horizontal mesh takes advantage of the unstructured-mesh capabilities of our nonhydrostatic solver. The variable-resolution meshes are generated such that there is a gradual change in mesh density from the coarse- to the high-resolution regions (Ringler et al. 2008), and they allow much more flexible mesh-refinement capabilities than approaches using a remapping of a structured mesh [e.g., the clustering of latitude and longitude grid lines used in Yeh et al. (2002) or the remapping of the global icosahedral mesh employed in a shallow-water (SW) equation solver by Tomita (2008)]. The smooth mesh transitions we use stand in contrast to the abrupt mesh transitions used in traditional two-way nested models (e.g., Skamarock et al. 2008) or in mesh refinement achieved directly through cell division (e.g., Walko and Avissar 2008). We believe the smooth mesh transition will ameliorate many of the difficulties associated with traditional nesting approaches, as indicated in the results of Ringler et al. (2011) using the shallow-water equations. Thorough testing of the variable-mesh capabilities of the MPAS nonhydrostatic solver is not a focus of this paper; we expect to report test results for this important capability in future publications.

The approach we are using to solve the nonhydrostatic equations can be considered an extension of existing techniques used in nonhydrostatic models, such as the Advanced Research Weather Research and Forecasting model (ARW-WRF; Skamarock and Klemp 2008), to the horizontal Voronoi mesh, and this includes the use of C-grid staggering. Using established techniques allows us to directly compare our Voronoi mesh results with those from state-of-the-art cloud–mesoscale models and to establish the accuracy and efficiency of the solver at the smallest resolved scales. Our formulation differs from other nonhydrostatic icosahedral mesh models, such as NICAM, in the use of C-grid staggering and the centroidal Voronoi mesh. C-grid staggering had not previously been incorporated in Voronoi or icosahedral hexagonal mesh models because of numerical problems recently resolved by Thuburn et al. (2009) and Ringler et al. (2010).

In this paper, we focus on describing the nonhydrostatic fluid-flow solver and demonstrate its accuracy in applications on the sphere and in nonhydrostatic test cases on Cartesian domains. We begin in section 2 by presenting the equations and the finite-volume split-explicit discretization used in the unstructured-mesh solver, and we discuss some aspects of the variable unstructured mesh and considerations involving the discretization. In section 3 we present nonhydrostatic test results on a Cartesian plane and compare them with results from existing nonhydrostatic solvers. Tests indicate that the solver is producing solutions comparable to those of existing cloud models. Test simulation results for large-scale flow on the sphere are presented in section 4, and include both idealized tests and example full-physics forecasts on uniform and variable-resolution meshes. These test results indicate that the solver is performing as well or better than existing Voronoi and icosahedral mesh models based on the results presented here and on shallow-water tests and theoretical considerations presented in the previous work of Ringler et al. (2011). We conclude in section 5 with a summary.

## 2. Solver formulation

### a. Continuous equations

We cast the governing equations in terms of a new height-based terrain-following vertical coordinate *ζ* following Klemp (2011). In this formulation, illustrated for simplicity on a Cartesian grid, the height of the coordinate surface is defined by

where *ζ* represents the nominal heights (ignoring terrain) of the coordinate surfaces, *A*(*ζ*) defines the relative weighting between the terrain-following coordinate and the pure height coordinate with 0 ≤ *A* ≤ 1 − *ζ*/*z _{t}*,

*z*is the height of the model top, and the array

_{t}*h*is specified to produce increased smoothing of the terrain influence with height with the requirement that

_{s}*h*(

_{s}*x*,

*y*, 0) =

*h*(

*x*,

*y*). For

*A*(

*ζ*) = 1 −

*ζ*/

*z*and

_{t}*h*(

_{s}*x*,

*y*,

*ζ*) =

*h*(

*x*,

*y*) (the terrain height), we recover the traditional terrain-following height coordinate, and for

*A*(

*ζ*) = 0 we recover a pure height coordinate. As described in Klemp (2011),

*A*(

*ζ*) and

*h*can be chosen to control the amount and scale of terrain influence on the vertical coordinate. The metric terms associated with the vertical coordinate transformation are represented by

_{s}**∇**

*ζ*= (

*ζ*,

_{x}*ζ*,

_{y}*ζ*), with

_{z}

*ζ*_{H}= (

*ζ*,

_{x}*ζ*) and

_{y}**z**

_{H}= −

*ζ*_{H}/

*ζ*. We use the traditional terrain-following height-coordinate option in all the test results presented in this paper.

_{z}We cast the continuous prognostic equations using the flux variables

where *ρ _{d}* is the density of dry air, ,

*q*represents the mixing ratio of the respective water species, and

_{j}is a modified moist potential temperature with *q _{υ}* representing the water vapor mixing ratio. The velocities (

*u*,

*υ*,

*w*) represent two orthogonal horizontal velocities and a vertical velocity radially outward from the center of the earth; the velocities obey the right-hand rule, , where

**i**×

**j**=

**k**. In the transformed coordinate Ω =

**V**·

**∇**

*ζ*is the component of the mass flux normal to the coordinate surfaces, where

**V**= (

*U*,

*V*,

*W*). The full equation set can be expressed as

Pressure is obtained diagnostically from the equation of state,

with *γ* = *c _{p}*/

*c*. Here,

_{υ}*ρ*is the density of moist air, and

_{m}where *q _{υ}*,

*q*, and

_{c}*q*are mixing ratios of vapor, cloud water, rainwater, etc. In (3),

_{r}*η*=

**k**·

**∇**×

*υ*_{H}+

*f*is the absolute vertical vorticity, and

*K*= |

*υ*_{H}|

^{2}/2 is the horizontal kinetic energy. In the curvature and Coriolis terms in (3) and (4),

*f*= 2Ω

_{e}sin

*ψ*,

*e*= 2Ω

_{e}cos

*ψ*,

*ψ*is the latitude, Ω

_{e}is the angular rotation rate of the earth,

*r*is the earth’s radius, and

_{e}*α*is the rotation angle between the line normal to the horizontal velocity and the meridians. The terms ,

_{r}*F*, , and represent sources and sinks from physics, subgrid models, and filters. Following the notation of Dutton (1986), we define the equation (

_{W}**∇**·

**V**

*b*)

_{ζ}=

**∇**

_{ζ}· (

**V**

_{H}

*b*) + ∂(Ω

*b*)/∂

*ζ*for any scalar

*b*, where

**∇**

_{ζ}refers to the divergence operator along a constant

*ζ*surface and

**V**

_{H}= (

*U*,

*V*).

The compressible nonhydrostatic equations (3)–(7) are written in flux form, with the horizontal momentum equations expressed in vector-invariant form, to help achieve the desired conservation properties in the discrete model. The use of the vector-invariant form of the horizontal momentum equation (3) also allows us to bypass the potentially problematic discretization of nonlinear horizontal momentum transport. Our discrete formulation, described in the following subsections, exactly conserves dry air mass, scalar mass (e.g., moisture), and Θ_{m}. The horizontal discretization is largely taken from the shallow-water model discretization of Ringler et al. (2010), with exceptions noted in section 2c below. The 3D discrete nonhydrostatic model, however, does not inherit the exact conservation of potential vorticity and energy (to time truncation) of the SW model.

### b. Temporal discretization

The fully compressible nonhydrostatic equations are solved using the time-split integration technique described in Klemp et al. (2007) for the height-coordinate equations. As in Klemp et al., we recast the governing equations (3)–(7) in terms of thermodynamic variables that are perturbations from a hydrostatically balanced reference state that is only a function of the geometric height *z*. In this way we reduce the truncation error associated with the horizontal pressure gradient terms and the roundoff error associated with the right-hand-side terms in the vertical momentum equation. We use the third-order Runge–Kutta scheme and explicit time splitting as described in Wicker and Skamarock (2002). The time-splitting technique integrates gravity waves and horizontally propagating acoustic waves on smaller explicit substeps within the three-step Runge–Kutta time integration. Vertically propagating acoustic waves are integrated implicitly. Thus, the Runge–Kutta time step is limited by the maximum advective velocity, and the acoustic time step is limited by the horizontal acoustic wave speed. The use of the vector-invariant form of the horizontal momentum equation (3) does not complicate this solution procedure, but note that we cast this equation in flux form to facilitate the acoustic integration. This gives rise to the additional term **υ**_{H }**∇**_{ζ} · **V** that does not appear in nonflux-form equations that are usually used in this context (e.g., Ringler et al. 2010). Again, following Klemp et al. (2007), we cast the acoustic-step equations in terms of perturbation variables from the values at time step *t* for the integration from *t* to *t* + Δ*t*.

### c. Spatial discretization

In the MPAS atmospheric solver, the continuous equations are spatially discretized on an unstructured C-grid centroidal Voronoi mesh following Thuburn et al. (2009) and Ringler et al. (2010). Quasi-uniform meshes as well as variable-resolution meshes can be generated for the sphere, for regional domains on the sphere, and for Cartesian planes using techniques described in Ringler et al. (2008). It is our use of the C-grid centroidal Voronoi mesh that distinguishes MPAS from other models such as NICAM.

A horizontal C-grid Voronoi mesh is depicted in Fig. 1. The horizontal momentum normal to the cell edge (*u* in Fig. 1) is prognosed at the cell edges. The coupled potential temperature Θ_{m}, density , and moisture *Q _{j}* are prognosed at the cell centers where they represent cell-averaged values in the finite-volume formulation. The vertical momentum is prognosed on the vertical cell faces located half a grid level above and below the cell center, consistent with a 3D C-grid discretization. All other quantities are diagnosed from the prognostic variables [e.g., pressure (8)]. The tangential component of the horizontal momentum needed at the cell edges,

**k**×

**V**

_{H}in (3), is diagnosed following Thuburn et al. (2009). The horizontal momentum components (

*U*,

*V*) needed at the cell center in the vertical momentum equation (4) are computed from a radial basis function reconstruction (Bonaventura et al. 2011) using the prognosed horizontal velocities at the cell edges. The absolute vertical vorticity

*η*is diagnosed following Ringler et al. (2010). Using the tangential velocity reconstruction described by Thuburn et al. (2009), this horizontal discretization for the C grid does not suffer from the problems of the nonstationary geostrophic mode (Ničković et al. 2002). The horizontal discretization conserves mass to machine roundoff.

We have two options for kinetic energy *K* defined at the cell centers and used in (3). The approach described in Ringler et al. (2010) defines the kinetic energy in the cell, here named *K _{c}*, as a sum over the edges of the cell:

where *u _{e}* is the edge normal velocity,

*A*is the area of the cell,

_{c}*l*is the edge length, and

_{e}*d*is the distance between cell centers sharing edge

_{e}*e*(see Fig. 1). The second option, developed by A. Gassmann (2011, personal communication), begins with a definition of the kinetic energy at cell vertices, which we name

*K*:

_{υ}where the three edges are those edges that share vertex *υ*, and *A _{υ}* is the area of the triangle containing the vertex (the triangle depicted by the dashed lines in Fig. 1). Using

*K*we construct a second cell-centered value of KE, named , by computing an area-weighted sum of kinetic energies

_{υ}*K*at the vertices of a given cell:

_{υ}where is the shaded area associated with a vertex of a given cell, as depicted in Fig. 1. The new value of the cell-centered kinetic energy used in (3) is given by the weighted sum of *K _{c}* [(10)] and [(12)]:

where *α* is a weighting coefficient; *α* = 1 recovers the Ringler et al. (2010) formulation (10) and *α* = ⅜ is used in the Gassmann formulation. The Ringler et al. (2010) formulation (10) guarantees energy conservation to the time-truncation error within their fully nonlinear SW equations solver. The formulation from Gassmann (*α* ≠ 0) will not conserve energy in the SW solver, but we have found that it removes a computational instability we have encountered in large-scale simulations of baroclinic waves. Further discussion of these kinetic energy formulations and simulation results is given in section 4.

The transport scheme used for the horizontal flux divergence calculations in (4), (5), and (7) on the irregular Voronoi mesh is described in Skamarock and Gassmann (2011). It uses nominally third- and fourth-order approximations to the scalar gradient to construct the scalar values on the cell edges in the conservative flux divergence calculation. The mass flux divergence in (6) is approximated by averaging the cell-averaged values of to the cell face from the two cells sharing the face (see Thuburn et al. 2009). We use the mass flux that is time averaged over the acoustic steps for scalar transport, thus maintaining consistency between the scalar transport and the continuity equation. The vertical flux divergence terms in (3), (4), (5), and (7) are calculated using the third- and fourth-order schemes as described in Wicker and Skamarock (2002) but using the reduced dissipation in the third-order scheme as described in Skamarock and Gassmann (2011). The monotonic option available in the Skamarock and Gassmann (2011) scheme is used for moisture in the NWP forecasts presented in section 4c. The monotonicity constraint introduces dissipation, and we do not use any additional explicit dissipation for moisture when the monotonic constraint is used.

As explained in Klemp et al. (2003), it is important to maintain consistency between the transport operator and the metric terms associated with the terrain transformation. In MPAS this requires that the diagnosis of the transformed vertical velocity Ω is consistent with the advection operator used in the thermodynamic equation (5). To satisfy this requirement, we formulate the transformed vertical velocity for cell *i* as

where *z _{i}* is the coordinate-surface height at cell

*i*for which Ω

_{i}is being computed. We apply the same Skamarock and Gassmann (2011) scheme to the right-hand side of (14), thus satisfying the consistency requirement. The expensive part of the Skamarock and Gassmann scheme, evaluating (

*z*−

*z*) at the cell faces, can be computed and stored before integration begins; thus, there is little additional cost in implementing the Klemp et al. (2003) consistency requirement.

_{i}Horizontal filters using the Laplacian are evaluated using a standard finite-volume discretization on the Voronoi mesh. The discrete version of the Laplacian ∇^{2}*ψ* = **∇** · **∇***ψ* of the scalar *ψ* at cell *i* is expressed as

where *n _{e}* is the number of edges making up cell

*i*,

*l*is the length of edge

_{e}*e*on cell

*i*,

*ψ*is the scalar value in cell

_{e}*j*sharing edge

*e*with cell

*i*, and

*d*is the distance between the centers of cells

_{e}*i*and

*j*. The formulation (15) follows directly from the definitions of the divergence and gradient operators in Thuburn et al. (2009) and Ringler et al. (2010), and experimental results show that (15) exhibits second-order convergence on smooth Voronoi meshes. A fourth-order hyperdiffusion is enabled by taking the discrete Laplacian of the discrete Laplacian (15).

The Laplacian of the momentum *u* is evaluated using the vector identity

where *u _{i}* is the edge normal velocity defined on cell edge

*i*,

*η*is the relative vertical vorticity, and

**∇**

_{ζ}·

**υ**is the horizontal divergence. The first term on the right-hand side of (16) is the difference between the horizontal divergences of the two cells sharing edge

*i*, and the second term is the difference of the vorticities defined at the vertices on either side of edge

*i*. The evaluation of the vertical vorticity at the vertices and the divergence for each cell is described in Ringler et al. (2010). As with scalar diffusion, a fourth-order hyperdiffusion for momentum can be evaluated by applying this discrete operator twice.

The horizontal filtering formulation of Smagorinsky (1963) uses the second-order Laplacians (15) and (16) along with an eddy viscosity *K _{h}* defined on the Cartesian plane using the Cartesian velocities (

*u*,

*υ*):

We have implemented this approach by evaluating *K _{h}* in (17) at cell centers on the sphere. The eddy viscosities are averaged to the cell faces such that the diffusion operator is in conservation form:

The evaluation of the eddy viscosity *K _{h}* is accomplished by projecting the velocities onto a tangent plane, integrating the squared terms involving the velocities in (17) over the cell, and applying Stokes theorem to transform the cell integrals to discrete line integrals around the cell edge. The evaluation is inexpensive because the prognostic normal and diagnostic tangential velocities exist, and the time-independent coefficients for the line integrals are precomputed and stored before the time integration begins.

### d. Discretization considerations

Compared with other Voronoi and icosahedral mesh models using finite-volume or finite-difference formulations, there are a number of new aspects in the MPAS solver formulation that benefit global multiscale (hydrostatic and nonhydrostatic scale) atmospheric simulations.

The first advance in MPAS is the formulation for the C-grid staggering of the prognostic variables. Ničković et al. (2002) examined hexagonal C-grid formulations for the shallow-water equations on an *f* plane. The tangential velocity needs to be reconstructed from the prognosed normal velocities on the cell faces of the C grid to evaluate the Coriolis term, and Ničković et al. showed that the stationary geostrophic mode in the linearized SW equations will be nonstationary using the most obvious reconstruction of the tangential velocity. The nonstationary geostrophic mode produces distortion in rotational modes that renders the scheme useless for most applications. Thuburn (2008) and Thuburn et al. (2009) developed a method to reconstruct the tangential velocity such that the geostrophic mode remains stationary for the linear SW equations. Ringler et al. (2010) extended this formulation to the nonlinear SW equations such that it conserves potential vorticity and potential enstrophy, conserves energy to the time truncation error, and allows for the dissipation of potential enstrophy following Sadourny and Basdevant (1985). These advances are based on the constraint that the horizontal mass divergence on the triangular mesh (the dashed triangular cell in Fig. 1; the dual of the hexagonal mesh) that is computed using the reconstructed tangential velocities is equal to the area-weighted sum of the divergences in the hexagons underlying the triangular cell. These conservation properties can be considered a generalization of the conservative rectangular-mesh discretizations of Sadourny (1975) and Arakawa and Lamb (1981).

Our use of C-grid staggering is guided by theory and practical experience. Mesoscale and cloud-scale motions are dominated by horizontally divergent gravity wave motions, and C-grid staggering provides twice the resolution of divergent modes compared to the unstaggered (A) grid; it does not require any averaging of the velocities or pressures in the pressure gradient and divergence terms as is required in the A-, B-, D-, and E-grid staggerings. Pressure and velocity averaging lead to stationary grid-scale modes (often referred to as *parasitic* modes) that must be filtered, and the parasitic modes associated with the divergence and pressure gradient terms will require a higher level of filtering on these meshes. Our experience is that the level of filtering needed on these other meshes is considerably higher than that needed to provide sinks for the downscale energy and enstrophy cascades in our full nonlinear simulations using the C-grid staggering. In practice, we find that solvers not using C-grid staggering need finer meshes to produce similarly resolved features (such as clouds). Our experience also indicates that higher-order differencing of the mass divergence and pressure gradient terms does not appreciably change the behavior for the non-C-grid meshes; parasitic modes remain, increased filtering is still needed, and higher mesh densities are still needed to produce comparably resolved solutions. Generally speaking, all grid staggerings can be made to work with some level of filtering, and the choices affect scheme efficiency (accuracy versus cost). We have found that the C-grid discretization results in the highest efficiency.

Randall’s (1994) analysis of geostrophic adjustment indicates that the C-grid staggering is not optimal for large-scale flows. Our intended applications for MPAS are cell spacings of the order 100 km and less. Our tests indicate that MPAS produces solutions similar in accuracy to other finite-volume and finite-difference models for large-scale flows [based on results in Ringler et al. (2010) for shallow-water tests, and our results for the Jablonowski and Williamson (2006a) baroclinic wave tests]. We have not identified any problems arising from the C-grid discretization or with the computational Rossby modes addressed in Thuburn et al. (2009). For higher-resolution applications (to resolve mesoscale circulations), the rotational modes will be well resolved regardless of the specific numerics. In this regime, numerical accuracy will be most strongly influenced by the treatment of the gravity wave modes.

The second advance incorporated into MPAS is the use of higher-order operators within the advection scheme. Atmospheric models using unstructured, irregular meshes often employ simple low-order operators within the advection algorithms. For example, the flux divergence term in the scalar conservation equation requires that mass flux be evaluated at the cell edge, and a common second-order formulation uses the average of the scalar values in the two cells sharing that edge. However, these schemes do not produce solutions on hexagonal meshes that are as accurate as those produced in state-of-the-art rectangular-mesh cloud and mesoscale models. Within MPAS we have implemented a transport scheme, described in Skamarock and Gassmann (2011), that uses a least squares fit polynomial to evaluate higher-order scalar derivatives used in the third- and fourth-order transport schemes [similar to that employed in ARW-WRF (Skamarock and Klemp 2008)]. A significant increase in accuracy is demonstrated with the Skamarock and Gassmann (2011) scheme over that obtained using other second-order formulations.

Finally, unstructured-mesh solvers such as MPAS make use of indirect addressing when building the horizontal operators during a time step. We keep vertical columns contiguous in memory in our FORTRAN implementations, and we find that our solvers have computational efficiencies similar to our rectangular (structured) grid solvers (e.g., ARW-WRF). MacDonald et al. (2011) have examined the question of efficiency in 3D atmospheric solvers and find that, on cache-based computer architectures, unstructured-mesh solvers having vertical columns contiguous in memory are as efficient as their structured-mesh counterparts.

## 3. Test cases: Nonhydrostatic flows on Cartesian planes

It is very costly and difficult to assess the nonhydrostatic response of global atmospheric solvers because very high resolutions are needed to resolve nonhydrostatic scales, particularly for convection, where cell spacings of a few kilometers or less are needed. The MPAS unstructured Voronoi mesh and solver can be used on Cartesian planes in addition to the sphere, and we have performed extensive tests of the solver for nonhydrostatic-scale motions on Cartesian planes using both 2D (*x*, *y*) and 3D test cases, including mountain waves, 2D and 3D squall lines, and 3D supercell thunderstorms. For the mountain wave cases we have exact linear and nonlinear solutions, and these tests have helped us verify the correctness of our coding. We present results from the Schär et al. (2002) test case that has been used to examine the numerical treatment of terrain and some aspects of the consistency of the discretization. A strongly nonlinear 2D density current test from Straka et al. (1993) is presented to demonstrate the nonlinear response of the numerics. In our experience, we have found that deep moist convection provides the most challenging tests of nonhydrostatic solver robustness (stability) and accuracy because of the significant latent heating occurring near the grid scale. Thus, we also present results from 3D idealized supercell thunderstorm simulations to demonstrate the robustness of the MPAS nonhydrostatic solver.

### a. 2D Schär test case

Schär et al. (2002) proposed a test case using flow over terrain containing small-scale structure that has been used to uncover some problems within terrain-following coordinate models. In the Schär et al. study, smoothing the coordinate surfaces helped remove the spurious motions generated using particular model formulations. Klemp et al. (2003) showed that the spurious motions examined by Schär et al. were associated with an inconsistency between the transport terms and the diagnosis of the vertical velocity in terrain-following coordinate model formulations.

The test case uses a terrain height,

where *H* = 250 m, *λ* = 4000 m, and *a* = 5000 m. A constant mean-state buoyancy frequency *N* = 0.01 s^{−1} is prescribed along with a constant horizontal environmental wind *U* = 10 m s^{−1}. Schär et al. (2002) and Klemp et al. (2003) use horizontal mesh spacings of 500 m and a vertical mesh spacing of 300 m nominally in their tests, and we use these same mesh spacings. The horizontal mesh in the 2D (*x*, *y*) configuration for MPAS is constructed using two rows of perfect hexagons in *x* along with periodic boundary conditions in *y*.

Figure 2 depicts the vertical velocity in the steady-state solutions for three test configurations; the first (Fig. 2a) uses a fully second-order model configuration, the second (Fig. 2b) uses the fourth-order transport scheme for the potential temperature (see Skamarock and Gassmann 2011) but only a second-order evaluation of the transformed vertical velocity using (14), and the third (Fig. 2c) uses the fourth-order scheme for both transport and the transformed vertical velocity diagnosis. The results can be compared with those in Klemp et al. (2003). The solution depicted in Fig. 2c, using the fourth-order scheme for both transport and the transformed vertical velocity diagnosis, best reproduces the analytic solution, and the second-order solution depicted in Fig. 2a reproduces the analytic solution fairly well. In contrast, the discretization using the fourth-order transport and second-order vertical velocity diagnosis in Fig. 2b has a very large error in the vertical velocity field above the mountain. This error is the same as that found in Schär et al. (2002) and Klemp et al. (2003). Our model formulation follows Klemp et al. by employing a consistent diagnosis of the vertical velocity using (14) to remove the source of the error.

### b. 2D density current

The density current test case is described in Straka et al. (1993). This 2D simulation of a highly nonlinear density current uses a fixed physical viscosity of 75 m^{2} s^{−1}; hence, converged solutions can be computed. The domain extends from −25.6 to +25.6 km in *x* and from 0 to 6.4 km in *z*, and the boundary conditions are periodic in *x* and upper solid and lower free-slip surfaces in *z*. There is no flow at the initial time and the atmosphere is neutrally stratified. The density current is produced by specifying an initial thermal perturbation of the form

where *R*^{2} = [(*x* − *x _{c}*)/

*x*]

_{r}^{2}+ [(

*z*−

*z*)/

_{c}*z*]

_{r}^{2},

*x*= 0 km,

_{c}*x*= 4 km,

_{r}*z*= 3 km, and

_{c}*z*= 2 km. This cold bubble descends to the surface and spreads out to produce left- and right-moving density currents.

_{r}For these tests, MPAS is configured to use the third-order upwind scheme from Skamarock and Gassmann (2011) for potential temperature transport with an upwinding coefficient of 0.25. Figure 3 depicts the simulation results at 900 s and these results can be directly compared with those from Straka et al. (1993). While the 25-m mesh solution is essentially converged, the 100-m mesh solution is beginning to exhibit noticeable differences from the converged solution (none of the published model results are visibly converged on the 100-m mesh). For example, the maximum temperature in the second eddy behind the leading edge is colder, and there is some evidence of overshoots in the temperatures closest to the surface. These errors are significantly larger on the 200-m mesh where many features can no longer be resolved. Using a standard second-order advection formulation (bottom panel of Fig. 3), the solution error increases dramatically because of the less-accurate transport scheme. Overall, the MPAS solutions compare well with other published cloud-model solutions for this test, including the model solutions presented in Straka et al. (1993) and more recently those for the Advanced Regional Prediction System (ARPS) model (Xue et al. 2000) and the ARW-WRF model (Skamarock and Klemp 2008).

### c. Supercell simulation

We have performed 3D supercell simulations with MPAS and a rectangular-mesh model. The horizontal domain (*x*, *y*) is square and doubly periodic with a length of 84 km. The domain height is 20 km and 40 levels are used in the vertical decomposition (Δ*z* = 500 m). Both models use a 4-s time step and 28 224 thermodynamic points on a horizontal plane, with cell-center spacings of 500 and ~539 m for the rectangular and Voronoi meshes, respectively. The horizontally homogenous environment is initialized using the sounding from Weisman and Klemp (1982) and a unidirectional-shear hodograph with a horizontal wind speed of zero at the surface, increasing linearly to 25 m s^{−1} at 5 km and constant above 5 km. A positive thermal perturbation in *θ* is used to initiate convection as described in Weisman and Klemp. Mirror-image right- and left-moving supercells are produced in these simulations, and depictions of the storms from the rectangular-mesh model and MPAS are given in Fig. 4, which shows the vertical velocity and rainwater fields at *z* = 5 km and the surface cold pool. The storms are very similar on both meshes, and perfect symmetry is maintained in both solvers. The solutions compare well with those published in the literature (e.g., Weisman and Klemp 1982). Moreover, the solutions are much more similar to each other than those generated using other discrete model formulations [e.g., Weisman and Klemp (1982) use the advective form of the governing equations], indicating that the differences in mesh configuration lead to solution differences that are very small relative to other possible model formulation differences.

The maximum vertical velocity as a function of time for these two simulations is given in Fig. 5. The maximum velocities are very similar in the two simulations, especially for the first half-hour of the simulation. They are sensitive to the configuration of the models, and we are using the same discretizations and filtering (second-order dissipation with a constant eddy viscosity *ν* = 500 m s^{−1}). The similarity of these two solutions to each other and to those in the cloud-modeling literature leads us to conclude that MPAS produces convective solutions of similar accuracy, quality, and robustness to our state-of-the-art cloud models, such as ARW-WRF (Skamarock and Klemp 2008). We have also found that MPAS has similar computational efficiency compared to rectangular-mesh cloud models such as ARW-WRF; the run times are within a few percent of each other and exhibit variability of a similar magnitude for different compilers, compiler options, and even runs on the same machines. Importantly, the cost of the additional horizontal momentum equation in MPAS [there are nominally three horizontal velocities on the Voronoi mesh for each thermodynamic point (cell) compared to two horizontal velocities for rectangular meshes] is offset by the ability to take longer time steps; experience, and a linear stability analysis, indicate that we can take a time step 1.25 times greater on the Voronoi mesh compared to the time step allowed on a rectangular mesh.

## 4. Test cases: Large-scale flow on the sphere

A standard test for 3D solvers on the sphere is the Jablonowski and Williamson (2006a) baroclinic wave test. We have simulated the baroclinic wave as described in the reference using horizontal mesh resolutions from approximately 480-km cell-center spacing (2562-cell mesh) to approximately 30-km cell spacing (655 362-cell mesh). The initial state consists of identical zonally symmetric unstable jets in both the northern and southern Hemispheres and represents an unstable but steady-state solution. Tests include simulations without any initial perturbation, in which case the steady solution should be maintained, and simulations of a baroclinic wave train that is triggered with a small perturbation in the Northern Hemisphere zonal wind.

To initialize the model state for these tests, we begin by computing the initial state on a 2D (*y*, *z*) mesh on the sphere where the horizontal dimension *y* extends from the south pole to the north pole and the vertical mesh spacing is that used for the full 3D mesh. We iteratively solve for the hydrostatically balanced thermodynamic state in each column in a manner similar to that described in Jablonowski and Williamson (2006a). We have found that there can be small but significant geostrophic imbalances using the analytic zonal velocities; thus we set the zonal velocities on the 2D mesh to the discrete geostrophic velocity *u _{g}*. We use a moist version of Jablonowski and Williamson’s (2006b) Eq. (A.2) to compute this balanced velocity:

We initialize the horizontal velocities by determining the zonal mass fluxes across each 3D cell edge from the 2D solution. We compute the thermodynamic state on the full 3D mesh using the same iterative procedure used for the 2D mesh. The initial state computed for the jet using this procedure is depicted in Fig. 6; it is in approximate geostrophic balance and is nondivergent.

In this section we also present examples of full-physics NWP forecasts on a uniform mesh and a variable-resolution mesh. We present these forecasts to demonstrate the robustness of the solver at large scales with full physics and earth’s terrain; evaluation of the full-physics MPAS forecast capabilities on uniform and nonuniform meshes will be the subject of future reports.

### a. Dry baroclinic wave simulations

Figure 7 depicts the surface pressure, 850-hPa vorticity, and temperature for the perturbed jet at day 9 using the 60-km (163 842 cell) mesh and can be directly compared with the solutions from Jablonowski and Williamson (2006a) in their Figs. 7, 8. MPAS used a time step of 450 s in the Runge–Kutta solver and a 75-s time step for the acoustic mode integration. A second-order Smagorinsky filter was used with a filter coefficient *c _{s}* = 0.125 and a constant length scale

*l*= 6.0 × 10

^{4}m. Within the acoustic substeps of the time-split scheme we used a 3D divergence damping coefficient of

*β*= 0.1 [see Klemp et al. (2007), their Eq. (19)] and a vertically implicit off-centering parameter

_{d}*β*= 0.1 [see Klemp et al. (2007), their Eq. (17)]. We use the third-order transport scheme described in Skamarock and Gassmann (2011) for the thermodynamic equation (5). The 60-km MPAS mesh is of similar mesh resolution to those used in the Jablonowski and Williamson study where most of the results are presented for the nominally ½° meshes (~55 km). The MPAS results in Fig. 7 are qualitatively very similar to those reported in Jablonowski and Williamson, both with regard to the structure in the fields and the maximum and minimum surface pressures and vorticities. We are using a time step that is significantly larger than those used by the fully explicit models in the Jablonowski and Williamson study; the larger MPAS time step is made possible by the split-explicit time-integration scheme that handles the gravity waves on the acoustic substep. We have also used this test case to examine model behavior using different advection schemes and the results are reported in Skamarock and Gassmann (2011). The test results reported in Skamarock and Gassmann show that the phase errors are significantly reduced using the third-order scheme compared with the second-order scheme for potential temperature transport.

_{s}We have also conducted the Jablonowski and Williamson (2006a) test case without the wind field perturbation to test the ability of the solver to maintain the initial steady-state balanced jet, and the results are given in Fig. 8 for the 240-km (10 242 cell) mesh. These results can be directly compared with the results depicted in Fig. 8 of Lauritzen et al. (2010), where six different dynamical cores were tested. Results for three different mesh orientations are plotted in Fig. 8, one with pentagons of the icosahedral-based mesh set at the poles, and two with these pentagons rotated latitudinally 45° and 90° from the poles. The tests with no mesh rotation maintain the steady state best (smallest perturbations in surface pressure at day 9) because the pentagons on the icosahedral-based mesh are not located in the jet region. The 45° and 90° mesh rotations place the pentagons within the jet, and the perturbations in the surface pressure field are larger because the truncation errors in these model formulations are largest around the pentagons and excite the unstable modes of the jet more quickly. Comparing these MPAS simulations to those presented in Lauritzen et al. (2010), the perturbations in the zonal flow in the MPAS surface pressure field are somewhat smaller in magnitude than those of the hydrostatic Colorado State University (CSU) model results [using both the sigma and hybrid coordinates; see Lauritzen et al. (2010), their Fig. 8] and much lower in amplitude compared to the Icosahedral Nonhydrostatic model (ICON) results [which use a triangular primal mesh; also plotted in Lauritzen et al. (2010), their Fig. 8]. Additionally, the surface pressure fields exhibit symmetry about the equator in the MPAS results whereas there is no apparent symmetry in the CSU and ICON model results, suggesting either asymmetries in the CSU and ICON meshes or initialization, or errors in the solvers. We also wish to point out that models based on latitude–longitude meshes do not introduce any zonal perturbations, and no unstable modes can be excited, except when the computational poles are shifted relative to the geographic poles as shown in Lauritzen et al. (2010).

Based on the comparisons of these test results with published results from other models, we find that the nonhydrostatic MPAS produces results with accuracy equal to or greater than that of other hydrostatic icosahedral finite-volume formulations on the sphere.

### b. Unfiltered baroclinic wave simulations

The baroclinic wave simulations presented in the previous section used a spatial filter with characteristics similar to those employed by the models tested in Jablonowski and Williamson (2006a). Frontal collapse occurs between days 9 and 10 in the simulations, and solutions degenerate into grid-scale noise in the frontal regions shortly thereafter if no spatial filtering is used. The nonhydrostatic MPAS will, however, run stably through day 10 in the simulation without spatial filtering. Figure 9 shows the solution for vorticity at day 8 on model level 5 (approximately 850 hPa) for two simulations where all spatial filtering is disabled in MPAS. One simulation uses the kinetic energy formulation (13) with *α* = 1, corresponding to Ringler et al. (2010), and the second simulation uses this formulation with *α* = ⅜. The noise evident in the simulation using *α* = 1 causes the unfiltered model to be unstable; the model aborts shortly after day 9. Any of the spatial filters available in MPAS (the second- and fourth-order filters using constant eddy viscosities and hyperviscosities, and the Smagorinsky-based filter) are sufficient to remove this instability. It is preferable to remove this instability in the basic solver formulation, and the kinetic energy evaluation (13) with *α* = ⅜ appears to accomplish this for this test. We do not have a theoretical analysis explaining this instability. While it appears to share some of the characteristics of the instability described by Hollingsworth et al. (1983), the formulation (13) does not lead to a cancellation of terms in the nonlinear momentum transport described by Hollingsworth et al. Furthermore, we do not have a theoretical justification for specifying *α* = ⅜. We have found that the instability is not apparent for 0.20 ≤ *α* ≤ 0.45 in (13).

### c. Full-physics global forecast examples

We have performed full-physics multiday forecasts on both quasi-uniform and variable-resolution meshes to assess the robustness of the solver. We employ the following model physics taken directly from the Weather Research and Forecasting model physics suite (see Skamarock and Klemp 2008; Skamarock et al. 2008): WRF Single-Moment 6-Class Microphysics scheme (WSM6), Kain–Fritsch convective parameterization, Yonsei University (YSU) PBL parameterization, Monin–Obukhov surface layer parameterization, the Noah land surface model, and the Rapid Radiative Transfer Model–GCM applications (RRTMG) longwave and shortwave radiation parameterizations. Descriptions and references for these parameterizations can be found in the ARW-WRF Technical Note (Skamarock et al. 2008). We use the same time step on the entire mesh for the variable-resolution mesh simulations, and this time step is chosen so that it is stable for the fine-mesh region. We use the Smagorinsky filter in these forecasts with the coefficient *c _{s}* = 0.125. The length scale in the Smagorinsky filter is scaled with the local grid spacing.

We have computed 5-day forecasts valid at 0000 UTC 28 October 2010 using a uniform-resolution mesh with a mean cell-center spacing of 60 km and using a nonuniform mesh with a mean cell-center spacing varying between approximately 162 km in the coarse-mesh region and 21 km in the fine-mesh region. Both meshes use the same number of cells (163 842). The variable-resolution mesh is the same as that depicted in Fig. 10 except the cell mesh spacing is ¼ of that in the figure; it is the 163 842-cell × 8 mesh used in Ringler et al. (2011) and it is constructed using the iterative technique described in Ringler et al. (2008). We initialized the forecasts using the National Centers for Environmental Prediction (NCEP) Climate System Forecast Reanalysis (CSFR) (Saha et al. 2010), and the CSFR at day 5 is also presented in Fig. 11. We chose this time period because of the very strong surface low pressure region that formed over the north-central United States during this forecast and the strong North Pacific jet preceding it that extended through the western mesh-transition region in the variable-resolution forecast. The time steps are 225 s for the uniform nominally 60-km mesh and 90 s for the variable-resolution 162–21-km mesh.

The two forecasts and the analysis shown in Fig. 11 are similar; the midlatitude baroclinic waves evident in the height field have similar locations and amplitudes with some differences in structure, and the general precipitation patterns are similar. The forecasts are slightly warmer than the analysis at 500 hPa in the tropics and show slightly more precipitation at the higher thresholds in the tropics. The higher tropical precipitation amounts, in addition to the lighter grid-scale precipitation amounts in the tropics, are consistent with our experiences using the Kain–Fritsch convective parameterization in WRF at these resolutions.

The fine-mesh region of the nonuniform-mesh forecast has some finer-scale structure compared to the uniform-resolution forecast, although this is difficult to discern in these global plots. Figure 12 depicts the 500-hPa relative vorticity at day 4 for the uniform- and variable-resolution-mesh forecasts over a portion of the Northern Hemisphere that includes the fine-mesh and mesh-transition regions. The finer structure over North America is evident in the variable-mesh forecast, for example, in the increased wave activity over the western United States and the stronger wrapping of vorticity south of Alaska. There is a smooth transition of the relative vorticity field in the variable-mesh-transition regions (south and west of Alaska and southeast of Greenland; see Fig. 10). As expected, the coarse-mesh region of the variable-mesh solution is less well resolved compared to the uniform-mesh solution, as is evident in the vorticity filaments over Europe and northern Africa.

These preliminary results suggest that the uniform-mesh configurations appear viable for traditional global NWP purposes and that the variable-resolution meshes may be suitable for high-resolution NWP and regional climate. Further evaluation of the solutions produced in the variable-resolution regions is needed; we expect that model physics and dissipation operators will need tuning or perhaps significant reformulations for the variable-resolution mesh. We also note here that MPAS can be run as a regional model with time-dependent lateral boundary conditions provided from previous MPAS simulations (i.e., one-way nesting) or from other sources.

## 5. Summary

We have described the formulation of a fully compressible nonhydrostatic model whose discretization uses a horizontal (spherical) centroidal Voronoi mesh with a terrain-following geometric-height vertical coordinate and C-grid staggering for momentum. The Voronoi meshes are unstructured and permit variable horizontal resolution, and our nonhydrostatic model solves the equations of motion directly on these unstructured meshes. We employ a vector-invariant form of the horizontal momentum equation to avoid discretization difficulties on the unstructured mesh associated with the nonlinear momentum transport and to allow for potential vorticity conservation in the horizontal discretization (Ringler et al. 2010). The temporal discretization uses the explicit time-split Runge–Kutta technique from Wicker and Skamarock (2002) and Klemp et al. (2007). The potential benefits of this formulation for the compressible flow solver are made possible by three advances. First, we are making use of the C-grid discretization techniques for Voronoi (nominally hexagonal) meshes described by Thuburn et al. (2009) and Ringler et al. (2010) that solve the problems associated with the nonstationary geostrophic mode analyzed by Ničković et al. (2002). Second, we are using higher-order transport operators as described in Skamarock and Gassmann (2011); the higher-order transport scheme allows MPAS to produce solutions of similar accuracy to that of present-day state-of-the-art cloud and mesoscale models, and also improves the large-scale response in early test simulations of baroclinic waves. Finally, MPAS employs an unstructured mesh that permits continuous grid refinement and demonstrates computational efficiency similar to our rectangular-mesh formulations on existing cache-based supercomputer architectures.

We have presented idealized test results from the MPAS nonhydrostatic solver for large-scale (hydrostatic) flows [e.g., the Jablonowski and Williamson (2006a) baroclinic wave on the sphere] and for nonhydrostatic-scale flows (e.g., mountain waves, density currents, and supercell thunderstorms). MPAS has produced results comparable to other state-of-the-art models in these tests. Specifically, the results show that our use of the C-grid Voronoi mesh and the vector-invariant form of the horizontal momentum equation does not compromise the robustness or accuracy of the nonhydrostatic solver. We have also presented preliminary results using a full-physics NWP configuration for MPAS that demonstrates its potential for NWP and regional climate applications on both quasi-uniform and variable-resolution meshes. The horizontal meshes are designed such that there is a smooth transition in mesh density between the coarse-and fine-resolution regions of the mesh. We have not observed any significant deleterious effects associated with the mesh transition in the solutions we have presented here or in other solutions we have computed with MPAS. As expected, the fine-mesh regions show finer-scale structure compared to the structure in the coarse-mesh region. Further research is needed to explore mesh refinement limitations and optimal mesh-transition characteristics (e.g., transition zone width and mesh resolution rate of change) for these Voronoi meshes and our nonhydrostatic solver. We will also be exploring the use of different time steps in regions of different mesh resolutions.

In our example variable-resolution forecast, the same model physics and subgrid-scale parameterizations were used over the entire mesh. The model physics and subgrid-scale parameterizations will need to be scale-aware in applications with widely varying mesh densities, particularly when mesh densities transition between resolving hydrostatic-scale and nonhydrostatic-scale motions.

## REFERENCES

*The Ceaseless Wind: An Introduction to the Theory of Atmospheric Motion.*Dover Publications, 617 pp.

*Mon. Wea. Rev.,*

**139,**2163–2169.

*J. Adv. Model. Earth Syst.,*

**2,**

**25,**392–403,

**139,**2962–2975.

## Footnotes

The National Center for Atmospheric Research is sponsored by the National Science Foundation.