An Adaptive Nonhydrostatic Atmospheric Dynamical Core Using a Multi-Moment Constrained Finite Volume Method

During the last few decades, increased attention has been paid toward improving numerical weather prediction (NWP) models for better resolving small-scale processes at high resolution. In terms of a dynamical core, much effort has been put forth to develop new numerical schemes and computational meshes with the intent of attaining higher accuracy and efficiency. Aside from the traditional methods, the high-order numerical schemes that are based on “local” or single-cell based spatial reconstruction are becoming more popular recently in developing the atmospheric dynamical cores because of their high-order accuracy, geometric flexibility, and their scalability on massive, parallel computer systems. The representative schemes are the spectral element method (Thomas and Loft, 2000; Iskandarani et al., 2002; Giraldo and Rosmond, 2004) and the discontinuous Galerkin method (Levy et al., 2007; Giraldo and Restelli, 2008; Nair et al., 2009; Blaise and St-Cyr, 2012) among others.

The multi-moment constrained finite-volume (MCV) method (Ii and Xiao, 2009) was recently proposed as a new framework to develop high-order schemes. Several point-wise values at equidistantly distributed solution points are defined as predicted variables (unknowns) within each computational cell. The formulations for the spatial discretization of unknowns are derived through the constraint conditions provided by introducing different types of moments (Yabe and Aoki, 1991; Yabe et al., 2001; Ii et al., 2005; Xiao et al., 2006; Ii and Xiao, 2007; Chen and Xiao, 2008), such as point values, volume-integrated averages, derivatives of different orders, and so on. Using the MCV method, 2D and 3D nonhydrostatic dynamical cores have been proposed in Li et al. (2013) and Qin et al. (2019).

Though computational power has made rapid gains recently, atmospheric modeling is still highly time-consuming and computationally expensive work. Considering the multi-scale phenomena at operating in the atmosphere, the development of adaptive numerical models, which are capable of running on the computational meshes with variable spatial resolution, is a very effective way to save on computational costs. Some existing methods can conduct multi-resolution simulations in atmospheric modeling, such as the nested grid (Pielke et al., 1992), the stretched grid (Yessad and Bénard, 1996), and the AMR grid (Berger and Oliger, 1984; Behrens, 1996; Stevens and Bretherton, 1996; Tomlin et al., 1997; Kessler, 1999). Using an AMR grid is very attractive among these techniques for its flexibility and efficiency in dynamic adjusting a computational mesh according to the time-evolution of predicted variables. AMR techniques have been applied in some atmospheric models to improve the computational efficiency, examples include the studies described in Skamarock et al. (1989), Skamarock and Klemp (1993), Blayo and Debreu (1999), Bacon et al. (2000), Hubbard and Nikiforakis (2003), Jablonowski et al. (2006), St-Cyr et al. (2008), and Chen et al. (2011). In this study, we will develop an adaptive MCV nonhydrostatic model with the applications of Berger and Oliger’s AMR algorithm based on the 2D model proposed in Li et al. (2013).

The remainder of this paper is organized as follows. A 2D nonhydrostatic dynamical core using a 3-point MCV scheme is briefly described in section 2. The setup and dynamic adjustment of an AMR mesh, along with numerical operations to exchange the solution information among different blocks under the framework of the MCV scheme, are described in section 3. Numerical results of the widely used benchmark tests are given in section 4, and a short summary is given in section 5.

A 2D compressible nonhydrostatic atmospheric dynamical core was proposed in Li et al. (2013) using a multi-moment constrained finite volume method. We implement this MCV dynamical core on the AMR grid in this study. The governing equations for the compressible nonhydrostatic dynamical core are written in conservative flux-form as

where the dependent variables are the density ($\rho $), the momentum vector $\left( {\rho {\boldsymbol{u}},\rho {\boldsymbol{w}}} \right)$, and the product of density and potential temperature $\rho \theta $; p is pressure and g is the acceleration of gravity. The effects of rotation are not considered in this implementation.

