HTML
--> --> -->Recently, another type of high-order scheme with local reconstructions, the so-called multi-moment constrained finite-volume method (MCV), has been developed (Ii and Xiao, 2009). The basic idea of this method is the multi-moment concept (Yabe and Aoki, 1991; Yabe et al., 2001; Xiao, 2004; Ii et al., 2005; Xiao et al., 2005, 2006; Xiao and Ii, 2007), where we make use of different kinds of discrete quantities, such as the point-wise value (PV), the volume-integrated average (VIA), and the derivatives of different orders, which are called “moments” in our context, to describe physical fields. With these moments locally defined, high-order reconstructions can be built in compact stencils. Instead of updating the moments directly (Chen and Xiao, 2008; Li et al., 2008; Akoh et al., 2010), the MCV formulation updates values at the points as the basic model variables located in each mesh element (Li et al., 2013b). The values at solution points are predicted by time evolution equations obtained from imposing a set of constrained conditions on different types of moments. In the MCV method, some constraint conditions are applied on the boundaries of the computational element, where the derivatives of the numerical flux are solved as derivative Riemann problems.
In comparison with other existing high-order schemes, the MCV method allows a less restrictive CFL condition for computational stability (Ii and Xiao, 2009). Besides, an MCV model has exact numerical conservation through a constraint on the VIA moment and shows superiority in computational efficiency and algorithmic simplicity. So far, this method has been successfully implemented to develop 2D global advection and shallow models on different grids (Ii and Xiao, 2010; Chen et al., 2012; Li et al., 2012; Li et al., 2013a; Chen et al., 2014a, b; Xie and Xiao, 2016), as well as a nonhydrostatic dynamic core on a Cartesian grid (Li et al., 2013b, 2016), which has been demonstrated as a robust high-order method for a new type of numerical framework. Zhang et al. (2017) also showed that the MCV method has good scalability on high-performance computers, thus allowing massively parallel computing with thousands of process cores.
In this paper, given the attractive results of the previous 2D three-point MCV nonhydrostatic atmospheric dynamics (Li et al., 2013b), we further extend it to a 3D nonhydrostatic dynamic core and verify its performance, which will set the stage for the development of practical global atmospheric modeling. The remainder of the paper is organized as follows: In section 2, the 3D governing equations with the effect of topography are described in detail. The basic numerical formulation of the three-point MCV scheme follows in section 3. Numerical results from benchmark tests are presented in section 4. And finally, a short conclusion is given in section 5.
where ρ is the air density; V = (u, v, w)T is the 3D velocity field in the Cartesian coordinate; f is the Coriolis parameter; and θ is the potential temperature, which is related to the air temperature T and pressure p by



2
2.1. Splitting of reference profile
As generally used in atmospheric modeling, thermodynamic variables are written as the sum of the mean state and perturbation (Skamarock and Klemp, 2008):where the mean state (denoted by the overbar) satisfies local hydrostatic balance:
Under this splitting, Eqs. (1)?(5) take the form
where

2
2.2. Incorporation of topography
Owing to the advantages of accuracy and flexibility, a terrain-following height-based coordinate is preferred in many nonhydrostatic models (Prusa et al., 2008; Skamarock et al., 2012), as originally proposed by Gal-Chen and Somerville (1975). Here, we use a more sophisticated terrain-following coordinate system (Sch?r et al., 2002) to map the Cartesian space (x, y, z) to computational space (x, y, ζ)where h is the topographic height, H is the altitude of the model top, and s is the scale height parameter.
The coordinate transformation enters the governing equations through the metric coefficients G1,3 and G2,3, and the Jacobian of the transformation

For an arbitrary field variable ?, the following chain rules (Clark, 1977) connecting the (x, y, z) with the (x, y, ζ) coordinate system are used:
Thus, the divergence operator of a flux F = (Fx, Fy, Fz) takes the following form:
The vertical velocity in the transformed coordinate is
2
2.3. Governing equations
Substituting the relations of Eqs. (16)?(21) into Eqs. (10)?(14), the governing equations in the transformed coordinate are written as follows:We rewrite Eqs. (22)?(26) into the following vector form:
where


