-
We define the horizontal and vertical grid structures used in this study before introducing the solvers in detail. A regular latitude–longitude grid is used in the horizontal discretization, where winds and tracers are staggered with the Arakawa C-grid (Arakawa and Lamb, 1977). The wind component is defined at the middle of the cell interface, while the tracer concentration is defined at the cell center and the Earth’s poles. We use a stationary, non-uniform, mass-based, and terrain-following vertical coordinate (Kasahara, 1974; Simmons and Burridge, 1981; Laprise, 1992), with the positive direction pointing downward, and top-down numbering. A Lorenz-type staggering is adopted in the vertical, which means that the vertical coordinate velocity and vertical mass flux are defined at half levels, and the tracer mass-mixing ratio is defined at full levels.
With a general vertical coordinate
$ \eta $ , the tracer conservation equation without sources or sinks can be written aswhere
$ \pi $ is the hydrostatic pressure (Laprise, 1992),${\partial \pi }/{\partial \eta }$ is the pseudo-density, and$ q $ is the tracer mixing ratio.$ {\nabla }_{\eta }\cdot \left(\,\right) $ represents the horizontal divergence operator, which can be discretized based on the Gauss theorem.$ {\mathbf{V}}_{h} $ is the two-dimensional horizontal velocity vector, and$ \dot{\eta } $ is the vertical coordinate velocity. For the purpose of spatial discretization, the layer-averaged formulation [Eq. (12) in Zhang et al. (2020)] is used in this study.$ {\mathbf{V}}_{h} $ and$ q $ represent layer-averaged states; consequently, the tracer equation can be further expressed aswhere
$ \delta \pi $ denotes the layer mass in a full level (k),$ \delta {\pi }_{k}={\pi }_{k+1/2}-{\pi }_{k-1/2} $ . The$ \delta $ operator in the third term denotes the difference between two half levels, which is defined at the full levels:To implement the adaptively implicit vertical transport scheme, the mass-weighted vertical velocity
$ \left(\frac{\partial \pi }{\partial \eta }\dot{\eta }\right) $ is split into explicit and implicit velocities for the explicit and implicit transport parts, respectively. Thus, the equation can be semi-discretized aswhere the explicit and implicit vertical velocities are partitioned from the full vertical velocity:
where the splitting parameter
$\, \beta $ is a transition function responsible for the splitting of the vertical velocity between a fully explicit solver ($ \,\beta =1 $ ) and a fully implicit solver ($\, \beta =0 $ ), which is determined by the local Courant number. We do not describe the calculating formulation of$ \,\beta $ and related parameters since these can be found in Wicker and Skamarock (2020). This approach has the advantage of rendering vertical advection stable and maintaining high accuracy in locations with low Courant numbers.The fractional step, also known as operator-splitting, is attractive for solving the three-dimensional equation because of its efficiency and high accuracy in each spatial dimension (e.g., Lin and Rood, 1996, 1997; LeVeque, 2002; Putman and Lin, 2007). The operator-splitting is applied between horizontal and vertical operators to make it simpler to implement an implicit integration scheme in the vertical, although the operator-splitting is first-order accurate in time. A form of Strang splitting (Strang, 1968) can also be used to achieve second-order accuracy. After calculating the horizontal operator, the updated tracer mass-mixing ratio is used for the vertical operator.
The computational procedures are as follows:
(i) Calculate the horizontal flux for the intermediate variable using
where the superscripts n and * denote the state at time level n and intermediate time level, respectively. In our model, the Lin and Rood (1996) FFSL approach is used, which is based on forward Euler explicit time differencing (a two time-level scheme). The essence of the FFSL method is that the two-dimensional flux divergence is split into two orthogonal one-dimensional flux-form transport operators. The flux is approximately equal to the air-mass flux times the tracer mixing ratio on the cell interface in a one-dimensional direction. Either the piecewise linear van Leer method or the Piecewise Parabolic Method (PPM, Colella and Woodward, 1984; Carpenter et al., 1990) can be used to reconstruct the tracer in a mesh cell by using mean values of surrounding cells. We choose the PPM for our model because of its accuracy and computational efficiency. In addition, a monotonicity constraint (or slope-limiter) can be imposed on the one-dimensional discrete solution, but, due to dimensional splitting, this does not ensure exact monotonicity in the overall two-dimensional solution (Lin and Rood, 1996; Kent et al., 2012). The slope-limiter would damp the numerical oscillations but introduce diffusion in the numerical solutions. We used the modified version of the PPM limiter from Appendix B of Lin (2004) for both the horizontal and vertical advections. A minor correction regarding the limiter is noted in the appendix of this paper.
A semi-Lagrangian approach is also used in the longitudinal (east–west) direction on the latitude–longitude grid to alleviate the pole problem. The flux across cell interfaces is calculated by integrating the swept volume over a large number of upstream zonal cells. This approach enlarges the upwind-biased computation stencil that relaxes the CFL linear stability constraint, especially near the poles, as stated in section 3 of Lin and Rood (1996). Since the method is in the flux-form (the FF abbreviation in FFSL), it guarantees the total tracer mass conservation, whereas the conventional method that tracks computational grid points does not make this guarantee. The pure Eulerian counterpart is used in the meridional (north–south) direction, so the time-step size is only constrained by the meridional Courant number.
For comparison, we also implement a third-order upwind flux formulation combined with a third-order Runge–Kutta (RK3) time integrator (Wicker and Skamarock, 2002; hereafter RK3O3) in the horizontal solver. Unlike FFSL, RK3O3 is a method of lines (Hyman, 1979; Shchepetkin, 2015), where the spatial and temporal discretizations are decoupled and treated separately. The FCT algorithm (Zalesak, 1979) is applied in RK3O3 to restore positive definiteness or monotonicity. It should be noted that the FCT operator can be applied either on each sub-step or the last sub-step of the RK3 scheme. In our tests, the former guarantees strict monotonicity, while the latter cannot (see section 3.1).
(ii) Update the mass-weighted mixing ratio in the vertical explicit residual with
where
${q}^{*}={{\left(\delta \pi q\right)}^{*}}/{\delta {\pi }^{n+1}}$ , the superscript ** denotes the state at intermediate time level, and the layer mass$ \delta {\pi }^{n+1} $ can be obtained by the continuity equation in the dynamical core of a full model. The layer mass is unchanged with time in the passive tracer advection tests. The explicit vertical velocity$ {\left(\frac{\partial \pi }{\partial \eta }\dot{\eta }\right)}_{\mathrm{E}} $ is used in the one-dimensional finite-volume scheme. As in the horizontal solver, we can obtain the vertical flux at half levels by using the third-order upwind flux operator or the PPM finite-volume method. The third-order upwind flux operator combines with the RK3 time-integrator method (RK3O3), whereas the PPM flux operator is equipped with the forward Eulerian time integration. Similarly, the FCT and slope-limiter can be imposed on the RK3O3 and PPM flux operators, respectively. The performance of these schemes will be evaluated by the following test cases.(iii) At the last step, update the tracer using an implicit algorithm with
Because it is the only inherently monotonic linear scheme (Godunov and Bohachevsky, 1959), the first-order upwind scheme is chosen to calculate the vertical mass flux at the half interface. While other higher-order nonlinear schemes are theoretically possible, the first-order upwind method is typically more cost-effective. Then, the mixing ratio at a new time level can be obtained by
${q}^{n+1}= {{\left(\delta \pi q\right)}^{n+1}}/{\delta {\pi }^{n+1}}$ . Since each unknown$ {q}_{k}^{n+1} $ depends on its upper$ {q}_{k-1}^{n+1} $ and lower$ {q}_{k+1}^{n+1} $ neighbors (subscript$ k $ is used for the variable located in the full layer in the vertical direction), the above set of equations form a tridiagonal system in a vertical column. The system of linear equations can be solved with a tridiagonal matrix algorithm.The summary of flux operators and monotonic limiters used in the model for the idealized tests is provided in Table 1.
The horizontal flux operator/limiter The vertical explicit flux operator/limiter The vertical implicit flux operator FFSL/slope-limiter PPM/slope-limiter 1st-order upwind RK3O3/FCT RK3O3/FCT Table 1. List of the horizontal and vertical flux operators, along with the corresponding limiters in the tracer transport model.
The horizontal flux operator/limiter | The vertical explicit flux operator/limiter | The vertical implicit flux operator |
FFSL/slope-limiter | PPM/slope-limiter | 1st-order upwind |
RK3O3/FCT | RK3O3/FCT |