Conservative Semi-lagrangian Scheme Applied to One-dimensional Shallow Water Equations

A forward-in-time semi-Lagrangian scheme developed by Sun et al. (1996) and Sun and Yeh (1997) has been applied to one-dimensional shallow water equations in both rotational and irrotational systems. After obtaining numerical results, we employ variation formulations (Sun and Sun 2004) with minimum correction to adjust both total mass and total energy so that they are conserved. Therefore, the scheme produces accurate, positive-definite solutions while conserving both mass and total energy. Comparing among different resolutions, the improvement on total energy is significant but less significant for a mass field in a coarse resolution model when it simulates the sharp discontinuities of surface waves, because the mass field calculation is quite accurate even without correction. The variation method proposed here can also be easily applied to multi-dimensional flows.


INTRODUCTION
Since 1959, semi-Lagrangian schemes have been studied by Wiin-Nielsen (1959) and many other scientists.Most semi-Lagrangian schemes employing backward trajectories have been reviewed by Staniforth and Côté (1991).Iterations are often required to solve the departure points of the backward trajectories if the velocity field is not constant, as discussed by Bates (1984), Kuo and Williams (1990), Bermejo and Staniforth (1992), and Huang (1994).On the other hand, forward trajectories can easily be obtained with great accuracy (Purser and Leslie 1994).Thus, Purser and Leslie (1991) have introduced the 'cascade interpolation method'.Sun et al. (1996) have presented the 'split interpolation method with clipping filter', and Sun and Yeh (1997) have developed a more general 'internet interpolation method' with unconditional stability, which has also been successfully applied to the NTU-Purdue nonhydrostatic atmospheric model (Hsu and Sun 2001;Sun and Hsu 2005) by Shieh et al. (2006) The use of semi-Lagrangian (SL) advection methods allows for relatively long time steps.It gives minimal phase error, minimizes the computational dispersion, and can handle sharp discontinuities (Nair and Machenhauer 2002).However, semi-Lagrangian does not guarantee the conservation of mass (Rancic 1992;Priestley 1993), which is important for a numerical model to simulate passive-scalar transport in the atmosphere or the properties of climate studies and ocean circulation (Priestley 1993).Hence, Rancic (1992), Priestley (1993), Gravel and Staniforth (1994), Lin and Rood (1996), Bermejo and Conde (2002), and others have developed different mass conserving, semi-Lagrangian schemes.The detailed mass conserving, semi-Lagrangian scheme derived from the piecewise parabolic method (PPM) based on the finite-volume scheme is referred to Rancic (1992).Priestley (1993) and Gravel and Staniforth (1994) applied a localized correction on the quasi-monotone semi-Lagrangian (QMSL) scheme of Bermejo and Staniforth (1992) to achieve mass conservation for a shallow water equation.On the other hand, Bermejo and Conde (2002) applied the variation formulation to achieve a global conservation of the QMSL.Based on Sun et al. (1996) and Sun and Yeh (1997), Sun and Sun (2004) (will be referred as SS hereafter) also applied the variation formulation with mass correction to obtain accurate, positive-definite, mass conserved solutions with simpler calculations.Furthermore, SS's results are as good as or more accurate than those obtained by Carpenter et al. (1990), Priestley (1993), and Bermejo and Conde's (2002) compared with analytical solutions for scalar, linear advection.Sun et al. (1996) have also proved that the accuracy of their semi-Lagrangian scheme is comparable or better than Kuo and Williams (1990) and others in solving nonlinear advection.Here, we apply the method developed by SS to one-dimensional shallow water equations.The simulations show that the scheme is capable of conserving both mass and total energy in rotational and irrotational systems.We are also applying this method to the two-dimensional shallow water equations to study interactions of vortices over topography based on the characteristic approach (Hirsch 1988(Hirsch , 1990;;Wang and Yeh 2005) and the cases presented by Wang and Yeh (2005).They will be presented in other papers.It is also noted that Arakawa and Lamb (1981) and Thuburn and Staniforth (2004) have developed an elegant finite difference scheme to conserve mass, energy, and absolute potential vorticity in shallow water equations.

Basic Equations and Semi-Lagrangian Scheme
The one-dimensional shallow water equations with the Coriolis force are: where h is the height.
Following Stoker (1957) and Yih (1979) and introducing c = (gh) 1/2 , we can replace (2.4) and (2.5) by (2.10) or in a Lagrangian form, (2.12) and from (2.6), we obtain: (2.13) Inspection of (2.11) -(2.12) reveals that the quantity u ± 2c is conserved along dx dt u c / = ± , or they mean that "message" u ± 2c is propagated through the fluid with speed ±c as well as carried with speed u (Yih 1979).On the other hand, v is conserved and carried by speed u alone.
Let us define the solutions obtained from (2.11) -(2.12) as p u c * * , and from (2.13) as v *, where c * = (gh *) 1/2 .Then, u * and h * (before including the effect of rotation) can be determined as follows (Yih 1979;Erbes 1993): (2.15) Then, apply u * and v * in equations (2.7) -(2.8), we obtain the analytic solution: (2.17) According to Sun et al. (1996) and Sun and Yeh (1997), the procedure to solve (2.11) -(2.13) includes: (a) Constructing the Lagrangian network induced by the motion of the fluid from the Eulerian network and finding the intersections of the networks by a general interpolation from the irregularly distributed Lagrangian grid to the regularly distributed Eulerian grid; (b) Applying the spatial filter to remove the unwanted short waves and the values beyond the constraints.

Conservation of Mass
This section follows the procedures discussed in SS closely.Suppose the height at time t n is h j n , then the summation of h j n over the entire domain of jm grids is: where H n is the total mass at n th -time step, and δx j is the space interval at j th grid and can be a variable.The summation of h * in (2.15) is: where the referenced values h max and h min are the extreme values of h j n and can be functions of time and space also, as discussed in SS, we propose a mass correction function, δ h j , at the j th grid, given by the polynomial: is the net mass flux between t n t = δ and t n t = + ( ) 1 δ , which is ignored in this paper.Equation (2.21) guarantees the conservation of mass.Equation (2.22) implies a minimal modification of the results obtained by the original Semi-Lagrangian scheme.Equation (2.20) is one of the simplest equations with two undetermined parameters required to satisfy both (2.21) and (2.22).It is noted that (2.20) approaches zero when h j * approaches h max or h min to ensure no overshooting of the results, which cannot be achieved by a lower-order polynomial.From (2.21) and (2.22), we can obtain: Hence, after calculating b from (2.24), and a from (2.23), we obtain from δ h j (2.20).Thus, the height of the j-th grid at time t n+1 is: which guarantees that: (2.28) Details are discussed in SS.SS obtained almost identical results when the correction function consists of two sine modes: where α π = − /( ) h h max min , a and b are constants to be determined.Hence, only (2.20) is used in this paper.

Conservation of Total Energy
If we multiply (2.1) by gh, (2.2) by u and (2.3) by v, we obtain: (2.30) , where h b is the basic height and h' is the perturbation, and integrate (2.30) over the entire domain between n t ∆ and ( ) n t +1 ∆ , we obtain: (2.31a) For a closed domain, or domain with a periodic boundary or far away boundaries, we obtain: δ is a constant according to the conservation of mass.Let us define the total Kinetic Energy, K h u v x j j j δ ; and the total potential energy, P g h The total energy at the n th time step, E K P n n n = + , should be conserved and equal to its initial value E 0 .
After solving û j , v j from (2.16) -(2.17) and h j n +1 from (2.27), we will apply the variation formulation so that the total energy can be conserved with a minimum modification to those obtained by (2.16) and (2.17).If we define: the perturbation of δκ ˆj for a close or periodic system should satisfy the following constraint: Here δκ ˆj is similar to δ h j in (2.20).We can follow the variation formulation (2.20 -2.22) to calculate δκ ˆj and the kinetic energy: (2.34) Then, we assume that the change of kinetic energy does not affect the wind direction, that is: The variation technique discussed here can be easily applied to multi-dimensional problems, as shown in SS, which conserves the total mass in two-dimensional flow.

NUMERICAL MODEL
The model based on A-grids consists of jm grids with a uniform space interval δx , which is integrated for N-time steps.The basic height h b = 200 m.The initial perturbed height h' is given by: where h a ' is the amplitude of height perturbation, x j and x c are the location at j th -grid and the central point, respectively; λ = = − ( ) / .
. gh b 1 0 447 2 4 km, is the radius of deformation if f = 10 -4 s -1 in a rotational fluid and g = 10 m s -2 ; where A is a parameter to control the length scale Lx of the initial disturbance, Lx A = λ .We also calculate the error of the total mass (Em) and error of the total energy (Ee) at the N th time step (last time step) as follows: Table 1 provides the parameters used in this study.The model is integrated with double precision.

NUMERICAL RESULTS IN AN IRROTATIONAL SYSTEM
In addition to pure gravity waves (i.e., f = 0), the results will show the numerical simulations of geostrophic adjustment from the initial shortwave (A < 1) and longwave (A > 1) disturbances of either velocity or mass field.

Case 4A
The system starts from rest (u = v = 0 everywhere) and an initial perturbation h' given by (3.1) with h a ' = 0.1 m, h b = 200 m, A = 0.2, δ t = 320 s, δx = 3.125 km, jm = 1601, and f = 0, as indicated by curve A in Fig. 1a and Case 4A-a in Table 1.As time increases, the perturbation propagates as surface gravity waves, as shown by B, C, and D at t = 1.07 × 10 4 , 2.13 × 10 4 , and 3.2 × 10 4 s, respectively.Because the perturbation is very small, the simulation is almost identical to the analytic solution.The waves propagate with a phase of 44.72 m s -1 (= gh ).The maximum Courant number CN, which is define as CN = (max ) / u gh t x + δ δ , = 4.58 with max .u ≈ 0 012 m s -1 .The total kinetic energy K is identical to the total potential energy P, as shown in Fig. 1b.This is called the principle of equipartition of energy and is valid in conservative dynamical systems undergoing small oscillations that are unaffected by planetary rotation (Kundu 1990).The error in mass field is very small in our simulations, as shown in Table 1.
Without applying the variation principle, the error of mass, Em = 0.5696 × 10 -11 and error of the total energy, Ee = -.2791× 10 -6 are very small for Case 4A-a shown in Table 1, because the scheme is very accurate for advection equations (Sun et al. 1996;Sun and Yeh 1997;Sun and Sun 2004;Shieh et al. 2006).
Case 4A-b has the same initial perturbation as Case 4A-a, except for a coarse resolution δx = 12.5 km and jm = 401.The results of Case 4A-b are very close to Case 4A-a except for the Courant number = 1.145, and slightly larger errors, Em = -0.90× 10 -11 and Ee = -0.15× 10 -3 , as shown in Table 1.

Case 4B
Case 4B-a is the same as Case 4A-a except h a ' = 20 m, as shown in Table 1 and Fig. 2. Because the height is always symmetric and velocity is asymmetric with respect to the central point, x c in our simulations, hereafter, we will show the results for x x c ≥ =0 only.The Courant number = 4.915 because of a stronger velocity ( max .u ≈ 2 2 m s -1 ) generated by a larger surface gradient.The phase speed is about 47 m s -1 , which is carried by u + c according to (2.9).The solutions depart from the linear solution as the tilt of both height field and velocity field increases with time.Similar results were obtained by Erbes (1993).Eventually, the wave collapses when the slope approaches infinity (i.e., ∂ ∂ h x / → ∞), which is similar to collapse of the Burger equation discussed in Sun et al. (1996).The error in total mass is still very small, although it is larger than Cases of 4A.The error in total energy is comparable between 4A-a and 4B-a.It is also noted that the ratio of P/E 0 in 4 B-a is slightly less than 0.5.
The solutions of Case 4B-b with coarse resolution δx = 12.5 km and jm = 401 are also included in Table 1.The error in mass is still negligible, Em = -.13 × 10 -5 but the error in the total energy, Ee = -.33 × 10 -2 , which is much larger than the fine resolution simulation of Case 4B-a.On the other hand, Case 4B-c in table 1 shows that both mass and total energy can be conserved when the variation principle is applied.It is also noted that the difference in either height or velocity fields is very small for most points between Cases 4B-b and 4B-c.

RESULTS WITH THE CORIOLIS FORCE
The propagation of waves becomes more complicated in a rotating system due to the existence of the Coriolis force (Rossby 1938;Gill 1982).The dispersion relation in linearized shallow water equations without mean flow becomes: where σ is frequency, f = 1.0 -4 s -1 is the Coriolis parameter, and k is the wave number.The phase speed of the dispersive wave is: (5.2) Gravity (or buoyancy) effect is more important than the Coriolis force for short wave (ck >> f ).
On the other hand, the Coriolis force becomes dominant for a long wave (ck << f).

Case 5A
The results for case 5A-a with f = 1.0 -4 s -1 , h a ' = 20 m, A = 0.2, δ t = 320 s, δx = 3.125 km, and jm = 1601 but without conservations of mass and energy are shown in Table 1 and Fig. 3, where only half of the domain ( x x c ≥ =0) will be presented because h is symmetric; while u and v are asymmetric with respect to x c .Figures 3a -c show the simulated h, u, and v at t = 0, 1.333 × 10 4 , 2.667 × 10 4 , and 4 × 10 4 s, which are indicated by A, B, C, and D, respectively.Figure 3d reveals the time sequences of the total energy, potential energy, and kinetic energy.At the beginning, the unbalanced pressure gradient force in x < λ generates the x-component wind and the fast propagating gravity waves, the amplitude of the waves is smaller than that in Case 4B, because part of pressure gradient is counterbalanced by the y-component velocity, as we can see at x = 1 × 10 2 km, the Coriolis force quickly turns a westerly wind to a northerly wind.Overshooting occurs in x < λ , thus, h' at x = 0 decreases from 20 to 1.8 m within 13000 s, then increases to 3.5 m at t = 40000 s.The leading wave propagates away with a phase speed of 47.5 m s -1 , which is similar to Case 4B but 0.5 m s -1 faster here due to the effect of rotation.The amplitude also gradually decreases with time.
The results reveal that for x < 100 km, the y-component momentum reaches steady state within 20000 s while the height fields still changes with time because Lx is smaller than the radius of deformation λ as shown in Figs.3a -c.
The time sequence of energy (Fig. 3d) reveals that half of the potential energy quickly converts to the kinetic energy within 3000 s.Then the energy conversion between potential energy and kinetic energy becomes much slower due to the inertia-gravity wave oscillation of the large scale waves shown in Figs.3a -c.The amplitude of oscillation decreases slowly with time as the major portion of the system gradually approaches geostrophic balance.
When a coarse resolution, δx = 12.5 km, is used, the simulated h at t = 0, 1.33 × 10 4 , 2.66 × 10 4 , and 4 × 10 4 s and the time evolution of energy are shown in Figs.4a -b, which  show a sharp peak in height field and 2.5% error in total energy at t = 40000 s as shown in Case 5A-b in Table 1.On the other hand, mass and total energy can be conserved when the variation principle is applied, as shown Case 5A-c and 5A-d in Table 1.

Case 5B
First, we introduce an initial perturbation h' as in Case 5A-a, then, a geostrophic wind balance is applied to obtain the y-component wind, which is shown as curve A in Fig. 5a for the right hand side of the domain ( x x c ≥ =0).After introducing y-component wind, we reset to h' = 0 initially, as curve A in Fig. 5b.Hence, there is a y-component wind but no height pertur- bation initially.Furthermore, conservation of mass and total energy is also applied in Case 5B-a.The results are shown in Table 1 and Figs.5a -d, in which we can see that the y-component wind does not change much with time within x < 200 km.On the other hand, h' increases rapidly between t = 0 and t = 1.33 × 10 4 s (Fig. 5b).The waves generated by an initial unbal- anced force propagate away as inertia-gravity waves (Figs.5b -c).Because of a little adjustment in the momentum field, about 95% of the total energy remains as kinetic energy, and only 5% of the initial kinetic energy is converted to potential energy (Fig. 5d).Without applying the variation principle to control the total mass and energy, we obtain 3% error in total energy, as shown in Fig. 5e and Case 5B-b in Table 1, although the difference in results between 5B-a and 5B-b small in both velocity and height fields.The error in total mass is also small as shown in Table 1.
From Cases 5A-B, we can see that the mass field adjusts to the momentum field for an initial short wave perturbation (A << 1) according to geostrophic adjustment (Gill 1982).

Case 5C
The initial condition is the same as Case 5A, except A = 3.6 and δx = 12.5 km and δ t = 1600 s (CN = 6.02).The simulated heights for Case 5C-a (without controlling for mass or energy) at t = 0, 1.333 × 10 4 , 2.667 × 10 4 , and 4 × 10 4 s are shown in Fig. 6a, and the time sequence of velocity vector (u, v) at x = 500, 1300, and 3825 km in Fig. 6b .Results show that the velocity adjusts to the mass field while around 90% of energy retains as potential energy, as shown in the time sequence of energy in Fig. 6c.Inertia-gravity oscillations exist in all fields with a period of about 57500 s with a small oscillation in mass (Fig. 6a) but a much larger in velocity field (Fig. 6b), although the waves are rather smooth compared with the disturbances generated by smaller scale perturbations in Cases 5A and 5B.The results remain about the same with δ t = 2666.7 s and CN = 10.05 (referred as Case 5C-b in Table 1).
It is also noted that the errors in mass and total energy remain very small even without applying the variation principle.

Case 5D
The initial condition of Case 5D-a is the same as Case 5B-a, except A = 3.6 here.Although about 75% of the energy remains as kinetic energy (Fig. 7c), the oscillation of both potential and kinetic energies is much larger than for Case 5C-a.This indicates that it is slower to approach a geostrophic balance with an unbalanced large-scale momentum field as initial condition.Table 1 also shows that the errors for Case 5D-a are quite small even without controlling for mass or energy.The errors remain very small with a larger time interval δ t = 2666.7 s as shown in Case 5D-b in Table 1.This paper reveals that we can apply the semi-Lagrangian scheme presented by Sun et al. (1996), Sun and Yeh (1997) to obtain accurate, positive-definite solutions from one-dimensional shallow water equations with or without the Coriolis force.According to SS, we can also add the variation formulations to guarantee the conservation of mass and total energy of the solutions.
Although in the 1-D case only left or right characteristic directions (i.e., ±x directions) exist, in the multi-dimensional cases there are an infinite number of characteristic directions that could be used for characteristic wave tracking, as discussed in Hirsch (1988Hirsch ( , 1990) ) and Wang and Yeh (2005).Hence, the characteristic equations of the semi-Lagrangian system in 2-D or 3-D are quite complicated.However, the semi-Lagrangian method is much more accurate than the finite difference method when they are applied to simulate flow with a sharp gradient.
After solving u, v, and h, we can easily apply the variation principle to achieve the conservation of mass and total energy in a 2-D or 3-D flow.Simulations for two-dimensional shallow water equations will be presented in future papers.

SUMMARY
A forward-in-time semi-Lagrangian scheme developed by Sun et al. (1996) and Sun and Yeh (1997) has been applied to one-dimensional shallow water equations in both rotational and irrotational systems.After obtaining the numerical results, we employ the variation principle (Sun and Sun 2004) to find the minimum corrections needed to adjust both height and total energy so that they are preserved.Because the minimum correction based on variation principle is applied, diffusion or dispersion introduced here is rather localized and is smaller than that applied by Bermejo and Conde (2002), as discussed in SS.The scheme produces accurate, positive-definite solutions while conserving both mass and total energy.The variation principles can significantly reduce the error in the total energy for coarse-resolution models; however, the improvement is less significant in the mass field because the original schemes produce quite accurate mass even without any adjustment.It is noted that most people are concerned with the conservation of mass, because total energy decays in a viscous flow.Whilst, an accurate calculation of the total energy should be important, it is more difficult than the mass field.
The simulations show that (non-dispersive) gravity waves are generated by an initial perturbation in an irrotational system.On the other hand, the flow in a rotational system goes through geostrophic adjustment; the waves generated by unbalanced forcing propagate away as inertia-gravity waves.Currently we are working on two-dimensional flows to study the interactions of vortices over topography and the problems discussed by Wang and Yeh (2005) based on the characteristic-based semi-Lagrangian method presented here.The characteristic equations and semi-Lagrangian scheme for multi-dimensional flows are quite complicated and will be presented in future papers.
under Grant Number 93-2745-P492 of Taiwan and National Center for High Performance Computing.

Table 1 .
Summary of the simulations for different cases with f