Development of a high-order global dynamical core using the discontinuous Galerkin method for an atmospheric large-eddy simulation (LES) and proposal of test cases: SCALE-DG v0.8.0

Focusing on the future global atmospheric simulations with a grid spacing of O(10–100 m), we developed a global nonhydrostatic atmospheric dynamical core with high-order accuracy by applying a discontinuous Galerkin method (DGM) for horizontal and vertical discretization. Furthermore, considering a global large-eddy simulation (LES), a Smagorinsky–Lilly turbulence model was introduced to the proposed global dynamical core in the DGM framework. By conducting several tests with various polynomial (p) orders, the impact of the high-order DGM on the accuracy of the numerical simulations of atmospheric flows was investigated. To show high-order numerical convergence, a few modifications were made in the experimental setup of existing test cases. In addition, we proposed an idealized test case to verify global-LES models, which is a global extension of an idealized planetary boundary layer (PBL) turbulence experiment performed in our previous studies. The error norms from the deterministic test cases, such as the linear-advection and gravity-wave tests, showed an optimal convergence rate achieved by an approximately p+1-order spatial accuracy when the temporal and round-off errors were sufficiently small. In the climatic test cases, such as the Held–Suarez test, the kinetic energy spectra indicated the advantage of effective resolution when large polynomial orders were used. In the LES experiment, the global model provided a reasonable vertical structure of the PBL and energy spectra because the results under shallow-atmosphere approximation reproduced those obtained in the plane computational domain well.

1 Introduction

Recently developed supercomputers have enabled us to conduct high-resolution global atmospheric simulations using a sub-kilometer horizontal grid spacing. For example, Miyamoto et al. (2013) conducted a global simulation at a horizontal grid spacing of 870 m and discussed the numerical convergence of statistical properties of deep moist convections. In the near future, this continuous development of computer technology is expected to enable us to perform global simulations using O(10–100 m) grid spacing (Satoh et al., 2019), which begin to explicitly represent turbulence in the inertial sub-range. Then, large-eddy simulation (LES) is a promising strategy, since in LES, the turbulence in a spatial scale larger than a spatial filter is explicitly calculated, whereas the effect of turbulence in a smaller spatial scale is parameterized using eddy viscosity and diffusion terms. By explicitly representing the large-scale eddies in boundary layers and the low-level clouds such as shallow cumuli, we expect to reduce a source of uncertainty associated with the parameterizations and improve representation of the global radiation budget in a realistic representation of Earth's atmosphere.

Considering future high-resolution atmospheric simulations such as global LES, we discussed the problem of low numerical accuracy of conventional atmospheric dynamical cores in Kawai and Tomita (2021, 2023) (hereafter referred to as KT2021 and KT2023, respectively). To perform LES precisely, we must ensure that the discretization errors do not dominate over sub-grid-scale (SGS) terms of turbulent models. Otherwise, the SGS terms might lose their physical significance. KT2021 investigated the order of accuracy necessary for advection schemes in the framework of conventional grid-point methods. In particular, the study derived two ratios associated with numerical diffusion and numerical dispersion: the ratio of decay time with the SGS terms to that of the numerical diffusion error terms and the ratio of phase speed due to the error in advection terms to that of the SGS terms. Moreover, we pointed out that the advection scheme requires at least seventh- or eighth-order accuracy to ensure that both ratios are less than 10⁻¹ at wavelengths longer than eight grid lengths for grid spacing simulations of O(10) m. However, in conventional grid-point methods, the required stencil grows larger as the order of accuracy increases. This can degrade the computational performance of the high-order methods in recent massively parallel computers. Thus, we focused on applicability of the discontinuous Galerkin method (DGM). At element boundaries, the representation of the flow field is allowed to be discontinuous, and the flux shared by two neighboring elements is calculated using approximated Riemann solvers. These computational features provide a straightforward strategy to achieve high-order discretization and computational compactness. KT2023 extended the discussion presented in KT2021 to the DGM framework and investigated a polynomial order necessary for precisely conducting LES. It indicated that the polynomial order needs to be higher than or equal to 4 when the upwind numerical fluxes and sufficiently scale-selective modal filters are used to ensure numerical stability.

The basis of the state-of-the-art global nonhydrostatic atmosphere dynamical cores was mainly developed during the 2000s–2010s. In these dynamical cores, low-order grid-point methods are often adopted. For example, the Nonhydrostatic ICosahedral Atmospheric Model (NICAM; Tomita and Satoh, 2004; Satoh et al., 2014), the Model for Prediction Across Scales (MPAS; Skamarock et al., 2012), and the Icosahedral Nonhydrostatic (ICON) model (Zängl et al., 2015) are based on either a totally first- or totally second-order scheme. The discretization accuracy has not always been a primary factor in the performance of atmospheric models because physical processes have various uncertain parameters. In situations where the grid spacing is coarser than the gray zone of turbulence, the totally second-order scheme may be appropriate in terms of computational cost and numerical robustness. However, as described above, it is crucial to increase numerical accuracy to precisely conduct the atmospheric LES. Furthermore, even in spatial resolutions coarser than that required by LES, it is undesirable that the shortest wavelength fully resolved by discretization methods, the so-called effective resolution, is significantly different from the grid spacing in terms of the physics–dynamics coupling. The low-order spatial scheme typically leads to significant discretization errors at wavelengths shorter than eight grid lengths (Kent et al., 2014). To decrease the gap between the effective resolution and grid scale, it is sensible to use high-order discretization methods in addition to designing better numerical filters for controlling the effective resolution.

Constructing high-order grid-point methods tends to be more complex for spherical geometries than in plane domains with structured grids. To achieve the high-order discretization accuracy horizontally, a conventional straightforward approach is the spectral transformation method based on the spherical-harmonics expansion. This approach provides desirable accuracy for numerical solutions in the wavelength range up to the truncated wave number while avoiding the problem of a restrictive time step near the poles due to the convergence of meridians. However, high-resolution global simulations can suffer from the large cost of data communication between all computational nodes in massively parallel supercomputers. On the other hand, some researchers have successfully developed global atmospheric dynamical cores based on high-order grid-point and element-based methods. The essence of the numerical methods can be found in horizontal discretization of the global shallow-water equations, for example, Ullrich et al. (2010) for a high-order finite-volume method (FVM) and Nair et al. (2005a) and Ullrich (2014) for high-order element-based methods. Ullrich and Jablonowski (2012b) proposed a global nonhydrostatic dynamical core (MCore) based on an FVM with a fourth-order horizontal reconstruction strategy. The Tempest model (Guerra and Ullrich, 2016) uses a high-order spectral-element method (SEM) horizontally. In the Non-hydrostatic Unified Model of the Atmosphere (NUMA; Kelly and Giraldo, 2012; Giraldo et al., 2013), which is applicable for both limited-area and global atmospheric simulations, the continuous and discontinuous Galerkin methods are adopted for the spatial discretization. The numerical-method prototype used in NUMA is utilized and extended to a global spectral-element dynamical core in the Navy Environmental Prediction System Utilizing a Nonhydrostatic Engine (NEPTUNE) for both horizontal and vertical discretization (e.g., Zaron et al., 2022). SEM is also used for the nonhydrostatic High-Order Method Modeling Environment (HOMME-NH; Dennis et al., 2005, 2012; Taylor et al., 2020) included in the Energy Exascale Earth System Model (E3SM) and for the nonhydrostatic dynamical core in the Korean Integrated Model (KIM) system (Hong et al., 2018). ClimateMachine uses a nodal discontinuous Galerkin method for both horizontal and vertical discretization. The corresponding regional dynamical core is described in Sridhar et al. (2022). In the case of classical high-order element-based methods, it is cumbersome to control the numerical instability caused by the aliasing errors with nonlinear terms (Winters et al., 2018). To overcome this problem, a split-form nodal DGM (e.g., Gassner et al., 2016) is a theoretical and computationally efficient approach. A similar method was successfully applied to a global shallow-water model in Ricardo et al. (2024) and to a global nonhydrostatic dynamical core for horizontal and vertical discretization in Souza et al. (2023). While conventional dynamical cores adopt a vertical discretization based on a low-order finite-difference method or FVM, some studies investigated the potential for the use of high-order vertical discretization (e.g., Guerra and Ullrich, 2016; Yi and Giraldo, 2020; Ishioka et al., 2022). For example, Guerra and Ullrich (2016) introduced an arbitrary-order vertical discretization using a staggered nodal finite-element method (FEM). They reported that high-order vertical discretization improves the representation of vertical dynamics at a relatively low vertical resolution.

By building on progress from the previous studies showing the applicability of the element-based methods to atmospheric-flow simulations, the current study attempted to develop a high-order global dynamical core using a nodal DGM both horizontally and vertically for future global atmospheric simulations with O(10–100 m) grid spacing. For a quasi-uniform spherical grid, a cubed-sphere projection was adopted. Moreover, a terrain-following coordinate was used to treat the topography. Although the numerical methods used in our dynamical core are similar to those used in previous studies that developed global DG dynamical cores such as NUMA and ClimateMachine, we consider the following points to be the unique contributions of the current study:

1. Considering global LES, we formulated SGS eddy viscous and diffusion terms with a Smagorinsky–Lilly type turbulent model in the DGM framework on the cubed-sphere coordinates. A discretization strategy for the scalar Laplacian operator on the cubed-sphere coordinates with the DGM is reported in previous studies (e.g., Nair, 2009). However, they did not treat the vector Laplacian operator for the vector quantities (for example, momentum). This might be because the rigorous form is so complex that it may not be worth the computational cost required for numerical stabilization. On the other hand, Ullrich (2014) presented a discretization strategy for the vector Laplacian operator with the continuous and discontinuous Galerkin methods. This approach can distinguish the divergence damping and vorticity damping with constant viscous coefficients. Guba et al. (2014) proposed a strategy of hyperviscosity with variable viscous coefficients in SEM where the vector Laplacian operator is applied to the Cartesian component of the vector fields. However, we did not directly use these approaches to the vector Laplacian when introducing the turbulent model. This is because we need to treat eddy viscous and diffusion coefficients as being dependent on local wind shear and stratification. In addition, we consider that the vector Laplacian operator applied to the vector component on the cubed-sphere coordinates can be convenient for the straightforward distinction between horizontal and vertical directions. Using tensor analysis, we systematically derived the eddy viscosity and diffusion terms. Then, we represented the corresponding semi-discretized equations using the DGM. Subsequently, a quantitative check was performed by conducting an idealized LES experiment of planetary boundary layer (PBL) turbulence. In particular, we extended a numerical experiment of idealized planetary boundary turbulence used in regional plane models (KT2021 and KT2023) to spherical geometry by slightly changing the initial condition.

2. We modified experimental settings of idealized test cases to demonstrate the numerical convergence with high-order dynamical cores. When using totally second-order dynamical cores, relatively large discretization errors may occur, which can overshadow the problems of ill-posed experimental settings. Fast numerical convergence achieved using high-order schemes is expected to enable detection of such problems in standard tests. Even when research interests do not include the dynamics, we consider an evaluation framework using a high-order dynamical core to be useful. For example, when new physical models are included, the physical performance can be directly evaluated by reducing numerical errors with the dynamical processes.

3. We attempted quantitative evaluations in a series of test cases for the global dynamical cores to investigate the impact of a high-order DGM on the numerical accuracy of atmospheric-flow simulations. Although the numerical-convergence characteristics of the DGM was closely investigated in the case of regional dynamical cores (e.g., Giraldo and Restelli, 2008; Brdar et al., 2013; Blaise et al., 2016), few studies have been conducted to demonstrate it for global nonhydrostatic dynamical cores.

The rest of this paper is organized as follows. In Sect. 2, the governing equations using the general curvilinear coordinates are formulated. Then, we introduce a cubed-sphere projection and a general vertical coordinate. We represent eddy viscous and diffusion terms associated with the turbulent model in the general curvilinear coordinates. Furthermore, we explain the spatial and temporal discretization adopted for the governing equations. In Sect. 3, we verify the proposed dynamical core through several idealized numerical experiments. Finally, the findings of this study and our future plans are summarized.

2 Model description

2.1 Governing equations

As governing equations for dry atmospheric flows, we used the three-dimensional, fully compressible nonhydrostatic equations based on the flux form (e.g., Ullrich and Jablonowski, 2012b). Following Li et al. (2020), a non-orthogonal curvilinear horizontal coordinate (ξ,η) was introduced. Subsequently, a general vertical coordinate ζ was introduced. For the horizontal coordinate transformation, the contravariant form of the metric tensor is represented by $G_{\mathrm{h}}^{ij}$ for $i,j=\mathrm{1},\mathrm{2}$. A three-dimensional metric tensor with the horizontal coordinate transformation is defined as

$$\begin{array}{}\text{(1)}& {\displaystyle }{G}^{ij}=\left(\begin{array}{ccc}{G}_{\mathrm{h}}^{\mathrm{11}}& G_{\mathrm{h}}^{\mathrm{12}}& \mathrm{0}\\ G_{\mathrm{h}}^{\mathrm{21}}& G_{\mathrm{h}}^{\mathrm{22}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right).\end{array}$$

The horizontal Jacobian is defined as $\sqrt{G_{\mathrm{h}}}=|G_{\mathrm{h}}^{ij}{|}^{-\frac{\mathrm{1}}{\mathrm{2}}}$. For the vertical coordinate transformation, the metric tensor is defined as $G_{\mathrm{v}}^{\mathrm{13}}=\partial \mathit{\zeta }/\partial \mathit{\xi },G_{\mathrm{v}}^{\mathrm{23}}=\partial \mathit{\zeta }/\partial \mathit{\eta }$ and the vertical Jacobian is defined as $\sqrt{G_{\mathrm{v}}}=\partial z/\partial \mathit{\zeta }$. The vertical velocity in the transformed vertical coordinate can be written using contravariant components of wind vector $(u^{\mathit{\xi }},u^{\mathit{\eta }},u^{\mathit{\zeta }})$ as

$$\begin{array}{}\text{(2)}& {\displaystyle }\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}\equiv {\displaystyle \frac{\mathrm{d}\mathit{\zeta }}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{\sqrt{G_{\mathrm{v}}}}}\left(u^{\mathit{\zeta }}+\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{13}}{u}^{\mathit{\xi }}+\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{23}}{u}^{\mathit{\eta }}\right).\end{array}$$

The final Jacobian composed of horizontal and vertical coordinate transformations can be represented as $\sqrt{G}=\sqrt{G_{\mathrm{h}}}\sqrt{G_{\mathrm{v}}}$. Hereafter, to briefly describe the formulations, the coordinate variables are sometimes expressed using $({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}},{\mathit{\xi }}^{\mathrm{3}})=(\mathit{\xi },\mathit{\eta },\mathit{\zeta })$. In addition, the Einstein summation notation will be applied for repeated indices when representing the geometric relations.

The compact form of the governing equations can be written as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{q}}{\partial t}}& +{\displaystyle \frac{\partial \left[\mathit{f}\left(\mathit{q}\right)+{\mathit{f}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})\right]}{\partial \mathit{\xi }}}\\ & +{\displaystyle \frac{\partial \left[\mathit{g}\left(\mathit{q}\right)+{\mathit{g}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})\right]}{\partial \mathit{\eta }}}\\ & +{\displaystyle \frac{\partial \left[\mathit{h}\left(\mathit{q}\right)+{\mathit{h}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})\right]}{\partial \mathit{\zeta }}}\\ & =\mathit{S}\left(\mathit{q}\right)+{\mathit{S}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q}).\end{array}\end{array}$$

Here, q is the solution vector defined as

$$\begin{array}{}\text{(4)}& \mathit{q}={\left(\sqrt{G}{\mathit{\rho }}^{\prime },\sqrt{G}\mathit{\rho }{u}^{\mathit{\xi }},\sqrt{G}\mathit{\rho }{u}^{\mathit{\eta }},\sqrt{G}\mathit{\rho }{u}^{\mathit{\zeta }},\sqrt{G}(\mathit{\rho }\mathit{\theta }{)}^{\prime }\right)}^T,\end{array}$$

