## Abstract

A standard nominally third-order upwind-biased spatial discretization of the flux-divergence operator was extended to a spherical icosahedral grid. The method can be used with multistage time-stepping schemes such as the Runge–Kutta method to compute the transport of variables on both hexagonal–pentagonal and triangular meshes. Two algorithms can be used to determine mesh cell face values: 1) interpolation using a quadratic function reconstructed subject to an integral constraint, or 2) calculation of the weighted mean of two linearly interpolated and extrapolated values. The first approach was adopted for a triangular mesh because the second approach depends on the mesh having a hexagonal or pentagonal shape. Both approaches were tested on the hexagonal–pentagonal mesh.

These schemes were subjected to standard transport tests on a spherical icosahedral grid. A three-stage Runge–Kutta time stepping method was used, and if necessary a flux limiter was applied to maintain monotonicity. The two different methods produced very similar solutions on a hexagonal–pentagonal mesh. Their accuracy was very close to the accuracy of a preexisting method designed for a Voronoi mesh only. When compared to another method that uses a quadratic polynomial interpolation, the phase error of the solutions was reduced, and their accuracy was much improved. The accuracies of the solutions were comparable on triangular and hexagonal–pentagonal meshes.

## 1. Introduction

A transport equation for atmospheric flow can be written in flux form as

where is the density, is the mixing ratio of a tracer, and is the flow velocity. The right-hand side of Eq. (1) is often called the flux-divergence term. The flux-divergence term appears in the Eulerian form of most of the prognostic equations for atmospheric flow if the equations are written in flux form. It is sometimes empirically evident that local temporal changes of temperature, moisture, and other atmospheric variables are greatly influenced by air masses transported from upstream. Accurate computation of the flux-divergence term is thus desirable for enhancing the accuracy and credibility of numerical simulations of atmospheric phenomena. The purpose of this report is to propose a new method to discretize the flux-divergence term on a spherical icosahedral grid.

A spherical icosahedral grid is generated by dividing an icosahedron and projecting it onto a sphere (e.g., Heikes and Randall 1995). These grids are beginning to be used in some atmospheric and oceanic models (e.g., Ringler et al. 2000; Tomita and Satoh 2004; Walko and Avissar 2008b; Skamarock et al. 2012; Gassmann 2013; Wan et al. 2013). One advantage of an icosahedral grid is its characteristic quasi uniformity, a feature that is suitable for maintaining the consistency of scale-dependent parameters in a model. Another advantage is that the structure of a spherical icosahedral grid fits naturally into massively parallel computer architectures. One example is the Nonhydrostatic Icosahedral Atmosphere Model (NICAM; Satoh et al. 2008), which was developed for the Earth Simulator to enable global “cloud resolving” simulations (e.g., Tomita et al. 2005; Miura et al. 2007). The NICAM is now being tested with more than 10 000 processors on the K computer, a Japanese new supercomputer in Kobe, Japan. As computer speeds increase, spherical icosahedral grids will probably be used more widely to study atmospheric phenomena and the earth's climate.

There are a number of options for spatially discretizing the flux-divergence term on a spherical icosahedral grid if single-stage forward time stepping is assumed. Lipscomb and Ringler (2005), Yeh (2007), and Miura (2007) used a linear reconstruction of tracer profiles to develop upwind-biased finite-volume schemes. Skamarock and Menchaca (2010) upgraded the scheme of Miura (2007) by replacing a linear reconstruction with a quadratic or quartic reconstruction. Miura and Skamarock (2013) proposed another quadratic reconstruction method to enable more accurate simulations with a small increase in computational costs. For triangular meshes, Walko and Avissar (2008a) used a second-order Crowley method, and Ii and Xiao (2010) developed third- and fourth-order multimoment transport schemes. In contrast, our options for transport schemes for multistage time-stepping schemes were relatively limited. A standard second-order centered scheme has been used with third-order Adams–Bashforth time stepping by Heikes and Randall (1995) or with the three-stage Runge–Kutta time stepping (RK3) of Wicker and Skamarock (2002) by Tomita and Satoh (2004). A quadratic reconstruction method was used with third-order Adams–Bashforth time stepping by Lee and MacDonald (2009, hereafter LM09). Weller et al. (2009) used another quadratic or cubic method with Crank–Nicholson time stepping.

