Verification of a non-hydrostatic dynamical core using horizontally spectral element vertically finite difference method: 2D Aspects

Abstract : The non-hydrostatic (NH) compressible Euler equations of dry atmosphere are solved in a simplified two dimensional (2D) slice (X-Z) framework employing a spectral element method (SEM) for the horizontal discretization and a finite difference method (FDM) for the vertical discretization. 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 terms and quadrature. The Euler equations used here are in a flux form based on the hydrostatic pressure vertical coordinate, which are the same as those used in the Weather Research and Forecasting (WRF) model, but a hybrid sigma-pressure vertical coordinate is implemented in this model. We verified the model by conducting widely used standard benchmark tests: the inertia-gravity wave, rising thermal bubble, density current wave, and linear hydrostatic mountain wave. The numerical results demonstrate that the horizontally spectral element vertically finite difference model is accurate and robust. By using the 2D slice model, we effectively show that the combined spatial discretization method of the spectral element and finite difference method in the horizontal and vertical directions, respectively, offers a viable method for the development of a NH dynamical core. The present core provides a practical framework for further development of three-dimensional (3D) non-hydrostatic compressible atmospheric models.


Introduction
There is a 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). The spectral element method (SEM) has been known as one of the most promising methods with high efficiency and accuracy (Taylor et al., 1997;Giraldo, 2001;Thomas and Loft, 2002). SEM is local in nature because 25 of having a large on-processor operation count (Kelly and Graldo, 2012 achieves this high level of scalability by decomposing the physical domain into smaller pieces with a small communication stencil. Also SEM has been shown to be very attractive in achieving high-order accuracy and geometrical flexibility on the sphere (Taylor et al., 1997;Giraldo, 2001;Giraldo et al., 2004). To date, the SEM has been successfully implemented in atmospheric modeling such 5 as in the Community Atmosphere Model -Spectral Element dynamical core (CAM-SE) (Thomas 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 grid meshes such as a cubed-sphere tiled with quadrilateral elements using SEM in the horizontal discretization and the finite differ-10 ence method (FDM) in the vertical. The robustness of the SEM has been illustrated through three-dimensional dry dynamical test cases (Thomas and Loft, 2005;Giraldo and Rosmond, 2004;Giraldo, 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 the combined horizontally 15 SEM and vertically FDM. Since testing a 3-D NH model requires a huge amount of computing resources, studying the feasibility of our approach in 2-D is an attractive alternatively to the development of a fully 3-D model. This is the case because a 2-D slice model effectively can test the practical issues resulting from the vertical discretization and time integration, prior to the construction of a full 3-D model. Although we could 20 also discretize the vertical direction with SEM (as is proposed in Kelly and Giraldo, 2012;Giraldo et al., 2013), we choose to use a conservative flux-form finite-difference method for discretization in the vertical direction, which provides an easy way for coupling the dynamics and existing physics packages.
We have developed a dry 2-D NH compressible Euler model based on SEM along 25 the x direction and FDM along the z direction for this purpose. Hereafter, this is simply referred to as the 2DNH model. We adopt the governing equation formulation proposed by Skamarock and Klemp (2008)  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | on the hydrostatic pressure vertical coordinate. In SK08, the terrain-following sigmapressure coordinate is used, but here we employ a hybrid sigma-pressure vertical coordinate. Park et al. (2013) (hereafter, PK13) provides 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 5 the sigma-pressure coordinate at the level of the actual coding implementation. Also, we built the 2DNH model using a time-split third-order Runge-Kutta (RK3) for the time discretization, which has been shown to work effectively in the WRF model. We keep 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 2DNH model. In this paper, we show the feasibility of the 2DNH model by conducting conventional benchmark test cases as well as focusing on the description of the numerical scheme for the spatial discretization. We verify the 2DNH by analyzing four test cases: the 15 inertia-gravity wave, rising thermal bubble, density current wave, and linear hydrostatic mountain wave.
The organization of this paper is as follows. In the next section we describe the governing equations with a definition of the prognostic and diagnostic variables used in our model, in which we present essential changes from SK08. In Sect. 3 we explain tempo-20 ral and spatial discretization including the spectral element formulation. In Sect. 4, we present the results of the 2DNH model using all four test cases. Finally we summarize the paper and propose future directions in Sect. 5.

Governing equations
We adopt the formulation of the governing equation set of SK08.
where p d is the hydrostatic pressure of dry air, B(η) is the relative weighting of the terrain-following coordinate vs. the normalized pressure coordinate, p s , p t , and p 0 are 5 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,η ≡ ∂η ∂t 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) 15 The flux-form Euler equations for dry atmosphere 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 the Coriolis and curvature which we ignore for simplicity. In

5
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 = p s (x, y) + p s where the variables denoted by overbars are 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 as where ∂p s ∂t is obtained by integrating vertically Eq. (6) from the surface (η = 1) to the top (η = 0) using a no-flux boundary condition such as Ω| η=0 or 1 = 0 and the specification of the vertical coordinate such as B(η = 1) = 1 and B(η = 0) = 0 as The above equation allows 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 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 5 equations that are used in our model.

Horizontal direction
For a given η level, we discretize the horizontal operators using the SEM. Therefore in 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 where x k represents N +1 grid points that correspond to the Gauss-Lobatto-Legendre

15
(GLL) points and ψ k (x) are the Nth-order Lagrange polynomials based on the GLL points. It is noteworthy that the ψ k have the cardinal property, i.e., they can be represented as Kronecker delta functions where ψ k are zero at all nodal points except x k (but are allowed to oscillate between nodal points). The GLL points ξ k in a reference coordinate system ξ ∈ [−1, +1] and the associated 20 quadrature weights ω(ξ k ), where P N (ξ) are the Nth-order Legendre polynomials, J = ∂x ∂ξ is the transformation Jacobian, and Ω e represents the non-overlapping elements. 5 We now introduce the polynomial expansions into our governing equations in the form of multiply by the basis function as a test function, and integrate to yield a system of 10 ordinary differential equations as such where k = 1, 2, · · · , N + 1, M e nk is the element based mass matrix given as

15
The right-hand sides of Eqs. (17) and (18) is evaluated using 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 nk , 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 5 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 , Θ, µ, α, p are at layer midpoints denoted by k = 1, 2, . . . , K where K is the total number of layers, and the variables W , Ω, ϕ live 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. 15 Here, we evaluate the vertical advection terms ( ∂η , 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 20 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. Likewise the vertical discretization quadrature rules for the calculations of Eqs. (9) and (10)

Temporal discretization
For integrating the equations, we use the time-split RK3 integration technique following the strategy of SK08, in which low-frequency modes due to advective forcings are explicitly advanced using a large time step of the RK3 scheme, but high-frequency modes are integrated over smaller time steps using an explicit forward-backward time integra-5 tion scheme for the horizontally propagating acoustic/gravity waves and a fully implicit scheme for vertically propagating acoustic waves and buoyancy oscillations (Klemp et al., 2007). This technique has been shown to work effectively within numerous nonhydrostatic models including the WRF model ( It is noted that in the procedure of the time-split RK3 integration, the difference between the approach used in this paper and 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 the forward-backward differencing on the small 15 time steps, then µ d can be computed directly from the specification of the vertical coordinate in Eq. (9) and Ω is obtained from the vertical integration.

Test cases
We validate the 2DNH model on four test cases of the linear hydrostatic mountain wave, density current, inertia-gravity wave, and rising thermal bubble experiments. All cases 20 but the mountain wave experiment do not have analytic solutions. Therefore, for the mountain wave experiment, numerical results of the 2DNH model are compared to analytic solutions (Durran and Klemp, 1983), and for the other experiments, we compare our results to the results of other published papers.
It should be mentioned that the horizontal SEM formulation is able to utilize arbitrary 25 order polynomials per element to represent the discrete spatial operators, but in this paper all the results presented use either 5th or 8th 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 5 in the domain and N is the number of the interval associated with irregularly spaced GLL quadrature points which is equivalent to the order of the basis polynomials. The average vertical grid spacing is defined in the same way of Eq. (24). Below, we use this convention to define the grid resolution.

Linear hydrostatic mountain wave test 10
In order to verify the 2DNH's feasibility to treat surface elevations associated with the vertical terrain-following coordinate, we simulate the linear hydrostatic mountain wave test introduced by Durran and Klemp (1983) (hereafter, DK83) in which the analytic steady-state solution is provided by using a single-peaked mountain with uniform zonal wind. To compare our results with the analytic and numerical solution shown in DK83, 15 the 2DNH is initialized using the same initial conditions and mountain profile in DK83 and we analyze our results using the same metrics of DK83. The mountain profile is given by 20 where the half-length of the mountain a m is 10 km, the height h m is 1 m, and the prescribed center of the profile is 0 km. The initial temperature is T 0 = 250 K for an isothermal atmosphere with the uniform zonal wind u = 20 m s −1 . In the isothermal case, the  our result in the horizontal velocity is in good agreement with DK83 and GR08. Figure 3 shows the normalized momentum flux values at various times to check vertical transport of horizontal momentum. It is observed that the flux is developing well and the simulations have reached steady-state after ut a = 60. It is noted that the mean momentum flux at that time is 97 % of its analytic value. It agrees well with DK83 20 as well as GR08; it is important to point out that the 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 the terrain-following vertical coordinate is well suited for the combination of the horizontal SEM and vertical FDM for spatial discretization even though we consider a small mountain.

2-D density current test
In order to verify the 2DNH's feasibility to control oscillations with numerical viscosity and evaluate numerical schemes in the 2DNH, we conduct 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 5 density current wave starts to spread symmetrically in the horizontal direction forming Kelvin-Helmholtz rotors. Following Straka et al. (1993) we employ a dynamic viscosity of ν = 75 m 2 s −1 to obtain converged numerical solutions.
For an initial cold bubble, the potential temperature perturbation is given as   Figure 6 shows the profiles of the potential temperature perturbation at the height of 1200 m. The results from the highest grid resolution of the simulations using 5th and 8th order basis polynomials are indistinguishable and well converged (Fig. 6a). Three minima corresponding to the three rotors agree well with other published solutions. 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 agree well with each other (Table 1), of which the numbers are comparable to those of GR08. In the numerical results from the different grid resolutions simulated by using 5th order basis polynomials, the potential temperature profile at the coarsest resolution of 400 m grid shows significant fluctuations (Fig. 6b). That of 8th order poly-15 nomials, however, tends to be relieved from the deviation from the converged solution (Fig. 6c). The above results suggest that the numerical solution can be converged more rapidly by using higher order of basis polynomial. Furthermore, the results in this paper show that an adequate convergence can be reached at grid resolutions finer than 200 m.

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 perturbation diverges to the left and right symmetrically in a channel with periodic lateral boundary conditions. The inertia-gravity wave test introduced by Skamarock and Klemp (1994)  Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | frequency of N = 0.01 s −1 with surface potential temperature of θ 0 = 300 K and a uniform zonal wind u = 20 m s −1 . In order to trigger the wave, the initial potential temperature perturbation θ is overlaid to the above initial state and is given as  It is noted that all experiments give almost the same values for potential temperature perturbation where these values in the range θ ∈ [−1.52 × 10 −3 , 2.83 × 10 −3 ] are comparable to other studies (e.g., GR08 and Li et al., 2013). GR08 give the ranges of θ ∈ [−1.51 × 10 −3 , 2.78 × 10 −3 ] from the model based on the spectral element and discontinuous Galerkin method. Also Li et al. (2013) show θ ∈ [−1.53 × 10 −3 , 2.80 × 10 −3 ] 5 using the high-order conservative finite volume model which are similar to our results.

Rising thermal bubble test
We also conduct 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 envi-10 ronment for an atmosphere at rest. The air that is warmer than the 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 are evaluated qualitatively.

15
The initial conditions we use follow those of GR08 in which the domain for the case is defined as (x, z) ∈ [0, 1] 2 km 2 . We consider no-flux boundary conditions for all four boundaries. The domain is initialized for a neutral atmosphere at rest with θ 0 = 300 K in hydrostatic balance. A 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. The model was run for a time of 700 s. It should be noted that an explicit second-order diffusion on coordinate surfaces is used with a viscosity coefficient of ν = 1 m 2 s −1 for all simulations 25 of this test. The numerical diffusion is applied for momentum and potential temperature 3732 Introduction

Tables Figures
Back Close

Full Screen / Esc
Printer-friendly Version

Interactive Discussion
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | along the horizontal and vertical directions so that it eliminates the erroneous oscillations at the small scale -while this amount of diffusion might seem excessive, it has been chosen because it allows the model to remain stable even after the bubble hits the top boundary. Figure 9 shows the potential temperature perturbation, horizontal wind, and vertical 5 wind fields for the simulations of two resolutions of 20 and 5 m horizontal and vertical grid spacings (∆x and ∆z), respectively, employing 5th order basis polynomials. In both simulations, the fine structures in the numerical solutions are well depicted with a perfectly symmetric distribution at the midpoint and sharp discontinuities of the fields along boundary lines of the bubble. At lower resolution, however, degradations in the 10 solution are visible in the potential temperature perturbation and vertical wind which are illustrated by fluctuations in the values as well as the concaving contours at the top of the bubble. It is noted that while 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 simulations here of 2DNH using the horizontally 15 SEM and vertically FDM 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. 10. Simulations were run with various resolutions of 5, 10, and 20 m, where the resolutions given are defined for both the horizontal and vertical directions. The results of 10 and 5 m resolutions are almost identical to each other. 20 The result of the lowest resolution of 20 m, however, shows a somewhat unresolved solution, in which the maximum value is underestimated and the phase shift is depicted. The time series for maximum potential temperature perturbation and maximum vertical velocity are shown in Fig. 11. In all simulations, the maximum vertical velocity increases as the maximum theta perturbation decreases. This shows that the thermal energy 25 of the theta perturbation leads to the acceleration of the vertical velocity. This result agrees well with the study of Ahmad and Lindeman (2007). 7,2014 Verification of a non-hydrostatic dynamical core

Summary and conclusions
The non-hydrostatic compressible Euler equations for a dry atmosphere are solved in a simplified 2-D slice (X-Z) framework by using the spectral element discretization (SEM) in the horizontal and the third-order finite difference scheme for the vertical discretization. The form of the Euler equations used here are the same as those used 5 in the Weather Research and Forecasting (WRF) model. We employ a hybrid sigmapressure 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 are separated into their horizontal and vertical components. In the horizontal components, the operators are discretized 10 using the 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 are discretized using the third-order upwind biased finite difference scheme for the vertical fluxes and 15 centered differences for the vertical derivatives. The time discretization relies on the time-split third-order Runge-Kutta technique.
We have presented results from idealized standard benchmark tests for large-scale flows (e.g., linear hydrostatic mountain wave) and for nonhydrostatic-scale flows (e.g., inertia-gravity wave, rising thermal bubble, and density currents). The numerical results 20 show that the present dynamical core is able to produce solutions of good quality comparable to other published solutions. These tests effectively reveal that the combined spatial discretization method of the spectral element and finite difference method in the horizontal and vertical directions, respectively, offers a viable method for the development of a NH dynamical core. Further research will be continued to couple the present Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Klemp, J. B., Skamarock, W. C., and Dudhia, J.: Conservative split-explicit time integration methods for the compressible nonhydrostatic equations, Mon. Weather Rev., 135, 2897-2913, 2007: Rotated versions of the Jablonowski steady-state and baroclinic wave test cases: a dynamical core intercomparison, J. Adv.