High-Order Semi-Lagrangian Schemes for the Transport Equation on Icosahedron Spherical Grids

Abstract

The transport process is an important part of the research of fluid dynamics, especially when it comes to tracer advection in the atmosphere or ocean dynamics. In this paper, a series of high-order semi-Lagrangian methods for the transport process on the sphere are considered. The methods are formulated entirely in three-dimensional Cartesian coordinates, thus avoiding any apparent artificial singularities associated with surface-based coordinate systems. The underlying idea of the semi-Lagrangian method is to find the value of the field/tracer at the departure point through interpolating the values of its surrounding grid points to the departure point. The implementation of the semi-Lagrangian method is divided into the following two main procedures: finding the departure point by integrating the characteristic equation backward and then interpolate on the departure point. In the first procedure, three methods are utilized to solve the characteristic equation for the locations of departure points, including the commonly used midpoint-rule method and two explicit high-order Runge–Kutta (RK) methods. In the second one, for interpolation, four new methods are presented, including (1) linear interpolation; (2) polynomial fitting based on the least square method; (3) global radial basis function stencils (RBFs), and (4) local RBFs. For the latter two interpolation methods, we find that it is crucial to select an optimal value for the shape parameter of the basis function. A Gauss hill advection case is used to compare and contrast the methods in terms of their accuracy, and conservation properties. In addition, the proposed method is applied to standard test cases, which include solid body rotation, shear deformation of twin slotted cylinders, and the evolution of a moving vortex. It demonstrates that the proposed method could simulate all test cases with reasonable accuracy and efficiency.

1. Introduction

Advection transport is a basic process in fluid dynamics, especially when it comes to the motion of various passive tracers in the atmosphere or ocean. It plays an important role in geographic modeling and simulation for kinds of geographical phenomena and processes [1,2,3]. The global advection transport model is also important in developing general circulation models (GCMs). The traditional latitude–longitude grid is easy for application but has singularities at the two poles [4]. Moreover, its non-uniform grid system would also seriously affect computational efficiency [5]. To address these issues, quasi-uniform grid systems without singularities or with weak singularities, such as Yin-Yang overset grids, block-structured quadrilateral meshes based on the cubed sphere, triangle or their dual meshes based on the regular icosahedron sphere, are becoming more and more popular in developing global transport models [6,7]. Each of these grids has its own merits and demerits. The Yin-Yang grids are composed of two patches of the latitude–longitude grid and belong to the chimera grid [8], which needs a special method to handle the internal boundaries for data exchange. The cubed sphere is composed of six patches. The coordinate system is not continuous across boundaries shared by two neighboring patches [9]. The triangle mesh or its dual hexagon-pentagon mesh based on a regular icosahedron sphere has relatively more uniform space resolution and has been utilized as a basis in many applications, especially in global atmosphere or ocean modeling [10,11]. In this study, the triangle mesh based on the icosahedron sphere is adopted for our transport model.

There are mainly two classes of numerical methods for the approximation of the advection equation on the sphere. The first class of methods belongs to the finite volume (FV) methods [12,13], which consider the fluid as a set of cells and then solve the cell-integrated flux-form advection equation. The solution is just the mean value of each cell, which is computed from the net flux into a cell over a time step. An FV method typically computes these fluxes by assuming subgrid reconstruction within the cell and then integrating it over the area of the cell flowing outward during a time step. It needs delicate techniques to construct stable, high-order numerical flux on the cell interface, and the time step is constrained by the CFL number if explicit time integration is used. The other one belongs to semi-Lagrangian (SL) methods [14,15], which form the advection equation in its “Lagrangian” or flow-following form. For a specific point, the SL methods use the winds (either prescribed or predicted) to trace along the trajectory backward to this point’s original location (called the departure point) on the previous time level. The value of the variable at that upwind point, usually interpolated from the surrounding grid points, is the advected value of the variable at the destination grid point.

The SL methods are designed based on the fixed set of the computational mesh with information propagating along characteristic curves, taking advantage of high spatial resolution in an Eulerian approach and computational efficiency with large time stepping in a Lagrangian method for transport problems [16]. Over the past few decades, the SL methods have been under great development and are widely used in fluid and kinetic simulations, weather forecasting, interface tracking, etc.

The SL method is also a popular choice for developing a global transport model since it allows a large time step without much loss of accuracy. In this paper, we consider the following horizontal transport equation on the spherical surface [17]:

$$\frac{\partial \varphi }{\partial t}+\nabla ·\left(\overrightarrow{V}\varphi \right)=0$$

(1)

where $\varphi$ is the fluid thickness or density of a specific tracer, $\overrightarrow{V}$ is the velocity vector tangent to the sphere for all t ≥ 0. The operator ∇ is assumed to be the horizontal divergence on the sphere.

Equation (1) could be expanded as follows:

$$\frac{\partial \varphi }{\partial t}+\nabla ·\left(\overrightarrow{V}\right)\varphi +\overrightarrow{V}·\nabla \varphi =0$$

(2)

In the case of a divergence-free flow, the Lagrangian form of Equation (2) is as follows:

$$\frac{\partial \varphi }{\partial t}+\overrightarrow{V}·\nabla \varphi =\frac{d\varphi }{dt}=0$$

(3)

where ∇ denotes the horizontal gradient operator on the sphere. The corresponding characteristic equation is as follows:

$$\frac{d\overrightarrow{x}}{dt}=\overrightarrow{V}$$

(4)

where $\overrightarrow{x}$ denotes the position vector of a point on the sphere.

In Section 2, based on the regular icosahedron, a triangle mesh for the generation of uniformly distributed node sets on the sphere surface would be described.

In Section 3, the detailed steps for implementation of the SL method are presented. It contains mainly of the following two procedures: firstly, find the location for the departure point by solving the trajectory equation; second, approximate the value at the departure point from its neighbor nodes by interpolation or fitting methods.

In Section 4, three different methods for the location of departure points are reported, which include the mid-rule method, the explicit classical four-order and five-order Runge–Kutta method. Additionally, comparison results of the accuracy of departure points’ location for these methods are evaluated. Once the departure point is obtained, then four new methods are based on interpolation using (1) linear interpolation; (2) polynomial fitting based on the least square method; (3) global radial basis function stencils (RBFs); (4) local RBFs. A test case, that is, the passive advection of Gauss hill, is utilized to evaluate the accuracy of the location of departure points and the mass conservation of different methods.

In Section 5, a series of numerical experiments that include the solid rotation of a cosine bell, the deformation of twin slotted cylinders, and the moving vortex test cases are utilized to evaluate the performance of the proposed global transport model. Additionally, the summary is given in Section 6.

2. Triangle Mesh Generation on a Sphere

