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

Authors: Yafei Nie; Kai Fu; Xianqing Lv
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

The CIP-MOC (Constrained Interpolation Profile/Method of Characteristics) is proposed to solve the tide wave equations with large time step size. The bottom topography and bottom friction, which are very important factor for the tidal wave model, are included to the equation of Riemann invariants as the source term. Numerical experiments demonstrate the good performance of the scheme. Compared to traditional semi-implicit (SI) finite difference scheme which is widely used in tidal wave simulation, CIP-MOC has better stability in simulating large gradient water surface change and has the ability to use much longer time step size under the premise of maintaining accuracy. Besides, numerical tests with reflective boundary conditions are carried out by CIP-MOC with large time step size and good results are obtained.


---

## Body

## 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) [1]. 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 [2] and Casulli [3] used semi-implicit scheme (hereafter SI) for the solution of shallow water equations; Lv and Zhang [4] used the semi-implicit scheme to solve tide wave equations and their computational format was used to study bottom friction coefficients [5] and tidal open boundary conditions [6]; 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 [9].There are growing interests in finding methods that enable using large time steps to solve shallow water equations. Erbes [9] 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 [10] 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. [11] 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 [9] 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. Section2 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 [12]. 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: (1)∂f∂t+u∂f∂x=0,the solution is [13](2)fn+1=Fx-uΔt,∂xfn+1=dFx-uΔtdx.Iff and ∂xf are known at each grid point, the profile between two adjacent points can be approximated by cubic polynomial(3)FiX=aiX3+biX2+∂xfiX+fi,where X=x-xi. Consider four constraints(4)Fxi=fin,Fxiup=fiupn,∂xFxi=∂xfin,∂xFxiup=∂xfiupn.then coefficients can be given as (5)ai=∂xfi+∂xfiupD2+2fi-fiupD3,(6)bi=3fiup-fiD2-2∂xfi+∂xfiupD,where sgn(ui)=1(ui≥0),-1(ui<0), iup=i-sgn(ui), and D=-Δx·sgn(ui).Thus, the profile off and ∂xf at next time step can be obtained as(7)fin+1=Fixi-uiΔt=aiξ3+biξ2+∂xfin·ξ+fin,(8)∂xfin+1=dFixi-uiΔtdx=3aiξ2+2biξ+∂xfin,where ξ=-uiΔ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 [14](9)∂ζ∂t+∂h0+ζu∂x=0,∂u∂t+u∂u∂x+g∂ζ∂x+ku=0,where t is time, x is Cartesian coordinate, h0x is undisturbed water depth, ζ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 [15] is used here. Replace the variable ζ with the water level H=ζ+h0+z, where zx represents the bottom topographic height (see Figure 1), we get the equivalent of (9) which can be written as (10)∂H∂t+∂Hu∂x=∂uz∂x,∂u∂t+u∂u∂x+g∂H∂x+ku=0.Figure 1
Shallow water variables.In a vector-matrix form, (10) can be rewritten as (11)∂∂tHu+uHgu∂∂xHu≡∂W∂t+AW∂W∂x=∂uz∂x-ku.The eigenvalue matrix Λ and eigenvector matrix L of coefficient matrix A are (12)Λ=λ+00λ-,L=ΓHg-ΓHg11,where ΓH=gH and λ±=u±ΓH. Thus, matrix A can be diagonalized as L-1AL, where (13)
L
-
1
=
1
2
g
Γ
H
1
-
g
Γ
H
1is the inverse of L. Then (11) can be written as (14)L-1∂W∂t+L-1ALL-1∂W∂x=L-1S,where S represents the right hand side of (11).Two characteristics about Riemann invariantsR±=ΓH±u/2 can be obtained from (14) as (15)∂R±∂t+λ±∂R±∂x=g2ΓH∂uz∂x∓12ku.By using CIP method, we have(16)ΓHn+1±un+12=ΓH±±u±2+g2∫tntn+11ΓH∂uz∂xdτ±∓12k·Δtu±+un+12,where ΓH± and u± are calculated by(17)ΓH±=ΓHCIPxi-λ±Δt,u±=uCIPxi-λ±Δt,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. ΓH and u at time step tn+1 can be obtained by solving linear equations (16) (see Figure 2):(18)ΓHn+1=12ΓH++ΓH-+12u+-u-+g2∫tntn+11ΓH∂uz∂xdτ++∫tntn+11ΓH∂uz∂xdτ--12kΔtu+-u-2,(19)un+1=1/2u++u-+2ΓH+-ΓH-+g∫tntn+11/ΓH∂uz/∂xdτ+-∫tntn+11/ΓH∂uz/∂xdτ--1/2kΔtu++u-/21+1/2kΔt.Figure 2
Space-time diagram that shows the characteristicsC±. The Riemann invariants consist of ΓH± and u±, λ± are characteristic speeds.Time integration of the bottom topographic term in (18) and (19) can be approximated using dx±/dτ±=u±ΓH as follows [16]:(20)∫tntn+11ΓH∂uz∂xdτ±≈-uz±-uznΓH±u±ΓH±.When applying the CIP method, spatial gradient at the next time step is also needed. Take spatial derivative of (15) (21)∂∂xR±∂t+λ±∂∂xR±∂x=G±,G± includes the terms related to spatial derivatives of λ±, ΓH, and u. Considering that the influence of G± is not significant in this problem, we neglect G± here. Then (21) can be solved as(22)∂xΓHn+1=12∂xΓH++∂xΓH-+12∂xu+-∂xu-,(23)∂xun+1=12∂xu++∂xu-+2∂xΓH+-∂xΓH-,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 [17]:(24)∂ζ∂t+∂h0+ζu∂x+∂h0+ζv∂y=0,∂u∂t+u∂u∂x+v∂u∂y-fv+kuu2+v2h0+ζ+g∂ζ∂x=0,∂v∂t+u∂v∂x+v∂v∂y+fu+kvu2+v2h0+ζ+g∂ζ∂y=0,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ζ with water level H=ζ+h0+z; we get the equivalent of (24) which can be written in vector-matrix form as (25)∂W∂t+AW∂W∂x+BW∂W∂y+S=0,where (26)
W
=
H
u
v
,
A
W
=
u
H
0
g
u
0
0
0
u
,
B
W
=
v
0
H
0
v
0
g
0
v
,
S
=
-
∂
u
z
∂
x
+
∂
v
z
∂
y
-
f
v
+
k
u
u
2
+
v
2
H
-
z
f
u
+
k
v
u
2
+
v
2
H
-
z
.Since matricesA 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:(27)∂W∂t+A∂W∂x+Sx=0Wn⟶W∗,(28)∂W∗∂t+B∂W∗∂y+Sy=0W∗⟶Wn+1,where(29)
S
x
=
-
∂
u
z
∂
x
-
1
2
f
v
+
1
2
k
u
u
2
+
v
2
H
-
z
1
2
f
u
+
1
2
k
v
u
2
+
v
2
H
-
z
,
S
y
=
-
∂
v
z
∂
y
-
1
2
f
v
+
1
2
k
u
u
2
+
v
2
H
-
z
1
2
f
u
+
1
2
k
v
u
2
+
v
2
H
-
z
.In order to get high accuracy in time, we calculate alternately in the x- and y-directions as follows:(30)Wn+2=LxLyLyLxWn,where Lx and Ly 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:(31)L-1∂W∂t+L-1ALL-1∂W∂x+L-1Sx≡L-1∂W∂t+λx-000λx+000uL-1∂W∂x+L-1Sx=0,where (32)
λ
x
±
=
u
±
Γ
H
,
L
=
-
Γ
H
g
Γ
H
g
0
1
1
0
0
0
1
,
L
-
1
=
1
2
-
g
Γ
H
1
0
g
Γ
H
1
0
0
0
2
.Equation (31) leads to the following three equations:(33)∂Rx±∂t+λx±∂Rx±∂x∓14fv±14kuu2+v2H-z-g2ΓH∂uz∂x=0,(34)∂v∂t+u∂v∂x+12fu+12kvu2+v2H-z=0,where Rx±=ΓH±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:(35)ΓH∗±12u∗=ΓH±±12u±±14fΔtv±+v∗2∓14kΔtu±+u∗2un2+vn2Hn-z+g2∫tntn+11ΓH∂uz∂xdτ±,(36)v∗=v0-12fΔtu0+u∗2-12kΔtv0+v∗2un2+vn2Hn-z.So we have(37)ΓH∗=12ΓH++ΓH-+k2′2u+-u-+f′2v+-v-+g2sx++sx-,(38)u∗=1/2k1′u++u-+2ΓH+-ΓH-+f′v++v-+2f′/k1′k2′v0-f′u0+gsx+-sx-k1′+f′2/k1,(39)v∗=k2′v0-f′u0+u∗k1′,where k1′=1+1/4kΔtun2+vn2/Hn-z, k2′=1-1/4kΔtun2+vn2/Hn-z, f′=1/4fΔt,(40)sx+=g2∫tntn+11ΓH∂uz∂xdτ+,sx-=g2∫tntn+11ΓH∂uz∂xdτ-.Similar to the derivation in one dimension, neglecting the influence of G±, the derivatives of each value at time step (n+1) are obtained(41)∂xΓH∗=12∂xΓH++∂xΓH-+12∂xu+-∂xu-,(42)∂xu∗=12∂xu++∂xu-+2∂xΓH+-∂xΓH-,(43)∂xv∗=∂v0∂x.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(44)∂yfi,j∗=∂yfi,j-sgnλx·∂yfi,j-∂yfi-sgnλx,jΔxλxΔt,where sgn(λx)=1(λx≥0),-1(λx<0).


