# A radial basis function method for the shallow water equations on a sphere

Natasha Flyer, Grady B. Wright

## Abstract

The paper derives the first known numerical shallow water model on the sphere using radial basis function (RBF) spatial discretization, a novel numerical methodology that does not require any grid or mesh. In order to perform a study with regard to its spatial and temporal errors, two nonlinear test cases with known analytical solutions are considered. The first is a global steady-state flow with a compactly supported velocity field, while the second is an unsteady flow where features in the flow must be kept intact without dispersion. This behaviour is achieved by introducing forcing terms in the shallow water equations. Error and time stability studies are performed, both as the number of nodes are uniformly increased and the shape parameter of the RBF is varied, especially in the flat basis function limit. Results show that the RBF method is spectral, giving exceptionally high accuracy for low number of basis functions while being able to take unusually large time steps. In order to put it in the context of other commonly used global spectral methods on a sphere, comparisons are given with respect to spherical harmonics, double Fourier series and spectral element methods.

## 1. Introduction

Radial basis functions (RBFs) have the advantage of achieving spectral accuracy in multidimensions for arbitrary node layouts with extreme algorithmic simplicity. For the purposes of interpolating multidimensional surfaces, the methodology has been around for nearly 40 years. However, it is only in the last 15 years that it has been applied to solving mixed partial differential equations (PDEs) containing parabolic and/or elliptic operators (cf. Kansa 1990; Hon & Mao 1998; Hon et al. 1999; Larsson & Fornberg 2003; Wright & Fornberg 2006; Driscoll & Heryundono 2007). Furthermore, with regard to spherical domains, it has only been considered for solving a scalar PDE for strictly linear operators in the last 5 years (Lê Gia 2004, 2005), with the first application to a purely hyperbolic operator done in 2007 (Flyer & Wright 2007). Given this perspective, the next obvious step would be to solve a system of coupled nonlinear PDEs on the sphere. Since the shallow water equations occur in a multitude of applications, they provide any ideal test bed to establish the viability of the RBF methodology in terms of numerical accuracy and stability in this setting.

Since no precursor in the RBF literature exists with regard to solving a system of coupled nonlinear PDEs on a sphere, it should be emphasized that the goal of this paper is to establish its performance and viability (both for larger values of the shape parameter and in the flat basis function limit). As a result, the paper (a) develops the first known RBF method for the shallow water equations on a sphere and, in doing so, gives the first application of RBFs to a system of coupled nonlinear PDEs on a sphere; (b) demonstrates and analyses why RBFs can take unusually long time steps; and (c) shows that RBFs give high accuracy when compared with the other spectrally accurate methods, when the same degrees of freedom are used.

RBFs can be a competitive methodology, giving high numerical accuracy when compared with other high-order spectrally accurate methods.

The methodology is tested on two cases in which the analytical solutions are known in order to be able to perform an accurate error analysis. The first test is a global steady-state flow where the initial (and for all time) velocity field is compactly supported (i.e. non-zero in a limited band region), admitting a solution that is represented as an infinite spherical harmonic (SH) expansion. The second test is more challenging as it represents unsteady flow, modelling nonlinear advection of the initial condition that must be kept intact over time with minimal dispersion. Physically, it can be described as a forced nonlinear system with a translating low-pressure centre that is superimposed on a jet stream. An overview of the paper is as follows: §2 gives an introduction to RBFs; §3 provides a simple derivation of the shallow water equations for the sphere in Cartesian coordinates; §4 derives the discrete RBF formulation of the shallow water equations that is used in the test cases; §5 derives the linearized shallow water equations that will be used in the time stability analysis; §6 are the numerical studies; §7 gives timed benchmarks and §8 summarizes the paper with future prospects.

## 2. Introduction to RBFs

The motivation of the RBF methodology originated with Hardy (1971) asking the question, ‘Given a set of sparse scattered data, , at the node locations , can an interpolant be constructed that adequately represents the unknown surface?’. It was first shown by Mairhuber (1956) that, in more than one dimension, interpolation by an expansion of basis functions, , , which are independent of the node locations, is not well posed. That is, there exists an infinite number of node configurations that will yield a singular interpolation problem. Hardy bypassed this singularity problem with a novel approach in which the interpolant is constructed from linear combinations of a single basis function that is radially symmetric about its centre and whose argument is dependent on the node locations. By giving up orthogonality, well-posedness of the interpolant and its derivatives for any set of distinct scattered nodes in any dimension is gained.

Piecewise smooth RBFs feature a jump in some derivative and thus can only lead to algebraic convergence. For example, the radial cubic |r|3, where is the Euclidian norm, has a jump in the third derivative, at x=xj, leading to the fourth-order convergence in one dimension, with the order of convergence increasing as the dimension increases (cf. Powell 1992). On the other hand, the evidence strongly suggests that infinitely smooth RBFs, such as , and , will lead to spectral convergence (Madych & Nelson 1992; Yoon 2001), as is demonstrated in this paper. Note that the infinitely smooth RBFs depend on a shape parameter ϵ. It was first shown by Driscoll & Fornberg (2002) that, in one dimension, in the limit of ϵ→0 (i.e. flat RBFs), the RBF methodology reproduces pseudospectral (PS) methods if the nodes are accordingly placed (i.e. equispaced nodes for Fourier methods, Gauss–Chebyshev nodes for Chebyshev methods, etc.). Similarly, on the surface of a sphere, Fornberg & Piret (2007) showed that in the limit of ϵ→0, RBFs reproduce SH in the sense that they span an equivalent space for any scattered node set.

