Abstract

The arbitrarily structured C grid, Thuburn–Ringler–Skamarock–Klemp (TRiSK), is being used in the Model for Prediction Across Scales (MPAS) and is being considered by the Met Office for their next dynamical core. However, the hexagonal C grid supports a branch of spurious Rossby modes, which lead to erroneous grid-scale oscillations of potential vorticity (PV). It is shown how these modes can be harmlessly controlled by using upwind-biased interpolation schemes for PV. A number of existing advection schemes for PV are tested, including that used in MPAS, and none are found to give adequate results for all grids and all cases. Therefore a new scheme is proposed; continuous, linear-upwind stabilized transport (CLUST), a blend between centered and linear-upwind with the blend dependent on the flow direction with respect to the cell edge.

A diagnostic of grid-scale oscillations is proposed that gives further discrimination between schemes than using potential enstrophy alone. Indeed, some schemes are found to destroy potential enstrophy while grid-scale oscillations grow.

CLUST performs well on hexagonal-icosahedral grids and unrotated skipped latitude–longitude grids of the sphere for various shallow-water test cases. Despite the computational modes, the hexagonal icosahedral grid performs well since these modes are easy and harmless to filter. As a result, TRiSK appears to perform better than a spectral shallow-water model.

1. Introduction

The formulation of the C grid on arbitrarily structured, two-dimensional orthogonal grids [Thuburn–Ringler–Skamarock–Klemp (TRiSK); Thuburn et al. 2009; Ringler et al. 2010] was an important breakthrough for weather and climate forecasting as this enables some of the desirable properties of the latitude–longitude C grid to hold on quasi-uniform grids of the sphere. (Orthogonal means that the line between cell centers is at right angles to the cell edges.) TRiSK ensures local conservation of mass and potential vorticity, steady geostrophic states, and conservation of energy as the time step tends to zero (but not potential enstrophy). These important properties on a compact and very cheap scheme are enabling models using quasi-uniform, Voronoi grids to compete with the accuracy of latitude–longitude grid (Weller et al. 2012), while removing time-step restrictions due to the convergence of the meridians toward the North and South Poles and thus enabling improved parallel scaling.

However no quasi-uniform C grid of the sphere seems ideal. If a 2D C grid is not constructed of quadrilateral cells and/or does not have exactly twice the number of velocity degrees of freedom (d.o.f.) as mass variable d.o.f., then computational modes will be present (Le Roux et al. 2005; Staniforth and Thuburn 2012). The C grid on hexagons has too many velocity d.o.f. and therefore suffers from computational modes consisting of a spurious branch of very slow Rossby waves in the dispersion relation between wavenumber and frequency (Thuburn 2008). The triangular C grid has too few velocity d.o.f. and therefore has two spurious branches of slow inertio-gravity modes (Danilov 2010). The cubed sphere has the correct ratio of d.o.f., but is usually nonorthogonal, apart from the conformal cubed sphere, which is not quasi-uniform since cells cluster toward the cube corners, increasingly with increasing resolution (Rančić et al. 1996; Putman and Lin 2007). Other quasi-uniform, orthogonal grids of the sphere have been tested using TRiSK (Weller et al. 2012), and using collocated finite volumes (Weller et al. 2009), but none are as accurate as the hexagonal icosahedron.

The computational modes on the hexagonal icosahedron consist of grid-scale oscillations of potential vorticity (PV) on the vertices or, equivalently, on the cells of the dual, triangular grid. An example of the impact of this mode can be seen in Fig. 1c, which shows the relative vorticity after 50 days for the Williamson et al. (1992) flow over a midlatitude mountain using a hexagonal grid of 40 962 cells and midpoint interpolation of PV from vertices to edges. This interpolation means that grid-scale PV oscillations can be present on the vertices, but not on the edges and so they are not present in the discretized solution of the momentum equation.

Fig. 1.

Relative vorticity after 50 days for the Williamson et al. (1992) flow over a mountain using various interpolation schemes for PV. Hexagonal grid of 40 962 cells. Time step is 100 s.

Fig. 1.

Relative vorticity after 50 days for the Williamson et al. (1992) flow over a mountain using various interpolation schemes for PV. Hexagonal grid of 40 962 cells. Time step is 100 s.