Skamarock and Gassmann (2011, hereafter SG11) have recently extended an upwind-biased transport scheme that is accurate to the third-order under uniform flow to a spherical icosahedral grid. SG11 used RK3 for the time integration. The same type of scheme has frequently been used in atmospheric models (e.g., Skamarock and Klemp 2008; Ishida et al. 2010), but was not available on a spherical icosahedral grid until SG11.They showed that a midlatitude cyclonic system was better simulated with the SG11 scheme than with the standard second-order centered scheme. Here we follow the terminology of SG11 with respect to the order of accuracy of the scheme. In other words, we refer to this scheme that is third-order accurate under uniform flow simply as “the (standard) third-order scheme.”

NICAM, for example, uses the same RK3 for temporal discretization and a standard second-order centered scheme for the spatial discretization of the flux-divergence term in all cases except passive tracers. It might therefore be expected that replacement of the second-order centered scheme with the SG11 scheme would improve simulations by NICAM. Unfortunately, however, the SG11 scheme is not applicable to NICAM because NICAM uses a non-Voronoi mesh (Fig. 1b), whereas the SG11 scheme assumes a Voronoi mesh (Fig. 1a). NICAM as well as other models that use non-Voronoi meshes cannot use the SG11 scheme.

The Voronoi hexagonal–pentagonal mesh can be configured from the triangular mesh of a spherical icosahedral grid by connecting neighboring orthocenters of the triangular cells (e.g., Heikes and Randall 1995). Alternatively, another type of hexagonal–pentagonal mesh may be configured by connecting the neighboring barycenters of triangular cells (Tomita et al. 2001). This latter mesh is the hexagonal–pentagonal mesh used by NICAM. Miura and Kimoto (2005) characterized this kind of mesh as barycentric. One of the advantages of a barycentric mesh over the Voronoi mesh is that the distortions of hexagonal cells are smaller (Miura and Kimoto 2005). In contrast, the major disadvantage is that the line connecting two neighboring nodes is not orthogonal to the intersecting cell face. As will be reviewed in the next section, the SG11 scheme requires a parameter, , a length scale that represents the distance between two cell nodes. For the Voronoi mesh, can be specified directly as the distance between two cell nodes. For a barycentric mesh, however, the determination of is not straightforward because of the lack of orthogonality. Using the distance between two cell nodes for may lead to an overestimation or underestimation. To overcome this limitation, we have introduced a new algorithm that does not require orthogonality at the cell face and can be applied to hexagonal–pentagonal meshes in general, regardless of whether the configuration consists of a Voronoi or non-Voronoi mesh. In addition, combined with an interpolation procedure, this algorithm can also be applied to a triangular mesh.

Section 2 uses initial consideration of a standard, third-order scheme on a uniform, one-dimensional grid to introduce two different perspectives for the determination of the cell face tracer value. We introduce two new transport schemes on a hexagonal–pentagonal mesh that are analogous to transport schemes in the one-dimensional case. Application to a triangular mesh then follows. In section 3, the new schemes are subjected to a cosine-bell advection test (Williamson et al. 1992), a slotted-cylinder advection test (Lipscomb and Ringler 2005), and a deformational flow test (Nair and Lauritzen 2010). Section 4 provides a summary.

## 2. Transport schemes

### a. A three-stage Runge–Kutta time stepping

We used the same RK3 scheme as SG11 for the temporal integration of the transport equation in this work. We assumed the density to be constant and rewrote Eq. (1) as follows:

This simplification leaves the problem essentially unchanged because we can determine in the right-hand side of Eq. (1) by using the same algorithm for or using another algorithm such as the second-order centered scheme. Alternatively, we can recover Eq. (1) by substituting into in the right-hand side of Eq. (2) and for in the left-hand side of Eq. (2), as suggested by SG11. Using the RK3 scheme for the temporal discretization, we can update from time to by the following three-stage calculations:

if the flux-divergence term inside the brackets is discretized.

### b. One-dimensional case

We start by considering a one-dimensional transport problem along the *x* coordinate with a uniform velocity *c* > 0 (Fig. 2) as the basis of our new algorithms on a spherical icosahedral grid. Given this scenario, Eq. (2) is rewritten as

The advection term on the right-hand side can be approximated at the point *x*_{0} by standard third-order upwind-biased finite differencing (e.g., Durran 2010) as

