A Semi-Lagrangian Type Solver for Two-Dimensional Quasi-Geostrophic Model on a Sphere

In this paper, we propose a numerical method based on semi-Lagrangian approach for solving quasi-geostrophic (QG) equations on a sphere. Using potential vorticity and stream-function as prognostic variables, two-order centered difference is suggested on the latitude-longitude grid. In our proposed numerical scheme, advection terms are expressed in a Lagrangian frame of reference to circumvent the CFL restriction. The pole singularity associated with the latitude-longitude grid is eliminated by a smoothing technique for the initial flow. Error analysis is provided for the numerical scheme.


Introduction
The quasi-geostrophic (QG) equations on a sphere are the major system of interest in weather forecasting and climate prediction [1], where the horizontal velocities are approximately geostrophic for extratropical synoptic-scale motions.It is important to investigate the QG equations because of both its intrinsic mathematical interest and its potential applications in dynamic meteorology and oceanography.The two-dimensional QG equations are originally established in modeling rotating fluids on the earth surface [2].Over the years, the study of two-dimensional QG equations has been an active research field; see, e.g., [3] [4].However, in spite of quite a number of contributions dealing with two-dimensional QG equations, there is little research on the evolution of two-dimensional QG equations on a sphere.Analytic solutions are rarely available due to the relatively complex nature of QG equations on a sphere; numerical simulations play an important role in the exploration of the QG equations, for example, [5] [6] [7].The accuracy of numerical investigation on QG models depends on many factors, including the knowledge of the initial state, the numerical methods employed, the resolution and so on.Some numerical methods based on finite difference methods can be applied to solve the QG equations on a spherical coordinate that may suffer from the pole singularity.Some efforts have made to alleviate the difficulties of the pole singularity.Kurihara [8] proposed a spherical grid system whose grid density on the globe is almost homogeneous, and the number of grid points per line of latitude near the poles is reduced.If the equations written in such coordinates are directly transformed into finite differences, excessive errors are committed in grids [9] [10].Some slight alteration was made in Kurihara's grid which alleviated the troubles [10] [11].On the other hand, in many global atmospheric applications, spatial discretization schemes are based on the spectral transform methods, in which solution fields are expressed as spherical harmonic expansions.Since the spherical harmonics are the natural representation of a two-dimensional field on the surface of a sphere, the spectral approach may provide better accuracy for the pole problem.The spectral transform method seems ideal for the spherical domain; however, it is too expensive for long time simulations.Especially at high spatial resolutions, since it requires associated Legendre transforms.
A framework of semi-Lagrangian semi-implicit methods was developed by Robert [12].It offers another possibility of developing fast numerical schemes for weather forecast models in complex systems.Over the years, researchers have devoted much attention to the further development and application of semi-Lagrangian methods.In [13] [14], Robert combined the semi-implicit integration scheme with a semi-Lagrangian treatment of the advection terms in a barotropic model.Robert et al. [15] extended this scheme to a multilevel model.They found that the time step could be increased by a further factor over an Eulerian semi-implicit scheme.For improving accuracy, McDonald [16] has incorporated a semi-Lagrangian treatment of advection in an efficient two-time-level integration scheme.The accuracy and stability of the semi-Lagrangian semi-implicit methods were investigated by McDonald [17].They showed that there is no restriction on the time step when considering the linear advection and following certain rules for the interpolation of functions at the trajectory points.Compared with Eulerian schemes, they have advantages of their own, such as high efficiency and easy to use in problems with nonuniform grids.In addition, semi-Lagrangian methods avoid the leading source of nonlinear instability in most wave-propagation problems since the nonlinear advection terms appearing in the Eulerian form of the momentum equations are eliminated when those equations are expressed in a Lagrangian approach [18].
Semi-Lagrangian methods and semi-implicit approach have become one of the most popular architectures used in numerical weather forecast and lots of other computational fluid dynamical applications; see e.g., [19] [20] [21].
In this paper, we consider the following non-dimensional 2D QG equations on a sphere: where Ψ is the potential vorticity; Φ is the stream-function; f is an external force; Ω is the rotational speed; r F is the planetary Froude number.r is the radius of sphere.It is challenging when initial flow is related to longitude, the discrete equation contains the 1 cosθ terms.It gives rise to larger truncation errors and reducing the order of convergence by one near the poles when θ approaches π 2

±
. As the stimulation propagated, these larger errors propagated over the sphere and eventually contaminated the rest of the solution.We provide a smoothing technique for the initial flow to overcome the above difficulty.Numerical experiments demonstrate the performance of our proposed method and investigate behaviors of QG model with geostrophic implications.
The rest of paper is structured as follows.In Section 2, we present the semi-Lagrangian discrete scheme of QG equations on a sphere.In Section 3, we carry out the detailed accuracy analysis of the method and explain some details.The conclusions are summarized in Section 4.