In TRiSK the advection is compatible so that the PV diagnosed from the normal velocity on each edge is identical to the PV as if it were advected on the dual grid. This means that higher-order or monotonic interpolation of PV from vertices to edges can be used to control grid-scale oscillations. This technique was introduced by Ringler et al. (2010) using an Anticipated Potential Vorticity Method (APVM) for use in the Model for Prediction Across Scales (MPAS) (http://mpas.sourceforge.net/) and will be described in section 3a. APVM acts like a smooth, upwind-biased advection and so grid-scale oscillations are reduced. However, the energy conserving averaging in TRiSK means that the implied PV advection on the dual grid would include this additional averaging and so theoretical results about any monotonic advection scheme do not hold. Notwithstanding, Ringler et al. (2010) found that APVM reduces grid-scale noise in the PV on the dual-triangular grid. However, closer inspection of the relative vorticity on the triangles using APVM reveals significant grid-scale noise in some test cases (Fig. 1i). This could either be noise introduced by the energy conserving averaging or an inadequacy of APVM. A scale-aware version of of APVM has been developed by Chen et al. (2011) in order to replace the time-step dependence with a dependence on the local turbulence with the aim of reproducing a realistic turbulent enstrophy cascade rather than removing computational modes. However, the computational modes are not well controlled using scale-aware APVM 1(Fig. 1j).

Grid-scale oscillations of PV on triangles are not controlled by the hexagonal C-grid discretization of Gassmann (2011) (A. Gassmann 2011, personal communication). Gassmann (2011) does not define PV on triangles, but on edges using a larger stencil than for TRiSK. However, the velocity patterns that lead to PV oscillations on triangles can still exist. The additional averaging means that they cannot be seen by the momentum equation and so they cannot be controlled. Conversely, in TRiSK, PV is defined in the most compact way, on the dual triangles and then averaged to the edges. If the averaging is done in an uncentered way (e.g., including some upwinding) then the discretization of the momentum equation will be able to see the grid-scale oscillations and act to remove them. Ongoing work by Gassmann controls computational modes better.

There has been some discussion whether the grid-scale oscillations of PV when using TRiSK on hexagons is a result of the computational modes or of the realistic turbulent cascade to smaller scales resulting in a buildup of vorticity at the grid scale, which must be removed since it cannot cascade farther downward to viscosity. It was argued by Weller et al. (2012) that any stationary grid-scale mode of the shallow-water equations on the rotating sphere must be a spurious computational artifact since the only true stationary modes are zonally symmetric. The nature of any grid-scale oscillations is additionally addressed here by comparing results of TRiSK on a hexagonal-icosahedron with those on a skipped latitude–longitude grid. Apart from the grid reductions, the skipped grid does not suffer from computational modes (besides the Coriolis mode of the C grid when the Rossby radius is poorly resolved; Randall 1994). However, a turbulent cascade to small scales will still result in a buildup of vorticity at the grid scale unless it is removed by some sort of diffusion or dissipation. Therefore, if this process is active then it should also affect results on the latitude–longitude grid.

The grids used are described in section 2. The numerical method is outlined in section 3 and the advection schemes are defined and some of their numerical properties reported. Results of 3 shallow-water test cases are presented in section 4 and final conclusions are drawn in section 5.

2. Quasi-uniform C grids of the sphere

Results are presented on two types of Voronoi tessellation of the sphere: the Heikes and Randall (1995) optimized hexagonal icosahedron (without the twist; Fig. 2a) grid generator from J. Thuburn (2004, personal communication), and a Voronoi tessellation of the cell centers of a skipped latitude–longitude grid (Fig. 2b). The skipped grid has 2:1 reductions in resolution toward the Pole when the cell-center to cell-center distance reaches less than 2.5 times that distance at the equator. For both grids, the Voronoi generating points are used as the cell centers, which implies that the grids are orthogonal. For all Voronoi tessellations, the edges are midway between the cell centers but the edge-crossing points are not midway along the edges. These grids are used as C grids, which implies that the mass variable is stored at the cell center and the velocity is normal to each edge, at the edge-crossing point.

Fig. 2.

Coarse versions of the quasi-uniform grids.

Fig. 2.

Coarse versions of the quasi-uniform grids.

The hexagonal grid was found to give consistently best results of the quasi-uniform grids presented by Weller et al. (2012), but the vorticity on the dual grid of triangles was subject to grid-scale noise. The skipped latitude–longitude grid behaves like a perfect latitude–longitude grid when there is no flow variability near the 2:1 grid resolution reductions. This is therefore a useful benchmark for assessing the results of the hexagonal grids. Various resolutions of each grid are used as shown in Table 1.

Table 1.

Statistics of the grids of the sphere. The Courant numbers are the maximum over the grid for the initial conditions. They are based on the flow speed and the gravity wave (g.w.) speed.

Statistics of the grids of the sphere. The Courant numbers are the maximum over the grid for the initial conditions. They are based on the flow speed and the gravity wave (g.w.) speed.
Statistics of the grids of the sphere. The Courant numbers are the maximum over the grid for the initial conditions. They are based on the flow speed and the gravity wave (g.w.) speed.

3. Numerical method

We are solving the shallow-water equations (SWEs) on the rotating sphere with the momentum equation in vector-invariant form and the continuity equation in flux form:

 
formula
 
formula

where h is the fluid depth, u is the horizontal wind, b is the height of the bottom boundary, g is the acceleration due to gravity, is the PV, u = k × u, f is the Coriolis parameter, k is the local outward-pointing unit normal to the sphere, and is the kinetic energy. These are solved using TRiSK (Thuburn et al. 2009; Ringler et al. 2010) on the orthogonal C grids described in section 2. The prognostic variables are the mass variable at cell centers (the height, h) and the velocity component normal to each cell edge (un). Only the aspects of the numerical method that are different from Weller et al. (2012) are described here. The various interpolation techniques tested are described in section 3a and the semi-implicit time stepping is described in section 3b.

a. Interpolation of PV from vertices to edges

APVM is a forward-in-time advection scheme; the space and time discretization are treated together by approximating the value on an edge by the quantity that passes through the edge in one time step. To second order, this can be approximated by the value at a departure point. The other advection schemes are split between space and time so the interpolation onto the edge is instantaneous and the time stepping is done separately with the Crank–Nicolson ordinary differential equation (ODE) solver.

These interpolation schemes are described with reference to the dual cells (rather than the primal vertices) onto the edge-crossing point. We are therefore interpolating from the upwind and downwind dual cells, with centers at xu and xd, onto the point xe on the edge between them. These components are shown in Fig. 3a.

Fig. 3.

(a) Primal grid (gray) and dual grid (black) and some components on these grids. Stencils of triangles for (b) APVM and (c) LUST, CLUST, and van Leer. (d) Figure 12 reprinted from Thuburn (2008, with permission from Elsevier) showing the smallest-scale structure of the spurious branch Rossby mode. Filled circles indicate positive vorticity; open circles indicate negative vorticity.

Fig. 3.

(a) Primal grid (gray) and dual grid (black) and some components on these grids. Stencils of triangles for (b) APVM and (c) LUST, CLUST, and van Leer. (d) Figure 12 reprinted from Thuburn (2008, with permission from Elsevier) showing the smallest-scale structure of the spurious branch Rossby mode. Filled circles indicate positive vorticity; open circles indicate negative vorticity.

1) APVM