and
are the flux vectors in the x, y and ζ directions, respectively.
2
2.4. Eigenvalues and eigenvectors of the governing equations
As in section 3, the eigenvalues and eigenvectors of the hyperbolic parts of the governing equations are needed when we compute the derivative Riemann problems. For convenience, we briefly discuss how to calculate them here. Considering the homogeneous part of the governing equations, Eqs. (22)?(26), the set of hyperbolic equations in 3D is given as follows:The linearized form of Eq. (28) is written as
where A, B and C are the Jacobian matrices corresponding to flux function f, g and h, which are calculated, respectively, by
where

For hyperbolicity, the eigenvectors and eigenvalues can be obtained from the Jacobian matrix. The eigenvalues of Jacobian matrix A, B and C are computed by
where λx, λy and λζ are the eigenvalues and I is the identity matrix.
The solution of the above determinants gives the respective eigenvalues:
The corresponding right and left eigenvectors (the inverse matrix of right eigenvectors) are as follows:
where

where q is a dependent state variable and f(q) is the corresponding flux function.
The 1D computational domain,




where

To build the MCV scheme, two kinds of moments—namely, the VIA moment over the segment and the PV moment at the two ends—are defined as
where

where




Given qi,s (s = 1, 2, 3), a cell-wise Lagrangian interpolation function Qi(x) can be built,
where
is the Lagrangian basis function.
Given interpolation function (46), Eqs. (40)?(42) can be explicitly expressed in the unknowns,
Using the relations of Eqs. (48)?(50), the evolution equations of different moments (43)?(45) are obtained as
The constraint conditions (51)?(53) finally lead to the equations to update the unknowns as follows:
At boundary points



where








The same procedure also applies for cell boundary

We write Eqs. (54)?(56) into the following vector-matrix form:
Denoting the entries of the matrix
by Ml,α (l =1, 2, 3 and α =1,…, 4), and the components of
by Fα, Eq. (60) can be written into a component form,
The 3D computational domain D is decomposed into non-overlapping mesh elements Ci,j,k indexed by i, j and k in the x, y and ζ directions respectively. Thus, we have
where the cell element spans over
with
and
and I, J and K represent the total numbers of cell elements in the x, y and ζ directions. The solution points (shown in Fig. 1b) for cell Ci,j,k are denoted by qi,j,k,l,m,n (l, m, n = 1, 2, 3) and defined at points (xi,j,k,l, yi,j,k,m, ζi,j,k,n). For convenience, we drop the cell indices i, j and k from now on and focus on the local control volume. As addressed in Ii and Xiao (2009), multi-dimensional MCV formulation in regular Cartesian grids can be obtained by directly applying 1D formulation along each direction. The semi-discretized formulation for updating the unknowns within cell element Ci,j,k in the 3D MCV framework in Eq. (27) is then obtained as
The terms


A similar expression is obtained in the ζ direction:
and S(ql,m,n) stands for the source terms.
Analogously,




So far, we have described the spatial discretization procedure of the three-point MCV scheme, and the resulting semi-discrete system of Eq. (65) is then solved by the third-order Total Variation Diminishing Runge?Kutta method (Shu, 1988) for time integration.
We denote the right-hand side of Eq. (65) by R(ql,m,n):
The values


where

2
4.1. Boundary conditions
In the following test cases, two types of boundary conditions are used, including a no-flux boundary and nonreflecting boundary.The no-flux boundary condition aims to omit the normal component of the velocity and only keep the tangential component. For velocity vector u, we enforce
where n is the outward normal vector from the boundary. Specifically, the velocity vector on one side of the boundary is equal to its negative value on the other side, while the scalar variables are copied directly.
The commonly used method of imposing nonreflecting boundary conditions involves adding absorbing layers (damping layers) along the lateral and top boundaries, which guarantee the waves exiting the domain of interest without reflection. This damping is added to the momentum and potential temperature evolution equations and always takes the following form:
where τ is the coefficient of the absorbing layer and qb is the prescribed reference state. The strength of the damping is determined by parameter τ, which is defined differently in different works. Here, following Giraldo and Restelli (2008), it is denoted as
where τ0 = 2.0 × 10?2 s?1,