Although the SL method belongs to the mesh-free method essentially, it needs a spatial discretization of the spherical surface. In this paper, a triangle mesh based on a regular icosahedron sphere is utilized [18]. As shown in Figure 1, first, a regular icosahedron was embedded in the sphere. To refine this grid, each geodesic edge is split into two equal halves by adding a new vertex at the midpoint, and each triangular face is consequently split into four sub-triangles. Figure 1 illustrates the mesh refinement procedure. The left panel shows a regular icosahedron mapped on the spherical surface, and the middle panel and right panel show the mesh with refinement levels equal to 1 and 2, respectively. Statistically, a triangle mesh which is refined N times will have 10 × 4^N +2 vertices and 20 × 4^N elements. It is clear that this triangle mesh is semi-structured; there are five triangle elements shared commonly by each of the original twelve vertices on the icosahedron, and each of the reminder grid vertices is shared by six triangles.

3. High-Order Semi-Lagrangian Method

For solving the transport equation, the main idea behind the semi-Lagrangian method is to find the location of the departure point for a grid point and then interpolate to the departure point from the values of its surrounding grid points. As shown in Figure 2, the implementation of the SL method can be divided into the following two procedures:

(1) For each grid point $\overrightarrow{x}$, integrate the characteristic Equation (4) backward for the location of the departure point ${\overrightarrow{x}}_D$. Usually the departure point does not coincide with the original node.

$$\begin{gathered} \frac{d\overrightarrow{x}}{dt}=\overrightarrow{V}(\overrightarrow{x},t) \\ \overrightarrow{x}(t)-\overrightarrow{x}(t-\Delta t)={\displaystyle {\int }_{t-\Delta t}^t\overrightarrow{V}(\overrightarrow{x},t)}dt \\ {\overrightarrow{x}}_D=\overrightarrow{x}(t-\Delta t)=\overrightarrow{x}(t)-{\displaystyle {\int }_{t-\Delta t}^t\overrightarrow{V}(\overrightarrow{x},t)}dt \end{gathered}$$

(5)

(2) Then some interpolation method is utilized to interpolate or fit the value on the departure point from those on the original nodes.

$$\begin{gathered} \varphi (\overrightarrow{x},t)=\varphi (\overrightarrow{x}-{\displaystyle {\int }_{t-\Delta t}^t\overrightarrow{V}}dt,t-\Delta t)=\varphi ({\overrightarrow{x}}_D,t-\Delta t) \\ {\varphi }^t(\overrightarrow{x})={\varphi }^{t-\Delta t}({\overrightarrow{x}}_D)\approx \mathrm{I}[{\varphi }^{t-\Delta t}(\overrightarrow{x})] \end{gathered}$$

(6)

where $\Delta t$ denotes the discrete time step, I denotes a specific interpolation or fitting operator.

In the following subsections, the above two main procedures would be illustrated, respectively.

3.1. Approximation for Departure Point

In this subsection, two main kinds of methods are utilized to integrate the characteristic Equation (4) for the location of the departure point; those are the Ritchie’s midpoint method and high-order explicit Runge–Kutta methods.

3.1.1. Ritchie’s Method

There are many ways of computing the right-hand integral in Equation (5) approximately, but the most commonly used form has been the midpoint rule, which was first proposed by Ritchie [19] and is iterative and interpolating. Giraldo [20] proved its simplification to the midpoint integration rule on the surface of a sphere. Additionally, the midpoint method was used widespread for the approximation of the characteristic equation in the semi-Lagrangian method [21,22,23]. It is implemented mainly as follows:

Integrate the trajectory equation half of a time step backward, yielding the following:

$${\overrightarrow{x}}_M={\overrightarrow{x}}_M-\frac{dt}{2}\overrightarrow{V}\left({\overrightarrow{x}}_M,t-dt/2\right)$$

(7)

where ${\overrightarrow{x}}_M$ enotes the midpoint of the characteristic curve and ${\overrightarrow{x}}_A$ denotes the arrival point. Usually, it was just the mesh vertex. The equation defines a recursive scheme because ${\overrightarrow{x}}_M$ is given implicitly in the relation. Usually, several iterations (4 or 5 iteration loops) are enough to obtain a convergent solution, and then the departure point ${\overrightarrow{x}}_D$ is calculated by the following:

$${\overrightarrow{x}}_D={\overrightarrow{x}}_A-dt\overrightarrow{V}\left({\overrightarrow{x}}_M,t-dt/2\right)$$

(8)

However, we have said nothing about the order of the interpolations, but obviously, the higher the order of the quadrature formulas, the better the trajectory accuracy, but more computational cost as well. The midpoint rule yields second-order accuracy and has been used quite successfully in 2D planar space. On the sphere, the midpoint rule has to be modified such that the new departure points computed by Equation (8) remain on the surface of the sphere. In other words, after each iteration we must apply the following projection:

$${\overrightarrow{x}}_D=\frac{R}{\left|{\overrightarrow{x}}_D\right|}{\overrightarrow{x}}_D$$

(9)

where $R$ denotes the radius of the sphere.

3.1.2. High-Order Runge–Kutta Methods

Although the Ritchie’s midpoint method could be implemented easily, it gives only second-order accuracy. Even higher order methods, such as the classical explicit four-order Runge–Kutta method (RK4), could also be used to solve the characteristic equation; it yields the following:

$$\begin{gathered} d{\overrightarrow{x}}_{k1}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A,t\right) \\ d{\overrightarrow{x}}_{k2}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+\frac{d{\overrightarrow{x}}_{k1}}{2},t-\frac{dt}{2}\right) \\ d{\overrightarrow{x}}_{k3}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+\frac{d{\overrightarrow{x}}_{k2}}{2},t-\frac{dt}{2}\right) \\ d{\overrightarrow{x}}_{k4}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+d{\overrightarrow{x}}_{k3},t-dt\right) \\ {\overrightarrow{x}}_D={\overrightarrow{x}}_A+\frac{d{\overrightarrow{x}}_{k1}+2d{\overrightarrow{x}}_{k2}+2d{\overrightarrow{x}}_{k3}+d{\overrightarrow{x}}_{k4}}{6} \end{gathered}$$

(10)

Additionally, classical Butcher’s Runge–Kutta method of order five (RK5) as follows [24]:

$$\begin{gathered} d{\overrightarrow{x}}_{k1}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A,t\right) \\ d{\overrightarrow{x}}_{k2}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+\frac{d{\overrightarrow{x}}_{k1}}{4},t-\frac{dt}{4}\right) \\ d{\overrightarrow{x}}_{k3}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+\frac{d{\overrightarrow{x}}_{k1}+d{\overrightarrow{x}}_{k2}}{8},t-\frac{dt}{4}\right) \\ d{\overrightarrow{x}}_{k4}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A-\frac{d{\overrightarrow{x}}_{k2}-2d{\overrightarrow{x}}_{k3}}{2},t-\frac{dt}{2}\right) \\ d{\overrightarrow{x}}_{k5}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+\frac{3d{\overrightarrow{x}}_{k1}+9d{\overrightarrow{x}}_{k4}}{16},t-\frac{3dt}{4}\right) \\ d{\overrightarrow{x}}_{k6}=-dt\overrightarrow{V}\left({\overrightarrow{x}}_A+\frac{-3d{\overrightarrow{x}}_{k1}+2d{\overrightarrow{x}}_{k2}+12d{\overrightarrow{x}}_{k3}-12d{\overrightarrow{x}}_{k4}+8d{\overrightarrow{x}}_{k5}}{7},t-dt\right) \\ {\overrightarrow{x}}_D={\overrightarrow{x}}_A+(7d{\overrightarrow{x}}_{k1}+32d{\overrightarrow{x}}_{k3}+12d{\overrightarrow{x}}_{k4}+32d{\overrightarrow{x}}_{k5}+7d{\overrightarrow{x}}_{k6})/90 \end{gathered}$$

