# CIP Method of Characteristics for the Solution of Tide Wave Equations.

1. Introduction

Tide wave is the water movement caused by the gravitational force of the celestial body. The control equations of tidal wave model are the original shallow water equations including source terms (e.g., bottom friction, bottom topography and horizontal eddy viscosity term) . Therefore, getting the solution of shallow water equations with source terms is of great significance.

In the simulation of ocean tidal waves, Eulerian schemes are widely used, for example, Backhaus  and Casulli  used semi-implicit scheme (hereafter SI) for the solution of shallow water equations; Lv and Zhang  used the semi-implicit scheme to solve tide wave equations and their computational format was used to study bottom friction coefficients  and tidal open boundary conditions ; the finite volume method is also frequently used [7, 8]. However, since the phase speed of the gravity waves is so fast, time step is strictly limited by CFL condition when using Eulerian schemes .

There are growing interests in finding methods that enable using large time steps to solve shallow water equations. Erbes  proposed a semi-Lagrangian method of characteristics (MOC) with quadratic interpolation to solve the shallow water equations with a large CFL number, but there was dispersion error produced. Ogata and Yabe  applied the CIP method to shallow waters combined with the MOC and showed low dispersion error and low numerical damping even with large CFL number. Toda et al.  proposed a new scheme by adopting the Conservative Semi-Lagrangian (CSL) based on the CIP method which shows good conservation. However, applying this method in realistic simulation of ocean tidal wave will meet a difficulty treating source terms including bottom topography, bottom friction, and boundary conditions.

This paper is designated to apply the CIP-MOC  to the solution of tide wave equations which are based on shallow water equations and include source terms of bottom topography and friction. The ideal numerical experiment shows that the CIP-MOC method outperforms SI scheme in stability and dispersion errors when using large time step. Problem with reflection boundary condition is also discussed.

This paper is organized as follows. Section 2 presents the CIP-MOC method. In Section 3, numerical experiments are carried out and comparison is made between the CIP and SI. Conclusion is given in Section 4.

2. Model and Method

2.1. CIP Method. CIP method is a compact, robust, less diffusive, and high-order scheme in computational fluid dynamics . Compared with other traditional interpolation methods, the internal information of grid is better described by CIP which uses spatial derivative in the description of function distribution in grid cells. We consider the following advection problem:

[partial derivative]f/[partial derivative]t + u [partial derivative]f/[partial derivative]x = 0, (1)

the solution is 

[f.sup.n+1] = F(x - u[DELTA]t),

[[partial derivative].sub.x][f.sup.n+1] = dF (x - u[DELTA]t)/dx. (2)

If f and [[partial derivative].sub.x]f are known at each grid point, the profile between two adjacent points can be approximated by cubic polynomial

[F.sub.i] (X) = [a.sub.i][X.sup.3] + [b.sub.i][X.sup.2] + [[partial derivative].sub.x][f.sub.1]X + [f.sub.i], (3)

where X = x - [x.sub.i]. Consider four constraints

[mathematical expression not reproducible] (4)

then coefficients can be given as

[a.sub.i] = ([[partial derivative].sub.x][f.sub.i] + [[partial derivative].sub.x] [f.sub.iup])/[D.sup.2] + 2 ([f.sub.i] - [f.sub.iup])/[D.sup.3], (5)

[b.sub.i] = 3 ([f.sub.iup] - [f.sub.i])/[D.sup.2] - (2[[partial derivative].sub.x][f.sub.i] + [[partial derivative].sub.x] [f.sub.iup])/D, (6)

where sgn([u.sub.i]) = 1 ([u.sub.i] [greater than or equal to] 0), -1 ([u.sub.i] < 0), iup = i - sgn([u.sub.i]), and D = -[DELTA]x x sgn([u.sub.i])

Thus, the profile of f and [[partial derivative].sub.x] f at next time step can be obtained as

[f.sup.n+1.sub.i] = [F.sub.i] ([x.sub.i] - [u.sub.i][DELTA]t) = [a.sub.i][[xi].sup.3] + [b.sub.i][[xi].sup.2] + [[partial derivative].sub.x][f.sup.n.sub.i] x [xi] + [f.sup.n.sub.i], (7)

