Implementation of a Conservative Two-step Shape-Preserving Advection Scheme on a Spherical Icosahedral Hexagonal Geodesic Grid

With the rapid development of supercomputing power, global atmospheric models will have ultra-high horizontal resolutions that can include various scales of atmospheric motion (e.g., Miura et al., 2007; Smolarkiewicz et al., 2015; Heinzeller et al., 2016). Quasi-uniform grids are generally preferred for this kind of model rather than regular latitude-longitude grids. The latter suffers from the well-known polar problem because of the convergence of zonal distance near the pole. This problem not only limits the numerical stability of discretized dynamical equations, but also introduces an unnecessary mesh refinement that may affect the appropriateness of subgrid-scale formulations.

For quasi-uniform mesh systems, the grid configuration and numerical algorithm for solving global models have undergone extensive developments. (Williamson, 1968) and (Sadourny et al., 1968) introduced a geodesic grid based on the projection of an inscribed subdivided icosahedron onto the sphere [(Sadourny, 1972) also introduced the cubed-sphere mesh]. Since then, many studies have investigated the solving of atmospheric equations on icosahedral-like mesh systems. These studies include solving shallow-water equations (e.g., Cullen, 1974; Masuda and Ohnishi, 1987; Heikes and Randall, 1995a, 1995b; Thuburn, 1997; Stuhne and Peltier, 1999; Tomita et al., 2001; Ringler and Randall, 2002; Bonaventura and Ringler, 2005; Walko and Avissar, 2008a; Rípodas et al., 2009; Lee and MacDonald, 2009; Weller et al., 2009; Ii and Xiao, 2010), formulating atmospheric dynamical cores (e.g., Ringler et al., 2000; Tomita and Satoh, 2004; Walko and Avissar, 2008b; Gassmann, 2013; Wan et al., 2013; Zäangl et al., 2015), and building realistic global atmospheric models (e.g., Randall et al., 2002; Majewski et al., 2002; Satoh et al., 2008; Skamarock et al., 2012; Bleck et al., 2015). Based on the choice of the control volume, icosahedral mesh systems can be categorized into triangular meshes and hexagonal-pentagonal meshes. The latter type was chosen in this study as an example to demonstrate our application.

In addition to the mesh configuration, the numerical operator is another component of solving atmospheric equations. In an atmospheric model, the transport/advection equation may be the simplest prognostic equation and requires approximating only a single operator. Hence, solving the advection equation is frequently used to test numerical methods (e.g., Lauritzen et al., 2014). A desired transport scheme needs to meet several key properties, including conservation of the tracer mass, the shape-preserving property that denies the occurrence of new maxima and minima (under non-divergent flows), preserving the correlation of tracers, and being computationally efficient in a spherical geometry. Finite-difference and finite-volume methods are usually employed for flux-form schemes (e.g., van Leer, 1977; Zalesak, 1979; Smolarkiewicz, 1984; Leonard, 1991; Lin et al., 1994; Yu, 1994; Hundsdorfer et al., 1995; Li, 2008). Many flux-form schemes have also been adopted or developed on unstructured grids. (Löhner et al., 1988) extended the flux-corrected transport (FCT; Boris and Book, 1973; Zalesak, 1979) scheme to an unstructured triangular grid with a finite-element discretization. (Smolarkiewicz and Szmelter, 2005) presented an unstructured-grid formulation of the multidimensional positive definite advection transport algorithm (Smolarkiewicz, 2006) in an arbitrary finite-volume framework. (Putman and Lin, 2007) generalized the flux-form semi-Lagrangian (Lin and Rood, 1996) scheme for quasi-uniform cubed-sphere grids. On icosahedral-mesh applications, (Thuburn, 1997) implemented a conservative and shape-preserving advection scheme based on polynomial interpolations (Thuburn, 1995). (Lipscomb and Ringler, 2005) reformulated an incremental remapping-style advection scheme on a spherical geodesic grid. (Miura, 2007) designed a simplified incremental remapping-style transport scheme for icosahedral grids. This scheme was further extended in (Skamarock and Menchaca, 2010) using higher-order reconstructions. (Steppeler et al., 2008) proposed a method to achieve high-order finite-difference schemes on icosahedral-type grids. (Skamarock and Gassmann, 2011) generalized a set of finite-difference operators for spherical hexagonal meshes. Other schemes that have been implemented on icosahedral grids include those by (Niwa et al., 2011), (Miura and Skamarock, 2013), (Chen et al., 2014), (Dubey et al., 2014) and others.