(11)

When the departure point ${\overrightarrow{x}}_D$ is obtained, then it should be mapped onto the sphere surface by Equation (9). In theory, the high-order Runge–Kutta methods could generate more accurate locations of departure points than the midpoint method. In Section 4, a numerical experiment would be used to validate the results.

3.2. Interpolation

Once the departure point is calculated, then the value at the departure point could be obtained by interpolation or fitting methods. In this section, four methods would be utilized, which include linear interpolation, polynomial fitting based on the least square method, global radial basis function stencil (RBFs), and local RBFs interpolation; the main ideas of these methods would be illustrated, respectively.

3.2.1. Linear Interpolation

The main idea for linear interpolation is as follows. First find the triangle element that the arrival point traced back into, then interpolate its value by those on the triangle’s vertex [25]. For instance, in Figure 3, assuming the red point ${\overrightarrow{x}}_D$ denotes the departure point for a specific vertex and it is located in the triangle element $P_0{P}_1{P}_2$, then its value could be approximated by the values of the three vertexes of the triangle as follows:

$$\varphi \left({\overrightarrow{x}}_D\right)=\frac{{\alpha }_0\varphi \left(P_0\right)+{\alpha }_1\varphi \left(P_1\right)+{\alpha }_2\varphi \left(P_2\right)}{{\alpha }_0+{\alpha }_1+{\alpha }_2}$$

(12)

where ${\alpha }_i,i=0,1,2$ denotes the spherical area of the sub-triangle formed by the departure point and the edges of the triangle opposite to vertex $P_i$.

This method could be implemented with little effort, and it can be seen that the interpolated value is just an area-weighted average of the values on three vertexes of the triangle.

3.2.2. Local Polynomial Fitting by Least-Square Method (Poly-LSQ)

The linear interpolation method above-mentioned has two-order accuracy. In order to attain more accurate interpolation results, high-order polynomial fitting was utilized. It needs a local grid stencil in 3D Cartesian spaces, and then the grid points in this local stencil would be mapped onto a local two-dimensional (2D) coordinate system by a general stereography projection [26,27], and then a bivariate polynomial was utilized to fit the values at grid points on the stencil. In order to obtain an approximate value at the departure point, first, find the nearest vertex for each departure point, and then a polynomial least-squares fit is calculated using values at the vertex and its nearest first-layer neighbor points.

The transformation between the spherical latitude–longitude coordinates and two-dimensional general stereographic coordinates (GSTC) is given by the following:

$$\begin{gathered} x=mR\cos \left(\theta \right)\sin \left(\lambda -{\lambda }_0\right) \\ y=mR\left[\sin \left(\theta \right)\cos \left({\theta }_0\right)-\cos \left(\theta \right)\sin \left({\theta }_0\right)\cos \left(\lambda -{\lambda }_0\right)\right] \\ m=\frac{2}{1+\sin \left(\theta \right)\sin \left({\theta }_0\right)+\cos \left(\theta \right)\cos \left({\theta }_0\right)\cos \left(\lambda -{\lambda }_0\right)} \end{gathered}$$

(13)

where m is the map factor and R is the radius of the earth, ($\theta ,\lambda$) and (x,y) denote the latitude–longitude coordinates and GSTC respectively. The center of the GSTC is at the tangent point where (x,y) = (0,0). In the transformed system, the positive x-axis is directed toward the east of the origin along the latitudinal circle of $\theta ={\theta }_0$ and the y-axis is defined as positive poleward of the origin along the meridian $\lambda ={\lambda }_0$.

In Figure 4a, an icosahedron sphere triangle mesh with refinement level 1 is shown. The north pole and its five nearest neighbor points shown in red dots are projected on the 2D tangent plane, and the blue lines denote projection lines. As shown in Figure 4b, for a specific vertex, its departure point is denoted by the red point ${\overrightarrow{x}}_D$. Firstly, find the vertex that is closest to this departure point, and then the value at ${\overrightarrow{x}}_D$ could be approximated by the values at the vertex and its first-order neighbor points, which leads to a full bivariate quadratic polynomial as follows:

$$\widehat{\varphi }\left(x,y\right)=a_0+a_{1}x+a_{2}y+a_3{x}^2+a_{4}xy+a_5{y}^2$$

(14)

In order to minimize the interpolation error, the fitted polynomial should be constrained by the following condition:

$$\mathrm{Min}{\displaystyle \sum }_{i=1}^{NBP}{\left[\widehat{\varphi }\left(p_i\right)-\varphi \left(p_i\right)\right]}^2$$

(15)

where the NBP = 6 or 7, is the number of first-order neighbor points for a specific vertex, including the vertex itself. $p_i$ denotes the projected point on the local 2D tangent plane.

For Equation (15), there are seven or six values on each stencil but only six polynomial terms. In order to determine the coefficient $a_i,i=0,1,\cdots ,5$, a Least-squares approach is needed because the system of equations is over-constrained. For a specific vertex, the stencil geometry is expressed in a local coordinate system with the vertex that is nearest the departure point as the origin so that the value at the vertex is equal to the constant coefficient $a_0$. Once the coefficients are obtained, an approximate value at the departure point could be obtained by inserting its mapped coordinates into the quadratic polynomial (14).

3.2.3. Global RBFs

In recent two decades, the mesh-free methods have achieved rapid development, among which the radial basis function (RBF) method is one of the most popular numerical methods, and it has been vastly applied to solve various kinds of physical problems with theoretical proofs on solvability and convergence presented in [28]. Increasingly, the global RBF collocation method is applied to some problems in the sphere [29,30].

Given a set of nodes $X=\left\{{\overrightarrow{x}}_k\right\},k=1,\cdots ,N$ on the sphere surface, the standard global RBF (GRBF) interpolant $\widehat{\varphi }\left(\overrightarrow{x}\right)$ to the data has the following form:

$$\widehat{\varphi }\left(\overrightarrow{x}\right)={\displaystyle \sum }_{k=1}^N{\omega }_{k}B\left(\Vert \overrightarrow{x}-{\overrightarrow{x}}_k\Vert \right)$$

(16)

where B denotes the kernel function of RBF, the notation $\overrightarrow{x}=\left(x,y,z\right)$ and ${\overrightarrow{x}}_k=\left(x_k,y_k,z_k\right)$, the Equation (16) is the linear combination of the kernel function of RBF. In this paper, we will limit our attention to the kernel functions B(r) listed in

Table 1. Definitions of some infinitely differentiable radial functions, the shape parameter c controls their flatness.