[[partial derivative].sub.x][f.sup.n+1.sub.i] = d[F.sub.i] ([x.sub.i] - [u].sub.i][DELTA]t)/dx = 3[a.sub.i][[xi].sup.2] + 2[b.sub.i][xi] + [[partial derivative].sub.x][f.sup.n.sub.i], (8)

where [xi] = -[u.sub.i][DELTA]t.

2.2. One-Dimensional Tide Wave Equations. Let bottom friction term be the linear form, then the one-dimensional tide wave equations in a primitive form are written as 

[partial derivative][zeta]/[partial derivative]t + [partial derivative] [([h.sub.0] + [zeta]) u]/[partial derivative]x = 0,

[partial derivative]u/[partial derivative]t + u [partial derivative]u/[partial derivative]x + g [partial derivative][zeta]/[partial derivative]x + ku = 0, (9)

where t is time, x is Cartesian coordinate, [h.sub.0](x) is undisturbed water depth, [zeta](x, t) is sea surface elevation above the undisturbed sea level, u(x, t) is the flow velocity, k is the bottom friction coefficient, and g is the gravitational acceleration.

To solve (9) with the method of characteristics and preserve the balance between the depth gradient and the bottom effect, the surface gradient method  is used here. Replace the variable [zeta] with the water level H = [zeta] + [h.sub.0] + z, where z(x) represents the bottom topographic height (see Figure 1), we get the equivalent of (9) which can be written as

[partial derivative]H/[partial derivative]t + [partial derivative](Hu)/[partial derivative]x = [partial derivative](uz)/[partial derivative]x,

[partial derivative]u/[partial derivative]t + u [partial derivative]u/[partial derivative]x + g [partial derivative]H/[partial derivative]x + ku = 0. (10)

In a vector-matrix form, (10) can be rewritten as

[mathematical expression not reproducible] (11)

The eigenvalue matrix [LAMBDA] and eigenvector matrix L of coefficient matrix A are

[mathematical expression not reproducible] (12)

where [[GAMMA].sub.H] = [square root of (gH)] and [[lambda].sup.[+ or -]] = u [+ or -] [[GAMMA].sub.H]. Thus, matrix A can be diagonalized as [L.sup.-1] AL, where

[mathematical expression not reproducible] (13)

is the inverse of L. Then (11) can be written as

[L.sup.-1] [partial derivative]W/[partial derivative]t + ([L.sup.-1] AL) [L.sup.-1] [partial derivative]W/[partial derivative]x = [L.sup.-1] S, (14)

where S represents the right hand side of (11).

Two characteristics about Riemann invariants [R.sup.[+ or -]] = [[GAMMA].sub.H] [+ or -] u/2 can be obtained from (14) as