To deal with bottom topography, the height-based, terrain-following coordinates $\left( {x,\zeta } \right)$ proposed in Schär et al. (2002) are adopted. The Jacobians of the transformation are $\sqrt G = \dfrac{{\partial z}}{{\partial \zeta }}$ and ${G_{13}} = \dfrac{{\partial \zeta }}{{\partial x}}$. In computational space, the uniform grid $\left( {x,\zeta } \right)$ is related to the Cartesian coordinate as

where ${h_s}$ is the height of the bottom of the mountain, H is the model top, and the parameter S = 3000 m is adopted. $\widetilde w$ is the contravariant vertical velocity component in the transformed coordinate defined as

where u and w are velocity components for the Cartesian coordinate system $\left( {x,z} \right)$

Splitting of the hydrostatic reference state is usually applied in atmospheric modeling to improve accuracy. As a result, the thermodynamical variables are written as

A 3-point MCV scheme of third-order accuracy (Ii and Xiao, 2009) is adopted to solve the governing equations, Eqs. (1)–(4) in this study. Hereafter, we briefly describe the MCV formulations by considering the 1D hyperbolic conservation law having the form of

where ${\boldsymbol{q}}$ is the vector of dependent variables and ${\boldsymbol{f}}$ is the vector of flux function.

Within each cell, three local degrees of freedom (DOFs), in terms of point-wise values ${{\boldsymbol{q}}_{i,m}}$ (m = 1 to 3), are defined at two endpoints ${x_{i \pm \frac{1}{2}}}$ and the center ${x_i}$ as shown in Fig. 1a for the cell ${C_i}$.

Figure 1.  Distribution of DOFs for the (a) 1D case and (b) 2D Cartesian grid.

With known values of local DOFs, a quadratic Lagrangian polynomial ${{\boldsymbol{Q}}_i}\left( x \right)$ can be constructed as

where the Lagrangian basis function is

To build the MCV scheme, two kinds of moments including the PV (point value) moment ${\overline {{}^P{\boldsymbol{q}}} _{i - \frac{1}{2}}}$ at cell interface and the VIA (volume-integrated average) moment ${\overline {{}^V{\boldsymbol{q}}} _i}$ over each cell are defined to derive the spatial discretization formulations (Ii and Xiao, 2009). At the cell interface, the PV moment (same as the DOF defined at the cell’s endpoint) is updated using the differential-form governing equations through solving the derivative Riemann problem. The VIA moment is updated using flux-form finite volume formulations and assures the rigorous numerical conservation. The updating formulations for the DOF defined at the cell center are then derived through a constraint condition on the VIA moment. The details of the numerical procedure to derive the updating formulations for the 3-point MCV scheme can be referred to Ii and Xiao (2009). For the sake of brevity, we directly give the semi-discrete formulations for three different DOFs as

The numerical fluxes ${\widehat {\boldsymbol{f}}_{i \pm \frac{1}{2}}}$ at cell interfaces are C0-continuous (i.e., the zeroth derivative of numerical fluxes at cell interfaces are continuous) and can be directly computed from the values of known DOFs. The first-order derivatives of flux functions ${\widehat {\boldsymbol{f}}_{x,i \pm \frac{1}{2}}}$ are generally discontinuous and determined by solving derivative Riemann problems. In this study, a local Lax-Friedrich (LLF) solver (Shu, 1988) is adopted and the derivatives of flux functions at the cell interface $x = {x_{i + \frac{1}{2}}}$ are evaluated by

where a is the maximum absolute value of the eigenvalues of the Jacobian matrix ${{\partial {\boldsymbol{f}}}}/{{\partial {\boldsymbol{q}}}}$. The derivatives of dependent variables and flux functions are evaluated using the Lagrangian interpolation polynomials having the form of Eq. (10), and the subscripts L and R denote the cells ${C_i}$ and ${C_{i + 1}}$.

As shown in Ii and Xiao (2009), the 3-point MCV scheme shown above is of 3rd-order accuracy and preserves the rigorous numerical conservation in terms of the VIA moment given by