In this study, we implemented a conservative finite-difference advection scheme, called the Two-step Shape-Preserving Advection Scheme (TSPAS; Yu, 1994), onto the spherical icosahedral-hexagonal-pentagonal (IHP) grid. The original scheme has been applied in several numerical models (e.g., Yu, 1995; Yu and Xu, 2004; Wang et al., 2004; Xiao et al., 2008; Shi et al., 2009; Zhang et al., 2013; Yu et al., 2015; Zhang and Li, 2016), all of which are built on regular latitude-longitude grids. To apply TSPAS on a geodesic grid, the original scheme was further generalized in a control-volume framework. TSPAS is a compound scheme that combines a low-order monotone scheme [the upwind (UP) scheme] and a high-order scheme [the Lax-Wendroff (LW) scheme (Lax and Wendroff, 1960)] that typically cannot ensure monotonicity. This approach sounds similar to the FCT approach. The difference is that in the FCT approach, the flux across each edge is computed by blending low- and high-order schemes with an anti-diffusive procedure; whereas, in the TSPAS approach, the flux across each edge is computed by either low- or high-order schemes.

In this paper, we detail the computational procedures of the generalized TSPAS and utilize several numerical tests to check the accuracy, efficiency and robustness of the implemented scheme. Numerical errors that arise from the discretized scheme are also computed. The two-step flux-form integration conserves the tracer mass and preserves the shape in a simple and efficient way. The FCT approach is also generalized in Appendix A, and the two schemes are compared.

The remainder of the paper is organized as follows. Section 2 describes the implementation of TSPAS on the IHP grid, and section 3 presents the results from various numerical tests. Concluding remarks are given in section 4.

The construction of an icosahedral grid starts from a raw icosahedron (inscribed inside a sphere), which has 12 vertices, 30 edges and 20 triangular faces. One-level finer grids are then generated by bisecting the edges of the former coarser grids (e.g., Baumgardner and Frederickson, 1985; Heikes and Randall, 1995a; Satoh et al., 2014). The nth bisection of the icosahedron is typically called grid level (Glevel) n. The division generates more spherical triangular faces that are formed by vertices. The vertices are then projected onto the sphere.

The spherical hexagons/pentagons are constructed by connecting the circumcenters (Fig. 1) of six/five neighboring spherical triangles (that share a common vertex). Pentagons only occur at the 12 vertices of the raw icosahedron. A hexagon that is formed by the circumcenters of spherical triangles is also referred to as a Voronoi (e.g., Heikes and Randall, 1995a). This differs from the barycenter-type hexagon, which is formed by connecting the barycenters of neighboring spherical triangles (e.g., Tomita et al., 2001).

Because spherical hexagons are not regular hexagons and are sort of distorted, grid optimizations are typically used for IHP grids to reduce truncation errors (Miura and Kimoto, 2005; Heikes et al., 2013). Many methods have been proposed for optimizing IHP grids. (Heikes and Randall, 1995b) slightly adjusted the intersection points of triangular and hexagonal/pentagonal grids to improve the accuracy. (Tomita et al., 2001) designed a spring grid system to reduce the grid noise that originates from the use of the Arakawa-A (Arakawa and Lamb, 1977) grid. Based on methods proposed by (Du et al., 1999) and (Du et al., 2003), (Ringler et al., 2008) introduced the Spherical Centroid Voronoi Tessellation (SCVT) for use in earth system modeling. The SCVT method is a more generalized surface tessellation approach that does not necessarily require starting from a subdivided icosahedron. When this method is used to optimize an IHP grid, the resultant grid can achieve the collocation and orthogonality properties, while the bisection property is not constrained (Miura and Kimoto, 2005). The SCVT optimization is used in this study. Open source icosahedral mesh tools are adopted to construct the mesh infrastructure (Renka, 1997; Peixoto and Barros, 2013), including the mesh generation and optimization.

Figure 1.  Illustration of the C grid discretization with a spherical hexagonal grid. The red triangles are the vertices of the spherical triangles (formed by the dashed lines), and the solid circles denote the circumcenters of spherical triangles; the solid lines that connect the solid circles form a hexagon. The red triangles are also the centroids of the hexagons given the SCVT optimization. See section 2.2 for additional details.

For simplicity, only a two-dimensional case is considered in this study. The flux-form advection equation may be written in the following form: \begin{equation} \dfrac{\partial\rho q}{\partial t}=-\nabla\cdot(\rho q{V}) , (1) \end{equation} where ρ is a density-like quantity, q is a mixing-ratio-like quantity, and V is the velocity vector. The continuity equation can be achieved by simply setting q to 1, which gives \begin{equation} \dfrac{\partial\rho}{\partial t}=-\nabla\cdot(\rho {V}) . (2)\end{equation}