Name of RBF	Abbreviation	Definition
Gaussian	GA	$e^{-{\left(cr\right)}^2}$
Multiquadric	MQ	$\sqrt{1+{\left(cr\right)}^2}$
Inverse multiquadric	IMQ	$\frac{1}{\sqrt{1+{\left(cr\right)}^2}}$
Inverse quadratic	IQ	$\frac{1}{1+{\left(cr\right)}^2}$

. All these four kernel functions are global supported and strictly positive definite. Other types of kernel functions that are local-compactly supported could be found in [31]. Note that there is a parameter c in each kernel function, it is known as the shape parameter.

If the interpolant $\widehat{\varphi }\left(\overrightarrow{x}\right)$ can be utilized to approximate the true solution, then the expansion coefficients are determine ${\omega }_k,k=1,\cdots ,N$ by enforcing $\widehat{\varphi }\left({\overrightarrow{x}}_k\right)=\varphi \left({\overrightarrow{x}}_k\right), k=1,\cdots ,N$, which can be expressed by the following linear system:

$$\left[\begin{array}{ccccc}B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_3\Vert )& & B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_N\Vert )\\ B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_3\Vert )& \cdots & B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_N\Vert )\\ B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_3\Vert )& & B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_N\Vert )\\ & ⋮& & \ddots & ⋮\\ B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_3\Vert )& \cdots & B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_N\Vert )\end{array}\right]\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_3\\ ⋮\\ {\omega }_N\end{array}\right]=\left[\begin{array}{c}\varphi \left({\overrightarrow{x}}_1\right)\\ \varphi \left({\overrightarrow{x}}_2\right)\\ \varphi \left({\overrightarrow{x}}_3\right)\\ ⋮\\ \varphi \left({\overrightarrow{x}}_N\right)\end{array}\right]$$

(17)

If the Gauss type or IMQ type basis function was used for RBF interpolation, the coefficient matrix in Equation (17) is positive-definite and invertible, and then the linear equations admit a unique solution. If the MQ type basis equation was utilized, the resulting matrix for interpolation may not be positive-definite. Linear-independent polynomial functions can be added to make the matrix positive-definite and invertible. In this case, the RBF interpolant is expressed as follows:

$$\widehat{\varphi }\left(\overrightarrow{x}\right)={\displaystyle \sum }_{k=1}^N{\omega }_{k}B\left(\Vert \overrightarrow{x}-{\overrightarrow{x}}_k\Vert \right)+{\gamma }_0+{\gamma }_{1}x+{\gamma }_{2}y+{\gamma }_{3}z$$

(18)

with the following constraints:

$${\displaystyle \sum }_{k=1}^N{\omega }_k={\displaystyle \sum }_{k=1}^N{\omega }_k{x}_k={\displaystyle \sum }_{k=1}^N{\omega }_k{y}_k={\displaystyle \sum }_{k=1}^N{\omega }_k{z}_k=0$$

(19)

The resulting linear system:

$$\left[\begin{array}{ccccccccc}B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_3\Vert )& & B(\Vert {\overrightarrow{x}}_1-{\overrightarrow{x}}_N\Vert )& 1& x_1& y_1& z_1\\ B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_3\Vert )& \cdots & B(\Vert {\overrightarrow{x}}_2-{\overrightarrow{x}}_N\Vert )& 1& x_2& y_2& z_2\\ B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_3-{\overrightarrow{x}}_3\Vert )& & B(\Vert x_3-{\overrightarrow{x}}_N\Vert )& 1& x_3& y_3& z_3\\ & ⋮& & \ddots & & ⋮\\ B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_1\Vert )& B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_2\Vert )& B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_3\Vert )& & B(\Vert {\overrightarrow{x}}_N-{\overrightarrow{x}}_N\Vert )& 1& x_N& y_N& z_N\\ 1& 1& 1& & 1& 0& 0& 0& 0\\ x_1& x_2& x_3& \cdots & x_N& 0& 0& 0& 0\\ y_1& y_2& y_3& & y_N& 0& 0& 0& 0\\ z_1& z_2& z_3& & z_N& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\omega }_1\\ {\omega }_2\\ {\omega }_3\\ ⋮\\ {\omega }_N\\ {\gamma }_0\\ {\gamma }_1\\ {\gamma }_2\\ {\gamma }_3\end{array}\right]=\left[\begin{array}{c}\varphi \left({\overrightarrow{x}}_1\right)\\ \varphi \left({\overrightarrow{x}}_2\right)\\ \varphi \left({\overrightarrow{x}}_3\right)\\ ⋮\\ \varphi \left({\overrightarrow{x}}_N\right)\\ 0\\ 0\\ 0\\ 0\end{array}\right]$$

(20)

The choice of basis function and shape parameter has a significant impact on the accuracy of RBF interpolation. Decreasing c or increasing the number of data points has a severe effect on the stability of the linear system (17) or (20). For a fixed shape parameter c, the condition number of the matrix in the linear system grows exponentially as the number of data points is increased. For a fixed number of data points, as the shape parameter c→0, the basis function becomes increasingly flat; meanwhile, the condition number of the linear system grows exponentially [32]. Figure 5 shows the order of magnitude for the condition number Log₁₀(cond(A)) as a function of shape parameter c under three different mesh refinement levels, where A denotes the coefficient matrix in Equations (17) or (20).

It can be seen that the condition number of the coefficient matrix becomes very large when the value of the shape parameter is close to zero for all four types of basis functions. In general, the condition number of the coefficient matrix grows exponentially as the mesh refinement level increases; the numerical computation tends to become unstable. Meanwhile, the condition number of the coefficient matrix decreased drastically as the value of the shape parameter increased. This is particularly significant because the highest accuracy is often found at some “small” shape parameter, which may result in a large condition number of interpolation matrices and even cause numerical instability. This conflict between accuracy and stability is sometimes referred to as the “trade-off principle”. However, there is no exact analytical formula to obtain the optimal value of the shape parameter. In most cases, it mainly relies on statistical methods [33,34,35,36], such as leave one out cross validation (LOOCV), maximum likelihood estimator (MLE), etc.

One issue with global RBF interpolation is that the coefficient matrix in the linear system is dense, thereby leading to considerable computational costs to solve it directly. Another issue is that the matrices can become ill-conditioned. There are a series of methods to tackle the problem, such as the Contour-Padé, RBF-RA, and RBF-QR algorithms [37]; however, they are not pursued and developed here.

3.2.4. Local RBFs

In order to relieve the problems of ill-conditioning and shape parameter sensitivity of the global RBF, the local RBF method has evolved, and the main idea of the local RBF (LRBF) method applied to the sphere is detailedly presented by [38]. Subsequently, some researchers use LRBF to deal with different complicated problems [39,40,41,42,43,44], etc. Similar to GRBF, LRBF is meshless in character and it needs no polygonization. In contrast to the GRBF, the LRBF is made on the overlapping sub-domains. For approximating the objective function at a given node, only its influence nodes are used. The overlapping sub-domains significantly reduce the size of the linear system at the cost of solving lots of small matrices. In addition, the local RBF method is comparatively less sensitive to the shape parameter. In many cases, the local RBF method can be as accurate as the global RBF method.

