## Abstract

The logically rectangular finite volume grids for two-dimensional partial differential equations on a sphere and for three-dimensional problems in a spherical shell introduced recently have nearly uniform cell size, avoiding severe Courant number restrictions. We present recent results with adaptive mesh refinement using the GeoClaw software and demonstrate well-balanced methods that exactly maintain equilibrium solutions, such as shallow water equations for an ocean at rest over arbitrary bathymetry.

## 1. Introduction

Recently, Calhoun *et al.* (2008*b*) introduced a class of logically rectangular grids (quadrilateral in two dimensions or hexahedral in three dimensions) that can be used for solving partial differential equations on a sphere, illustrated in figure 1*a*. These grids have the virtue that the cell sizes are nearly uniform (unlike latitude–longitude grids, for example), while maintaining the simple (*i*,*j*,*k*) indexing of logically rectangular grids. On the sphere the ratio of largest to smallest grid cell area is less than 1.7. A potential disadvantage is that the four corners of the computational rectangle are mapped to quadrilaterals that are nearly triangular; the grid lines are far from orthogonal in these regions. We have demonstrated, however, that appropriate finite volume methods can be effectively used on these grids that achieve second-order accuracy on smooth solutions and that also capture shock waves and other discontinuities sharply (Calhoun *et al.*2008*a*,*b*; Calhoun & Helzel in press). In particular, we have presented results for the shallow water equations on the sphere for both shock waves and the Rossby–Haurwitz test problem of Williamson *et al.* (1992) on a single grid as well as some preliminary results with adaptive mesh refinement (AMR).

In this paper we present additional results for shallow water on the sphere and introduce a finite volume method for flow over arbitrary topography or bathymetry on the sphere. We show that the method is ‘well balanced’ for the case of an ocean at rest, meaning that it preserves the stationary steady state in this case and also accurately computes small amplitude wave propagation on top of this equilibrium solution.

The test problems are implemented in GeoClaw, a subset of the Clawpack software that has been developed for geophysical flows in work with David George (George & LeVeque 2006, 2008; LeVeque & George 2007; George 2008), initially in the context of tsunami modelling in the mid-latitudes on a latitude–longitude grid. The complete source code for all the examples presented here will be made available on the Web (www.clawpack.org/links/sphere).

The AMR code in GeoClaw is based on AMRClaw, originally developed by Berger & LeVeque (1998). Refinement is done on logically rectangular patches in the computational domain, an approach that has been widely used, including for three-dimensional flow on the sphere in the recent work of Hubbard & Nikiforakis (2003), for example. The grid mappings we use here have the advantage that the computational domain is a single rectangle, greatly simplifying the application of adaptive refinement.

## 2. The grid mappings

Figure 2 illustrates the mapping from the computational domain (figure 2*a*) to the quadrilateral grid we use on the sphere (figure 2*d*).

Details of this mapping are presented in Calhoun *et al.* (2008*b*) and a Python implementation is displayed in figure 3. Here (xc,yc) is a point in the computational domain [−3,1]×[−1,1] and (xp,yp,zp) is the image on the sphere. The square [−3,−1]×[−1,1] is mapped to the lower hemisphere, while [−1,1]×[−1,1] is mapped to the upper hemisphere.

Note one modification from the mapping presented in Calhoun *et al.* (2008*b*): the expression for *D* used here involves rather than *d*(2−*d*) as in the original mapping. This leads to slightly more uniform cells and is more natural since equally spaced points along a diagonal are then mapped to points that are equally spaced in latitude in the image. In practice it makes little difference in computed results. As written here, this mapping produces the grid shown in figure 1*a*.

We also mention a modified version of this mapping, which can be obtained by un-commenting the line `R = Rsphere * sin(pi*d/2)`. With this choice the grid shown in figure 1*b* is obtained. This grid has the possible advantage that one set of grid lines are lines of constant latitude. (The latitude lines are concentric squares in the computational domain.) However, this grid is not as smooth near the poles as the grid shown in figure 1*a*, and in practice has not been found to be advantageous, even for the test problems of §7 where the true solution is axi-symmetric (results not shown here). The grid may prove useful for specialized applications however, such as purely zonal advection.