Figure 1 shows the Arakawa-C control volume discretization in the context of an IHP mesh. The scalar values (e.g., ρ and q) are placed at the centroid of a control volume cell, which is also the vertex of a spherical triangle under the SCVT optimization. The velocity fields are defined at the midpoint of each hexagonal cell wall. Let Q=ρ q; then, the right side of Eq. (1) can be transformed into a line integral using Gauss's theorem: \begin{equation} \dfrac{\partial Q_0}{\partial t}=-\dfrac{1}{S_0}\sum_{i=1}^{i=N}{F}_i\cdot {n}_il_i , (3)\end{equation} where Q₀ is the scalar value at the centroid of a hexagonal/pentagonal cell 0, S₀ is the area of cell 0. A variable with a subscript denotes the target variable at a particular location (similar hereafter). n_i is a unit outward normal vector of the cell wall i, N is 6 for a hexagonal cell and 5 for a pentagonal cell, l_i is the spherical length of the cell wall i, and is located between the vertex 0 and i. F_i has a form of \(\tilde{Q}_iV_i\), which denotes the mass flux across the cell wall i (\(\tilde{Q}_i\) is an imaginary quantity at the cell wall). F_i is assumed to be located at the intersecting point between the Voronoi's and triangle's edges (Fig. 1).

With these configurations, the next step is to build the divergence operator with an appropriate spatial scheme to calculate the flux at the cell wall and an appropriate time scheme for the time integration. The transport equation can then be numerically solved. Combining TSPAS with the configuration of the IHP grid, it may be formulated as follows.

The first step is to calculate an intermediate flux at each cell wall using the high-order scheme but with a slightly increased flux (* denotes an intermediate stage): \begin{equation} F_i^*={F}_i\cdot {n}_i=\beta_0\left[0.5U_i(Q_0+Q_i)-0.5|U_i|(Q_i-Q_0)\dfrac{|U_i|\Delta t}{\Delta m_{i0}}\right] ,(4) \end{equation} where β₀ is a scalar that is placed at vertex 0 and is used to add a slight increment to the actual flux. The formula for calculating β₀ will be given in section 2.3, and β₀≥ 1 is a prerequisite. U_i=V_i· n_i is the wind speed normal to the cell wall i, (the brown arrow in Fig. 1; the sign is outward positive). Q₀ and Q_i (Fig. 1) are scalar values defined at vertices 0 and i. ∆ t is the time step, and ∆ m_i0 is the spherical length between vertex 0 and i. Next, we need to pre-update the scalar field at vertex 0 using Eqs. (3) and (4) and evaluate the updated \(Q_0^*[Q_0^*-Q_0-(∆ t/S)\sum_{i=1}^{i=N}F_i^*l_i\). An A value is then introduced and defined at each vertex as follows: \begin{eqnarray} A_0&=&(Q_0^*-Q_{\max})(Q_0^*-Q_{\min}) ,(5)\\ Q_{\max}&=&{\rm MAX}(Q_0,Q_1,Q_2,\ldots,Q_N) ,(6)\\ Q_{\min}&=&{\rm MIN}(Q_0,Q_1,Q_2,\ldots,Q_N) , (7)\end{eqnarray} where MAX and MIN denote taking the maximum and minimum values of vertex 0 and all of its neighboring vertices. By applying this operation globally, the entire A field will be evaluated.

The second step is to choose an appropriate scheme to perform the real/final integration. If A values at both vertex 0 and i are less than zero, the flux at the cell wall i will be calculated using the high-order LW scheme as follows: \begin{eqnarray} F_{{\rm LW},i}&=&{F}_{{\rm LW},i}\cdot {n}_i\nonumber\\ &=&0.5U_i(Q_0+Q_i)-0.5|U_i|(Q_i-Q_0)\dfrac{|U_i|\Delta t}{\Delta m_{i0}} ;(8) \end{eqnarray} otherwise, the flux will be calculated using the UP scheme as follows: \begin{equation} F_{{UP},i}={F}_{{\rm UP},i}\cdot {n}_i=0.5U_i(Q_0+Q_i)-0.5|U_i|(Q_i-Q_0) . (9)\end{equation} The choice between the two schemes can be achieved by using an "if statement" or Eqs. (9)-(12) in (Yu, 1994). This is the entire integration procedure of TSPAS on the IHP grid. Note that by using Eqs. (8) or (9) for each cell wall, we are assuming a one-dimensional finite-difference operator in each direction. The flux in each direction of a given cell wall is calculated using the same equation [Eq. (8) or Eq. (9)] but with opposite signs. Thus, the global sum of the fluxes will vanish, and the conservation is guaranteed (it has been checked that the scheme conserves the tracer mass to machine precision). A concise flowchart of the computational procedure is presented in Appendix B.

Up to this point, only β is unknown. Starting from Eq. (3), the discretized two-step integration with a forward-in-time scheme from time level n to n+1 may be written as \begin{eqnarray} Q_0^*&=&Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_i^*l_i ,(10)\\ Q_{n+1,0}&=&Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_il_i . (11)\end{eqnarray} Meanwhile, considering Eqs. (4), (8) and (9), we may obtain \begin{eqnarray} F_{{\rm UP},i}&=&F_{{\rm LW},i}+0.5|U_i|(Q_i-Q_0)\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right) ,(12)\\ F_i^*&=&\beta_0F_{{\rm LW},i} . (13)\end{eqnarray} Note that Q₀ and Q_n,0 are the same unless otherwise explicitly mentioned. As stated previously, TSPAS combines the fluxes that are computed by the LW and UP schemes. Hence, only three possible situations may occur:

(1) All of the fluxes are calculated by the LW scheme. Thus, \begin{equation} Q_{n+1,0}=Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_{{\rm LW},i}l_i .(14) \end{equation} Using Eqs. (10) and (13), we obtain \begin{equation} -\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_{{\rm LW},i}l_i=\dfrac{Q_0^*-Q_{n,0}}{\beta_0} . (15)\end{equation} Thus, \begin{equation} Q_{n+1,0}=Q_{n,0}+\dfrac{Q_0^*-Q_{n,0}}{\beta_0}=\dfrac{Q_0^*+(\beta_0-1)Q_{n,0}}{\beta_0} . (16)\end{equation}

Note that as long as one cell wall around vertex 0 chooses the LW scheme (A₀<0), there must be \(Q_\min\leq Q_0^*\leq Q_\max\). Because β₀≥ 1, we obtain \begin{eqnarray} \dfrac{Q_{\min}+(\beta_0-1)Q_{\min}}{\beta_0}&\leq& \dfrac{Q_0^*+(\beta_0-1)Q_{n,0}}{\beta_0}\nonumber\\ &\leq& \dfrac{Q_{\max}+(\beta_0-1)Q_{\max}}{\beta_0} . (17)\end{eqnarray} Therefore, \(Q_\min\leq Q_n+1,0\leq Q_\max\). Thus, the shape-preserving property is guaranteed when β₀≥ 1.

(2) All of the fluxes are calculated by the UP scheme. The shape-preserving property is guaranteed when the stability condition ([(|U_i|∆ t)/∆ m_i0]-1≤ 0) is satisfied.

(3) A total of k fluxes are calculated by the UP scheme, and N-k fluxes are calculated by the LW scheme. Using Eqs. (11) and (12), we obtain \begin{eqnarray} Q_{n+1,0}&\!=\!&Q_{n,0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=N}F_{{\rm LW},i}l_i-\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|(Q_i-Q_0)\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] ,(18)\\ Q_{n+1,0}&\!=\!&\dfrac{Q_0^*\!+\!(\beta_0\!-\!1)Q_{n,0}}{\beta_0}-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\! \left[0.5|U_i|Q_i\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\!+\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i| Q_0\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] ,(19)\\ Q_{n+1,0}&\!=\!&\dfrac{Q_0^*}{\beta_0}+\!\left\{1\!-\!\dfrac{1}{\beta_0}\!+\!\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\right\}Q_0\!-\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|Q_i\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] . (20)\end{eqnarray} Because 1/β₀≥ 0 and [(|U_i|∆ t)/∆ m_i0]-1≤ 0 when the stability condition is satisfied, it only requires \begin{equation} 1-\dfrac{1}{\beta_0}+\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\geq 0 \ \ (21)\end{equation} to achieve \begin{eqnarray} &&\dfrac{Q_{\min}}{\beta_0}+\left\{1-\dfrac{1}{\beta_0}+\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)\right]l_i\right\}Q_{\min}-\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|Q_{\min}\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right]\leq Q_{n+1,0}\leq\nonumber\\ &&\dfrac{Q_{\max}}{\beta_0}+\left\{1-\dfrac{1}{\beta_0}+\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k} \left[0.5|U_i|\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)\right]l_i\right\}Q_{\max}-\quad\nonumber\\ &&\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|Q_{\max}\left(\dfrac{|U_i|\Delta t}{\Delta m_{i0}}-1\right)l_i\right] . (22)\end{eqnarray} That is, \(Q_\min\leq Q_n+1,0\leq Q_\max\).