Ringler et al. (2010) described an APVM for controlling grid-scale PV oscillations:

 
formula

where qe is the PV on edge e; qu and qd are the PV in the upwind and downwind dual cells (triangles), respectively; ue is the velocity of edge e; eq is the PV gradient on edge e; and δt is the time step. In this implementation, eq is calculated as a difference between qu and qd divided by the distance between them and then the orthogonal component is reconstructed using the gradients of q normal to the other edges of the dual cells on either side of e using the TRiSK reconstruction weights (Thuburn et al. 2009) on the dual mesh. The reasoning for this scheme is that the PV is interpolated to a point upstream of the edge; an approximate departure point. When the flow is normal to the dual cell edge (parallel to the primal cell edge) APVM is a blend between midpoint differencing and first-order upwind. At other angles, information from lateral dual cells is used, through the eq term. Ringler et al. (2010) used two different methods of calculating the component of eq parallel to the dual edge; this was estimated using the difference between q in the primal cells at either end of the dual edge or was ignored (T. Ringler 2011, personal communication). Using q in the primal cells (i.e. the hexagons) effectively has a stencil of 10 triangles on a hexagonal grid whereas ignoring the component parallel to the dual edge has a stencil of just two dual cells (just qu and qd). The version with 10 triangles is less effective at removing grid-scale oscillations than the version of APVM described above, which has a stencil of 6 triangles (results not shown). The stencil of 6 triangles used for this implementation of APVM is shown in Fig. 3b.

2) Scale-aware APVM

The time step controls the amount of upwinding used in APVM, which is not necessarily proportional to the amount of upwinding needed to suppress grid-scale oscillations. This motivated the scale-aware APVM (Chen et al. 2011) in which the parameter is proportional to the strength of the grid-scale eddies rather than the time step. A derivation resulted in the following formula for the edge vorticity:

 
formula

where α is a constant independent of the scale and state of the flow with an optimal value found to be α = 0.0013, is the kinetic energy averaged over a region around edge e, is the height averaged in the same region, and ℓ is the distance between cell centers. Q. Chen (2011, personal communication) advised that the global averages for and were used for their experiments although more local values would be more relevant and more computationally efficient. Global values are used here to avoid more local variations in the upwinding of qe.

3) LUST

We are not searching for strict boundedness of PV since the consequences of unbounded PV advection are not catastrophic on the solution (unlike negative density or temperature). However, we are searching for a scheme that removes grid-scale PV oscillations more effectively than APVM. The majority of the schemes that are tested are blends between centered, linear interpolation (cℓ) and linear upwind (ℓu). These blends will be denoted linear-upwind stabilized transport (LUST). It is expected that these will be more accurate than APVM as they are a blend between two linear accuracy schemes rather than a linear and first-order scheme. However, neither scheme guarantee boundedness.

Centered, linear interpolation is defined as

 
formula

and linear upwind is a Taylor-series expansion about the upwind dual cell center:

 
formula

where the cell center gradient at the upwind cell, uq, is calculated using Gauss’s theorem with linear interpolation:

 
formula

where qec is the PV interpolated onto edge e′ using centered linear interpolation, de is the length of dual cell edge e′, and ie is the outward-pointing unit normal to e′. LUST with blend b (LUSTb) is thus defined as

 
formula

On a hexagonal grid, LUST has a stencil of four triangles (Fig. 3c). Some comments on the choice of b are made in section 3a(6) below and in section 4, results are presented using LUST15, LUST25, and LUST50 (i.e., 15%, 25%, and 50% linear upwind relative to centered linear, respectively).

4) CLUST

A problem with LUST, in common with monotonic advection schemes, is that there is switching between asymmetric stencils on one side or the other of an edge as the flow changes direction. This is particularly bad if the flow is parallel with the edge; hence, the upwind side of the edge is not well defined. The consequences on accuracy will be demonstrated in section 4c. APVM does not have this problem since the stencil is always symmetric. To overcome this problem with LUST, the blending with centered linear interpolation is dependent on the flow direction, yielding continuous LUST (CLUST):

 
formula

where ue is the velocity normal to edge e and ue is the total velocity at edge e. CLUST also has a stencil of four triangles.

5) Monotonic scheme with a van leer limiter

This standard van Leer scheme is a total variation diminishing (TVD) blend between centered linear and first-order upwind:

 
formula

This uses a stencil of four triangles.

6) Comments on the advection schemes

