Development of a non-hydrostatic atmospheric model using the Chimera grid method for a steep terrain

In a high-resolution atmospheric model, the terrain is resolved in detail, and thus steeper and more complex terrain can be represented. Thus, we developed an atmospheric model that represents the terrain by using the Chimera grid method, which represents the computational region as a composite of overlapping grids. To test this model, we simulated a lee wave over a semicircular mountain and one over a tall semi-elliptical mountain. The results show that the Chimera grid method can be used to represent a very steep terrain which terrain-following coordinates and numerically generated coordinates cannot represent appropriately


Introduction
Recent developments in computing have rapidly increased atmospheric models' resolution.In high-resolution models, the terrain is resolved in more detail, and thus steeper and more complex terrain can be resolved.Currently, terrain-following coordinates are usually used to represent the terrain in atmospheric models (Gal-Chen and Somerville, 1975;Clark, 1977).However, over steep and complex terrain, such coordinates are greatly deformed, and thus they induce serious errors (Mesinger and Janjić, 1985).Therefore, a good method for the representation of terrain is one of the challenges of high-resolution atmospheric models.
Previous research has handled this challenge by using one of two approaches.One approach is to use Cartesian coordinates and the other one is to use boundary-fitted coordinates.Cartesian coordinates do not require a transformation.Thus, the coordinates are not deformed, even when the terrain is steep and complex.Moreover, the governing equations can be represented in a simpler form.However, the terrain features are unlikely to coincide with the grid lines of the Cartesian coordinates.Thus, various ways have been proposed for representing the terrain: e.g. a cut-cell method (Adcroft et al., 1997;Yamazaki and Satomura, 2010) and a thin-wall approximation (Steppeler et al., 2002).These methods successfully simulated flow over steep slopes.Conversely, Cartesian coordinates require an additional scheme to efficiently increase the resolution near the ground level (Yamazaki and Satomura, 2012).
On the other hand, boundary-fitted coordinates are transformed following the topography.As the topography boundary coincides with the lower boundary of the model, it is easy to apply boundary conditions.Moreover, it is easy to increase the vertical resolution near the ground.However, with this coordinate system, the truncation error increases due to the deformation of the coordinates over steep terrain (Mesinger and Janjić, 1985).Klemp (2011) proposed smoothed terrain-following coordinates obtained by modifying the transformation to reduce these deformation at higher elevations.These transformed coordinates succeeded in reducing the error, but there are limitations near the ground as it is only a vertical modification.Satomura (1989) used numerically generated coordinates, which have high orthgonality because a vertical transformation and a horizontal transformation are used.This coordinate system successfully simulated flow over a semicircular mountain that was so steep that terrain-following coordinates induced serious errors.However, this system cannot be generated over complex terrain in which the slope angle changes abruptly: e.g. a cliff.More improvements are necessary if such complex terrain is to be represented adequately.
In the field of computational fluid dynamics, complex boundaries such as those of an airplane wing are treated, and the Chimera grid method is commonly used (Benek et al., 1986).This method is also referred to overlapping grids, overset grids or composite overlaid grids.In this method, the computational region is represented by a composite of overlapping grids, and one grid is used to represent the entire computational region.Other grids are used to represent local boundaries that are too complex to be represented by the main grid.In the field of atmospheric sciences, the Chimera grid is used for global models to represent the entire Earth (Dudhia and Bresch, 2002;Peng et al., 2006;Staniforth and Thuburn, 2012).However, this method has not yet been used to represent complex terrain.
In this article, we present an atmospheric model that uses the Chimera grid to represent steep and complex terrain.To assess this method, we simulated flow over a semicircular mountain and over a tall semi-elliptical mountain.

Model description
We consider the two-dimensional (x: horizontal, z: vertical) full-compressible atmosphere for a resolution from tens of meters to hundreds of meters.The governing equations of this model were transformed from Cartesian coordinates (x, z) = (x 1 , x 2 ) to general coordinates (, ) = ( 1 ,  2 ).The transformed governing equations are (1) where (U 1 , U 2 ) are the (x 1 , x 2 ) components of the flow velocity and (v 1 , v 2 ) are the contravariant component of the flow velocity in the ( 1 ,  2 ) coordinate, respectively, p is the total pressure,  is the total density,  is the total potential temperature, (g 1 , g 2 ) = (0, − g) is the gravitational acceleration, R d is the gas constant, p 0 is a reference pressure of 10 5 Pa, C p is the specific heat at constant pressure, 1 is the Jacobian.Note that some of these variables can be divided into basic components p (z) , p (z) and  (z), and perturbation components p ' (x, z, t),  ' (x, z, t) and  ' (x, z, t), respectively.Basic components p (z), p (z), and  (z) fulfill: The last terms in the right-hand side of Equations ( 1) and (3), diff .U i and diff . are source terms due to mixing and diffusion, for which we employed the first-order gradient transport theory as: where K m and K h are the coefficients of mixing and diffusion for the momentum and the heat, respectively.The value of K m is constant value of 0.1 m 2 /s and K h = 3K m (Deardorff, 1972).To damp upward propagating waves before they are reflected from the upper boundary, we added a Rayleigh damping term − (z)U i to the right-hand side of the Equation ( 1).The value of  is determined as following.
where D is a height scale, which is set to 15 km.In addition, to simulate flows stably by removing high-frequency waves, we also added fourth-order diffusion terms in the  1 and  2 directions to the right hand side of the Equations ( 1)-( 3).At points next to a boundary point, the fourth-order diffusion terms were replaced by the second-order diffusion terms.We adopted the second-order central difference for the space discretization and the fourth-order Runge-Kutta scheme for time integration.For generating the grid, we used the variational method proposed by Brackbill and Saltzman (1982).As shown in Figure 1(a), scalar variables and velocity components were arranged on the centers of cells and the corners of cells, respectively, which was used in Yamazaki and Satomura (2010).

