A positivity-preserving conservative Semi-Lagrangian Multi-moment Global Transport Model on the Cubed Sphere

A positivity-preserving conservative semi-Lagrangian transport model by multi-moment finite volume method has been developed on the cubed-sphere grid. In this paper, two kinds of moments, i.e. point values (PV moment) at cell boundaries and volume integrated average (VIA) value, are defined within a single cell. The PV moment is updated by a conventional semi-Lagrangian method, while the VIA moment is cast by the flux form formulation that assures the exact numerical conservation. Different from the spatial approximation used in CSL2 (conservative semi-Lagrangian scheme with second order polynomial function) scheme, a monotonic rational function which can effectively remove non-physical oscillations and preserve the shape, is reconstructed in a single cell by the PV moment and VIA moment. The resulting scheme is inherently conservative and can allow a CFL number larger than one. Moreover, the scheme uses only one cell for spatial reconstruction, which is very easy for practical implementation. The proposed model is evaluated by several widely used benchmark tests on cubed-sphere geometry. Numerical results show that the proposed transport model can effectively remove unphysical oscillations compared with the CSL2 scheme and preserve the numerical non-negativity, and it has the potential to transport the tracers accurately in real atmospheric model.


Introduction
Global advection transport describes the motion of various tracers, which is a basic process in atmospheric dynamics. The numerical result of advection transport is significant in developing general circulation models (GCMs). Traditional latitude-longitude grid is very easy for application but has singularities at poles. Moreover, its nonuniform grid system would seriously affect computational efficiency. To address these issues, quasi-uniform grid systems without singularities or with weak singularities, such as the cubed-sphere grid, yin-yang grid and icosahedral grid, are becoming more and more popular in developing global transport model. Among them, cubed-sphere grid is much more attractive due to its computational merits, such as locally structured grid, local mass conservation and quasiuniform grid. Some transport models based on the cubed-sphere grid can be seen in [2,3,5,6,7,11,13,15]. In this paper, we consider the cubed-sphere with gnomonic projection for our transport model.
Semi-Lagrangian method [14] is widely used in transport model for it permit a large CFL number without reducing accuracy. But the traditional semi-Lagrangian method has a serious shortcoming, that is, it is not mass conservative. To address this problem, many efforts have been made to develop the conservative semi-Lagrangian scheme. Nakamura et al. [9] proposed such a scheme based on their previous Constrained Interpolation Profile (CIP) method [19], and they called it CIP-CSL. In their method, the point values at cell boundaries and the cell-averaged value are used to construct the piecewise interpolation profile. The point values are updated by the semi-Lagrangian approach, while the cell average or volume average values are calculated by the flux-form formulation. The semi-Lagrangian approach permits a large CFL number, and the flux-form formulation of updating cell averaged values makes the scheme inherently conservative in terms of cell averaged values. Xiao and Yabe [17] introduced a slope limiter in CIP-CSL scheme to suppress oscillations around discontinuities, but the stencil for spatial reconstruction extended, from one cell to three cells. Instead of the cubic polynomial function used in CIP-CSL2, Xiao et al. [18] used a rational interpolation as a substitution, they called it CSLR, which used only one cell as stencil to construct interpolation function and remove non-physical oscillations simultaneously, but this scheme cant completely preserve positivity. In this paper, we make some modification on the CSLR scheme to get a non-negative scheme and extend it to the cubed-sphere grid to develop a global transport model. The paper is organized as follows. In section 2, we will review the 1D algorithm of CSLR scheme its modification. In section 3, we extend this formula to the cubed-sphere grid. Section 4 presents several kind of benchmark tests to evaluate the performance of the proposed global transport model. A brief summary is given in section 5.

Spatial reconstruction
To reconstruct the spatial profile, two kinds of moments are introduced in each cell, as illustrated in Fig. 1, PV moments at cell boundaries and VIA moment in C i (i = 1, 2, . . . , N) are defined as • The PV moments P q i±1/2 (t) = q x i±1/2 , t (1) where q(x, t) is the transport quantity, ∆x is the grid spacing. Using these three moments, a rational function can be reconstructed in cell i and the coefficients are determined by R i x i+1/2 = P q i+1/2 (5) 1 ∆x An alternative rational function can be expressed as β i is the same as in Eq. (3). Using Eq. (4) (5) and (6), other coefficients can be determined. Details can be found in [17].