where Δ*x* is the grid interval. This equation can be rewritten in flux form as

where the cell interface values and are determined as follows:

Equation (3) means that the flux-divergence term can be approximated by dividing the difference between outflow from the 0th cell through the cell interface at *x*_{1/2} and inflow to the 0th cell through the interface at *x*_{−1/2} by the grid interval.

where *β* = 1 gives the same form as Eq. (4), and *β* = 0 leads to the interface value that constitutes a part of fourth-order centered differencing (Wicker and Skamarock 2002). SG11 used the fact that the second-order spatial derivative of *q* can be discretized as

to rewrite Eq. (5) in the following form:

An alternative way to obtain the same result as Eq. (4) is to consider an interpolation that uses a piecewise parabolic reconstruction. The tracer profile inside the cell *x*_{0} is approximated by a quadratic profile (the gray line in Fig. 2b):

Using the constraints *f*(−Δ*x*) = *q*_{−1}, *f*(0) = *q*_{0} and *f*(Δ*x*) = *q*_{1}, the coefficients *a*_{0}, *a*_{1}, and *a*_{2} are determined, and Eq. (7) becomes

Substituting *x* = Δ*x*/2 into Eq. (8) enables the cell interface value to be determined as

This value does not equal Eq. (4), and a transport scheme based on this algorithm is thus different from the standard third-order scheme.

We now add a correction term to the quadratic profile in Eq. (8) as follows:

This adjustment produces the shift shown schematically in Fig. 2b. The form of is determined from the following constraint:

This constraint causes the shifted tracer profile integrated over the *x*_{0} cell to be equal to the tracer amount inside the cell (if *q*_{0} is the mean in the cell). Substituting Eq. (10) into Eq. (11) causes the form of to become

and the shifted quadratic profile is specified as follows:

Substituting *x* = Δ*x*/2 into Eq. (12) allows the cell interface value to be determined as follows:

It was shown that Eq. (4) is recovered by use of this upwind-biased interpolation to the cell interface with the shifted piecewise parabolic profile in Eq. (12) constrained by condition in Eq. (11). The transport scheme using this algorithm to determine and in Eq. (3) is therefore equivalent to the standard third-order scheme. We refer to this method as “type I” in this work and have applied it to both the hexagonal–pentagonal and triangular meshes.

At the end of this subsection, we introduce another algorithm, called “type II,” that recovers Eq. (4) through a combination of linear interpolation and extrapolation. In Fig. 2c, is linearly interpolated from adjacent tracer values on both sides of the interface as

and is linearly extrapolated from two tracer values located on the upwind side of the interface as

It is straightforward to show that

and thus the weighted mean of and recovers Eq. (4). This type-II algorithm is simpler and its implementation is easier than the type-I scheme. We will apply this algorithm to the hexagonal–pentagonal mesh, but its application to the triangular mesh is difficult as will be discussed in section 2d.

### c. Application to a hexagonal–pentagonal mesh

We used Gauss's theorem to discretize the flux-divergence term of Eq. (2) about the 0th cell on either the hexagonal–pentagonal mesh (Fig. 3a) or the triangular mesh (Fig. 3b) as follows:

where *q _{i}*(

*i*= 0, 1, … ,

*N*

_{0}) is the mixing ratio at node point

**H**

_{i}(Fig. 3a) or

**T**

_{i}(Fig. 3b),

*A*

_{0}is the area of the 0th cell,

*N*

_{0}is the number of the cells surrounding the 0th cell (

*N*

_{0}= 6 or 5 for the hexagonal–pentagonal mesh and

*N*

_{0}= 3 for the triangular mesh),

*l*(

_{i}*i*= 1 , … ,

*N*

_{0}) is the length of the

*i*th cell face shared by the 0th cell and the

*i*th surrounding cell,

**n**

_{i}is the unit vector normal to the

*i*th cell face and pointing outward from the 0th cell, and and are the mixing ratio and velocity vector, respectively, representing those quantities for the

*i*th cell face. If density is not constant, the right-hand side of Eq. (1) might be discretized as follows:

where is the density of the *i*th cell face and might be determined by using the same algorithm for or by using another algorithm such as the second-order centered scheme of Tomita et al. (2001). In the following, we explain how to determine if the flow direction is outward from the 0th cell across the *i*th cell face ).