The quadrilateral grid we use on the sphere is easily extended to a three-dimensional hexahedral grid in the normal direction. The computational domain is again purely Cartesian [−3,1]×[−1,1]×[*R*_{1},*R*_{2}] for a shell with inner radius *R*_{1} and outer radius *R*_{2}. We are currently working on the development of a well-balanced AMR method for solving atmospheric flow problems on this grid, for situations in which the flow is a small perturbation of an equilibrium with a hydrostatic pressure gradient in the normal direction. This work will be presented elsewhere.

## 3. Shallow water equations on the sphere

Since there is no smooth coordinate system for our grids on the sphere, we keep track of the fluid velocity in terms of a full three-dimensional vector in each finite volume grid cell. These define a velocity vector that should be tangent to the sphere. Similar approaches are used, for example, in Swarztrauber *et al.* (1997) and Giraldo *et al.* (2002).

We use the geopotential height and solve the shallow water equations on the sphere in the form
3.1Here *q*=[*ϕ*,*ϕu*,*ϕv*,*ϕw*]^{T} and is the surface geopotential.
3.2is the momentum flux. The symbol denotes the projection matrix that projects a vector at to the tangent plane,
3.3The divergence of the flux and the last term on the right in equation (3.1), which appears for flow over variable topography specified by a surface geopotential , are projected to the tangent plane. The Coriolis term involving the rotation rate Ω has already been projected to the tangent plane. Giraldo *et al.* (2002) provide justification that these equations are equivalent to standard formulations of shallow water on the sphere.

Note that for shallow incompressible water of depth *h* on the sphere the geopotential is *ϕ*=*gh*, where *g* varies with latitude to incorporate the centrifugal force on a rotating sphere.

## 4. Topography and well-balanced methods

In applications such as tsunami modelling, it is very important that the numerical method maintains the equilibrium solution consisting of an ocean at rest, for which the depth *h* and the bathymetry *B* both vary greatly, but the surface elevation *η* ≡ *h*+*B* is constant and the velocities are identically 0. Numerical methods often do not maintain this steady state unless care is taken, since the divergence of the flux and the source term are both large but cancel out. In particular, fractional step methods in which one alternates between solving the homogeneous conservation laws and the source term equations do not work well at all (e.g. LeVeque 1998). In terms of the geopotential *ϕ*, the ocean at rest equilibrium solution has *ϕ*+*gB* constant (where again *g* varies with latitude if the sphere is rotating).

The high-resolution wave-propagation method that we present in §6 exactly preserves such an equilibrium solution and also accurately captures small amplitude waves on top of such a steady state. Much research has recently been devoted to developing such ‘well-balanced’ methods, particularly for shallow water equations (e.g. Gosse 2000, 2001; Audusse *et al.* 2004; Bouchut 2004; Castro *et al.* 2008).

## 5. Boundary conditions

For finite volume methods of the type we use, boundary conditions at the edges of the computational rectangle are generally implemented using the standard procedure of introducing two rows of ‘ghost cells’ at each edge of the computational domain. Values in these cells are set at the start of each time step by copying data from elsewhere in the domain, since each edge of the domain is adjacent to another edge after performing the mapping illustrated in figure 2. After setting the ghost cells in this manner, the wave-propagation algorithm is applied at all points in the computational domain. We fill two rows of ghost cells because this method has a stencil that extends two grid cells out in each direction (the adjacent cell value is used in solving the Riemann problem at the cell interface and the value in the next proximate grid cell is required in applying limiters to the resulting waves for the high-resolution method).

Let (*x*_{c},*y*_{c}) represent a point in the computational domain. After the mapping, the boundaries *x*_{c}=−3 and *x*_{c}=+1 meet up along the equator and these boundaries simply have standard periodic boundary conditions:
5.1for −1≤*y*_{c}≤1. At the boundary *ε*=0, but to fill the ghost cells these expressions are used for *ε*>0. These are easily implemented on a 2*m*×*m* grid by setting
5.2for *j*=1,2, … ,*m* at the start of each time step.

