A Positivity-preserving Conservative Semi-Lagrangian Multi-moment Global Transport Model on the Cubed Sphere

Global advection transport describes the motion of various passive tracers in the atmosphere, which is a basic process in atmospheric dynamics. The advection transport model is important in developing general circulation models (GCMs). The traditional latitude-longitude grid is very easy for application but has singularities at the poles. Moreover, its nonuniform grid system would also 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 models. Among those grids, the cubed-sphere grid is usually preferred due to its computational merits, such as locally structured grid and quasi-uniform grid. Recently, many transport models have been developed on the cubed-sphere grid, such as the discontinuous Galerkin transport models (Nair et al., 2005; Guo et al., 2014, 2016), the conservative semi-Lagrangian multi-tracer (CSLAM) model (Lauritzen et al., 2010), the finite volume transport model (Norman and Nair, 2018), and the multi-moment transport models (Chen et al., 2011; Tang et al., 2018). In this study, the cubed-sphere grid with gnomonic projection is adopted for our transport model.

The semi-Lagrangian method is a popular choice for developing a global transport model, since it allows a large time step without reducing accuracy. The traditional semi-Lagrangian method defines a set of parcels that arrive at the Euler computational grid at every time step, and then these parcels are traced back to find their departure locations at the previous time step. A review of the semi-Lagrangian method can be seen in Staniforth and Côté (1991). However, the traditional semi-Lagrangian method has a serious shortcoming regarding mass conservation. To deal with this, many efforts have been made to develop the conservative semi-Lagrangian method. The finite-volume semi-Lagrangian (FVSL) method is a popular one which can be probably separated into two categories (departure volume based and flux form based). The departure volume based FVSL method initially finds the departure volume and then remaps the departure volume from the given Euler computational grid. Examples of using this method can be seen in Nair and Machenhauer (2002), Nair et al. (2002), Zerroukat et al. (2002), and Lauritzen et al. (2010). The flux form based FVSL method calculates the flux traveling across the interfaces of a cell and uses a flux form formulation to update. Lin and Rood (1996) is a typical example of this kind.

Nakamura et al. (2001) proposed a flux-form FVSL method based on their previous Constrained Interpolation Profile (CIP) scheme (Yabe and Aoki, 1991), calling it CIP-CSL. In their method, the point values at cell boundaries and the cell-averaged value are used to reconstruct 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 time step, and the flux-form formulation of updating cell-average values makes the scheme inherently conservative in terms of cell-integrated average values. Xiao and Yabe (2001) introduced a slope limiter in the 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. (2002) utilized a rational interpolation as an alternative, calling it the CSLR scheme, which used only one cell as stencil to reconstruct the interpolation function and could remove nonphysical oscillations simultaneously. However, this scheme can’t completely preserve positivity. In this paper, we make some modifications on the CSLR scheme to make it non-negative and extend it to the cubed-sphere grid to develop a global transport model.

The paper is organized as follows. In section 2, we introduce the algorithm of the CSLR method and its modifications on the Cartesian geometry. In section 3, we extend this formula to the cubed-sphere grid. Section 4 presents several kinds of benchmark tests to evaluate the performance of the proposed global transport model. And a brief summary is given in section 5.

To reconstruct the spatial approximation profile, two kinds of moments are introduced in each cell, as illustrated in Fig. 1. Point value (PV) moments at cell boundaries and the volume integrated average (VIA) moment in C_i ($i = 1,2,...,N$) are defined as:

Figure 1.  Illustration of moments in one dimension.

● The PV moments

● The VIA moment

where $q\left({x,t} \right)$ is the transport quantity and $\Delta x = {x_{i + \frac{1}{2}}} - {x_{i - \frac{1}{2}}}$ is the grid spacing.

In the CSLR1 method, a rational function is reconstructed within a single cell $i$

and three constraint conditions are applied as

after which the coefficients can be determined as

