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.
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. 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 [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: the solution is [13]
If and are known at each grid point, the profile between two adjacent points can be approximated by cubic polynomialwhere . Consider four constraintsthen coefficients can be given as where , , and .
Thus, the profile of and at next time step can be obtained aswhere .
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]where is time, is Cartesian coordinate, is undisturbed water depth, is sea surface elevation above the undisturbed sea level, is the flow velocity, is the bottom friction coefficient, and 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 , where represents the bottom topographic height (see Figure 1), we get the equivalent of (9) which can be written as
In a vector-matrix form, (10) can be rewritten as The eigenvalue matrix and eigenvector matrix of coefficient matrix are where and . Thus, matrix can be diagonalized as , where is the inverse of . Then (11) can be written as where represents the right hand side of (11).
Two characteristics about Riemann invariants can be obtained from (14) as By using CIP method, we havewhere and are calculated bythe 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. and at time step can be obtained by solving linear equations (16) (see Figure 2):
Time integration of the bottom topographic term in (18) and (19) can be approximated using as follows [16]:
When applying the CIP method, spatial gradient at the next time step is also needed. Take spatial derivative of (15) includes the terms related to spatial derivatives of , , and . Considering that the influence of is not significant in this problem, we neglect here. Then (21) can be solved asthen values and derivatives of each value at time step 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]:where and are and components of flow velocity and represents the Coriolis parameter.
In the same way as one dimension, replace variable with water level ; we get the equivalent of (24) which can be written in vector-matrix form as where
Since matrices and 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:where
In order to get high accuracy in time, we calculate alternately in the x- and y-directions as follows:where and 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 :where Equation (31) leads to the following three equations:where 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 can be discretized as follows:So we havewhere , , ,Similar to the derivation in one dimension, neglecting the influence of , the derivatives of each value at time step are obtainedThe 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 interpolationwhere .
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 The initial conditions arewhere 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 cells at the final time is shown in Figure 3. We keep as the SI scheme is strictly restricted by the CFL condition, and when using the CIP.
(a)
(b)
(c)
(d)
Since flow velocity and are far smaller than wave velocity , the CFL number is defined as
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 . 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
Since the spatial gradient is also required when using the CIP method, the spatial gradient condition at the boundary is obtained according to (48):
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:where . Let the bottom friction coefficient and the constant Coriolis parameter be .
The computation domain is , and and 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 , 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(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.
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.