The accuracy of an RBF approximation can be improved by increasing the number of terms in the expansion and/or decreasing the shape parameter ϵ (Larsson & Fornberg 2003). In both cases, the shifted RBFs in the expansion approach the constant 1 as is clearly seen in the case for ϵ=0.1 in figure 1, which leads to ill-conditioning. The Contour–Padé algorithm (Fornberg & Wright 2004) can be used to bypass the RBF ill-conditioning mentioned above for the case of a fixed (relatively small) number of terms and increasingly small values of ϵ (even ϵ=0). Furthermore, Fornberg & Piret (2007) have recently developed an algorithm (RBF-QR) for bypassing the ill-conditioning for RBF interpolation on the surface of the sphere both as the number of terms is increased and ϵ is decreased right to 0. We will implement the RBF-QR to study how the error and time stability vary as a function of the shape parameter ϵ and what are the optimal choices for it.

Figure 1

GA RBF for various values of the shape parameter ϵ. ϵ=10, solid line; ϵ=1, dot-dashed line; ϵ=0.1, dashed line.

A good way to introduce the RBF methodology is through interpolation, since the differentiation matrices can be obtained by applying the exact differential operator to the interpolant (2.1) and then evaluating it at the data locations. As mentioned above, RBFs approximate a function f(x) sampled at some set of N distinct node locations by translates of a single radially symmetric function ϕ(r). For example, given the data values at the node locations , the RBF interpolant s(x) to the data is defined by(2.1)where the expansion coefficients, , are found by enforcing the collocation conditions such that the residual is 0 at the data locations. This is equivalent to solving the symmetric linear system of equations(2.2)where A is the interpolation matrix. An example of RBF interpolation in two dimensions is illustrated in figure 2. For RBFs such as the Gaussian (GA), inverse multiquadric (MQ) and inverse quadratic, the matrix in (2.2) is positive definite regardless of the distinct node locations and the dimension. For details on the well-posedness of (2.2) see Cheney & Light (2000, chs 12–16).

Figure 2

(a) Data values , (b) the RBF collocation functions and (c) the resulting interpolant.

Studies have shown that if the shape parameter, ϵ, is kept fixed throughout the domain (as will be done in the current study—variable shape parameter is needed when implementing local mesh refinement; Wertz et al. 2006; Driscoll & Heryundono 2007; Fornberg & Zuev 2007; Flyer & Lehto submitted), best results are achieved with roughly evenly distributed nodes (Iske 2004). Since only a maximum of 20 nodes can be evenly distributed on a sphere, there are a multitude of algorithms to define ‘even’ distribution for larger numbers of nodes, such as equal partitioned area, convex hull approaches and electrostatic repulsion (Sherwood 2007). Although any of these will suffice, we have decided to use an electrostatic repulsion or minimum energy (ME) approach since the nodes do not line up along any vertices or lines, emphasizing the arbitrary node layout and coordinate-free nature of a RBF methodology.

The ME node sets have the property that h, a measure of the spacing of the nodes, decays approximately uniformly similar to the inverse of the square root of the number of nodes N (Womersley & Sloan 2003), i.e. . Thus, they are similar to a uniform discretization of the unit square. In figure 3, the distribution for 1849 nodes on the unit sphere is displayed. For infinitely smooth RBFs, Jetter et al. (1999) show that, provided the underlying function being interpolated is sufficiently smooth, RBF interpolants converge in the L norm such as h−1/2 ec/4h for some constant c>0 that depends on the RBF. For the ME node sets, convergence will thus proceed as .

Figure 3

Minimum energy (ME) nodes on the sphere, N=1849.

## 3. Cartesian form of the shallow water equations on the sphere

As was discussed in Flyer & Wright (2007), the RBF differentiation matrix approximating a continuous operator is completely independent of the orientation and type of coordinate system in which the operator was originally posed. For example, the gradient operator has no singularities on the surface of the sphere. However, if we choose to represent it in spherical coordinates (λ, θ), where λ and θ are the respective longitudinal and latitudinal variables, singularities will occur at the poles (θ=±(π/2)) simply due to the coordinate system. Unlike orthogonal expansions on a latitude–longitude grid, we can apply the surface gradient operator to a RBF and there will be no evidence of the pole singularities as the operator, itself, is inherently smooth on the sphere. For PDEs, where the scalar variable is always acted on by a smooth spatial operator, one can easily pose the equation in spherical coordinates and then apply the RBF methodology (as was done in Flyer & Wright (2007) and Fornberg & Piret (2008)). However, when solving the shallow water equations this is not the case. If a spherical coordinate system were to be used, the directional velocity vector components, u (latitudinal) and v (longitudinal), will inherently carry the pole singularities in their solution since the unit vector is discontinuous at the poles. As a result, the Cartesian form of the shallow water equations will be used since such a coordinate system has no singularities.

The full shallow water equations in a three-dimensional Cartesian coordinate system for a rotating fluid are as follows:(3.1)(3.2)where f is the Coriolis force; and are the velocity vectors; h is the geopotential height; and x={x, y, z}T represents the position vector.