where ${\beta _i}$ is predetermined in Eq. (3) [see Xiao et al. (2002) for details], ${\gamma _i} = 1 + {\beta _i}\Delta x$, and $\varepsilon $ is a very small number, such as $\varepsilon = 1 \times {10^{ - 20}}$, for avoiding a zero denominator in Eq. (10).

Consider the following one-dimensional transport equation,

where $u$ is the velocity.

● Updating the PV moments:

The PV moments are updated by the traditional semi-Lagrangian approach. Rewriting Eq. (11) in an advection form gives

and it can be viewed as an advection equation plus a source term, $- q {{\partial u}}/{{\partial x}}$. The advection part is calculated by the semi-Lagrangian concept

where ${\xi _{i - \frac{1}{2}}}$ is the departure point at previous time step $t = n\Delta t$ corresponding to the arrival point ${x_{i - {\frac{1}{2}}}}$ at next time step $t = \left({n + 1} \right)\Delta t$, where $\Delta t$ is the time interval, the subscript $I$ is the index of the cell which contains the departure point ${\xi _{i -\frac{1}{2}}}$, and departure point is simply calculated by

where $u\left({\xi _{i -\frac{1}{2}}^ * } \right)$ is the velocity at the first guess point $\xi _{i - \frac{1}{2}}^ * = {x_{i - \frac{1}{2}}} - {u_{i - \frac{1}{2}}}\Delta t$. In general, $\xi _{i - \frac{1}{2}}^ *$ would not be identical with the point at the cell interface, so the velocity at the first guess point $\xi _{i - \frac{1}{2}}^ *$ is calculated by linear interpolation using known velocity at two interfaces of the cell which contains $\xi _{i - \frac{1}{2}}^ *$.

After the “source term” in Eq. (12) is simply approximated by a central difference formulation, the semi-discretized form of the transport equation can be written as

● Updating the VIA moment:

The VIA moment is updated by the flux-form concept

where ${g_{i + {\frac{1}{2}}}}$ is the flux as $q$ goes through the boundary ${x_{i + {\frac{1}{2}}}}$ during time interval $\left[ {n\Delta t,\,\left({n + 1} \right)\Delta t} \right]$, which is calculated by analytically integrating the interpolation function along the trajectory of ${x_{i + {\frac{1}{2}}}}$

Preserving the positivity of certain physical quantities requires that the minimum value ${q_{\min }}$ should not be less than zero. However, the point values calculated by Eq. (15) may produce negative values. Since the conservation of the PV moment is not required in the context of the multi-moment finite volume scheme, an easy and effective modification for the PV moments is used:

Despite this modification, in the specific case when a “valley” shape near the lower boundary is transported, the negative values may still appear. As illustrated in Fig. 2, if the PV moments at the cell boundary are bigger than the VIA moment, the reconstructed rational function would produce “undershoots”. Thus, a further modification of the approximation profile is needed:

Figure 2.  Illustration of the rational reconstruction when a “valley” is advected.

where ${q_{\max }}$ is the maximum value of transport quantity and $\varepsilon $ is a small parameter, such as $\varepsilon = {10^{ - 3}}$. It should be noted that the modification of Eq. (19) can guarantee the spatial approximation profile is above zero, and by using the flux-form formula of VIA moment we can obtain an absolutely positive result. Therefore, after utilizing these two modifications the numerical result can strictly preserve positivity.

In this paper, the scheme using Eq. (3) for spatial reconstruction is called CSLR1, and the scheme with two-step modifications is called CSLR1-M hereafter. When $\beta = 0$ in Eq. (3), the scheme reduces to CSL2 (Yabe et al., 2001).

Given the known PVs and VIAs at the previous time step, the CSLR1-M algorithm updating procedure can be summarized as follows:

1) Using Eq. (3) and the modification of Eq. (19), the reconstructed profile within each cell can be determined.

2) Point values are updated by Eq. (13) and Eq. (15).

3) Cell-averaged values are updated by Eq. (16).