SG11 developed a method to estimate by using an equation similar to Eq. (6) on a Voronoi-type hexagonal–pentagonal mesh. For example, the mixing ratio of the *i*th cell face is estimated as

where Δ*x _{i}* is the distance between

**H**

_{0}and

**H**

_{i}(Fig. 1a), and the second-order derivative term is estimated for the upwind-side cell (the 0th cell at this time) by a polynomial fit along the line connecting

**H**

_{0}and

**H**

_{i }. SG11 showed that this scheme worked very well on a Voronoi mesh. For non-Voronoi meshes, however, an ambiguity in specifying Δ

*x*makes it difficult to use the SG11 scheme directly. In addition, the first term on the right-hand side of Eq. (15) benefits from a property of the Voronoi mesh, namely that the cell face bisects exactly. Without this property, the second-order accuracy of this term may be compromised.

_{i}LM09 used another method to estimate , a piecewise parabolic reconstruction on a hexagonal–pentagonal mesh, but neither LM09 nor Lee et al. (2010) explained the details of the algorithm. The algorithm may be analogous to the one used to derive Eq. (9) in one-dimension and can be implemented as follows. We assume a quadratic profile in the neighborhood of **H**_{0}, that is,

We use least squares on the local two-dimensional *x*–*y* plane about **H**_{0} to determine the unknown coefficients *a _{j}* (

*j*= 1, … , 5) (Stuhne and Peltier 1996; Majewski et al. 2002). The coordinates of

**H**

_{i}on the local plane may be specified by , where

**n**

_{x}and

**n**

_{y}constitute an orthonormal basis set. By assuming

*q*=

_{i}*f*(

*x*,

_{i}*y*) (

_{i}*i*= 1, … ,

*N*

_{0}) for the surrounding cells, we obtain the following matrix equation:

where