The MCV scheme can be easily extended to the multi-dimensional models by applying the above one-dimensional formulations in different directions one-by-one (Ii and Xiao, 2009). As shown in Fig. 1b, nine local DOFs in terms of point-wise values are defined within the cell ${C_{i,j}}$ to build the two-dimensional numerical models using the 3-point MCV scheme. Similarly, the numerical conservation is assured by the VIA moment ${\overline {{}^V{\boldsymbol{q}}} _{i,j}}$, which is determined by the two-dimensional Simpson’s rule as

where the indices i and j are omitted on the right-hand side of Eq. (17) for the sake of brevity.

For time-stepping, a third-order Runge-Kutta scheme is adopted to update the following semi-discrete ordinary differential equation (ODE) obtained from the above spatial discretization:

The solution is advanced from ${{\boldsymbol{q}}_{(n)}}$ at time ${t_{(n)}}$ to ${{\boldsymbol{q}}_{(n + 1)}}$ at time ${t_{(n + 1)}}$ by

with $\Delta t = {t_{(n + 1)}} - {t_{(n)}}$ as the time step. The time integration strategy for the AMR grids will be demonstrated in later contents.

The details of the numerical formulations for a two-dimensional MCV nonhydrostatic dynamical core are described in Li et al. (2013).

The AMR algorithm proposed by Berger and Oliger (1984) uses a set of blocks to cover the whole computational domain. On each block, the structured computational mesh is constructed to implement the spatial discretization. These blocks are grouped into different levels according to their spatial resolutions and the blocks with finer resolution always overlay the blocks belonging to the next coarser levels. The governing equations can be solved on each block independently if enough ghost cells, which are either copied from the adjacent blocks of the same level or evaluated by coarse-fine interpolations from blocks of the next coarser level, are supplemented. Additionally, the AMR algorithm also implements the dynamic adjustment of blocks according to the time evolution of the predicted fields. Since the high-order discretization is implemented using the single-cell based stencils in MCV models, the coarse-fine interpolations are easily accomplished and extra numerical errors due to the non-uniform resolution can be more effectively suppressed in comparison with existing variable-resolution models, as shown in our previous study (Chen et al., 2011). The operations to generate and adjust the AMR grid along with some additional numerical procedures to implement the MCV model on an AMR grid are briefly introduced in this section. Without losing generality, here, we consider the uniform Cartesian grid $\left( {x,y} \right)$ on each block, i.e., the grid spacing is ${h_k} = \Delta {x_k} = \Delta {y_k}$ on level k.

In the beginning, a base block (level 0) is first built, which covers the whole computational domain and has the coarsest resolution. Then, finer levels are generated one by one following a similar procedure. Given the kth-level blocks, the next finer level is generated as follows.

All cells belonging to level k are checked by some refinement criterion to find those cells where the physical fields should be solved using finer resolution to improve the accuracy. Berger and Oliver (1984) introduced a criterion based on the error estimation using Richardson extrapolation. Though it is a general methodology for different types of applications, a more efficient criterion is to check the known flow structures of some representative physical variable in atmospheric modeling. In this study, we flag the cell to be refined if the variation of some indicative variable across the cell is larger than a prescribed threshold.

The simplest way is to use just one block, which is the smallest rectangle (in 2D computational space) is to cover all flagged cells. However, it is inefficient since, in general, the domain consisting of the flagged cells is highly irregular, and using a single structured block will result in many un-flagged cells being involved at the finer level $k{\text{ + 1}}{\text{.}}$ To improve the efficiency, we divide one big block into two small ones and it is expected that the area of each small one then can be reduced to cover the flagged cells. The division of an existing block will be conducted continuously until the ratio between the area covered by un-flagged cells and the total area on this block is smaller than a prescribed value. In this study, we adopted a division method based on pattern recognition (Berger and Rigoutsos, 1991) instead of a simple bisection to speed up this operation.

The above procedure is repeated until the finest level has been constructed. During the simulations, the computational meshes are completely or partially adjusted every several time steps to track the evolution of predicted fields. This adjustment starts from the fine level and most of the above operations can be applied with some modifications. When the grid adjustment is finished, values of predicted variables in the newly generated fine cells are interpolated from the coarser level. Coarse-fine interpolation operations based on a single-cell stencil are described later in this section.

