Verification of a non-hydrostatic dynamical core using the horizontal spectral element method and vertical finite difference method: 2-D aspects

Abstract. The non-hydrostatic (NH) compressible Euler equations for dry atmosphere were solved in a simplified two-dimensional (2-D) slice framework employing a spectral element method (SEM) for the horizontal discretization and a finite difference method (FDM) for the vertical discretization. By using horizontal SEM, which decomposes the physical domain into smaller pieces with a small communication stencil, a high level of scalability can be achieved. By using vertical FDM, an easy method for coupling the dynamics and existing physics packages can be provided. The SEM uses high-order nodal basis functions associated with Lagrange polynomials based on Gauss–Lobatto–Legendre (GLL) quadrature points. The FDM employs a third-order upwind-biased scheme for the vertical flux terms and a centered finite difference scheme for the vertical derivative and integral terms. For temporal integration, a time-split, third-order Runge–Kutta (RK3) integration technique was applied. The Euler equations that were used here are in flux form based on the hydrostatic pressure vertical coordinate. The equations are the same as those used in the Weather Research and Forecasting (WRF) model, but a hybrid sigma–pressure vertical coordinate was implemented in this model. We validated the model by conducting the widely used standard tests: linear hydrostatic mountain wave, tracer advection, and gravity wave over the Schar-type mountain, as well as density current, inertia–gravity wave, and rising thermal bubble. The results from these tests demonstrated that the model using the horizontal SEM and the vertical FDM is accurate and robust provided sufficient diffusion is applied. The results with various horizontal resolutions also showed convergence of second-order accuracy due to the accuracy of the time integration scheme and that of the vertical direction, although high-order basis functions were used in the horizontal. By using the 2-D slice model, we effectively showed that the combined spatial discretization method of the spectral element and finite difference methods in the horizontal and vertical directions, respectively, offers a viable method for development of an NH dynamical core.


Introduction
There is growing interest in developing highly scalable dynamical cores using numerical algorithms under petascale computers with many cores (with the goal of exascale computing just around the corner), and the spectral element method (SEM), with high efficiency and accuracy, is known to be one of the most promising methods (Taylor et al., 1997;Giraldo, 2001;Thomas and Loft, 2002).SEM is local in nature because it has a large on-processor operation count (Kelly and Giraldo, 2012).SEM achieves this high level of scalability by decomposing the physical domain into smaller pieces with a small communication stencil.Additionally, SEM has been shown to be very attractive for achieving highorder accuracy and geometrical flexibility on the sphere (Taylor et al., 1997;Giraldo, 2001;Giraldo and Rosmond, 2004).
To date, SEM has been implemented successfully in atmospheric modeling, such as in the community atmosphere model-spectral element (CAM-SE) dynamical core (Thomas Published by Copernicus Publications on behalf of the European Geosciences Union. S.-J.Choi et al.: Verification of a non-hydrostatic dynamical core and Loft, 2005) and the scalable spectral element Eulerian atmospheric model (SEE-AM) (Giraldo and Rosmond, 2004).These models consider the primitive hydrostatic equations on global grids, such as a cubed sphere tiled with quadrilateral elements using SEM in the horizontal discretization and the finite difference method (FDM) in the vertical.The robustness of SEM has been illustrated through three-dimensional dry dynamical test cases (Giraldo and Rosmond, 2004;Giraldo, 2005;Thomas and Loft, 2005;Taylor et al., 2007;Lauritzen et al., 2010).
The ultimate objective of our study is to build a 3-D non-hydrostatic (NH) model based on the compressible Navier-Stokes equations using SEM in the horizontal discretization and FDM in the vertical.Because testing a 3-D NH model requires a large amount of computing resources, studying the feasibility of our approach in 2-D is an attractive alternative to the development of a fully 3-D model.This is the case because a 2-D slice model can effectively test the practical issues resulting from the vertical discretization and time integration prior to construction of a full 3-D model.Although we could discretize the vertical direction using SEM (as proposed in Kelly andGiraldo, 2012, andGiraldo et al., 2013), we chose to use a finite difference method for discretization in the vertical direction because it provides an easy way to couple the dynamics and existing physics packages.
For this objective, we developed a dry 2-D NH compressible Euler model based on SEM along the x direction and FDM along the z direction, which we hereafter refer to as the 2-D NH model.We adopted the governing equation formulation proposed by Skamarock and Klemp (2008) (SK08 hereafter), which is used in the Weather Research and Forecasting (WRF) model.The Euler equations are in flux form based on the hydrostatic pressure vertical coordinate.In SK08, the terrain-following sigma-pressure coordinate is used, but we here employed a hybrid sigma-pressure vertical coordinate.Park et al. (2013) (PK13 hereafter) provided a clue for the equation set in the hybrid sigma-pressure in their Appendix, in which the hybrid sigma-pressure coordinate is applied to the hydrostatic primitive equations and can be modified exactly to the sigma-pressure coordinate at the level of the actual coding implementation.We also built the 2-D NH model using a time-split, third-order Runge-Kutta (RK3) for the time discretization, which has been shown to be effective in the WRF model.We kept the temporal discretization of the model as similar as possible to the WRF model in order to more directly discern the differences related to the discrete spatial operators between the two models.This provides robust tools for development and verification of the 2-D NH model.
In this paper, we demonstrate the feasibility of the 2-D NH model by conducting conventional benchmark test cases and by focusing on the description of the numerical scheme for the spatial discretization.We verify the 2-D NH by analyzing six test cases: inertia-gravity wave, rising thermal bub-ble, density current wave, linear hydrostatic mountain wave, and tracer advection and gravity wave over the Schär-type mountain.
The organization of this paper is as follows.In the next section, we describe the governing equations with definitions of the prognostic and diagnostic variables used in our model.In Sect.3, we explain the temporal and spatial discretization including the spectral element formulation.In Sect.4, we present the results of the 2-D NH model using all four test cases, and finally, in Sect.5, we summarize the paper and propose future directions.