Next, we confine the flow to the surface of the unit sphere, which first requires constructing a linear operator, P, for projecting vectors in a three-dimensional Cartesian space onto a plane tangent to the sphere. The construction of P is quite simple since, if (x, y, z) is a point on the unit sphere, then the normal at (x, y, z) is just x. Thus, if u is a vector at (x, y, z), then xxTu gives the projection of u onto x, and uxxTu gives the projection of u onto the plane tangent to the sphere at (x, y, z). As a result, we define the projection operator P as,(3.3)where I is the 3×3 identity matrix. The px, py and pz represent the projection operators in the x, y and z directions, respectively.

Next, the gradient operator appearing in (3.1) must be constrained so that when it is applied to a scalar, it produces a vector that is tangent to the unit sphere. This is done by replacing all occurrences of ∇ in (3.1) with the operatorsince P projects arbitrary three-dimensional Cartesian vectors onto a plane tangent to the unit sphere. The divergence operator appearing in (3.2) must also be restricted, so that it produces the surface divergence of a vector field. This is also accomplished with P by simply taking the dot product of the vector field with the projected gradient operator P∇.

The final step in confining the motion to the unit sphere is to project the right-hand side (r.h.s.) of (3.1), with the modified differential operators, onto the corresponding , and directions. For example, in the case of the u momentum equation, this is done by taking the dot product of the modified r.h.s. of (3.1) with the px projection vector.

Putting all these pieces together, the final shallow water equations on the surface of the unit sphere in Cartesian coordinates is given by(3.4a)(3.4b)Note that the only spatial operator that needs to be discretized is the projected gradient, P∇, and its components, px.∇, py.∇, pz.∇.

## 4. The RBF method on the sphere

In §4ac, we develop the discrete operators that are necessary to spatially discretize the shallow water equations using RBFs. In §4a, we define how to formulate the projected gradient operator in RBFs. In §4b, we do a step-by-step construction of the three needed differentiation matrices, , and . Finally, §4c demonstrates the entire RBF formulation of the shallow water equations to be implemented computationally.

### (a) Formulation of the projected gradient operator in RBFs

Let x={x, y, z} and be points on the unit sphere, then the Euclidean distance from x to xk is(4.1)It is important to note that the distances are not great circle arcs measured along the surface, but are the Euclidean distances measured straight through the sphere. The reason being is that RBFs do not ‘feel’ the geometry of the domain in which they are applied nor the dimension, but only the scalar distances between the nodes and the locations at which the RBFs are centred (in our study the RBFs centres and the nodes coincide, but this need not be the case).

Let ϕk(r(x)) be an RBF centred at xk. Using the chain rule, the gradient of ϕk(r(x)) is given by(4.2)where we have used ′ to denote differentiation with respect to r. The projected RBF gradient operator is then simply obtained by applying the projection matrix P=IxxT to (4.2)Now, we have all the components necessary to build the action of on an RBF representation of a geophysical quantity.

### (b) Constructing the RBF differentiation matrices

Given a geophysical quantity f(x)=f(x, y, z) known at the node locations on the surface of a unit sphere, we first represent f(x) as an RBF expansion(4.3)where the coefficients ck are determined by collocation. Applying the projected gradient operator in the x direction to (4.3) and evaluating it at gives(4.4)where A is the interpolation matrix in (2.2) and we have used the property that c=A−1f. The differentiation matrix is the discrete RBF approximation to the x-component of the projected gradient operator. Using the same concept, and are derived as follows:(4.5)(4.6)where the entries in the matrices By and Bz are respectively given byNote that all three differentiation operators (4.4)–(4.6) are well defined with no pole singularities.

### (c) The discretized RBF shallow water equations

Let be the geopotential height and be the components of the velocity field u sampled at the node locations . Then, the r.h.s. of the Cartesian form of the shallow water equations in §3 can now be written in discretized form r.h.s.D aswhere ∘ denotes element by element multiplication of the vectors. The full discretized equations are given by(4.7)(4.8)where , , are the vectors of px, py, pz evaluated at the nodes . The method-of-lines technique is used to advance the system in time. As a note, for flows on spheres of radius a, the modification to the method is simple—leave x and xj as points on the unit sphere and divide , and by a.

## 5. Eigenvalue stability: the linearized equations

For each test case, we will perform an eigenvalue stability analysis in order to understand: (i) why the RBF method can take long time steps and (ii) the impact of ϵ-refinement on stability. However, in order to do so, we must linearize the shallow water equations about a given state. In all cases, we linearize about the initial condition, since this is the solution for all time for the steady-state case and for unsteady flow the shape of the initial condition is advected intact.

Assume and are the approximate solutions to (3.4a) and (3.4b) and can be written aswhere u0 and h0 are the initial conditions. Substituting this expansion into (3.4a) and (3.4b) and only considering perturbations of O(δ) gives the linearized equations(5.1a)

(5.1b)

(5.1c)

(5.1d)

## 6. Numerical studies

For each test case, we evaluate the convergence properties of the RBF method both with respect to uniformly increasing the number of nodes (h-refinement) and decreasing ϵ (ϵ-refinement). We also perform a stability/eigenvalue study to explore how the time errors behave as a function of ϵ and the time step. Also, how well linear stability analysis predicts the growth of unstable modes in the full RBF solution is explored. Finally, the performance of the RBF method is put into context with the following commonly used spectral models in the literature: a de-aliased SH method (Jakob-Chien et al. 1995), a double Fourier with a SH filter (DF/SHF; Spotz et al. 1998) and a spectral element (SE) method with filtering (Taylor et al. 1997). The RBF method does not need to use any spatial filtering to achieve the numerical accuracies reported.

### (a) Steady-state test case: compactly supported wind field