4) Modifying the PV moments by Eq. (18) ensures positive PV moments at next time step.

It is noted that given the monotonicity of rational function and the PVs at cell boundaries as predicted variables, the CSLR1-M scheme can easily facilitate a positive-preserving property, as shown in this paper.

A second order Strang dimension-splitting time-stepping (Strang, 1968) technique is adopted to extend the 1D algorithm to the two-dimensional Cartesian case. For the sake of simplicity, we collectively define the 1D CSLR1 and CSLR1-M algorithm as

which means given the known point values ${\overline {{}^Pq} ^n}$ at cell interfaces, volume integrated values ${\overline {{}^Vq} ^n}$, and the velocity $u$ at time $t = n\Delta t$, we can use the 1D CSLR algorithm to update the PVs ${\overline {{}^Pq} ^{n + 1}}$ and VIAs ${\overline {{}^Vq} ^{n + 1}}$.

In the two-dimensional case (see Fig. 3), four kinds of moments are introduced within cell ${C_{ij}}$:

Figure 3.  Illustration of moments defined in a two-dimensional case.

● Volume integrated average (VIA):

where $\Delta x$ and $\Delta y$ are grid spacing in the $x$- and $y$-directions, respectively.

● Point value (PV): four point-values located at vertices

● Line-integrated average values along $x$-direction

● Line-integrated average value along $y$-direction

Consider the two-dimensional transport equation in the Cartesian coordinates:

where $u$ and $v$ are the velocity in the $x$- and $y$-directions, respectively.

By using the dimension-splitting technique, the transport Eq. (25) is split into two 1D equations:

The updating procedures for a time step $\Delta t$ are as follows:

1) Update Eq. (26) for ${{\Delta t} / 2}$. By using the ${\rm{CSLR}} $$ \left({{{\overline {{}^Pq} }^n},{{\overline {{}^{Lx}q} }^n},u,{{\overline {{}^Pq} }^{\left(1 \right)}},{{\overline {{}^{Lx}q} }^{\left(1 \right)}}} \right)$ algorithm [where the superscript (1) means the result of step 1 and hereafter the superscript $\left(x \right)$ means the result of step $x$ similarly], PV moments and line integrated values along the $x$-direction are updated. And by ${\rm{CSLR}}\left({{{\overline {{}^{Ly}q} }^n},{{\overline {{}^Vq} }^n},u,{{\overline {{}^{Ly}q} }^{\left(1 \right)}},{{\overline {{}^Vq} }^{\left(1 \right)}}} \right)$, in which ${\overline {{}^{Ly}q} ^n}$ is viewed as point values in the 1D algorithm, the line-integrated values along the $y$-direction and VIA moments are updated. The four kinds of moments defined in our method are all updated in the $x$-direction for ${{\Delta t} / 2}$.

2) Similarly, update Eq. (27) for $\Delta t$ by using the ${\rm{CSLR}}\left({{{\overline {{}^Pq} }^{\left(1 \right)}},{{\overline {{}^{Ly}q} }^{\left(1 \right)}},v,{{\overline {{}^Pq} }^{\left(2 \right)}},{{\overline {{}^{Ly}q} }^{\left(2 \right)}}} \right)$ and ${\rm{CSLR}}\left({{\overline {{}^{Lx}q} }^{\left(1 \right)}},{{\overline {{}^Vq} }^{\left(1 \right)}}, $$ v,{{\overline {{}^{Lx}q} }^{\left(2 \right)}},{{\overline {{}^Vq} }^{\left(2 \right)}} \right)$ algorithms.

3) Update Eq. (26) for another ${{\Delta t} / 2}$ by using the ${\rm{CSLR}}\left({{{\overline {{}^Pq} }^{\left(2 \right)}},{{\overline {{}^{Lx}q} }^{\left(2 \right)}},u,{{\overline {{}^Pq} }^{n + 1}},{{\overline {{}^{Lx}q} }^{n + 1}}} \right)$ and ${\rm{CSLR}}\left({{\overline {{}^{Ly}q} }^{\left(2 \right)}},{{\overline {{}^Vq} }^{\left(2 \right)}}, $$ u,{{\overline {{}^{Ly}q} }^{n + 1}},{{\overline {{}^Vq} }^{n + 1}} \right)$ algorithms.