where ρ and θ denote the density and potential temperature, respectively. To accurately treat nearly balanced flows, the density ρ and pressure P (thus ρθ) are decomposed as $q(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)=q_{\mathrm{hyd}}(\mathit{\xi },\mathit{\eta },\mathit{\zeta })+q^{\prime }(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)$, where q_hyd denotes a variable satisfying the hydrostatic balance and q^′ denotes the deviation. In Eq. (3), f(q), g(q), and h(q) are inviscid fluxes in the ξ, η, and ζ directions, respectively. The horizontal inviscid fluxes are represented as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathit{f}\left(\mathit{q}\right)=\left(\begin{array}{c}\sqrt{G}\mathit{\rho }{u}^{\mathit{\xi }}\\ \sqrt{G}(\mathit{\rho }{u}^{\mathit{\xi }}{u}^{\mathit{\xi }}+G_{\mathrm{h}}^{\mathrm{11}}{P}^{\prime })\\ \sqrt{G}(\mathit{\rho }{u}^{\mathit{\eta }}{u}^{\mathit{\xi }}+G_{\mathrm{h}}^{\mathrm{21}}{P}^{\prime })\\ \sqrt{G}\mathit{\rho }{u}^{\mathit{\zeta }}{u}^{\mathit{\xi }}\\ \sqrt{G}\mathit{\rho }\mathit{\theta }{u}^{\mathit{\xi }}\end{array}\right),\\ & \mathit{g}\left(\mathit{q}\right)=\left(\begin{array}{c}\sqrt{G}\mathit{\rho }{u}^{\mathit{\eta }}\\ \sqrt{G}(\mathit{\rho }{u}^{\mathit{\xi }}{u}^{\mathit{\eta }}+G_{\mathrm{h}}^{\mathrm{12}}{P}^{\prime })\\ \sqrt{G}(\mathit{\rho }{u}^{\mathit{\eta }}{u}^{\mathit{\eta }}+G_{\mathrm{h}}^{\mathrm{22}}{P}^{\prime })\\ \sqrt{G}\mathit{\rho }{u}^{\mathit{\zeta }}{u}^{\mathit{\eta }}\\ \sqrt{G}\mathit{\rho }\mathit{\theta }{u}^{\mathit{\eta }}\end{array}\right),\end{array}\end{array}$$

and the vertical inviscid fluxes are represented as

$$\begin{array}{}\text{(6)}& \mathit{h}\left(\mathit{q}\right)=\left(\begin{array}{c}\sqrt{G}\mathit{\rho }\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}\\ \sqrt{G}[\mathit{\rho }{u}^{\mathit{\xi }}\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}+(G_{\mathrm{v}}^{\mathrm{13}}{G}_{\mathrm{h}}^{\mathrm{11}}+G_{\mathrm{v}}^{\mathrm{23}}{G}_{\mathrm{h}}^{\mathrm{12}}\left)P^{\prime }\right]\\ \sqrt{G}[\mathit{\rho }{u}^{\mathit{\eta }}\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}+(G_{\mathrm{v}}^{\mathrm{13}}{G}_{\mathrm{h}}^{\mathrm{21}}+G_{\mathrm{v}}^{\mathrm{23}}{G}_{\mathrm{h}}^{\mathrm{22}}\left)P^{\prime }\right]\\ \sqrt{G}\mathit{\rho }{u}^{\mathit{\zeta }}\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}+{\sqrt{G}}_{\mathrm{h}}{P}^{\prime }\\ \sqrt{G}\mathit{\rho }\mathit{\theta }\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}\end{array}\right).\end{array}$$

Furthermore, S(q) represents the source terms as

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{S}\left(\mathit{q}\right)=\left(\begin{array}{c}\mathrm{0}\\ \sqrt{G}(F_{\mathrm{H}}^{\mathrm{1}}+F_M^{\mathrm{1}}+F_{\mathrm{C}}^{\mathrm{1}})\\ \sqrt{G}(F_{\mathrm{H}}^{\mathrm{2}}+F_M^{\mathrm{2}}+F_{\mathrm{C}}^{\mathrm{2}})\\ \sqrt{G}(F_{\mathrm{buo}}+F_{\mathrm{C}}^{\mathrm{3}})\\ \mathrm{0}\end{array}\right),\end{array}$$

where $F_{\mathrm{H}}^i$ for i=1 and 2 is the horizontal pressure gradient term with hydrostatic balance and can be expressed as

$$\begin{array}{}\text{(8)}& {\displaystyle }{F}_{\mathrm{H}}^i=-{\displaystyle \frac{G_{\mathrm{h}}^{i{m}^{\prime }}}{\sqrt{G_{\mathrm{v}}}}}\left[{\displaystyle \frac{\partial \left(\sqrt{G_{\mathrm{v}}}{P}_{\mathrm{hyd}}\right)}{\partial {\mathit{\xi }}^{m^{\prime }}}}+{\displaystyle \frac{\partial \left(G_{\mathrm{v}}^{m^{\prime }\mathrm{3}}\sqrt{G_{\mathrm{v}}}{P}_{\mathrm{hyd}}\right)}{\partial {\mathit{\xi }}^{\mathrm{3}}}}\right].\end{array}$$

Here, note that $m^{\prime }=\mathrm{1}$ and 2. $F_M^i=-{\mathrm{\Gamma }}_{ml}^i(\mathit{\rho }{u}^m{u}^l+G^{ml}{P}^{\prime })$ is the source terms due to the horizontal curvilinear coordinate, where m and l take values of 1, 2, and 3 and ${\mathrm{\Gamma }}_{ml}^i$ is the Christoffel symbol of the second kind, which means the spatial variation in the basis vector. $F_{\mathrm{C}}^i=-G^{im}{\mathit{ϵ}}_{jml}\mathrm{2}{\mathrm{\Omega }}^m\mathit{\rho }{u}^l$ is the Coriolis terms, where ϵ_jkl is the three rank Levi-Civita tensor and Ω^m is the components of the angular-velocity vector. $F_{\mathrm{buo}}=-{\mathit{\rho }}^{\prime }(a/r{)}^{\mathrm{2}}g$ is the buoyancy term, where r is the radial coordinate, a is the planetary radius, and g is the standard gravitational acceleration. To close the equation systems, the pressure P is calculated using the state equation for the ideal gas as

$$\begin{array}{}\text{(9)}& {\displaystyle }P=P_{\mathrm{0}}{\left({\displaystyle \frac{R}{P_{\mathrm{0}}}}\mathit{\rho }\mathit{\theta }\right)}^{\frac{C_p}{C_{\mathrm{v}}}},\end{array}$$

where P₀ is a constant pressure; R is the gas constant; and C_v and C_p are the specific heat at constant volume and constant pressure, respectively. The actual values for these constants are provided in Table 1. In Eq. (3), the terms with an SGS subscript are associated with a turbulent model; f_SGS(q,∇q), g_SGS(q,∇q), and h_SGS(q,∇q) are the parameterized eddy fluxes, while S_SGS(q,∇q) is the source terms with the curvilinear coordinates. The terms associated with the turbulent model are detailed in Sect. 2.2.

Table 1 Values of parameters.

As a horizontal curvilinear coordinate, we adopted an equiangular gnomonic cubed-sphere projection (Sadourny, 1972; Ronchi et al., 1996) to map a cube onto a sphere. Compared to a conformal projection (Rančić et al., 1996), we preferred this projection for generating more uniform grids in high spatial resolutions, although the non-orthogonal basis needs to be treated. In each panel of the cube, a local coordinate using the central angles (α,β) ($\in [-\mathit{\pi }/\mathrm{4},\mathit{\pi }/\mathrm{4}]$) was introduced and related to the horizontal coordinates (ξ,η) by $\mathit{\xi }=\mathit{\alpha },\mathit{\eta }=\mathit{\beta }$. Based on the derivation with the coordinate transformation in previous studies (e.g., Nair et al., 2005b; Ullrich et al., 2012; Li et al., 2020), the horizontal contravariant metric tensor and the horizontal Jacobian for the equiangular gnomonic cubed-sphere projection can be written as, respectively,

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{G}_{\mathrm{c}}^{ij}& ={\displaystyle \frac{{\mathit{\delta }}^{\mathrm{2}}}{r^{\mathrm{2}}(\mathrm{1}+X^{\mathrm{2}})(\mathrm{1}+Y^{\mathrm{2}})}}\left(\begin{array}{cc}\mathrm{1}+Y^{\mathrm{2}}& XY\\ XY& \mathrm{1}+X^{\mathrm{2}}\end{array}\right)\text{ and}\\ & \sqrt{G_{\mathrm{c}}}={\displaystyle \frac{r^{\mathrm{2}}(\mathrm{1}+X^{\mathrm{2}})(\mathrm{1}+Y^{\mathrm{2}})}{{\mathit{\delta }}^{\mathrm{3}}}},\end{array}\end{array}$$

where X=tan α, Y=tan β, $\mathit{\delta }=\sqrt{\mathrm{1}+X^{\mathrm{2}}+Y^{\mathrm{2}}}$, and r is the radial coordinate. The Christoffel symbol of the second kind ${\mathrm{\Gamma }}_{ml}^i$ is represented as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\mathrm{\Gamma }}_{ml}^{\mathrm{1}}=\left(\begin{array}{ccc}\frac{\mathrm{2}X{Y}^{\mathrm{2}}}{{\mathit{\delta }}^{\mathrm{2}}}& \frac{-Y(\mathrm{1}+Y^{\mathrm{2}})}{{\mathit{\delta }}^{\mathrm{2}}}& \frac{{\mathit{\delta }}_{\mathrm{S}}}{r}\\ \frac{-Y(\mathrm{1}+Y^{\mathrm{2}})}{{\mathit{\delta }}^{\mathrm{2}}}& \mathrm{0}& \mathrm{0}\\ \frac{{\mathit{\delta }}_{\mathrm{S}}}{r}& \mathrm{0}& \mathrm{0}\end{array}\right),\\ & {\mathrm{\Gamma }}_{ml}^{\mathrm{2}}=\left(\begin{array}{ccc}\mathrm{0}& \frac{-X(\mathrm{1}+X^{\mathrm{2}})}{{\mathit{\delta }}^{\mathrm{2}}}& \mathrm{0}\\ \frac{-X(\mathrm{1}+X^{\mathrm{2}})}{{\mathit{\delta }}^{\mathrm{2}}}& \frac{\mathrm{2}{X}^{\mathrm{2}}Y}{{\mathit{\delta }}^{\mathrm{2}}}& \frac{{\mathit{\delta }}_{\mathrm{S}}}{r}\\ \mathrm{0}& \frac{{\mathit{\delta }}_{\mathrm{S}}}{r}& \mathrm{0}\end{array}\right),\\ & {\mathrm{\Gamma }}_{ml}^{\mathrm{3}}={\mathit{\delta }}_{\mathrm{S}}{\displaystyle \frac{r(\mathrm{1}+X^{\mathrm{2}})(\mathrm{1}+Y^{\mathrm{2}})}{{\mathit{\delta }}^{\mathrm{4}}}}\\ & \left(\begin{array}{ccc}-(\mathrm{1}+X^{\mathrm{2}})& XY& \mathrm{0}\\ XY& -(\mathrm{1}+Y^{\mathrm{2}})& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right),\end{array}\end{array}$$

where δ_S is a switch for shallow-atmosphere approximation. The components of the angular-velocity vector included in the Coriolis terms $F_{\mathrm{C}}^i$ are given as

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\mathrm{\Omega }}^{\mathrm{1}}=\mathrm{0},{\mathrm{\Omega }}^{\mathrm{2}}={\mathit{\delta }}_{\mathrm{S}}{\displaystyle \frac{\mathit{\omega }\mathit{\delta }}{r(\mathrm{1}+Y^{\mathrm{2}})}},\\ & {\mathrm{\Omega }}^{\mathrm{3}}=\mathit{\omega }{\displaystyle \frac{Y}{\mathit{\delta }}},\text{for the equatorial panels, and}\\ & {\mathrm{\Omega }}^{\mathrm{1}}=-{\mathit{\delta }}_{\mathrm{S}}{\displaystyle \frac{s\mathit{\omega }X\mathit{\delta }}{r(\mathrm{1}+X^{\mathrm{2}})}},{\mathrm{\Omega }}^{\mathrm{2}}=-{\mathit{\delta }}_{\mathrm{S}}{\displaystyle \frac{s\mathit{\omega }Y\mathit{\delta }}{r(\mathrm{1}+Y^{\mathrm{2}})}},\\ & {\mathrm{\Omega }}^{\mathrm{3}}={\displaystyle \frac{s\mathit{\omega }}{\mathit{\delta }}},\text{for the polar panels},\end{array}\end{array}$$

where ω is the angular velocity of the planet and an index s has a value of 1 and −1 for the northern and southern polar panels, respectively. In the numerical experiments in Sect. 3, shallow-atmosphere approximation was applied. Then, r and δ_S are treated as follows: the radial coordinate r in Eqs. (10)–(12) and the buoyancy term in Eq. (6) are replaced by the planetary radius a. In Eqs. (11) and (12), the terms with δ_S are ignored. In addition, the pressure contribution in $F_M^i$ disappears because the relation of ${\mathrm{\Gamma }}_{ml}^i{G}^{ml}=\mathrm{0}$ is satisfied in shallow-atmosphere approximation.

To treat the topography, the traditional terrain-following coordinate (Phillips, 1957; Gal-Chen and Somerville, 1975) was adopted as a general vertical coordinate. The vertical coordinate transformation can be expressed as

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\zeta }=z_{\mathrm{T}}{\displaystyle \frac{z-h}{z_{\mathrm{T}}-h}},\end{array}$$

where z is the height coordinate, z_T is the top height of the computational domain (we assume it is a constant value), and h is the surface height. The corresponding Jacobian and metric tensor can be represented as

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \sqrt{G_{\mathrm{v}}}=\mathrm{1}-{\displaystyle \frac{h}{z_{\mathrm{T}}}},\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{13}}=\left({\displaystyle \frac{\mathit{\zeta }}{z_{\mathrm{T}}}}-\mathrm{1}\right){\displaystyle \frac{\partial h}{\partial \mathit{\xi }}}\text{ and}\\ & \sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{23}}=\left({\displaystyle \frac{\mathit{\zeta }}{z_{\mathrm{T}}}}-\mathrm{1}\right){\displaystyle \frac{\partial h}{\partial \mathit{\eta }}},\end{array}\end{array}$$

respectively.

2.2  Formulation of eddy viscous and diffusion terms in general curvilinear coordinates

Considering global LES in future high-resolution simulations, this subsection describes eddy viscous and diffusion terms in the general curvilinear coordinates. We utilized a Smagorinsky–Lilly type model (Smagorinsky, 1963; Lilly, 1962) that considered the stratification effect (Brown et al., 1994). This turbulent model was also used in KT2021 and KT2023. As a spatial filter, Favre filtering (Favre, 1983) was used. We do not explicitly denote the symbol representing the spatial filter because the filtering approach is essentially the same as that explained in Appendix A of KT2023. The difficulties in the derivation of viscous and diffusion terms are caused by the gradient of vector quantities and the spatial divergence with the non-orthogonal basis because the manipulations grow increasingly complex. However, previous studies that utilized tensor analysis help us provide a systematic derivation (e.g., Ullrich, 2014; Rančić et al., 2017). In the absence of a vertical coordinate transformation, the parameterized fluxes with the turbulent model can be represented in the general curvilinear coordinates as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& {\mathit{f}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})=\left(\begin{array}{c}\mathrm{0}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{11}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{12}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{13}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}_{*}^{\mathrm{1}}\end{array}\right),\\ & {\mathit{g}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})=\left(\begin{array}{c}\mathrm{0}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{21}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{22}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{23}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}_{*}^{\mathrm{2}}\end{array}\right),\\ & {\mathit{h}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})=\left(\begin{array}{c}\mathrm{0}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{31}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{32}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}^{\mathrm{33}}\\ -\sqrt{G}\mathit{\rho }{\mathit{\tau }}_{*}^{\mathrm{3}}\end{array}\right),\end{array}\end{array}$$

and the source term can be given by

$$\begin{array}{}\text{(16)}& {\mathit{S}}_{\mathrm{SGS}}(\mathit{q},\mathrm{\nabla }\mathit{q})=\left(\begin{array}{c}\mathrm{0}\\ -\sqrt{G}{\mathrm{\Gamma }}_{ml}^{\mathrm{1}}\mathit{\rho }{\mathit{\tau }}^{ml}\\ -\sqrt{G}{\mathrm{\Gamma }}_{ml}^{\mathrm{2}}\mathit{\rho }{\mathit{\tau }}^{ml}\\ -\sqrt{G}{\mathrm{\Gamma }}_{ml}^{\mathrm{3}}\mathit{\rho }{\mathit{\tau }}^{ml}\\ \mathrm{0}\end{array}\right).\end{array}$$

In the equations, τ^ij is the contravariant components of parameterized eddy viscous flux tensor (i=1, 2, and 3 and j=1, 2, and 3) and can be written as