The computational mode at the grid scale of the spurious Rossby mode branch for the C grid on hexagons has uniform height and a velocity field as shown in Fig. 3d (from Thuburn 2008). When this mode is present, some upwinding is needed so as to use more of qu than qd in the approximation of qe, so that the momentum equation can see the mode and act to remove it. However, for this mode, ueeq is exactly equal to zero, so APVM, scale-aware APVM, and CLUST will not produce any upwinding in the presence of this mode alone, so it will not be damped. The situation is worst for scale-aware APVM since the upwinding factor is |ueeq| ueeq rather than just ueeq as in APVM. This potential problem must be balanced against the need not to include upwinding when the mean flow is close to perpendicular to an edge of the primal grid (near parallel to an edge of the dual grid), as shown in Fig. 3a. If there is no clear upwind dual cell relative to an edge, then midpoint should be used. Otherwise round-off errors in a perfectly aligned wind could lead to a bias toward qu over qd in the estimate of qe. This is a potential problem for LUST and van Leer.

7) Numerical analysis of LUST and choice of parameter b

LUST can be analyzed on a uniform, 1D grid considering interpolation of q from points xj−1, xj, and xj+1 onto xj+1/2½, where xj = x0 + jΔx and qj = q(xj) for constant, positive wind, u. LUST can then be written as

 
formula

The value b = ½ on a uniform 1D grid gives third-order accuracy interpolation (the well-known quadratic upwind or QUICK scheme). However, using LUST to calculate a divergence using Gauss’s theorem requires b = ⅔ for third-order accuracy. In release 2.1.0 of OpenFOAM (OpenCFD Foundation 2011) b = 0.25 is recommended.

To guide the choice of b we can also use von Newmann stability analysis to find the amplification factor and phase speed for different values of b. If we assume that the linear advection equation, has solutions at time step n and for wavenumber k and if we assume that there are no time-stepping errors, then the amplification factor, A, has the following magnitude:

 
formula

where the Courant number, c = uΔtx and the phase speed, ub, is given by

 
formula

where u is the exact phase speed and the true amplification factor should have |A| = 1. Errors in the amplification factor magnitude and phase speed for LUST with various values of b are shown in Figs. 4a,b. Midpoint (b = 0) gives zero amplification magnitude errors (no diffusion) but artificially slows the smallest wavelengths more than any other value of b (dispersion errors), whereas larger values of b (more linear upwinding) give larger amplitude errors (numerical diffusion damps grid-scale waves) but smaller phase speed errors. For values of b greater than ⅔, the phase speed can be faster than u. The balance between diffusion errors (high b) and dispersion errors (low b) is not clear from this analysis. Numerical experiments are needed.

Fig. 4.

(a) Amplification factor and (b) phase speed errors for LUST on a uniform grid in 1D. (c) Advection of a square wave and a cosine hill using LUST.

Fig. 4.

(a) Amplification factor and (b) phase speed errors for LUST on a uniform grid in 1D. (c) Advection of a square wave and a cosine hill using LUST.

The 1D advection of a square wave and a cosine bell (Fig. 4c) show that values of b > 0 suppress grid-scale oscillations, but do not prevent spurious undershoots and overshoots. This simulation advects the initial profile once around the domain using 40 grid points, a Courant number of 0.2, and Heun time stepping (for consistency with the semi-implicit Crank–Nicolson time stepping used with TRiSK).

b. Crank–Nicolson semi-implicit time stepping

The momentum equation (1) and the continuity equation (2) are discretized using second-order, two time level, Crank–Nicolson time stepping:

 
formula
 
formula

where superscripts 0 and 1 represent values at the old and new time steps and superscript ½ is the arithmetic mean between the old and new time step:

 
formula

Here u1 from (14) is substituted into (15) and, using (16) and after some rearrangement we get the Helmholtz equation for h1:

 
formula

In order that this is a linear equation, the first h1 inside each ⋅ has been replaced by a lagged value, h; the value of h1 from the previous iteration (during the first iteration it is simply h0). Two outer iterations are performed each time step in order that lagged values are sufficiently up to date. All values on the right-hand side of (17) are also lagged.

TRiSK energy-conserving operators are used for all the spatial discretizations in (17). The implementation of these operators in OpenFOAM is described in Weller et al. (2012). The matrix resulting from (17) is solved using a conjugate gradient linear equation solver with incomplete Cholesky preconditioning (from OpenFOAM) to a tolerance of 10−2 times the initial residual during the first iteration and to 10−12 during the second iteration. After each iteration, u1 is updated using (14).

4. Results

Results from test cases 2 and 5 of Williamson et al. (1992) (steady-state solid-body rotation and the flow over a mountain) and the Galewsky et al. (2004) barotropically unstable jet are presented on different resolutions of hexagonal icosahedral grid and, for the jet, on a skipped latitude–longitude grid. Grid resolutions and time steps are given in Table 1. The simulations use a variety of advection schemes for PV. The time steps used are generally very short; far shorter than the stability limit. This is because grid-scale oscillations are stronger when the time step is shorter, so using a short time step is more challenging.

a. Williamson et al. (1992): Steady-state, solid-body rotation