The Chimera grid method
In the Chimera grid, the computational region is represented by a composite of overlapping grids.In this article, we refer to the grid used to represent the entire computational region as the global grid, and we refer to a grid used to represent a detail of the terrain as a local grid.These grids are partially overlapping around each grid boundary, as shown in Figure 1(b).At each time step, the governing equations are solved at each grid point, and then the boundary point u I is obtained by interpolating between surrounding grid points u ij , using them with a boundary condition: where S ij is interpolation coefficient and N is the order of the interpolation stencil.We used a first-order bilinear interpolation with stencil order N = 2 and a third-order Lagrange interpolation with stencil order N = 4.However, we used the bilinear interpolation instead of the Lagrange interpolation near the ground where the number of points used for interpolation was not sufficient.The interpolation was not performed in the physical space (x, z) but in the computational space (, ), as the physical space coordinates were deformed, and this caused the interpolation to be complex.In the computational space (, ), the grid points were ordered on a rectangle.The interpolation coefficients S ij were Local grid Global grid represented as a product of the interpolation coefficients s i , s j for each direction: where ( I ,  I ), ( i ,  i ) are the coordinates of the interpolation point u I and the surrounding grid points u ij , respectively.The coordinates ( I ,  I ) are needed to evaluate the interpolation coefficients S ij , but, in general, they are unknown.Sherer and Scott (2005) proposed high-order offsets, and with these, ( I ,  I ) are evaluated by using the two-dimensional Newton method: In Equation ( 11), (F x , F z ) is the error of the interpolated physical space coordinates, defined as the difference between ( ) and the true physical space coordinates (x I , z I ).Iteration continues until the error (F x , F z ) is smaller than a predetermined tolerance.The grid interface was set to fulfill the following condition.The ratio of the size of an interpolation point cell to that of the surrounding grid points should be close to unity.The area of the overlapping region must be large enough for each interpolation point to have N stencil points but should be as small as possible.
When introducing the Chimera grid to an atmospheric model, it is necessary to consider the stratification of the atmosphere.As thermodynamic variables change vertically in a stratified atmosphere, vertical interpolation may induce errors.Therefore, only the perturbation components were interpolated.

Simulation and results
To test the Chimera grid, we simulated a lee wave over a semicircular mountain and one over a tall semi-elliptical mountain.At the foot of each of these mountains, the slope angle changes rapidly and it becomes very steep.Thus, terrain-following coordinates induce serious errors (Satomura, 1989).We compared the numerical solutions with the analytical solution for Bousinessq system by Miles (1968).To mimic Bousinessq system within our fully compressible system model, we set the value of surface potential temperature  0 at 4000 K, where the sound speed increased about four times of that at 300 K.

Semicircular mountain
We simulated a lee wave over a semicircular mountain of height H = 1000 m.For the initial condition, we assumed that the atmosphere had a constant buoyant frequency of N = 0.01 s − 1 , and a constant horizontal velocity of U = 10 m s − 1 .The horizontal and vertical grid intervals were Δx = Δz = 250 m, respectively.The domain consisted of 2000 horizontal and 120 vertical cells, respectively.We performed the simulation using the terrain-following coordinates, the numerically generated coordinate, and the Chimera grid. Figure 2 shows the grid for each coordinate system.In the Chimera grid, the global grid uses Cartesian coordinates, and around the mountain, the local grid uses polar coordinates.The terrain-following coordinates were greatly deformed above the mountain, but there was less deformation with the numerically generated coordinates and the Chimera grid. Figure 3 shows the analytical solution of the vertical velocity over the semicircular mountain and the differences of the numerical solutions at t = 10  We performed an quantitative comparison of the Chimera grid and the numerically generated coordinates, and we evaluated the accuracy of our proposed model over the semicircular mountain.We used Equations ( 14) and ( 15) to compute the horizontal  velocity error  U and the vertical velocity error  W .
where N x and N z are the numbers of grid points in the x and z directions, respectively; and U c i,k and W c i,k are the analytical solution of the horizontal and vertical velocities, respectively.The domain for the error evaluation was a rectangular domain of − 3 km ≤ x ≤ 3 km and 0 km ≤ z ≤ 3 km.Three values of the grid intervals (Δx, Δz) were tested, 400 m, 200 m and 100 m. Figure 4 shows how the error changed with changes in the grid interval.For each coordinate system, the order of the convergence is less than second-order although we adopted the second-order finite difference scheme.This may be partially due to the difference between the compressible system and the Bousinessq system.However, the error seems to be converging for each coordinate system.Moreover, for each grid interval, the error for the Chimera grid is almost equal to or smaller than that for the numerically generated coordinates.Thus, it can be concluded that the Chimera grid yields an equally or more accurate numerical solution than the numerically generated coordinates in the semicircular mountain case.This is likely due to the fact that the grid of the Chimera grid was less deformed than that of the numerically generated coordinates at the foot of the mountain.Moreover, a comparison of the results of the bilinear interpolation with those of the Lagrange interpolation shows that the error of the bilinear interpolation is almost equal to or smaller than that of the Lagrange interpolation.Thus, it can also be concluded that the decrease in accuracy is negligibly small when using the first-order bilinear interpolation.