In this study, a time marching step is chosen to keep the same CFL (Courant-Friedrichs-Lewy) number on different levels. In two-dimensional case, the CFL number on level k is $\frac{{u\Delta {t_k}}}{{\Delta {x_k}}} + \frac{{v\Delta {t_k}}}{{\Delta {y_k}}}$. Considering that the refinement ratio of grid resolution is an integer, $\lambda $, between levels $k$and $k + 1$, the corresponding grid spacings and time steps on these two adjacent levels are ${h_k} = \lambda {h_{k + 1}}$ and $\Delta {t_k} = \lambda \Delta {t_{k + 1}}$. With variable time steps on different levels, the governing equations are advanced for $\lambda $ steps on level $k + 1$ during one step on level $k$. In Fig. 2, the time marching procedure from ${t_n}$ to ${t_{n + 1}} = {t_n} + \Delta t$ is illustrated for a 3-level AMR grid. Each level is advanced through Eq. (19) with its own time integration interval. For the sake of simplicity, the refinement ratio is chosen as $\lambda = 2$ hereafter in all illustrative examples. A special recursive procedure is adopted for time marching on an AMR grid, i.e., any coarse level cannot be advanced to the next time step until all finer levels have reached the current time instant. Additionally, the coarse-fine interpolations are conducted both in time and space to evaluate the ghost cells. In this study, the temporal interpolation is accomplished by a linear reconstruction of the predicted fields at two adjacent time steps.

Figure 2.  Diagram depicting the subcycling of AMR levels in time.

To guarantee the numerical conservation, flux correction should be conducted along the interfaces between the coarse and fine levels (Berger, 1989, 1998). As shown in Fig. 3, $F_{i + 1,j,k}^{}$ represents the time-averaged numerical flux during one time step on level $ k $, and is used to update $ \overline {{}^Vq} _{i,j,k}^{} $, the VIA moment of the coarse cell $C_{i,j,k}^{}$. During the same time interval, VIA moments of the adjacent fine cells $ C_{2i + 1,2j,k + 1}^{} $ and $C_{2i + 1,2j - 1,k + 1}^{}$ are advanced by time-averaged numerical fluxes $f_{2i + 1,2j,k + 1}^{}$ and $f_{2i + 1,2j - 1,k + 1}^{}$. Since the relation $F_{i + 1,j,k}^{} = {(f_{2i + 1,2j - 1,k + 1}^{} + f_{2i + 1,2j,k + 1})}/{2}$ is not assured by the raw scheme, the AMR model is not conservative without a flux correction operation.

Figure 3.  Flux correction along the interface between the coarse and fine levels.

In the proposed model, the VIA moment of a coarse cell is corrected to assure the numerical conservation by:

where $A_{i,j,k}^{}$ is the area of the coarse cell, $C_{i,j,k}^{}$.

The solution transfer between coarse and fine levels plays an important role to accomplish the numerical modeling on an AMR grid. In the proposed model, numerical solutions on level $k$ are only exchanged with the two adjacent levels, $k - 1$ and $k + 1$.

Since the numerical solutions with finer resolution are of higher accuracy, they are used to replace the solutions in the overlap regions on the next coarser level through a fine-to-coarse solution transfer. As shown in Fig. 4, at the solution points defined along the edges of the coarse cell (denoted by solid circles), the values of predicted variables are directly copied from the corresponding solution points (denoted by hollow circles) on the finer level. At the center of the coarse cell, the point-wise values (denoted by solid triangles) are evaluated to preserve the numerical conservation, i.e., it is calculated through Eq. (17) with the VIA obtained through the relation

Figure 4.  Fine-to-coarse solution transfer in the overlap regions.

The coarse-to-fine solution transfer is adopted to supplement the ghost cells and to evaluate the solutions within the newly generated fine cells.