Williamson et al.’s (1992) test case 2 is steady state, so errors are calculated with respect to the analytic initial conditions sampled at the cell and edge centers. Convergence with resolution of the ℓ2 and ℓ error norms of height after 5 days are shown for the hexagonal grid in the Figs. 5a,c (The ℓ1 error norms behave very similarly to the ℓ2 error norms and are not shown). For sufficient resolution, ℓ2(h) converges with second-order accuracy, as expected using TRiSK on the near-centroidal Heikes and Randall (1995) hexagonal grid. Here ℓ(h) converges with first-order accuracy due to errors at the grid distortions along the base icosahedron edges where the grid is less centroidal than elsewhere (the computational points are not at the cell centers and the edge crossover points are not at the edge midpoints). The height error norms do not clearly show the differences between the PV advection schemes since the PV only influences the height directly through the nonlinear terms in the momentum equation and the nonlinear terms are small for this test case. Therefore, we also present the ℓ2 and ℓ error norms of the PV in Figs. 5b,d. For sufficient resolution, these all converge with first-order accuracy since vorticity is a spatial derivative of velocity and height and velocity are treated with close to second-order accuracy. The PV errors are much higher using midpoint and scale-aware APVM since these schemes do not dampen the computational mode.

Fig. 5.

(a)–(d) Convergence with resolution of the ℓ2 and ℓ error norms of height and PV for the Williamson et al. (1992) test case 2 on hexagonal grids with 642, 2562, 10 242, and 40 962 cells and for various advection schemes for PV. (e),(f) Diagnostics for the 40 962 cell grid using a time step of 450 s.

Fig. 5.

(a)–(d) Convergence with resolution of the ℓ2 and ℓ error norms of height and PV for the Williamson et al. (1992) test case 2 on hexagonal grids with 642, 2562, 10 242, and 40 962 cells and for various advection schemes for PV. (e),(f) Diagnostics for the 40 962 cell grid using a time step of 450 s.

TRiSK conserves mass and PV to within machine precision and energy dependent on the time step and time-stepping scheme. We have verified that these conservation properties are independent of the advection scheme for PV (not shown). However, conservation of potential enstrophy is dependent on the PV advection and it is desirable that this should decay slowly rather than increase. The change in normalized potential enstrophy (as defined in Williamson et al. 1992) is shown in Fig. 5e for the solid-body rotation test case using the hexagonal grid of 40 962 cells and a time step of 450 s. The growing computational modes lead to increasing potential enstrophy using midpoint and scale-aware APVM. For this test case, scale-aware APVM is behaving exactly as midpoint since uq is equal to zero on the large scales and for the growing mode and the damping in scale-aware APVM is heavily dependent on uq. The other upwinded schemes have the desired effect of slowly removing potential enstrophy, but does this suppress the computational mode?

The shortcoming of potential enstrophy as a diagnostic of the effectiveness of a scheme at suppressing computational modes is that it is dependent on the large and small scales, so enstrophy can be destroyed overall while grid-scale oscillations grow. To diagnose solely the presence of grid-scale checkerboard patterns, we propose the diagnostic:

 
formula

where i is the index over all dual cells and Ii is an indicator function for the presence of a checkerboard pattern:

 
formula

where the js are the neighbor dual cells to i and nj is the number of neighbors of cell i. This implies that dual cell i only makes a contribution if has one sign while all the neighbors have the opposite sign, which only occurs at grid-scale oscillations.

Midpoint and scale-aware APVM do not suppress the computational modes and so G grows using these schemes (Fig. 5f), but not with the others.

b. Williamson et al. (1992): Flow over a midlatitude mountain

If Williamson et al.’s (1992) test case 5 is run for 50 days rather than just 15 days, the grid-scale oscillations due to the computational modes of the hexagonal grid begin to overwhelm the large-scale solution unless they are damped by an upwinded advection scheme for PV (Fig. 1). This simulation originally motivated the search for a better advection scheme than APVM since grid-scale PV oscillations are present (e.g., northwest of the mountain using APVM). All of the LUST schemes control these oscillations with LUST50 being the smoothest LUST and van Leer being smoother still. It is also interesting to compare these with results from the National Centre for Atmospheric Science Spectral Transform Model (STSWM; Jakob-Chien et al. 1995), updated by Ripodas (2012) both with and without diffusion. At T213 (320 × 640 grid points), the solution is very noisy without diffusion and mostly smooth with a diffusion coefficient of 8 × 1012 m2 s−1, but neither of the solutions capture the detail captured by TRiSK, despite TRiSK being a low-order scheme at lower resolution (~120 km as opposed to ~60 km). It is therefore not considered useful to compare TRiSK with STSWM more quantitatively for this test case. The scale-aware APVM removes some of the PV oscillations for this test case but not as many as APVM.

The proximity of the anticyclone (negative vorticity) north of the mountain to the Pole varies between the simulations. The schemes that can maintain tighter vorticity gradients and PV filaments without introducing PV oscillations (LUST25, LUST15, and CLUST25) have the anticyclone closer to the Pole. The others, whether they are very smooth (LUST50, van Leer, and STSWM with diffusion) or excessively oscillatory (STSWM without diffusion, midpoint, and both APVMs) have the anticyclone farther south. It is therefore possible that the enstophy dissipation caused by smooth schemes or the nonlinear interactions caused by grid-scale oscillations erroneously move the anticyclone south.

The normalized potential enstrophy and the grid-scale oscillation diagnostic G in (18) are shown in Fig. 6. For this long test case, the scale-aware APVM only generates half the amount of potential enstrophy as midpoint, but all of the other schemes slowly remove potential enstrophy with van Leer removing the most and APVM removing the least. The overall removal of potential enstrophy by APVM hides the fact that grid-scale oscillations are maintained by APVM (Fig. 1i). However the extent of grid-scale oscillations using both APVMs is confirmed by G in Fig. 6b. The value of G is exactly as expected using the other schemes: LUST15 > CLUST25 > LUST25 > LUST50. Using van Leer almost no grid-scale oscillations are generated, consistent with this being a monotonic advection scheme.

