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

 
formula

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.

Fig. 1.

Schematic illustrations of a hexagonal cell in (a) a Voronoi mesh and (b) a non-Voronoi mesh. The triangular mesh shown by the dotted lines is the same.

Fig. 1.

Schematic illustrations of a hexagonal cell in (a) a Voronoi mesh and (b) a non-Voronoi mesh. The triangular mesh shown by the dotted lines is the same.

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:

 
formula

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:

 
formula

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

 
formula

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

 
formula

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

 
formula

where the cell interface values and are determined as follows:

 
formula

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 x1/2 and inflow to the 0th cell through the interface at x−1/2 by the grid interval.

Fig. 2.

(a) A tracer transport problem on a one-dimensional grid. (b) The determination of the cell interface value by using the quadratic profile , which is obtained by shifting the quadratic profile , that is reconstructed inside the 0th cell so that it passes through , , and . (c) The determination of the cell interface value from the weighted mean of the interpolated value and the extrapolated value .

Fig. 2.

(a) A tracer transport problem on a one-dimensional grid. (b) The determination of the cell interface value by using the quadratic profile , which is obtained by shifting the quadratic profile , that is reconstructed inside the 0th cell so that it passes through , , and . (c) The determination of the cell interface value from the weighted mean of the interpolated value and the extrapolated value .

SG11 modified Eq. (4) as follows:

 
formula

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

 
formula

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

 
formula

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 x0 is approximated by a quadratic profile (the gray line in Fig. 2b):

 
formula

Using the constraints f(−Δx) = q−1, f(0) = q0 and fx) = q1, the coefficients a0, a1, and a2 are determined, and Eq. (7) becomes

 
formula

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

 
formula

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:

 
formula

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

 
formula

This constraint causes the shifted tracer profile integrated over the x0 cell to be equal to the tracer amount inside the cell (if q0 is the mean in the cell). Substituting Eq. (10) into Eq. (11) causes the form of to become

 
formula

and the shifted quadratic profile is specified as follows:

 
formula

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

 
formula

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

 
formula

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

 
formula

It is straightforward to show that

 
formula

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:

 
formula

where qi(i = 0, 1, … , N0) is the mixing ratio at node point Hi (Fig. 3a) or Ti (Fig. 3b), A0 is the area of the 0th cell, N0 is the number of the cells surrounding the 0th cell (N0 = 6 or 5 for the hexagonal–pentagonal mesh and N0 = 3 for the triangular mesh), li (i = 1 , … , N0) is the length of the ith cell face shared by the 0th cell and the ith surrounding cell, ni is the unit vector normal to the ith cell face and pointing outward from the 0th cell, and and are the mixing ratio and velocity vector, respectively, representing those quantities for the ith cell face. If density is not constant, the right-hand side of Eq. (1) might be discretized as follows:

 
formula

where is the density of the ith 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 ith cell face ).

Fig. 3.

Schematic illustrations of (a) a hexagonal cell and (b) a triangular cell. (c) A hexagonal cell and (d) a pentagonal cell where the gray-shaded triangles are used to determine gradients for interpolation and extrapolation.

Fig. 3.

Schematic illustrations of (a) a hexagonal cell and (b) a triangular cell. (c) A hexagonal cell and (d) a pentagonal cell where the gray-shaded triangles are used to determine gradients for interpolation and extrapolation.

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 ith cell face is estimated as

 
formula