$$\begin{array}{}\text{(17)}& {\mathit{\tau }}^{ij}=-\mathrm{2}{\mathit{\nu }}_{\mathrm{SGS}}\left(S^{ij}-{\displaystyle \frac{G^{ij}}{\mathrm{3}}}D\right)-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{G}^{ij}{K}_{\mathrm{SGS}},\end{array}$$

where S^ij is the strain velocity tensor, ν_SGS is the eddy viscosity, D is the divergence of the three-dimensional velocity, and K_SGS is the SGS kinetic energy. The strain velocity tensor is represented as

$$\begin{array}{}\text{(18)}& S^{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(G^{im}{\displaystyle \frac{\partial u_{,m}^j}{\partial {\mathit{\xi }}^m}}+G^{jm}{\displaystyle \frac{\partial u_{,m}^i}{\partial {\mathit{\xi }}^m}}\right),\end{array}$$

using the covariant derivative of the contravariant velocity component

$$\begin{array}{}\text{(19)}& u_{,j}^i={\displaystyle \frac{\partial u^i}{\partial {\mathit{\xi }}^j}}+u^m{\mathrm{\Gamma }}_{jm}^i.\end{array}$$

The eddy viscosity is written as

$$\begin{array}{}\text{(20)}& {\mathit{\nu }}_{\mathrm{SGS}}=C_{\mathrm{s}}{\mathrm{\Delta }}_{\mathrm{SGS}}\left|S\right|,\end{array}$$

where C_s, Δ_SGS, and $\left|S\right|$ represent the Smagorinsky constant, the filter length, and the norm of strain tensor defined as $\sqrt{\mathrm{2}{G}_{im}{G}_{jn}{S}^{ij}{S}^{mn}}$, respectively. The parameterized eddy diffusive flux can be written as

$$\begin{array}{}\text{(21)}& {\mathit{\tau }}_{*}^i=-{\mathit{\nu }}_{\mathrm{SGS}}^{*}{G}^{ij}{\displaystyle \frac{\partial \mathit{\theta }}{\partial {\mathit{\xi }}^j}},\end{array}$$

where ${\mathit{\nu }}_{\mathrm{SGS}}^{*}$ is the eddy diffusion coefficient. For further details of the turbulent model, refer to Sect. 2.2 of Nishizawa et al. (2015).

2.3 Spatial discretization

We performed the spatial discretization for Eq. (3) based on a nodal DGM (e.g., Hesthaven and Warburton, 2007). In each cubed-sphere panel, the three-dimensional computational domain Ω was divided using non-overlapping hexahedral elements. To relate the coordinates $({\mathit{\xi }}^{\mathrm{1}},{\mathit{\xi }}^{\mathrm{2}},{\mathit{\xi }}^{\mathrm{3}})=(\mathit{\alpha },\mathit{\beta },\mathit{\zeta })$ with the local coordinates $\stackrel{\mathrm{̃}}{\mathit{x}}\equiv ({\stackrel{\mathrm{̃}}{x}}^{\mathrm{1}},{\stackrel{\mathrm{̃}}{x}}^{\mathrm{2}},{\stackrel{\mathrm{̃}}{x}}^{\mathrm{3}})$ in a reference element Ω_e, we adopted a linear mapping defined as

$$\begin{array}{}\text{(22)}& {\stackrel{\mathrm{̃}}{x}}^i=\mathrm{2}{\displaystyle \frac{{\mathit{\xi }}^i-{\mathit{\xi }}_{\mathrm{e}}^i}{h_{\mathrm{e}}^i}},\end{array}$$

where ${\mathit{\xi }}_{\mathrm{e}}^i$ and $h_{\mathrm{e}}^i$ represent the center position and width of the element, respectively, in the ξⁱ direction. By equally dividing the (α,β) plane, we generated a horizontal mesh including the $N_{\mathrm{e},\mathrm{h}}\times N_{\mathrm{e},\mathrm{h}}$ finite elements. The center horizontal position of $(i^{\prime },j^{\prime })$th element can be written as

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{\mathit{\alpha }}_{i^{\prime }}& =-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}{N}_{\mathrm{e},\mathrm{h}}}}\left(i^{\prime }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right),\\ & {\mathit{\beta }}_{j^{\prime }}=-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}{N}_{\mathrm{e},\mathrm{h}}}}\left(j^{\prime }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right).\end{array}\end{array}$$

Using the tensor product of one-dimensional Lagrange polynomials $l_m\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)=l_{m_{\mathrm{1}}}\left({\stackrel{\mathrm{̃}}{x}}^{\mathrm{1}}\right)l_{m_{\mathrm{2}}}\left({\stackrel{\mathrm{̃}}{x}}^{\mathrm{2}}\right)l_{m_{\mathrm{3}}}\left({\stackrel{\mathrm{̃}}{x}}^{\mathrm{3}}\right)$, a local approximated solution within each element Ω_e can be represented as

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{\mathit{q}}^{\mathrm{e}}{|}_{{\mathrm{\Omega }}_{\mathrm{e}}}(\stackrel{\mathrm{̃}}{\mathit{x}},t)& =\sum _{m_{\mathrm{1}}=\mathrm{1}}^{p+\mathrm{1}}\sum _{m_{\mathrm{2}}=\mathrm{1}}^{p+\mathrm{1}}\sum _{m_{\mathrm{3}}=\mathrm{1}}^{p+\mathrm{1}}{\mathit{Q}}_{m_{\mathrm{1}},m_{\mathrm{2}},m_{\mathrm{3}}}^{\mathrm{e}}\left(t\right)\\ & l_{m_{\mathrm{1}}}\left({\stackrel{\mathrm{̃}}{x}}^{\mathrm{1}}\right)l_{m_{\mathrm{2}}}\left({\stackrel{\mathrm{̃}}{x}}^{\mathrm{2}}\right)l_{m_{\mathrm{3}}}\left({\stackrel{\mathrm{̃}}{x}}^{\mathrm{3}}\right).\end{array}\end{array}$$

In Eq. (24), the coefficients $Q_{m_{\mathrm{1}},m_{\mathrm{2}},m_{\mathrm{3}}}^{\mathrm{e}}$ are the unknown degrees of freedom (DOFs) and p is the polynomial order. In this study, the Legendre–Gauss–Lobatto (LGL) points were used for interpolation and integration nodes. We defined a representative grid spacing at the equator, which approximately corresponds to that in the grid-point methods as

$$\begin{array}{}\text{(25)}& {\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}={\displaystyle \frac{\mathit{\pi }a}{\mathrm{2}{N}_{\mathrm{e},\mathrm{h}}(p+\mathrm{1})}}.\end{array}$$

Similar to the case of the horizontal direction, we defined a representative vertical grid spacing Δ_v. For the uniform vertical-element size, ${\mathrm{\Delta }}_{\mathrm{v}}=z_{\mathrm{T}}/\left(N_{\mathrm{e},\mathrm{v}}\right(p+\mathrm{1}\left)\right)$, where N_e,v is the number of vertical elements. Hereinafter, we simply refer to Δ_h,eq and Δ_v as the horizontal and vertical grid spacing, respectively.

By applying the Galerkin approximation to Eq. (3), a strong form of the semi-discretized equations can be obtained as

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\mathit{q}}^{\mathrm{e}}(\stackrel{\mathrm{̃}}{\mathit{x}},t)l_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}=\\ & -\sum _{j=\mathrm{1}}^{\mathrm{3}}\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\displaystyle \frac{\partial {\mathit{F}}_j({\mathit{q}}^{\mathrm{e}},\mathit{G})}{\partial {\mathit{\xi }}^j}}{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}\\ & -\underset{\partial {\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left[\widehat{\mathit{F}}({\mathit{q}}^{\mathrm{e}},\mathit{G})-\mathit{F}({\mathit{q}}^{\mathrm{e}},\mathit{G})\right]\cdot \mathit{n}{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\partial \mathrm{E}}\mathrm{d}S\\ & +\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left[\mathit{S}\left({\mathit{q}}^{\mathrm{e}}\right)+{\mathit{S}}_{\mathrm{SGS}}({\mathit{q}}^{\mathrm{e}},\mathit{G})\right]l_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}},\end{array}\end{array}$$

where $({\mathit{F}}_{\mathrm{1}},{\mathit{F}}_{\mathrm{2}},{\mathit{F}}_{\mathrm{3}})=(\mathit{f}+{\mathit{f}}_{\mathrm{SGS}},\mathit{g}+{\mathit{g}}_{\mathrm{SGS}},\mathit{h}+{\mathit{h}}_{\mathrm{SGS}})$ is the flux vector tensor, $\widehat{\mathit{F}}$ is the numerical flux at the element boundary ∂Ω_E, and n is the outward unit vector normal to ∂Ω_E. In the volume and surface integrals, J^E and J^∂E represent the transformation Jacobian with the general curvilinear coordinates and local coordinates within each element. Note that, because of the linear mapping in Eq. (22), the associated geometric factors such as J^E and J^∂E have constant values when the volume and surface integrals are calculated. For the turbulent model, we need to evaluate the eddy viscous flux tensor and diffusion flux, which include a few gradient terms with quantities such as $\mathit{\chi }=(u^{\mathit{\xi }},u^{\mathit{\eta }},u^{\mathit{\zeta }},\mathit{\theta })$, denoted by $\mathit{G}=(\partial \mathit{\chi }/\partial {\mathit{\xi }}^{\mathrm{1}},\partial \mathit{\chi }/\partial {\mathit{\xi }}^{\mathrm{2}},\partial \mathit{\chi }/\partial {\mathit{\xi }}^{\mathrm{3}})$ in Eq. (26). The gradient discretization in the ξ^j direction is given by

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\mathit{\rho }{\mathit{G}}_j{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}=\\ & \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left[{\displaystyle \frac{\partial {\mathit{\rho }}^{\mathrm{e}}{\mathit{\chi }}^{\mathrm{e}}}{\partial {\mathit{\xi }}^j}}-{\mathit{\chi }}^{\mathrm{e}}{\left({\displaystyle \frac{\partial \mathit{\rho }}{\partial {\mathit{\xi }}^j}}\right)}^{\mathrm{e}}\right]l_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}\\ & +\underset{\partial {\mathrm{\Omega }}_{\mathrm{e}}}{\int }\left(\widehat{\mathit{\rho }\mathit{\chi }}-{\mathit{\rho }}^{\mathrm{e}}{\mathit{\chi }}^{\mathrm{e}}\right){\mathit{n}}_{{\stackrel{\mathrm{̃}}{x}}^j}\cdot \mathit{n}{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\partial \mathrm{E}}\mathrm{d}S,\end{array}\end{array}$$

where ${\mathit{n}}_{{\stackrel{\mathrm{̃}}{\mathit{x}}}^j}$ is the unit vector in the ${\stackrel{\mathrm{̃}}{x}}^j$ direction and the density gradient is calculated by

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& \underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\left({\displaystyle \frac{\partial \mathit{\rho }}{\partial {\mathit{\xi }}^j}}\right)}^{\mathrm{e}}{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}=\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{\displaystyle \frac{\partial {\mathit{\rho }}^{\mathrm{e}}}{\partial {\mathit{\xi }}^j}}{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\mathrm{E}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}\\ & +\underset{\partial {\mathrm{\Omega }}_{\mathrm{e}}}{\int }(\widehat{\mathit{\rho }}-{\mathit{\rho }}^{\mathrm{e}}){\mathit{n}}_{{\stackrel{\mathrm{̃}}{x}}^j}\cdot \mathit{n}{l}_{\mathit{m}}\left(\stackrel{\mathrm{̃}}{\mathit{x}}\right)J^{\partial \mathrm{E}}\mathrm{d}S.\end{array}\end{array}$$

For the numerical flux of the inviscid terms, the Rusanov flux (Rusanov, 1961) was used as a simple choice of the approximated Riemann solvers. Its numerical dissipation is provided based on the maximum absolute eigenvalue of the Jacobian matrix on the left and right sides of the element boundary. Previous studies (e.g., Li et al., 2020) formulated the Rusanov flux, taking into account the horizontal and vertical coordinate transformations, as

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{\widehat{\mathit{F}}}_{\mathrm{invis}}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left\{\left[{\mathit{F}}_{\mathrm{invis}}\left({\mathit{q}}^{+}\right)+{\mathit{F}}_{\mathrm{invis}}\left({\mathit{q}}^{-}\right)\right]\right.\\ & \left.\cdot \mathit{n}-{\mathit{\lambda }}_{\max }\left[{\mathit{q}}^{+}-{\mathit{q}}^{-}\right]\right\},\end{array}\end{array}$$

where λ_max is the maximum of the absolute value of eigenvalues of the flux Jacobian in the direction n and q⁻ and q⁺ represent the interior and exterior values at ∂Ω_e. At the element boundaries in the horizontal directions (ξ and η), λ_max can be represented as

$$\begin{array}{}\text{(30)}& {\mathit{\lambda }}_{\max ,\mathit{\xi }}=\left|u^{\mathit{\xi }}\right|+\sqrt{G_{\mathrm{c}}^{\mathrm{11}}}{c}_{\mathrm{s}},{\mathit{\lambda }}_{\max ,\mathit{\eta }}=\left|u^{\mathit{\eta }}\right|+\sqrt{G_{\mathrm{c}}^{\mathrm{22}}}{c}_{\mathrm{s}},\end{array}$$

where $c_{\mathrm{s}}=\left[\right(C_p/C_{\mathrm{v}})RT{]}^{\mathrm{1}/\mathrm{2}}$ is the speed of a sound wave. For the vertical direction ζ, λ_max can be represented as

$$\begin{array}{}\text{(31)}& {\mathit{\lambda }}_{\max ,\mathit{\zeta }}=\left|\stackrel{\mathrm{̃}}{u^{\mathit{\zeta }}}\right|+{\left[\mathrm{1}/\sqrt{G_{\mathrm{v}}}+G_{\mathrm{v}}^{\mathrm{13}}{G}_X+G_{\mathrm{v}}^{\mathrm{23}}{G}_Y\right]}^{\mathrm{1}/\mathrm{2}}{c}_{\mathrm{s}},\end{array}$$

where $G_X=G_{\mathrm{v}}^{\mathrm{13}}{G}_{\mathrm{c}}^{\mathrm{11}}+G_{\mathrm{v}}^{\mathrm{23}}{G}_{\mathrm{c}}^{\mathrm{12}}$ and $G_Y=G_{\mathrm{v}}^{\mathrm{13}}{G}_{\mathrm{c}}^{\mathrm{21}}+G_{\mathrm{v}}^{\mathrm{23}}{G}_{\mathrm{c}}^{\mathrm{22}}$. The central flux was adopted as the numerical flux of the gradient G and the SGS fluxes $({\mathit{f}}_{\mathrm{SGS}},{\mathit{g}}_{\mathrm{SGS}},{\mathit{h}}_{\mathrm{SGS}})$ with the turbulent model.

When the same nodes are used for interpolation and integration (i.e., collocation), a matrix form of Eqs. (26) and (27) can be obtained as

$$\begin{array}{}\text{(32)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{d{\mathit{q}}^{\mathrm{e}}}{\mathrm{d}t}}=& -\sum _{j=\mathrm{1}}^{\mathrm{3}}{d}_j{D}_{{\stackrel{\mathrm{̃}}{x}}^j}{\mathit{F}}_j({\mathit{q}}^{\mathrm{e}},\mathit{G})-\sum _{f=\mathrm{1}}^{\mathrm{6}}{s}_{\partial {\mathrm{\Omega }}_{\mathrm{e},f}}{L}_{\partial {\mathrm{\Omega }}_{\mathrm{e},f}}\\ & \left[\widehat{\mathit{F}}({\mathit{q}}^{\mathrm{e}},\mathit{G})-\mathit{F}({\mathit{q}}^{\mathrm{e}},\mathit{G})\right]\cdot \mathit{n}\\ & +\mathit{S}\left({\mathit{q}}^{\mathrm{e}}\right)+{\mathit{S}}_{\mathrm{SGS}}({\mathit{q}}^{\mathrm{e}},\mathit{G}),\end{array}\text{(33)}& {\displaystyle }\begin{array}{rl}\mathit{\rho }{\mathit{G}}_j=& d_j{D}_{{\stackrel{\mathrm{̃}}{x}}^j}\left({\mathit{\rho }}^{\mathrm{e}}{\mathit{\chi }}^{\mathrm{e}}\right)-{\mathit{\chi }}^{\mathrm{e}}{\left({\displaystyle \frac{\partial \mathit{\rho }}{\partial {\mathit{\xi }}^j}}\right)}^{\mathrm{e}}\\ & +\sum _{f^{\prime }=\mathrm{1}}^{\mathrm{2}}{s}_{\partial {\mathrm{\Omega }}_{\mathrm{e},f^{\prime }}}{L}_{\partial {\mathrm{\Omega }}_{\mathrm{e},f^{\prime }}}\left(\widehat{\mathit{\rho }\mathit{\chi }}-{\mathit{\rho }}^{\mathrm{e}}{\mathit{\chi }}^{\mathrm{e}}\right){\mathit{n}}_{{\stackrel{\mathrm{̃}}{x}}^j}\cdot \mathit{n},\end{array}\end{array}$$