The coefficient vector **a** is obtained from the equation where ^{#} means the inverse matrix (^{−1}) for the pentagonal cell and the pseudoinverse matrix (*) for the hexagonal cell. We used the singular vector decomposition (SVD) algorithm (e.g., Golub and Reinsch 1970) to compute ^{−1} and *. By substituting the coordinates of the center of the cell face, we estimated the mixing ratio of the *i*th cell face as follows:

where are the coordinates of **F**_{i} = (**T**_{i−1} + **T**_{i})/2 on the local coordinate system. Note that **T**_{i−1} means for *i* = 1. The LM09 scheme produces better results than the standard second-order-centered scheme (LM09; Lee et al. 2010). As discussed previously, however, the upwind-biased interpolation using the piecewise parabolic profile in Eq. (16) does not recover the standard third-order scheme. As a result, the LM09 scheme is less accurate than the SG11 scheme.

The LM09 scheme can be improved by applying the type-I algorithm. A shift term is added to Eq. (16) as follows:

and the integral constraint

is used to obtain

Following Miura and Skamarock (2013), we used the fact that the integral of any quadratic profile over the area of a triangle can be computed by using the three Gauss quadrature points at the center of each edge. Equation (21) was then discretized as

where *α _{i}* is the area of the subtriangle Δ

*H*

_{0}

*T*

_{i}_{−1}

*T*, and and are the coordinates of

_{i}**F**

_{i}= (

**T**

_{i−1}+

**T**

_{i})/2 and

**E**

_{i}= (

**H**

_{0}+

**T**

_{i})/2, respectively, in the local coordinate system (Fig. 3a).

Substitution of into Eq. (19) allows one to estimate the mixing ratio of the *i*th cell face as in a manner similar to Eq. (18). This method, however, is not accurate enough because *f _{c}* is quadratic, and thus, , with the exception of special cases. Alternatively, we could use two Gauss quadrature points and to compute the line integral of

*f*. Substituting the coordinates of into Eq. (19), the mixing ratio of the

_{c}*i*th cell face is determined as follows:

Use of Eq. (22) in Eq. (14) results in an extension of the standard third-order scheme on a hexagonal–pentagonal mesh. Notably, this type-I transport scheme does not depend on specific features of hexagons or pentagons, and it is thus usable on both Voronoi and non-Voronoi meshes.

We next apply the type-II algorithm to a hexagonal–pentagonal mesh. The type-II scheme is somewhat simpler than its type-I counterpart. By comparing results from the type-I and type-II schemes with each other, we can verify that those two algorithms are valid extensions of the standard third-order scheme. The linear interpolation and extrapolation to the *i*th cell face can be written in a two-dimensional *x–y* coordinate system about **H**_{0} as

With these two values, we construct the cell face value as

The gradient terms, and , are defined in terms of the interpolation to **F**_{1}, for example, as

for a hexagonal cell (Fig. 3c) and as

for a pentagonal cell (Fig. 3d). The gradient of the *i*th triangular cell, Δ*H*_{0}*H _{i}H_{i}*

_{+1}, is given by

where

Here, **n**_{x} and **n**_{y} are the orthonormal basis vectors of the local *x–y* coordinate system about **T**_{i}, and (*x*_{0}, *y*_{0}), (*x _{i}*,

*y*), and (

_{i}*x*

_{i+}_{1},

*y*

_{i+}_{1}) are the coordinates of

**H**

_{0},

**H**

_{i}, and

**H**

_{i+1}, respectively.

### d. Application to a triangular mesh

As demonstrated in the next section, the type-II scheme is as accurate as the type-I scheme, despite its greater ease of implementation. However, the type-I algorithm is superior to the type-II scheme in two respects. First, the type-II scheme depends on the axisymmetric structure of hexagons. As a result, its accuracy may be compromised if the distortion of a hexagonal cell is large and not an exact extension of the standard third-order scheme about the pentagons. In contrast, the type-I scheme makes no assumptions about the shape of hexagons and can treat pentagons in a consistent manner. Second, the application of the type-I scheme to a triangular mesh is straightforward, whereas application of the type-II scheme is not. It is difficult to determine the form of gradients that mimic Eq. (13) on a triangular mesh. Therefore, we have applied only the type-I algorithm to the triangular mesh.

Before constructing a quadratic profile, we interpolated tracer values from triangular cell nodes to hexagonal cell nodes by following the method of Holmes and Connel (1989) as improved by Kim et al. (2003) to reduce the degrees of freedom from the number of the triangular nodes to the number of the hexagonal or pentagonal nodes. For example, we interpolated the tracer values from to in Fig. 3b as follows. We wrote the interpolation in the form of a weighted mean as

where the were the weighting factors for the mixing ratio at , and the summation was taken for . The weighting factors were chosen to satisfy the condition of a zero pseudo-Laplacian:

where are the coordinates of in the local *x–y* coordinate system, the origin of which is . The weighting factors were written as

and were determined by solving an optimization problem in which the cost function

was minimized subject to the constraint in Eqs. (26), where . With the use of Lagrangian multipliers, this problem is equivalent to the problem of minimizing the following function:

leading to

The two Lagrangian multipliers are obtained from

where , , , , and .

With the interpolated values at , we reconstructed a piecewise parabolic profile,

inside the 0th triangular cell (Fig. 3b). Assuming , where are the coordinates of in the local *x–y* coordinate system about , we obtained the following matrix equation:

where

Note that Eq. (28) is slightly changed from Eq. (17), and the coefficient vector , which has six components, is determined by . The reconstructed profile satisfies for the six surrounding nodes and is an exact fit.

Similarly to the case of a hexagonal–pentagonal mesh, we adjusted Eq. (27) by adding a shift term as follows:

Application of the integral constraint,

where is the area of the 0th triangular cell, leads to

where are the coordinates of the center of the *i*th cell face, and . Note that means for . With the shifted quadratic profile, the mixing ratio of the *i*th cell face is determined as follows:

where and are the coordinates of and , respectively, in the local coordinate system.

### e. Computational cost and implicit diffusion

The type-I scheme might appear more complicated than the LM09 scheme and the type-II scheme. For all three of the schemes, however, the determination of in Fig. 4, for example, can be modified into the same form as follows:

where and

Here, wind direction across the cell face is not assumed. The values determined from the left-hand side and right-hand side of the cell face can be written, respectively, as

We can compute and store the factors and before the time loop, and the computational costs of the LM09 scheme, the type-I scheme, and the type-II scheme are thus the same.

In a manner similar to SG11, we were able to reduce the influence of the implicit numerical diffusion of the type-I and type-II schemes by changing the parameter *β* in Eq. (33). The type-I and type-II schemes are equivalent to the standard third-order scheme for *β* = 1 and to the nominally fourth-order scheme for *β* = 0.

## 3. Test results

A cosine-bell advection test (Williamson et al. 1992), a slotted-cylinder advection test (Lipscomb and Ringler 2005), and a deformational flow test (Nair and Lauritzen 2010) were performed for a hexagonal–pentagonal mesh to examine the properties of the type-I and the type-II schemes. We also conducted the same set of tests for the type-I scheme on a triangular mesh. The spherical icosahedral grid was generated by the iterative splitting method of Heikes and Randall (1995), and its shape was modified slightly by the spring method (Tomita et al. 2001; Miura and Kimoto 2005). This generated grid was composed of triangular cells and was used as a triangular mesh. From the triangular mesh, a hexagonal–pentagonal mesh was generated by defining the vertices of hexagons or pentagons at the barycenter of each triangular cell. We used grids having 2562, 10 242, 40 962, 16 382, and 655 362 node points over the globe. Following Miura and Skamarock (2013), we designated them as H16, H32, H64, H128, and H256 based on the number of node points on each edge of the original icosahedron. The RK3 scheme was adopted for temporal integration, and the flux limiter of Thuburn (1995, 1996) was used, with a slight modification (the appendix) to ensure monotonicity of the solutions where not explicitly noted. Accuracy of the schemes was measured by the and norms, the definitions of which are

where and are the computed and exact solutions on the *i*th cell node, respectively, and is the area of the *i*th cell.

### a. Advection of a cosine-shaped bell

The cosine-bell advection test proposed by Williamson et al. (1992) has been widely used to evaluate transport schemes on icosahedral grids (e.g., Heikes and Randall 1995; Thuburn 1997; Tomita et al. 2001; Lipscomb and Ringler 2005; Yeh 2007; Miura 2007; Mittal et al. 2007; LM09; Skamarock and Menchaca 2010; SG11; Miura and Skamarock 2013). Because our computed solutions were only weakly sensitive to the flow angle, we set the angle between the axis of the solid body rotation and the pole axis of the sphere to zero. For the H16, H32, H64, H128, and H256 grids, the time intervals were 7200, 3600, 1800, 900, and 450 s, respectively, for the hexagonal–pentagonal mesh and 3600, 1800, 900, 450, and 225 s, respectively, for the triangular mesh. The Courant number is 0.5 approximately in each case. We confirmed that the schemes were stable up to the Courant number of about unity although solutions were more diffusive than the solutions of the Courant number of 0.5 (not shown).

Figures 5a and 5d appear almost identical to Figs. 3b and 3a, respectively, of Lee et al. (2010) except for the difference in flow angle. There were nonnegligible phase errors in both cases, both with and without the flux limiter. This result confirms that our implementation of the LM09 scheme is consistent with LM09 and Lee et al. (2010). Figure 6 indicates that the *l*_{2} and *l*_{∞} norms of the LM09 scheme were also very close to those shown in Fig. 2 of Lee et al. (2010).

The solutions and error norms of the type-I and type-II schemes were almost equivalent and were better than the LM09 scheme counterparts. The phase errors were significantly smaller than the LM09 scheme phase errors (Fig. 5). The convergence rates of the *l*_{2} and *l*_{∞} norms were close to second order without the flux limiter (Fig. 6a). In contrast, with the flux limiter the *l*_{∞} norm converged slower than second order, whereas the *l*_{2} norm converged nearly to third order (Fig. 6b). For this test several researchers have likewise reported convergence rates faster than second order with flux limiters (e.g., Miura 2007; SG11). The solutions and error norms of the type-I and type-II schemes appeared very similar to those shown in Figs. 3 and 4 for the third-order scheme of SG11. This result confirms that the three different algorithms (SG11, the type-I, and the type-II) for extending the standard third-order scheme to the spherical icosahedral grid can reproduce almost equivalent results. In Eq. (22), we used the two Gauss quadrature points, **G**_{i±}, instead of the midpoint of the cell face, **F**_{i}, to compute the line integral. If we determined the face value simply from the equation (the type-I´ scheme in Figs. 6b and 7), the phase error enlarged, and accuracy was compromised. The two Gauss quadrature points are therefore desirable, although this requirement increases the complexity of preprocessing. Figure 8 shows that the results of the type-I scheme were as good on a triangular mesh as on a hexagonal–pentagonal mesh. The phase error was smaller, and the convergence of the error norms was significantly faster, than the analogous results in Figs. 4 and 5, respectively, of Walko and Avissar (2008a).

### b. Advection of a slotted cylinder

We performed the slotted-cylinder advection test of Lipscomb and Ringler (2005). The cosine bell in the previous test was replaced with a slotted cylinder. This test is harder than the cosine-bell advection test because it includes sharp discontinuities around the cylinder. The grid used and the time intervals were the same as those used in the cosine-bell advection test.

In Fig. 9, the solution of the LM09 scheme for an H64 grid was very similar to the analogous results in Fig. 4 of Lee et al. (2010). The solutions for the type-I and type-II schemes were almost equivalent again. The distortion of the solution was smaller for both the type-I and the type-II schemes than for the LM09 scheme. The differences in the error norms, however, were minor between the type-I or type-II scheme and the LM09 scheme, and the convergence rates of the *l*_{2} norm were less than first order for all of the schemes (Fig. 10). These results reflect the fact that most of the error comes from insufficient representation of the discontinuities around the cylinder. The type-I scheme performs equally well on a triangular mesh and hexagonal mesh (Fig. 11).

### c. Advection by deformational flow

We performed a test with a nondivergent deformational flow as proposed by Nair and Lauritzen (2010). The initial scalar and velocity fields were given by Eqs. (14), (31), and (32) in Nair and Lauritzen (2010). The parameter in the equations of the flow velocities was defined as , following Harris et al. (2011), where is the radius of the sphere and the period is 12 days. For the H16, H32, H64, H128, and H256 grids, the time intervals were 3600, 1800, 900, 450, and 225 s, respectively, for the hexagonal–pentagonal mesh, and 1800, 900, 450, 225, and 100 s, respectively, for the triangular mesh.

The initial pair of Gaussian hills (Fig. 12a) was strongly deformed until day 6 (Fig. 12b), a result that appears similar to the reference solution of Nair and Lauritzen (2010). Because the solutions for the type-I and type-II schemes were very similar on day 12, that comparison is not included in Fig. 12. The solution was less distorted for the type-I scheme than the LM09 scheme after the 12-day integrations. Without the flux limiter, the convergence rates of the error norms approached third order for the type-I and the type-II schemes, whereas the convergence rates were nearly second order for the LM09 scheme (Fig. 13a). With the flux limiter, the convergence of the *l*_{∞} norm was slower than second order, but the convergence of the *l*_{2} norm was still nearly third order (Fig. 13b). We found that the type-I and the type-II schemes were more accurate than the LM09 scheme for this test. The fact that the type-I scheme produced better results on a triangular mesh than a hexagonal–pentagonal mesh (Fig. 14) reflects the fact that the time intervals were almost halved and that the triangular mesh had almost twice the number of cells as the hexagonal–pentagonal mesh for the same grid.

### d. Reducing implicit numerical diffusion

By changing *β* in Eq. (33), we were able to reduce the influence of implicit numerical diffusion in the type-I and the type-II schemes, as was demonstrated by SG11. It is obvious from a comparison of Fig. 15 to Figs. 5 and 9 that the *β* = 0.25 case is significantly less diffusive than the original (*β* = 1.0). In compensation for the smaller implicit diffusion, the symmetry of the solution in Fig. 15d is weakly violated relative to Fig. 9b. For *β* = 0.0, the solution is distorted in the rear side of the cosine bell (Fig. 15b), or else the solution becomes noisy (Fig. 15e). These defects with *β* = 0.0 have already been shown in Figs. 3 and 12 of SG11. We can mitigate this negative effect by making the time interval smaller. If the time interval is halved from 3600 to 1800 s, the benefit from the nominally fourth-order scheme emerges (Figs. 15c,f). The solutions are even less diffusive than the *β* = 0.25 solutions, but the asymmetry about the 90°W line is larger than is the case for *β* = 1.0 and *β* = 0.25. As suggested by SG11, *β* = 0.25 seems to be a reasonable choice if a weaker implicit diffusion is required.

## 4. Summary

In this work we proposed two new algorithms to extend the standard third-order upwind-biased transport scheme to a spherical icosahedral grid. In the first algorithm, the quadratic profile of the tracer was reconstructed inside the upwind-side cell [Eq. (16) for hexagonal or pentagonal cells or Eq. (27) for triangular cells]. The constructed profile was then shifted to satisfy the integral constraint with respect to the amount of tracer [Eq. (20) for hexagonal or pentagonal cells or Eq. (30) for triangular cells]. By using the shifted quadratic profile [Eq. (19) for hexagonal or pentagonal cells or Eq. (29) for triangular cells], we determined the cell face value with Eq. (22) for hexagonal or pentagonal cells or with Eq. (31) for triangular cells. We characterized this algorithm as type I and tested it on both hexagonal–pentagonal and triangular meshes. The second algorithm involved determination of the cell face value from a weighted mean in Eq. (25) of a linearly interpolated value in Eq. (23) and a linearly extrapolated value in Eq. (24). We characterized this algorithm as type II and tested it only on the hexagonal–pentagonal mesh. The type-II algorithm was not applied to a triangular mesh because the extrapolation from the upwind side was difficult to define. Both the type-I and type-II schemes can be used with multistage time integration schemes like RK3 in a manner similar to that used for the SG11 scheme, which is another extension of the standard third-order scheme.

The type-I and type-II schemes were subjected to the cosine-bell advection test, the slotted-cylinder advection test, and the deformational flow test and were compared with the LM09 scheme and the SG11 scheme. On a hexagonal–pentagonal mesh, the results of the type-I scheme, the type-II scheme, and the SG11 scheme were equivalent. The solutions were very similar, and the error norms were almost the same. Those schemes were obviously superior to the LM09 scheme, at least in the tests conducted in this study. The type-I scheme worked on a triangular mesh as well as on a hexagonal–pentagonal mesh. On a triangular mesh the type-I scheme produced a more accurate solution with a smaller phase error in the cosine-bell advection test compared to the results of Walko and Avissar (2008a). By changing the parameter β in Eq. (33), we could reduce the influence of implicit numerical diffusion. With respect to the cosine-bell advection test and the slotted-cylinder advection test, *β* = 0.25 seemed a reasonable choice, as was suggested by SG11, because decreasing *β* led to increasing distortion of the solution. The computational costs of the temporal integration were the same for the type-I, type-II, and LM09 schemes because the determination of the cell face value could be summarized by the same formulations in Eqs. (32)–(34).

The type-I scheme can be more broadly applied than the type-II scheme and the SG11 scheme because it assumes only that the shape of a cell is hexagonal, pentagonal, or triangular on a spherical icosahedral grid. In contrast, the type-II scheme assumes symmetry about the cell center, and the SG11 scheme assumes a Voronoi mesh. As a result, the type-I scheme can be used on non-Voronoi meshes and, in addition, may be preferable when a local mesh refinement is adopted (e.g., Tomita 2008). The algorithm for the type-I scheme, which is based on a quadratic profile reconstruction subject to an integral constraint on the amount of tracer, will possibly be useful when we extend the standard third-order scheme to three-dimensional space.

## Acknowledgments

Hiroaki Miura thanks Prof. David Randall for supporting his visiting to Colorado State University; a part of this work was done during that visit. This work was supported by the Grant-in-Aid for Young Scientists (B) of Ministry of Education, Culture, Sports, Science, and Technology of Japan (Grant 22740310).

### APPENDIX

#### Flux Limiter

We used a flux limiter to ensure the monotonicity of solutions. The original method of Thuburn (1995) and Thuburn (1996) was slightly modified to allow less diffusive solutions. For the cell of Fig. 3a, for example, the first guess of the cell face mixing ratio, , which is determined by Eq. (32), was adjusted in the following procedure. Note that means for and means for .

- For each cell, outflow bounds were calculated: where When there were no outflow boundaries, and were set to some very large negative and positive numbers, respectively. If the density was not constant, we modified as follows: where is the density of the cell, and is the density estimated at the
*i*th cell face.

## REFERENCES

*Numerical Methods for Fluid Dynamics: With Applications to Geophysics.*2nd ed. Springer, 516 pp.

*Proc. AIAA Ninth Computational Fluid Dynamics Conf.,*Buffalo, NY, AIAA, 25–39.

*Proc. 16th AIAA Computational Fluid Dynamics Conf.,*Orlando, FL, AIAA, 1–11.