In this section, we extend the proposed scheme to the cubed-sphere grid to develop a global transport model. The quasi-uniform cubed-sphere grid (Sadourny, 1972) with equiangular central projection is adopted in this paper, as shown in Fig. 4, which has six identical cube faces with local coordinate $\left({\alpha,\,\beta } \right) = \left[ { - {\pi /4},\,{\pi / 4}} \right]$. It is worth mentioning that the conventional tropic-belt arrangement (Nair et al., 2005) is used in this paper although the staircase arrangement (Chen, 2021) is a good interlock pattern which has better symmetry for patch information exchange. The two-dimensional transport equation in local coordinates can be written as

Figure 4.  Schematic of a cubed-sphere grid with $12 \times 12 \times 6$ meshes.

where $\sqrt G $ is the Jacobian of transformation and $\left({{u^1},{u^2}} \right)$ are the contravariant components on the local coordinate $\left({\alpha, \beta } \right)$. Details of the transformation calculation can be found in (Nair et al., 2005). Note that the transport equation in the computational space of the cubed-sphere grid is the same as the two-dimensional Cartesian coordinate’s transport equation, besides substituting $\left({x,\,y} \right)$ with $\left({\alpha,\,\beta } \right)$, which makes it convenient to update the solutions as in the 2D Cartesian case.

As shown in Fig. 5, we divided the cubed-sphere grid into three directions $\left({\xi,\,\eta,\,\zeta } \right)$, and Eq. (28) is split into three sequential 1D equations along three directions (Guo et al., 2014):

Figure 5.  Schematic for three directions on the cubed-sphere grid. Top left is the ξ-direction along the α-direction on Patch 1, Patch 2, Patch 3, and Patch 4; Top right is the η-direction along the β-direction on Patch 1, Patch 3, Patch 5, and Patch 6; bottom is the ζ-direction along the β-direction on Patch 2 and Patch 4 and along the α-direction on Patch 5 and Patch 6.

where ${U_\xi }$, ${U_\eta }$, and ${U_\zeta }$ are the velocity along $\xi $-, $\eta $-, and $\zeta $-directions, respectively.

Then, the numerical solutions are updated by the splitting algorithm. In each direction, the moments are updated by the one-dimensional algorithm, similar to the case in the two-dimensional Cartesian geometry. Given the known point values, line-integrated values, and cell-integrated values, the final updating procedures on a sphere for a time step $\Delta t$ are summarized as follows:

1) Update in $\xi $-direction for ${{\Delta t} / 2}$ by using ${\rm{CSLR}}\left({{{\overline {{}^Pq} }^n},{{\overline {{}^{L\xi }q} }^n},{U_\xi },{{\overline {{}^Pq} }^{\left(1 \right)}},{{\overline {{}^{L\xi }q} }^{\left(1 \right)}}} \right)$ (where $\overline {{}^{L\xi }q} $ is the line-integrated average value along the $\xi $-direction) and ${\rm{CSLR}}\left({{\overline {{}^{L{\xi _{\bot} }}q} }^n}, $$ {{\overline {{}^Vq} }^n},{U_\xi },{{\overline {{}^{L{\xi _{\bot} }}q} }^{\left(1 \right)}},{{\overline {{}^Vq} }^{\left(1 \right)}} \right)$(where the mark ${\xi _ \bot }$ means the direction perpendicular to $\xi $, i.e. $\overline {{}^{L{\xi _ \bot }}q} $ is the line-integrated average values along the direction perpendicular to $\xi $, and $\overline {{}^{L{\eta _ \bot }}q} $ and $\overline {{}^{L{\zeta _ \bot }}q} $ are the line-integrated average values along the direction perpendicular to $\eta $ and $\zeta $, respectively);

