Validation of a Global Nonhydrostatic Atmospheric Dynamical Core Using Discontinuous Galerkin Method

1. Introduction

In future high-resolution atmospheric simulations with grid spacing of O(10 m), such as global large eddy simulations (LES), we consider that discretization accuracy for atmospheric dynamical cores is one of issues. It is possible that the numerical errors dominate over the effect of sub-grid scale parametrizations. In the framework of conventional grid-point methods, Kawai and Tomita (2021, MWR) discussed this issue and indicated the eighth-order accuracy is required for advection schemes to precisely conduct LES of planetary boundary layer turbulence. On the other hands, the state-of-the-art global nonhydrostatic atmospheric models typically adopt low-order grid-point methods such as totally second-order accuracy (e.g., Tomita and Satoh, 2004). To achieve significant high-order accuracy, Kawai & Tomita (2023, MWR) focuses on discontinuous Galerkin method (DGM) because of straightforward strategy for high-order discretization and computational compactness compared to the conventional grid-point methods. In addition, in terms of the dynamics-physics coupling, DGM with many degrees of freedom per element may give us novel findings for utilizing accurate dynamical fields to calculate the physical processes.

To investigate the applicability of DGM to global atmospheric simulations, we are now developing a global dynamical core based on DGM and validating it by various idealized numerical experiments. In this presentation, we will show the formulation of our dynamical core and preliminary results of validation experiments.

2. Formulation of global dynamical core

The governing equation is the fully compressible nonhydrostatic equations of dry atmosphere. The spherical geometry and topography are treated by an equiangular cubed-sphere projection and a basic terrain-following coordinate, respectively. The spatial discretization is based on nodal DGM (e.g., Hesthaven and Warburton, 2007). The computational domain is divided by hexahedral elements. Each element has (p+1)³ degree of freedoms where p is the order of expansion polynomial. The fluxes at element boundaries are evaluated using Rusanov flux. As the temporal scheme, we use a horizontally explicit and vertically implicit (HEVI) scheme with the third-order accuracy.

3. Validation experiments

To validate our developing dynamical core for global atmospheric model, we perform a series of numerical experiments of linear advection, internal gravity wave, baroclinic instability, mountain wave, and Held-Suarez tests. We confirmed that the simulation results for all test cases qualitatively reproduce the corresponding reference solutions obtained from previous studies. Further, we investigate the behavior of numerical convergence for various p to evaluate the impact of high-order DGM on atmospheric fields. This presentation focuses on two test cases of gravity wave and baroclinic instability.

The experimental setup for gravity wave test is based on Baldauf and Brdar (2013). The horizontal and vertical effective grid spacing are changed as (Δ_h, Δ_z)〜(313 km, 417 m), (156 km, 208 m), (78 km, 104 m) using p=1, 3, 7. We investigate the self-convergence of numerical solutions using a reference solutions obtained by an experiment with (Δ_h, Δ_z)〜(39 km, 52 m) with p=7. Figure (a) shows a dependence of L2 error norms of winds on the spatial resolution. This indicates the convergence rate for the zonal wind follows almost p+1-order accuracy. On the other hand, the convergence rate for the vertical wind related to fast modes is smaller than the p+1 order accuracy for p>4. The reason is related to the third-order temporal accuracy of HEVI scheme because the convergence rate approaches to p+1 order spatial accuracy with the decrease of timestep (see the red and magenta lines in the figure).

The experimental setup for baroclinic instability test is based on Jablonowski and Williamson (2006) (hereinafter, referred to as JW06). The horizontal effective grid spacing is changed as Δ_h〜250 km, 125 km, 63 km using p=7. Following JW06, we fix the number of vertical degrees of freedom, which is set to 24 in this study. The reference solution is obtained by an experiment with Δ_h〜31 km with p=7. Figure (b) shows a temporal evolution of L2 error norms of surface pressure and the spatial resolution dependence. The dashed lines represent the results from Mcore (Ullrich and Jablonowski, 2012), in which a fourth-order finite volume method is used as the horizontal discretization. The results of both models become comparable when the baroclinic wave starts to develop significantly. The decrease in the error slows down when Δ_h<250 km even during the first five days. One of the reasons may be related to the experimental setup of fixed vertical resolution. However, the L2 error norms for Δ_h<250 km are within a range of uncertainty suggested by JW06 (represented by shade in the figure), and we consider that the numerical solutions from our model based on DGM are reasonable.