For a specific vertex ${\overrightarrow{x}}_i$, assuming the node set $\left\{{\overrightarrow{x}}_{ik}\right\},k=1,\cdots ,n$ is its n nearest neighbor points including the ${\overrightarrow{x}}_i$ itself ($n\ll N$). The local RBF interpolant ${\widehat{\varphi }}_i\left(\overrightarrow{x}\right)$ for the vertex ${\overrightarrow{x}}_i$ takes the form:

$${\widehat{\varphi }}_i\left(\overrightarrow{x}\right)={\displaystyle \sum }_{k=1}^n{\omega }_{ik}B(\Vert \overrightarrow{x}-{\overrightarrow{x}}_{ik}\Vert )$$

(21)

In the local region influenced by ${\overrightarrow{x}}_i$, in order to approximate the true solution by the above Equation (21), we could obtain a linear system written in matrix form as follows:

$$\left[\begin{array}{ccccc}B(\Vert {\overrightarrow{x}}_{i1}-{\overrightarrow{x}}_{i1}\Vert )& B(\Vert {\overrightarrow{x}}_{i1}-{\overrightarrow{x}}_{i2}\Vert )& B(\Vert {\overrightarrow{x}}_{i1}-{\overrightarrow{x}}_{i3}\Vert )& & B(\Vert {\overrightarrow{x}}_{i1}-{\overrightarrow{x}}_{in}\Vert )\\ B(\Vert {\overrightarrow{x}}_{i2}-{\overrightarrow{x}}_{i1}\Vert )& B(\Vert {\overrightarrow{x}}_{i2}-{\overrightarrow{x}}_{i2}\Vert )& B(\Vert {\overrightarrow{x}}_{i2}-{\overrightarrow{x}}_{i3}\Vert )& \cdots & B(\Vert {\overrightarrow{x}}_{i2}-{\overrightarrow{x}}_{in}\Vert )\\ B(\Vert {\overrightarrow{x}}_{i3}-{\overrightarrow{x}}_{i1}\Vert )& B(\Vert {\overrightarrow{x}}_{i3}-{\overrightarrow{x}}_{i2}\Vert )& B(\Vert {\overrightarrow{x}}_{i3}-{\overrightarrow{x}}_{i3}\Vert )& & B(\Vert {\overrightarrow{x}}_{i3}-{\overrightarrow{x}}_{in}\Vert )\\ & ⋮& & \ddots & ⋮\\ B(\Vert {\overrightarrow{x}}_{in}\cdots {\overrightarrow{x}}_{i1}\Vert )& B(\Vert {\overrightarrow{x}}_{in}-{\overrightarrow{x}}_{i2}\Vert )& B(\Vert {\overrightarrow{x}}_{in}-{\overrightarrow{x}}_{i3}\Vert )& \cdots & B(\Vert {\overrightarrow{x}}_{in}-{\overrightarrow{x}}_{in}\Vert )\end{array}\right]\left[\begin{array}{c}{\omega }_{i1}\\ {\omega }_{i2}\\ {\omega }_{i3}\\ ⋮\\ {\omega }_{in}\end{array}\right]=\left[\begin{array}{c}\varphi \left({\overrightarrow{x}}_{i1}\right)\\ \varphi \left({\overrightarrow{x}}_{i2}\right)\\ \varphi \left({\overrightarrow{x}}_{i3}\right)\\ ⋮\\ \varphi \left({\overrightarrow{x}}_{in}\right)\end{array}\right]$$

(22)

In the above equation $i1,i2,\cdots ,in$, is a vector associated with the center point ${\overrightarrow{x}}_i$ that contains the index of itself and the indices of its $n-1$ nearest neighboring vertexes, ${\omega }_{ik},k=1,\cdots ,n$ is a vector of expansion coefficients.

For a specific departure point, firstly, find the nearest vertex close to it, and then search for the n nearest neighboring vertexes (including the vertex itself) around this vertex. The procedure could be implemented efficiently by the KDtree search method, at last, a corresponding linear system (20) is formed, when the expansion coefficients are obtained, then Equation (19) could be utilized to interpolate the value at the departure point.

The local RBF method suffers from the same problem as the global RBF, yet less severe. Another issue is that there is no exact theory for the optimal size of the local stencil for a specific vertex. Several values are commonly used, such as 31, 50, 74, and 101 are very popular, when the shape parameter of kernel function is fixed, the condition number of the matrix in the linear system grows quickly as the size of local stencil is increased. Therefore, the local stencil size is chosen as 31 in the following numerical test cases if the local RBF method is involved. In Figure 6, a local stencil of size 31 on the icosahedron spherical triangle mesh with refinement level 4 was shown, the center vertex is shown in blue asterisk, and its nearest neighbor points are shown in the red dot.

For the above-mentioned four interpolation methods, the linear interpolation method and the Poly-LSQ method could be implemented without much effort; however, we need to take extra care in selecting the value of the shape parameter in the basis function if the RBF methods are utilized. If no specific explanations, for all the numerical experiments in Section 4 and Section 5, the density of tracer or height fields are placed at the nodes of triangle elements.

4. Trajectory Accuracy and Mass Conservation Evaluation

4.1. Trajectory Accuracy

In a similar manner described in McGregor [45], the accuracy of the trajectories is defined by the following normalized L₂ norm:

$${\widehat{\varphi }}_i\left(\overrightarrow{x}\right)=\sqrt{\frac{{∯}_S^{}{\left({\overrightarrow{x}}_D-{\overrightarrow{x}}_D^{exact}\right)}^{2}ds}{{∯}_S^{}{\left({\overrightarrow{x}}_D-{\overrightarrow{x}}_A\right)}^{2}ds}}$$

(23)

where ${\overrightarrow{x}}_A$ denotes the arrival point, and ${\overrightarrow{x}}_D$, ${\overrightarrow{x}}_D^{exact}$ denote the calculated departure point and the exact one respectively.${∯}_S^{}ds$ denotes the surface integral on the sphere.

This test case is advection of a Gaussian hill around the sphere [46]. The Gaussian hill is prescribed as follows:

$${\widehat{\varphi }}_i=h_0{e}^{-b{r}^2}$$

(24)

where $r=\Vert \overrightarrow{x}-{\overrightarrow{x}}_0\Vert$ denotes the Euclidean distance between the two points $\overrightarrow{x}$, ${\overrightarrow{x}}_0$, $h_0=0.95$ and $b=5$ respectively. The 3D Cartesian coordinates $\overrightarrow{x}=\left(x,y,z\right)$ and the latitude–longitude coordinates $\left(\theta ,\lambda \right)$ on spherical surface are related through

$$\left(x,y,z\right)=\left(R\cos \theta \cos \lambda ,R\cos \theta \sin \lambda ,R\sin \theta \right)$$

(25)

where R is the radius of the sphere. The coordinates for the center of the Gaussian distribution ${\overrightarrow{x}}_0$ are computed by inserting the tracer center $\left({\theta }_0,{\lambda }_0\right)$ into Equation (25) and evaluating the right-hand side.