In Fig. 5, the coarse-to-fine solution transfer to supplement the ghost cells is illustrated. A coarse cell $C_{i,j,k}^{}$ is adjacent to the fine cells $C_{2i + 1,2j,k + 1}^{}$ and $C_{2i + 1,2j - 1,k + 1}^{}$ on level $k + 1$. Fine cells $C_{2i,2j,k + 1}^{}$ and $C_{2i,2j - 1,k + 1}^{}$ are ghost cells, where the solutions should be interpolated from the coarse cell $C_{i,j,k}^{}$ to provide the boundary conditions. Ghost points denoted by solid circles coincide with solution points within $C_{i,j,k}^{}$, and the values of predicted variables can be directly copied from the known solutions. Ghost points, denoted by solid squares, are not solution points of level $k$, and values of physical fields have to be evaluated through single-cell ($C_{i,j,k}^{}$) based Lagrangian interpolation polynomials, which are 2D quadratic polynomials. For some variable, $Q$, these take the form of

Figure 5.  Coarse-fine interpolations based on multi-moments for evaluating the ghost cells.

where the coefficients are decided by nine constraint conditions provided by known DOFs.

The coarse-fine interpolations to evaluate predicted variables in newly generated fine cells are more complicated, as the rigorous numerical conservation should be preserved. As shown in Fig. 6, the coarse cell $C_{i,j,k}^{}$ is divided into four fine cells $C_{2i - 1,2j - 1,k + 1}^{}$, $C_{2i - 1,2j,k + 1}^{}$, $C_{2i,2j - 1,k + 1}^{}$, and $C_{2i,2j,k + 1}^{}$. The predicted variables should be evaluated at all solution points on level $k + 1$ (denoted by solid shapes). At solution points denoted by solid circles and squares, the same algorithm is adopted either by copy or interpolation, as mentioned above, to set up the fine ghost cells. However, at internal points denoted by solid triangles, the point-wise values should be determined to preserve the numerical conservation. First, the VIAs of some variable $q$ within four fine cells are evaluated by

Figure 6.  Coarse-fine interpolations based on multi-moments for evaluating the flow variables in newly generated fine cells.

Then the point-wise value at the cell center denoted by the solid triangle is obtained through the two-dimensional Simpson’s rule [see Eq. (17)] in each newly generated fine cell.

To check the performance of the adaptive model in improving the computational efficiency, the widely used benchmark test cases are carried out in this section. The performance of the MCV dynamical core on a fixed grid has been validated in Li et al. (2013) and the numerical results are comparable to those of existing advanced models. In this study, we run the MCV dynamical core on adaptive grids with different configurations. The adaptive grid is denoted by $N \times M \times l \times \lambda $. N and M represent the number of cells along the horizontal and vertical directions respectively, which means the base block (the coarsest level) consists of $N \times M$ computational cells, the finest level is level $l$ and the refinement ratio between two adjacent levels is $\lambda $. It is noted that the minimum grid spacing is allowed by the finest level l. The normalized ${l_1}$, ${l_2}$ and ${l_\infty }$ errors (Williamson et al., 1992) and computational costs are examined. Definitions of the normalized errors are given as follows:

where $\varOmega$ is the computational domain, $\overline q $ and $ {\overline q _{\text{e}}} $ are the numerical and reference solutions in terms of cell-integrated averages. All numerical tests in this study are carried out on the AMD EYPC 7702 CPU (single processor).

The initial hydrostatic state is specified as

where the Exner pressure is given by $\Pi = {\left( {\frac{p}{{{p_0}}}} \right)^{\frac{{{R_d}}}{{{c_p}}}}}$. The constants are $ {c_p} = 1004.5{\text{ J k}}{{\text{g}}^{ - 1}}{\text{ }}{{\text{K}}^{ - 1}} $ and ${R_d} = 287{\text{ J k}}{{\text{g}}^{ - 1}}{\text{ }}{{\text{K}}^{ - 1}}$.

The basic state of potential temperature, $\theta $, is chosen to be a constant

or specified with a constant Brunt-Väisälä frequency as

where $N_0^2 = g {{{\rm{d}}{\text{ln}}\overline \theta }}/{{{\rm{d}}z}}$. The hydrostatic density $\overline \rho $ is then computed by:

where $\overline T = \frac{{\overline \theta }}{{\overline \Pi }}$.