Governing equations
We adopted the formulation of the governing-equation set of SK08.Here, we implemented the hybrid sigma-pressure coordinate introduced in PK13, which only considers the hydrostatic primitive equation.The hybrid sigma pressure coordinate is defined with η ∈ [0, 1] as where p d is the hydrostatic pressure of dry air; B(η) is the relative weighting of the terrain-following coordinate versus the normalized pressure coordinate; and p s , p t , and p 0 are the hydrostatic surface pressure of dry air, the top-level pressure, and a reference sea level pressure, respectively.A more detailed description of the hybrid sigma-pressure coordinate can be found in the Appendix of PK13.The definition of the flux variables are where v H = (u, v) and w are the velocities in the horizontal and vertical directions, respectively; η ≡ dη dt is the ηcoordinate (contravariant) vertical velocity; θ is the potential temperature; and µ d is the mass of the dry air in the layers, defined as (3) The flux-form Euler equations for dry atmosphere to be recast using perturbation variables are expressed as where φ is the geopotential; α d is the inverse density for dry air; and F V H and F W represent forcing terms of Coriolis and curvature, which we ignore for simplicity.In Eqs. ( 4)-( 8), the governing equations are described with perturbation variables, such as p = p(z)+p , φ = φ(z)+φ , α d = ᾱd (z)+α d , and p s = ps (x, y) + p s , where the variables denoted by an overbar are the reference state variables that satisfy hydrostatic balance.
For completeness, the diagnostic relation for is given by integrating Eq. ( 6) vertically from the surface (η = 1) to the material surface: where ∂t is obtained by integrating Eq. ( 6) vertically from the surface (η = 1) to the top (η = 0) using a no-flux boundary condition, such as | η=0 or 1 = 0 .∂p s ∂t is defined as The above equation allows for p s to be evolved forward in time, where we then compute µ d directly from Eq. ( 5).The diagnostic relation for the dry inverse density is given as and the full pressure for dry atmosphere is This concludes the description of the governing equations used in our model; in the next section, we describe the discretization of the continuous form of the governing equations that are used in our model.

Horizontal direction
For a given η level, we discretized the horizontal operators using SEM.Therefore, in the 2-D (x − z) slice framework, we focus on the SEM discrete gradient operator for 1-D (x).In SEM, we approximate the solution in non-overlapping elements e as q(x, t) where x i represents the N + 1 grid points that correspond to the Gauss-Lobatto-Legendre (GLL) points and ψ i (x) are the N th-order Lagrange polynomials based on the GLL points.
It is worth noting that the ψ i have the cardinal property, i.e., they can be represented as Kronecker delta functions where ψ i are zero at all nodal points except x i .
The GLL points ξ i in a reference coordinate system ξ ∈ [−1, +1] and the associated quadrature weights ω(ξ i ), are introduced for the Gaussian quadrature: e where P N (ξ ) are the Nth-order Legendre polynomials, J = ∂x ∂ξ is the transformation Jacobian, and e represents the non-overlapping elements.
We now introduce the polynomial expansions into our governing equations in the form of multiply by the basis function ψ i as a test function, and integrate to yield a system of ordinary differential equations, such as where i = 1, 2, • • •, N +1, M e j i is the element-based mass matrix given as The right-hand sides of Eqs. ( 17) and (18) are evaluated using the Gaussian quadrature of Eq. (15).It is noted that using GLL points for both interpolation and integration results in a diagonal mass matrix M e j i , which means that the inversion of the mass matrix is trivial.
The horizontal derivatives included in the right-hand side of Eq. ( 17) are evaluated using the analytic derivatives of the basis functions as follows: Note that the non-differential operations, such as cross products, are computed directly at grid points since we use nodal basis functions associated with Lagrange polynomials based on the GLL points.In order to satisfy the equations globally, we use the direct stiffness summation (DSS) operation.For a more detailed description of the SEM, see Giraldo and Rosmond (2004), Giraldo and Restelli (2008), and Kelly and Giraldo (2012).