2) Update in $\eta $-direction for ${{\Delta t} / 2}$ by using ${\rm{CSLR}}\left({{{\overline {{}^Pq} }^{\left(1 \right)}},{{\overline {{}^{L\eta }q} }^{\left(1 \right)}},{U_\eta },{{\overline {{}^Pq} }^{\left(2 \right)}},{{\overline {{}^{L\eta }q} }^{\left(2 \right)}}} \right)$ and ${\rm{CSLR}}\left({{\overline {{}^{L{\eta _{\bot} }}q} }^{\left(1 \right)}},{{\overline {{}^Vq} }^{\left(1 \right)}}, $$ {U_\eta },{{\overline {{}^{L{\eta _{\bot} }}q} }^{\left(2 \right)}},{{\overline {{}^Vq} }^{\left(2 \right)}} \right)$;

3) Update in $\zeta $-direction for $\Delta t$ by using ${\rm{CSLR}}\left({{{\overline {{}^Pq} }^{\left(2 \right)}},{{\overline {{}^{L\zeta }q} }^{\left(2 \right)}},{U_\zeta },{{\overline {{}^Pq} }^{\left(3 \right)}},{{\overline {{}^{L\zeta }q} }^{\left(3 \right)}}} \right)$ and ${\rm{CSLR}}\left({{\overline {{}^{L{\zeta _{\bot} }}q} }^{\left(2 \right)}},{{\overline {{}^Vq} }^{\left(2 \right)}}, $$ {U_\zeta },{{\overline {{}^{L{\zeta _{\bot} }}q} }^{\left(3 \right)}},{{\overline {{}^Vq} }^{\left(3 \right)}} \right)$;

4) Update in $\eta $-direction for another ${{\Delta t} / 2}$ as in step 2;

5) Update in $\xi $-direction for ${{\Delta t} / 2}$ as in step 1.

Note that the cubed-sphere grid is not continuous across the cube patch boundaries, and some special treatments are needed. In Eq. (14), if the first guess point $\xi _{i - \frac{1}{2}}^ * $ moves across the cube patch boundary, we first calculate the time that the arrival point reaches the cube patch boundary:

where ${x_{\rm{b}}}$ is the coordinate of cube patch boundary. Then, the departure point is calculated by:

where ${u_{\rm{b}}}$ is the velocity at ${x_{\rm{b}}}$. At cube patch boundaries, the ‘source term’ is calculated by one-side difference instead of the central difference in Eq. (15).

We should note that the flux across the cube patch boundaries is calculated only once in this study. As shown in Fig. 6, A is the arrival point on a 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 the A_d−A line segment. If A_d is on Patch 1, the same process is executed for Patch 1, and so on for each patch boundary.

Figure 6.  Illustration of departure points along patch boundary.

To verify the performance of the proposed transport model, several widely used benchmark tests, including solid body rotation, moving vortices, and deformational flow tests are performed on the spherical mesh.

The normalized errors proposed by Williamson et al. (1992) are used:

where $\Omega $ is the whole computational domain and $q$ and ${q_t}$ refer to numerical solutions (volume-integrated average in our paper) and exact solutions, respectively.

The solid-body rotation test (Williamson et al., 1992) is widely used in two-dimensional spherical transport modeling to evaluate the performance of a transport model. The wind components in the latitude-longitude coordinates $ \left(\lambda,\theta \right) $ are defined as:

where $\left({{u_{\rm{s}}},{v_{\rm{s}}}} \right)$ is the velocity vector, ${u_0} = {{2\pi R} / {1036800}}$ (1036800 s equals 12 days), which means it takes 12 days to complete a full revolution on the sphere, R is the radius of the sphere, and $ \alpha $ 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.

The initial condition of a cosine bell test is specified as:

where ${r_{\rm{d}}}$ is the great circle distance between $\left({\lambda,\,\theta } \right)$ and the center of the cosine bell, located at $\left({{{3\pi }/ 2},\,0} \right)$, ${r_0} = {{7\pi R} / {64}}$ is the radius of the cosine bell, and ${h_0} = 1$.

The normalized errors on $32 \times 32 \times 6$ meshes and with 256 time steps compared with other existing published semi-Lagrangian schemes, the PPM-M scheme (Zerroukat et al., 2007) and CSLAM-M (Lauritzen et al., 2010), are presented in Table 1. The result shows that CSLR1 and CSLR1-M get almost the same result. And our scheme is comparable to the PPM-M scheme, and the result in the near-pole flow direction ($\alpha = {\pi / 2}$ and $\alpha = {\pi / 2} - 0.05$) is better than the CSLAM-M scheme.

To check the influence of the weak singularities at the eight vertices of the cubed-sphere gird, this test is conducted with $\alpha = {\pi / 4}$ to pass through four vertices. The history of normalized errors (CSLR1 and CSLR1-M are almost the same, so we only present the result of CSLR1-M here) are shown in Fig. 7. We can see that the normalized errors have little fluctuations (except the ${l_\infty }$ errors at around day 4 and day 10) when the flow passes four weak singularities.

Figure 7.  History of normalized errors of the solid body rotation of a cosine bell for one revolution on grid N = 32 (number of cells in one direction on each cell), 256 time steps and with α = π/4.

To demonstrate the ability of the CSLR1-M scheme using a large Courant number to transport, we use 72 time steps (local maximum Courant number is about 1.78) with rotation angle $\alpha = {\pi / 2}$ to complete one revolution. The normalized errors are ${l_1} = 0.052$, ${l_2} = 0.046$, and ${l_\infty } = 0.061$.

A non-smooth step cylinder is calculated to evaluate the non-oscillatory property. The initial distribution is specified as

where ${r_{\rm{d}}}$ is the great circle distance between $\left({\lambda,\theta } \right)$ and $\left({{{3\pi } / 2},0} \right)$, which is the center of the step cylinder, ${r_1} = {2 / 3}R$ and ${r_2} = {1 / 3}R$.

In this test, we set $\alpha = {\pi / 4}$, which is the most challenging case of the rotation test where the step cylinder moves through four vertices and along two boundary edges of the cubed-sphere grid to complete a full revolution. Here, we use 90×90×6 meshes and 720 time steps to conduct this test. The numerical results after 12 days are shown in Fig. 8, and we can see that the CSL2 scheme will generate obvious oscillations around the discontinuities. By using the CLSR1 and CSLR1-M approaches, these nonphysical oscillations are effectively removed. The maximum and minimum value of CSL2 are ${q_{\max }} = 1034.23$ and ${q_{\min }} = - 2.45$, and for CLSR1 and CSLR1-M they are ${q_{\max }} = 1001.85$ and ${q_{\min }} = 0$. The history of relative mass errors is given in Fig. 9, which shows that the relative mass errors are up to the tolerance of machine precision, therefore the proposed global transport model is exactly mass conservative during the simulation procedure.

Figure 8.  Numerical results of solid body rotation of the step cylinder after one revolution (12 days). (a) is the result of CSL2, (b) is the result of CSLR1, and (c) is the result of CSLR1-M.

Figure 9.  The time history of relative mass error for solid body rotation of the step cylinder test case by the CSLR1-M scheme.

The second benchmark test we used is the moving vortices test proposed by Nair and Jablonowski (2008). The wind component of this test is a combination of the solid body rotation test and two vortices, and it is much more complicated than the solid body rotation test. The velocity fields on the sphere are specified as:

where ${u_{\rm{s}}}$ and ${v_{\rm{s}}}$ are calculated by Eqs. (37) and (38), and the rotation angle of this test is set to be $\alpha = {\pi / 4}$. ${\rho _0} = 3$, ${\lambda _{\rm{c}}}\left(t \right)$ and ${\theta _{\rm{c}}}\left(t \right)$ are the center of the moving vortex at time $ t $, and the calculation procedure of ${\lambda _{\rm{c}}}\left(t \right)$ and ${\theta _{\rm{c}}}\left(t \right)$ can be found in (Nair and Jablonowski, 2008).