In this test case, the flow field is nonlinear and compactly supported but still infinitely differentiable (see Williamson et al. (1992), test case 3). For simplicity, the initial condition and solution for all time for the wind field in spherical coordinates is given by(6.1)where , θb=−π/6, θe=π/2 and . usph(θ) has a maximum at , corresponding to a maximum wind positioned at 30° N (π/6) and is 0 outside of θb<θ<θe. The latitudinal component of the velocity field is 0. The analytic solution is then defined by a rotation of the coordinate system to tilt the wind field an angle α relative to the polar axis. For α=π/3, the wind field is graphed in figure 4a as an orthographic projection centred at 60° S and 0° E. The geopotential height field is defined by numerically integrating to machine precision(6.2)where f=2Ω sin θ, Ω=7.292×10−5 s−1 and a=6.37122×106 m. The geopotential height field is graphed in figure 4b. The test case is run for 5 days.

Figure 4

Analytical solution for the steady-state test case in the rotated coordinate system. (a) The velocity field sampled at the N=3136 ME nodes; (b) the contours of the height field from 2100 to 3000 m in increments of 100 m. Both plots are orthographic projections centred at 60° S (−π/3) and 0° E, which is what the solution would look as in non-rotated coordinates.

#### (i) Error study with regard to h-refinement

Figure 5 displays the error (exact minus numerical) in the RBF solution of the geopotential height field on an unrolled sphere for N=1849 (figure 5a) and N=3136 (figure 5b), respectively. The contour lines represent the exact solution in 100 m intervals, beginning at 3000 m. In figure 5a, all deviations in the solution of less than 10−4 m are shown in white, while in figure 5b deviations from the exact solution of less than 10−6 m are shown in white. As expected, both figures demonstrate that the dominant error is concentrated in the area with the steepest gradients. This test would be a prime candidate for local node refinement, as clustering the nodes in such areas would have a high likelihood of severely decreasing the error as has been demonstrated in Flyer & Lehto (submitted). Furthermore, note that there is little evidence of any dispersive wave trains in either plots, indicating that the r.h.s. of the PDEs in (4.7) are being adequately resolved.

Figure 5

The error (exact-numerical) for the height field from the steady-state test case at t=5 days in ‘unrolled’ spheres (a) N=1849 nodes and (b) N=3136 nodes. Solid black lines correspond to contours of the exact solution plotted from 2100 to 3000 m in increments of 100 m. ϵ=3.25 for all N.

Figure 6a shows the time traces of the relative 2 error for the 5 day simulation for various values of N and fixed ϵ=3.25. Note that the error for all values of N barely grows as a function of time during the simulation period. To determine spectral convergence from the time traces, we plot the relative 2 error as a function of spatial resolution (i.e. ) in figure 6b and see that indeed it is spectral. The figure also demonstrates that the 2 error is relatively insensitive to the value of ϵ used between 3.25 and 4. For larger nodes sets, if ϵ goes below 3.25 ill-conditioning will set in. More on the sensitivity of the calculation to ϵ will be addressed in §6a(ii) on ϵ-refinement.

Figure 6

(a) Relative 2 error in the height field for the steady-state test case as a function of time and N. ϵ=3.25 for all N. (b) Relative 2 error in the height field for the steady-state test case at t=5 days as a function and ϵ. ϵ=4, cross; ϵ=3.75, circles; ϵ=3.5, squares; ϵ=3.25, diamonds.

Figure 7 shows the relative 1, 2 and norms of the error for the RBF method as a function of time for N=3136 and Δt=10 min. This spatial resolution was selected since it results in a relative 2 error of O(10−10), which is also the order of the lowest relative 2 error reported by many of the other spectral methods (table 1). For this time step and node set, figure 8 shows the time traces as a function of ϵ. One can see that as ϵ decreases the error decreases. However, the rate at which the error decreases slows down as the value of ϵ=3.25 is approached. One can then ask, is this the effect of ill-conditioning as the shape of the RBFs become flat? The next logical question would then be, ‘what if we implemented the RBF-QR algorithm (Fornberg & Piret 2007), which bypasses this ill-conditioning problem and allows for calculations in the limit as ϵ→0. Will we see the error drastically decrease?’

Figure 7

Relative 1, 2 and errors in the height field for the steady-state test case as a function of time for N=3136, ϵ=3.25 and Δt=10 min.

View this table:
Table 1

Comparison of commonly used spectral methods for steady-state flow with α=π/3. (The number in parentheses in the RBF section corresponds to the square root of N, which is inversely proportional to the node spacing. The number in parentheses in the SH section corresponds to the number of SH coefficients updated in time. RBF and DF use the same time-stepping scheme. SE uses a third-order Adams–Bashforth method. SH uses a semi-implicit time-stepping scheme denoted by the asterisks. For the SH 1849 case (Spotz et al. 1998), gives Δt=3 min when using a leapfrog scheme as is done in RBF and DF.)

Figure 8

Relative 2 errors in the height field for the steady-state test case as a function of time and different values of ϵ.

#### (ii) Error and time stability study with regard to ϵ-refinement

Our objective is to explore the accuracy and stability of the RBF shallow water method as the shape parameter ϵ decreases to zero. Since the solution to this test case can be accurately resolved with a low-resolution model, as demonstrated in §6a(i), we fix the node set at N=484 for the study. To compute the RBF differentiation matrices for this resolution in the low ϵ regime (i.e. 0≤ϵ≤0.7 for N=484), we use the RBF-QR algorithm.