Vertical direction
Using a Lorenz staggering, the variables V H , , µ, α, and p are at layer midpoints denoted by k = 1, 2, . .., K, where K is the total number of layers, and the variables W , , and φ are at layer interfaces defined by k + 1 2 , k = 0, 1, . .., K, so that η K+1/2 = η top and η 1/2 = η Bottom = 1. Figure 1 describes the grid points and the allocation of the variables.Here, we evaluate the vertical advection terms ∂( v H ) ∂η , ∂( w)  ∂η , and ∂( θ ) ∂η and vertical derivative terms ∂p ∂η and ∂φ ∂η .The former is discretized using the third-order upwind-biased discretization in Hundsdorfer et al. (1995), which is given as where f corresponds to the flux, such as v H , and η = η k+1/2 − η k−1/2 is the thickness of the layer.The latter is discretized by the centered finite difference, which is given as where g corresponds to the variables p and ϕ.Likewise, the vertical discretization integration rules for the calculations of Eqs. ( 9) and (10) follow the finite difference naturally as

Explicit diffusion
In addition to the governing equations, a viscous term might be needed to conduct some tests.The viscosity used here is an explicit Laplacian (∇ 2 ) diffusion operator on coordinate surfaces.In order to implement the Laplacian operator for a model flux variable µ d a, we multiply by the basis function ψ as a test function and integrate using the divergence theorem to yield the weak form equation where ν h denotes the eddy viscosity coefficient and the term with e is a boundary integral that accounts for internal faces (neighboring elements share faces).Because we used the periodic boundary condition in this study, the boundary integral term of the right-hand side can be ignored in all elements, which allows us to rewrite the equations as After introducing the polynomial expansions, such as a(x, t) = N+1 i=1 ψ i (x)a N (x i , t), the integrals of the above equation can be approximated using SEM.A description of the Laplacian operator using SEM can also be found in Denis et al. (2011).The vertical Laplacian operator for a model flux variable µ d a is added to a governing equation as follows: where ν v denotes the vertical eddy viscosity coefficient and α is the inverse density.It is noted that the above term is not more than ν v ∂z 2 .The vertical derivative term ∂ ∂η is discretized by the centered finite difference.

Temporal discretization
To integrate the equations, we used the time-split RK3 integration technique following the strategy of SK08.In the time-split RK3 integration, low-frequency modes due to advective forcings are explicitly advanced using a large time step in the RK3 scheme, but high-frequency modes are integrated over smaller time steps.Among the high-frequency modes, horizontally propagating acoustic/gravity waves are advanced using an explicit forward-backward time integration scheme and vertically propagating acoustic waves and buoyancy oscillations are advanced using a fully implicit scheme (Klemp et al., 2007).For numeric stability, acoustic-mode filterings of the forward centering of the vertically implicit portion and divergence damping of the horizontal momentum equation are used, which is the same as in the WRF model (Skamarock et al., 2008).It is notable that the timesplit RK3 integration scheme is third-order accurate for linear equations and second-order accurate for nonlinear equations (SK08).
This technique has been shown to work effectively within numerous non-hydrostatic models, including the WRF model (Skamarock et al., 2008), the Model for Prediction Across Scales (MPAS) (Skamarock et al., 2012), and the Non-hydrostatic Icosahedral Atmospheric Model (NICAM) (Satoh et al., 2008).It is also noted that, in the procedure of the time-split RK3 integration, the difference between the approach used in this paper and that in SK08 comes from the vertical coordinate.Since we use the hybrid sigma-pressure coordinate, the equation for p s (Eq.6) should be first stepped forward in time using forward-backward differencing on the small time steps, then µ d can be computed directly from the specification of the vertical coordinate in Eq. ( 9) and can be obtained from the vertical integration.