The refinement criterion is chosen as the variation of the potential temperature perturbation for the rising thermal bubble, density current flow, and internal gravity waves, as well the variation of horizontal velocity for the Schär mountain case. Considering a two-dimensional cell $C_{i,j,k}^{}$ on level $k$, it is flagged to be refined if

where $\delta $ is a prescribed threshold, and $q$ can be the potential temperature perturbation, $\theta '$, or the horizontal velocity, $u$.

Atmospheric motions caused by thermodynamic effects (Carpenter et al., 1990; Wicker and Skamarock, 1998; Ahmad and Lindeman, 2007; Norman et al., 2011) are common phenomena in nature. They are extensively used to verify the performance of dynamical cores. In this test, a thermal bubble that is warmer than the ambient air is initialized by specifying a positive potential temperature perturbation as

where $r = \sqrt {{{\left( {x - {x_0}} \right)}^2} + {{\left( {z - {z_0}} \right)}^2}} $ and $ R = 2000 $ m (the radius of the bubble). The initial center of the thermal bubble is located at $ \left( {{x_0},{\text{ }}{z_0}} \right) = \left( {10000{\text{ m}},{\text{ }}2000{\text{ m}}} \right) $.

In this case, the uniform background potential temperature is specified as $\overline \theta = 300{\text{ }}K$, the maximal perturbation is $\Delta \theta = 2{\text{ }}K$ and the computational domain is [0, 20] km $ \times $ [0, 10] km. The simulation runs for 1000 s and all boundaries are slippery walls. The explicit dissipation filter (Li et al., 2013) with a viscosity coefficient, $\mu = 10{\text{ }}{{\text{m}}^2}{\text{ }}{{\text{s}}^{ - 1}}$, is used to eliminate numerical noise. The threshold for refinement is set as $\delta = 0.04$. During the simulation, the thermal bubble rises and finally attains a mushroom-like shape. The numerical results on a uniform grid, with a resolution of 25 m $ \times $ 25 m, are used as a reference solution to calculate the normalized errors.