Moments updating
After getting the piecewise spatial reconstruction on the entire domain, the updating procedure can be conducted. Consider the following one-dimensional transport equation • Updating the PV moments: The PV moments are updated by the traditional semi-Lagrangian approach. Rewrite Eq. (8) in an advection form it can be viewed as an advection equation plus a source term −q ∂u ∂x . The advection part is calculated by the semi-Lagragian concept P q n+1 q i−1/2 = R n I x ip (10) where x ip is the departure point at previous time step t = n∆t corresponding to the arrival point x i−1/2 at next time step t = (n + 1)∆t, the subscript I is the index of the cell which contains the departure point x ip , and departure point is simply calculated by where u x 1 ip is the velocity at predict point x 1 ip = x i−1/2 − u i−1/2 ∆t . In general, x 1 ip would not be identical with the point at cell interface, and the velocity at predict point x 1 ip is calculated by linear interpolation using known velocity at two interface of the cell which contains x 1 ip . On the cubed-sphere grid, if predict point is outside of patch boundary, the departure point is calculated by The source term in Eq. (9) is simply approximated by • Updating the VIA moment: The VIA moment is updated by the flux-form concept where g i+1/2 is the flux of q going through the boundary x i+1/2 during [n∆t, (n + 1)∆t], which is calculated by analytically integrating the interpolation function along the trajectory of x i+1/2

Modifications for positivity preserving and reducing overshoots
For a problem with a lower boundary q min = 0 and upper boundary q max . The point values calculated by Eq.(13) may produce negative values and values exceed q max . An easy and effective modification for the PV moments is used In addition, even the PV moments are all no less than 0, negative values may also appear when a valley near lower boundary is advected. As illustrated in Fig. 2, when the PV moments at cell boundary are bigger than the VIA moment, the reconstructed rational function would produce undershoots. In this condition, if the VIA is around lower boundary, negative values may appear. Similarly, when a peak is advected, overshooting of PV moments may also appear. Thus, a further modification is needed where εis a small parameter, such as ε = 10 −3 , to avoid the valley or peak to be advected. We should note that, the modification of Eq.(17) can guarantee the spatial approximation profile is above 0, and by using the flux-form formula of VIA moment we can get an absolutely positive value. But around the upper boundary, by using Eq. (18) we can only get the spatial approximation profile below q max , but cannot ensure the flux out is no less than the flux in, so overshooting of VIA moment around upper boundary may still exist. Thus by using these modification of PV moments, the numerical result can strictly preserve positive and reduce overshooting.
In this paper, the scheme using Eq.(7) as spatial reconstruction is called CSLR1, the scheme with two steps modification is called CSLR1-M. When β = 0 in Eq. (7), the scheme becomes CSL2 [20]. Given the known point values and cell average values at previous time step, the updating procedure can be summarized as follows: 1) Using Eq. (7), the spatial approximation of transport property on the whole domain can be determined.  (17) and (18) to ensure a positive value and suppress overshoots at next time step.

Extending to the cubed-sphere grid
In this section, we extend the CSLR1 and its modification CSLR1-M scheme to the cubed-sphere grid to develop a global transport model. The cubed-sphere grid we used in this paper is constructed by equiangular central projection, by this way six identical local coordinate (α, β) = [−π/4, π/4] are constructed.
The two-dimensional transport equation in local coordinate can be written as where √ G is the Jacobian of transformation, and u 1 , u 2 is the contravariant components on the local coordinate, details can be found in [11].
In the two-dimensional case, as shown in Fig.3, four kinds of moments are introduced within C i j : • Volume integrated average (VIA): where ∆α and ∆β are grid spacing in the α and β direction respectively.
• Point value (PV): four point-values are located at the vertices • Line integrated average values along α direction: • Line integrated average values along β direction: To ensure conservation property, we followed [5] to divide the cube-sphere grid into three direction, see in Fig.4. Firstly, we conduct the update procedure in ξ direction, i.e. α direction on patch 1,2,3,4 for ∆t/2. It should be noted that the moments along patch boundaries are only updated for once. As shown in Fig.5, A is the arrival point on patch boundary and A d is the corresponding departure point on Patch 4, the point value of point A is calculated on Patch 4 and the flux across A is calculated by integrating the spatial approximation profile on Patch 4 along A d A. If A d is on Patch 1, vice versa. By the same way, we update the moments in η direction, i.e. β direction on patch 1,3,5,6 for ∆t/2 ; update in ζ direction, i.e. α direction on patch 5,6 and β direction on patch 2,4 for ∆t ; then another ∆t/2 in direction η and ξ direction separately to complete a full update procedure on the sphere geometry.