Therefore, the shape-preserving property is guaranteed when Eq. (21) is satisfied. Thus, we require \begin{equation} \beta_0\geq \dfrac{1}{1-\frac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|\left(1-\frac{|U_i|\Delta t}{\Delta m_{i0}}\right)l_i\right]} \ \ (23)\end{equation} under the condition \begin{equation} 1-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}\left[0.5|U_i|\left(1-\dfrac{|U_i|\Delta t}{\Delta m_{i0}}\right)l_i\right]>0 . (24)\end{equation} Let γ_i=|U_i|[1-(|U_i|∆ t)/∆ m_i0]l_i and note that γ_i≥ 0 under the stability condition. We can find a value of γ_max that satisfies \begin{equation} \gamma_{\max}={\rm MAX}(\gamma_1,\gamma_2,\gamma_3\cdots \gamma_N) . (25)\end{equation} Then, \begin{equation} 1-\dfrac{\Delta t}{S_0}\sum_{i=1}^{i=k}0.5\gamma_i>1-\dfrac{k\Delta t}{S_0}0.5\gamma_{\max} ;(26) \end{equation} namely, \begin{equation} \dfrac{1}{1-\frac{k\Delta t}{S_0}0.5\gamma_{\max}}>\dfrac{1}{1-\frac{\Delta t}{S_0}\sum_{i=1}^{i=k}0.5\gamma_i} . (27)\end{equation} Thus, we let \begin{equation} \beta_0={\rm MAX}\left(1,\dfrac{2}{2-\frac{k\Delta t}{S_0}\gamma_{\max}}\right) . (28)\end{equation}

Eq. (28) is a generalized form of β that was used by (Yu, 1994), and it is a formula that deals with all possible situations. Ideally, the value of k may be determined for each vertex at each time step using iterative routines. However, this will substantially increase the computational burden. For simplicity, we use k=3 globally throughout the integration. Normally, a larger (smaller) k value leads to a larger (smaller) β value, and the scheme may choose the monotone scheme more (less).

The numerical tests that are conducted in this study are given in (Williamson et al., 1992) and (Nair and Lauritzen, 2010). We use the solid body rotation test of (Williamson et al., 1992) and the deformational flow cases (1-4) of (Nair and Lauritzen, 2010). The numerical errors are calculated following (Williamson et al., 1992), including \begin{eqnarray} L_1&=&\dfrac{I(|Q_{\rm C}-Q_{\rm T}|)}{I(|Q_{\rm T}|)} ,(29)\\[1mm] L_2&=&\sqrt{\dfrac{I[(Q_{\rm C}-Q_{\rm T})^2]}{I[(Q_{\rm T})^2]}} ,(30)\\[1mm] L_{\infty}&=&\dfrac{\max(|Q_{\rm C}-Q_{\rm T}|)}{\max(|Q_{\rm T}|)} ,(31)\\[1mm] h_{\max}&=&\dfrac{\max(Q_{\rm C})-\max(Q_{\rm T})}{\Delta h} , (32)\end{eqnarray} \begin{eqnarray} h_{\min}&=&\dfrac{\min(Q_{\rm C})-\min(Q_{\rm T})}{\Delta h} ,(33)\\[1mm] \Delta h&=&\max(Q_{\rm T})-\min(Q_{\rm T}) , (34)\end{eqnarray} where I denotes the global integration, Q _C denotes the computational solution, and Q _T denotes the true/analytical solution, \(h_\max,h_\min\) and ∆ h represent deviations. The operators max and \(\min\) denote taking the maximum and minimum values on the globe, respectively. All of the numerical tests are performed on a unit sphere (radius = 1) for the normalization. The duration of integration is T=5, as suggested by (Nair and Lauritzen, 2010). We solve Eqs. (1) and (2) with a constant initial density field (ρ=1). Results are presented and diagnosed for Q in this paper. The initial spatial distributions of the tracer mixing ratio q are given in each subsection below. The resolutions of the icosahedral mesh are Glevel 4, 5 and 6, which have 2562, 10242 and 40962 points, respectively. The resolutions correspond to approximately 4^°, 2^° and 1^° radians. The integration is divided into 600, 1200 and 2400 steps for Glevel 4, 5 and 6, such that the maximum Courant number has an approximately constant value at all resolutions. All of the results are post-interpolated to a uniform rectangular grid. We run each test with TSPAS and FCT. The error statistics from the two schemes are given in Tables 1-7. Figures that show the FCT's results are provided in a supplementary document.