Normalized errors and elapsed CPU time for different grids are given in Table 1, where the results are grouped according to the finest resolution on an adaptive grid. The normalized CPU time is also computed by dividing the CPU time on the coarsest uniform grid. From Table 1, it can be found that the normalized errors of the numerical results on the uniform and AMR grids within each group are quite similar, noting that much less CPU time is consumed by the AMR model. Contour plots of the potential temperature perturbation $\left( {\theta '} \right)$ and the absolute errors on a $100 \times 50 \times 3 \times 2$ grid at different times are shown in Fig. 7. The solid thick lines represent the interfaces between different levels. The symmetric distribution of $\theta '$ is perfectly reproduced on the adaptive grid as in our previous studies (Li et al., 2013). It is observed that the fine levels are dynamically adjusted to fit the change of flow field. AMR grids properly locate the disturbed region and put the fine blocks there to assure numerical accuracy. In undisturbed regions, the coarse blocks are applied to save computational costs. The relative total mass error along the simulation is also recorded, which supports that the numerical conservation property is well achieved through the flux correction procedure.

Figure 7.  Contour plots of a potential temperature perturbation ($\theta '$) and its errors for a rising thermal bubble on $100 \times 50 \times 3 \times 2$ grid at (a) t = 0 s, (b) t = 250 s, (c) t = 500 s, (d) t = 750 s, (e) t = 1000 s, and (f) absolute errors.

In the density current test case (Straka et al., 1993; Giraldo and Restelli, 2008), a cold bubble is put in a neutrally stratified atmosphere by specifying a negative potential temperature perturbation as

where

$ \overline \theta = 300{\text{ K}} $, $\Delta \theta = - 15{\text{ K}}$, $\left( {{x_0},{\text{ }}{z_0}} \right) = \left( {0 , 3000 } \right)\;{\rm{m}}$ and $\left( {{x_r},{\text{ }}{z_r}} \right) = \left( {4000 , 2000 } \right)\;{\rm{m}}$.

During the simulation, the cold bubble drops to the ground and spreads out in the horizontal direction to form three Kelvin-Helmholtz shear instability rotors along the cold frontal surface. The computational domain of this case is [–26.5, 26.5] km $ \times $ [0, 6.4] km and the simulation time is 900 s. All of the boundaries are slippery walls. A dissipation filter with a coefficient of $\mu = 10{\text{ }}{{\text{m}}^2}{\text{ }}{{\text{s}}^{ - 1}}$ is used to meet the requirement of a physical process. The threshold of the refinement criterion of this case is $\delta = 0.2$. Again, the numerical results on a uniform grid with a resolution of 25 m × 25 m are adopted as reference solutions.

Normalized errors of $\theta '$ and elapsed CPU time on different grids are given in Table 2. It can be observed that the AMR model consumes much less computational costs at the price of slightly larger normalized errors. Contour plots of $\theta '$ and the absolute errors on a $132 \times 32 \times 2 \times 4$ grid are shown in Fig. 8, which have the finest resolution of 50 m $ \times $ 50 m. Only half of the computational domain is depicted due to the symmetry of the flow field. The numerical result agrees well with those on the uniform grid with the same resolution given in Li et al., (2013).

Figure 8.  Contour plots of the potential temperature perturbation ($\theta '$) and its errors for the density current on $132 \times 32 \times 2 \times 4$ grid at (a) t = 0 s, (b) t = 225 s, (c) t = 450 s, (d) t = 675 s, (e) t = 900 s, and (f) absolute errors.

The internal gravity waves test involves the evolution of a potential temperature perturbation in a channel with periodic boundary conditions in the horizontal direction. Initial conditions used in Skamarock and Klemp (1994) are adopted. The potential temperature field is initialized as

where $ {\theta _0} = 300{\text{ K}} $, $H = 10{\text{ km}}$, $\Delta \theta = 0.01{\text{ K}}$, ${x_0} = 100{\text{ km}}$, and $a = 5{\text{ km}}$.

The background state of potential temperature $ \overline \theta \left( z \right) $ is obtained by using a constant Brunt-Väisälä frequency. A constant mean flow of $\overline u = 20{\text{ m }}{{\text{s}}^{ - 1}}$ for the advection of the internal gravity waves is adopted. The bottom and top boundaries are slippery walls. The computational domain of this case is [0, 300] km $ \times $ [0, 10] km, and the simulation time is 3000 s. $\delta = 1.8 \times {10^{ - 4}}$ is chosen for grid refinement. Numerical results on a 250 m $ \times $ 25 m grid are taken as the reference solutions. It is noted that the computational cells are no longer square, an aspect ratio of grid spacing $\Delta x = 10\Delta z$ is adopted in this case.

Normalized errors of $\theta '$ and elapsed CPU times on various grids are given in Table 3. Comparing to the first two test cases, the effect of AMR grids in saving computational costs is less obvious. The reason is that the internal gravity waves are spreading in the horizontal direction during the simulation and the disturbed regions are continuously growing as shown in Fig. 9. Maximum and minimum values of vertical velocity and potential temperature perturbation on the AMR grids with the finest resolution of $\Delta z = 100{\text{ m}}$after 3000 s are given in Table 4, which is the same as that obtained by Li et al (2013). The distribution of the absolute errors of $\theta '$ on a uniform grid and a 75 × 25 × 3 × 2 grid at $t = 3000{\text{ s}}$ is depicted in Figs. 10a, b. No obvious increase of errors is found around the boundaries between different levels, which reveals that the computational mode due to the change of grid resolution is effectively suppressed in the proposed AMR model. The general errors in the AMR model are also affected by the prescribed refinement threshold. More accurate solutions are obtained when a larger area of the computational domain is refined with a smaller refinement threshold.

Figure 9.  Contour plots of potential temperature perturbation ($\theta '$) for the internal gravity waves on $75 \times 25 \times 3 \times 2$ grid at (a) t = 0 s, (b) t = 750 s, (c) t =1500 s, (d) t = 2250 s, and (e) t = 3000 s.

Figure 10.  Error distribution of the potential temperature perturbation ($\theta '$) for the internal gravity waves on the (a) uniform grid and (b) $75 \times 25 \times 3 \times 2$ grid at t = 3000 s.

The Schär mountain case (Schär et al., 2002) is used to test the ability of a model to deal with the effects of complex terrain. A mountain with five peaks is used as topography which is specifically defined as

where ${h_0} = 250{\text{ m}}$, ${a_0} = 5000{\text{ m}}$, and ${\lambda _0} = 4000{\text{ m}}$. An initial background state of potential temperature is obtained by using a constant Brunt-Väisälä frequency of ${N_0} = {10^{ - 2}}{\text{ }}{{\text{s}}^{ - 1}}$ and ${\theta _0} = 280{\text{ K}}$. A constant mean flow of $\overline u = 10{\text{ m }}{{\text{s}}^{ - 1}}$along the horizontal is also initialized. This case is running on a domain of [–25, 25] km $ \times $ [0, 21] km for 10 hours. Non-reflecting boundary conditions are implemented by putting sponge layers in the area of $z \geqslant 9$ km for the top boundary and $ \left| x \right| \geqslant 15 $ km for the lateral inflow and outflow boundaries. The bottom boundary is a slippery wall. A refinement threshold, $\delta = 0.3$, is applied to this case.

The profile of the Schär mountain is depicted in Fig. 11e, noting that only part of the computational domain is displayed. Contour plots of the horizontal velocity, $u$, at various stages in the simulations are displayed in Figs. 11a–d. Since the disturbed regions quickly extended to cover the entire computational domain, it is not expected that the computational costs can be significantly saved by AMR models. Numerical results on a uniform fine grid with a resolution of 125 m × 105 m are used as reference solutions and the corresponding normalized errors are given in Table 5. The elapsed CPU time of the AMR model is about 44% of that on the uniform grid. The AMR grid does not affect accuracy since the numerical errors after 10 hours are almost the same as those on the uniform grid with the same resolution.

Figure 11.  Contour plots of horizontal velocity ( u ) of the Schär mountain case at (a) t = 7200 s, (b) t = 16800 s, (c) t = 26400 s, (d) t = 36000 s, and (e) the shape and size of the Schär mountain.

In this paper, an adaptive nonhydrostatic atmospheric dynamical core is proposed by using the multi-moment constrained, finite-volume method and Berger-Oliger’s adaptive mesh refinement algorithm. The high-order scheme is constructed based on two kinds of local degrees of freedom. As a result, the spatial reconstruction for a 3rd-order scheme is accomplished based on a single-cell-based stencil. The compact spatial stencil is beneficial for the efficient implementation of coarse-fine interpolation and reduction of the extra numerical errors due to the non-uniform resolution of an AMR grid. The proposed model runs on the adaptive grid consisting of blocks with different resolutions. The blocks with fine resolutions are placed in those regions that contain large variations of physical fields in order to improve the numerical accuracy, while the blocks with coarse resolutions are placed in the remainder of the area to save computational costs. With the evolution of flow fields, the distributions of blocks are dynamically adjusted to guarantee that the regions with complex flow structures are always resolved with fine resolutions. Numerical results of widely used benchmark tests reveal the effectiveness of the AMR technique in saving computational costs. It is observed that the use of an adaptive grid has only a very slight negative impact on computational accuracy in comparison with results on fixed grids with the same resolution.

Today, finer predictions of small-scale weather and climate phenomena and their interactions with large-scale ones gain more and more attention. Though the supercomputing capacity has been greatly increased recently, the use of an AMR grid would be beneficial in reducing unnecessary computational cells and make it more affordable to accurately track spatially complex and time-dependent small-scale phenomena. However, using an AMR grid has the undesirable side-effect of making it more difficult to achieve high scalability, since it is not easy to keep satisfying load balance among different CPUs on time-varying computational meshes. It is a major challenge to construct a practical atmospheric model using adaptive mesh.

Promising results by the proposed adaptive model have shown great potential in applying the proposed numerical framework to develop a more efficient atmospheric dynamical core. Future research will mainly focus on the extension of an adaptive model to the three-dimensional regional and global dynamical cores.

Acknowledgements. This work is supported by The National Key Research and Development Program of China (Grants Nos. 2017YFA0603901 and 2017YFC1501901) and The National Natural Science Foundation of China (Grant No. 41522504).