Numerical experiments
In this section, we use several widely used benchmark tests to verify the performance of the proposed transport model. Including solid body rotation, moving vortices and deformational flow test on the spherical mesh.
The normalized errors and relative maximum and minimum errors proposed by Williamson et al. [16] are used: where Ω is the whole computational domain, q and q t refer to numerical solutions (volume integrated average in our paper) and exact solutions respectively.

Solid-body rotation tests
Solid-body rotation test [16] is widely used in two-dimensional spherical transport model to evaluate the performance of a transport model. The wind components in the latitude-longitude coordinates (λ, θ) are defined as where (u s , v s ) is the velocity vector, u 0 = 2πR/12days, which means it takes 12 days to complicate a full rotate on the sphere, R is the radius of the sphere, and α is a parameter which controls the rotation angle. In this test, two kinds of initial conditions are used, including a cosine bell and a step cylinder. (a) Solid body rotation of a cosine bell The initial condition of a cosine bell test is specified as where r d is the great circle distance between (λ, θ) and the center of the cosine bell, located at (3π/2, 0), r 0 = 7πR/64 is the radius of the cosine bell, h 0 = 1. The normalized errors on 30 × 30 × 6 meshes and with 256 time-steps compared with other existing published Semi-Lagrangian scheme, the PPM-M scheme [21] and CSLAM-M [7], are presented in Table 1. With a qusi-smooth initial condition and quiet simple velocity field, the modification procedure has little effect on the numerical result. The result shows that CSLR1 and CSLR1-M has almost the same result. And our scheme is comparable to PPM-M scheme, the result in near pole flow direction (α = π/2 and α = π/2 − 0.05) is better than CSLAM-M scheme.
To check the influence of the weak singularities at 8 the vertices of the cubed-sphere gird, this test is conducted with α = π/4 to pass through four vertices. The history of normalized errors (CSLR1 and CSLR1-M are almost the same, here we only present the result of CSLR1-M) are shown in Fig.6. We can see that the normalized errors have no visible fluctuations when the flow passes four weak singularities.
To demonstrate the ability of using large CFL number to transport, we use 72 time-steps to complete one revolution. The normalized errors compared with CSLAM-M are shown in Table 2, the proposed scheme in this paper has larger l 1 and l 2 error, while l ∞ error is smaller than CSLAM-M. Compare with the result using 256 time-steps with the same grid resolution and flow direction, a larger CFL number condition even get a better result in terms of norm errors.
(b) Solid body rotation of a step cylinder   A non-smooth step cylinder is calculated to evaluate the non-oscillatory property. The initial distribution is specified as where r d is the great circle distance between (λ, θ) and (3π/2, 0), which is the center of the step cylinder, r 1 = 2/3R and r 2 = 1/3R. In this test, we set α = π/4, which is the most challenging case, the step cylinder moves through four vertices and two edges to complete a full rotation. Here we use 90 × 90 × 6 meshes and with 720 steps to conduct this test. The numerical results after 12 days are shown in Fig. 7, we can see that the CSL2 scheme will generate oscillations around the discontinuities. But by using the CLSR1 and CSLR1-M approach, these unphysical oscillations are effectively removed. The relative maximum and minimum of CSL2 are q max = 3.77 × 10 −2 and q min = −2.80 × 10 −2 , for CLSR1 it is q max = 1.36 × 10 −3 and q min = 0, for CSLR1-M is q max = −9.23 × 10 −5 and q min = 0. The result indicates that CSLR1-M can effectively remove overshoots in this test. The history of normalized mass errors is given in Fig.  8, which shows that the normalized mass error is within the tolerance of machine precision, therefore the proposed global transport model is exactly mass conservative during simulation procedure.