where $D_{{\stackrel{\mathrm{̃}}{x}}^j}$ represents the differential matrix for the ${\stackrel{\mathrm{̃}}{x}}^j$ direction, $L_{\partial {\mathrm{\Omega }}_{\mathrm{e},f}}$ represents the lifting matrix with the surface integral for the fth element surface, and $L_{\partial {\mathrm{\Omega }}_{\mathrm{e},f^{\prime }}}$ represents the same for the f^′th element surface in the gradient operator for the ${\stackrel{\mathrm{̃}}{x}}^j$ direction. The components of these matrices are given as

$$\begin{array}{}\text{(34)}& \begin{array}{rl}& {\left(D_{{\stackrel{\mathrm{̃}}{x}}^j}\right)}_{\mathit{m},{\mathit{m}}^{\prime }}=M^{-\mathrm{1}}\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{l}_{\mathit{m}}{\displaystyle \frac{\partial l_{{\mathit{m}}^{\prime }}}{\partial {\stackrel{\mathrm{̃}}{x}}_j}}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}},\\ & {\left(L_{\partial {\mathrm{\Omega }}_{\mathrm{e}},j}\right)}_{\mathit{m},{\mathit{m}}^{\prime }}=M^{-\mathrm{1}}\underset{\partial {\mathrm{\Omega }}_{\mathrm{e},j}}{\int }{l}_{\mathit{m}}{l}_{{\mathit{m}}^{\prime }}^{\partial {\mathrm{\Omega }}_{\mathrm{e},j}}\mathrm{d}S,\end{array}\end{array}$$

where M denotes the mass matrix and is given by

$$\begin{array}{}\text{(35)}& M_{\mathit{m},{\mathit{m}}^{\prime }}=\underset{{\mathrm{\Omega }}_{\mathrm{e}}}{\int }{l}_{\mathit{m}}{l}_{{\mathit{m}}^{\prime }}\mathrm{d}\stackrel{\mathrm{̃}}{\mathit{x}}.\end{array}$$

The density gradient term is calculated as

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{\left({\displaystyle \frac{\partial \mathit{\rho }}{\partial {\mathit{\xi }}^j}}\right)}^{\mathrm{e}}& =d_j{D}_{{\stackrel{\mathrm{̃}}{x}}^j}{\mathit{\rho }}^{\mathrm{e}}\\ & -\sum _{f^{\prime }=\mathrm{1}}^{\mathrm{2}}{s}_{\partial {\mathrm{\Omega }}_{\mathrm{e},f^{\prime }}}{L}_{\partial {\mathrm{\Omega }}_{\mathrm{e},f^{\prime }}}\left(\widehat{\mathit{\rho }}-{\mathit{\rho }}^{\mathrm{e}}\right){\mathit{n}}_{{\stackrel{\mathrm{̃}}{x}}^j}\cdot \mathit{n}.\end{array}\end{array}$$

Note that, in Eqs. (32), (33), and (36), $d_j=\partial {\stackrel{\mathrm{̃}}{x}}^j/\partial {\mathit{\xi }}^j$ and $s_{\partial {\mathrm{\Omega }}_{\mathrm{e}},f^{\prime }}=J_{\partial {\mathrm{\Omega }}_{\mathrm{e},f^{\prime }}}/J^{\mathrm{E}}$ are constant values in the volume and surface integrals, respectively. We changed the calculation method of mass and lifting matrices depending on temporal discretization. This is detailed in Sect. 2.4.

The balance between the pressure gradient and buoyancy terms should be carefully treated in the discrete momentum equation (e.g., Blaise et al., 2016; Orgis et al., 2017). In the above formulation, because a different discretization space is used between the terms, a numerical imbalance is possible and may cause spurious oscillations, which can destabilize the simulations. To avoid this incompatibility, the vertical polynomial order for the density in the buoyancy term was reduced by one following Blaise et al. (2016).

2.4 Temporal discretization

The semi-discretized equations in Eq. (26) can be represented as the following ordinary differential equation (ODE) system as

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\mathrm{d}\mathit{q}}{\mathrm{d}t}}=\mathcal{S}(\mathit{q},\mathrm{\nabla }\mathit{q})+\mathcal{F}(\mathit{q},\mathrm{\nabla }\mathit{q}),\end{array}$$

where 𝒮(q,∇q) and ℱ(q,∇q) represent the tendencies with slow and fast contributions, respectively. This study adopted Runge–Kutta (RK) schemes to solve the ODE system from t=nΔt to $t=(n+\mathrm{1})\mathrm{\Delta }t$, where Δt is the time step and n is a natural number. In this subsection, we describe two approaches for temporal discretization, namely, the horizontally explicit and vertically implicit (HEVI) and horizontally explicit and vertically explicit (HEVE) approaches.

We introduce two types of Courant number, which are used to explain time step settings for the numerical experiments in Sect. 3. For the horizontal advection test, the advective Courant number associated with the horizontal wind is defined as $C_{\mathrm{r},\mathrm{adv}}=U_{\mathrm{0}}\mathrm{\Delta }t/{\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}$, where U₀ is the representative wind speed. For other numerical experiments, the acoustic Courant number associated with the sound wave propagation is defined as $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=c_{\mathrm{s}}\mathrm{\Delta }t/\mathrm{\Delta }$, where Δ is the grid spacing. In particular, for the HEVI approach, $\mathrm{\Delta }={\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}$.

2.4.1 HEVI approach

If the aspect ratio of horizontal grid spacing to its vertical counterpart is large, it is impractical to use fully explicit temporal schemes because the vertically propagating sound waves severely restrict the time step. A strategy to avoid computational cost in such a case is the HEVI approach. The terms corresponding to vertical dynamics with a fast timescale are evaluated using an implicit temporal scheme, while the remaining terms are evaluated using an explicit temporal scheme. This procedure is regarded as a framework of an implicit–explicit (IMEX) time integration scheme (Bao et al., 2015; Gardner et al., 2018). General formulation of the IMEX RK scheme (e.g., Ascher et al., 1997) with ν stages can be represented as

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{\mathit{q}}^{\left(s\right)}& ={\mathit{q}}^n+\mathrm{\Delta }t\sum _{s^{\prime }=\mathrm{1}}^{s-\mathrm{1}}{a}_{s{s}^{\prime }}\mathcal{S}(t+c_{s^{\prime }}\mathrm{\Delta }t,{\mathit{q}}^{\left(s^{\prime }\right)})\\ & +\mathrm{\Delta }t\sum _{s^{\prime }=\mathrm{1}}^s{\stackrel{\mathrm{̃}}{a}}_{s{s}^{\prime }}\mathcal{F}(t+{\stackrel{\mathrm{̃}}{c}}_{s^{\prime }}\mathrm{\Delta }t,{\mathit{q}}^{\left(s^{\prime }\right)})\mathrm{for}s=\mathrm{1},\mathrm{\dots },\mathit{\nu },\\ {\mathit{q}}^{n+\mathrm{1}}& ={\mathit{q}}^n+\mathrm{\Delta }t\sum _{s=\mathrm{1}}^{\mathit{\nu }}{b}_s\mathcal{S}(t+c_s\mathrm{\Delta }t,{\mathit{q}}^{\left(s\right)})\\ & +\mathrm{\Delta }t\sum _{s=\mathrm{1}}^{\mathit{\nu }}{\stackrel{\mathrm{̃}}{b}}_s\mathcal{F}(t+{\stackrel{\mathrm{̃}}{c}}_s\mathrm{\Delta }t,{\mathit{q}}^{\left(s\right)}),\end{array}\end{array}$$

where $a_{s{s}^{\prime }}$, b_s, and c_s define the explicit temporal integrator; ${\stackrel{\mathrm{̃}}{a}}_{s{s}^{\prime }}$, ${\stackrel{\mathrm{̃}}{b}}_s$, and ${\stackrel{\mathrm{̃}}{c}}_{s^{\prime }}$ define the implicit temporal integrator; and $c_{\mathrm{s}}={\sum }_{s^{\prime }=\mathrm{1}}^{s-\mathrm{1}}{a}_{s{s}^{\prime }}$ and ${\stackrel{\mathrm{̃}}{c}}_s={\sum }_{s^{\prime }=\mathrm{1}}^{s-\mathrm{1}}{\stackrel{\mathrm{̃}}{a}}_{s{s}^{\prime }}$ represent time when slow and fast terms are evaluated, respectively. These coefficients are compactly represented using a “double Butcher tableau”, as shown in Table 2. Note that, in the table of the explicit part, $\mathcal{A}=\mathit{\left\{}{a}_{s{s}^{\prime }}\mathit{\right\}}$ with $a_{s{s}^{\prime }}=\mathrm{0}$ for $s^{\prime }\ge s$. On the other hand, for the implicit part, $\stackrel{\mathrm{̃}}{\mathcal{A}}=\mathit{\left\{}{\stackrel{\mathrm{̃}}{a}}_{s{s}^{\prime }}\mathit{\right\}}$ with ${\stackrel{\mathrm{̃}}{a}}_{s{s}^{\prime }}=\mathrm{0}$ for $s^{\prime }>s$ in the case of the diagonally implicit RK scheme.

Table 2 Double Butcher tableau for a third-order IMEX RK scheme proposed by Kennedy and Carpenter (2003).

In this study, the terms associated with vertical mass flux, vertical pressure gradient, vertical flux of potential temperature, and buoyancy in Eq. (3) were treated as fast terms, whereas the other terms were treated as slow terms. To minimize contaminating the spatial accuracy of a high-order DGM by temporal errors present in a low-order HEVI scheme, a third-order scheme proposed by Kennedy and Carpenter (2003) was adopted; it includes four explicit and three implicit evaluations. The corresponding double Butcher tableau is given in Table 2. In the implicit part of each stage, the corresponding nonlinear equation system is solved using Newton's method. In each iteration, the linearized equation system is solved. Obtaining accurate solutions of the nonlinear equation system generally requires numerous iterations. However, this study performed a single iteration in Newton's method (i.e., Rosenbrock approach), significantly reducing the computational cost. A similar approach has been used in previous studies (Ullrich and Jablonowski, 2012a). In the case of the collocation approach, because the horizontal dependency between all nodes within an element vanishes, the vertically implicit evaluation can be parallelly performed at each horizontal node.

For the case of HEVI, the volume and surface integrations in Eqs. (34) and (35) were evaluated using inexact integration with the LGL nodes. Consequently, M and ${\mathit{L}}_{\partial {\mathrm{\Omega }}_{\mathrm{e},\mathrm{3}}}$ became diagonal matrices, which further simplified the matrix structure associated with the vertical spatial operator.

2.4.2 HEVE approach

When we consider a horizontal grid spacing with O(10 m) such as in LES, the ratio of horizontal to vertical grid spacing approaches unity. The advantages of the HEVI approach decrease. Thus, it is suitable to adopt a fully explicit temporal approach, referred to as the HEVE approach. In such cases, RK schemes with a strong-stability-preserving (SSP) property (Gottlieb et al., 2001) are often used in combination with the DGM. Similar to KT2023, this study adopted a 10-stage RK scheme with the fourth-order accuracy proposed by Ketcheson (2008). The corresponding Butcher tableau is given in Table 3. When using the HEVE approach, entries of the matrices in Eqs. (34) and (35) were directly calculated following Sect. 3.2 in Hesthaven and Warburton (2007).

Table 3 Butcher tableau for a fourth-order fully explicit RK scheme with 10 stages proposed by Ketcheson (2008).

2.5 Modal filtering

For a high-order DGM, numerical instability is likely to occur in advection-dominated flows because the numerical dissipations with the upwind numerical fluxes weaken. Furthermore, we adopted a collocation approach due to its computational efficiency, as described in Sect. 2.3. One drawback is that the aliasing errors with evaluations of the nonlinear terms can drive numerical instability. To suppress this numerical instability, a modal filter was used as an additional stabilization mechanism. The filter matrix for the three-dimensional problem can be obtained as

$$\begin{array}{}\text{(39)}& \mathcal{F}=V^{\mathrm{3}\mathrm{D}}{C}^{\mathrm{3}\mathrm{D}}{V}^{\mathrm{3}\mathrm{D}},\end{array}$$

where V^3D represents the Vandermonde matrix associated with the LGL interpolation nodes (in Eq. 24) and C^3D represents the diagonal cutoff matrix. The entries of C^3D are defined as

$$\begin{array}{}\text{(40)}& C_{(m_{\mathrm{1}},m_{\mathrm{2}},m_{\mathrm{3}}),(m_{\mathrm{1}}^{\prime },m_{\mathrm{2}}^{\prime },m_{\mathrm{3}}^{\prime })}^{\mathrm{3}\mathrm{D}}={\mathit{\delta }}_{m_{\mathrm{1}},m_{\mathrm{1}}^{\prime }}{\mathit{\sigma }}_{m_{\mathrm{1}}}^{\mathrm{h}}{\mathit{\delta }}_{m_{\mathrm{2}},m_{\mathrm{2}}^{\prime }}{\mathit{\sigma }}_{m_{\mathrm{2}}}^{\mathrm{h}}{\mathit{\delta }}_{m_{\mathrm{3}},m_{\mathrm{3}}^{\prime }}{\mathit{\sigma }}_{m_{\mathrm{3}}}^{\mathrm{v}},\end{array}$$

where ${\mathit{\sigma }}_i^{\mathrm{h}}$ and ${\mathit{\sigma }}_i^{\mathrm{v}}$ represent the decay coefficient for the one-dimensional horizontal and vertical modes i, respectively. Based on Hesthaven and Warburton (2007), a typical choice of the coefficient for mode i is provided with an exponential function as