Test cases
We validated the 2-D NH model with six test cases: linear hydrostatic mountain-wave, tracer-advection, and gravity-wave tests over Schär Mountain, as well as density current, inertia-gravity wave, and rising thermal bubble experiments.The last three cases do not have analytic solutions.Therefore, for the mountain experiments, the numerical results of the 2-D NH model were compared with analytic solutions (Durran and Klemp, 1983;Schär et al., 2002); for the other experiments, we compared our results with the results of other published papers.
It should be mentioned that the horizontal SEM formulation is able to utilize arbitrary-order polynomials per element to represent the discrete spatial operators, but in this paper all the results presented use either fifth-or eighth-order polynomials.The averaged horizontal grid spacing is defined as where x n is the internal grid spacing within the element, which is regularly spaced in the domain, and N is the number of intervals associated with irregularly spaced GLL quadrature points, which is equivalent to the order of the basis polynomials.The average vertical grid spacing is defined as in Eq. ( 26).Below, we use this convention to define the grid resolution.The resolutions and time-step sizes used for all the cases are summarized in Table 1.

Linear hydrostatic mountain-wave test
We simulated the linear hydrostatic mountain-wave test introduced by Durran and Klemp (1983) (DK83 hereafter) in which the analytic steady-state solution is provided by using a single-peak mountain with uniform zonal wind.To compare our results with the analytic and numerical solutions shown in DK83, the 2-D NH was initialized using the same initial conditions and mountain profile as in DK83, and we analyzed our results using the same metrics as DK83.
The mountain profile is given by where the half-length of the mountain a m is 10 km, the height h m is 1 m, and the prescribed center x c of the profile is 0 km.The initial temperature is T 0 = 250 K for an isothermal atmosphere with the uniform zonal wind ū = 20 m s −1 .In the isothermal case, the Brunt-Väisälä frequency N 2 = g d(ln θ) dz ≈ g 2 c p T 0 yields the potential temperature as which is one of the prognostic variables in our model.The domain is defined as The bottom boundary uses a no-flux boundary condition, whereas the lateral and top boundaries use sponge layers.The sponged zone is 10 km deep from the top and 50 km wide from the lateral boundaries.Over the sponge-layer zone, the prognostic variables are relaxed to the basic initial hydrostatic state.The model is integrated with a grid resolution of x = 2 km using fifth-order basis polynomials per element and z = 375 m for a nondimensional time of ūt a = 60, which corresponds to 8.33 h.Additionally, the model is run without diffusion or viscosity.
Figure 2 shows the numerical and analytic solutions at steady state for the horizontal and vertical velocities, which agree reasonably well.The vertical velocity fields match very closely, although the extrema in the horizontal velocity field are underestimated by the numerical model.The underestimated extrema in the horizontal velocity were also shown in both models of DK83, which used x = 2 km and z = 200 m, and in Giraldo and Restelli (2008) (GR08 hereafter), which used x = 1.2 km and z = 240 m with 10thorder basis polynomials.Our result in the horizontal velocity is in good agreement with DK83 and GR08.
To check the vertical transport of horizontal momentum, Fig. 3 shows the normalized momentum flux values at various times.It is observed that the flux has developed well and that the simulations reach steady state after ūt a = 60.It is also noted that the mean momentum flux at this time is 97 % of its analytic value.The result agrees well with DK83 as well as GR08; however, it is important to point out that the  Rising thermal bubble (fifth-order basis function) x = 20 and z = 20 x = 10 and z = 10 x = 5 and z = 5 Durran-Klemp model is based on the FD method in both directions, while the Giraldo-Restelli model is based on SEM in both directions.The mountain test shows that the terrainfollowing vertical coordinate is well suited for the combination of horizontal SEM and vertical FDM for spatial discretization, even though we considered a small mountain.

