An MCV Nonhydrostatic Atmospheric Model with Height-Based Terrain following Coordinate : Tests of Waves over Steep Mountains

1National Meteorological Center, China Meteorological Administration, Beijing 10086, China 2Center of Numerical Weather Prediction, China Meteorological Administration, Beijing 10086, China 3State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, Beijing 100086, China 4Department of Energy Sciences, Tokyo Institute of Technology, Yokohama 226-8502, Japan 5School of Human Settlement and Civil Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China


Introduction
Mountain weather processes, such as lee waves, rotors, and downslope windstorms, which have great influence on the air quality over complex mountainous terrains [1,2], involve a wide range scales of air motions and present a challenge to numerical modeling.With the rapid development of computer hardware, it is now possible for the atmospheric models to represent the complexity in topography with the increasing horizontal resolutions and thus provide more accurate predictions for the mountain weather processes.The ability to simulate the atmospheric flows over complex mountainous areas becomes highly demanded for the nonhydrostatic atmospheric numerical models, such as MM5 [3], MC2 [4], COAMPS [5], LM [6], ARPS [7], WRF [8], and GRAPES [9].Though the significant advancements have been achieved during the past several decades, adequate simulations and predictions of the complex terrain-forced weather processes, for example, mountain waves, still remain an issue unsatisfactorily resolved.Further efforts are still required to develop more reliable dynamic cores with accurate representations of the topographic effects to improve the simulation of atmospheric flows over complex terrains.
Recently, a new nonhydrostatic model was developed by using the multimoment constrained finite volume (MCV) method [10].Different from the conventional finite volume method, the unknowns (or the degree of freedom (DOFs)) are defined as the values at the solution points distributed within each mesh cell.In contrast to the direct multimoment method [11][12][13][14][15][16][17][18] where the moments are directly used as the predicted unknowns, the MCV formulation [19] updates the DOFs as the point values at the solution points whose time evolution equations are derived by applying a set of constrained conditions imposing on different types of 2 Advances in Meteorology moments, such as volume-integrated average (VIA), point value (PV), and spatial derivative values (DV) and is thus simple, efficient, and easy to implement.Being a high order scheme, more accurate numerical results can be obtained in terms of the equivalent DOF resolution in comparison with the traditional finite volume method even with relatively coarse grid spacing.The rigorous numerical conservation in MCV model is exactly guaranteed by a constraint on the VIA through a finite volume formulation of flux form.
Being a new nodal-type high order conservative method, it is much beneficial to compute the metric and source terms in an MCV model, which are always involved in the treatments of spherical geometry in the horizontal direction and coordinate transformation in the vertical direction for topographic effect.The MCV model has some appealing features for atmospheric modeling, such as the rigorous numerical conservation, good computational efficiency, and flexible configuration for solution points, and thus can be expected as a practical framework of the dynamic cores, which well balances the numerical accuracy and algorithmic complexity [10,20].The competitive results of the widely adopted benchmark tests can be referred to [10,[20][21][22][23][24].
To deal with the bottom topography, the height-based terrain following coordinate is adopted in our MCV model.Since Phillips' pioneering work [25], the terrain following coordinates, mainly classified into the pressure-based coordinate and the height-based coordinate, have been widely adopted to represent the underlying mountainous surface in atmospheric models.During the past several decades, the height-based coordinate [26] has got an increasing popularity mainly due to its applicability for both hydrostatic and nonhydrostatic models and computational simplicity.A variational grid generation technique [27] was adapted to mountain wave simulation as well.Modified version of terrain following coordinate has been also proposed [28] to circumvent to some extent the drawbacks of the coordinate in representing steep topography.In this study, a set of benchmark tests of mountain waves generated by a constant background flow over mountains with different steepness, which essentially represent the complex mechanisms involved in mountain weather processes, are examined by the MCV nonhydrostatic model using the height-based terrain following coordinate in order to verify the performance of the MCV model in simulating the mountain weather processes and its potential for the further numerical investigations on the flows in the atmosphere boundary layer and the air quality over complex mountainous terrain.The numerical results are evaluated in comparison with the semianalytical solutions obtained from the linear theory [29], as well as the numerical solutions from other representative models.
The remainder of this paper is organized as follows.In Section 2 the compressible nonhydrostatic atmospheric model using the MCV scheme and the height-based terrain following coordinate is briefly introduced.Section 3 describes the mountain tests and the semianalytical solutions through the linear theory.Section 4 discusses numerical results of various mountain waves by the MCV model.Finally, a short conclusion is given in Section 5.

2D MCV Nonhydrostatic Atmospheric Model
In order to consider complex topography as the bottom boundary of the atmospheric model, a height-based terrain following coordinate is used in this study to map the physical domain (, ) to the computational domain (, ) through the transformation where   () is the elevation of topography,   the altitude of the model top, and  ∈ [0,   ].
Using a height-based terrain following coordinate, 2D compressible and nonhydrostatic governing equations for atmospheric dynamics are written in flux form as where where  is density, (, ) are velocity vector in the physical domain, w = / is the vertical velocity in the transformed coordinates,  13 = /, and √  = / is the Jacobian of transformation.
The thermodynamic variables are split into a reference state and the deviations to improve the accuracy of the numerical model as  (x, ) =  () +   (x, ) ,  (x, ) =  () +   (x, ) , where the reference pressure () and the density () satisfy the local hydrostatic balance, x is the position vector,   =  0 ()  , and  0 =  0 () −1 .The constants used in the simulations are specified as follows.Gravitational acceleration  = 9.80616 ms −2 , ideal gas constant for dry air   = 287 Jkg −1 K −1 , specific heat at constant pressure   = 1004.5Jkg −1 K −1 , specific heat at constant volume  V = 717.5 Jkg −1 K −1 ,  =   / V = 1.4,reference pressure at the surface  0 = 10 5 Pa, and constant  0 =    −  / V 0 .The MCV scheme is adopted in this model to solve the governing equations (2).The MCV scheme is a general numerical framework for developing high order numerical models to solve the hyperbolic systems.A major feature, which distinguishes MCV scheme from other conventional numerical schemes, is the local high order spatial reconstruction.For the sake of brevity, we omit the details of the numerical formulation of MCV nonhydrostatic model in the present paper.The fourth-order MCV scheme and the 3rd TVD Runge-Kutta time scheme [30] are adopted in this study.A local Lax-Friedrich approximate Riemann solver [31] is used for computational efficiency.The interested readers are referred to [10] for details.
The boundary conditions on the bottom surface are of crucial importance in atmospheric models, especially for simulations of the waves generated by complex topography.For the test cases studied in this paper, no-flux condition is imposed along the bottom boundary and nonreflecting condition is used for the lateral and the top boundaries.
The no-flux boundary condition requires the velocity field to satisfy the relation where n is the outward unit normal vector of the bottom surface and u = (, )  is velocity vector on the bottom boundary.
The nonreflecting boundary conditions are realized by a sponge layer along the lateral and top boundaries that relaxes the numerical solution to the prescribed reference.The damping terms are added to the momentum and potential temperature equations as where  is the relaxation coefficients and q  is the specified reference state.More details about the strength of the Rayleigh damping  and the formulations of semidiscrete noflux condition are described in [10].

Mountain Wave Test Cases
Satomura et al. [32] suggested a set of mountain wave tests to evaluate the capability of atmospheric models in reproducing topographic effects.The isolated bell-shaped bottom mountain to trigger the waves is specified as where ℎ 0 is the maximum height of the mountain,  0 is the center of physical domain, and  is the half width of the mountain.
Six test cases with different mountain height and horizontal width adopted in this study are given in Table 1.In these cases, the mountain slopes (measured by averaged inclination angles ) vary from 0.006 to 45 degrees.The cases of large inclination angles indicate steep mountains and are very challenging for the height-based terrain following coordinate.The initial hydrostatic conditions in these tests are specified in terms of Exner pressure Π = (/ 0 )   /  and potential temperature  via the hydrostatic relation Π/ = −/(  ).Setup of each case will be described in detail in Section 4.
The semianalytical solutions to mountain waves are derived from the linear theory.We first briefly introduce the linear theory before discussing the numerical results.
Using the linear theory [29], the small-amplitude 2D mountain waves can be described by a single partial differential equation as where  is the vertical displacement of a parcel at a steady state and  =  0 / is Scorer parameter which is constant for an elastic, isothermal, and nonshear uniform flow. 0 is Brunt-Väisälä frequency.For a bell-shaped mountain studied in this paper, the analytic solution of (8) can be obtained by using the Fourier transform method.The vertical displacement of a parcel (, ) is obtained by where  is the wave number, h() is Fourier transform of the mountain shape,  = √ −1, and Re{ } indicates the real part of a complex number.From (9), it is found that the solution for (, ) is proportional to the maximum height of mountain ℎ 0 since h() = ℎ 0  exp(−).The flow pattern of the mountain wave depends on the dimensionless quantity .When  ≪ 1, the wave disappears and a simple uniform flow is retrieved.Forced by the lower topographic surface, the flow climbs up and goes down along the mountain and the highest speed and lowest pressure occur at the top of the mountain.When  ≫ 1, the buoyancy-dominated hydrostatic flow over an isolated ridge develops and the disturbance energy is propagating upward away from the mountain.When  = 1, the flow is in the nonhydrostatic regime, and the nonhydrostatic mountain waves will appear.It is also necessary to consider another nondimensional parameter ℎ 0 , that is, mountain height scaled by Scorer parameter, which suggests the applicability of linear theory.Generally, the small value of ℎ 0  (ℎ 0  ≪ 1) in test cases 1 to 3 meets the linear theory assumption since the mountain is relatively low in comparison with the vertical stratification.While ℎ 0  ≈ 1 in test case 4, the nonlinear effects have to be considered for a high mountain.
The vertical flux of horizontal momentum is a measure of the performance of a numerical model to simulate the wave energy.Following Eliassen and Palm [33] and Smith [29], the vertical flux of horizontal momentum (called momentum flux hereafter) is defined as where () is the reference density and   and   are the wind disturbances from the basic state.() is constant for linear mountain waves in a constant mean wind [33].Furthermore, when the linear mountain waves are irrotational and hydrostatic, the analytic momentum flux is given as where  0 and  0 are the reference density and horizontal velocity values at the ground level.hours.The MCV model can accurately capture the mountain waves triggered by the gentle slope, except the numerical maxima and minima contours in the vertical velocity field are slightly weaker than the analytic ones.Compared with the existing numerical results obtained by high order schemes such as spectral element (SE) and discontinuous Galerkin (DG), the present results look quite similar to those in [34].

Numerical Results
The momentum fluxes at different instants are plotted in Figure 1(c).Values normalized by   given in (11) are shown.As the linear hydrostatic mountain waves exist in this test, the vertical variation of the normalized momentum flux is usually used to examine the numerical dissipation of a model.As shown in Figure 1(c), as the mountain wave develops into a mature state, the momentum flux gradually approaches a vertically uniform distribution with a value close to unity of the analytical solution.At  = 10 (i.e., nondimensional time / = 72), the numerical flux is about 0.99 at the surface and approaches 0.9645 at the height of double vertical wavelength ( = 2/ ≈ 6.4 km).Present results are competitive to those shown in [7,35], where the flux at  = 6.4 km is 94% of its steady value at nondimensional time of 60 in [35] and the flux reaches 0.96 at later time at a height just below their Rayleigh damping layer (12 km) in [7].The results reveal that the high order accuracy of MCV model is much beneficial to improve the simulation of nondissipative hydrostatic mountain waves.

Case 2:
Linear Nonhydrostatic Mountain Waves.The initial state of the atmosphere is specified as a flow of uniform horizontal velocity  = 10 ms −1 with Brunt-Väisälä frequency  0 = 0.01 s −1 .The reference potential temperature and Exner pressure are respectively, with  0 = 280 K.The computational domain is [0, 144000] m × [0, 30000] m with a resolution of Δ = 360 m and Δ = 300 m.The model runs for 5 hours.The flow is in the nonhydrostatic regime since the dimensionless parameter  = 1.The averaged inclination angle is about 0.06 degrees, which is larger than that of case 1.The nonreflecting boundary conditions in this case are imposed by placing the damping term in the computational domain  ≥ 15 km for the top boundary and the thickness of 30 km for the lateral boundary.Figures 2(a) and 2(b) show the numerical and the semianalytical solutions of vertical velocity  after 5 hours (nondimensional time of 180) for linear nonhydrostatic mountain waves.It is observed that the linear nonhydrostatic mountain waves are distinguished from the linear hydrostatic ones by the dispersive character of wave trains behind the mountain peak.The simulated vertical velocity agrees well with the analytical solution and the numerical result of other existing high order schemes, for example, DG3 model in [34].

Case 3: Simple
Flows.Two types of waves are checked in case 3, which are denoted as A3 and A4 cases following Satomura et al. [32].In these cases, the wave structures diminish, and the topographic forcing is limited only to the lower atmosphere, which is referred to as the simple flow pattern in this paper.The reference constant wind at the initial time is  = 10 ms −1 with Brunt-Väisälä frequency  0 = 0.02 s −1 for both cases.The reference potential temperature The model runs for 20 minutes (nondimensional time of 120) until it reaches the steady state.Compared with cases 1 and 2, the mountain slope is steeper, that is, roughly 26.5 degrees in regard to averaged inclination angles.The damping term is imposed in the computational domain  ≥ 4 km for the top boundary and the thickness of 3 km for the lateral boundary to assure the nonreflecting boundary conditions.
Figures 3(a) and 3(b) depict the numerical and the semianalytical solutions of vertical velocity  at nondimensional time of 120.The numerical results by MCV model are visually identical to those by the linear theory.In case of the equivalent DOF resolution, that is, Δ DOF = Δ/3 = 20 m, the present results agree well with those in Satomura et al. [32].The dimensionless parameter  = 0.2 indicates that a simple irrotational flow will appear in the simulation.It is manifested by the streamlines shown in Figure 3(c), where the steady flows climb up to the peak of mountain and then go down, and the largest speed and the lowest pressure are observed at the top of the mountain.The momentum flux scaled by   given in ( 11) is presented in Figure 3(d).It should be Advances in Meteorology pointed out here that we use ℎ 0 = 27 m, instead of the real mountain height, to calculate   to make the normalized values of momentum flux distributed around unity for better illustration.The values of normalized flux in our model remain between 0.2 and 1.2.Similar treatment is also adopted in following test cases.The profile below the nonreflecting absorbing layer at 4 km shows small fluctuations, which are also observed in [32,  The model is integrated for 10 minutes (nondimensional time of 120).The nonreflecting boundary conditions in this case are imposed by placing damping term in the computational domain  ≥ 1 km for the top boundary and | − 1.5 km| ≥ 1.2 km for the lateral boundary.
Figures 4(a) and 4(b) show that the numerical and the semianalytical solutions of vertical velocity are at nondimensional time of 120, which agree quite well.Our results also agree well with those in [32] when using the equivalent DOF resolution; that is, Δ DOF = 5m.The concavity of vertical velocity near the ground is observed in Advances in Meteorology the present results, which also appears in MRI/NPD-NHM model with horizontally explicit and vertically implicit (HE-VI) method and TSO model in [32].Compared with A3 case, the dimensionless parameter  is smaller and equals 0.1.As a result, the wave structure predicted by the linear theory is the simple pattern flow and confirmed by streamlines shown in Figure 4(c).
The mountain slope in this case is steeper than the A3 case due to a narrower half mountain width , and the averaged inclination angle is about 45 degrees.This case verifies the performance of the proposed numerical model to simulate the waves generated by steep mountains.Furthermore, the normalized flux profile scaled by formulation (11) with ℎ 0 = 21 m is shown in Figure 4(d).It is observed that the flux profile below the nonreflecting absorbing layer at 1 km is nearly unity, which is better than those in Satomura et al. [32].

Case 4: Nonlinear Mountain
Waves.Similar to case 3, two types of mountain waves are tested in case 4 and denoted by D1 and D2 cases in Satomura et al. [32].The initial constant mean flow is  = 10 ms −1 with Brunt-Väisälä frequency  0 = 0.01 s −1 for both cases.The reference potential temperature and the Exner pressure are computed the same as in case 2. The model runs for 100 minutes.In case 4 the parameter ℎ 0  ≈ 1, and the mountain is specified to be high enough to take the nonlinear effects into account [29], which is different from other test cases in this study.Compared with the results shown in [32], our results agree well with those of TSO model when using the equivalent DOF resolution; that is, Δ DOF = 50 m.It is noticed that there is slight difference in the amplitudes between our result and the analytical solution.It may owe to the nonlinear effects of the high mountain.Different from case 3, the dimensionless parameter  equals 0.5 so that the waves will propagate upward.The streamlines are plotted in Figure 5(c) and it is observed that the mountain waves propagate upward behind the mountain peak.The normalized flux profile scaled by formulation (11) with ℎ 0 = 250 m is presented in Figure 5(d).It can be seen that the flux profile below 12 km is nearly unity, whereas the flux fluctuation is visible under the height of 8 km.Although the mountain is high enough to take into account the nonlinear effects, the flux profile from the numerical model is close to that of the linear theory except the fluctuation structure.

D2 Case.
In this case all mountain parameters are the same as in the D1 case except the half width of mountain  = 250 m, which leads to the increase of the averaged inclination angle of mountain to about 45 degrees.The dimensionless parameter  = 0.25, and in this case the vertically propagating wave will not be noticeable.The boundary conditions are imposed the same as D1 case.
The numerical and semianalytical solutions of vertical velocity  at nondimensional time of 240 are plotted in Figures 6(a) and 6(b).It is seen that the wave propagates with the similar structure compared with that by the linear theory.However, the vertical velocities are larger than the analytical solution due to the significant nonlinear effects of the high mountain, which are also observed in [32].The dimensionless parameter  equals 0.25 and is smaller than D1 case, which means the upward propagating waves are relatively weak.This is confirmed by the numerical results shown in Figure 6(c) for the streamlines.The normalized flux profile scaled by formulation (11) with ℎ 0 = 180 m is shown in Figure 6(d).
The values of normalized flux by the MCV model are nearly unity and fluctuate with the height.The flux profile structure at nondimensional time of 240 is similar to those of MRI/NPD-NHM model (HE-VI) and TSO model in [32] when using the equivalent DOF resolution Δ DOF = 50 m.

Conclusions
The mountain waves triggered by a constant flow over a bell-shaped mountain are investigated by MCV compressible nonhydrostatic atmospheric model with the application of the height-based terrain following coordinate.A set of mountain wave test cases in which the mountain slopes range from 0.006 to 45 degrees are systematically checked.Comparison with the semianalytical solution from the linear theory reveals that the present model can accurately reproduce various mountain waves in the linear regime.The numerical results are also compared with other existing models, including high order models using SE and DG schemes as well as other operational and research models, which verifies the numerical accuracy of the present model in representing complex topography of different steepness.
Verified in this work, the MCV model can accurately capture various mountain-triggered waves by using the heightbased terrain following coordinate, which thus provides a tool for better understanding the mountainous weather processes and the related air pollution in the boundary layers over mountainous terrains.
Being nodal-type high order scheme, the MCV method has significant advantage in treating geometrical metric terms.It can be concluded from the present study that the present numerical framework, which adopts the MCV scheme and the height-based terrain following coordinate, is promising for further developing 3D nonhydrostatic atmospheric model.

4. 1 .
Case 1: Linear Hydrostatic Mountain Waves.The initial state is an isothermal atmosphere with the constant temperature of  = 250 K moving at a constant velocity  = 20 ms −1 .The Brunt-Väisälä frequency  0 = /√   in this case.As a result, the reference Exner pressure is Π = exp(−(/  )).The computational domain is [0, 240000] m × [0, 30000] m with the grid spacing of Δ = 1200 m and Δ = 240 m.The MCV model is integrated up to 10 hours.According to the linear theory, the parameter  determines the hydrostatically balanced mountain waves in this test.Figures1(a) and 1(b) show the numerical and the semianalytical solutions of vertical velocity component  at  = 10

Figure 1 :
Figure 1: Numerical results of linear hydrostatic mountain waves.Numerical solution (a), semianalytical solution (b), and normalized momentum flux (c) are shown.The contour values vary from −0.005 to 0.005 with an interval of 0.0005.Negative values are denoted by dashed lines.

Figure 2 :
Figure 2: Same as Figure 1, but for linear nonhydrostatic case.

Figure 3 :
Figure 3: Numerical results of case A3.Numerical solution (a), semianalytic solution (b), streamlines from 0 to half vertical wavelength height over an isolated mountain (c), and normalized momentum flux (d) are shown.Contour interval is 0.25 ms −1 .All negative values are denoted by dashed lines.

4. 4 . 1 .
D1 Case.The computational domain is [0, 18000] m × [0, 30000] m with a grid spacing of Δ = 150 m and Δ = 150 m.Similar with A3 case, the averaged inclination angle is about 26.5 degrees.The damping term is placed in the computational domain  ≥ 10 km for the top boundary and || ≥ 10 km for the lateral boundary.Figures 5(a) and 5(b) give the numerical and the semianalytical solutions of vertical velocity  at nondimensional time of 120.The pattern of vertical velocity simulated by MCV model compares quite well with that of the linear theory.

Figure 6 :
Figure 6: Same as Figure 3, but for Case D2.Contour interval is 0.5 ms −1 in this case.

Table 1 :
Configuration of mountains in different test cases.