$$\begin{array}{}\text{(41)}& {\displaystyle }{\mathit{\sigma }}_i=\left\{\begin{array}{ll}\mathrm{1}& \mathrm{if}\mathrm{0}\le i\le p_{\mathrm{c}}\\ \exp \left[-{\mathit{\alpha }}_m{\left({\displaystyle \frac{i-p_{\mathrm{c}}}{p-p_{\mathrm{c}}}}\right)}^{p_m}\right]& \mathrm{if}{p}_{\mathrm{c}}\le i\le p,\end{array}\right.\end{array}$$

where p_c, p_m, and α_m represent the cutoff parameter, the order of the filter, and the non-dimensional decay strength, respectively. In this study, p_c was considered 0. We applied the filter ℱ to the solution vector q (in Eq. 4) at the final stage of the RK scheme with a time step Δt. Then, the decay timescale for the highest mode can be regarded as approximately equal to $\mathrm{\Delta }t/{\mathit{\alpha }}_m$. We set the order p_m and decay coefficient α_m such that the strength of the filter should ensure numerical stability while being as weak as possible. We checked the impact of the modal filter on the convergence rate in Sect. 3.1. In addition, the investigation into how much the modal filters can contaminate the eddy viscosity with the turbulent model and the energy spectra was performed in KT2023.

3 Verification of the dynamical core

We conducted several tests to verify our dynamical core. These tests are summarized in Table 4. For investigating the behavior of numerical convergence, the number of elements and the polynomial order were changed as detailed in Table 5. The dissipation mechanisms used in the numerical experiments are summarized in Table 6. When evaluating numerical errors for the deterministic experiments such as linear-advection, gravity-wave, mountain-wave, and baroclinic-wave tests, we used the following error norms as

$$\begin{array}{}\text{(42)}& \begin{array}{rl}& L_{\mathrm{1},\mathrm{error}}={\displaystyle \frac{{\sum }_{\mathrm{E}}{\int }_{{\mathrm{\Omega }}_{\mathrm{E}}}\left|\mathit{\psi }(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)-{\mathit{\psi }}_{\mathrm{ref}}(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)\right|\mathrm{d}\mathit{x}}{{\sum }_{\mathrm{E}}{\int }_{{\mathrm{\Omega }}_{\mathrm{E}}}\mathrm{d}\mathit{x}}},\\ & L_{\mathrm{2},\mathrm{error}}\\ & ={\left[{\displaystyle \frac{{\sum }_{\mathrm{E}}{\int }_{{\mathrm{\Omega }}_{\mathrm{E}}}{\left[\mathit{\psi }(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)-{\mathit{\psi }}_{\mathrm{ref}}(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)\right]}^{\mathrm{2}}\mathrm{d}\mathit{x}}{{\sum }_{\mathrm{E}}{\int }_{{\mathrm{\Omega }}_{\mathrm{E}}}\mathrm{d}\mathit{x}}}\right]}^{\mathrm{1}/\mathrm{2}},\\ & L_{\inf ,\mathrm{error}}=\max \left[\mathit{\psi }(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)-{\mathit{\psi }}_{\mathrm{ref}}(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)\right],\end{array}\end{array}$$

where $\mathit{\psi }(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)$ and ${\mathit{\psi }}_{\mathrm{ref}}(\mathit{\xi },\mathit{\eta },\mathit{\zeta },t)$ denote the numerical and reference solutions, respectively, and ∑_E represents the summation over all elements. Except for the linear-advection test case, the results obtained from a sufficiently high-resolution experiment were used as the reference solution because the exact solution is unknown. In such a case, the numerical solution was interpolated into the computational grid with the highest-resolution experiment when evaluating the error norms.

Table 4 Summary of numerical experiments for validating the proposed dynamical core.

For idealized climatological or turbulent-flow simulations, such as the Held–Suarez and global-LES tests, it is difficult to directly evaluate the numerical convergence using the error norms defined in Eq. (42). In the long-term integration, the chaotic behavior of the nonlinear systems can diverge from the numerical solutions. In the turbulent-flow simulations, a smaller-scale structure becomes more apparent as the grid spacing decreases until the spatial resolution reaches the physical-dissipation scale. Thus, for the test cases, we mainly investigated the impact of the polynomial order on the shortest wavelength at which the energy spectra began to separate from that in the reference solution.

3.1 Linear advection

To verify the spatial discretization with the cubed-sphere geometry, we conducted a two-dimensional linear-advection test for a scalar quantity q. The experimental setup is similar to test case 1 of Williamson et al. (1992). The longitudinal and latitudinal components of horizontal wind were prescribed by a solid-body rotation as

$$\begin{array}{}\text{(43)}& \begin{array}{rl}u& =u_{\mathrm{0}}(\cos \mathit{\varphi }\cos {\mathit{\varphi }}_{\mathrm{0}}+\sin \mathit{\varphi }\cos \mathit{\lambda }\sin {\mathit{\varphi }}_{\mathrm{0}}),\\ & v=-u_{\mathrm{0}}\sin \mathit{\lambda }\sin {\mathit{\varphi }}_{\mathrm{0}},\end{array}\end{array}$$

where λ and ϕ are the longitude and latitude coordinates, respectively; $u_{\mathrm{0}}=\mathrm{2}\mathit{\pi }a/\left(\mathrm{12}\right[\mathrm{d}\left]\right)$; and ϕ₀ denotes the angle between the axis of solid-body rotation and the north pole. We considered three values of ϕ₀=0, $\mathit{\pi }/\mathrm{4}$, and $\mathit{\pi }/\mathrm{2}$ rad to investigate the impact of singularity with four corners of each panel in the cubed sphere. Although a cosine-bell profile is often given as an initial profile of the advected field, we used a Gaussian profile to confirm the order of accuracy higher than 2. The profile is defined as

$$\begin{array}{}\text{(44)}& q(\mathit{\lambda },\mathit{\varphi })=\exp \left(-{\displaystyle \frac{d(\mathit{\lambda },\mathit{\varphi })}{D}}\right),\end{array}$$

where D is the characteristic horizontal scale and d is the circle distance between a position on the sphere (λ,ϕ) and the center position of Gaussian profile (λ_c,ϕ_c), which is calculated by

$$\begin{array}{}\text{(45)}& d(\mathit{\lambda },\mathit{\varphi })=a\arccos \left[\sin {\mathit{\varphi }}_{\mathrm{c}}\sin \mathit{\varphi }+\cos {\mathit{\varphi }}_{\mathrm{c}}\cos \mathit{\varphi }\cos (\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{c}})\right].\end{array}$$

In this experiment, $D=a/\mathrm{5}$ and $({\mathit{\lambda }}_{\mathrm{c}},{\mathit{\varphi }}_{\mathrm{c}})=(\mathrm{3}\mathit{\pi }/\mathrm{2},\mathrm{0})$.

To investigate a convergence rate, we changed the horizontal grid spacing Δ_h,eq from 313 to 39 km for p=1, 3, and 7 and from 417 to 52 km for p=11. As a temporal scheme, we adopted a fully explicit fourth-order RK scheme described in Sect. 2.4.2. To focus on the spatial errors, we set sufficiently small time steps such that $C_{\mathrm{r},\mathrm{adv}}=\mathrm{7.41}\times {\mathrm{10}}^{-\mathrm{2}}$. In this experiment, the modal filter was not required because the upwind numerical flux provided sufficient numerical stabilization.

Figure 1 Dependence of (a) L₁, (b) L₂, and (c) L_inf errors at t=12 d on horizontal resolution in a two-dimensional linear-advection test using p=1, 3, 7, and 11. The solid, dashed, and dotted colored lines represent the results for ϕ₀=0, $\mathit{\pi }/\mathrm{4}$, and $\mathit{\pi }/\mathrm{2}$, respectively. A dashed black line labeled “On” indicates the slope with nth-order accuracy.

Figure 1 shows the dependence of the L₁, L₂, and L_inf errors at 12 d on the horizontal resolution. As theoretically expected, we obtained p+1-order spatial accuracy for $p=\mathrm{1},\mathrm{3}$, and 7. For p=11 in the high spatial resolutions, the discretization error with the fourth-order temporal scheme, or the round-off error, degrades the convergence rate of 12th-order spatial accuracy. In the figure, the dashed lines represent the error norms in the case of ${\mathit{\varphi }}_{\mathrm{0}}=\mathit{\pi }/\mathrm{4}$ rad when the Gaussian profile passes over the singular points on cubed-sphere mesh. Their magnitudes were similar to those obtained for ϕ₀=0 and $\mathit{\pi }/\mathrm{2}$ rad, which are represented by solid and dashed lines, respectively. For p=1, the numerical errors were almost independent of the angle of the rotation axis. The errors for ${\mathit{\varphi }}_{\mathrm{0}}=\mathit{\pi }/\mathrm{4}$ rad can be smaller than those observed for ϕ₀=0 and $\mathit{\pi }/\mathrm{2}$ rad (e.g., p=3). The reason has not been confirmed, but we have found similar results in previous studies (Ullrich et al., 2010). In summary, when applying the DGM to the advection problem, we consider the influence of singularity with the cubed-sphere coordinate to be quite small.

Figure 2 The impact of the modal filters on the horizontal-resolution dependence of the L₁, L₂, and L_inf errors at t=12 d in a two-dimensional linear-advection test for the case ϕ_m=0 using p=1, 3, 7, and 11: (a) p_m=64, (b) p_m=32, (c) p_m=16, and (d) p_m=8. In each p_m, we changed α_m to 0 (without the filter), 10⁻³, 10⁻¹, and 10¹.

Because the modal filters are used in other test cases, we checked the impact of the modal filters on the numerical convergence. We conducted the linear-advection test in the case of ϕ₀=0, where the order and decay coefficient of the filter changed as p_m=64, 32, 16, and 8 and ${\mathit{\alpha }}_m={\mathrm{10}}^{-\mathrm{3}}$, 10⁻¹, and 10¹. Figure 2 shows the impact of the modal filters on the horizontal-resolution dependence of the error norms. The results indicate that the filters can degrade the original convergence rate and increase the numerical errors because the modal filter diminishes the high modes in the polynomial expansion. If the case of using high-order modal filters such as p_m≥32, the degradation of the convergence rate was in the range of 1–3 for p=7 and 11, even if we set sufficiently large values of α_m such that the highest mode was immediately decayed after one time step. For p=3, the degradation of the convergence rate appeared less obvious. However, it should be noted that the errors without the modal filter were much larger compared to p=7 and 11. Thus, for p=3, the effect of the increased error due to the filters may be more pronounced in the representation of the flow fields.

3.2 Internal gravity wave

To check wave propagation with pressure gradient and buoyancy terms, test cases of gravity waves are often utilized. For example, Tomita and Satoh (2004) also performed an internal-gravity-wave test. However, the basic state and initial perturbation produce vertically high modes, and nonlinear terms can develop small flow structures. This is inconvenient for investigating numerical convergence. On the other hand, the experimental setting based on Baldauf and Brdar (2013), which originally assumed a two-dimensional computational domain, can focus on a single mode. This study presents a global domain version of the gravity-wave test in Baldauf and Brdar (2013). The initial condition was a resting isothermal atmosphere of T₀=300 K, which corresponds to a constant Brunt–Väisälä frequency of $\sqrt{g^{\mathrm{2}}/\left(C_p{T}_{\mathrm{0}}\right)}\sim \mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{2}}$ s⁻¹. Furthermore, we added a small temperature perturbation with a Gaussian profile as

$$\begin{array}{}\text{(46)}& T^{\prime }=\mathrm{\Delta }T\exp \left(-{\displaystyle \frac{d}{D}}\right)\sin \left(n_{\mathrm{v}}\mathit{\pi }{\displaystyle \frac{z}{z_{\mathrm{T}}}}\right)\exp \left(-{\displaystyle \frac{g}{\mathrm{2}R{T}_{\mathrm{0}}}}z\right),\end{array}$$

where ΔT is the amplitude, D is the characteristic horizontal scale, n_v is the index with the vertical mode, and d was calculated using Eq. (45). In this experiment, we set z_T=10 km, ΔT=0.01 K, $D=a/\mathrm{5}$, n_v=1, and $({\mathit{\lambda }}_{\mathrm{c}},{\mathit{\varphi }}_{\mathrm{c}})=(\mathrm{0},\mathit{\pi })$. The Coriolis force and topography were not considered.

The horizontal and vertical grid spacing were changed from (313 km, 417 m) to (39 km, 52 m) using p=1, 3, and 7, whereas, for p=11, they were changed from (208 km, 417 m) to (102 km, 104 m). As the temporal scheme, we adopted an IMEX Runge–Kutta scheme with the third-order accuracy, as described in Sect. 2.4.1. For the HEVI scheme, we set the time step such that $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{1.34}\times {\mathrm{10}}^{-\mathrm{1}}$ for p=1, 3, and 7 and $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{1.26}\times {\mathrm{10}}^{-\mathrm{1}}$ for p=11. To investigate the impact of temporal errors, we also conducted additional experiments with smaller time steps for p=7 and p=11 where the above Courant number was reduced by factors of $\mathrm{1}/\mathrm{2}$ and $\mathrm{1}/\mathrm{4}$. In the absence of a modal filter, the self-convergence of numerical solutions was investigated. The reference solution was obtained from a high-resolution experiment where the horizontal and vertical grid spacing were (20 km, 26 m) with p=7 and $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{6.30}\times {\mathrm{10}}^{-\mathrm{2}}$.

Figure 3 The spatial distribution of potential temperature, zonal wind, and vertical wind at the equator after (a, b, c) t=0.5 d and (d, e, f) t=2 d obtained from a gravity-wave test case with $({\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}},{\mathrm{\Delta }}_{\mathrm{v}})=$(78 km, 104 m) using p=7.

To present the temporal evolution of gravity waves, Fig. 3 shows the spatial distribution of potential temperature, zonal wind, and vertical wind after t=0.5 and t=2 d. Based on this result, the horizontal phase speed was estimated to be approximately 58 m s⁻¹. This result corresponds well to the value based on a linear theory under the hydrostatic approximation of $N{z}_{\mathrm{T}}/\left(\mathit{\pi }{n}_{\mathrm{v}}\right)\sim \mathrm{57}$ m s⁻¹.

Figure 4 Dependence of the L₁, L₂, and L_inf errors on spatial resolution for the (a) density perturbation (ρ^′), (b) horizontal wind (u^ξ), (c) vertical wind (w), and (d) perturbation of potential-temperature-weighted density ((ρθ)^′) after t=2 d in a gravity-wave test case.

Figure 4 shows the dependence of error norms on spatial resolutions for the density perturbation ρ^′, horizontal wind u^ξ, vertical wind w, and perturbation of potential-temperature-weighted density (ρθ)^′. For relatively low-order p values, such as p=1 and p=3, nearly a p+1-order accuracy was observed for the four variables in sufficiently high spatial resolutions. However, due to the fast wave modes, temporal errors for the third-order HEVI scheme can dominate over the spatial errors in the cases of larger p values and higher resolutions. This behavior was evident for the error norms of all variables except for the horizontal wind. For the case of p=7 with $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{1.34}\times {\mathrm{10}}^{-\mathrm{1}}$, the convergence rate had approximately third-order slopes for ρ^′ and (ρθ)^′ and between second- and third-order slopes for w. As the time step decreased, the convergence rate for p=7 approached the p+1 order. Even for the cases of a small Courant number, due to an increase in round-off errors, the reduction in the error norms for ρ^′ and u^ξ stopped as the grid spacing decreased. For p=11, the problem of round-off errors was worse. Stagnating error reduction appeared in the spatial resolutions coarser than that in p=7, and the errors increased with the spatial resolution. Note that the influence of round-off errors might be overestimated because the amplitude of initial perturbation was significantly small and no modal filter was used in this experiment. Thus, the problem is considered to be not critical in practical simulations that include the modal filtering or turbulent schemes.

3.3 Mountain wave

Adopting the basic terrain-following coordinate introduced in Sect. 2.1 together with low-order schemes is known to produce large numerical errors with pressure gradient terms and to develop spurious flows (e.g., Zängl, 2012). However, such issues can be avoided using a high-order DGM. To check the numerical behavior of the basic terrain-following coordinate in a high-order DGM, we performed a mountain-wave test on a reduced planet radius based on Klemp et al. (2015) (referred to as KSP2015) and test case 2-1 in a Dynamical Core Model Intercomparison Project (DCMIP) test case document (Ullrich et al., 2012). Here, the planetary radius was set to $a/X_{\mathrm{r}}$, where X_r=166.7 is the scaling factor. In this experiment, the rotation was not considered. KSP2015 considered a topography profile in the form of

$$\begin{array}{}\text{(47)}& h^{\mathrm{KSP}\mathrm{2015}}(\mathit{\lambda },\mathit{\varphi })=h_{\mathrm{0}}\exp \left(-{\displaystyle \frac{{\stackrel{\mathrm{̃}}{d}}^{\mathrm{2}}}{d_{\mathrm{0}}^{\mathrm{2}}}}\right){\cos }^{\mathrm{2}}\left(\mathit{\pi }{\displaystyle \frac{\stackrel{\mathrm{̃}}{d}}{d_{\mathrm{1}}}}\right)\cos \mathit{\varphi },\end{array}$$

where $\stackrel{\mathrm{̃}}{d}=a/X_{\mathrm{r}}(\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{c}})$ (here, λ_c=π), d₀=5000 m, and d₁=4000 m. The maximum height of mountain h₀ was set to 25 m. In this topography profile, the mountain-wave structure along the equator is comparable to the results with a two-dimensional Schär type mountain (Schär et al., 2002). On the other hand, from the perspective of investigating the numerical convergence, it is undesirable for the zonal scale of topography to decrease with the latitudes and eventually become zero at the poles. To ensure that the minimum horizontal scale is sufficiently resolved in high-resolution simulations, we eliminated the undulation of the mountain at the high latitudes using a tapering function as

$$\begin{array}{}\text{(48)}& h(\mathit{\lambda },\mathit{\varphi })=h^{\mathrm{KSP}\mathrm{2015}}(\mathit{\lambda },\mathit{\varphi }){\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}+\tanh \left({\displaystyle \frac{\left|\mathit{\varphi }\right|-\mathit{\pi }/\mathrm{3}}{\mathrm{8}\mathit{\pi }/\mathrm{180}}}\right)\right].\end{array}$$

As the initial condition, we assumed a resting isothermal atmosphere of 300 K. KSP2015 considered an impulsive start where a zonal wind in solid-body rotation (u=U₀cos ϕ, where u₀=20 m s⁻¹) and the corresponding balanced state were initially given. However, such an impulsive start produces initial shocks with small spatial scales, which complicates the discussion on the numerical convergence. To mitigate the influence of an impulsive start, we gradually accelerated the wind using relaxation terms with the timescale of 60 s. Further details are provided in Appendix B1.