Tracer-advection and gravity-wave tests over the Schär-type mountain
In order to verify the feasibility of 2-D NH to treat steep surface elevations associated with the vertical terrain-following coordinate, we performed the tracer-advection and gravitywave experiments introduced by Schär et al. ( 2002) (SC02 hereafter), in which the mountain is defined by a five-peak mountain chain, over the Schär-type mountain.To compare our results with the numerical solution shown in SC02, the initial conditions and mountain profiles are the same as those of SC02.
For the tracer-advection test, the mountain profile is given by where h 0 = 3 km, a = 25 km, and λ = 8 km.The prescribed wind profile is given by where u 0 = 10 m s −1 , z 1 = 4 km, and z 2 = 5 km, and the initial tracer is assigned as , where amplitude q 0 = 1, location (x 0 , z 0 ) = (−50, 9) km, and the half-width (A x , A z ) = (25, 3) km.Since the domain is defined as (x, z) ∈ [−150, 150] × [0, 25] km 2 , the tracer is centered directly over the mountain at time 2500 s.The model is integrated with a grid resolution of x = 300 m using fifth-order basis polynomials per element and z = 250 m using 100 levels for 5000 s.The model is run without any diffusion, filter, or limiter.It should be noted that the advection equation used in this study is the advective form defined as  The numerical solutions and the error field are shown in Fig. 4. The figure uses the same contouring interval as in SC02.Even at t = 2500 s, when the center of the tracer is located over the center of the mountain, the distribution of the initial tracer is generally maintained (Fig. 4a), which means that 2-D NH using the horizontal spectral element method and vertical finite difference method can produce numerical solutions of good quality in response to the strong vertical gradient in the coordinate deformation.It is worth noting that the error in Fig. 4b at t = 5000 s gives ranges of −2.71 × 10 −2 , 2.35 × 10 −2 , which are substantially small, and that the error is distributed mainly over the mountain where distortion of the computational grid is significant.The Schär-type mountain gravity-wave test was initialized in a stratified atmosphere with the Brunt-Väisälä frequency of N = 0.01 s −1 , the constant mean flow of ū = 10 m s −1 , and the initial temperature of T 0 = 288 K.In the Schär-type mountain gravity wave, the highest mountain peak was h 0 = 250 m, which is relatively lower than that in the advection test.The mountain profile is given by where a = 5 km and λ = 4 km, and the domain is defined as (x, z) ∈ [−30, 30] × [0, 21] km 2 .The model was integrated with a grid resolution of x = 300 m using fifth-order basis polynomials per element and z = 250 m using 80 levels for 10 h without any diffusion or viscosity.The bottom boundary had a no-flux boundary condition, while the lateral and top boundaries had sponge layers.The sponged zone was 10 km deep from the top and 5 km wide from the lateral boundaries.
Over the sponge layer zone, the prognostic variables were relaxed to the initial state.
Figure 5 shows the simulated results of the perturbed horizontal and vertical wind speeds after 10 h.In comparison with the analytic solution, the numerical solutions match quite well.The results of the present model are also very similar to the results of other numerical models (Giraldo and Restelli, 2008;Li et al., 2013).For a quantitative comparison, we present the root-mean-square errors for u , w , and θ in Table 2.These values are very comparable with those of other numerical models (Giraldo and Restelli, 2008;Li et al., 2013).Table 2. Root-mean-square errors (RMSEs) of the Schär-type mountain wave after 10 h for x = 300 m using fifth-order polynomials per element and z = 250 m using 80 levels.
Variable RMSE

2-D density current test
In order to verify the feasibility of 2-D NH to control oscillations with numerical viscosity and evaluate numerical schemes in 2-D NH, we conducted the density current test suggested by Straka et al. (1993).The density current test is initialized using a cold bubble in a neutrally stratified atmosphere.When the bubble touches the ground, the density current wave starts to spread symmetrically in the horizon- tal direction, forming Kelvin-Helmholtz rotors.Following Straka et al. (1993), we employed a dynamic viscosity of ν = 75 m 2 s −1 to obtain converged numerical solutions.The viscosity used here is an explicit Laplacian (∇ 2 ) diffusion operator on coordinate surfaces.For an initial cold bubble, the potential temperature perturbation is given as where θ c = −15 K and r = Figure 6 shows the potential temperature perturbation after 900 s for 400, 200, 100, and 50 m grid spacings ( x) using fifth-order basis polynomials per element.All simulations used z = 64 m grid spacing vertically.As expected, the higher-resolution experiments produced better solutions than the lower-resolution experiments.At the very lowest resolution of 400 m, only two of the three Kelvin-Helmholtz rotors were generated with somewhat coarsened frontal surfaces.In the experiment with a resolution of 200 m, the three rotors appeared, but the numerical solution still suffered from the coarsening of frontal surfaces.The solutions on grids finer than 100 m converged with the three rotor structures adequately simulated.The converged solution was almost identical to other published solutions (e.g., Straka et al., 1993;Skamarock and Klemp, 2008;GR08).
In order to examine the effect of higher order of the basis polynomials than fifth-order, we show profiles of the potential temperature perturbation at the height of 1200 m in the simulations using fifth-order polynomials together with the simulations using eighth-order polynomials (Fig. 7).Note that the simulations using eighth-order polynomials have the same number of GLL grid points as the simulations using fifth-order basis polynomials.This was achieved by using a lower number of elements in the eighth-order experiment than in the fifth-order experiment as the number of grid points at a given level becomes n e × n p , in which n e refers to the number of elements and n p denotes the polynomial order of the elements.It is also noted that we arbitrarily choose eighth-order as the higher order.The results from the highest grid resolution of the simulations using fifth-and eighthorder polynomials are indistinguishable and well converged, with three minima corresponding to the three rotors, which agree well with other published solutions (Fig. 7a).In addition to the profiles, the front location (−1 K of potential temperature perturbation at the surface) and the extrema of the pressure perturbation and potential temperature perturbation agreed well with each other (Table 3).The numbers in Table 3 are comparable to those of GR08.Although the potential temperature profiles of the simulations using fifth-order polynomials tend to have more fluctuations than those using eighth-order basis polynomials, the simulations do not show a large difference between using eighth-order and fifth-order basis polynomials (Fig. 7b and c).
In order to investigate the characteristics of the convergence more clearly, a self-convergence test was carried out.were interpolated to the equidistant grid of x = 400 and z = 50 and then used to evaluate errors.Here, we evaluated the error by using the relative L2 error defined by where q simulation and q ref represent the model solution and reference solution, respectively.The resulting L2 norm of the error in the potential temperature perturbation θ is plotted in Fig. 8.At the highest resolution of x = 50 m, it is noted that the experimental convergence rate reaches the convergence rate 2, which depends on the accuracy of the time-split RK3 integration scheme in relatively uniform spacing in the vertical direction.Note that it could be theoretically first-order accuracy with resolution if fully non-uniform vertical spacing were used, since the centered difference scheme in the vertical direction is implemented.Additionally, it is shown that the error of the solutions of the eighth-order basis function is slightly smaller than that of the fifth-order basis function.

