## Introduction

Numerical weather prediction (NWP) model uses various discretization method and grid structure (Haltiner and Williams, 1980; Gates, 1992; Klemp, 2011; Kang and Cheong, 2017). As the computing resources to run the model becomes powerful, the spatial resolution of the NWP model has been increased steadily (Tomita and Satoh, 2004). In general, the vertical resolution is much higher than the horizontal resolution because the vertical variation of atmospheric motion and state is much stronger than the horizontal due to a large horizontal-to-vertical ratio of atmospheric layer (Haltiner and Williams, 1980).

In traditional NWP models the grid-size of the horizontal grid is not uniform, as is characterized by the latitude-longitude grid system (Smagorinsky, 1963; Arakawa, 1972; Hoskins and Simmons, 1975; Arakawa and Lamb, 1977). Although the lat-lon grid system is still used in some models and operational NWP centers (Saito et al., 2006; Skamarock and Klemp, 2008), the quasi-uniform grid system has gained increasing popularity (Rivier et al., 2002; Tomita and Satoh, 2004; Nair et al., 2005; Zhang and Rancic, 2007; Marras et al., 2015; Choi and Hong, 2016; Kang and Cheong, 2017, 2018). In the vertical direction, the grid-size varies substantially with height, typically decreasing from the middle level towards lower and upper boundaries (Gal-Chen and Somerville, 1975; Haltiner and Williams, 1980; Simmons and Burridge, 1981; Saito et al., 2006). This is because, unlike the horizontal domain, the atmosphere is bounded by the Earth’s surface and upper boundary, and it is needed to resolve the atmospheric motion and state with finer grid near the boundaries where the flow speed changes steeply (Haltiner and Williams, 1980; Leuenberger et al., 2010; Klemp, 2011). Since there is no clear physical boundary in the upper atmosphere, not as the surface boundary, the NWP model usually sets a constant pressure level or height level as the upper boundary.

The nonuniform grid in the vertical direction is regarded as an essential component of the NWP model as is convinced in numerical weather predictions (Haltiner and Williams, 1980; Saito et al., 2006; Klemp 2011). However, the nonuniform grid, in comparison with the uniform grid, is known to provide less discretization accuracy for the governing equations of the atmospheric motions and state (Arakawa and Lamb, 1977; Durran, 1999). In spite of this inherent deficiency of the nonuniform grid, little attention has been paid to the systematic accuracy assessment so far.

In this study, the accuracy of the nonuniform grid for the two-dimensional transport equation is examined. For simplicity and clarity, the grid structure is given to be symmetric with respect to the middle point of the vertical domain. The nonuniform grid is established in such a way to maintain the ratio of larger-to-smaller grid-size of neighboring vertical grid-points to be the same. The paper is organized as follows. Next section will present the two-dimensional transport equation and boundary conditions. Section 3 shows the discretization method and theoretical error analysis with a focus on the nonuniformity of the grid-size. Initial conditions and analytic solution are provided in section 4. Results of numerical experiment and discussion in comparison with the theoretical analysis will be given in the subsequent section. Concluding remark is presented in the final section.

## Transport Equation and Boundary Conditions

Two dimensional transport (advection) equation for an arbitrary scalar variable *h* is written as:

where *u* and *w* are the advection velocity of *x* and *z* direction, respectively. Velocity field is given to be nondivergent, that is, $\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0$
. The nondivergence of the flow makes it feasible to introduce the streamfunction such that $u=\frac{\partial \psi}{\partial z}$
and $w=\frac{\partial \psi}{\partial x}$
.

Model domain is set as $0\le x\le L$ and $0\le z\le H$ , and the boundary condition of Eq. (1) is set as follows:

In the case where the scalar variable at boundaries is given initially zero, it remains as zero-valued at any time during time evolution. This is expressed as the following relationship:

## Discretization Method