Discretization Schemes
In this section, we introduce the semi-Lagrangian scheme for QG equations.We will focus on the time discretization and ignore the space variable for the moment.The first equation of QG model is expressed as where Then the semi-Lagrangian scheme for above equation is ( ) Here subscripts a and d refer to evaluation at an arrival and departure point, respectively.We firstly iteratively calculate departure point using some first guess and an interpolation formula, the order of the interpolation is much less important.So we use linear interpolation here.Secondly, an cubic interpolation formula is adopted to evaluate Ψ at upstream point [22].Finally, evaluate Ψ at arrival points at time t t + ∆ .
By doing so, we have Using the two-stage trajectory calculation, Equation (3) can be computed as follows, ( ) ( ) . where , λ θ denote the longitude and latitude of the departure point.Comparing to the classical Runge-Kutta midpoint method, a slight difference is that the first stage uses ( ) ( ) , , , , , , , , , , ; the latter is more convenient if , dicted at the same time as Φ .After we located the grid cell containing the departure point of the characteristic, we adopt the four (eight) surrounding points in the interpolation scheme ; for grid cells close to the poles this means that we resort to points located across the pole.
Let us introduce, for example, cubic Lagrange interpolation.Set , i a ii q j b jj p − = + − = + , where ( ) , ii jj indicates the lower left corner of the cell containing the departure point and 0 , 1 p q < < .We then define ( ) , , , , P p p p p = ( )( ) and similarly for q.When 1 1 jj NY < < − , let ϕ indicates the variable to be interpo- lated, we denote .
With the propositions above, we have T T , , ii q jj p Q P W ϕ + + = alternatively, it can be expanded as follows, p q q q q p q q q q p q q q q p q q q ϕ ϕ ϕ ϕ ϕ ii jj ii jj In the case 1 jj = , only the first term on the right hand side of Equation ( 5) needs to be changed as .
ii NX jj ii NX jj ii NX jj ii NX jj p q q q q ϕ ϕ ϕ ϕ − , a similar procedure can be adopted.When 0 jj = (and analo- gously jj NY = ), the first two terms on the right hand side of Equation ( 5) should be substituted as p q q q q p q q q q ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ and similarly for the case jj NY = .
For spatial approximation, we adopt unstaggered grid that is uniform in longitude and latitude.Stream-function Φ and potential vorticity Ψ are collocated on nodes ( ) , i j that are the corners of grid cells with spatial increments ( ) the Equator and 1 2 j NY = + the North Pole) [22].From the second equation of (1) we have For the second term of right hand side of Equation ( 6), we use the following approximation With the definition of the discrete operators above, the spatial discretization for ( 6) is defined as follows.When 2 1 j NY ≤ ≤ − , the second-order centered difference scheme for ( 6) is When 1 j = , the grid points corresponding to 1 j − located across the pole, the modified ( 8) is Similarly when j NY = , the scheme is cos 2 cos sin .

Error Analysis
In this section, we will carry out the error estimate of our new algorithm.We let 0 f = for simplicity.Suppose that ( ) is the sufficiently smooth and continuous solution.The truncation error E satifies , 2 cos sin cos , .
The semi-Lagrangian scheme of the first equation in (1) reads Here subscript d refers to evaluation at an departure point.By expanding ( ) , ,t λ θ Ψ in a Taylor series about the estimated departure point and evaluating an expression of the form ( ) where ( ) and the summation represents an ( ) The first term of right hand side of Equation ( 18) is determined by the error in the trajectory calculation, and the second term of right hand side of equation ( 18) is determined by the error in the interpolation of n Ψ to the departure point.
We first calculate the error of the Runge-Kutta discretization 2 ( ) ( ) ( ) . 2 Similarly we can get ( ) ( ) ( ) Employ Taylor expansions of Substituting ( 22) into (23) gives Which that the Runge-Kutta scheme (4) generates an ( ) tion toward the total error in the semi-Lagrangian approximation.Now suppose that the velocity data , u v are available only at discrete locations on the space-time mesh.Ideally the velocity at time 1 2 n t + would be computed by interpolation between times n t and 1 n t + .Such interpolation cannot, however, be performed when semi-Lagrangian methods are used to solve prognostic equations for the velocity itself, because the velocity at 1 n t + will be needed for the trajectory calculations before it has been computed.This problem is generally avoided by extrapolating the velocity field forward in time using data from the two previous time levels such that ( ) ( ) Suppose that the extrapolated velocity field , , , , , Substituting the preceding equation into (20) shows that the use of ,

Conclusion
In this paper, we present a numerical method which combines semi-Lagrangian me-thod with two-order centered difference scheme for solving two-dimensional quasigeostrophic equations on a sphere.In our approach, potential vorticity and streamfunction are used as prognostic variables.Advection terms are expressed in a Lagrangian frame of reference to avoid the necessity of stable constraint.The pole singularity is eliminated by means of employ a smoothing technique.An error analysis is presented.
In view of the spherical coordinates, the grid is nonuniform in ,x y .The size of the cells reduces when we move toward the poles.No functional variables are collocated on the poles.More specifically, for fixed NX and NY , we pose 2π )