Space-time transformation in flux-form semi-Lagrangian schemes

With a finite volume approach, a flux-form semi-Lagrangian (TFSL) scheme with space-time transformation was developed to provide stable and accurate algorithm in solving the advection-diffusion equation. Different from the existing fluxform semi-Lagrangian schemes, the temporal integration of the flux from the present to the next time step is transformed into a spatial integration of the flux at the side of a grid cell (space) for the present time step using the characteristic-line concept. The TFSL scheme not only keeps the good features of the semi-Lagrangian schemes (no Courant number limitation), but also has higher accuracy (of a second order in both time and space). The capability of the TFSL scheme is demonstrated by the simulation of the equatorial Rossby-soliton propagation. Computational stability and high accuracy makes this scheme useful in ocean modeling, computational fluid dynamics, and numerical weather prediction.


InTroducTIon
From a physical point of view, advection of a passive tracer is the simple transition of a quantity without diffusion and dispersion. Numerical approaches in atmospheric and oceanic modeling inevitably introduce diffusion (or dissipation) and dispersion into the approximate solution. The numerical diffusion and dispersion are aliens to the process that is being modeled Fan 1998, 1999). As applied to a constituent advection problem, these numerical artifacts manifest themselves as nonphysical mixing by numerical diffusion, nonphysical highs and lows in the constituent field caused by dispersion, and nonphysical tracer spectra caused by trapping in nonpropagating small spatial scales (Rood 1987). For example, the commonly used upwind scheme is conditionally stable (with the Courant number being much smaller than 1) and some artificial viscosity is introduced. Hence, less the numerical diffusion and dispersion errors equates to better model performance.
Many numerical algorithms have been proposed to reduce numerical diffusion and dispersion errors and to keep the numerical stability. The flux-form semi-Lagrangian scheme is among them. Using the flux-form semi-Lagrangian schemes, artificial viscosity is reduced and stability is kept without the limitation of a Courant number (Casulli 1990(Casulli , 1999. In this study, we use a finite volume approach to develop time-space transformed flux-form semi-Lagrangian (TFSL) scheme. This scheme has an explicit form and much less diffusion and dispersion errors.
The stability and accuracy of numerical schemes for ocean models are usually verified using the propagation of a Rossby soliton on an equatorial beta-plane. In principle, the soliton propagates to the west at a fixed phase speed, without a change of shape. Since the uniform propagation and shape preservation of the soliton are achieved through a delicate balance between linear wave dynamics and nonlinearity. In other words, the Rossby soliton is non-diffusive and non-dispersive (Boyd 1980), which makes it a perfect test case for verification of numerical schemes in ocean models since any diffusion and dispersion in the numerical solution of the Rossby soliton are computational errors. Interested readers are referred to the website: http://marine. rutgers.edu/po/index.php?model=test-problems. Terr. Atmos. Ocean. Sci., Vol. 21, No. 1, 17-26, February 2010 To show the benefit of using the TFSL scheme, we first show instability and large diffusion and dispersion errors in numerical solution of the Rossby soliton using the existing schemes such as the flux-form upwind, flux-form central, Lax-Wendroff, and flux-form semi-Lagrangian schemes. Then, we will describe the procedures of the TSFL scheme development and verification. The rest of paper is organized as follows. Section 2 describes the equatorial Rossby soliton and its usefulness for an ocean model verification. Section 3 shows the failure of the three existing schemes (upwind, central, Lax-Wendroff, semi-Largangian) in simulating the equatorial Rossby soliton. Section 4 introduces the TFSL scheme. Section 5 derives the analytical form of the amplification factor of the TFSL-scheme. Section 6 shows the capability of the TFSL-scheme in simulating the equatorial Rossby soliton. Section 7 presents our conclusions.

roSSby SoLITon
Let X be the angular frequency of earth's rotation and R be the earth radius, and let (x, y) be the spatial coordinates with unit vectors (i, j) and t be the time. Consider a single layer of homogeneous ocean layer with depth of H. Lamb's parameter is defined by where g is the gravitational acceleration. The length and time are non-dimensionalized by For the mean ocean depth H = 4 km, the earth radius R = 6370 km, and X = 2π / (86400 s), the length and time scale are L = 543 km, T = 22.39 hr. Let (x, y) be the non-dimensional Cartesian coordinates, (u, v) be the non-dimensional velocity components in the meridional and latitudinal directions, and φ be the non-dimensional surface elevation. After defining and transforming the nonlinear shallow water wave equations into a frame of reference moving with the linear wave, the flow variables (u, v, φ) for the mode-1 wave can be represented by (Boyd 1980) ( , , ) where S is treated as a source/sink term. Evidently Eq. (7) has the analytical solution Eq. (6). Since the analytical solution Eq. (6) exists, the Rossby soliton Eq. (7) is a perfect test case for verifying the stability and accuracy of numerical schemes since the diffusion term has been changed into the given source/sink term. To do so, the fluid is assumed to occupy the equatorial region surrounding the earth. The zonal direction is discretized into 120 cells (i.e., resolution at 3° longitude). The increment Δs is given by The independent variables (s, t) are discretized by s i = s i -1 + Δs, t n = t n -1 + Δt, i = 1, 2, ....; n = 1, 2, ...., with Δt the time step. The dependent variable (η) at (s i , t n ) is represented by

SeverAL exISTIng SchemeS
Equation (7) can be discretized using the flux-form upwind scheme, and the semi-Lagrangian scheme, In order to compare the difference between numerical and exact solutions (a westward propagating Rossby soliton), the zonal equatorial strip is assumed to be infinitely long. When the Rossby soliton travels over n × 120 cells, it goes around the earth once n times (called n cycles). The exact solution at t = 0 is taken as the initial condition, with s = 0 denoting 0° longitude. Three difference Eqs. (10) -(12) are solved numerically from the initial condition Eq. (15) representing the upwind, central, and Lax-Wendroff schemes (Lax and Wendroff 1960) with varying Δt at each time step for a given Courant is due to the fact that the proposed TFSL scheme will be reduced to the Lax-Wendroff scheme for C ≤ 0.5 [see Eq. (36)]. After obtaining the numerical solution, ( , ) x t i n h , substituting it into Eq. (4c) yields φ(s i , y j , t n ). The accuracy of the schemes can be verified through their capability in predicting the westward propagation of the Rossby soliton. To do so, the surface elevation φ(s i , y j , t n ) is plotted with contour values of 2. 13, 4.26, 6.4, 8.53, 10.66, 12.79, 14.93, and 17.06 cm. All the numerical schemes greatly distort the Rossby soliton (Fig. 1). When the ratio of the root-mean square error versus the root-mean of the analytical solution is greater than 100%, the numerical solution is considered divergence. Figure 1 shows that the numerical solution diverges at 27°45'W using the flux-form central scheme, at 54°45'W using the Lax-Wendroff scheme, and at 30°W after one cycle around the earth using the flux-form upwind scheme. Comparing Figs.1b -d to Fig. 1a, the numerical solutions are totally different from the analytical solution.

Semi-Lagrangian method
Consider the advection of a passive scalar φ(x, t) by the velocity u(x, t). The Eulerian formulation is given by where x is the position vector, D/Dt denotes the material derivative, while the Lagrangian counterpart is where the subscript 'p' shows the fluid particle in Lagrangian sense. Although Eqs. (16) and (17) carry the same physical information, their discretization and numerical implementation is different: Eq. (16) is discretized on an Eulerian grid with a finite number of grid points and then time-advanced, while Eq. (17) is integrated for a finite number of fluid particles. Semi-Lagrangian methods combine both Eulerian and Lagrangian points of view; the scalar field is discretized on an Eulerian grid, but is advanced in time using Eq. (17). The key element in accomplishing this is the identification of each grid point x i as the arrival point, for instance, at t + ∆t, of a particle originating from x * i at time t. The algorithm has three steps: (a) The particle associated with each grid point x i at time t + ∆t is traced back to its location x * i at time t, is obtained by interpolating the known values at neighboring grid points, where P is any interpolation operator and ( ) x ik t denotes the set of interpolation points associated with x * i , for example, the nodes of the cell containing x * i ; (c) Finally, the scalar is updated, Thus, the main issues of the semi-Lagrangian method are the backward integration in step (a) and the interpolation in step (b).

Flux Form
Equation (16) can be rewritten in the flux form with inclusion of diffusion, we obtain the finite difference equation of the flux-averaged transport, (a) i j k i j k n n i j k n n 1 2 1 1 2 1 where (F, G, H) are components of the vector F, and represents the temporal average (from t n to t n + 1 ). The tilde represents the volume average over ijk X , The hat represents the combined volume ( ijk X ) and temporal average (from t n to t n + 1 ), For the finite volume ijk DX , the flux at x = x i -1/2 and t = t n is calculated by To solve Eq. (22) numerically, we need to compute the temporally integrated fluxes, , , , If these fluxes are computed using the semi-Lagrangian method, it is called the flux-form semi-Lagrangian scheme (Casulli 1990(Casulli , 1999Lin and Rood 1996).

Transformation of Temporal Integration into Spatial mean
For simplicity and no loss of generality, we consider one dimensional problem of Eq. (22) without source/sink term (i.e., S ijk t = 0), From the semi-Lagrangain consideration, we have In the existing flux-form semi-Lagrangian schemes, the temporally integrated flux F / ( ) , i n n 1 2 1 -+ [similar for F / ( ) , i n n 1 2 1 + + ] is given by the mean value at the two time steps t n and t n + 1 (e.g., Casulli 1990;Tanguay et al. 1990), Using the characteristic-line concept, the flux at time step t n + 1 and location x i -1/2 can be transformed into the flux at time step t n and location x i -1/2 -C (Fig. 2), The bracket [ ] represents the round-off integer. Similarly, the temporally averaged flux at the right boundary ( (from t n to t n + ∆t) are transformed into the spatially averaged fluxes over multiple grids at time step t n with weights of δ 1/2 , δ 1 , ...., δ m . If the characteristic line at t n is beyond the boundary, the boundary condition can be used to calculate F / ( ) , i n n 1 2 1 -+ (Fig. 3), where F b n is the boundary value between F n 1 and F n 1 1 + , and is interpolated by for C ≤ 1/2. This term can be regarded as the numerical (positive) diffusion which leads to computational stability.
Different schemes have different algorithms to compute the temporally averaged fluxes F / ( ) , i n n 1 2 1 -+ and F / ( ) , i n n 1 2 1 + + (from t n to t n + ∆t). The TFSL scheme has second order accuracy in both time and space.

STAbILITy oF The TFSL Scheme
The stability of numerical schemes is an important issue in solving the advection Eq. (16). In section 3, we showed the instability of the existing schemes (upwind, central, and Lax-Wendroff). To determine the stability of the TFSL scheme Eq. (36), the Fourier series expansion is used. Decay or growth of an amplification factor indicates whether or not the numerical algorithm is stable (von Neumann and Richtmyer 1950). Assuming that at any time step t n , the compute solution i n z is the sum of the exact solution

SImuLATIng The roSSby SoLITon uSIng The TFSL Scheme
The TFSL-scheme Eq. (36) is only for a spatially variant and temporally invariant u. When u [or -f 1 η in Eq. (7)] at x i -1/2 varies with time from t n to t n + 1 , concept of variant characteristic lines can be used to determine u(x i -1/2 , t) with sub time-steps (δt 1/2 , δt 1, ..., δt m ) (between t n and t n + 1 ) from u(x, t n ) at grid points (x i -1 ,..., x im , x i* ), and for u > 0 the time from the left neighboring grid x i -[k + 1] to x i -k is given by (Fig. 5), Equation (7) for the Rossby soliton is discretized using the flux form, where S i t is the temporally-spatially averaged source term After the numerical solution η(x i , t n ) is obtained, substituting it into Eq. (4c) yields φ(s i , y i , t n ) as shown in Fig. 6. Note that the flux-form semi-Lagrangian scheme is highly distorted (Fig. 6) with the numerical solution diverging at 40°30'W, which may be caused by the error accumulation. However, the TFSL-scheme is quite stable and accurate. After propagating westward around the earth the numerical Rossby soliton (using the TFSL scheme) appears to be almost non-diffusive and non-dispersive.
To show the quality of the TFSL-scheme, the difference Eq. (47) is integrated for C = 1.5 for a long time period corresponding to the Rossby soliton propagates westward around the earth 5 times. The solution φ(s i , y i , t n ) is stable all the time (Fig. 7). The relative root-mean-square error is calculated to illustrate the accuracy of the TFSL scheme. Table 1 shows RRMSE at the end of first five cycles around the earth. The error varies from 2.66% for the first cycle to 3.53% for the fifth cycle.

concLuSIonS
(1) This study shows that the TFSL scheme is a promising stable and accurate method for solving the advectiondiffusion equation. The Fourier analysis shows that the TFSL scheme has second-order accuracy in time and space. This scheme retains the good features of semi-Lagrangian schemes (no Courant number limitation) with higher accuracy. Computational stability and higher accuracy than the widely used schemes (central, upwind, Lax-Wendroff, semi-Lagragian) makes this technique useful in ocean modeling, computational fluid dynamics, and numerical weather prediction.
(2) Several major features distinguish the TFSL scheme from existing schemes, both Eulerian and semi-Lagrangian. First, the flux (F) at the side of each grid cell is computed not from a single time step (present or next) but from a temporal integration from the present time step to the next time step. Second, this temporal integration is transformed into a spatial integration at the present time step using the characteristic line method.
(3) The equatorial Rossby soliton is used to test the capability of the TFSL scheme since it has exact solution. The equation is solved numerically from the soliton initially located at the equator and 0° longitude with an overall Courant number of 0.75. The upwind, central, and Lax-Wendroff schemes greatly distort the Rossby soliton and diverge as it propagates. However, the TFSL scheme does not distort the Rossby soliton and converges as it propagates many cycles around the earth. With an overall Courant number of 1.5, the numerical Rossby soliton can still propagate many cycles around the earth using the TFSL scheme, but diverges at 40°30'W using the flux-form semi-Lagrangian scheme.
(4) Application of the TFSL scheme to the atmospheric and oceanic models needs more research. This is because that the highly accurate treatment of the source term enabled good performance of the TFSL scheme. However, this is a very special case because the analytical solution of S (source term) is known for this system. In ocean modeling, computational fluid dynamics, or numerical weather prediction, source term analytical solutions are usually unknown and thus several iteration processes will be required for the departure/arrival point estimation as well as the source term estimation (for spatial and temporal averaging), which causes the loss of efficiency and accuracy.
(5) The TFSL scheme was developed on the basis of a finite volume approach. It is relatively easy to extend one-dimensional space-time transformation Eqs. (28b) and (28c) into three-dimensional transformations. The space integration of the flux is conducted over the twodimensional surface of the finite volume, and the time integration is for that volume (see section 4.2). This will be reported in a separate paper in the near future.