Figure 9 shows time traces for the 5 day run of the relative 2 error for ϵ varying from 0 to 0.7. Leapfrog with a Robert filter of γ=0.07 and a time step of 18 min was used. Clearly, time instability sets in as ϵ tends to zero, with the worse case being ϵ=0 corresponding to SH as previously mentioned. There are an infinite set of node layouts, such as a latitude–longitude grids or any perturbation to them, that fall on or close to the zeros of the SH. The result is singular or close to singular SH interpolation matrices. As a result, a least-squares approach or over-sampling is needed in order to ‘pin down’ a solution when using this basis. Thus, one conclusion is that the optimal ϵ is not in the ϵ→0 limit. However, one possibility as to why instability sets in as ϵ becomes small is that we are not using a small enough time step. However, figure 10 easily throws that notion aside as we see that whether the time step is 2, 6 or 18 min, the blow-up of the 2 error is identical.

Figure 9

Relative 2 error in the height field for steady-state test case. The different curves correspond to different values of ϵ used in the RBF solution, with the top curve corresponding to ϵ=0 and the bottom curve to ϵ=0.7 in increments of 0.05.

Figure 10

Relative 2 error in the height field for steady-state test case. ϵ is fixed at 0.2 and the time step is varied as indicated in the legend. Δt=120 s, cross; Δt=360 s, circles; Δt=1080 s, squares.

A necessary condition for the RBF method to be stable for all time is for the spectrum of the RBF operator for the entire r.h.s. of the linearized system (5.1a)–(5.1d) to be contained within the stability domain of the time-stepping scheme. Figure 11 plots the stability domain for the leapfrog scheme with a Robert's filter of γ=0.07 along with the spectrum of the RBF operator for the linearized equations for ϵ=0, 0.2 and 0.7. The most important trend to note is that as ϵ decreases the real part of the eigenvalues spread off the imaginary axis and into the right half plane both near the centre and at the ends of the spectrum. As was noted in figure 9, the worst case is for ϵ=0 when RBFs reproduce SH. However, if we look at the eigenvector associated with the eigenvalue having the largest real part in this case, as displayed in figure 12, we see such modes are spurious and not physical. As is well known in the SH literature, one needs to impose filters to take care of such eigenmodes. However, for the N=484 case, the eigenmodes associated with the eigenvalues that are slightly off the imaginary axis in figure 11c do not come into play for time integration periods in the order of a month.

Figure 11

Spectrums of the RBF approximation operator of the r.h.s. of the linearized system in equations (5.1a) and (5.1b) for different values of ϵ. Each spectrum from (a) to (c) has been scaled by the time step Δt=1080 s (18 min), where ξt×eigenvalues. The stability domain of time-stepping scheme, leapfrog scheme with a Robert filter of γ=0.07, is plotted as a grey line. (a) ϵ=0; (b) ϵ=0.2; and (c) ϵ=0.7. For all figures N=484.

Figure 12

The unit eigenvector for the height field corresponding the eigenvalue with the largest real part in figure 11a.

We would also like to see how well the linear eigenvalue stability analysis predicts the growth of the unstable modes in the full RBF solution for this test case. If it is a good predictor, then the error in the RBF solution for a given ϵ and Δt should grow as(6.3)where β denotes the eigenvalue of the linearized RBF operator with the largest real part, scaled by Δt, which is not inside the time stability domain and k is the number of time steps taken. In figure 13a–c, we plot this predicted growth rate together with the error in the RBF solution for the height field for the same three values of ϵ used in the previous figures (ϵ=0, 0.2, 0.7). Recall that for the first two of these ϵ values the solution went unstable well before t=5 days, while for the last case the solution was stable and accurate. In figure 13a, we see the predicted growth matches the observed growth very nicely. In figures 13b,c, the linear stability theory actually predicts a slightly more elevated result than what is observed.

Figure 13

Relative 2 error in the RBF solution of the geopotential height for the steady-state test case together with the predicted growth rate (solid line) from equation (4.5). The value of ϵ (cross) is as follows: (a) ϵ=0; (b) ϵ=0.2; and (c) ϵ=0.7. For all figures Δt=1080 s and N=484.

#### (iii) Discussion of spectral models and comparison

The methods presented in tables 1 and 2 are described below (Jakob-Chien et al. 1995; Taylor et al. 1997; Spotz et al. 1998).

1. SH requires twice as many grid points as basis functions. Owing to a severe Courant–Friedrichs–Lewy (CFL) condition, a semi-implicit time stepping is commonly used (Jakob-Chien et al. 1995). Also, SH are notorious for incorrectly increasing the energy in low modes through aliasing. As a result, SH are never run without de-aliasing, normally using Orszag's 2/3 rule (Gottlieb & Orszag 1977). So, for example, a SH method that uses 1849 basis functions is in fact using 4096 basis functions with a grid of 8192 nodes. However, only 1849 bases are updated in time in spectral space due to de-aliasing. This leads to the forward SH transform (which is O(M3) operations, where M is the number of associated Legendre roots in the latitudinal direction), being performed on 4096 bases using 8192 grid point values, and the reverse transform on 1849 spectral coefficients since the others are set to zero. Fast Fourier transforms (FFTs) can be used in the longitudinal direction, resulting in an operation count of O(N3/2) per time step (N is the total number of grid points). Lastly, they involve high computational complexity.

View this table:
Table 2