2
4.2. Test cases
34.2.1. Rising thermal bubble
This test is taken from Marras et al. (2013). The domain is bounded within [0, 3200] m × [0, 3200] m × [0, 4000] m. The mesh resolution is set to 80 m. This test case uses a hydrostatically balanced reference state based on a uniform potential temperature θ0 = 300 K. No-flux boundary conditions are applied along all six boundaries. The initial u velocity, v velocity and w velocity are set to zero. The potential temperature perturbation has the following form:where

The initial potential temperature perturbation field and the simulated results after 480 s are shown in Fig. 2. The thermal bubble (Fig. 2a) is located in the center of the domain in the beginning. The warm perturbation generates acceleration in the inner region of the bubble, accompanied by downdrafts on either side of the bubble. Since the highest temperature inside the bubble is located in the center, the center of the bubble rises faster. This creates sharp temperature gradients in the upper part of the bubble. After 480 s, it evolves into a classical mushroom shape. Similar to previous studies such as Chen et al. (2013) and Marras et al. (2013), the present model captures the evolution of the bubble in time, and the flow structures are reproduced with high solution quality. Figures 2b?d display the xz cross section of potential temperature perturbation, horizontal wind and vertical wind at 480 s. As expected, the evolution of potential temperature perturbations remains axisymmetric at all times and is very competitive when compared to that of Marras et al. (2013) (see their Fig. 19), which applies the variational multiscale stabilization (VMS) method to suppress the instabilities of dominated convection in the centered-nature finite element framework. Different from the VMS method, the MCV scheme combined with the approximate Roe solver presents small oscillations of potential temperature perturbations in the larger gradient of the mushroom edge. The horizontal wind and vertical wind at 480 s look closely like those of Marras et al. (2013) as well. In a quantitative manner, the front of the bubble finally reaches the height of approximately 2400 m, which agrees well with Marras et al. (2013). By comparing the thermal bubble test results with previous studies, the results of the 3D MCV nonhydrostatic model look very competitive.

3
4.2.2. 3D linear hydrostatic mountain
In this test, the initial state consists of an isothermal atmosphere with

where hc = 1 m is the mountain height, ac = 10 km is the half-width of the mountain, and (xc, yc) = (120, 120) km is the center position of the mountain. The domain is [0, 240] km × [0, 240] km × [0, 24] km, with grid spacing of ?x = ?y =1500 m and ?ζ = 300 m. While no-flux boundary conditions are imposed along the bottom boundary, nonreflecting boundary conditions are used near the upper and lateral boundaries. The sponge layers are placed in the region of ζ ≥18 km for the upper boundary and the thickness of 60 km for the lateral boundaries. As mentioned in section 2.2, the hybrid vertical coordinate (Sch?r et al., 2002) is employed. It can be verified that the Froude number U/(Nh0) > 1, so the flow is in the hydrostatic range. The model runs for five hours.
Figures 3a and b show the numerical results of the u velocity perturbation and vertical velocity after one hour, respectively. The results present the features of vertically propagating, with a small tilt component, which are in good agreement with the results reported in Giraldo (2011). The density perturbation field after five hours is plotted in Fig. 4a. Thanks to the theory of Smith (1988), the analytic solution of this problem has been made available. Our results are well matched with the analytic solution (Smith, 1988) near the mountain, but deviate away from the mountain due to the influence of the Rayleigh damping layer. We should not expect to agree exactly with the analytic solution, for it is solved under the linear Boussinesq approximation. Therefore, the results of previously published work are also used for comparison. We observe agreement with the results of Blaise et al. (2016) (see their Fig. 16a), suggesting our model is correctly capturing the mountain waves triggered by the gentle slope. It is worth mentioning that, compared with the results of Blaise et al. (2016) obtained by the discontinuous Galerkin method, our simulated density perturbation field shows small oscillations, especially at the lower levels. Additional verification is performed by comparing the density perturbations along the line x = y = 120 km. It is noted that the simulated vertical profile in Fig. 4b matches closely with the reference work (Blaise et al., 2016, Fig. 16b). In general, our 3D MCV nonhydrostatic model simulates the linear hydrostatic mountain waves quite well.