Equation (1) is solved numerically for a given initial field of *h* and velocity field (*u, w*) using the third-order Runge Kutta scheme (Nair et al., 2005) and the FDM for the time and the space, respectively. The grid system used in this study is uniform in xdirection but nonuniform in the vertical direction, as shown in Fig. 1. Such a grid system is a typical one adopted in the numerical weather prediction models (Haltiner and Williams, 1980; Tomita and Satoh, 2004; Skamarock and Klemp, 2008). The grid interval in *z*-direction decreases monotonically from the midlevel toward the boundaries. In the present study, the ratio of two adjacent grid-intervals is kept at a constant value γ in the lower half domain and γ ^{−1} in the upper half domain, which is controlled as a model parameter. A typical value of γ in operational model is 1.0-1.2 (e.g., Saito et al., 2006; Haiden et al., 2017). Location of variables is depicted in Fig. 2 together with the vertical nonuniform grids. Non-staggered grid is used in the horizontal direction, while staggered grid is adopted in the vertical direction: All variables are located in the same vertical line, and the scalar variable and horizontal velocity are placed in the middle of vertical grid-points. By choosing this gridsystem, it is possible to focus on the effect of the vertical nonuniform grid on the transport equation.

Third-order Runge Kutta scheme is written as follows (Nair et al., 2005):

where the superscript *n* means time-level of timestepping procedure, that is, ${t}_{n}\left(=n\delta t\right)$
with δt being the time-step size, and *T* implies the time tendency of the equation.

Spatial discretization is performed with FDM: For the horizontal discretization the second-order centered FDM is used, while for the vertical discretization the staggered-grid 3-points stencil FDM is chosen. The discrete form of these methods can be derived based on the Taylor series. The procedure to derive the finite difference form related to these methods is explained briefly. Taylor series of an arbitrary function $f\left(z\right)$ at z + α is given as follows

where ${f}^{\left(n\right)}=\frac{{\partial}^{n}f}{\partial {z}^{n}}$ and $n!=n\cdot \left(n-1\right)\cdots 2\cdot 1$ . Let the function $f\left(z\right)$ be defined on discrete grid-points ${z}_{k}$ , which is denoted by $f\left({z}_{k}\right)$ . The grid interval is defined as ${\Delta}_{k}={z}_{k}-{z}_{k-1}$ . Then, the above equation is written with grid-point indices for an arbitrary grid-point ${z}_{p}$ as:

where the subscript *p* is an integer, ${\overline{\Delta}}_{p,k}={z}_{p}-{z}_{k}$
, and *R* means higher-order terms than the fourth.

Second-order FDM in *x*-direction is approximated with a centered difference as following

where (*i, k*) means the grid-index corresponding to (*x _{i}*,

*z*) and Ξ implies the gridsize in the

_{k}*x*-direction. Eq. (7) was obtained by neglecting the second-order or higher order terms from the Taylor expansions. Staggered-grid FDM on 3-points stencil for the vertical advection term is approximated in terms of an average of finite difference around the target point (Arakawa and Lamb, 1977; Haltiner and Williams, 1980; Tomita and Satoh, 2004; Cheong, 2006; Haiden et al., 2017):

As is expected from the Taylor series, this gives the second order accuracy for a uniform grid system and constant vertical velocity, while it is only *first-order**accuracy* for nonuniform grid system regardless of the specific value of vertical velocity. Since the discretization error for the nonuniform grid is first-order, this scheme invites numerical diffusion in the vertical direction depending on the sign of advection speed. For uniform grid and constant speed (*c*), the finite difference form of the vertical advection term is written as

When the above identity is applied to the advection equation, it is clear that the second-order derivative term acts as diffusion for *c* < 0 and upward increasing grid-size (i.e., the case of ${\Delta}_{k+1}>{\Delta}_{k-1}$
) and *c* > 0 and downward increasing grid-size (i.e., the case of ${\Delta}_{k+1}<{\Delta}_{k-1}$
). Negative diffusion is unavoidable in this scheme because three points including the target point at the center should be incorporated unlike the upwind differencing which uses only two points selectively.