## 2.1. CIP Method

CIP method is a compact, robust, less diffusive, and high-order scheme in computational fluid dynamics [12]. 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: (1)∂f∂t+u∂f∂x=0,the solution is [13](2)fn+1=Fx-uΔt,∂xfn+1=dFx-uΔtdx.Iff and ∂xf are known at each grid point, the profile between two adjacent points can be approximated by cubic polynomial(3)FiX=aiX3+biX2+∂xfiX+fi,where X=x-xi. Consider four constraints(4)Fxi=fin,Fxiup=fiupn,∂xFxi=∂xfin,∂xFxiup=∂xfiupn.then coefficients can be given as (5)ai=∂xfi+∂xfiupD2+2fi-fiupD3,(6)bi=3fiup-fiD2-2∂xfi+∂xfiupD,where sgn(ui)=1(ui≥0),-1(ui<0), iup=i-sgn(ui), and D=-Δx·sgn(ui).Thus, the profile off and ∂xf at next time step can be obtained as(7)fin+1=Fixi-uiΔt=aiξ3+biξ2+∂xfin·ξ+fin,(8)∂xfin+1=dFixi-uiΔtdx=3aiξ2+2biξ+∂xfin,where ξ=-uiΔ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 [14](9)∂ζ∂t+∂h0+ζu∂x=0,∂u∂t+u∂u∂x+g∂ζ∂x+ku=0,where t is time, x is Cartesian coordinate, h0x is undisturbed water depth, ζ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 [15] is used here. Replace the variable ζ with the water level H=ζ+h0+z, where zx represents the bottom topographic height (see Figure 1), we get the equivalent of (9) which can be written as (10)∂H∂t+∂Hu∂x=∂uz∂x,∂u∂t+u∂u∂x+g∂H∂x+ku=0.Figure 1
Shallow water variables.In a vector-matrix form, (10) can be rewritten as (11)∂∂tHu+uHgu∂∂xHu≡∂W∂t+AW∂W∂x=∂uz∂x-ku.The eigenvalue matrix Λ and eigenvector matrix L of coefficient matrix A are (12)Λ=λ+00λ-,L=ΓHg-ΓHg11,where ΓH=gH and λ±=u±ΓH. Thus, matrix A can be diagonalized as L-1AL, where (13)
L
-
1
=
1
2
g
Γ
H
1
-
g
Γ
H
1is the inverse of L. Then (11) can be written as (14)L-1∂W∂t+L-1ALL-1∂W∂x=L-1S,where S represents the right hand side of (11).Two characteristics about Riemann invariantsR±=ΓH±u/2 can be obtained from (14) as (15)∂R±∂t+λ±∂R±∂x=g2ΓH∂uz∂x∓12ku.By using CIP method, we have(16)ΓHn+1±un+12=ΓH±±u±2+g2∫tntn+11ΓH∂uz∂xdτ±∓12k·Δtu±+un+12,where ΓH± and u± are calculated by(17)ΓH±=ΓHCIPxi-λ±Δt,u±=uCIPxi-λ±Δt,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. ΓH and u at time step tn+1 can be obtained by solving linear equations (16) (see Figure 2):(18)ΓHn+1=12ΓH++ΓH-+12u+-u-+g2∫tntn+11ΓH∂uz∂xdτ++∫tntn+11ΓH∂uz∂xdτ--12kΔtu+-u-2,(19)un+1=1/2u++u-+2ΓH+-ΓH-+g∫tntn+11/ΓH∂uz/∂xdτ+-∫tntn+11/ΓH∂uz/∂xdτ--1/2kΔtu++u-/21+1/2kΔt.Figure 2
Space-time diagram that shows the characteristicsC±. The Riemann invariants consist of ΓH± and u±, λ± are characteristic speeds.Time integration of the bottom topographic term in (18) and (19) can be approximated using dx±/dτ±=u±ΓH as follows [16]:(20)∫tntn+11ΓH∂uz∂xdτ±≈-uz±-uznΓH±u±ΓH±.When applying the CIP method, spatial gradient at the next time step is also needed. Take spatial derivative of (15) (21)∂∂xR±∂t+λ±∂∂xR±∂x=G±,G± includes the terms related to spatial derivatives of λ±, ΓH, and u. Considering that the influence of G± is not significant in this problem, we neglect G± here. Then (21) can be solved as(22)∂xΓHn+1=12∂xΓH++∂xΓH-+12∂xu+-∂xu-,(23)∂xun+1=12∂xu++∂xu-+2∂xΓH+-∂xΓH-,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 [17]:(24)∂ζ∂t+∂h0+ζu∂x+∂h0+ζv∂y=0,∂u∂t+u∂u∂x+v∂u∂y-fv+kuu2+v2h0+ζ+g∂ζ∂x=0,∂v∂t+u∂v∂x+v∂v∂y+fu+kvu2+v2h0+ζ+g∂ζ∂y=0,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ζ with water level H=ζ+h0+z; we get the equivalent of (24) which can be written in vector-matrix form as (25)∂W∂t+AW∂W∂x+BW∂W∂y+S=0,where (26)
W
=
H
u
v
,
A
W
=
u
H
0
g
u
0
0
0
u
,
B
W
=
v
0
H
0
v
0
g
0
v
,
S
=
-
∂
u
z
∂
x
+
∂
v
z
∂
y
-
f
v
+
k
u
u
2
+
v
2
H
-
z
f
u
+
k
v
u
2
+
v
2
H
-
z
.Since matricesA 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:(27)∂W∂t+A∂W∂x+Sx=0Wn⟶W∗,(28)∂W∗∂t+B∂W∗∂y+Sy=0W∗⟶Wn+1,where(29)
S
x
=
-
∂
u
z
∂
x
-
1
2
f
v
+
1
2
k
u
u
2
+
v
2
H
-
z
1
2
f
u
+
1
2
k
v
u
2
+
v
2
H
-
z
,
S
y
=
-
∂
v
z
∂
y
-
1
2
f
v
+
1
2
k
u
u
2
+
v
2
H
-
z
1
2
f
u
+
1
2
k
v
u
2
+
v
2
H
-
z
.In order to get high accuracy in time, we calculate alternately in the x- and y-directions as follows:(30)Wn+2=LxLyLyLxWn,where Lx and Ly 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:(31)L-1∂W∂t+L-1ALL-1∂W∂x+L-1Sx≡L-1∂W∂t+λx-000λx+000uL-1∂W∂x+L-1Sx=0,where (32)
λ
x
±
=
u
±
Γ
H
,
L
=
-
Γ
H
g
Γ
H
g
0
1
1
0
0
0
1
,
L
-
1
=
1
2
-
g
Γ
H
1
0
g
Γ
H
1
0
0
0
2
.Equation (31) leads to the following three equations:(33)∂Rx±∂t+λx±∂Rx±∂x∓14fv±14kuu2+v2H-z-g2ΓH∂uz∂x=0,(34)∂v∂t+u∂v∂x+12fu+12kvu2+v2H-z=0,where Rx±=ΓH±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:(35)ΓH∗±12u∗=ΓH±±12u±±14fΔtv±+v∗2∓14kΔtu±+u∗2un2+vn2Hn-z+g2∫tntn+11ΓH∂uz∂xdτ±,(36)v∗=v0-12fΔtu0+u∗2-12kΔtv0+v∗2un2+vn2Hn-z.So we have(37)ΓH∗=12ΓH++ΓH-+k2′2u+-u-+f′2v+-v-+g2sx++sx-,(38)u∗=1/2k1′u++u-+2ΓH+-ΓH-+f′v++v-+2f′/k1′k2′v0-f′u0+gsx+-sx-k1′+f′2/k1,(39)v∗=k2′v0-f′u0+u∗k1′,where k1′=1+1/4kΔtun2+vn2/Hn-z, k2′=1-1/4kΔtun2+vn2/Hn-z, f′=1/4fΔt,(40)sx+=g2∫tntn+11ΓH∂uz∂xdτ+,sx-=g2∫tntn+11ΓH∂uz∂xdτ-.Similar to the derivation in one dimension, neglecting the influence of G±, the derivatives of each value at time step (n+1) are obtained(41)∂xΓH∗=12∂xΓH++∂xΓH-+12∂xu+-∂xu-,(42)∂xu∗=12∂xu++∂xu-+2∂xΓH+-∂xΓH-,(43)∂xv∗=∂v0∂x.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(44)∂yfi,j∗=∂yfi,j-sgnλx·∂yfi,j-∂yfi-sgnλx,jΔxλxΔt,where sgn(λx)=1(λx≥0),-1(λ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 [18]. The difference is that we add the effect of the bottom friction here.The bottom topography is given by(45)bx=0.25cos⁡10πx-1.5+1,if1.4≤x≤1.6,0,otherwise.The initial conditions are(46)Hx,0=1+ζ,if1.1≤x≤1.2,1,otherwise.ux,0=0,k=2.0,where ζ=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 [19]. 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 Δt=T/1000 as the SI scheme is strictly restricted by the CFL condition, and Δt=T/1000,T/400,T/240 when using the CIP.Figure 3
A small perturbation of one-dimensional steady state water with a pulseζ=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.
(a) (b) (c) (d)Since flow velocityu and v are far smaller than wave velocity ΓH, the CFL number is defined as (47)CFL≈ΓHmax·ΔtΔx.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.Figure 4
Effect of bottom friction. Results computed by CIP with CFL=2.61 for k=0, 1, and 2.


### 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(48)∂H∂Ω=0,∂v∂Ω=0,uΩ=0.Since the spatial gradient is also required when using the CIP method, the spatial gradient condition at the boundary is obtained according to (48): (49)∂xHΩ=0,∂xvΩ=0,∂xu∂Ω=0.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 Figure5).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.The initial conditions are given as follows:(50)ζx,y,0=0.001·exp⁡-x-Xmax/2σx2-y-Ymax/2σy2,zx,y=0.5·exp⁡-x-Xmax/3σx2-y-Ymax/3σy2,Hx,y,0=1+ζx,y,0-zx,y,ux,y,0=0,vx,y,0=0,where σx=σy=25. Let the bottom friction coefficient k=2×10-3 and the constant Coriolis parameter f be 1×10-3.The computation domain isΩ=0,200×0,200, and Δx and Δ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.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.
(a) (b) (c) (d)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.
(a) (b) (c) (d)Figure 8
The comparison of cross-sections between topographic (a) and nontopographic conditions (b) with different CFL numbers.
(a) (b)


## 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 [18]. The difference is that we add the effect of the bottom friction here.The bottom topography is given by(45)bx=0.25cos⁡10πx-1.5+1,if1.4≤x≤1.6,0,otherwise.The initial conditions are(46)Hx,0=1+ζ,if1.1≤x≤1.2,1,otherwise.ux,0=0,k=2.0,where ζ=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 [19]. 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 Δt=T/1000 as the SI scheme is strictly restricted by the CFL condition, and Δt=T/1000,T/400,T/240 when using the CIP.Figure 3
A small perturbation of one-dimensional steady state water with a pulseζ=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.
(a) (b) (c) (d)Since flow velocityu and v are far smaller than wave velocity ΓH, the CFL number is defined as (47)CFL≈ΓHmax·ΔtΔx.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.Figure 4
Effect of bottom friction. Results computed by CIP with CFL=2.61 for k=0, 1, and 2.


## 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(48)∂H∂Ω=0,∂v∂Ω=0,uΩ=0.Since the spatial gradient is also required when using the CIP method, the spatial gradient condition at the boundary is obtained according to (48): (49)∂xHΩ=0,∂xvΩ=0,∂xu∂Ω=0.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 Figure5).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.The initial conditions are given as follows:(50)ζx,y,0=0.001·exp⁡-x-Xmax/2σx2-y-Ymax/2σy2,zx,y=0.5·exp⁡-x-Xmax/3σx2-y-Ymax/3σy2,Hx,y,0=1+ζx,y,0-zx,y,ux,y,0=0,vx,y,0=0,where σx=σy=25. Let the bottom friction coefficient k=2×10-3 and the constant Coriolis parameter f be 1×10-3.The computation domain isΩ=0,200×0,200, and Δx and Δ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.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.
(a) (b) (c) (d)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.
(a) (b) (c) (d)Figure 8
The comparison of cross-sections between topographic (a) and nontopographic conditions (b) with different CFL numbers.
(a) (b)


## 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.



---

*Source: 3469534-2018-07-02.xml*