Comparison between commonly used spectral methods for the unsteady flow test case. (The number in parentheses in the RBF section corresponds to the square root of N, which is inversely proportional to the node spacing. The number in parentheses in the SH section corresponds to the number of SH coefficients updated in time. RBF and DF use the same time-stepping scheme. SE uses a third-order Adams–Bashforth method. SH uses a semi-implicit time-stepping scheme denoted by the asterisks. For the SH 1849 case (Spotz et al. 1998), gives Δt=3 min when using a leapfrog scheme as is done in RBF and DF.)

2. DF. Unlike SH, DF has a one-to-one correspondence between grid points and spectral coefficients. For comparison, a SH method that uses 1849 bases with the necessary 8192 grid points is approximately equivalent to a DF using 8192 grid points (table 1). DF for the same reasons as SH also requires the use of some type of de-aliasing filter. Neither SH nor DF have the option of local mesh refinement. Owing to FFTs, DF has only an operation count of O(N log N) per time step.

3. SE. Grid generation is a large overhead cost for SE and include a variety of approaches such as icosahedral or the projection of a cube onto the surface of a sphere. SE allows for local node refinement but implementation is algorithmically highly complex. The method considered here (a cubed sphere; Taylor et al. 1997) uses a tensor product of Legendre polynomials on each element. Owing to the clustering of Legendre grid points near the ends of each element boundary, the CFL restriction is severe and increases as the order of polynomials used increases. SE models generally use some type of filtering or instability will set in. The operation count for SE is O(kP2), where k is the number of elements and P is the order of the polynomial expansion on each element.

4. RBFs. RBFs are algorithmically very simple. The RBF method presented in this paper can be coded in less than 100 lines of Matlab, using only level 1–3 BLAS operations. Owing to close clustering of the eigenvalues of the linearized RBF operator about the imaginary axis as seen in §6a(ii), RBFs can take comparatively much longer time steps (table 1). Unlike the other methods, the RBF method can go without de-aliasing for longer time scales. Since aliasing is a spectral phenomena where all information is global, it is hypothesized that the locality properties of RBF coefficients (Fornberg et al. 2008) and the fact that they only pick up information locally (Fornberg & Flyer 2005; Fornberg et al. in press) contributes to this phenomena. Furthermore, for a given spatial resolution, RBFs achieve a higher numerical accuracy (table 1). To implement global RBFs, as done in this paper, O(N2) operations are used. However, reducing this cost is discussed at the end of the paper and is currently under implementation.

For all RBF results listed in table 1, MQ is used with ϵ=3.25 (although any infinitely smooth RBFs can be used, but GA tends to be more sensitive to the value of ϵ). Leapfrog with a Robert's filter of γ=0.07 is used to advance the model in time. All results are for α=π/3, tilting the wind field in (6.1) so that it flows directly over the poles. For all N in table 1, RBFs can easily take longer time steps with regard to stability but the time step was chosen so that temporal discretization errors matched spatial discretization errors, which is the optimal scenario. The RBF case is the only method that does not use a de-aliasing/spatial filter.

### (b) Unsteady flow test case

The test case is a low-pressure system, initially centred at (λ0, θ0)=(0, π/4), superimposed on a jet stream that is symmetrical about the equator (see test case 4 in Williamson et al. (1992)). Forcing terms are added to the shallow water equations to constrain the motion of the system, so that the initial condition is nonlinearly advected intact.

#### (i) Derivation of analytical solution in Cartesian form

Since the Cartesian description of the solution does not appear to have been given before in the literature, we first discuss this topic.

Let xc(t) denote the centre of the translating low-pressure system in Cartesian coordinates at time t. Then, xc(t) is given bywhere a is the radius of the Earth and (λ0, θ0) is the initial position of the low-pressure system. Following the notation of Williamson et al. (1992), let denote the stream function for the low-pressure system at time t without the superimposed jet streamwhere ψ0 and σ are constants defined below. With this definition, the velocity field in Cartesian coordinates is given byThe first term in this expression is the velocity field for the jet stream with maximum velocity u0, while the second term is the velocity field for the low-pressure system. The geopotential height is given by(6.4)where f=2Ωz is the Coriolis force andThe constants are specified as follows: (θ0, λ0)=(π/4, 0); ψ0=−0.03(gh0/f0); σ=(12.74244)2; gh0=105 m2 s−2; Ω=7.292×10−5 s−1; ; and the maximum velocity, u0 is 20 m s−1. Figure 14a,b display the initial (t=0) velocity and height fields, respectively. The solution at any time t looks identical to these figures, but with the centre of the low-pressure system shifted to λ=(u0/a)t.

Figure 14

Initial (a) velocity field and (b) height field with N=3136 for the unsteady flow test case plotted as orthographic projections centred at 45° N and 0° E. The contours in (a) range from 10 600 to 10 100 m in intervals of 50 m.

These values of and are not analytical solutions of the shallow water equations (3.4a) and (3.4b), but can be made analytical solutions by adding forcing terms Fu, Fv, Fw and Fh to the respective r.h.s.s of (3.4a) and (3.4b). The Fu forcing term for the u momentum equation (3.4a) is given by substituting and into the r.h.s. of (3.4a) and then subtracting this from the time derivative of , i.e.The forcing terms Fv, Fw and Fh are similarly computed through the respective substitution of the analytical solutions and into (3.4a) and (3.4b). As specified in Williamson et al. (1992), the complete forced shallow water system is to be simulated for 5 days. The standard fourth-order Runge–Kutta (RK4) scheme is used to advance the RBF method in time. This was chosen over leapfrog with a Robert filter since the unsteady test case is much more sensitive to time truncation errors.