The initial wind fields are as follows:

$$\begin{gathered} u=u_0\left(\cos \left(\theta \right)\cos \left(\alpha \right)+\sin \left(\theta \right)\cos \left(\lambda -\frac{\pi }{2}\right)\sin \left(\alpha \right)\right) \\ V=-u_0\sin \left(\lambda -\frac{\pi }{2}\right)\sin \left(\alpha \right) \end{gathered}$$

(26)

where the flow constant $u_0=2\pi R/\left(64 \mathrm{hours}\right)$, $\alpha$ is the rotation angular between the wind flow direction and the equatorial plane, the prescribed velocity field is that of solid-body rotation, such that the tracer moves around the earth once with a rotational period of 64 h.

The tracer center was set as $\left({\theta }_0,{\lambda }_0\right)=\left(0,0\right)$, the rotation angle of flow-vectors $\alpha =0$. Then the mid-rule method and high-order Runge–Kutta methods were utilized to solve the trajectory equation for the departure points’ location respectively. The exact location of the test case is given, and comparison of calculated departure points’ location for the midpoint rule method, RK4 and RK5 methods under different time steps are given in

Table 2. The comparison for departure points’ location error between exact solution and these of three different methods.

Refinement Level	Node Number	Time-Step	Midpoint Rule	RK4	RK5
Level 3	642	2 h	0.0012	$5.4257\times {10}^{-6}$	$2.3382\times {10}^{-8}$
		4 h	0.0049	$8.6429\times {10}^{-5}$	$8.1214\times {10}^{-7}$
		8 h	0.0205	1.4 × 10−3	$3.2846\times {10}^{-5}$
Level 4	2562	2 h	0.0012	$5.4263\times {10}^{-6}$	$2.3368\times {10}^{-8}$
		4 h	0.0049	$8.6440\times {10}^{-5}$	$8.1174\times {10}^{-7}$
		8 h	0.0205	1.4 × 10−3	$3.2840\times {10}^{-5}$
Level 5	10,242	2 h	0.0012	$5.4264\times {10}^{-6}$	$2.3364\times {10}^{-8}$
		4 h	0.0049	$8.6442\times {10}^{-5}$	8.1164 $\times {10}^{-7}$
		8 h	0.0205	1.4 × 10−3	$3.2838\times {10}^{-5}$

.

It can be seen from Table 2 that when the number of nodes is fixed, as the time step decreased to half of the precedent one, the error of departure points’ location decreased to about 1/4 of the last one; therefore, the midpoint rule method has two-order accuracy. Similarly, it can be seen that the RK4 method achieved four-order accuracy while the RK5 method has at least a five-order convergence rate. Obviously, the RK5 method owns the most accurate results. Therefore, the RK5 method is utilized for computing the location of the departure point in all the following numerical experiments. We also note that an increase in the spatial resolution has little impact on the location errors of departure points when the time step is fixed.

4.2. Mass Conservations

In order to check the conservation of the presented scheme, the mass conservations is commonly used, that is, the total mass should be kept as a constant as follows:

$$\begin{gathered} \mathrm{Mass}: I_1\left(t\right)={\iint }_{S^2}^{}\varphi \left(\overrightarrow{x},t\right)ds \\ \frac{d{I}_1\left(t\right)}{dt}=0 \end{gathered}$$

(27)

where ${\iint }_{S^2}^{}\varphi \left(\overrightarrow{x},t\right)ds$ denotes the integral of $\varphi$ on the spherical surface S². The numerical relative error for $I_1\left(t\right)$ is reported using the following:

$$(I_1\left(t\right)-I_1\left(0\right))/I_1\left(t\right)$$

(28)

In order to facilitate comparison between different methods for the following test case simulations, the L₁, L₂ and $L_{\infty }$ norms of mass error are specified as follows [47]:

$$\begin{array}{c}{L}_1=\frac{{\displaystyle {∯}_S\left|{\varphi }_n(j)-{\varphi }_r(j)\right|ds}}{{\displaystyle {∯}_S\left|{\varphi }_r(j)\right|ds}}\approx \frac{{\displaystyle \sum \left|{\varphi }_n(j)-{\varphi }_r(j)\right|A_j}}{{\displaystyle \sum \left|{\varphi }_r(j)\right|A_j}}\\ L_2=\sqrt{\frac{{\displaystyle {∯}_S{({\varphi }_n(j)-{\varphi }_r(j))}^{2}ds}}{{\displaystyle {∯}_S{\varphi }_r{(j)}^{2}ds}}}\approx \sqrt{\frac{{\displaystyle \sum {({\varphi }_n(j)-{\varphi }_r(j))}^2{A}_j}}{{\displaystyle \sum {\varphi }_r{(j)}^2{A}_j}}}\\ L_{\infty }=\frac{\underset{j}{\mathrm{Max}}\left|{\varphi }_n(j)-{\varphi }_r(j)\right|}{\underset{j}{\mathrm{Max}}\left|{\varphi }_r(j)\right|}\end{array}$$

(29)

The function ${\varphi }_n\left(j\right)$ denotes the numerical solution defined at the ${\overrightarrow{x}}_j$ positions on the numerical mesh. The index j represents the specific cell on the pentagon–hexagon mesh dual to the icosahedron spherical triangle mesh, and $A_j$ denotes the spherical area associated with the Voronoi cell controlled by vertex ${\overrightarrow{x}}_j$. The function ${\varphi }_r\left(j\right)$ is the reference solution that has been calculated at or interpolated to the same ${\overrightarrow{x}}_j$ positions. The reference solution represents a high-resolution solution that is sufficiently accurate for the computation of the error norms if an analytic solution is not available. The function $\underset{j}{\max }\left|\varphi \right|$ finds the global maximum of the $\left|\varphi \right|$ evaluated at the ${\overrightarrow{x}}_j$ positions of the numerical mesh.

In order to test the interpolation accuracy of four different methods and check the mass conservation, the above-mentioned Gauss hill advection is utilized as the benchmark again, and the mesh refinement level is equal to 4. The relative error norm of the tracer calculated by four different methods for interpolation is shown in Figure 7. It can be seen that when the time step is too large, the error norm increases. Obviously, the linear interpolation method gives results with relatively large errors in Figure 7a. If the time step is more than 4 h, the method would generate large errors. In comparison, the polynomial fitting by the least-square method could generate more accurate results, as shown in Figure 7b. With an appropriate value of shape parameter, the global RBF and local RBF could also give more accurate results in Figure 7c,d, for the global RBF method, the Gauss basis function is utilized, and the shape parameter c =1.5; for the local RBF method, the Gauss basis function is utilized, the size of the local stencil is 31, and the shape parameter c = 1.5. It is obvious the last three methods generate results that are in the same order of magnitude. We also note that the global RBF and local RBF methods give much like results.

5. Numerical Experiments

In this section, the semi-Lagrangian method will be utilized to model the evolution process of the advection of a passive tracer; three test cases are considered. They are referred to as the solid body rotation, the shear deformation of twin slotted cylinders, and the moving vortex. In the first case, the mass conservation properties of our scheme would be numerically evaluated. Additionally, the modeling results would be presented in the other two cases.