Moving vortices on the sphere
The second benchmark test we used is the moving vortices proposed by Nair et al. [12], the wind component of which is a combination of rotation test and two vortices, it is much more complicated than the solid-body rotation test. The velocity field on the sphere is specified as where u s and v s are calculated by Eq.(29) and (30), the rotation angle of this test is set to be α = π/4. ρ 0 = 3, λ c (t) and θ c (t) are the center of the moving vortex at time t, the calculation procedure of λ c (t) and θ c (t) can be found in [12].
The tracer field is defined as where γ is a parameter to control the smoothness of tracer field, (λ , θ ) is the rotated spherical coordinates, which can be calculated by θ (λ, θ) = sin −1 sin θ sin θ p + cos θ cos θ p cos λ − λ p and λ p , θ p = (π, π/2 − α) is the North Pole of the rotated spherical coordinate. In this test, we followed [13] to set γ = 10 −2 to conduct a large gradient in tracer distribution. When t = 0 in Eq. (38), we get the initial condition. The test is conducted on 80 × 80 × 6 grid and uses 400 steps to move for 12 days. Contour plot is shown in Fig.  9, compared with the exact solution, the proposed scheme can correctly simulate this complicated procedure. The errors of this test are presented in Table 3, we can see that the normalized errors l 1 ,l 2 , l ∞ and q max of CSLR1 and CSLR1-M have little difference, CSLR1 scheme have small undershoots, while CSLR1-M can completely preserve non-negativity.

Deformational flow test
The last benchmark test used in our paper is deformational flow test proposed by Nair and Lauritzen [10], which is a very challenging test case. The flow field is nondivergent and time-dependent: where κ = 2, T = 5, and λ = λ − (2πt/T ). Two kinds of initial conditions are checked here, including twin slotted cylinder case to evaluate the positivity preserving property and correlated cosine bells to evaluate the nonlinear correlations between tracers [8]. By the given flow field, the initial distributions will change into thin bars in the first half period, then return to its initial state in the second half period.
(a) Deformation of twin slotted cylinder The initial condition is defined as where r 0 = 0.5 and r i (i = 1, 2) means the great circle distance between the center of the two slotted cylinder and a given point. The center of the two slotted cylinders are located at (λ 1 , θ 1 ) = (5π/6, 0) and (λ 2 , θ 2 ) = (7π/6, 0), respectively. The numerical result of deformational flow of CSLR1-M scheme on 90 × 90 × 6 grid and with 390 time-steps (local CFL max is about 3) is shown in Fig. 10. The result at first half period is shown in Fig. 10(b), the two slotted cylinders are deformed into two thin filaments by the background flow field. Fig. 10(c) is result at final time, it is clearly that the proposed scheme can correctly reproduce this complicate deformational flow. Normalized errors are shown in Table 4. Both schemes are non-negative at final step, but the CSLR1 scheme will generate negative values at some time step, while CSLR1-M can always preserve positive along the simulation procedure. And the CSLR1-M scheme reduces the overshooting slightly.
where h i = 1 2 1 + cos πr i r 0 for i = 1, 2. To compare with other schemes, we set this test in 60 × 60 × 6 meshes and with 600 time steps. The result is shown in Table 5, the error by our scheme is slightly larger than those by CSLAM and SLDG scheme.
(c) Deformation of correlated cosine bells To check the ability of preserve nonlinear correlated relations between two tracers, we used two kinds of tracers. One is the twin cosine bells defined by Eq.(44), and the other one is the correlated cosine bells:  where Ψ(q) = −0.8q 2 + 0.9. This test is conducted on 90 × 90 × 6 meshes with 1800 time-steps. The scatter plot of numerical result at t = T/2 is shown in Fig. 11. The solution of cosine bells is in the X-direction, and the correlated cosine bells is in the Ydirection. The result shows that CSL2 scheme has visible undershoots, while CSLR1 scheme can effectively remove these unphysical undershoots. Table 6 lists mixing diagnostics. CSLR1 scheme produce more real mixing diagnostic and less range-preserving unmixing. Because of these two schemes are not shape-preserving, some overshooting appears as expected. But the CSLR1 scheme has less overshooting which is consist with the scatter plot.

Summary
In this paper, a non-negative and conservative semi-Lagrangian transport scheme based on multi-moment concept has been developed on the cubed-sphere grid. By using two kind of moments such as point value and volume integrated average, a rational function is constructed as spatial approximation function within a single cell. To get a nonnegative scheme, two kinds of modifications are conducted on the original CSLR1 scheme. The benchmark tests above demonstrate that CSLR1 scheme is non-oscillatory and can preserve the non-linear correlations between tracers. In case of valley of the transported field, CSLR1 scheme is unable to preserve positivity due to the rational function properity. The definite positivity can be achieved by the easy and effective modifications. In addition, the semi-Lagrangian approach allows the proposed scheme to use large time step, which can greatly improve computational efficiency. Unlike many other non-oscillatory scheme, for example, the WENO limiter used in DG and finite volume method, the proposed model in this paper uses only one cell as stencil for spatial reconstruction, which makes this model very efficient and easy for practical application.