Fig. 6.

Flow diagnostics for the Williamson et al. (1992) flow over a midlatitude mountain using various interpolation schemes for PV. Hexagonal grid of 40 962 cells. A time step of 100 s.

Fig. 6.

Flow diagnostics for the Williamson et al. (1992) flow over a midlatitude mountain using various interpolation schemes for PV. Hexagonal grid of 40 962 cells. A time step of 100 s.

The flow over a mountain test case is not run for 50 days on the skipped latitude–longitude grid because errors from the grid reductions would be too damaging to the global solution after such a long simulation.

c. Galewsky et al. (2004): Barotropically unstable jet

Before presenting results of the barotropically unstable jet with an initial perturbation that initiates the instability, the jet is simulated without the initial perturbation in order to assess the model’s ability to maintain the geostrophically balanced but unstable jet. Maintaining steady but unstable flow with a grid that is not zonally symmetric is known to be very difficult (e.g., Ii and Xiao 2010). The ℓ1(h) and ℓ1(q) error norms for this jet for 5 days are presented in Fig. 7 using the hexagonal grid of 163 782 cells, using various advection schemes for PV and a time step of 100 s and are compared with results of the same test case by Ii and Xiao (2010). Since this is a steady-state case the errors are calculated as differences from the initial conditions. The TRiSK results have 655 122 d.o.f., whereas Ii and Xiao (2010) have 60 × N2 + 2 d.o.f. for each variable (1 for height and 2 for velocity) making 648 002 d.o.f. altogether for N = 60. The early times of slow error growth for both norms are due to the ability of TRiSK to maintain geostrophic balance. The errors initially grow fast while they are moving from the initial conditions to close to discrete balance but then grow more slowly for a few days while they are close to discrete balance. This contrasts with the results of Ii and Xiao (2010) using a higher-order model that does not have stationary geostrophic modes. Their results grow close to exponentially throughout. Once the near-steady-state breaks down (2–3 days) the TRiSK errors grow more rapidly than Ii and Xiao (2010) due to the low-order accuracy of TRiSK.

Fig. 7.

The ℓ1 error norms of (a) height and (c) potential vorticity for the steady-state (but unstable) jet of Galewsky et al. (2004) without the initial perturbation on a hexagonal grid of 163 782 cells using various advection schemes for PV and a time step of 100 s. (b) Reprinted from Ii and Xiao (2010) with permission from Elsevier. The TRiSK results have d.o.f. equivalent to N = 60 for Ii and Xiao (2010).

Fig. 7.

The ℓ1 error norms of (a) height and (c) potential vorticity for the steady-state (but unstable) jet of Galewsky et al. (2004) without the initial perturbation on a hexagonal grid of 163 782 cells using various advection schemes for PV and a time step of 100 s. (b) Reprinted from Ii and Xiao (2010) with permission from Elsevier. The TRiSK results have d.o.f. equivalent to N = 60 for Ii and Xiao (2010).

Midpoint interpolation and scale-aware APVM do not control the computational modes of TRiSK and so the ℓ1(q) errors are higher using these schemes (Fig. 7c). However, the computational modes do not have a direct or immediate impact on the height errors (Fig. 7a). In fact, most of the schemes that control the computational mode in PV initially increase the height errors, apart from CLUST, which has the lowest errors for both norms.

The relative vorticity of the inviscid barotropically unstable jet with the initial perturbation using all the interpolation schemes for PV on the hexagonal-icosahedral grid of 163 782 cells at day 6 is shown in Fig. 8 and is compared with results of STSWM at T213 (320 × 640) without diffusion. Both these resolutions correspond to about 60 km. Without diffusion, STSWM clearly generates considerable spectral ringing but the large scales are very similar for both models. A time step of 100 s is used for both models, which resolves gravity waves (Table 1) and enhances the generation of spurious grid-scale oscillations when using TRiSK.

Fig. 8.

Relative vorticity after 6 days for the Galewsky et al. (2004) barotropically unstable jet on a hexagonal grid of 163 782 cells, using various advection schemes for PV and time step of 100 s in comparison to STSWM at T213 (320 × 640). Plots from 10° to 80°N with a tick at 45°N. Zoomed region from 57°–68°N, 145°–122°W.

Fig. 8.

Relative vorticity after 6 days for the Galewsky et al. (2004) barotropically unstable jet on a hexagonal grid of 163 782 cells, using various advection schemes for PV and time step of 100 s in comparison to STSWM at T213 (320 × 640). Plots from 10° to 80°N with a tick at 45°N. Zoomed region from 57°–68°N, 145°–122°W.

All of the schemes remove more grid-scale oscillations than the scale-aware APVM. On close inspection of the APVM solution, particularly in the zoomed region, spurious grid-scale oscillations are present. These are removed when using any version of LUST or van Leer, with van Leer and LUST50 being smoother than LUST25, LUST15, and CLUST25.

Using van Leer, the wave at 160°E has grown more than all of the other solutions and additional comparisons with the higher-resolution results presented in Galewsky et al. (2004) reveal that this additional growth is erroneous. The low-order accuracy when the limiter is active in van Leer can spuriously trigger barotropic instability since truncation errors can upset the balance. Initially, the jet has meridional velocity, υ = 0. A small error in υ (say υ < 0) combined with upwind advection will lead to positive PV anomaly being advected from north of the jet center toward the jet center. Since ∂υ/∂t is controlled by −qu, the positive PV will decrease υ farther. Thus, the small error will grow.