The horizontal and vertical grid spacing changed from (625 m, 500 m) to (156 m, 125 m) using p=3, 7, and 11. The model top was set to 30 km. As the temporal scheme, we used a fully explicit fourth-order RK scheme described in Sect. 2.4.2. For the HEVE scheme, we set the time steps such that $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{2.63}\times {\mathrm{10}}^{-\mathrm{1}}$. The reflection of waves at the model top was suppressed by introducing a sponge layer at z>15 km, where the vertical-element size linearly increased with the altitude. Moreover, a lateral sponge layer was placed to reduce the disruption of targeting mountain-wave structures by initial shocks propagating globally. For the details on the sponge layer, refer to Appendix B2. The reference solution was obtained from a high-resolution experiment where the horizontal and vertical grid spacing were (78 m, 62.5 m) with p=7. In this test case, to ensure numerical stability, we used a weak modal filter, which is summarized in Table 7.

Figure 5 The spatial distribution of vertical wind at the equator obtained from a mountain-wave test case with a Schär-like mountain: (a) numerical solution at t=2 h obtained from $({\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}},{\mathrm{\Delta }}_{\mathrm{v}})=$ (625 m, 500 m) using p=7 and (b) two-dimensional linear analytic solution on a flat plane (shown for comparison).

Figure 5a shows the spatial distribution of vertical wind after 2 h. For comparison, a linear analytic solution on a flat plane in the two-dimensional Cartesian coordinates is shown in Fig. 5b (the derivation can be found in Appendix A of KSP2015). Because the characteristic wavelength of the mountain scaled by the Scorer parameter is $d_{\mathrm{0}}N/U_{\mathrm{0}}\sim$ O(1), this setting corresponds to a nonhydrostatic regime of mountain waves. In such a regime, the waves with small-scale wavelengths are trapped near the surface, while large-scale waves propagate upward. The obtained wind pattern reproduced the results shown in Fig. 2a of KSP2015 well. On the other hand, the numerical solution and the linear analytic solution of a flat plane were slightly different. For example, the vertical wavelength of the large-scale waves in Fig. 5a was shorter compared to that in Fig. 5b. Based on the consideration using our regional dynamical core, the difference might be caused by the spherical experimental setup. Thus, we expect this difference to decrease as the planetary radius increases, while the spatial scale of the mountain remains unchanged.

Figure 6 Dependence of the L₁, L₂, and L_inf errors on spatial resolution for the (a) perturbation of density (ρ^′), (b) horizontal wind (u^ξ), (c) vertical wind (w), and (d) perturbation of potential-temperature-weighted density ((ρθ)^′) after 2 h in a mountain-wave test case.

Figure 6 shows the dependence of error norms on spatial resolution for the density perturbation ρ^′, horizontal wind u^ξ, vertical wind w, and perturbation of potential-temperature-weighted density (ρθ)^′. A comparison performed at the fixed DOF shows that the numerical errors decreased with the increase in polynomial order, although the convergence rate was smaller than p+1-order accuracy. For example, the slope of the L₂ error norm was approximately $\mathrm{3}/\mathrm{4}$ of that with p+1-order accuracy. Based on additional experiments with the corresponding two-dimensional setup, the sub-optimal convergence can be related to several factors such as the modal filter and the spatial discretization for Jacobian cofactors ($\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{13}}$ and $\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{23}}$). For further details, refer to Appendix B3. Because the modal filters shave off the high modes in the polynomial expansion, the convergence rate can be degraded. When requiring a convergence rate with a certain order of the accuracy, we need to increase the polynomial order according to the filter intensity.

Table 5 Summary of the polynomial order p, the number of elements, and the resulting equatorial resolution Δ_h,eq in the numerical experiments. For the number of elements, we denote the number of horizontally one-dimensional elements in a panel of the cubed sphere as N_e,h and the number of vertical elements as N_e,v.

Table 6 Summary of dissipation mechanism in the numerical experiments. In the table, ∘ means it is included, while − means it is not included. For the modal filtering in the mountain-wave, baroclinic-wave, and Held–Suarez tests, “scale-selective” means a high-order filter with p_m≥8 and “strong” means a large decay coefficient of α_m≥ O(1). The parameters of the filters are detailed in Tables 7, 8, and 9.

Table 7 Modal-filter orders and decay coefficients used in the mountain-wave test.

Table 8 Modal-filter orders and decay coefficients used in the baroclinic-instability test. Because the vertical resolution was fixed when increasing the horizontal resolution, the decay coefficient for the vertical filter α_m,v was reduced in proportion to the time step.

3.4 Baroclinic instability

Baroclinic instability is a typical phenomenon in the mid-latitudes. It includes small-scale structures such as front and filament formations. We conducted an idealized numerical experiment based on Jablonowski and Williamson (2006) (referred to as JW2006). In the aspect of numerical method, dynamical cores must accurately represent the wave growth process. In addition, it is necessary to treat the developing small-scale flow structures while ensuring numerical stability. For the experimental setup, the initial zonally symmetric fields were expressed using the analytic expressions of a steady-state solution of the adiabatic inviscid primitive equations. To trigger baroclinic instability, a perturbation with a Gaussian profile was added to the zonal wind in the northern hemisphere. For further details on parameter values, refer to Sect. 2a of JW2006.

We investigated the dependence of numerical solutions on the horizontal resolution as performed in JW2006. The horizontal grid spacing Δ_h,eq changed from 250 to 32 km with a fixed total vertical DOF of 24 for p=3 and 7, whereas it changed from 208 to 52 km with a fixed total vertical DOF of 36 for p=11. We used a stretched vertical grid based on Eq. (102) in Ullrich et al. (2012). The stretching parameter was set such that the vertical grid spacing Δ_v near the surface took the values of 305, 523, and 426 m for $p=\mathrm{3},\mathrm{7}$, and 11, respectively. The stretching is further detailed in Appendix C. For the third-order HEVI scheme, we set the time steps such that $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{1.68}\times {\mathrm{10}}^{-\mathrm{1}}$ for p=3 and 7 and $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{1.26}\times {\mathrm{10}}^{-\mathrm{1}}$ for p=11. Furthermore, the modal filter was utilized to maintain the numerical stability. Its parameters are summarized in Table 8. When calculating the L₂ error in surface pressure, we used the results obtained from the corresponding highest-resolution experiment for each p as the reference solution to directly discuss the behavior of numerical convergence associated with the horizontal spatial or temporal accuracy.

Figure 7 Spatial distribution of the surface pressure and temperature at 850 hPa after $t=\mathrm{4},\mathrm{6},\mathrm{8}$, and 10 d in a baroclinic-instability test. We present the results obtained from the experiment using p=7 with ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{63}$ km and a vertical DOF value of 24.

Figure 8 Horizontal-resolution dependence of (a) surface pressure and (b) temperature at 850 hPa after t=9 d in a baroclinic-instability test. We present the results obtained from the experiments using p=7 with a vertical DOF value of 24.

Figure 7 shows the temporal evolution of the baroclinic wave for the case of ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{63}$ km using p=7. The obtained horizontal distributions of surface pressure and temperature at 850 hPa were similar to those reported in the previous studies. For example, see Fig. 5 in JW2006, which was obtained from the finite-volume (FV) dynamical core (Lin and Rood, 1996, 1997). The wave grew very slowly for 4 d. After that, the highs and lows deepened significantly and the wave began to break on the eighth day. Figure 8 shows dependence of the surface pressure and temperature at 850 hPa (after 9 d) on the horizontal spatial resolution for p=7. The same figure obtained from the FV dynamical core can be seen in Fig. 6 of JW2006. Our dynamical core provided reasonably accurate numerical solutions for experiments performed at high spatial resolution. These solutions were comparable to the reference solutions reported in the previous studies. In addition, the effective resolution was apparently higher than that of the low-order global dynamical core. For example, in the marginally resolved simulation setting, ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}\sim \mathrm{250}$ km, the amplitude and phase errors were small.

Figure 9 Dependence of the L₂ error norm for surface pressure on horizontal spatial resolution in a baroclinic-instability test using (a) p=3, (b) p=7, and (c) p=11. Note that the reference solution for each p is the result of the corresponding highest-resolution experiment.

For a quantitative evaluation of the horizontal-resolution dependence, Fig. 9 shows the temporal evolution of the L₂ error norm of the surface pressure for $p=\mathrm{3},\mathrm{7}$, and 11. In the figure, the gray shading represents an uncertainty range of reference solutions estimated by various dynamical cores in JW2006. In our evaluation strategy, we successfully captured the numerical convergence with horizontal discretization and temporal errors. This is because the vertical spatial errors have similar values among different horizontal resolution cases with the same p, and these errors virtually cancel out when the L₂ error is evaluated. Until about 6 d (except the initial adjustment stage), the L₂ error decreased with the horizontal resolution. The magnitude was significantly small compared to that reported in previous studies. For example, in the horizontal grid spacing of 50 km (0.5°), the L₂ error was $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{2}}$ for the FV dynamical core and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ for MCore (Ullrich and Jablonowski, 2012b). After 6 d when the baroclinic wave started to develop significantly, the L₂ errors rapidly grew and the difference between horizontal resolutions decreased. For p=3, the feature of numerical convergence at the mature stage was similar to that obtained from MCore. In summary, the L₂ errors for ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}<\mathrm{250}$ km are within the uncertainty range suggested by JW2006. Thus, we consider the numerical solutions obtained from the proposed model to be reasonable.

3.5 Held–Suarez test

As a long-term idealized benchmark for real climate simulations, we conducted the Held–Suarez test (Held and Suarez, 1994), which used a prescribed forcing that mimics complex physics parameterizations. The boundary layer friction was represented in the form of Rayleigh damping. The diabatic heating–cooling effect was represented using a Newtonian relaxation term to a prescribed temperature in radiative equilibrium T_e. For further details on these terms, see p. 1826 in Held and Suarez (1994). In this study, a resting atmospheric field in hydrostatic balance with T_e was given as the initial condition.

To investigate the spatial-resolution dependence, the horizontal grid spacing Δ_h,eq and vertical DOF changed from (208 km, 32) to (52 km, 128) for p=3 and 7 and from (208 km, 36) to (52 km, 144) for p=11. The vertical grid spacing was stretched using the strategy in Appendix C. For example, when the vertical DOF value is 32, the vertical grid spacing near the surface becomes approximately 350 m. In the cases of ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{208},\mathrm{104}$, and 52 km, the temporal integration was performed for 1200 d. Data from the first 200 d were discarded during the statistical analysis. For high-resolution cases of ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{52}$ and 26 km, the results after the spin-up period with coarser resolutions were used as the initial data, and a temporal integration was conducted for 1000 d. As the temporal scheme, we adopted the third-order HEVI scheme with the acoustic Courant number of $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{1.26}\times {\mathrm{10}}^{-\mathrm{1}}$ for p=3 and 7 and $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{7.56}\times {\mathrm{10}}^{-\mathrm{2}}$ for p=11. Moreover, we used the modal filters summarized in Table 9. Note that the large decay coefficients were set to stabilize long temporal integrations with nonlinear flow processes. The reference solution was obtained from a high-resolution experiment where Δ_h,eq and vertical DOF were (26 km, 256) with p=7.

Table 9 Modal-filter orders and decay coefficients used in the Held–Suarez test after 250 d of spin-up experiments. Note that we temporarily increased α_m,h to 6.0×10⁰ for 20 d during the 1000 d integration for the case of ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{26}$ km using p=7.

Figure 10 Zonally and temporally averaged atmospheric fields in a statistical equilibrium state: (a) zonal wind, (b) temperature, (c) meridional eddy flux of zonal momentum, (d) meridional eddy flux of temperature, (e) eddy kinetic energy, and (f) eddy temperature variance obtained from a Held–Suarez test with ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{208}$ km using p=7. As is typically done in previous studies, when taking the zonal and temporal average, we used the 1000 d of data after the spin-up calculation.

Figure 10 shows the zonally and temporally averaged atmospheric fields in a statistical equilibrium state for ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{208}$ km using p=7. The obtained pattern and strength of general circulations were similar to the results obtained using nearly the same horizontal spatial resolution used in previous studies (e.g., Wan et al., 2008). For a single westerly jet in each hemisphere, a maximum velocity of 32 m s⁻¹ was obtained at P=250 hPa. Easterlies existed in the equatorial and polar lower atmosphere and near the model top at low latitude. As shown in Fig. 10c–f, the baroclinic-wave activity in the proposed DG model was similar to that reported in previous studies. At P=250 hPa in the mid-latitudes, the magnitude of eddy momentum flux reached approximately 70 m² s⁻². The maximum of the poleward eddy heat flux was located at P=850 hPa in the mid-latitudes, and its value reached approximately 22 K m s⁻¹. The eddy kinetic energy and temperature variance reached maximum values of approximately 430 m² s⁻² at P=250 hPa and 45 K² at P=800 hPa in the mid-latitudes, respectively.

Figure 11 Dependence of absolute peak values on the spatial resolution: (a) eddy temperature variance, (b) eddy heat flux, (c) eddy kinetic energy, and (d) eddy momentum flux. The time-averaging period is the same as that in Fig. 11. In each panel, the difference between the northern and southern hemispheres is represented using the error bars. The colored boxes labeled by T04, W08, JW12, and W13 denote the corresponding peak values indicated by the results reported in Tomita and Satoh (2004), Wan et al. (2008), Ullrich and Jablonowski (2012b), and Wan et al. (2013), respectively. Because the peak values were estimated from the contour figures, note that the uncertainty is large, and its range is roughly represented by the box height.

As discussed in a previous study (Wan et al., 2008), these eddy quantities such as the eddy kinetic energy and temperature variance are sensitive to the spatial resolution. As shown in Fig. 11a and c, the absolute peak values increased with the spatial resolution and began to converge when the horizontal grid spacing was less than 50 km. The convergence of peak values with p=7 and 11 was faster than that in the case of p=3. For comparison, the corresponding peak values indicated from previous studies are shown by the colored boxes in the figure. The obtained trend of spatial-resolution dependence for p=3 was similar to that reported by studies using conventional low-order grid-point methods (Tomita and Satoh, 2004; Wan et al., 2013). On the other hand, the peak values from p=7 and 11 at the horizontal grid spacing of approximately 200 km were similar to the results obtained using the horizontal spectral method (Wan et al., 2008).