Test case 1: solid body rotation

This test case describes the advection of a cosine bell once around the sphere [48,49,50]. The initial shape of the bell is formulated as follows:

$$\varphi =\{\begin{array}{c}\frac{h_0}{2}\left(1+\cos \frac{\pi r}{R_b}\right), if r<R_b \\ 0, if r\ge R_b \end{array}$$

(30)

where $h_0=1000 \mathrm{m}$ denotes the initial bell maximize height. $R_b=R/3$ is the radius of the cosine bell,$r$ denotes the great circle distance between $\left(\lambda ,\theta \right)$ and the center of the cosine bell $\left({\lambda }_0,{\theta }_0\right)=\left(0,0\right)$, given by the following:

$$r=R{\cos }^{-1}(\sin \left(\theta \right)\sin \left({\theta }_0\right)+\cos \left(\theta \right)\cos ({\theta }_0)\cos (\lambda -{\lambda }_0))$$

(31)

The numerical results of the SL scheme are implemented on Tri_Icosa_L5 mesh with the time step equal to 12 hours (h). The initial wind fields are the same as Equation (26), while the parameter $u_0=2\pi R/\left(12 \mathrm{days}\right)$, days = 24 h.

To check the influence of the wind flow direction $\alpha$ upon the icosahedron sphere gird, this test is conducted with $\alpha =0$ (eastward along the equator) and $\alpha =\pi /2$ (across the north and south poles) respectively. In Figure 8a, it shows L₁, L₂, and $L_{\infty }$ norm error for the cosine bell test case by the global RBF method with Gauss kernel function; the shape parameter is 16, and the wind direction angle $\alpha =0$; Figure 8b shows the same content but the wind direction angle $\alpha =\pi /2$. In Figure 8c, it shows L₁, L₂ and $L_{\infty }$ norm error for the cosine bell test case by the local RBF method with Gauss kernel function; the shape parameter is 20, the wind direction angle $\alpha =0$; Figure 8d shows the same content but the wind direction angle $\alpha =\pi /2$. It can be seen that the history of normalized errors of mass by the local RBF method with different flow angles is almost the same; four panels show the L₁, L₂ error norms are in the same order of magnitude, which means that the wind direction angle has little effect on the total mass. This accords with the facts. We also note that the normalized errors by the global RBF method are slightly more accurate than those by the local RBF method.

In order to evaluate the efficiency of the proposed SL method, another three methods for the approximation of the transport equation on the sphere were utilized for comparison, those are, the RBF method, pseudo-spectral method (PS), and two-order finite difference (FD2) method [51]. The SL method and RBF method are implemented on an icosahedron sphere triangle mesh and the Gauss basis function is used. The PS method and FD2 method are on latitude–longitude grids. The explicit four-order classical Runge–Kutta method was used in all methods for time integration except the SL method.

Table 3. Comparison results between the proposed SL method, the PS method, FD2 method, and RBF method under different space resolutions.

Method	Basis Function and c	Node Number	Time-Step	${\mathit{L}}_{\mathit{2}}$ Norm	${\mathit{L}}_{\mathit{\infty }}$ Norm
SL method (Global RBFs)	GA, c = 1.5	642	4 h	0.0443	0.0371
	GA, c = 6.0	2562	2 h	0.0046	0.0030
	GA, c = 16.0	10,242	1 h	0.0011	0.0011
PS method	-	$32\times 16$	9/40 h	0.2243	0.1874
	-	$64\times 32$	9/80 h	0.0558	0.0460
	-	$128\times 64$	9/160 h	0.0556	0.0483
FD2 method	-	$32\times 16$	9/40 h	1.0790	0.6150
	-	$64\times 32$	9/80 h	1.0211	0.6056
	-	$128\times 64$	9/160 h	0.7574	0.5983
RBF method	GA, c = 3.2	642	1/20 h	0.0987	0.0471
	GA, c = 6.0	2562	1/80 h	0.0242	0.0061
	GA, c = 16.0	10,242	1/320 h	0.0043	0.0021

shows the comparison results of mass error between the proposed SL method, the PS method, the FD2 method, and the RBF method under different spatial resolutions.

From the $L_2$ and $L_{\infty }$ norm of mass error shown in Table 3, it can be seen that the results obtained by the SL method based on global RBF interpolation own the most accuracy, and then the RBF method follows; at last, the PS method and FD2 method give least accurate results. It could also be noted that the SL method admits a relatively larger time step than other methods.

Test case 2: Deformation of twin slotted cylinders

In this test case, the tracer’s shape is not smooth, and it is more challenging and designed to assess the shape-preserving properties of the proposed method, the tracer has a distribution shape of twin slotted cylinders defined as follows [52,53,54]:

$${\varphi (\lambda ,\theta ,t)|}_{t=0}=\{\begin{array}{lr}1,& if r_i\le R_b and |\lambda -{\lambda }_i|\ge \frac{R_b}{6} for i=0,1\\ 1,& if r_0\le R_b and |\lambda -{\lambda }_0|<\frac{R_b}{6} and \theta -{\theta }_0<-\frac{5R_b}{12}\\ 1,& if r_1\le R_b and |\lambda -{\lambda }_1|<\frac{R_b}{6} and \theta -{\theta }_1>\frac{5R_b}{12}\\ 0,& otherwise\end{array}$$

(32)

where $R_b$ denotes the radius of the cylinder,$R_b=0.5$. $r_i=r_i\left(\lambda ,\theta \right)$ is the great-circle distance from the center of the cylinder located at $\left({\lambda }_i,{\theta }_i\right),i=0,1$. The initial positions of the centers of two cylinders are at $\left({\lambda }_0,{\theta }_0\right)=\left(-\frac{\pi }{6},0\right)$ and $\left({\lambda }_1,{\theta }_1\right)=\left(\frac{\pi }{6},0\right)$ respectively. The slots are oriented in opposite directions for the two cylinders so that they are symmetric with respect to the flow.

The wind field is non-divergent but highly deformational. The initial distributions are deformed into thin filaments halfway through the simulation while they are being transported along the zonal direction by the solid-body component of the flow. Note that an exact solution for this test is only available at the final time t = T, and it is just the same as the initial distributions. The time-dependent non-divergent wind field is defined as follows:

$$\begin{gathered} u\left(\lambda ,\theta ,t\right)=\kappa {\sin }^2\left(\widehat{\lambda }\right)\sin \left(2\theta \right)\cos \left(\pi t/T\right)+2\pi \cos \left(\theta \right)/T \\ V\left(\lambda ,\theta ,t\right)=\kappa \sin \left(2\widehat{\lambda }\right)\cos \left(\theta \right)\cos \left(\pi t/T\right) \end{gathered}$$

(33)

where $\widehat{\lambda }=\lambda -\frac{2\pi t}{T},\kappa =2.0$, and the final time T = 5 in non-dimensional time units.

An icosahedron sphere triangle mesh with refinement level 5 is used. The initial distribution of the tracer is shown in Figure 9a, and Figure 9b shows the tracer’s restoration to the original shape at the final time by the Poly-LSQ method with $\mathsf{\Delta }t=0.1$. Figure 9c shows the tracer’s shape at the final time by the global RBF interpolation based on the Gauss kernel function with shape parameter c = 18.0, time step $\mathsf{\Delta }t=0.1$; Figure 9d shows the tracer’s shape at the final time by the local RBF method based on the Gauss kernel function with shape parameter c = 19.0, time step $\mathsf{\Delta }t=0.1$.

It can be seen that the above three interpolation methods could reproduce this complicated deformational flow correctly, while the shape edges of the twin cylinders are smeared, and small overshoots or undershoots are observed near the sharp edges of the twin cylinders. For the Poly-LSQ method and the local RBF method, the oscillations are more obvious in the upwind direction near the edge of the twin cylinders than in the global RBF method. In general, the global RBF method performs the best.

Test case 3: moving vortex

In order to further test the advection scheme representing tracer evolution in a non-divergent but deforming flow, this test case is about the evolution of idealized cyclogenesis on the sphere, which was first described by Nair et al. [55], and then was used as a benchmark for transport processes on a spherical surface. The initial condition was spun up by an angular velocity field, resulting in two diametrically opposed vortices. Let $\left(\acute{\lambda },\acute{\theta }\right)$ be the rotated coordinate system with the north pole at $\left({\lambda }_p,{\theta }_p\right)$ with respect to the regular spherical coordinates $\left(\lambda ,\theta \right)$. It follows:

$$\begin{gathered} \acute{\lambda }=\tan 2^{-1}\left(\sin \left(\lambda -{\lambda }_p\right),\sin ({\theta }_p)\cos \left(\lambda -{\lambda }_p\right)-\cos ({\theta }_p)\tan (\theta )\right) \\ \acute{\theta }={\sin }^{-1}\left(\sin \left(\theta \right)\sin ({\theta }_p)+\cos \left(\theta \right)\cos ({\theta }_p)\cos \left(\lambda -{\lambda }_p\right)\right) \end{gathered}$$

(34)

The wind fields are a combination of wind vectors of the solid rotation and that of the deformational flow, which are specified as follows:

$$\begin{gathered} u=u_0\left(\cos \left(\theta \right)\cos \left(\alpha \right)+\cos \left(\lambda \right)\sin \left(\theta \right)\sin \left(\alpha \right)\right) \\ +\mathsf{\omega }\left(\acute{\theta }\right)\left[\sin \left({\theta }_c\left(t\right)\right)\cos \left(\theta \right)-\cos \left({\theta }_c\left(t\right)\right)\cos \left(\lambda -{\lambda }_c\left(t\right)\right)\sin \left(\theta \right)\right] \\ V=-u_0\sin \left(\lambda \right)\sin \left(\alpha \right)+\mathsf{\omega }\left(\acute{\theta }\right)\left[\cos \left({\theta }_c\left(t\right)\right)\sin \left(\lambda -{\lambda }_c\left(t\right)\right)\right] \end{gathered}$$

(35)

where $\left({\lambda }_c\left(t\right),{\theta }_c\left(t\right)\right)$ denotes the center of the vortex, which is free to move anywhere on the sphere. The scaled tangential velocity is as follows:

$$\mathsf{\omega }\left(\acute{\theta }\right)=\{\begin{array}{c}{u}_0\frac{3\sqrt{3}}{2\rho \left(\acute{\theta }\right)}{\sec \mathrm{h}}^2\left(\rho \left(\acute{\theta }\right)\right)\tanh \left(\rho \left(\acute{\theta }\right)\right) if \rho \left(\acute{\theta }\right)\ne 0\\ 0 if \rho \left(\acute{\theta }\right)=0 \end{array}$$

(36)

where $\rho \left(\acute{\theta }\right)={\rho }_0\cos \left(\acute{\theta }\right)$ is the radial distance of the vortex. The background flow constant $u_0=2\pi R/\left(12 \mathrm{days}\right)$. The exact solution in non- dimensional units at time t is given by the following:

$$\varphi \left(\acute{\lambda },\acute{\theta },0\right)=1-\tan \mathrm{h}\left(\frac{\rho \left(\acute{\theta }\right)}{\gamma }\sin \left(\acute{\lambda }-\omega \left(\acute{\theta }\right)t\right)\right)$$

(37)

where $\gamma$ is a parameter defining the characteristic width of the frontal zone. We set ${\rho }_0=3$ and $\gamma =5$.With these parameters, the initial condition $\varphi \left(\acute{\lambda },\acute{\theta },0\right)$ is displayed in Figure 10a.

In this case, an icosahedron sphere triangle mesh with refinement level 5 is used, global RBF interpolation based on the Gauss kernel function was implemented and shape parameter c = 16, time step $\mathsf{\Delta }t=0.5$ h. Initially a vortex center is located at the following:

$$\left({\lambda }_c\left(t\right),{\theta }_c\left(t\right)\right){|}_{t=0}=\left(\frac{\pi }{2},0\right)$$

(38)

with this setup the other vortex center will be placed at the diametrically opposite point $\left(-\frac{\pi }{2},0\right)$. Thus the vortices have the initial and final positions (after one complete revolution) geographically located at 90° E and 90° W on the equator, as shown in Figure 10b,c.

In Figure 10c, it can be seen that a couple of vortices with correct rotation structures are evolving and there are no obvious oscillations around large gradients. The proposed scheme can simulate this complicated procedure successfully.

6. Conclusions

A computationally efficient and high-order backward trajectory semi-Lagrangian approach has been proposed for the solution of the transport equation on the sphere. We test several kinds of interpolation methods for the values at departure points.

The main conclusions are as follows:

(1)

Although the semi-Lagrangian method is not constrained by the CFL number, the time step for backward integration could not be too large arbitrarily, especially for flows of complex structures;

(2)

A high-order explicit integration method could be utilized to solve the characteristic equation in order to obtain more accurate locations of departure points. The SL method does not guarantee mass conservation; however, the mass is conserved well in test cases for numerical experiments when the solution is smooth;

(3)

Among different interpolation methods, results obtained by the linear interpolation method show too much diffusion, especially in the great gradient region of the solution. Compared to linear interpolation, more accurate results could be obtained by the Poly-LSQ method and the global or local RBF interpolation method. However, the value of the shape parameter for the basic function is key to the RBF interpolation. AN inappropriate value would cause a great loss of accuracy in the interpolation results, even numerical instability. It costs a lot of effort to select an optimal value. While this problem does not occur for the linear and Poly-LSQ methods. When the solution is not smooth, it tends to generate small ripples in the steep gradient region of the solution if the Poly-LSQ method, global RBFs, or local RBFs would be utilized. Therefore, the techniques to suppress numerical oscillations should be developed in further research;

(4)

Essentially, the proposed SL method is meshless and implemented in 3D Cartesian coordinates; it could be extended to other kinds of grids on the sphere directly.