All of the TRiSK solutions suffer from spurious stripes in the vorticity field that are larger than the dual grid size, apart from the very smooth van Leer solutions. These are thought to be due to the energy-conserving averaging used in TRiSK, which modifies the implied advection scheme for PV. However, the origins, consequences, and solutions of this are the subject of ongoing research.

The barotropically unstable jet is also simulated on a skipped latitude–longitude grid of 59 954 cells (192 × 384) using a time step of 100 s and compared with solutions on a skipped latitude–longitude grid of 241 970 cells (384 × 768) in Fig. 9. The jet is initially centered at 45°N and the grid reductions are at 60°N so the waves in the jet should not reach the grid reductions until day 5. Therefore, unless there is spurious early barotropic instability, this grid largely acts like a full latitude–longitude grid for this test case run for 6 days.

Fig. 9.

Relative vorticity after 6 days for the Galewsky et al. (2004) barotropically unstable jet on a reduced latitude–longitude grid of 59 954 cells, using various advection schemes for PV and a time step of 100 s in comparison to a reduced latitude–longitude grid of 241 970 cells. Plots from 10° to 80°N with a tick at 45°N. Zoomed region from 57°–68°N, 145°–122°W.

Fig. 9.

Relative vorticity after 6 days for the Galewsky et al. (2004) barotropically unstable jet on a reduced latitude–longitude grid of 59 954 cells, using various advection schemes for PV and a time step of 100 s in comparison to a reduced latitude–longitude grid of 241 970 cells. Plots from 10° to 80°N with a tick at 45°N. Zoomed region from 57°–68°N, 145°–122°W.

LUST15, LUST25, and LUST50 all have spurious release of barotropic instability and erroneous waves at 160°E for this test case, although jets aligned with the grid should be easy to maintain at coarse resolution. The spurious waves have grown so much using van Leer that they have reached the grid reductions sooner and thus spurious solutions have had longer to grow. LUSTb and van Leer are all based on upwind stencils that switch suddenly between upwind sides when the flow is along a dual cell edge (upwind-biased schemes are derived assuming the flow is normal to the edge in question). Again, considering the initial jet, which should have υ = 0, a round-off error in υ leading to υ < 0 will lead to a stencil for interpolating PV, which is north of an east–west edge. This will lead to a positive bias in the estimation of PV at the edge that will act to further decrease υ, thus triggering the barotropic instability. The poor behavior of LUST on this grid and test case in comparison to APVM was the original motivation for CLUST, which does not use upwinding when the flow is aligned with an edge. The CLUST results do not have any spurious waves for this test case.

The normalized potential enstrophy and the grid-scale oscillation diagnostic G in (18) for the unstable jet are shown in Fig. 10. On the hexagonal grid, midpoint and scale-aware APVM generate enstrophy (top left) due to the growth of the computational mode (the grid-scale oscillations of PV) whereas all the other schemes lead to losses of enstrophy, most strongly for APVM. While APVM leads to the largest loss of enstrophy, it does not do well at minimizing G (bottom left of Fig. 10), implying that it is destroying enstrophy at large scales without preventing oscillations at the grid scale. Conversely, LUST and CLUST remove less enstrophy while producing fewer grid-scale oscillations.

Fig. 10.

Flow diagnostics for the Galewsky et al. (2004) barotropically unstable jet for various advection schemes for PV and a time step of 100 s.

Fig. 10.

Flow diagnostics for the Galewsky et al. (2004) barotropically unstable jet for various advection schemes for PV and a time step of 100 s.

On the skipped grid, midpoint and scale-aware APVM do not lead to such big increases in enstrophy (top right of Fig. 10) since TRiSK on this grid does not suffer from as many computational modes. Therefore, the reductions in enstrophy due to the other advection schemes are not desirable.

5. Summary and conclusions

It has been shown that grid-scale oscillations associated with the computational modes of the hexagonal C grid are easy to control using an upwind-biased advection scheme for PV when using the TRiSK algorithm for solving the vector invariant shallow-water equations. APVM (Ringler et al. 2010) and the scale-aware version (Chen et al. 2011) do not always control the computational mode but a number of proposed split space/time schemes perform better. Linear upwind stabilized transport (LUST) is a fixed blend between linear upwind and centered linear. Values of the blend of 0.5, 0.25, and 0.15 linear upwind have been found to be sufficient to control the mode without destroying the balance, conservation, or accuracy of TRiSK. A blend of 0.25 is recommended. A monotonic, van Leer blend between first-order upwind and centered linear controls the mode but is more dissipative of potential enstrophy.

However, LUST does not work well for all test cases. Like van Leer and other upwind-biased schemes, there is switching between stencils as the flow goes from upwind on one side to upwind on the other side of an edge. This can lead to, for example, spurious release of barotropic instability. We therefore propose CLUST, a linear combination of centered-linear and linear-upwind interpolation, with the blend being dependent on the flow direction so that when flow is parallel to an edge, CLUST becomes centered linear and when it is normal to an edge it is LUST. The flow direction dependence removes the abrupt switching between upwind-biased stencils on one side or the other of an edge, which can introduce errors, and CLUST25 works well on all the test cases.