Inertia-gravity-wave test
This test examines the evolution of a potential temperature perturbation θ in a constant mean flow with a stratified atmosphere.This initial potential temperature perturbation θ radiates symmetrically to the left and right in a channel with periodic lateral boundary conditions.The inertia-gravity-wave test introduced by Skamarock and Klemp (1994) (SK94 hereafter) serves as a tool to investigate the accuracy for NH dynamics.We also used this experiment to check the consistency of the results at various resolutions.The parameters for the test were the same as those of SK94.The initial state was a constant Brunt-Väisälä frequency of N = 0.01 s −1 with a surface potential temperature of θ 0 = 300 K and a uniform zonal wind of ū = 20 m s −1 .In order to trigger the wave, the initial potential temperature perturbation θ was overlaid the above initial state and is given as where θ c = 0.01 K, z c = 10 km, x c = 100 km, and a c = 5 km.The domain was defined as (x, z) ∈ [0, 300] × [0, 10] km 2 .We used periodic lateral boundary conditions and no-flux boundary conditions for both the bottom and top boundaries.
The simulation was performed for 3000 s with no viscosity.
Figure 9 shows the solution θ at the initial time and at time 3000 s with horizontal resolution x = 250 m and vertical resolution z = 250 m.For comparison, the figure uses the same contouring interval as in SK94 and Giraldo and Restelli (2008).The results were produced with eighth-order polynomials per element.We conducted the 2-D NH model with various basis polynomial orders at the same resolution, and the simulated results were found to be very comparable.SK94 provides an analytic solution for the case of the Boussinesq equations; however, it is only valid for the Boussinesq equations and we used the fully compressible equations in our model.Using the analytic solution for only qualitative comparisons, we found that the extrema of our results are comparable to the analytic values.Compared with the results of Giraldo and Restelli (2008), for which the fully compressible equations were also used, our results appear very similar.
Figure 10 shows profiles along 5000 m for various horizontal resolutions.All models show consistently identical solutions with symmetric distribution about the midpoint (x = 160 km), which is the location to which the initial perturbation moved by the horizontal flow of 20 m s −1 after 3000 s.Even in coarser-resolution experiments, it does not exhibit phase errors, although the maxima and minima near the midpoint (x = 160 km) are slightly damped.Table 4 shows the extrema of vertical velocities and potential temperature perturbations for the results of various horizontal resolutions after 3000 s.All the experiments give almost the same values for potential temperature perturbation, which is  in the range θ ∈ −1.52 × 10 −3 , 2.83 × 10 −3 .These values are comparable to those of other studies.For example, GR08 gave the ranges of θ ∈ −1.51 × 10 −3 , 2.78 × 10 −3 from the model based on the spectral element and discontinuous Galerkin methods.Additionally, Li et al. (2013), using the high-order conservative finite volume model, showed θ ∈ −1.53 × 10 −3 , 2.80 × 10 −3 .