Along *y*_{c}=±1 the boundary conditions are only slightly more complicated. At each *y*_{c}=±1 boundary, the segment −3≤*x*_{c}≤−1 matches up with the segment −1≤*x*_{c}≤1 with the orientation reversed, as seen in figures 2 and 4, and hence
5.3for −3≤*x*_{c}≤1. This is implemented by copying data according to
5.4for *i*=1,2, … ,2*m*.

These boundary conditions were implemented for single-grid calculations in the example presented in Calhoun *et al.* (2008*b*) and Fortran code is available at the webpage for that paper. These boundary conditions were recently implemented in the AMR version of GeoClaw as one possible choice of boundary conditions, facilitating the use of AMR for flow problems on the sphere. With AMR, the non-standard boundary conditions (5.3) require some work to implement, since a grid that is adjacent to the top or bottom boundary will generally have to fill ghost cell values from a different grid at the same level or by interpolating from coarser grids, and logic must be incorporated to search for the relevant data at the appropriate interior point as defined by equation (5.3). Similar issues arise for standard periodic boundary conditions, which were implemented in Berger’s original AMR code that formed the basis for the AMR routines in Clawpack.

Some preliminary results using AMR on the sphere were presented in Calhoun *et al.* (2008*a*), but with a different approach to specifying boundary conditions. We mention this alternative since it may be useful to those wishing to apply other AMR software to the sphere using our grid mapping. The computational domain can be doubled in size to a full square [−3,1]×[−3,1] by defining
5.5for −3≤*x*_{c}≤1. This simply reflects all the data through the point (−1,−1). The mapping can be extended in an obvious way so that the new computational rectangle [−3,1]×[−3,−1] is mapped to a second copy of the sphere and each copy will contain the same value of *q* at each point on the sphere. The advantage of introducing this second copy is that the boundary conditions at *y*_{c}=−3 and *y*_{c}=1 are now the standard periodic boundary conditions in this direction,
5.6for −3≤*x*_{c}≤1. Hence, any AMR code that can handle periodic boundary conditions on a rectangular domain can be applied on the sphere, at double the computer time but without any need to reprogram the AMR boundary conditions.

## 6. The finite volume method

The cell average in cell (*i*,*j*) at time *t*_{n} is updated by waves that propagate into the cell from the neighbouring cell interfaces, which are determined by solving one-dimensional Riemann problems at each interface. The approach we use is the ‘*f*-wave formulation’ of the wave-propagation algorithm introduced in Bale *et al.* (2002) and also discussed in ch. 17 of LeVeque (2002). This means that the waves are determined by computing the jump in the flux plus a contribution from the source terms, and then decomposing this vector into eigenvectors of a linearized Jacobian matrix defined at each cell interface, as described in more detail below. These waves are then also used to define correction terms that yield a method that is second-order accurate on smooth solutions. Limiters are applied to these waves (by comparing with waves from adjacent Riemann problems) in order to avoid non-physical oscillations and produce high-resolution results. ‘Transverse Riemann problems’ must also be solved to determine the portion of waves that move transversely into corner-adjacent cells and this appears to be an important aspect of these algorithms in successfully solving multi-dimensional problems on the highly skewed grids of figure 1. Complete details of this approach can be found in LeVeque (2002).

Here we concentrate on describing the unique aspects of the algorithm for flow over topography on the sphere, which is primarily concerned with the algorithm for solving the one-dimensional Riemann problem at a cell edge to obtain a suitable set of *f*-waves.

Consider, for example, the interface (*i*−1/2,*j*) between cells (*i*−1,*j*) and (*i*,*j*). Let be the momentum components of the cell averages on the (*i*,*j*) grid cell at time *t*_{n}, which we suppose to lie in the tangent plane at the cell centre. We will solve a one-dimensional Riemann problem between states
6.1where and . To obtain the velocities we first project the momentum components of the data in the two neighbouring cells, and , onto a common direction. We choose the direction that is normal to the cell edge and that lies in the tangent plane at the midpoint of the cell edge. This defines the normal momentum components of the Riemann problem data, *ϕ*_{l}*u*_{l} and *ϕ*_{r}*u*_{r}. We also compute tangential momentum components *ϕ*_{l}*v*_{l} and *ϕ*_{r}*v*_{r} by projecting the cell momenta onto the direction determined by the cell edge.