The tracer field is defined as:

where $\gamma $ is a parameter to control the smoothness of the tracer field, $\left({\lambda ',\theta '} \right)$ is the rotated spherical coordinates, which can be calculated by:

and $\left({{\lambda _{\rm{p}}},\,{\theta _{\rm{p}}}} \right) = \left({\pi,\,{\pi / 2} - \alpha } \right)$ is the North Pole of the rotated spherical coordinate. In this test, we followed Norman and Nair (2018) to set $\gamma = {10^{ - 2}}$ to conduct a large gradient in tracer distribution to check the non-oscillatory property and the performance of positivity preserving. When $t = 0$ in Eq. (46), we get the initial condition.

This test is conducted on 80 × 80 × 6 meshes and uses 400 time steps to move forward 12 days. The contour plots in Fig. 10 show that compared with the exact solution, our proposed scheme can simulate this complicated procedure well. The plot along the equator is presented in Fig. 11, and it shows that there are no obvious oscillations around large gradients. The normalized errors of CSLR1 and CSLR1-M are almost the same, being ${l_1} = 5.295 \times {10^{ - 2}}$, ${l_2} = 0.1295$, and ${l_\infty } = 0.5667$, respectively. The histories of minimum values are shown in Fig. 12, where we can see that the CSLR1 scheme would produce negative values during the simulation procedure, while the minimum values of CSLR1-M can completely preserve positivity (the minimum values are within the machine precision).

Figure 10.  Contour plot of moving vortices after 12 days. (a) is the exact solution, (b) is the result of the CSLR1-M scheme.

Figure 11.  Plot along the equator for the moving vortices test at 12 days.

Figure 12.  The histories of minimum values q_min of the moving vortices test. (a) is the result of the CSLR1 scheme, (b) is the result of the CSLR1-M scheme.

The last benchmark test used in our paper is the deformational flow test proposed by Nair and Lauritzen (2010), which is the most challenging test case. The nondivergent and time-dependent flow fields are defined as:

where $\kappa = 2$, $T = 5$, and $\lambda ' = \lambda - \left({{{2\pi t} /T}} \right)$.

Two kinds of initial conditions are checked here, including the twin slotted cylinders case to evaluate the positivity preserving property and correlated cosine bells to evaluate the nonlinear correlations between tracers (Lauritzen and Thuburn, 2012). By the given flow fields, the initial distributions will be deformed into thin bars during the first half period, then return to its initial state during the second half period.

The initial condition is defined as:

where ${r_0} = 0.5$ and ${r_i}\left({i = 1,2} \right)$ represent the great circle distances between the center of the two slotted cylinders and a given point. The centers of the two slotted cylinders are located at $\left({{\lambda _1},{\theta _1}} \right) = \left({{{5\pi } / 6},0} \right)$ and $\left({{\lambda _2},{\theta _2}} \right) = \left({{{7\pi } /6},0} \right)$, respectively.

The numerical results of deformational flow of the CSLR1-M scheme with 90 × 90 × 6 meshes and with 390 time steps (local maximum Courant number is about 3) are shown in Fig. 13. As shown in Fig. 13b, the two slotted cylinders are deformed into two thin filaments by the background flow field during the first half period. Figure 13c gives the counters of the slotted cylinders at the final time, and it is indicated that the proposed scheme can correctly reproduce this complicated deformational flow and does not produce oscillations. The histories of minimum values are shown in Fig. 14, which indicates that the CSLR1 scheme would produce negative values, while the CSLR1-M scheme keeps minimum values within the tolerance of machine precision, which can be viewed as non-negativity. The Normalized errors are ${l_1} = 0.3287$, ${l_2} = 0.3321$, and ${l_\infty } = 0.9415$ for both the CSLR1 and CSLR1-M schemes.