In the solid body rotation test, a cosine bell is advected by prescribed ambient winds on the sphere. The wind field and initial scalar values are as follows: \begin{eqnarray} u&=&u_0(\cos\theta\cos\alpha+\sin\theta \cos\lambda\sin\alpha) ,(35)\\ v&=&-u_0\sin\lambda\sin\alpha ,(36)\\ u_0&=&\dfrac{2\pi}{T} ,(37)\\ q&=&\left\{ \begin{array}{l@{\quad}l} 0.5\left[1+\cos\left(\dfrac{\pi r}{R}\right)\right] & if\ r<R\\ 0 & if\ r\geq R \end{array} \right.,\quad (38)\end{eqnarray} where θ is the latitude, Λ is the longitude, α is the angle between the axis of the solid body rotation and the polar axis of the spherical coordinate system, and r is the great circle distance between vertex i and the center of the initial cosine bell. This can be calculated as follows: \begin{equation} r=\arccos(\sin\theta_i\sin\theta_0+\cos\theta_i\cos\theta_0\cos(\lambda_0-\lambda_i)) ,(39)\end{equation} where i and 0 denote the positions of the vertex and the center of the initial cosine bell, respectively. R=1/3 is the radius of the cosine bell. After the duration T, the cosine bell returns to its initial position. Tests are performed for α=0 and π/2.

Figure 2 shows the results from the test with α=0, in which the cosine bell is located at the equator (in the sense of earth). The results are shown at t=0, T/2, T for three resolutions. The error measures for these tests are given in Table 1. The numerical precision increases with increasing resolution. No overshooting or undershooting values occur in this case for TSPAS. A comparison of the results of TSPAS and FCT shows that TSPAS generally has slightly smaller L₁,L₂, and L_\infty errors at the different resolutions. FCT maintains larger maximum values, while TSPAS more effectively prevents the undershooting values.

Figure 2.  Solid body rotation tests with $\alpha=0$ for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

Figure 3 shows the results from a test with α=π/2, in which the cosine bell starts moving from the South Pole (in the sense of earth). The numerical precision (Table 2) is quite similar to those of the test with α=0 (Table 1), owing to the quasi-uniform mesh in the polar and equator regions. Note that although most of the results in Tables 1 and 2 appear to be the same, they are actually different if more decimal places are considered. The shape-preserving feature of TSPAS is also guaranteed. The differences between TSPAS and FCT are similar to those in Table 1.

Figure 3.  Solid body rotation tests with $\alpha=\pi/2$ for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

Figure 4.  Deformational flow test Case 1 for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T/2$ (bottom).

(Nair and Lauritzen, 2010) proposed a class of benchmark deformational flow test cases for the two-dimensional horizontal linear transport problems on a sphere. The scalar field undergoes severe deformation during the simulation, and the flow reverses its course at the half time and the scalar field returns to its initial position and shape at the final time step. This process makes the analytical solution available at the end of the simulation, and facilitates an assessment of the accuracy of the underlying transport scheme. We conducted four cases, as described in (Nair and Lauritzen, 2010). The selected initial scalar fields are two symmetrically located cosine bells that are defined as follows: \begin{equation} h_i(\lambda,\theta)=\dfrac{1}{2}\left[1+\cos\left(\dfrac{\pi r_i}{R}\right)\right]if\ r_i<R , (40)\end{equation} where R=0.5 is the base radius of the bells, and r_i=r_i(Λ,θ) is the great circle distance between (Λ,θ) and a specified center (Λ_i,θ_i) of the initial cosine bell. The initial conditions consist of a background value b and two cosine bells with centers at position 1 and 2: \begin{equation} q(\lambda,\theta)=\left\{ \begin{array}{l@{\quad}l} b+ch_1(\lambda,\theta) & if\ r_1<R\\ b+ch_2(\lambda,\theta) & if\ r_2<R\\ b & {\rm otherwise} \end{array} \right., (41)\end{equation} where b=0.1 and c=0.9. The centers of the cosine bells depend on the choice of different test cases as follows.

3.2.1. Case 1

This is a non-divergent flow test case. The wind field is given as follows: \begin{eqnarray} u(\lambda,\theta,t)&=&k\sin^2\left(\dfrac{\lambda}{2}\right)\sin(2\theta)\cos\left(\dfrac{\pi t}{T}\right) ,(42)\\ v(\lambda,\theta,t)&=&\dfrac{k}{2}\sin\lambda\cos\theta\cos\left(\dfrac{\pi t}{T}\right) . (43)\end{eqnarray}

Figure 5.  Deformational flow test Case 2 for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

The flow parameter k is 2.4, and the centers of the initial cosine bells are kept at (Λ₁,θ₁)=(π,π/3) and (Λ₂,θ₂)=(π,-π/3). Figure 4 shows the results from this test case, and Table 3 shows the results from the error statistics. The numerical errors generally decrease as resolution increases. TSPAS shows slight undershooting values, which are generally smaller (absolute values, similar hereafter) than those produced by FCT. The L₁, L₂, and L_\infty errors from FCT are generally smaller than those from TSPAS at different resolutions, and FCT maintains larger maximum values.

Figure 6.  Deformational flow test Case 3 for Glevel 4 (left column), 5 (center column), and 6 (right column) at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

3.2.2. Case 2

This is also a non-divergent case. The wind field is given as follows: \begin{eqnarray} u(\lambda,\theta,t)=k\sin^2\lambda\sin(2\theta)\cos\left(\dfrac{\pi t}{T}\right) ,(44)\\ v(\lambda,\theta,t)=k\sin(2\lambda)\cos\theta\cos\left(\dfrac{\pi t}{T}\right) . (45)\end{eqnarray} The flow parameter k is 2. The centers of the initial cosine bells are kept at (5π/6,0) and (Λ₂,θ₂)=(7π/6,0). Figure 5 shows the results for this test case, and Table 4 shows the error statistics. For both schemes, the numerical errors for this case are overall larger than those of Case 1, but the undershooting values are generally smaller than those of Case 1. The differences between TSPAS and FCT are similar to those in Case 1. TSPAS generally produces smaller undershooting values than FCT, while FCT has relatively smaller errors and maintains larger maximum values.

3.2.3. Case 3

This is a divergent case. The wind field is given as follows: \begin{eqnarray} u(\lambda,\theta,t)&=&-k\sin^2(\lambda/2)\sin(2\theta)\cos^2\theta\cos\left(\dfrac{\pi t}{T}\right),(46)\\ v(\lambda,\theta,t)&=&\dfrac{k}{2}\sin\lambda\cos^3\theta\cos\left(\dfrac{\pi t}{T}\right) . (47)\end{eqnarray} The flow parameter k is 1. The centers of the initial cosine bells are kept at (Λ₁,θ₁)=(3π/4,0) and (Λ₂,θ₂)=(5π/4,0). Figure 6 shows the results from this test case, and Table 5 shows the error statistics. This case produces the smallest numerical errors among all four cases with cosine bells, but generates the largest undershooting values. The differences between the two schemes are similar to those in Case 1 and Case 2.

Figure 7.  Deformational flow test Case 4 for Glevel 4 (left column), 5 (center column), and 6 (right column), at time $t=0$ (top), $t=T/2$ (middle row) and $t=T$ (bottom).

3.2.4. Case 4

\begin{eqnarray} u(\lambda,\theta,t)&=&k\sin^2\lambda'\sin 2\theta\cos\left(\dfrac{\pi t}{T}\right)+2\pi\cos\theta/T ,(48)\\[-0.5mm] v(\lambda,\theta,t)&=&k\sin(2\lambda')\cos\theta\cos\left(\dfrac{\pi t}{T}\right) ,(49) \end{eqnarray} \begin{eqnarray} \lambda'&=&\lambda-\dfrac{2\pi t}{T}. (50)\end{eqnarray}

Figure 8. 

This is a non-divergent case. The flow parameter k is 2. The centers of the initial cosine bells are kept at (Λ₁,θ₁)=(5π/6,0) and (Λ₂,θ₂)=(7π/6,0). Figure 7 shows the results from this test case, and Table 6 shows the error statistics. This case generates the largest errors of the four cases. The undershooting values for this case are quite small. TSPAS prevents undershooting values more effectively, while FCT still maintains larger maximum values. In contrast to the previous three deformational tests, TSPAS produces smaller L₁, L₂, and L_\infty errors in this case.

3.2.5. Case 3 with slotted cylinders

To test the shape-preserving properties of the schemes in a more challenging way, we rerun Case 3 but with discontinuous slotted cylinders (Fig. 8). The errors statistics are presented in Table 7. Both schemes produce overshooting and undershooting values in this case. Similar to Case 3 with cosine bells, TSPAS generally produces smaller overshooting and undershooting values at various resolutions but with slightly larger errors.

3.2.6. Rates of convergence

To better check the order of accuracy, rates of convergence were estimated for each case, as shown in Table 8. The rate of convergence is computed by regressing the logarithms of error norms (L₁ and L₂) against the logarithms of the maximum angular distance between vertices at different resolutions. The regression coefficient is taken as the convergence rate. As shown in the table, the computed rates of convergence are around first-order accuracy, which is generally consistent with the fact that two schemes are composed by first-order and second-order schemes in practice.

In the above text, the generalized TSPAS is described. Several numerical tests were conducted to examine the simulation results, and a comparison between TSPAS and FCT was made. The table constituting Appendix C summarizes the overall performance of TSPAS and FCT. One may notice that these two schemes share a similar philosophy to satisfy the shape-preserving property; namely, find an optimal combination of schemes to maintain accuracy under the shape-preserving condition. The similar error norms and order of accuracy are consistent with the fact that both methods are composed of two same schemes. Nevertheless, they generate combinations in quite different ways.

The FCT approach blends fluxes computed by two schemes at each edge. The blending coefficient is computed by explicitly prohibiting those conditions when maximum net anti-diffusive fluxes are higher than total inflow/outflow anti-diffusive fluxes. This sets the range of the blending coefficient from zero to one. In the TSPAS approach, the value of β is used to modify the LW-updated scalar field. Only when the pre-updated scalar values at both sides of a cell edge satisfy the shape-preserving rule will the scheme choose the LW scheme; otherwise, it will choose the monotone scheme. Therefore, the combination generated by TSPAS may contain more monotone and diffused solutions than FCT. This explains why TSPAS generally produces smaller undershooting/overshooting values, while FCT maintains larger maximum values of the scalar fields. Meanwhile, the procedures needed by TSPAS to generate a combination are more simplified, which helps TSPAS to be less time-consuming (FCT costs ∼1.3-1.4 times more than TSPAS).

One may also notice that both schemes sometimes generate undershooting and overshooting. Reasons for this problem are complex, including the challenging flows and distributions of scalar values. From the aspect of the algorithm, the multidimensional application is a key reason. In the case of FCT, (Zalesak, 1979) discussed the impact of multidimensional applications on the monotonicity of the algorithm. We report its impact on TSPAS in the following text.

As formulated in section 2.3, the value of k is a consequence of multi-dimensional applications. This parameter will not be introduced in a one-dimensional case. k is used to quantify how many edges will be computed by UP. This is a precondition because the scheme will always choose a combination of UP and LW. The value of k does not set a limitation on the number of fluxes that are computed by the UP scheme. Once k is set, it determines how likely a cell wall will choose UP (or not choose). But this possibility is not entirely associated with k, because β is the parameter that balances the options, and β is only partially contributed by k. The problem of using a fixed k value is that the options between two schemes, sometimes, may not rigidly satisfy the shape-preserving rule. Hence, the scheme may generate under/overshooting values, especially under the environment of severe deformational flows and sharp gradients. Tests on different k values suggest that the sensitivity is not significant. We use k=3 because it is the median of the number of edges of a hexagonal mesh.

In this study, a finite-difference advection scheme, TSPAS, was generalized and tested on an icosahedral hexagonal-pentagonal geodesic grid. The generalized scheme conserves the tracer mass and preserves monotonicity in a simple and efficient approach. Several transport tests were used to check the accuracy, efficiency and robustness of the generalized scheme when it is applied on an optimized IHP grid. The results demonstrate the feasibility of the generalized two-step scheme on an unstructured geodesic grid. Tests at various resolutions show that the numerical precision increases as the resolution increases.

A comparison between TSPAS and FCT suggests that, generally, TSPAS produces smaller undershooting values, whereas FCT maintains larger maximum values. The errors derived from the L₁, L₂ and L\(\infty\) metrics show that the two schemes generate comparable errors. TSPAS generates smaller errors in two solid-body rotation tests and in the deformational test Case 4, whereas FCT produces smaller errors in the other four cases. Both schemes show an approximate first-order accuracy in idealized tests. Meanwhile, TSPAS is relatively more efficient in terms of the computational expense due to its simpler formulations, and it can be conveniently implemented on other unstructured grids. Moreover, for a practical global model application, the scheme can be easily adapted for a massively parallel environment given that it only requires a very limited node communication.