From these Riemann data we can define an approximate Jacobian matrix, using for example the Roe-averaged Jacobian . The eigenvalues of this matrix define the wave speeds and the eigenvectors are used to determine the *f*-waves. This is done by decomposing
6.2as a linear combination of the eigenvectors of ,
6.3The vectors *𝒵*^{p} are the *f*-waves used in the wave-propagation algorithm.

The quantity Δ in equation (6.2) is the flux difference modified by the source term from the varying bathymetry. For the special case of an ocean at rest, all terms in equation (6.2) involving velocities drop out and all that remains are the following terms in the second element of Δ: the pressure gradient and the source term . For the ocean at rest, the sea surface is constant and so it can be easily verified that these two terms cancel. Hence, Δ is the zero vector in this case, with the result that the waves *𝒵*^{p} all have zero strength. For an ocean near rest, it is only the deviation from steady state that is split into waves and propagated with our wave-propagation algorithm.

When AMR is used, some care must be exercised in interpolating values between coarse grids and finer level grids, as required for example when initializing a new fine grid in some region or interpolating ghost cell values from coarser levels. One crucial aspect is to do interpolations based on and then reconstruct *ϕ*, rather than interpolating *ϕ* itself. This is discussed in more detail in the context of tsunami modelling in George (2008).

The area of a quadrilateral grid cell on the sphere can exactly be computed as the sum of the area of two spherical triangles. Note that with our grid mapping, refined grid cells do not exactly match up with the underlying coarse grid cell in physical space. Therefore, we calculate the area of every grid cell as the sum of the areas of grid cells on the finest possible level that might be obtained by refining the grid cell. In an analogous way we also calculate the surface geopotential in each grid cell.

## 7. Numerical results

As a test case for the well-balanced method with adaptive refinement as described in this paper, we have chosen a problem for which the solution varies only with latitude so that a one-dimensional axisymmetric code can be used to generate a fine grid reference solution. A test of this nature was also presented in Calhoun *et al.* (2008*b*), but only on a single grid with no topography.

The one-dimensional reference solution (plotted as a solid line in figure 6) is obtained by approximating the system
7.1Here *ϕ*=*ϕ*(*θ*,*t*) is the geopotential height, is the surface geopotential. *U*=*U*(*θ*,*t*) and *V* =*V* (*θ*,*t*) are velocities in the latitudinal and the longitudinal direction, respectively. Initial value problems of system (7.1) are approximated for latitude *θ* ∈ [−*π*/2,*π*/2]. Reflecting boundary conditions are used on both sides of this interval.

For our test problem we specify axisymmetric topography by choosing the surface potential to be the following function of latitude *θ*:
7.2where *θ* varies from −*π*/2 at the south pole to *π*/2 at the north pole. Sea level is taken to be the surface of geopotential 0, so the surface potential (7.2) corresponds to an ocean of depth 40 000/*g*≈4000 m with an axisymmetric ridge about the sphere at a latitude of −*π*/6 or 30° south. At the peak of the ridge the ocean depth is reduced by half. The ocean at rest corresponds to so that .

As initial conditions we perturb this to
7.3corresponding to an axisymmetric elevation above sea level in a ring of latitude *π*/6. Note the decay factor −1000 in the Gaussian, giving a width of roughly 200 km on the Earth. This is a reasonable wave length for a tsunami wave.

We initialize the zonal velocity *V* to *V* (*θ*,*t*=0)=0 and the meridional momentum *ϕU* to
7.4Here is the approximate wave speed and so the vector (*ϕ*,*ϕU*) is approximately a 1-wave eigenvector, which results in most of the initial perturbation moving southwards, with a much weaker wave moving north. In the two-dimensional code on the sphere we must initialize the three-dimensional momentum components of the solution vector (*ϕu*,*ϕv*,*ϕw*) to agree with this in each grid cell by
7.5