where Δxi is the distance between H0 and Hi (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 H0 and Hi . SG11 showed that this scheme worked very well on a Voronoi mesh. For non-Voronoi meshes, however, an ambiguity in specifying Δxi 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.

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 H0, that is,

 
formula

We use least squares on the local two-dimensional xy plane about H0 to determine the unknown coefficients aj (j = 1, … , 5) (Stuhne and Peltier 1996; Majewski et al. 2002). The coordinates of Hi on the local plane may be specified by , where nx and ny constitute an orthonormal basis set. By assuming qi = f(xi, yi) (i = 1, … , N0) for the surrounding cells, we obtain the following matrix equation:

 
formula

where

 
formula

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 ith cell face as follows:

 
formula

where are the coordinates of Fi = (Ti−1 + Ti)/2 on the local coordinate system. Note that Ti−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:

 
formula

and the integral constraint

 
formula

is used to obtain

 
formula

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

 
formula

where αi is the area of the subtriangle ΔH0Ti−1Ti, and and are the coordinates of Fi = (Ti−1 + Ti)/2 and Ei = (H0 + Ti)/2, respectively, in the local coordinate system (Fig. 3a).

Substitution of into Eq. (19) allows one to estimate the mixing ratio of the ith cell face as in a manner similar to Eq. (18). This method, however, is not accurate enough because fc 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 fc. Substituting the coordinates of into Eq. (19), the mixing ratio of the ith cell face is determined as follows:

 
formula

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 ith cell face can be written in a two-dimensional x–y coordinate system about H0 as

 
formula
 
formula

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

 
formula

The gradient terms, and , are defined in terms of the interpolation to F1, for example, as

 
formula

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

 
formula

for a pentagonal cell (Fig. 3d). The gradient of the ith triangular cell, ΔH0HiHi+1, is given by

 
formula

where

 
formula

Here, nx and ny are the orthonormal basis vectors of the local x–y coordinate system about Ti, and (x0, y0), (xi, yi), and (xi+1, yi+1) are the coordinates of H0, Hi, and Hi+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

 
formula

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:

 
formula

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

 
formula

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

 
formula

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:

 
formula

leading to

 
formula

The two Lagrangian multipliers are obtained from

 
formula

where , , , , and .

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

 
formula

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:

 
formula

where

 
formula

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:

 
formula

Application of the integral constraint,

 
formula

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

 
formula

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

 
formula

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:

 
formula
 
formula

where and

 
formula

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

 
formula

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.

Fig. 4.

A schematic illustration of the determination of the cell face value from the left-hand side and from the right-hand side.

Fig. 4.

A schematic illustration of the determination of the cell face value from the left-hand side and from the right-hand side.

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

 
formula

where and are the computed and exact solutions on the ith cell node, respectively, and is the area of the ith 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 l2 and l norms of the LM09 scheme were also very close to those shown in Fig. 2 of Lee et al. (2010).

Fig. 5.

The cosine-bell advection test on the H32 hexagonal–pentagonal mesh. The solutions after the 12-day integration. (a) LM09, (b) type-I, and (c) type-II schemes without the flux limiter. (d) LM09, (e) type-I, and (f) type-II schemes with the flux limiter. The true solution is indicated by the dotted lines in (d)–(f).

Fig. 5.

The cosine-bell advection test on the H32 hexagonal–pentagonal mesh. The solutions after the 12-day integration. (a) LM09, (b) type-I, and (c) type-II schemes without the flux limiter. (d) LM09, (e) type-I, and (f) type-II schemes with the flux limiter. The true solution is indicated by the dotted lines in (d)–(f).

Fig. 6.

Dependence of the and norms on grid refinement in the cosine-bell advection test. These norms were evaluated after a 12-day integration. The numbers inside the parentheses on the abscissa are the corresponding glevels of Tomita et al. (2001). (a) Without the flux limiter and (b) with the flux limiter. The thick solid and solid lines are reference lines for third-order and the second-order convergences, respectively.

Fig. 6.

Dependence of the and norms on grid refinement in the cosine-bell advection test. These norms were evaluated after a 12-day integration. The numbers inside the parentheses on the abscissa are the corresponding glevels of Tomita et al. (2001). (a) Without the flux limiter and (b) with the flux limiter. The thick solid and solid lines are reference lines for third-order and the second-order convergences, respectively.

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 l2 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 l2 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, Gi±, instead of the midpoint of the cell face, Fi, 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).

Fig. 7.

Computed solution minus true solution on the H256 hexagonal–pentagonal mesh. The type-I scheme (a) using two Gauss quadrature points or (b) using the midpoint to estimate the line integral for each cell face.

Fig. 7.

Computed solution minus true solution on the H256 hexagonal–pentagonal mesh. The type-I scheme (a) using two Gauss quadrature points or (b) using the midpoint to estimate the line integral for each cell face.

Fig. 8.

The solutions of the cosine-bell advection test on the H32 triangular mesh (a) without the flux limiter, and (b) with the flux limiter. (c) The and norms. The norms of the type-I scheme with the flux limiter on the hexagonal mesh (HEX) have been added as references. The thick solid and the solid lines are reference lines for third-order and the second-order convergences, respectively.

Fig. 8.

The solutions of the cosine-bell advection test on the H32 triangular mesh (a) without the flux limiter, and (b) with the flux limiter. (c) The and norms. The norms of the type-I scheme with the flux limiter on the hexagonal mesh (HEX) have been added as references. The thick solid and the solid lines are reference lines for third-order and the second-order convergences, respectively.

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 l2 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).

Fig. 9.

The slotted-cylinder advection test. The solutions after a 12-day integration. (a) LM09, (b) type-I, and (c) type-II schemes on an H32 hexagonal–pentagonal mesh. (d) LM09, (e) type-I, and (f) type-II schemes on an H64 hexagonal–pentagonal mesh.