3
4.2.3. Stratified flow past a steep 3D hill
In this test case the simulations of the complex flow behind a steep isolated hill are challenging for the newly developed numerical model. The steep isolated hill is defined bywhere the peak height of the hill h0 = 1500 m, the half-width L = 3000 m, and r is the distance to the mountain center (x0, y0). A plot of the hill profile is shown in Fig. 5. The domain size is [0, 15 000] m × [?6000, 6000] m × [0, 6000] m, respectively, in the x, y and ζ direction. The mesh consists of 151 × 121 × 61 grid points with uniform grid spacing ?x = ?y = ?ζ = 100 m. The model is integrated up to 1200 s when the main features of the solution are already established. The sponge layer extends upward from a height of 4500 m and at the lateral outflow boundary with 1500 m thickness. Following Smolarkiewicz et al. (2013), the initial state is characterized by uniform horizontal wind (u = 5 m s?1) with constant buoyancy frequency (N = 10?2 s?1). The reference potential temperature is θ0 = 300 K. The specified environmental profile and the height of the hill result in a low Froude number flow. Here, the Froude number equals 1/3. The flow is in the nonhydrostatic range. This test case also complements the quantitative test of our 3D nonhydrostatic model. For convenience of comparison, the same contour intervals in the figures of this test case are used as in previous studies.

Figure 6 shows the vertical velocity along the lower surface of z = h(x, y). Wind vectors are overlaid on the contours of vertical velocity. It can be seen in Fig. 6 that the symmetric structures of both vertical velocity and wind vectors are perfectly reproduced, similar to previous studies (Smolarkiewicz et al., 2013), especially the two large eddies on the lee side and the lower flow splitting.

Figure 7a displays the vertical velocity in the central xz cross section (y = 0). This figure is generated by linearly interpolating the vertical velocity in the transformed coordinate ζ to the geometric height coordinate z. It can be seen that the MCV nonhydrostatic model captures fine-scale features of the turbulent wake on the lee side and the gravity-wave response aloft. The magnitude of vertical velocity also agrees quantitatively with the simulation of Smolarkiewicz et al. (2013). In comparison with the reference results (Smolarkiewicz et al., 2013; Szmelter et al., 2015), our model introduces small spurious noises, although numerical dissipation is not imposed. Additionally, our numerical results show a better simulation of the positive vortex on the lee side. Therefore, we can again confirm the competitive performance of our model.

The profiles of vertical velocity in the xy cross section at the height of z (= h0/3) = 500 m are plotted in Fig. 7b. This horizontal flow patterns also compare well with existing results [see Fig. 6 in Smolarkiewicz et al. (2013) and Fig. 7 in Szmelter et al. (2015)]. The fine-scale eddy structures are perfectly resolved by our model, especially the two evident eddies behind the hill. Figure 8 gives the xy cross section of vertical velocity at z = 2500 m. No visually distinguishable differences are observed between our results and those of Smolarkiewicz et al. (2013) (see their Fig. 9). Overall, it can be seen that, without any artificial dissipation terms included in the equation set, the 3D MCV nonhydrostatic model not only simulates a smoother numerical solution, but correctly reproduces all the small-scale characteristics generated by the steep hill.

The numerical results of various benchmark tests indicate that the present MCV nonhydrostatic model is very competitive when compared to most other existing advanced models. The mountain wave tests verify the capability of the present model in accurately resolving the small-scale characteristics generated by terrain effects, especially when the topographic inclination is steep. The present model is a practical and promising framework to simulate nonhydrostatic atmospheric flow.
The implementation of more complicated physical processes in the model is still in progress, and our preliminary tests with a simple wet physical process look very promising. The related results will be reported in another paper in the near future. Meanwhile, the extension of the present MCV nonhydrostatic dynamics into the spherical system is another step to be undertaken in future studies.
Acknowledgments. This work was supported by the National Key Research and Development Program of China (Grant Nos. 2017YFC1501901 and 2017YFA0603901) and the Beijing Natural Science Foundation (Grant No. JQ18001).