Rising thermal bubble test
We also conducted the rising thermal bubble test to verify the consistency of the scheme in the model to simulate thermodynamic motion (Wicker and Skamarock, 1998).This test considers the time evolution of warm air in a constant potential temperature environment for an atmosphere at rest.The air that is warmer than ambient air rises due to buoyant forcing, which then deforms due to the shearing motion caused by gradients of the velocity field and eventually shapes the thermal bubble into a mushroom cloud.Because the test case has no analytic solution, the simulation results were evaluated qualitatively.
The initial conditions we used follow those of GR08, in which the domain for the case is defined as (x, z) ∈ [0, 1] 2 km 2 .We used no-flux boundary conditions for all four boundaries.The domain was initialized for a neutral atmosphere at rest with θ 0 = 300 K in hydrostatic balance.The potential temperature perturbation to drive the motion is given as where θ c = 0.5 K, r = (x − x c ) 2 + (z − z c ) 2 with (x c , z c ) = (500, 350) m, and r c = 250 m.The model was run for a time of 700 s.It should be noted that an explicit Laplacian (∇ 2 ) diffusion on coordinate surfaces was used with a viscosity coefficient of ν = 1 m 2 s −1 for all simulations of this test.The numerical diffusion was applied for momentum and potential temperature along the horizontal and vertical directions to eliminate erroneous oscillations at the small scale.Although this amount of diffusion might seem excessive, it was chosen because it allows the model to remain stable even after the bubble reaches the top boundary.
Figure 11 shows the potential temperature perturbation, horizontal wind field, and vertical wind field for the simulations of the two resolutions of 20 and 5 m horizontal and vertical grid spacings ( x and z), respectively, employing fifth-order basis polynomials.In both simulations, the fine structures in the numerical solutions are well depicted, with a symmetric distribution at the midpoint and sharp discontinuities of the fields along the boundary lines of the bubble.At lower resolution, however, degradations in the solution are visible in the potential temperature perturbation and vertical wind, as illustrated by fluctuations in the values as well as the concave lines at the top of the bubble.It is noted that although the numerical solution of the model using the spatially centered FDM of Wicker and Skamarock (1998) shows spurious oscillations in the potential temperature field, the present simulations of 2-D NH using SEM horizontally and FDM vertically is devoid of these oscillations.We also show the vertical profiles of potential perturbation at x = 500 m after 700 s for various resolutions in Fig. 12. Simulations were run with the resolutions of 5, 10, and 20 m, where the resolutions given are defined for both the horizontal and vertical directions.The results of the 10 and 5 m resolutions are almost identical.The result of the lowest 20 m resolution, however, shows a somewhat unresolved solution, in which the maximum value is underestimated and the phase shift is depicted.Time series for maximum potential temperature perturbation and maximum vertical velocity are shown x z ∆ = m (thin dashed line), and , 5 x z ∆ ∆ =m (thick solid line). in Fig. 13.In all simulations, the maximum vertical velocity increases as the maximum theta perturbation decreases.This shows that the thermal energy of the theta perturbation leads to the acceleration of the vertical velocity.This result agrees well with the study of Ahmad and Lindeman (2007).

Summary and conclusions
The non-hydrostatic compressible Euler equations for a dry atmosphere were solved in a simplified 2-D slice (X-Z) framework by using spectral element method (SEM) for the horizontal discretization and finite difference method (FDM) for the vertical discretization.The form of the Euler equations used here is the same as those used in the Weather Research and Forecasting (WRF) model.We employed a hybrid sigma-pressure vertical coordinate, which can be converted exactly into a sigma-pressure coordinate at the level of the actual coding implementation.
For the spatial discretization, the spatial operators were separated into their horizontal and vertical components.In the horizontal components, the operators were discretized using SEM, in which high-order representations are constructed through the GLL grid points by Lagrange interpolations in elements.Using GLL points for both interpolation and integration results in a diagonal mass matrix, which means that the inversion of the mass matrix is trivial.In the vertical components, the operators were discretized using the third-order upwind-biased finite difference scheme for the vertical fluxes and centered differences for the ver- We presented results from idealized standard benchmark tests for large-scale flows (e.g., mountain-wave tests) and for non-hydrostatic-scale flows (e.g., inertia-gravity wave, rising thermal bubble, and density current).By varying the viscosity between test cases, the numerical results showed that the present dynamical core is able to produce high-quality solutions comparable to other published solutions.These tests effectively revealed that the combined spatial discretization method of the spectral element and finite difference methods in the horizontal and vertical directions, respectively, offers a viable method for the development of a NH dynamical core.Further work will be needed to achieve accurate solutions for a resting atmosphere over steep orography with minimal diffusion and to implement a horizontal diffusion operator in physical space, although horizontal diffusion on the coordinate surface was used in this study.Further research will also be conducted to couple the present core with the existing physics packages and extend the 2-D slice framework to develop a 3-D dynamical core for the global atmosphere in which the cubed-sphere grid is used for the spherical geometry.