Fig. 9.

The slotted-cylinder advection test. The solutions after a 12-day integration. (a) LM09, (b) type-I, and (c) type-II schemes on an H32 hexagonal–pentagonal mesh. (d) LM09, (e) type-I, and (f) type-II schemes on an H64 hexagonal–pentagonal mesh.

Fig. 10.

As in Fig. 6, but for the slotted-cylinder advection test. The solid line is a reference line for first-order convergence.

Fig. 10.

As in Fig. 6, but for the slotted-cylinder advection test. The solid line is a reference line for first-order convergence.

Fig. 11.

The solutions of the slotted-cylinder advection test on (a) H32 and (b) H64 triangular meshes. (c) As in Fig. 10, but for the triangular mesh. The norms of the type-I scheme on the hexagonal mesh (HEX) have been added as a reference.

Fig. 11.

The solutions of the slotted-cylinder advection test on (a) H32 and (b) H64 triangular meshes. (c) As in Fig. 10, but for the triangular mesh. The norms of the type-I scheme on the hexagonal mesh (HEX) have been added as a reference.

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 l2 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.

Fig. 12.

The deformational flow test. (a) The initial and true solutions. (b) The solution of the type-I scheme after a 6-day integration on the H64 hexagonal–pentagonal mesh. True (dotted lines) and computed (solid lines) solutions of the (c) LM09 and (d) type-I schemes after a 12-day integration on the H64 hexagonal–pentagonal mesh.

Fig. 12.

The deformational flow test. (a) The initial and true solutions. (b) The solution of the type-I scheme after a 6-day integration on the H64 hexagonal–pentagonal mesh. True (dotted lines) and computed (solid lines) solutions of the (c) LM09 and (d) type-I schemes after a 12-day integration on the H64 hexagonal–pentagonal mesh.

Fig. 13.

As in Fig. 6, but for the deformational flow test.

Fig. 13.

As in Fig. 6, but for the deformational flow test.

Fig. 14.

(a) As in Figs. 12c,d, but for the H64 triangular mesh. (b) As in Fig. 8c, but for the deformational flow test.

Fig. 14.

(a) As in Figs. 12c,d, but for the H64 triangular mesh. (b) As in Fig. 8c, but for the deformational flow test.

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.

Fig. 15.

The solutions of the type-I scheme on an H32 hexagonal–pentagonal mesh for (a)–(c) the cosine-bell advection test and for (d)–(f) the slotted-cylinder advection test. The parameter β was set to (a),(d) 0.25 and (b),(c),(e),(f) 0.0 . The time interval was (a),(b),(d),(e) 3600 and (c),(f) 1800 s.

Fig. 15.

The solutions of the type-I scheme on an H32 hexagonal–pentagonal mesh for (a)–(c) the cosine-bell advection test and for (d)–(f) the slotted-cylinder advection test. The parameter β was set to (a),(d) 0.25 and (b),(c),(e),(f) 0.0 . The time interval was (a),(b),(d),(e) 3600 and (c),(f) 1800 s.

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 .

  1. For each face, inflow bounds were calculated: 
    formula
     
    formula
  2. For each face, the first guess was adjusted so that it lay within the inflow bounds: 
    formula
  3. For each face, the inflow bounds were revised: 
    formula
     
    formula
     
    formula
  4. For each cell, the minimum and maximum permitted values at the next step are given by the smaller and larger values between and the minimum and maximum of the inflow bounds over the inflow faces: 
    formula
     
    formula
  5. For each cell, outflow bounds were calculated: 
    formula
     
    formula
    where 
    formula
    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: 
    formula
    where is the density of the cell, and is the density estimated at the ith cell face.
  6. For each face, the revised first guess was adjusted to lie within the outflow bounds: 
    formula

REFERENCES