Tall semi-elliptical mountain
Finally, we simulated a lee wave over a tall semi-elliptical mountain.As the slope changes sharply at the foot of the mountain, numerically generated coordinates cannot properly simulate this topography.Huppert and Miles (1969) obtained the analytical solution of a lee wave over a tall semi-elliptical mountain.They concluded that the behavior of such a wave depends on the ratio of the mountain's height H to its width B, on  = H/B, and on the nondimensional height  = NH/U, which is the reciprocal Froude number.When  >  c , the amplitude of a lee wave develops so large that the  field becomes convectively unstable.When  ≤  c , a stable lee wave develops and converges to a stationary solution.Here, the critical non-dimensional height  c monotonically increases with .We simulated a lee wave for this topography with a fixed ratio of  = 5, and with an average slope of 80 ∘ .According to Huppert and Miles (1969), when  = 5,  c is around 1.6.For the initial conditions, we assumed the atmosphere had a horizontal velocity of U = 10 ms − 1 , and a buoyant frequency of N = 0.01 s − 1 .The horizontal and vertical grid intervals were Δx = Δz = 150 m.We performed the simulation for three types of mountain: .Figure 5 shows the potential temperature field after 60 min of simulation for each mountain type in using the Chimera grid with the bilinear interpolation.In Figure 5(a) and (b), we can see that stable lee wave simulations are simulated when  ≤  c , and in Figure 5(c), we can see that a lee wave develops too largely and becomes unstable around 4 km above the mountain when  >  c .These results are consistent with the analysis of a lee wave over similar topography by Huppert and Miles (1969).It is evident that the Chimera grid can appropriately simulate time evolutions of flow fields complex terrain.

Summary
We have developed an atmospheric model using the Chimera grid method and have performed a high-resolution simulation of air flow over steep and complex terrain.The results show that the Chimera grid can reduce the errors that are produced by using terrain-following coordinates for steep terrain, and can simulate the flow appropriately over a very steep mountain for which the numerical coordinates cannot be generated.The Chimera grid would be useful for the simulation over such a very steep terrain, e.g. a cliff.However, the criterion for deciding which coordinate system should be used is not clear now.Moreover, the accuracy of the results highly depends on the grid configuration and shape of terrain.As future work, it is necessary to clarify the dependence of the accuracy of the numerical solutions on the grid configuration and the shape of topography in more detail.
The interpolation method used in this article did not consider the global conservation of physical quantities.To satisfy the global conservation, we should introduce a conservative interpolation, which interpolates fluxes of physical quantities around interpolation points (Chesshire and Henshaw, 1994).A conservative interpolation method was also developed for a global model (Peng et al., 2006), and would be able to be introduced into our model without extra difficulty.As for the extension of our model to three dimension, it will not be difficult as three dimensional grid generation was developed in the field of the computational fluid dynamics (Petersson, 1999).

Figure 1 .
Figure 1.(a) Arrangement of variables.Open circles and open squares represent scalar points and velocity points, respectively.(b) Grid interaction.Overlapping area of the global grid and the local grid.Point closed square and closed circle are the boundary point of the global grid and the local grid, respectively.These points are interpolated from around point of the other gird open circle and open square.

Figure 2 .
Figure 2. Representation of a semicircular mountain of height H = 1000 m: (a) the terrain-following coordinates; (b) the numerically generated coordinates; (c) Chimera grid with bilinear interpolation; and (d) Chimera grid with Lagrange interpolation.

Figure 3 .Figure 4 .
Figure 3. Vertical velocity field and the difference from the analytical solution after 10 h over a semicircular mountain.(a) analytical solution of Miles (1968) (Courtesy of S. Masuda) with contour interval of 1.0 m s −1 ; (b)-(e) the differences of numerical solution from the analytical solution with contour interval of 0.2 m s −1 .(b) terrain-following coordinates; (c) numerically generated coordinates; (d) Chimera grid with bilinear interpolation; and (e) Chimera grid with third-order Lagrange interpolation.The discontinuity of contour near the mountain in panel (d)-(e) is due to the interface of grids.