[partial derivative][R.sup.[+ or -]/[partial derivative]t + [[lambda].sup.[+ or -]] [partial derivative][R.sup.[+ or -]]/[partial derivative]x = g/2[[GAMMA].sub.H] [partial derivative](uz)/[partial derivative]x [- or +] 1/2 ku, (15)

By using CIP method, we have

[mathematical expression not reproducible] (16)

where [[GAMMA].sup.[+ or -].sub.H] and [u.sup.[+ or -]] are calculated by

[[GAMMA].sup.[+ or -].sub.H] = [[GAMMA].sub.HCIP] ([x.sub.i] - [[lambda].sup.[+ or -]] [DELTA]t),

[u.sup.[+ or -]] = [u.sub.CIP] ([x.sub.i] - [[lambda].sup.[+ or -]] [DELTA]t), (17)

the subscript CIP in variables means that they are interpolated with the CIP scheme, and the bottom friction term is discretized with the Crank-Nicholson method which is implicit in time. [[GAMMA].sub.H] and u at time step [t.sup.n+1] can be obtained by solving linear equations (16) (see Figure 2):

[mathematical expression not reproducible] (18)

[mathematical expression not reproducible] (19)

Time integration of the bottom topographic term in (18) and (19) can be approximated using [dx.sup.[+ or -]]/d[[tau].sup.[+ or -]] = u [+ or -] [[GAMMA].sub.H] as follows :

[mathematical expression not reproducible] (20)

When applying the CIP method, spatial gradient at the next time step is also needed. Take spatial derivative of (15)

[partial derivative] ([[partial derivative].sub.x] [R.sup.[+ or -]])/[partial derivative]t + [[lambda].sup.[+ or -]] [partial derivative] ([[partial derivative].sub.x] [R.sup.[+ or -]])/[partial derivative]x = [G.sup.[+ or -]], (21)

[G.sup.[+ or -]] includes the terms related to spatial derivatives of [[lambda].sub.[+ or -]], [[GAMMA].sub.H], and u. Considering that the influence of [G.sup.[+ or -]] is not significant in this problem, we neglect [G.sup.[+ or -]] here. Then (21) can be solved as

[mathematical expression not reproducible], (22)

[mathematical expression not reproducible], (23)

then values and derivatives of each value at time step (n + 1) are obtained.

2.3. Two-Dimensional Tide Wave Equations. Assuming that pressure is hydrostatic and density is constant, the depth averaged two-dimensional tidal model without horizontal eddy viscosity is given as follows :

[mathematical expression not reproducible], (24)

where u(x, y, t) and v(x, y, t) are x and y components of flow velocity and f represents the Coriolis parameter.

In the same way as one dimension, replace variable C with water level H = [zeta] + [h.sub.0] + z; we get the equivalent of (24) which can be written in vector-matrix form as

[mathematical expression not reproducible], (25)

Since matrices A and B are not commutative; thus, there is no matrix to diagonalize them at the same time. Split (25) into two sequential phases; the source term is divided into two terms added to each directional phase:

[partial derivative]W/[partial derivative]t + A [partial derivative]W/[partial derivative]x + [S.sub.x] = 0 [W.sup.n] [right arrow] [W.sup.*], (27)

[partial derivative][W.sup.*]/[partial derivative]t + B [partial derivative][W.sup.*]/[partial derivative]y + [S.sub.y] = 0 [W.sup.*] [right arrow] [W.sup.n+1], (28)

where

[mathematical expression not reproducible] (29)

In order to get high accuracy in time, we calculate alternately in the x- and y-directions as follows:

[W.sup.n+2] = [L.sub.x][L.sub.y] [L.sub.y][L.sub.x][W.sup.n], (30)

where [L.sub.x] and [L.sub.y] represent operation of (27) and (28), respectively.

Here we give the solution process in x-direction. Similar to the one-dimensional case, left multiply (27) by the eigenvector matrix L:

[mathematical expression not reproducible] (31)

where

[mathematical expression not reproducible] (32)

Equation (31) leads to the following three equations:

[mathematical expression not reproducible] (33)

[mathematical expression not reproducible] (34)

where [R.sup.[+ or -].sub.x] = [[GAMMA].sub.H] [+ or -] u/2 are the Riemann invariants. These three equations all have the form of convective equations that can be calculated by CIP method. The Coriolis force, bottom topography, and bottom friction term are added to each Riemann invariant along characteristic line. Thus, the Riemann invariants and V can be discretized as follows:

[mathematical expression not reproducible] (35)

[mathematical expression not reproducible] (36)

So we have

[mathematical expression not reproducible] (37)

[mathematical expression not reproducible] (38)

[mathematical expression not reproducible] (39)

where [k'.sub.1] = 1 + (1/4)k[DELTA]t([square root of ([([u.sup.n]).sup.2] + [([v.sup.n]).sup.2])]/([H.sup.n] - z)), [k'.sub.2] = 1 + (1/4)k[DELTA]t([square root of ([([u.sup.n]).sup.2] + [([v.sup.n]).sup.2])]/([H.sup.n] - z)), f' = (1/4) f [DELTA]t,

[mathematical expression not reproducible] (40)

Similar to the derivation in one dimension, neglecting the influence of [G.sup.[+ or -]], the derivatives of each value at time step (n + 1) are obtained

[mathematical expression not reproducible] (41)

[mathematical expression not reproducible] (42)

[[partial derivative].sub.x][v.sup.*] = [partial derivative][v.sup.0]/[partial derivative]x. (43)

The operation in the x-direction is completed.

The operation in y-direction is essentially the same as the x-direction and the y derivative is computed by linear interpolation

[mathematical expression not reproducible] (44)

where sgn([[lambda].sub.x]) = 1 ([[lambda].sub.x] [greater than or equal to] 0), -1 ([[lambda].sub.x] < 0).

3. Numerical Results

3.1. A Small Perturbation of One-Dimensional Steady State Water. To demonstrate that the CIP method can be applied to the computation of ocean tidal waves with bottom topography and bottom friction, we first consider the quasi-stationary test case given in . The difference is that we add the effect of the bottom friction here.

The bottom topography is given by

[mathematical expression not reproducible] (45)

The initial conditions are

[mathematical expression not reproducible] (46)

where [zeta] = 0.001 is constant amplitude of the perturbation. According to the propagation mechanism of wave, the small disturbance will be split into two smaller waves propagating to both sides. However, because of the Gibbs phenomenon, the small pulse of perturbation is a challenge for many numerical schemes . Numerical experiments are carried out with both the CIP and SI schemes; the surface level with Nx = 2000 cells at the final time T = 0.2 is shown in Figure 3. We keep [DELTA]t = T/1000 as the SI scheme is strictly restricted by the CFL condition, and [DELTA]t = T/1000, T/400, T/240 when using the CIP.

Since flow velocity u and v are far smaller than wave velocity [[GAMMA].sub.H], the CFL number is defined as

CFL [approximately equal to] [[GAMMA].sub.H max] x [DELTA]t/[DELTA]x. (47)

The experimental results are similar to those in [19, 20]. It can be seen that nonphysical oscillations occur when simulating marching wave using the SI scheme, waves are broken in the process of travel, and the amplitude of nonphysical oscillations increases gradually in time, which also means that the scheme does not have the ability to simulate large gradient surface change. However, the shape of the simulated wave is smooth and has no oscillations when using CIP until CFL = 2.61. Figure 4 gives the results computed by CIP with CFL=2.61 for k=0, 1, and 2. It can be seen that the separated wave height is less than half of the initial perturbation because of the bottom friction. Numerical experiments show that the CIP method can successfully simulate large gradient water surface with both bottom topography and bottom friction using large time step sizes.

3.2. Smooth Surface Wave Propagation Problem. Since boundary condition plays a significant role in the shallow water flow of an actual marine area, we consider a smooth surface wave propagation problem here with reflective boundary conditions. The conditions for the x-direction are

[partial derivative]H/[partial derivative][OMEGA] = 0,

[partial derivative]v/[partial derivative][OMEGA] = 0

u([OMEGA]) = 0. (48)

Since the spatial gradient is also required when using the CIP method, the spatial gradient condition at the boundary is obtained according to (48):

[[partial derivative].sub.x]H ([OMEGA]) = 0,

[[partial derivative].sub.x]v ([OMEGA]) = 0,

[[partial derivative].sub.x]u/[partial derivative][OMEGA] = 0. (49)

When the traced upstream departure point on the shore, the corresponding actual upstream departure point in the waters is approximated using reflective boundary conditions (see Figure 5).

The initial conditions are given as follows:

[mathematical expression not reproducible] (50)

where [[sigma].sub.x] = [[sigma].sub.y] = 25. Let the bottom friction coefficient k = 2 x [10.sup.-3] and the constant Coriolis parameter f be 1 x [10.sup.-3].

The computation domain is n = [0, 200] x [0, 200], and [DELTA]x and [DELTA]y are set to 1.0. The contour map of the water surface at t=20 and t=40 is shown in Figures 6 and 7, respectively. It can be seen that at t=20 the contour line is not a complete circle because of the impact of the bottom topography (similar to the one-dimensional case). At t = 40, the reflected waves are superimposed on each other and eventually form what is shown in Figure 7. The results show that when time step increases to three times to the SI method, CIP-MOC scheme still simulate the change process of water surface very well, even considering the bottom topography (see Figure 8), boundary, and Coriolis effect.

4. Conclusion

In this paper, we apply the CIP-MOC to the solution of tide wave equations using large time step size. Bottom topography, bottom friction, and Coriolis term are included to the equation of Riemann invariants as the source term. One-dimensional experiment shows that the CIP-MOC scheme can keep the shape very well and has no oscillations even when the time step is increased to four times of the SI scheme in the simulation of large gradient water surface. Bottom friction effect is also tested in the experiment. In two-dimensional case, the problem with reflected boundary conditions is considered. Numerical experiments show that the CIP-MOC scheme is succeeded in simulating the change process of water surface using time step that is three times of the SI method. All the results validated the good stability and low dispersion error of the CIP-MOC method.

https://doi.org/10.1155/2018/3469534

Data Availability

All results presented in the article were produced from model simulations. Therefore, there is no data to be made available. Researchers who wish to replicate the study will use the equations and parameters described in the article. With such equations and parameters, researchers can use modeling simulations to replicate the figures presented in the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

Partial support for this research was provided by the National Key Research and Development Plan through Grant 2016YFC1402304, the Key Research and Development Plan of Shandong Province through Grants 2016ZDJS09A02 and ZR2016AB16, and the National Natural Science Foundation of China through Grant 41606006.

References

 J. Pedlosky, Geophysical fluid dynamics, Springer Science & Business Media, 2013.

 J. O. Backhaus, "A semi-implicit scheme for the shallow water equations for application to shelf sea modelling," Continental Shelf Research, vol. 2, no. 4, pp. 243-254, 1982.

 V. Casulli, "Semi-implicit finite difference methods for the two-dimensional shallow water equations," Journal of Computational Physics, vol. 86, no. 1, pp. 56-74, 1990.

 X. Lu and J. Zhang, "Numerical study on spatially varying bottom friction coefficient of a 2D tidal model with adjoint method," Continental Shelf Research, vol. 26, no. 16, pp. 1905-1923, 2006.

 H. Pan, Z. Guo, and X. Lv, "Inversion of tidal open boundary conditions of the M2 constituent in the bohai and yellow seas," Journal of Atmospheric and Oceanic Technology, vol. 34, no. 8, pp. 1661-1672, 2017.

 Z. Guo, H. Pan, W. Fan, and X. Lv, "Application of surface spline interpolation in inversion of bottom friction coefficients," Journal of Atmospheric and Oceanic Technology, vol. 34, no. 9, pp. 2021-2028, 2017.

 R. H. Weisberg and L. Zheng, "A finite volume coastal ocean model simulation of the tide, buoyancy, and wind-driven circulation of Tampa Bay," Journal of Geophysical Research, vol. 111, article C01005, 2006.

 C. Chen, H. Huang, R. C. Beardsley, H. Liu, Q. Xu, and G. Cowles, "A finite volume numerical approach for coastal ocean circulation studies: comparisons with finite difference models," Journal of Geophysical Research: Oceans, vol. 112, no. 3, Article ID C03018, 2007.

 G. Erbes, "A Semi-Lagrangian Method of Characteristics for the Shallow-Water Equations," Monthly Weather Review, vol. 121, no. 12, pp. 3443-3452, 1993.

 Y. Ogata and T. Yabe, "Multi-dimensional semi-lagrangian characteristic approach to the shallow water equations by the CIP method," International Journal of Computational Engineering Science, vol. 5, no. 3, pp. 699-730, 2004.

 K. Toda, Y. Ogata, and T. Yabe, "Multi-dimensional conservative semi-Lagrangian method of characteristics CIP for the shallow water equations," Journal of Computational Physics, vol. 228, no. 13, pp. 4917-4944, 2009.

 T. Yabe and T. Aoki, "A universal solver for hyperbolic equations by cubic-polynomial interpolation. I. One-dimensional solver," Computer Physics Communications, vol. 66, no. 2-3, pp. 219-232, 1991.

 T. Yabe and E. Takei, "A new higher-order Godunov method for general hyperbolic equations," Journal of the Physical Society of Japan, vol. 57, no. 8, pp. 2598-2601, 1988.

 Y. Anle and L. Fengqi, Physical Oceanography, China Ocean University Press, 1992.

 J. G. Zhou, D. M. Causon, C. G. Mingham, and D. M. Ingram, "The surface gradient method for the treatment of source terms in the shallow-water equations," Journal of Computational Physics, vol. 168, no. 1, pp. 1-25, 2001.

 R. Akoh, S. Ii, and F. Xiao, "A CIP/multi-moment finite volume method for shallow water equations with source terms," International Journal for Numerical Methods in Fluids, vol. 56, no. 12, pp. 2245-2270, 2008.

 R. W. Lardner, "Optimal assimilation of current and surface elevation data in a two-dimensional numerical tidal model," Applied Mathematical Modelling, vol. 19, no. 1, pp. 30-38, 1995.

 R. J. LeVeque, "Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wave-propagation algorithm," Journal of Computational Physics, vol. 146, no. 1, pp. 346-365, 1998.

 Q. Zhu, Z. Gao, W. S. Don, and X. Lv, "Well-balanced hybrid compact-WENO scheme for shallow water equations," Applied Numerical Mathematics, vol. 112, pp. 65-78, 2017.

 Y. Xing and C.-W. Shu, "High order finite difference WENO schemes with the exact conservation property for the shallow water equations," Journal of Computational Physics, vol. 208, no. 1, pp. 206-227, 2005.

Yafei Nie, (1) Kai Fu (iD), (2) and Xianqing Lv (iD) (1)

(1) Physical Oceanography Laboratory, Ocean University of China, Qingdao 266100, China

(2) School of Mathematical Sciences, Ocean University of China, Qingdao 266100, China

Correspondence should be addressed to Kai Fu; kfu@ouc.edu.cn

Received 2 February 2018; Accepted 26 May 2018; Published 2 July 2018

Caption: Figure 1: Shallow water variables.

Caption: Figure 2: Space-time diagram that shows the characteristics [C.sup.[+ or -]]. The Riemann invariants consist of [[GAMMA].sup.[+ or -].sub.H] and [u.sup.[+ or -]], [[lambda].sup.[+ or -]] are characteristic speeds.

Caption: Figure 3: A small perturbation of one-dimensional steady state water with a pulse [zeta] = 0.001. (a) shows SI scheme with CFL=0.63; (b), (c), and (d) are the results of the CIP method with CFL=0.63, 1.57, and 2.61.

Caption: Figure 4: Effect of bottom friction. Results computed by CIP with CFL=2.61 for k=0, 1, and 2.

Caption: Figure 5: Space-time diagram that shows the reflecting boundary conditions. The filled circle represents the approximated upstream departure point, while the hollow circle represents the virtual upstream point on the shore.

Caption: Figure 6: Smooth surface wave propagation problem. The contour of the water level H at t=20. (a) is the result of SI with CFL=0.63; (b), (c), and (d) are the results of CIP-MOC with CFL=0.63, 1.42, and 1.90.

Caption: Figure 7: Smooth surface wave propagation problem. The contour of the water level H at t=40. (a) is the result of SI with CFL=0.63; (b), (c), and (d) are the results of CIP-MOC with CFL=0.63, 1.42, and 1.90.

Caption: Figure 8: The comparison of cross-sections between topographic (a) and nontopographic conditions (b) with different CFL numbers.
Title Annotation: Printer friendly Cite/link Email Feedback Research Article; Constrained Interpolation Profile Nie, Yafei; Fu, Kai; Lv, Xianqing Advances in Mathematical Physics Report 1USA Jan 1, 2018 3831 Quantum Walk in Terms of Quantum Bernoulli Noise and Quantum Central Limit Theorem for Quantum Bernoulli Noise. On the Convergence Ball and Error Analysis of the Modified Secant Method. Interpolation Tsunamis Wave equation Wave equations