In the case where the advection velocity is not a constant value, Eq. (9) is rewritten as

where

Equation (11) states that the coefficient of the discarded second-order term depends on the gridinterval and the advection velocity, which can be either positive (numerical adverse diffusion) or negative (numerical diffusion) as in the case of constant speed. This combined effect is expected to result in different behavior of the discretization error from the case of constant advection velocity.

## Initial Condition, Analytic Solution, and Model Parameters

In numerical modeling of the atmospheric circulation, analytic solutions are often used to assess the accuracy of discretization (Browning et al., 1989; Williamson et al., 1992; Cheong, 2006; Cheong and Park, 2007; Flyer and Wright, 2007). In this study, the velocity field is specified as a pair of circulation cells of different sign (cell L and cell R), of which streamfunction, *x*-direction and *z*-direction velocity components are given as

and

The flow pattern is very similar to the rotational flow on the sphere (Nair et al., 2005; Flyer and Wright, 2007; Cheong et al., 2015). The angular velocity of the advection flow is then written as

Scalar field is given in the form of a Gaussian bell, which is expressed as

where ${r}_{0}=3.0H$ and ${D}_{\theta}=\sqrt{{\left(0.3H\right)}^{2}-{\left(r-0.3H\right)}^{2}}+{10}^{-4}H$ , and ${\theta}_{c}\left(r,t\right)$ denotes the azimuthal location of the Gaussian bell. Due to the circularity of the advection flow, it can be written as ${\theta}_{c}\left(r,t\right)=\omega \left(r\right)t$ , which enables to specify the scalar field at any time (ω=angular speed).

The horizontal domain is set to be twice of the vertical domain so that two circular advection flows can be included in the model domain. The model equation is scaled is such a way that both the maximum advection flow and the vertical range are unity. In the control run, the number of the grid-point in *z*-direction is set *N _{z}*= 30 and that of

*x*-direction is set as

*N*= 2

_{x}*N*, with which the spatial resolutions in both directions are the same for the uniform grid spacing. Time integration is carried out until the maximum advection flow transports the scalar field along azimuthal direction by the nondimensional distance of 0.3 at the radial distance of

_{z}*r*= 0.25

*H*. The time-step size is given large enough to make the CFL number be unity for the smallest grid-size, which is smaller than the marginal value ($\sqrt{3}$ ) for the stability (Nair et al., 2005; Kang and Cheong, 2017).

## Results and Discussion

As depicted in section 3, the numerical error will depend on grid structure (that is, non-uniformity) as well as the advection flow configuration even though the discretization method is the same. This issue was taken into consideration by defining the initial scalar field at various locations. Figure 3 illustrates two examples of the initial distribution of the scalar field. In the first example, which corresponds to the control run, two scalar fields are located in the same position and the scalar variable is advected away from the boundary, as indicated by arrows. In the case of the nonuniform grid, this means that it is transported from fine vertical-resolution region to coarse verticalresolution region. In the lower-panel example, the left cell has the initial scalar field in the same position the upper panel. On the other hand, the right cell has the scalar field near the upper boundary so that it is transported from the coarse vertical-resolution region to the fine vertical-resolution region. It is expected that transport of the scalar field across different resolution grids will result in different error behavior. To assess the accuracy, the root mean squared error (RMSE) is measured for each circulation cell:

where *h _{s}* and

*h*mean the simulated field and reference field, respectively.