LUST and CLUST are dependent on a coefficient and exhaustive sensitivity to this coefficient has not been undertaken. However, it is expected that LUST and CLUST are less sensitive to this coefficient than schemes that blend centered linear with first-order upwind since LUST and CLUST blend two linear schemes. The disadvantage of this is that monotonicity cannot be guaranteed. The motivation for LUST and CLUST was to control computational modes and grid-scale oscillations, not to guarantee monotonicity and they do this well.

A diagnostic of grid-scale oscillations is proposed that is more discriminating than using potential enstrophy alone. Some schemes, in particular APVM, can lead to reductions in global potential enstrophy while grid-scale oscillations grow, as measured by this diagnostic.

A major objection raised by Staniforth and Thuburn (2012) to the hexagonal icosahedral grid of the sphere was that it supports a branch of spurious Rossby modes. Here we show that these can be cheaply and effectively controlled without upsetting the accuracy, balance, or conservation properties of TRiSK. The results of using TRiSK on a hexagonal grid often appear to be better than those of a spectral transform model (STSWM; Jakob-Chien et al. 1995; Ripodas 2012) and well in comparison to a higher-order model in a hexagonal grid (Ii and Xiao 2010). The hexagonal icosahedron with TRiSK appears to be a good choice for a quasi-uniform global atmosphere model.

Acknowledgments

H.W. acknowledges support from NCAS Climate and NERC Grants NE/H002774/1, NE/H001166/1, and NE/I022086/1. I am also grateful to John Thuburn, Colin Cotter, and Todd Ringler for useful discussions and for constructive comments from two anonymous reviewers.

REFERENCES

REFERENCES
Chen
,
Q.
,
M.
Gunzburger
, and
T.
Ringler
,
2011
:
A scale-invariant formulation of the anticipated potential vorticity method
.
Mon. Wea. Rev.
,
139
,
2614
2629
.
Danilov
,
S.
,
2010
:
On utility of triangular C-grid type discretization for numerical modeling of large-scale ocean flows
.
Ocean Dyn.
,
60
(
6
),
1361
1369
.
Galewsky
,
J.
,
R.
Scott
, and
L.
Polvani
,
2004
:
An initial-value problem for testing numerical models of the global shallow-water equations
.
Tellus
,
56A
,
429
440
.
Gassmann
,
A.
,
2011
:
Inspection of hexagonal and triangular C-grid discretizations of the shallow water equations
.
J. Comput. Phys.
,
230
,
2706
2721
.
Heikes
,
R.
, and
D.
Randall
,
1995
:
Numerical integration of the shallow-water equations on a twisted icosahedral grid. Part II: A detailed description of the grid and an analysis of numerical accuracy
.
Mon. Wea. Rev.
,
123
,
1881
1887
.
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
.
Jakob-Chien
,
R.
,
J.
Hack
, and
D.
Williamson
,
1995
:
Spectral transform solutions to the shallow water test set
.
J. Comput. Phys.
,
119
,
164
187
.
Le Roux
,
D.
,
A.
Sèneb
,
V.
Rostanda
, and
H.
Hanert
,
2005
:
On some spurious mode issues in shallow-water models using a linear algebra approach
.
Ocean Modell.
,
10
(
1–2
),
83
94
.
OpenCFD Foundation
, cited
2011
:
The opensource CFD toolbox. [Available online at http://www.openfoam.org.]
Putman
,
W.
, and
S.-J.
Lin
,
2007
:
Finite-volume transport on various cubed-sphere grids
.
J. Comput. Phys.
,
227
,
55
78
.
Rančić
,
M.
,
R.
Purser
, and
F.
Mesinger
,
1996
:
A global shallow-water model using an expanded spherical cube: Gnomonic versus conformal coordinates
.
Quart. J. Roy. Meteor. Soc.
,
122
(
532
),
959
982
.
Randall
,
D.
,
1994
:
Geostrophic adjustment and the finite-difference shallow-water equations
.
Mon. Wea. Rev.
,
122
,
1371
1377
.
Ringler
,
T.
,
J.
Thuburn
,
W.
Skamarock
, and
J.
Klemp
,
2010
:
A unified approach to energy conservation and potential vorticity dynamics for arbitrarily-structured C-grids
.
J. Comput. Phys.
,
229
,
3065
3090
.
Ripodas
,
P.
, cited
2012
:
Spectral transform shallow water model. [Available online at http://icon.enes.org/swm/stswm/.]
Staniforth
,
A.
, and
J.
Thuburn
,
2012
:
Horizontal grids for global weather and climate prediction models: A review
.
Quart. J. Roy. Meteor. Soc.
,
138
(
662
),
1
26
.
Thuburn
,
J.
,
2008
:
Numerical wave propagation on the hexagonal C-grid
.
J. Comput. Phys.
,
227
,
5836
5858
.
Thuburn
,
J.
,
T.
Ringler
,
W.
Skamarock
, and
J.
Klemp
,
2009
:
Numerical representation of geostrophic modes on arbitrarily structured C-grids
.
J. Comput. Phys.
,
228
(
22
),
8321
8335
.
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
.
Weller
,
H.
,
J.
Thuburn
, and
C.
Cotter
,
2012
:
Computational modes and grid imprinting on five quasi-uniform spherical C grids
.
Mon. Wea. Rev.
,
140
,
2734
2755
.
Williamson
,
D.
,
J.
Drake
,
J.
Hack
,
R.
Jakob
, and
P.
Swarztrauber
,
1992
:
A standard test set for numerical approximations to the shallow water equations in spherical geometry
.
J. Comput. Phys.
,
102
,
211
224
.