Figure 12 Energy spectra of horizontal wind in a statistical equilibrium state at (a) P=850 hPa and (b) P=250 hPa in a Held–Suarez test. As explained in the legend, the difference between spatial resolutions is represented by line types and the line color indicates the polynomial order. The results from the reference experiment are shown by solid black lines. Lower panels represent the compensated spectra, which is proportional to E(n)n³. The temporal average was calculated over 1000 d after the spin-up period. In the highest-resolution case (${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{26}$ km), it was performed over 300 d. VDOF: vertical degree of freedom.

Figure 12 shows kinetic energy spectra of horizontal winds at p=850 and P=250 hPa. As reported in previous studies (e.g., Malardel and Wedi, 2016; Tolstykh et al., 2017), the obtained spectra had the n⁻³ slope at the spherical-harmonic degrees between 10–100. The steeper slope compared to −3 reflects the influence of a numerical dissipation mechanism with the upwind numerical flux and the modal filter. For the cases of p=7 and 11, the obtained energy spectra followed that for the reference experiment at the wavelengths longer than eight grid lengths well. In the spatial-resolution dependence of peak values with the eddy quantities shown in Fig. 11, there was no significant difference between p=7 and p=11, whereas the improvement of effective resolution by the higher polynomial order of p=11 can be observed in the energy spectra as the grid spacing decreases. For p=3, the energy spectra overlapped with that of the reference experiment at a wavelength range longer than 10–20 grid lengths. Furthermore, for ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{208}$ km with p=3, the entire spectra were smaller than the reference solution. Thus, using strong modal filters to ensure numerical stability for long-term integration has a significant effect on the spectra and effective resolutions when using p≤3. These results indicate that there is room for improving our treatment of the nonlinear terms to weaken the modal filters.

3.6 Planetary boundary layer turbulence experiment on a small planet

As a first step toward future global LES with O(10 m) grid spacing, we performed a global extension of the LES experiment of idealized planetary boundary layer turbulence in Nishizawa et al. (2015), KT2021, and KT2023. Currently, it is not feasible to conduct a global LES for a planetary size of Earth using a uniform spatial resolution of O(10 m). To save the required computational resources, the planetary radius was set to 3.4 km. Although this value is significantly different from that in realistic planets such as Earth and the effect of spherical geometry may affect the convection structure, we consider this test to be useful to verify the turbulent model described in Sect. 2.2. We focused on the case of applying shallow-atmosphere approximation because we expected the results to be comparable to those reported in our previous studies. This approximation is obviously unsuitable for discussing the physical aspect in this experimental setting. For the case without approximation, refer to Appendix D. The experimental setup is as follows. The altitude of model top was set to 3 km. There were no rotation and topography. Initially, we set a stable stratification with a vertical gradient of potential temperature of 4 K km⁻¹ and added random perturbations with an amplitude of 1 K. Because it is difficult to consider a uniform wind in the global situation, there was no initial motion in contrast to that in our previous studies. To drive thermal convections, a constant heat flux of 200 W m⁻² was imposed at the surface. To focus on the turbulent model, radiation and moist processes were not considered. In the turbulent model, we set the filter length to double that of the local grid spacing, which follows our previous studies. A reflection of waves at the model top was prevented using a sponge layer, where the vertical wind was decayed by the Rayleigh damping. The e-folding time varied as the half-cosine function from 0 at z=2 km to 10 s at the model top.

To check the impact of polynomial order on the energy spectra as in KT2023, we changed p to 3, 4, and 7 while setting the horizontal and vertical grid spacing to be approximately 10 m. Numerical stability was ensured using a modal filter with parameters $p_{m,\mathrm{h}}=p_{m,\mathrm{v}}=\mathrm{32}$ and ${\mathit{\alpha }}_{m,\mathrm{h}}={\mathit{\alpha }}_{m,\mathrm{v}}={\mathrm{10}}^{-\mathrm{3}}$. As the temporal scheme, a fully explicit fourth-order RK scheme described in Sect. 2.4.2 was used for the inviscid terms, whereas the forward Euler scheme was adopted for the SGS terms. We set the time step such that $C_{\mathrm{r},{\mathrm{c}}_{\mathrm{s}}}=\mathrm{4.38}\times {\mathrm{10}}^{-\mathrm{1}}$. The integration time was 4 h for the case of p=7. To reduce the computational cost, the output at 3 h was used as the initial condition of the other experiments for which the integration time was 1 h.

Figure 13 The horizontal distribution of vertical wind at z=500 m after t=4 h when shallow-atmosphere approximation is applied in the LES of an idealized planetary boundary layer turbulence for the case of ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}=\mathrm{10}$ m using p=7: (a) northern hemisphere (NH), (b) southern hemisphere (SH), and (c) their corresponding cross-sections along the equator.

Figure 14 The vertical structure of the PBL temporally averaged during the last 30 min: (a) potential temperature, (b) resolved eddy heat flux plus SGS heat flux, (c) variance of vertical wind, and (d) skewness of vertical wind for $p=\mathrm{3},\mathrm{4}$, and 7. The gray shading represents the results obtained from KT2023 using the plane model.

Figure 13 shows the horizontal distribution for vertical wind at z=500 m and cross-section along the equator after t=4 h in the case of p=7. The convective cells had polygonal structures with a horizontal scale of approximately 2–3 km. The height of the PBL reached between 1–1.5 km. To present the vertical structure of the PBL, Fig. 14 shows the vertical distribution of potential temperature, turbulent transport of heat and momentum, and skewness of vertical wind for p=7. In these panels, the gray shading represents the results obtained from KT2023 using the plane regional model. The results obtained in this study were similar to those reported in KT2023.

Figure 15 (a) Density-weighted energy spectra E(n) of three-dimensional wind at the height of 500 m for $p=\mathrm{3},\mathrm{4}$, and 7. The dash-dotted gray line represents aE(n), where $a=\mathrm{8.0}\times {\mathrm{10}}^{\mathrm{1}}$. (b) Spectra normalized by the result of p=7. (c) Partially expanded view of energy spectra in the short-wavelength range.

Figure 15 shows the kinetic energy spectra of three-dimensional wind at z=500 m, which were temporally averaged during the last 30 min. The features of the obtained energy spectra were similar to those reported in KT2023. At wavelengths longer than eight grid lengths, the slope of spectra was approximately $-\mathrm{5}/\mathrm{3}$. On the other hand, the slope of spectra at the shorter wavelength range deepened due to the turbulent model, numerical dissipation with the upwind numerical flux, and modal filter. KT2023 indicated that p>3 is required for the effect of numerical diffusion term to be sufficiently small compared to that of the SGS eddy viscosity term at the wavelength longer than eight grid lengths. This is true for global LES as shown in Fig. 15b.

4 Conclusions

For conducting future high-resolution atmospheric simulations precisely, our previous study (KT2021) indicated that conventional low-order discretization methods used in the state-of-the-art global nonhydrostatic dynamical cores have a problem with numerical errors because it is possible to contaminate the effect of physical parameterization schemes. To overcome this issue, we developed a global nonhydrostatic atmospheric dynamical core of dry atmosphere using the discontinuous Galerkin method (DGM) as the spatial discretization because the DGM has several advantages over grid-point methods, including the simple strategy of high-order discretization and the high count of floating-point operations per second (FLOPS) in recent parallel supercomputers. Furthermore, considering global LES, we formulated a Smagorinsky–Lilly type turbulent model in the cubed-sphere coordinates and discretized it in the DGM framework. To verify the proposed global dynamical core, several numerical experiments, from the linear-advection test to the Held–Suarez test, were conducted. To demonstrate the high-order numerical convergence, the experimental setup of existing test cases were slightly modified. In addition, an idealized test case was proposed to check the behavior of global dynamical cores including the turbulent model. Through the numerical experiments with various polynomial (p) orders and spatial resolutions, we discussed the impact of high-order spatial discretization on the quality of numerical solutions in the atmospheric-flow simulations.

For the deterministic test cases, such as the linear-advection and gravity-wave tests, p+1-order spatial accuracy was confirmed until the temporal discretization and round-off errors became significant compared to the spatial errors. In the gravity-wave test, it was observed that the temporal errors with the third-order HEVI scheme can contaminate the convergence rate of high-order spatial discretization even when using the horizontal acoustic Courant number of O(0.1). To investigate the numerical performance of the terrain-following coordinate with the DGM, we conducted a mountain-wave test case based on that used in Klemp et al. (2015). However, we made some modifications to investigate high-order numerical convergence. When comparing the results for a fixed DOF, the advantage of a large polynomial order was apparent in terms of the fast numerical convergence, although the resultant convergence rate was slightly smaller than the optimal order associated with the spatial discretization. The results of the baroclinic-instability test showed that, when p≥3 and ${\mathrm{\Delta }}_{\mathrm{h},\mathrm{eq}}<\mathrm{240}$ km, the obtained L₂ error norms of surface pressure entered the uncertainty range indicated by the previous studies. We confirmed the rapid numerical convergence over the second-order accuracy until the mature stage was reached. Subsequently, the sharp gradient with the front structure developed and the waves began to break, which made it difficult to identify the numerical convergence with the high-order schemes.

For test cases in which small-scale turbulent structures developed, such as the Held–Suarez test and the LES experiment of PBL turbulence, we mainly focused on the energy spectra in terms of the effective resolution. In the Held–Suarez test, where the turbulence model was not used, the extent of the dissipation effect with the numerical flux and modal filters was clearly visible. Based on the comparison with the energy spectra for the reference experiment, we confirmed that the effective resolution was improved as the polynomial order increased. When we used high-order modal filters with large decay coefficients to ensure numerical stability during long temporal integration, the effective resolution was estimated to be between 10–20 grid lengths for p=3 and eight grid lengths for p=7 and 11. To enhance the effective resolutions by weakening the modal filters, we consider the entropy-stable DGM adopted in Souza et al. (2023) to be a promising method, although this topic is beyond the scope of current study. To check the behavior of turbulent model included in global dynamical cores, we proposed an idealized test case of global LES considering a small planetary radius, which is an extension of the experimental setup used in KT2021 and KT2023 with regional plane models. In our numerical experiments with shallow-atmosphere approximation, the convective cell pattern and vertical structures of the PBL reproduced the results of the regional plane model well. We confirmed that the obtained energy spectra obeyed the Kolmogorov spectra of turbulence at the wavelength range longer than eight grid lengths when p>3 was used together with the Rusanov numerical flux and a scale-selective modal filter. This result was consistent with the numerical criteria discussed in KT2023.

This study demonstrated the applicability of a high-order DGM to global atmospheric dynamical cores via a series of numerical experiments. However, several tasks required to conduct realistic atmospheric simulations were not performed. To treat the effect of topography in LES, we must also consider the vertical coordinate transformation in the SGS terms of a turbulent model. Such a formulation can be achieved using the chain rule, as performed in the differential terms with inviscid fluxes. A related issue is the treatment of topography with steep slopes in high-order strategies. To investigate whether the DGM-based dynamical core is robust for realistic topography, a Held–Suarez experiment with realistic topography may be appropriate. Such work is expected to yield deep insights into the impact of effective resolutions of topography on large-scale flows when a high-order DGM is used. Furthermore, a severe time step restriction for explicit temporal schemes is one of the unresolved issues in a high-order DGM. We expect that the computational overhead would be ignored in several cases. A coarser spatial resolution can be used due to the high-order numerical convergence or when the small communication cost in the DGM is taken advantage of. However, to accelerate DG dynamical cores in all situations, developing sophisticated temporal treatments is an important future step. Finally, it is indispensable to perform a coupling between the physics (such as moist and radiation processes) and DGM-based dynamics. Recent studies begin to discuss the potential difficulties with the element-based methods. For example, in the context of SEM, Herrington et al. (2019) indicated that a straightforward evaluation for physics tendencies at irregular nodes within the element causes a grid imprinting along the element boundaries. To solve this problem, they introduced a physics grid with a quasi-equal volume coarser than the node intervals with the dynamics when calculating the physics tendencies. While taking care of the advantages associated with the effective resolutions of high-order dynamical cores, we will explore how to treat the coupling of physics and dynamics in the DGM framework. For actual operational runs including physical processes and other factors, such as realistic topography, the degree of numerical filters depends on situations and cannot be generalized now. It is an important issue to produce a kind of criterion for numerical stabilization.

Appendix A: Additional linear-advection test

Figure A1 Dependence of (a) L₁, (b) L₂, and (c) L_inf errors at t=12 d on the horizontal resolution in test case 4 in NL2010 with the Gaussian hills using p=1, 3, 7, and 11. The dashed green line represents the case of applying a modal filter (MF) for p=11.

In Sect. 3.1, by conducting a linear-advection test with a smooth profile in a solid-body rotation flow, we tested the spatial discretization with the cubed-sphere geometry and checked the convergence rate with the high-order DGM. In this section, we performed a linear-advection test using the Gaussian hills and slotted cylinder profiles in a deformation flow, which is test case 4 presented in Nair and Lauritzen (2010) (hereinafter, referred to as NL2010). The experimental setup with the spatial resolution, polynomial order, and time step was similar to that described in Sect. 3.1. To compare the errors reported in Guba et al. (2014) (hereinafter, referred to as G2014), we normalized the errors following Appendix C of NL2010.

Figure A1 shows the dependence of error norms on the horizontal resolution when the Gaussian hills were given as the initial condition of the tracer. Because of the infinitely smooth profile, we obtained p+1-order accuracy for $p=\mathrm{1},\mathrm{3},\mathrm{7}$, and 11. The behavior of numerical convergence and the magnitude of errors were comparable to those in G2014 (see the values of l₁, l₂, and l_∞ for “Gauss.” with the hyperviscosity (but without the limiter) in Tables 1–4 and Fig. 4 of G2014). In our results, the L_inf error was larger than the unity for p=11 in the coarsest spatial resolution. This reflects a numerical instability with the aliasing errors near the static stagnation point of the deformation flow. It occurred when the static stagnation point was located at the element boundaries and numerical dissipation was not sufficient. By introducing a very weak modal filter with ${\mathit{\alpha }}_{m,\mathrm{h}}=\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{2}}$ and $p_{m,\mathrm{h}}=\mathrm{64}$, we can control the numerical instability as shown by the dashed green line in Fig. A1.

Figure A2 Dependence of (a) L₁, (b) L₂, and (c) L_inf errors at t=12 d on the horizontal resolution in test case 4 in NL2010 with the slotted cylinders using p=1, 3, 7, and 11. The dashed lines represent the cases of applying a modal filter for p=7 and 11.

Figure A2 shows the dependence of error norms on the horizontal resolution when the slotted cylinders were given as the initial tracer profile. Because of the discontinuous field, we cannot expect a convergence rate higher than the first-order accuracy. Even when using the high polynomial orders of p>3, we obtained at most the first order for the L₁ error norm. For the L_inf error, the convergence rate was near the zeroth order. The behavior of slow numerical convergence and the magnitude of numerical errors were similar to those reported in G2014 (see the error norms for “Cyl.” with the hyperviscosity (but without the limiter) in Tables 1–4 of G2014). As seen in the case of Gaussian hills, due to the numerical instability near the static stagnation point, we observed very large L_inf error values for the cases of p=7 and 11 with the coarse spatial resolution. The modal filter used for the case of Gaussian hills can suppress the numerical instability as shown by the dashed lines in Fig. A2.

Appendix B: Additional information for the mountain-wave test

In this section, we detail the spin-up strategy and sponge layer, which were used in the mountain-wave test described in Sect. 3.3. In addition, we consider reasons why the obtained convergence rate was slightly less than the optimal order accuracy.

B1 Spin-up strategy

To mitigate the influence of impulsive start on numerical solutions, we gradually accelerated the wind as performed in previous studies with a regional experimental setup (e.g., Durran, 1986; Sachsperger et al., 2016). The initial condition was a resting isothermal atmosphere and was represented as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{u}^{\mathit{\xi }}=\mathrm{0},u^{\mathit{\zeta }}=\mathrm{0},u^{\mathit{\eta }}=\mathrm{0},\\ \text{(B1)}& {\displaystyle }& {\displaystyle }P=P_{\mathrm{0}}\exp \left(-{\displaystyle \frac{gz}{R{T}_{\mathrm{0}}}}\right),{\displaystyle }& {\displaystyle }\mathit{\rho }={\displaystyle \frac{P_{\mathrm{0}}}{R{T}_{\mathrm{0}}}}\exp \left(-{\displaystyle \frac{gz}{R{T}_{\mathrm{0}}}}\right),\end{array}$$

where P₀=10⁵ Pa and T₀=300 K. To accelerate a zonal wind, we added the relaxation terms in the right-hand side of governing equations as

$$\begin{array}{ll}{\displaystyle \frac{\partial \sqrt{G}{\mathit{\rho }}^{\prime }}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{f}}\sqrt{G}{\mathit{\rho }}^{\prime },\\ {\displaystyle \frac{\partial \sqrt{G}\mathit{\rho }{u}^{\mathit{\xi }}}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{f}}\sqrt{G}\left(\mathit{\rho }{u}^{\mathit{\xi }}-\mathit{\rho }{U}^{\mathit{\xi }}\right),\\ \text{(B2)}& {\displaystyle \frac{\partial \sqrt{G}\mathit{\rho }{u}^{\mathit{\eta }}}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{f}}\sqrt{G}\left(\mathit{\rho }{u}^{\mathit{\eta }}-\mathit{\rho }{U}^{\mathit{\eta }}\right),{\displaystyle \frac{\partial \sqrt{G}\mathit{\rho }{u}^{\mathit{\zeta }}}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{f}}\sqrt{G}\left(\mathit{\rho }{u}^{\mathit{\zeta }}-\mathit{\rho }{U}^{\mathit{\zeta }}\right),\\ {\displaystyle \frac{\partial \sqrt{G}(\mathit{\rho }\mathit{\theta }{)}^{\prime }}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{f}}\sqrt{G}(\mathit{\rho }\mathit{\theta }{)}^{\prime },\end{array}$$

where $(U^{\mathit{\xi }},U^{\mathit{\eta }},U^{\mathit{\zeta }})$ are the vector components of prescribed wind and α_f is a time-dependent coefficient with Rayleigh forcing terms, which is provided in this subsection. Note that we set the hydrostatic balance part of pressure and density to

$$\begin{array}{}\text{(B3)}& P_{\mathrm{hyd}}=P_{\mathrm{0}}\exp \left(-{\displaystyle \frac{u_{\mathrm{eq}}}{\mathrm{2}R{T}_{\mathrm{0}}}}{\sin }^{\mathrm{2}}\mathit{\varphi }-{\displaystyle \frac{gz}{R{T}_{\mathrm{0}}}}\right),{\mathit{\rho }}_{\mathrm{hyd}}={\displaystyle \frac{P_{\mathrm{hyd}}}{R{T}_{\mathrm{0}}}},\end{array}$$