REFERENCES
Durran
,
D. R.
,
2010
: Numerical Methods for Fluid Dynamics: With Applications to Geophysics. 2nd ed. Springer, 516 pp.
Gassmann
,
A.
,
2013
:
A global hexagonal C-grid non-hydrostatic dynamical core (ICON-IAP) designed for energetic consistency
.
Quart. J. Roy. Meteor. Soc.
,
139
,
152
175
.
Golub
,
G. H.
, and
C.
Reinsch
,
1970
:
Singular value decomposition and least squares solutions
.
Numer. Math.
,
14
,
403
420
.
Harris
,
L. M.
,
P. H.
Lauritzen
, and
R.
Mittal
,
2011
:
A flux-form version of the conservative semi-Lagrangian multi-tracer transport scheme (CSLAM) on the cubed sphere grid
.
J. Comput. Phys.
,
230
,
1215
1237
.
Heikes
,
R.
, and
D. A.
Randall
,
1995
:
Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part I: Basic design and results of tests
.
Mon. Wea. Rev.
,
123
,
1862
1880
.
Holmes
,
D. G.
, and
S. D.
Connel
,
1989
: Solution of the two-dimensional Navier–Stokes equations on unstructured adaptive grids. AIAA Paper 89-1932, Proc. AIAA Ninth Computational Fluid Dynamics Conf., Buffalo, NY, AIAA, 25–39.
Ii
,
S.
, and
F.
Xiao
,
2010
:
A global shallow water model using high order multi-moment constrained finite volume method and icosahedral grid
.
J. Comput. Phys.
,
229
,
1774
1796
.
Ishida
,
J.
,
C.
Muroi
,
K.
Kawano
, and
Y.
Kitamura
,
2010
: Development of a new nonhydrostatic model ASUCA at JMA. WGNE Blue Book: Research activities in atmosphere and ocean modeling, WMO, 5.11–5.12. [Available online at http://www.wcrp-climate.org/WGNE/BlueBook/2010/documents/author-list.html.]
Kim
,
S.-E.
,
B.
Makarov
, and
D.
Caraeni
,
2003
: A multi-dimensional linear reconstruction scheme for arbitrary unstructured grids. AIAA Paper 2003-3990, Proc. 16th AIAA Computational Fluid Dynamics Conf., Orlando, FL, AIAA, 1–11.
Lee
,
J.-L.
, and
A. E.
MacDonald
,
2009
:
A finite-volume icosahedral shallow-water model on a local coordinate
.
Mon. Wea. Rev.
,
137
,
1422
1437
.
Lee
,
J.-L.
,
R.
Bleck
, and
A. E.
MacDonald
,
2010
:
A multistep flux-corrected transport scheme
.
J. Comput. Phys.
,
229
,
9284
9298
.
Lipscomb
,
W. H.
, and
T. D.
Ringler
,
2005
:
An incremental remapping transport scheme on a spherical geodesic grid
.
Mon. Wea. Rev.
,
133
,
2335
2350
.
Majewski
,
D.
, and
Coauthors
,
2002
:
The operational global icosahedral-hexagonal gridpoint model GME: Description and high-resolution tests
.
Mon. Wea. Rev.
,
130
,
319
338
.
Mittal
,
R.
,
H. C.
Upadhyaya
, and
O. P.
Sharma
,
2007
:
On near-diffusion-free advection over spherical geodesic grids
.
Mon. Wea. Rev.
,
135
,
4214
4225
.
Miura
,
H.
,
2007
:
An upwind-biased conservative advection scheme for spherical hexagonal–pentagonal grids
.
Mon. Wea. Rev.
,
135
,
4038
4044
.
Miura
,
H.
, and
M.
Kimoto
,
2005
:
A comparison of grid quality of optimized spherical hexagonal–pentagonal geodesic grids
.
Mon. Wea. Rev.
,
133
,
2817
2833
.
Miura
,
H.
, and
W. C.
Skamarock
,
2013
:
An upwind-biased transport scheme using a quadratic reconstruction on spherical icosahedral grids
.
Mon. Wea. Rev.
,
141
,
832
847
.
Miura
,
H.
,
M.
Satoh
,
T.
Nasuno
,
A. T.
Noda
, and
K.
Oouchi
,
2007
:
A Madden-Julian oscillation event realistically simulated by a global cloud-resolving model
.
Science
,
318
,
1763
1765
.
Nair
,
R. D.
, and
P. H.
Lauritzen
,
2010
:
A class of deformational flow test cases for linear transport problems on the sphere
.
J. Comput. Phys.
,
229
,
8868
8887
.
Ringler
,
T. D.
,
R. P.
Heikes
, and
D. A.
Randall
,
2000
:
Modeling the atmospheric general circulation using a spherical geodesic grid: A new class of dynamical cores
.
Mon. Wea. Rev.
,
128
,
2471
2490
.
Satoh
,
M.
,
T.
Matsuno
,
H.
Tomita
,
H.
Miura
,
T.
Nasuno
, and
S.
Iga
,
2008
:
Nonhydrostatic icosahedral atmospheric model (NICAM) for global cloud resolving simulations
.
J. Comput. Phys.
,
227
,
3484
3514
.
Skamarock
,
W. C.
, and
J. B.
Klemp
,
2008
:
A time-split nonhydrostatic atmospheric model for weather research and forecasting applications
.
J. Comput. Phys.
,
227
,
3465
3485
.
Skamarock
,
W. C.
, and
M.
Menchaca
,
2010
:
Conservative transport schemes for spherical geodesic grids: High-order reconstructions for forward-in-time schemes
.
Mon. Wea. Rev.
,
138
,
4497
4508
.
Skamarock
,
W. C.
, and
A.
Gassmann
,
2011
:
Conservative transport schemes for spherical geodesic grids: High-order flux operators for ODE-based time integration
.
Mon. Wea. Rev.
,
139
,
2962
2975
.
Skamarock
,
W. C.
,
J. B.
Klemp
,
M. G.
Duda
,
L. D.
Fowler
,
S.-H.
Park
, and
T. D.
Ringler
,
2012
:
A multiscale nonhydrostatic atmospheric model using centroidal Voronoi tessellations and C-grid staggering
.
Mon. Wea. Rev.
,
140
,
3090
3105
.
Stuhne
,
G. R.
, and
W. R.
Peltier
,
1996
:
Vortex erosion and amalgamation in a new model of large scale flow on the sphere
.
J. Comput. Phys.
,
128
,
58
81
.
Thuburn
,
J.
,
1995
:
Dissipation and cascades to small scales in numerical models using a shape-preserving advection scheme
.
Mon. Wea. Rev.
,
123
,
1888
1903
.
Thuburn
,
J.
,
1996
:
Multidimensional flux-limited advection schemes
.
J. Comput. Phys.
,
123
,
74
83
.
Thuburn
,
J.
,
1997
:
A PV-based shallow-water model on a hexagonal–icosahedral grid
.
Mon. Wea. Rev.
,
125
,
2328
2347
.
Tomita
,
H.
,
2008
:
A stretched grid on a sphere by new grid transformation
.
J. Meteor. Soc. Japan
,
86A
,
107
119
.
Tomita
,
H.
, and
M.
Satoh
,
2004
:
A new dynamical framework of nonhydrostatic global model using the icosahedral grid
.
Fluid Dyn. Res.
,
34
,
357
400
.
Tomita
,
H.
,
M.
Tsugawa
,
M.
Satoh
, and
K.
Goto
,
2001
:
Shallow water model on a modified icosahedral geodesic grid by using spring dynamics
.
J. Comput. Phys.
,
174
,
579
613
.
Tomita
,
H.
,
H.
Miura
,
S.
Iga
,
T.
Nasuno
, and
M.
Satoh
,
2005
:
A global cloud-resolving simulation: Preliminary results from an aqua planet experiment
.
Geophys. Res. Lett.
,
32
, L08805, doi:10.1029/2005GL022459.
Walko
,
R. L.
, and
R.
Avissar
,
2008a
:
The Ocean–Land–Atmosphere Model (OLAM). Part I: Shallow-water tests
.
Mon. Wea. Rev.
,
136
,
4033
4044
.
Walko
,
R. L.
, and
R.
Avissar
,
2008b
:
The Ocean–Land–Atmosphere Model (OLAM). Part II: Formulation and tests of the nonhydrostatic dynamic core
.
Mon. Wea. Rev.
,
136
,
4045
4062
.
Wan
,
H.
, and
Coauthors
,
2013
:
The ICON-1.2 hydrostatic atmospheric dynamical core on triangular grids—Part 1: Formulation and performance of the baseline version
.
Geosci. Model Dev. Discuss.
,
6
,
59
119
.
Weller
,
H.
,
H. G.
Weller
, and
A.
Fournier
,
2009
:
Voronoi, Delaunay, and block-structured mesh refinement for solution of the shallow-water equations on the sphere
.
Mon. Wea. Rev.
,
137
,
4208
4224
.
Wicker
,
L. J.
, and
W. C.
Skamarock
,
2002
:
Time-splitting methods for elastic models using forward time schemes
.
Mon. Wea. Rev.
,
130
,
2088
2097
.
Williamson
,
D. L.
,
J. B.
Drake
,
J. J.
Hack
,
R.
Jakob
, and
P. N.
Swarztrauber
,
1992
:
A standard test set for numerical approximations to the shallow water equations in spherical geometry
.
J. Comput. Phys.
,
102
,
211
224
.
Yeh
,
K.-S.
,
2007
:
The streamline subgrid integration method: I. Quasi-monotonic second-order transport schemes
.
J. Comput. Phys.
,
225
,
1632
1652
.