Vertical discretization with finite elements for a global hydrostatic model on the cubed sphere

Abstract

A formulation of Galerkin finite element with basis-spline functions on a hybrid sigma-pressure coordinate is presented to discretize the vertical terms of global Eulerian hydrostatic equations employed in a numerical weather prediction system, which is horizontally discretized with high-order spectral elements on a cubed sphere grid. This replaces the vertical discretization of conventional central finite difference that is first-order accurate in non-uniform grids and causes numerical instability in advection-dominant flows. Therefore, a model remains in the framework of Galerkin finite elements for both the horizontal and vertical spatial terms. The basis-spline functions, obtained from the de-Boor algorithm, are employed to derive both the vertical derivative and integral operators, since Eulerian advection terms are involved. These operators are used to discretize the vertical terms of the prognostic and diagnostic equations. To verify the vertical discretization schemes and compare their performance, various two- and three-dimensional idealized cases and a hindcast case with full physics are performed in terms of accuracy and stability. It was shown that the vertical finite element with the cubic basis-spline function is more accurate and stable than that of the vertical finite difference, as indicated by faster residual convergence, fewer statistical errors, and reduction in computational mode. This leads to the general conclusion that the overall performance of a global hydrostatic model might be significantly improved with the vertical finite element.

Introduction

Developing a numerical scheme for the discretization of governing equations is one of the key issues in constructing a numerical model of the atmosphere. Since the spatial terms of the governing equations employed in most numerical weather prediction (NWP) models [1], [2] are completely decoupled in the horizontal and vertical directions, both terms might be discretized with different numerical schemes. Despite various spatial discretization schemes, such as finite volumes, finite elements, and their variants, have been adopted for the horizontal terms of the governing equations, the finite difference scheme of Simmons and Burridge [3] is still commonly used for the discretization of vertical terms. This conventional central finite difference is first-order accurate for non-uniform grids [4], which is too low accurate, compared with the horizontal discretization schemes used in recently developed dynamical cores [5], [6], [7]. Moreover, all central difference schemes for solving advection equations suffer from dispersion errors [8], so the vertical finite difference (VFD) scheme with the Lorenz staggering grids has been known to increase the computational mode in distributions of temperature through incorrect hydrostatic balance [9], [10]. This results in an increase in numerical instability.

To overcome the limits of the VFD in terms of accuracy and stability, vertical discretization with finite elements has been attempted in several atmospheric models [11], [12], [13], [4]. Staniforth et al. [11] first tried to develop a full Galerkin model, using spherical harmonics for the horizontal discretization and finite elements for the vertical discretization in the sigma coordinate. Using the vertical finite element (VFE) with a linear basis function for the vertical terms of baroclinic primitive equations, they obtained promising results with short-range forecasts. A comparison of the accuracy between the VFE and VFD was performed with the linearized model of baroclinic primitive equations in the sigma vertical coordinate by Béland et al. [12]. They reported that the VFD generated the computational mode near the top of the model, resulting in a convergence problem. Studies of the vertical discretization based on the Galerkin finite element were extended to higher-order basis functions and other vertical coordinate systems. The VFE with second-order quadratic basis functions in the sigma coordinate showed a systematic improvement over the 10-day control runs that were vertically discretized with the finite difference [13]. Their works on further developing the VFE scheme continued with a hybrid vertical coordinate system [14]. For the global hydrostatic spectral model with semi-Lagrangian advection, Untch et al. [4] employed both linear and cubic basis functions in the hybrid vertical coordinate system. Since semi-Lagrangian vertical advection was used in the model, only the vertical integral terms of the continuity and hydrostatic equations were discretized with the Galerkin finite element. They found that the VFE with the cubic basis function reduced vertically propagating computational modes arising from the VFD scheme with the Lorenz staggering grid by producing the phase speeds of gravity waves more accurately. This resulted in improved model performance, especially in the stratosphere. Recently, attempts to discretize with the VFE scheme in non-hydrostatic models have been achieved in limited-area models [15], [16]. Most previous studies of vertical discretization with the finite element have been involved in global and regional atmospheric models using horizontal spectral methods with spherical harmonics and Fourier series, respectively. For the conservation properties of the VFE scheme, the energy conservation with barotropic equations was proved by Stepperler [13], which was given independent of the order of the basis functions.

In this work, a VFE scheme applied to the global hydrostatic version of Korea Institute of Atmospheric Prediction Systems Integrated Model (KIM) is developed to replace the VFD scheme because of its poor accuracy and stability. The present version of KIM is based on the high-order spectral element method on a cubed sphere grid and the finite difference method of Simmons and Burridge [3] in Lorenz grids of hybrid sigma-pressure (or η) vertical coordinates, of which the hydrostatic model was originated from High-Order Method Modeling Environment (HOMME) [5]. However, it is difficult to extend the same high-order spectral element method for vertical discretization owing to coupling with the physics at the points interpolated within each element. Therefore, the main purpose of this development is not only to keep using the high-order discretization scheme for the vertical terms, but also to easily couple them with existing physics packages, which have been developed for a Lorentz staggering grid system. This results in the model remaining in the same framework of the Galerkin method as the horizontal discretization. The global hydrostatic equations used in this study are based on a fully Eulerian approach, while HOMME uses horizontally Eulerian and vertically Lagrangian discretizations [17]. Unlike the work of Untch et al. [4], the fully Eulerian hydrostatic equations to be discretized with the VFE scheme are required to derive both the vertical derivative and integral operators. The basis-spline function employed when computing these operators as a weight function of the finite elements is obtained from the de-Boor recursive algorithm [18], in which the knots are a set of hybrid η vertical grids.

The remainder of this paper is organized as follows. In Section 2, the three-dimensional global hydrostatic equations to be vertically discretized with the VFE are presented on the hybrid η coordinate. Section 3 describes how to derive the basis functions and vertical operators, and discretize the vertical terms of the governing equations. In Section 4, the accuracy and stability of the vertical discretization schemes are demonstrated through two- and three-dimensional test cases using several assessment tools such as the bias score (BIAS), root-mean-square-error (RMSE), normalized error norms, and qualitative comparison. These test cases include three-dimensional idealized and short-range hindcasting cases. Finally, this work is summarized and concluded with future work suggestions in Section 5.