Figure 1 .
Figure 1.Grid points of columns within an element having four GLL points.The hybrid sigma-pressure coordinates are illustrated, and the closed (open) circles on the solid (dashed) line indicate the location of the variables at layer midpoints (interfaces).
density current (fifth-and eighth-order basis function) x = 400 and z = 64 x = 200 and z = 64 x = 100 and z = 64 x = 50 and z = 64 t = 0.3 Inertia-gravity wave (eighth-order basis function) x = 1250 and z = 250 x = 500 and z = 250 x = 250 and z = 250 x = 125 and z = 250 t = 1 -state flow of (top) horizontal velocity (m/s) and (bottom) vertical velocity r 1-m high mountain at nondimensional time 60 ut a = with a grid resolution of m using 5th-order basis polynomials per element and 375 z ∆ = m.The numerical represented by solid lines and the analytic solution is represented by dashed lines.

Figure 2 .
Figure2.Steady-state flow of (top) horizontal velocity (m s −1 ) and (bottom) vertical velocity (m s −1 ) over 1 m high mountain at nondimensional time ūt a = 60 with a grid resolution of x = 2 km using fifth-order basis polynomials per element and z = 375 m.The numerical solution is represented by solid lines and the analytic solution is represented by dashed lines.

1FIG. 3 .Figure 3 .
FIG. 3. Vertical flux of horizontal momentum, normalized by its analytic value at several 2 nondimensional times ut a .M and M H are the momentum flux of the numerical and analytic 3 solutions.4 5

Figure 4 .
FIG. 4. Tracer advection test over the topography (red line).(a) Advective tracer at time 2 0 (black line), 2500 s (orange), and 5000 s (blue).The contour values are from −1.0 to 1.0 3 with an interval of 0.1.(b) Error at time 5000 s.The contour values are from −0.24 to 0.2 with 4an interval of 0.01.The numerical solutions were obtained with a grid resolution of 5 300 x ∆ = m using 5th-order basis polynomials per element and

Figure 6 .
Figure 6.Potential temperature perturbation after 900 s using grid spacing of (a) x = 400 m, (b) x = 200 m, (c) x = 100 m, and (d) x = 50 m, with fifth-order basis polynomials per element for the density current.All simulations use z = 64 m grid spacing.The contour values are from −14.5 to −0.5 with an interval of 1.0.
bubble at (x c , z c ) = (0, 3000) m and the size parameter (x r , z r ) = (4000, 2000) m.No-flux boundary conditions were used for all boundaries, and the model was integrated for 900 s on the domain [−25 600, 25 600] × [0, 6400] m 2 .In this study, the potential temperature perturbation of θ c = −15 K was adopted for comparison with GR08 andLi et al. (2013).Straka et al. (1993) originally used a −15 K temperature perturbation.The −15 K potential temperature corresponds to −13.53 K temperature.

Figure 7 . 6 Figure 8 .
FIG. 7. Profiles of (a) potential temperature perturbation after 900 s along 1200 m height ing grid spacing of

6 Figure 9 .
FIG. 9. Potential temperature perturbation at the initial time (top) and time 3000 s 2 (bottom) for

Figure 11 .
FIG. 11.Plots of (a, b) potential temperature perturbation (K), (c, d) horizontal wind 2 (m/s), and (e, f) vertical wind (m/s) for the rising thermal bubble test after 700 s with (left) 3 , 20 x z ∆ ∆ = m and (right) , 5 x z ∆ ∆ =m resolution for the rising thermal bubble test.All 4 simulations use 5th-order basis polynomials per element.Negative values are denoted by 5 dashed lines and positive values are denoted by solid lines.6 7

Table 4 .
FIG. 12. Vertical profiles of the potential temperature perturbation for the rising thermal bble test at x = 500 m after 700 s for various resolutions:, 20x z ∆ ∆ = m (thin solid line),

Figure 12 .
Figure 12.Vertical profiles of the potential temperature perturbation for the rising thermal bubble test at x = 500 m after 700 s for various resolutions: x, z = 20 m (thin solid line), x, z = 10 m (thin dashed line), and x, z = 5 m (thick solid line).

5 Figure 13 .
FIG. 13.Domain maximum potential temperature perturbation (top) and vertical wind 2 (bottom) for the rising thermal bubble test.All simulations use 5th-order basis polynomials 3 per element, and the vertical resolutions are the same as the horizontal resolutions.4 5

S
.-J.Choi et al.: Verification of a non-hydrostatic dynamical core tical derivatives.The time discretization relied on the timesplit, third-order Runge-Kutta technique.

Table 1 .
Summary of the resolutions and time-step sizes used for the tests.