#### (ii) Error study with regard to h-refinement

Figure 15 displays the error (exact minus numerical) for the RBF solution of the height field for N=1849 and 3136 nodes. It is clear that the largest errors are predominantly located where the gradients in the solution are the highest as would be expected. For the latter case (N=3136), the time traces of the 1, 2 and error in the height field over the full 5 day integration period with a Δt=15 min are given in figure 16. The trend in the growth of the 1 error observed in this figure is similar to that reported in the SE method (Taylor et al. 1997).

Figure 15

The error (exact − numerical) at t=5 days for the height field from the unsteady flow test case. The results are for (a) N=1849 and (b) N=3136 nodes and ϵ=3.25. Solid lines correspond to contours of the exact solution over the same values as figure 14b.

Figure 16

Relative 1, 2 and errors in the height field for the unsteady flow test as a function of time for the case of N=3136, ϵ=3.25 and Δt=15 min.

Figure 17a shows that the relative 2 error for various N barely grows as a function of time during the simulation period. To determine the convergence rate from these time traces, we plot the relative 2 error as a function of spatial resolution (i.e. ) in figure 17b and see that the RBF method is converging spectrally to the true solution of the height field. In contrast to figure 6b, figure 17b shows that the relative 2 error is more sensitive to the value of ϵ, decreasing almost an order of magnitude as ϵ decreases from ϵ=4 to 3.25 for N=5041. Even for smaller values of N, such as N=3136, we see this pattern, as is illustrated by the time traces of the relative 2 error for varying ϵ values in figure 18. As in the steady-state test case, we again ask the question, ‘Will the error decrease as ϵ keeps decreasing?’. We will again be employing the RBF-QR algorithm (Fornberg & Piret 2008) to answer this question.

Figure 17

(a) Relative 2 error in the approximate solution of the height field for the unsteady flow test case 3 as a function of time and different spatial resolutions. The shape parameter was fixed at ϵ=3.25 for all N values. (b) Relative 2 error in the approximate solution of the height field for the unsteady flow test case at t=5 days for different spatial resolutions and values of the shape parameter ϵ. ϵ=4, cross; ϵ=3.75, circle; ϵ=3.5, squares; ϵ=3.25, diamonds.

Figure 18

Relative 2 errors in the height field for the unsteady flow test case with N=3136 as a function of time and different values of the shape parameter ϵ.

#### (iii) Error and time stability study with regard to ϵ-refinement

To study the accuracy and stability with respect to the shape parameter ϵ, we restrict our attention to the N=1849 case. This is a larger node set than used in the previous ϵ-refinement study due to proper resolution of the translating low-pressure system.

As in the previous test, we study time stability by analysing the spectrum of the RBF operator for the linearized system (5.1a)–(5.1d). However, unlike the previous test case, the linearized system changes every time step since the analytical solution is not a steady-state solution. In the results that follow, we will thus look at the spectrum for the case of t=0. Although not demonstrated here, the spectrum for other values of t is very similar, which should be expected since the analytical solution at time t is just a translation of the t=0 solution.

Figure 19a–d displays the stability domain for RK4 along with the spectrum of the RBF operator for the linearized equations for ϵ=0, 0.2, 0.7 and 1.7, respectively. For the first three of these values of ϵ, the RBF-QR algorithm was used to compute the differentiation matrices, while for ϵ=1.7 the direct approach of inverting the interpolation matrix A was used. Note that the ϵ=0 spectrum has been scaled by Δt=360 s, whereas the other spectrums have been scaled by Δt=1200 s. This change was necessary since the eigenvalue with the largest negative real part for the ϵ=0 case would not fit in the stability domain with Δt=1200 s. We see a similar trend to the steady-state test case, in that there is a terrible degradation of the eigenvalues into the right half plane both near the centre and at the ends of the spectrum as ϵ is decreased. However, for the ϵ=1.7 case, the eigenvalues are tightly compacted around the imaginary axis and do not spread far up the imaginary axis. This is the reason why RBFs give high accuracy with much longer time steps when compared with SE, which use orthogonal basis such as Legendre polynomials. Such bases have eigenvalues that are high up and low down in the left half plane, severely restricting the time steps (Fornberg 1996, §4.4).

Figure 19

Spectrums of the RBF approximation operator of the r.h.s. of the linearized system in equations (5.1a) and (5.1b) for different values of ϵ. The spectrum of (a) has been scaled by the time step Δt=360 s (6 min), while the spectrums from (b) to (d) have been scaled by Δt=1200 s (20 min) with ξt×eigenvalues. The stability domain of the RK4 time-stepping scheme has also been plotted as a grey line in each figure. (a) ϵ=0; (b) ϵ=0.2; (c) ϵ=0.7; and (d) ϵ=1.7.

To see how well the linear eigenvalue theory predicts the error growth of the unstable modes in the model, we plot in figure 20a–c the predicted growth according to (6.3) together with the 2 error of the height field for ϵ=0.2, 0.7 and 1.7. We do not include the results for ϵ=0 since the solution blows up after only eight time steps with Δt=360 s. Figure 20a can only be run for 100 time steps before the method goes completely unstable. All the plots show the eigenvalue analysis that predicts the growth relatively accurately. Note that in the plots from figure 20b,c, the x-axis ranges from 0 to 720 time steps, which is actually a 10 day integration period instead of 5 days as specified in the test case. In the latter of these plots, we see the growth rate from 5 days (360 time steps) to 10 days is quite acceptable.