We take *A*=10 as the amplitude of the perturbation in the initial data. This gives *A*/*g*≈1 m, which is realistic for tsunami modelling and tests the ability of our approach to handle a wave that has very small amplitude relative to the variations in topography (about 0.05% in this case).

We include Coriolis effects and take Ω=7.292×10^{−5}, the rotation rate of the Earth. For the radius of the Earth, we use *R*=6.371×10^{6}. In both the one- and two-dimensional code we compute latitude relative to a perfect sphere, though this could be replaced by an appropriate geoid.

Figure 5 shows how the AMR refinement might look at time *t*=0 in the case of very coarse grids.

Results of this test are shown in figure 6. Figure 6*a* shows scatter plots of the solution when a single 800×400 grid is used on the sphere, at three different times. Figure 6*b* shows results where refinement is done near the main wave with the following refinement factors:

level 1: 50×25 grid,

level 2: refined by 4 (effective resolution 200×100), and

level 3: refined by 4 (effective resolution 800×400).

For this particular example, the uniform grid test took roughly 1300 s on a desktop machine, while the AMR calculation required 360 s. For this axisymmetric calculation much of the grid was refined since the hump extends around the sphere. For modelling a more localized region the benefits of AMR would be more pronounced.

The AMR results exhibit comparable accuracy to the uniform fine grid within the region of maximum refinement. This can be seen in the plots of max-norm error versus time in figure 7, which track the error in the Gaussian peak and are essentially the same with or without refinement. The peak is cut off in figure 6 in order to show the accuracy better in regions away from the peak.

In regions where the grid is less refined, such as the upper hemisphere at later times, the AMR calculation gives somewhat reduced accuracy but still quite acceptable given the coarseness of the grid. Note in particular that there is no spurious reflection of wave energy near the edges of the refinement patches.

Also note that in all cases the scatter plot looks as smooth near *θ*=0 as elsewhere, indicating that the highly skewed cells on the equator are not causing accuracy problems. Finally, note that no spurious waves are generated from the region near *θ*=−*π*/6≈−0.5 where the large topography feature is located. This test problem is analogous to the one introduced by LeVeque (1998) for one-dimensional shallow water equations, where it is demonstrated that all accuracy is lost if a fractional step procedure is used for the topography source term. The ability to capture this small amplitude wave depends on our use of the *f*-wave technique.

We have also tested the code with a wave of amplitude *A*=5000, in which case the Gaussian rapidly steepens into a shock wave. The high-resolution shock capturing algorithms employed in GeoClaw can handle this case as well. Results analogous to figure 6 have been obtained (not shown here).

## 8. Conclusions

We have demonstrated that the quadrilateral grids introduced in Calhoun *et al.* (2008*b*) can be successfully used to solve the shallow water equations on the surface of the sphere with adaptive mesh refinement. This has been accomplished using GeoClaw by enhancing this software package with a new choice of numerical boundary conditions needed for our quadrilateral geometry, as described in §5. This software uses high-resolution shock capturing finite volume methods that accurately capture shock wave propagation. We have also implemented the numerical method using the *f*-wave formulation described in §4, which ensures that the code is well balanced for an ocean at rest. Waves of small amplitude relative to the variation in the bottom topography are well captured, as illustrated in figure 6. Our quadrilateral grid has nearly constant cell areas, leading to efficient solution with explicit methods having Courant number between 0.5 and 1 in all cells. A potential difficulty is that some cells are highly skewed, but with the multi-dimensional wave propagation algorithms we employ, this does not appear to cause any loss of accuracy.

The computer code used to generate the numerical results and figures presented in this paper, along with links to the GeoClaw software, can be found at www.clawpack.org/links/sphere.

## Acknowledgements

This work was supported in part by DOE grant DE-FG02-88ER25053, AFOSR grant FA9550-06-1-0203, NSF grant DMS-0106511 and DFG grant HE 4858/1-1.

## Footnotes

One contribution of 10 to a Theme Issue ‘Mesh generation and mesh adaptation for large-scale Earth-system modelling’.

- © 2009 The Royal Society