_{r}Figure 4 presents the simulation results for the control run. As is expected from flow configuration, the error distribution is exactly symmetric with respect to the vertical line *x* = *L*/2 (=0.5 in non-dimensional unit). This is because the centered FDM is employed in *x*-direction, which, unlike the vertical discretization, has the same error for both directions (negative and positive *x*-direction). Summary of simulations with various initial location of tracer, spatial resolutions, and non-uniformity, is shown in Table 1. The results clearly exhibit the second-order error convergence for the uniform grid. For instance, the error of case III is reduced by a factor of about 4 in comparison with the case II. The same holds true for the cases of III and IV. It is noticeable that the error for the nonuniform grid is quite sensitive to the initial location of the tracer as can be seen the results of the cases from V to X. In the cases of the nonuniform grid, the error of downdraft motion in the stretching-grid region (i.e., the grid-size increases with height) is larger than the downdraft in the shrinking-grid region (i.e., the gridsize decreases with height). This tendency becomes more conspicuous for increased nonuniformity, e.g., the cases of from V to X. Such a behavior is well illustrated in Fig. 5 which corresponds to the case V. Error patterns appear in wave train with alternating sign, and more contour lines, being indicative of larger error, is clearly seen in the left circulation cell. The reason for this asymmetric error behavior may be that the initial condition of Gaussian bell is resolved not in an identical accuracy. It is certain that the initial tracer located in the fine-grid region better resolves the Gaussian bell than that located in the coarse-grid region. This is well illustrated in Fig. 6 where the results of case V are plotted. Obviously, the initial tracer over the left cell is defined on smaller number of grid-points than the right cell.

To identify the combined effect of nonuniformity of grid and the stretching or shrinking motion in the vertical direction more clearly, an experiment with circular (i.e., axisymmetric) pattern of tracer distribution was conducted. Since both the flow pattern and the tracer distribution are axisymmetric, the error pattern also should be axisymmetric if the error source associated with the spatial discretization is axisym- metric. Result of this case is shown in Fig. 7, where the error pattern is not found to be axisymmetric. This is attributed to the non-uniformity of the grid and stretching or shrinking motion in the vertical direction. Combined effect of these two factors is represented with Eq. (11b). The error source may be termed as either numerical diffusion (*E*_{2}< 0) or numerical adverse-diffusion (*E*_{2}> 0). The largest error, therefore, is directly associated with the largest value of *E*_{2}, which may be observed in the updraft or downdraft near the boundaries. The adverse diffusion causes the time integration unstable while the diffusion does not. Of the two factors in Eq. (11b), the stretching or shrinking of the flow in the vertical direction is more important than the non-uniformity of the grid as for the two dimensional circulation. This is because the ratio of stretching velocity, for adjacent two grid points, is well above 3 or 4 for in the level near the boundaries, being much larger than the value of γ. Since the error pattern is advected along the circular trajectory, it can be imagined easily that the error at the initial stage is distributed near the boundaries. That is, the vertical distribution of |*E*_{2}| is severely confined to the region of strong velocity-stretching. If the results are applied to the atmospheric circulation model, it may be stated that the updraft in the boundary layer accompanied by the horizontal convergence is vulnerable to large numerical error.

## Concluding Remarks

The effect of nonuniform grid on the two dimensional transport phenomenon was investigated based on the numerical simulations and error analysis. Main concern of the study was the vertical discretization and nonuniform staggered-grid which resembled the weather forecast model in a qualitative sense. Discretization method in time was the thirdorder Runge-Kutta scheme, while for spatial discretization the FDM on the staggered-grid 3-points stencil in the vertical direction and the centered FDM in the horizontal direction were used. Error analysis indicated that the vertical discretization produces numerical adverse-diffusion or diffusion when the vertical-flow slope (that is, vertical convergence or divergence) is large. This effect tends to be strengthened by the presence of non-uniformity of the grid. However, the effect of the latter was smaller than that of the effect of vertical-flow slope for the typical non-uniformity of the grid. This result implies that in the numerical model the rising or sinking motion, accompanied by horizontal convergence or divergence, respectively, near the boundary, may be vulnerable to the numerical instability due to the adverse diffusion. And the possibility of the instability becomes large in the case of non-uniformity of grid-size.