Figure 20

Relative 2 error of the height field for the unsteady test case together with the predicted growth rate (solid line) from equation (6.3). Δt=1200 s (20 min) and N=1849. The value of ϵ (crosses) used is indicated below each figure. Figures (b) and (c) are run for 720 time steps, corresponding to 10 days of integration. (a) ϵ=7; (b) ϵ=0.7; and ϵ=1.7.

#### (iv) Comparative results

As with the previous test case, we compare the performance of the RBF method to several other high-order methods in table 2. From each study, the best results were chosen in terms of numerical accuracy. The 2 error is normalized with regard to the norm of the true height field without the mean term in (6.4). The larger RBF time steps in table 2 as opposed to table 1 are due to using the RK4, for which the stability domain extends further up the imaginary axis than the leapfrog scheme. In table 2, the time step was chosen at the break point where the temporal discretization error matched the spatial discretization error. For example, for the N=4096 case, the relative 2 error at the end of the 5 day run is graphed as a function of the time step in figure 21. The RBF method could easily have taken a larger timestep than 8 min (as reported in table 2), but temporal discretization errors would begin to dominate as shown by the large growth in the error after an 8 min time step. Again, the RBF method uses no spatial filtering for the reported results.

Figure 21

Relative 2 error in the height field at 5 days for the unsteady flow test case as a function of the time step for ϵ=3.25 and N=4096.

## 7. Performance benchmarks

In order to give the reader a feeling of the time requirements for running the RBF method, we gave some timing benchmarks in table 3 for both test cases. Although RBF matrices are full and thus their inversion to calculate the three needed differentiation matrices is O(N3), this is a pre-processing step that is done only once. At every time step, a matrix–vector multiply is needed, requiring O(N2) operations. Table 3 gives both the runtime per time step as well as total runtime needed to achieve the given 2 error (or similarly as a function of N). All computations were performed in Matlab v. 7.5 R2007b running on Dell PowerEdge 2950 Server with two 2.66 GHz Intel Xeon X5335 Quad Core processors. BLAS multithreading in Matlab was enabled with a maximum of eight cores available.

View this table:
Table 3

Runtime results for the RBF method for a 5 day simulation of steady-state and unsteady forced flows. (Time is given in seconds (s). See §7 for details.)

The runtimes are greater for the unsteady flow test case than steady-state flow since (i) the RK4 method was used as opposed to a leapfrog scheme due to greater sensitivity in time truncation errors and (ii) for unsteady flow, the forcing functions have to be evaluated at each time step.

## 8. Summary and future prospects

The paper develops a stable spectrally convergent RBF method for solving the shallow water equations on a sphere and evaluates its convergence and time stability properties using test cases where analytical solutions are known so that an exact error study in time and space could be performed. Results are put into context with respect to those published in the literature for other commonly used spectral methods. The general findings with regard to the steady and unsteady flow tests considered are as follows.

1. h-Refinement. For both test cases, spectral convergence is easily achieved. The error (exact numerical) is concentrated in the region where the solution has the steepest gradients, outside of which little dispersion was seen. For all values of N tested, the 2 error barely grew as function of time over the 5 day simulation period.

2. ϵ-Refinement. Interestingly, in the limit as ϵ→0 (with ϵ=0 being SH), time instability sets in. Performing a linearized eigenvalue stability analysis about the steady-state solution of the RBF scheme, eigenvalues are noted to spread off the imaginary axis and into the right half plane. However, away from this limit, eigenvalues are tightly clustered around the imaginary axis and near the origin in contrast to PS methods. This allows for long time steps. In both test cases, the linear eigenvalue stability analysis of the RBF method proved to be a decent measure of error growth.

3. Comparative study. For the steady flow test case, the RBF method achieves a relative 2 error of 6.88×10−12 with N=5041 nodes and a time step of 6 min using a leapfrog scheme with a Robert's filter of 0.07. The highest accuracy noted in the literature is 2(10−13) for a DF method when used with a SH filter with N=32 768 nodes and a time step of 90 s (Spotz et al. 1998), also using a leapfrog scheme. For the unsteady test case, the RBF method achieved a relative 2 error of 1.02×10−8 with N=5041 nodes and a time step of 6 min using the standard RK4. The highest accuracy noted in the literature is 4×10−5 for a SE method with N=24 576 nodes and a time step of 45 s (Taylor et al. 1997).

The RBF node sets used in this study are for roughly uniform resolutions of approximately 500 km or greater on the surface of the Earth. However, to go to much higher resolutions (e.g. on the order of 10 km) and to do three-dimensional modelling requires many more nodes. Since global RBFs (as used in this paper) require full matrices, they are not a practical computing option for high-resolution models in either two or three dimensions, especially with regard to parallelization. In these cases, RBF methods, such as RBF finite differences, may give an alternative viable approach and are currently under development for a variety of applications.

## Acknowledgments

The National Center for Atmospheric Research is sponsored by the National Science Foundation. N.F. was supported by NSF grant ATM-0620100. G.B.W. was supported by NSF grant ATM-0801309. The authors would like to thank Dr Amik St-Cyr, Dr Ram Nair and Dr Aime Fournier.

## Footnotes

• Received January 20, 2009.
• Accepted February 25, 2009.

View Abstract