which satisfies a dynamically balanced state associated with a zonal wind in solid rotation (u_eqcos ϕ). Then, the perturbation at the initial time is given by $P^{\prime }=P-P_{\mathrm{hyd}}$, ${\mathit{\rho }}^{\prime }=\mathit{\rho }-{\mathit{\rho }}_{\mathrm{hyd}}$.

Figure B1 After 2 h in a mountain-wave test case with a global model and spatial distribution of the (a) L₂ error and (b) local convergence rate for the vertical wind at the equator. For the L₂ error, we show the result obtained from the experiment where the horizontal and vertical grid spacing (Δ_h,Δ_v) are set to 156 and 125 m, respectively, using p=3. The cell color in the figure corresponds to the average values within the finite element. When evaluating the local convergence rate, we used the element average of the L₂ error obtained from two experiments: a coarse-resolution experiment (Δ_h=625 m, Δ_v=500 m) and the highest-resolution experiment for p=3 (Δ_h=156 m, Δ_v=125 m). To see the large-scale structure of the local convergence rate, the moving average was taken across the five elements horizontally. In each figure, the white lines represent the vertical wind in the highest-resolution experiment for p=3.

As the horizontal component of prescribed wind, we considered a zonal wind in solid-body rotation, where u_eq=20 m s⁻¹. The corresponding (U^ξ,U^η) can be calculated by considering the coordinate conversion between the cubed-sphere and geographic coordinates. To improve the inconsistency with the no-flux boundary condition at the surface, the vertical component was added in the form of

$$\begin{array}{}\text{(B4)}& {\displaystyle }{U}^{\mathit{\zeta }}=-{\sqrt{G}}_{\mathrm{v}}(G_{\mathrm{v}}^{\mathrm{13}}{U}^{\mathit{\xi }}+G_{\mathrm{v}}^{\mathrm{23}}{U}^{\mathit{\eta }})\exp \left(-{\displaystyle \frac{\mathit{\zeta }}{H_{\mathrm{f}}}}\right),\end{array}$$

where H_f was set to 2 km in this study. This modification also reduces the influence of initial shock. On the other hand, the coefficient with the forcing terms was given as ${\mathit{\alpha }}_{\mathrm{f}}\left(t\right)=w\left(t\right){\mathit{\tau }}_{\mathrm{f}}^{-\mathrm{1}}$, where

$$\begin{array}{lll}{\displaystyle }w& {\displaystyle }=\mathrm{1}& {\displaystyle }\mathrm{for}\mathrm{0}\le t\le t_{\mathrm{1}},\\ \text{(B5)}& {\displaystyle }w& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}-\cos \left(\mathit{\pi }{\displaystyle \frac{t-t_{\mathrm{1}}}{t_{\mathrm{2}}-t_{\mathrm{1}}}}\right)\right]& {\displaystyle }\mathrm{for}{t}_{\mathrm{1}}\le t\le t_{\mathrm{2}},{\displaystyle }w& {\displaystyle }=\mathrm{0}& {\displaystyle }\mathrm{for}t\ge t_{\mathrm{2}},\end{array}$$

with τ_f as the forcing timescale. In this study, these parameters were set to τ_f=60 s, t₁=200 s, and t₂=1800 s.

B2 Sponge layer

To suppress a reflection of waves at the model top, we introduced a sponge layer at the upper computational domain. In addition, to reduce the disruption of our targeting structure of mountain waves due to the global propagation of initial shocks, a lateral sponge layer was placed. As in Eq. (B2), linear damping terms were added to the governing equations as follows:

$$\begin{array}{ll}{\displaystyle \frac{\partial \sqrt{G}{\mathit{\rho }}^{\prime }}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{s}}\sqrt{G}{\mathit{\rho }}^{\prime },\\ {\displaystyle \frac{\partial \sqrt{G}\mathit{\rho }{u}^{\mathit{\xi }}}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{s}}\sqrt{G}\left(\mathit{\rho }{u}^{\mathit{\xi }}-U^{\mathit{\xi }}\right),\\ \text{(B6)}& {\displaystyle \frac{\partial \sqrt{G}\mathit{\rho }{u}^{\mathit{\eta }}}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{s}}\sqrt{G}\left(\mathit{\rho }{u}^{\mathit{\eta }}-\mathit{\rho }{U}^{\mathit{\eta }}\right),{\displaystyle \frac{\partial \sqrt{G}\mathit{\rho }{u}^{\mathit{\zeta }}}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{s}}\sqrt{G}\left(\mathit{\rho }{u}^{\mathit{\zeta }}-\mathit{\rho }{U}^{\mathit{\zeta }}\right),\\ {\displaystyle \frac{\partial \sqrt{G}(\mathit{\rho }\mathit{\theta }{)}^{\prime }}{\partial t}}=& {\displaystyle }\mathrm{\dots }-{\mathit{\alpha }}_{\mathrm{s}}\sqrt{G}(\mathit{\rho }\mathit{\theta }{)}^{\prime }.\end{array}$$

The decay coefficient was given as ${\mathit{\alpha }}_{\mathrm{s}}=(\mathrm{1}-w(t\left)\right)({\mathit{\alpha }}_{\mathrm{s},\mathrm{h}}+{\mathit{\alpha }}_{\mathrm{s},\mathrm{v}})$, where α_s,h and α_s,v are the coefficients for lateral and upper sponge layers, respectively. To avoid the sponge layer interfering with the initial forcing in Eq. (B2), as the initial forcing weakens, we gradually activated the sponge layer using the coefficient (1−w(t)). The coefficient for the upper sponge layer was given as

$$\begin{array}{}\text{(B7)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{s},\mathrm{v}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{s},\mathrm{v}}^{-\mathrm{1}}}{\mathrm{2}}}\left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}+\tanh \left({\displaystyle \frac{z-(z_{\mathrm{T}}+z_{\mathrm{sp}})/\mathrm{2}}{{\mathit{\delta }}_{\mathrm{sp},\mathrm{v}}(z_{\mathrm{T}}-z_{\mathrm{sp}})}}\right)\right]\right\},\end{array}$$

whereas, for the lateral sponge layer,

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}{\mathit{\alpha }}_{\mathrm{s},\mathrm{h}}& ={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{s},\mathrm{h}}^{-\mathrm{1}}}{\mathrm{2}}}\left\{\left[\mathrm{1}-\tanh \left({\displaystyle \frac{\mathit{\lambda }-\mathit{\pi }/\mathrm{4}}{{\mathit{\delta }}_{\mathrm{sp},\mathrm{h}}\mathit{\pi }/\mathrm{2}}}\right)\right]\right.\\ & \left.+\left[\mathrm{1}+\tanh \left({\displaystyle \frac{\mathit{\lambda }-\mathrm{7}\mathit{\pi }/\mathrm{4}}{{\mathit{\delta }}_{\mathrm{sp},\mathrm{h}}\mathit{\pi }/\mathrm{2}}}\right)\right]\right\}\\ & \cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}+\tanh \left({\displaystyle \frac{\left|\mathit{\varphi }\right|-\mathit{\pi }/\mathrm{3}}{\mathrm{8}\mathit{\pi }/\mathrm{180}}}\right)\right],\end{array}\end{array}$$

where z_T is the height of model top and τ^s,v and τ^s,h are the decay timescales corresponding to the upper and lateral sponge layers. Note that the coefficient for the lateral sponge layer is multiplied by a tapering function in the latitudinal direction to avoid an infinite zonal scale near the poles, as performed in Eq. (48). In this study, z_sp=15 km, ${\mathit{\delta }}_{\mathrm{sp},\mathrm{v}}={\mathit{\delta }}_{\mathrm{sp},\mathrm{h}}=\mathrm{0.16}$, and ${\mathit{\tau }}_{\mathrm{s},\mathrm{v}}={\mathit{\tau }}_{\mathrm{s},\mathrm{h}}=\mathrm{100}$ s.

Figure B2 After 2 h in a two-dimensional mountain-wave test case, the spatial distribution of the (a) L₂ error and (b) local convergence rate for the vertical wind in the cases of the (a, b) numerically calculated Jacobian cofactor and (c, d) analytical Jacobian cofactor for p=3. In the L₂ error, we show the results obtained from the experiments with Δ_h=156 m and Δ_v=125 m. When evaluating the local convergence rate, we used the results obtained from two experiments: a coarse-resolution experiment (Δ_h=312 m, Δ_v=250 m) and the highest-resolution experiment for p=3 (Δ_h=39 m, Δ_v=31.25 m). In each figure, the white lines represent the vertical wind in the reference experiment.

B3 Investigation on the degradation of the optimal numerical convergence

In Fig. 6, the convergence rate obtained from the mountain test case was slightly smaller than that achieved for p+1-order accuracy. We consider the reasons behind this result to be as follows. First, to evaluate the differentials with the Jacobian cofactors ($\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{13}}$ and $\sqrt{G_{\mathrm{v}}}{G}_{\mathrm{v}}^{\mathrm{23}}$), we used same discretization operator, as described in Sect. 2.3. This strategy is beneficial to simply satisfy the geometric conservation law identity in the discretized equations. However, because the calculated geometric factors have the order p, it is possible to degrade the optimal convergence. Figure B1a and b show the spatial distribution of numerical errors for vertical wind and the local convergence rate, respectively, for p=3. The numerical error was large near the surface where the mountain exists. Furthermore, the relatively slow convergence rate appeared. The rate near the surface was between 2 and 3, while it approached the value of 4 at locations apart from the surface. Second, the modal filter can reduce the convergence rate during the long-term temporal integrations, even if we adopted a high-order modal filter with a relatively small decay coefficient.

To increase the certainty in our speculations, we conducted additional numerical experiments. To simplify the investigations and save computational resources, we treated the corresponding two-dimensional experimental setup. With respect to the Jacobian cofactors, we considered two cases: (i) the case where it is numerically given using the same discretization operator mentioned in Sect. 2.3 and (ii) the case when it is given by analytically evaluating the spatial derivatives at the node. In addition, to discuss the impact of modal filters on the convergence rate, we considered the case of no modal filter for p=3 because we found that the 2 h temporal integration without filters can be performed only for p=3. As performed with the global-model case, we conducted several numerical experiments, changing the spatial resolutions and polynomial orders. To evaluate the error norms, we used the results from the reference experiments with p=7, where horizontal and vertical grid spacing were Δ_h=78 m and Δ_v=62.5 m (z<15 km), respectively.

Figure B2a and b show the spatial distribution of numerical errors for vertical wind and the local convergence rate obtained from the two-dimensional experiments with p=3. As shown in the global experiment (see Fig. B1), the convergence rate near the mountain achieved only the third-order accuracy in the cases of a numerically calculated Jacobian cofactor. On the other hand, when the analytical Jacobian cofactor was used, the numerical errors near the mountain decreased and the convergence rate approached the fourth-order accuracy. Thus, the calculation strategy of a Jacobian cofactor is one of the reasons for sub-optimal convergence.

Figure B3 Dependence of the L₁, L₂, and L_inf errors on spatial resolution for the (a) density perturbation (ρ^′), (b) horizontal wind (u^ξ), (c) vertical wind (w), and (d) perturbation of potential-temperature-weighted density ((ρθ)^′) after 2 h in a mountain-wave test case with the two-dimensional experimental setup. Note that the cyan lines represent the results for the case of p=3 without the modal filter (MFoff).

Figure B3 shows that the dependence of the L₁, L₂, and L_inf errors on the spatial resolution. First, we focus on the results with p=3. When the metric cofactors were analytically evaluated and the modal filter was not used, the fourth-order accuracy was observed except for the density. In the case of a numerically calculated Jacobian cofactor, the convergence rates of the L₂ and L_inf errors were characterized by the sub-optimal order because the Jacobian cofactors have only pth-order accuracy, as mentioned above. Such behavior was observed for horizontal wind, vertical wind, and the perturbation of potential-temperature-weighted density based on the comparison between cases (i) and (ii). On the other hand, as indicated by the blue and cyan lines, the order reduction with the modal filters was obvious for the horizontal wind, while for other variables, there was little influence. This may be because the filters act not only on the perturbation part of horizontal wind but also on the mean flow part. For higher-order cases (p=7 and 11), the filters are unavoidable for ensuring numerical stability in a classical DGM. Then, the convergence rate can be limited by the modal filters, and the analytical Jacobian cofactor would have little impact. Even for case (ii), L₂ and L_inf errors in horizontal and vertical wind had a convergence rate slightly less than the optimal order. As for the density, note that the third-order accuracy was obtained for p=3 even when using the analytical Jacobian cofactor and removing the modal filter. It remains unclear why the density error does not decrease in the optimal order. We may need to pursue how to discretely deal with the hydrostatic balance (e.g., Li and Xing, 2018) and investigate the boundary errors with a not-normal flux condition near the surface.

Appendix C: Vertical grid stretching in the baroclinic-wave and Held–Suarez tests

In the baroclinic-wave and Held–Suarez tests, a function form for stretching the vertical-element size is similar to that in Ullrich and Jablonowski (2012b). We calculated the vertical coordinate ζ at the top element boundary of the k^′th element as

$$\begin{array}{}\text{(C1)}& {\mathit{\zeta }}_{k^{\prime }+\frac{\mathrm{1}}{\mathrm{2}}}=z_{\mathrm{T}}{\displaystyle \frac{\mathrm{1}}{\sqrt{b+\mathrm{1}}-\mathrm{1}}}\left[\sqrt{b{\left({\displaystyle \frac{k^{\prime }}{N_{\mathrm{e},\mathrm{v}}}}\right)}^{\mathrm{2}}+\mathrm{1}}-\mathrm{1}\right],\end{array}$$

where b is a positive parameter and N_e,z is the number of elements in the vertical direction. As b decreases, the vertical-element size near the surface reduces compared to the upper domain of the model. Here b=20 for p=3 and 7 and b=5 for p=11 in the baroclinic-wave test, while b=3 for p=3 and 7 and b=5 for p=11 in the Held–Suarez test.

Appendix D: The effect of not using shallow-atmosphere approximation in a global PBL turbulence experiment

In Sect. 3.6, we showed the results of the PBL turbulence experiment with shallow-atmosphere approximation. By applying this approximation, the obtained results become comparable with those reported in KT2021 and KT2023, which used the plane regional model. However, because the planetary radius is set to be 3.4 km, this approximation is not suitable for discussing physical aspects such as the impact of sphere geometry on the convective cells. This section shows the results when shallow-atmosphere approximation is not applied.

Figures D1 and D2 show the horizontal distributions of convective cells and vertical structure of the PBL when shallow-atmosphere approximation is not applied. In Fig. D2, all results with shallow-atmosphere approximation are represented by the gray shading. An open cell pattern with the characteristic horizontal scale of a few kilometers was observed irrespective of whether we apply shallow-atmosphere approximation. On the other hand, the winds became weaker and the PBL was shallower compared to that in shallow-atmosphere approximation. We consider such changes to be consistent with the effect of spherical geometry because horizontal area increases with altitude.

Figure D3 shows the energy spectra when shallow-atmosphere approximation is not applied. The major feature of energy spectra in the inertial sub-range remained mostly unchanged. For example, the wavelength range obeyed a power law of $-\mathrm{5}/\mathrm{3}$ and the relation of effective resolution with polynomial order.

Figure D1 Same as Fig. 13, but these results were obtained without shallow-atmosphere approximation. The horizontal distribution for vertical wind at z=500 m after t=4 h: (a) northern hemisphere (NH), (b) southern hemisphere (SH), and (c) the corresponding cross-sections along the equator.

Figure D2 Same as Fig. 14, but these figures show the vertical structure of the PBL temporally averaged during the last 30 min without shallow-atmosphere approximation: (a) potential temperature, (b) resolved eddy heat flux plus SGS heat flux, (c) variance of vertical wind, and (d) skewness of vertical wind. In these panels, the blue, red, and yellow lines represent the results for p=3, 4, and 7, respectively. We indicate the corresponding results with shallow-atmosphere approximation, shown in Fig. 14, by the gray shading.

Figure D3 Same as Fig. 15, but these figures show the results without shallow-atmosphere approximation: (a) density-weighted energy spectra E(n) of three-dimensional wind at a height of 500 m for p=3, 4, and 7. The dash-dotted gray line represents aE(n), where $a=\mathrm{8.0}\times {\mathrm{10}}^{\mathrm{1}}$. (b) Spectra normalized by the result of p=7. (c) Partial expanded view of energy spectra in the short-wavelength range.