Figure 13.  Numerical result of deformational flow of slotted cylinder after one period by the CSLR1-M scheme. (a) is the exact solution, (b) is the numerical solution at half cycle, (c) is the numerical solution after one cycle.

Figure 14.  Histories of minimum values q_min of deformation of the twin slotted cylinder test case. (a) is the result of the CSLR1 scheme, (b) is the result of the modified CSLR1-M scheme.

To check the ability of preserving nonlinearly correlated relations between two tracers, we used two kinds of tracers. One is the quasi-smooth twin cosine bells:

where ${h_i} = \dfrac{1}{2}\left[ {1 + \cos \left({\dfrac{{\pi {r_i}}}{{{r_0}}}} \right)} \right]$ for $i = 1,2$.

The other one is the correlated cosine bells:

where $ \psi \left(q \right) = - 0.8{q^2} + 0.9$.

This test is conducted on $90 \times 90 \times 6$ meshes with 1800 time steps. The scatter plot of numerical result at $t = {T / 2}$ is shown in Fig. 15. The solution of cosine bells is in the x-direction, and the correlated cosine bells is in the y-direction. The mixing diagnostics are ${l_r} = 1.05 \times {10^{ - 3}}$, ${l_u} = 2.40 \times {10^{ - 5}}$, and ${l_0} = 5.57 \times {10^{ - 4}}$, respectively (see Lauritzen and Thuburn, 2012) for the detail definition of these three parameters). The CSLR1-M scheme is built using a monotone rational polynomial with modest accuracy, which always overly flattens the maximum and minimum values, as shown in the bottom-right corner of Fig. 15. In the whole, the scattering points of the CSLR1-M scheme are almost located inside the convex hull which means that the CSLR1-M scheme can preserve nonlinearly correlated relations between tracers well.

Figure 15.  Scatter plot of nonlinearly correlated cosine bell at t = T/2.

In this paper, a non-negativity and conservative semi-Lagrangian transport scheme based on a multi-moment finite volume method has been developed on the cubed-sphere grid. By using the PV moment and VIA moment, a rational function is constructed as a spatial approximation function within a single cell instead of the non-monotonic CSL2 scheme to suppress the numerical oscillations and keep the monotonicity. In terms of multi-moment concepts, the VIA moment is cast by utilizing the flux form formulation to guarantee the exact numerical conservation. In the CSLR1 scheme, the semi-Lagrangian method is adopted to update the PV moments, which keeps good properties of the semi-Lagrangian scheme. To simplify the implementation in curvilinear (cubed-sphere) geometry, a dimension-splitting time stepping strategy is combined with the multi-moment finite volume method. In the case of a valley of the transported field, two kinds of modifications are conducted on the original CSLR1 scheme for exactly positive-definite preservation. Note that the improved CSLR1-M scheme dose not degrade the accuracy of the original CSLR1 scheme. The numerical results show that the CSLR1-M scheme is non-oscillatory and can preserve the non-linear correlations between tracers. In addition, the semi-Lagrangian approach permits a large time step, which can greatly improve computational efficiency. The quality of the present transport modelling has been demonstrated by the widely used benchmark tests on a cubed-sphere grid. The results reveal that the developed transport modelling not only can effectively remove nonphysical oscillations, but it can also preserve the non-negativity of numerical solutions, which indicates that it has the potential to simulate the various tracers accurately for real applications. When the mass sources such as evaporation, condensation, etc. are involved in real simulation, they can be added to each PV variable through a fractional step in the multi-moment model after the tracers are advected by the CSLR scheme. Furthermore, a positivity constraint is also imposed on the source terms.

Acknowledgments. This work was supported by the National Key Research and Development Program of China (Grant Nos. 2017YFC1501901 and 2017YFA0603901) and the Beijing Natural Science Foundation (Grant No. JQ18001).

Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.