A CDG-FE method for the two-dimensional Green-Naghdi model with the enhanced dispersive property

TL;DR

This paper develops a numerical method combining central discontinuous Galerkin and finite element techniques to solve the two-dimensional Green-Naghdi model with enhanced dispersive properties, ensuring accuracy, efficiency, and physical fidelity.

Contribution

It introduces a novel CDG-FE scheme for the 2D Green-Naghdi model that improves dispersive effects and maintains well-balanced and positivity-preserving features.

Findings

The proposed method accurately captures dispersive wave behavior.

Numerical tests demonstrate the scheme's efficiency and stability.

The approach effectively preserves physical properties like positivity.

Abstract

In this work, we investigate numerical solutions of the two-dimensional shallow water wave using a fully nonlinear Green-Naghdi model with an improved dispersive effect. For the purpose of numerics, the Green-Naghdi model is rewritten into a formulation coupling a pseudo-conservative system and a set of pseudo-elliptic equations. Since the pseudo-conservative system is no longer hyperbolic and its Riemann problem can only be approximately solved, we consider the utilization of the central discontinuous Galerkin method which possesses an important feature of needlessness of Riemann solvers. Meanwhile, the stationary elliptic part will be solved using the finite element method. Both the well-balanced and the positivity-preserving features which are highly desirable in the simulation of the shallow water wave will be embedded into the proposed numerical scheme. The accuracy and efficiency of the numerical model and method will be illustrated through numerical tests.

Tables

Table 1. Table 1: L 2 superscript 𝐿 2 L^{2} errors and orders of accuracy of ( h , u ) ℎ 𝑢 (h,u) .

		$h$						$u$
$\Delta x$		$P^1$			$P^2$			$P^1$			$P^2$
		$L^2$ error	Order		$L^2$ error	Order		$L^2$ error	Order		$L^2$ error	Order
1		2.28E-01	—		7.80E-02	—		5.16E-01	—		1.05E-01	—
0.5		6.00E-02	1.93		9.94E-03	2.97		1.27E-02	2.03		1.47E-02	2.84
0.25		1.53E-02	1.97		1.27E-03	2.97		2.94E-02	2.11		1.84E-03	2.99
0.125		3.53E-03	2.12		1.64E-04	2.96		7.03E-03	2.07		2.29E-04	3.01

Table 2. Table 2: L ∞ superscript 𝐿 L^{\infty} errors on ( h + b , u , v ) ℎ 𝑏 𝑢 𝑣 (h+b,u,v) for the stationary solution at t = 10 𝑡 10 t=10 .

Case		precision		$h+b$		$u$		$v$
A		single		5.96E-08		8.86E-10		9.73E-10
A		double		2.56E-16		2.81E-16		2.23E-16
B		single		5.96E-08		4.21E-09		2.45E-09
B		double		4.93E-16		5.32E-16		2.33E-16

Equations

$$\left\{\begin{array}[]{lclcl}h_{t}+(hu)_{x}+(hv)_{y}=0,\\ (hu)_{t}+\left(hu^{2}+\frac{1}{2}gh^{2}+\frac{1}{3}h^{3}\Phi+\frac{1}{2}h^{2}\Psi\right)_{x}+\left(huv\right)_{y}=-\left(gh+\frac{1}{2}h^{2}\Phi+h\Psi\right)b_{x},\\ (hv)_{t}+\left(huv\right)_{x}+\left(hv^{2}+\frac{1}{2}gh^{2}+\frac{1}{3}h^{3}\Phi+\frac{1}{2}h^{2}\Psi\right)_{y}=-\left(gh+\frac{1}{2}h^{2}\Phi+h\Psi\right)b_{y},\end{array}\right.$$

$$\Phi$$

$$\Psi$$

$$u_t$$

$$v_t$$

$$\left\{\begin{array}[]{lclcl}h_{t}+(hu)_{x}+(hv)_{y}=0,\\ (hu)_{t}+\left(hu^{2}+\frac{1}{2}gh^{2}+\frac{1}{3}h^{3}\Phi+\frac{1}{2}h^{2}\Psi\right)_{x}+\left(huv\right)_{y}=-\left(gh+\frac{1}{2}h^{2}\Phi+h\Psi\right)b_{x},\\ (hv)_{t}+\left(huv\right)_{x}+\left(hv^{2}+\frac{1}{2}gh^{2}+\frac{1}{3}h^{3}\Phi+\frac{1}{2}h^{2}\Psi\right)_{y}=-\left(gh+\frac{1}{2}h^{2}\Phi+h\Psi\right)b_{y}.\end{array}\right.$$

$$\Phi$$

$$\Psi$$

$$h=h_0+\tilde{h},u=\tilde{u},v=\tilde{v},$$

$$\left\{{\begin{array}[]{*{20}{l}}\tilde{h}_{t}+h_{0}\tilde{u}_{x}+h_{0}\tilde{v}_{y}=0~{},\\ \tilde{u}_{t}+g\tilde{h}_{x}-\frac{1}{3}h_{0}^{2}\left(\alpha\tilde{u}_{txx}+\alpha\tilde{v}_{tyx}+(\alpha-1)\tilde{h}_{xxx}+(\alpha-1)\tilde{h}_{yyx}\right)=0~{},\\ \tilde{v}_{t}+g\tilde{h}_{y}-\frac{1}{3}h_{0}^{2}\left(\alpha\tilde{u}_{txy}+\alpha\tilde{v}_{tyy}+(\alpha-1)\tilde{h}_{xxy}+(\alpha-1)\tilde{h}_{yyy}\right)=0~{}.\end{array}}\right.$$

$$\left(\begin{array}[]{l}\tilde{h}\\ \tilde{u}\\ \tilde{v}\end{array}\right)=\left(\begin{array}[]{l}\hat{h}\\ \hat{u}\\ \hat{v}\end{array}\right){e^{i(k_{1}x+k_{2}y-\omega t)}},$$

$$\left({\begin{array}[]{*{20}{c}}-\omega&h_{0}k_{1}&h_{0}k_{2}\\ gk_{1}+\frac{\alpha-1}{3}gh_{0}^{2}k_{1}(k_{1}^{2}+k_{2}^{2})&-\omega(1+\frac{\alpha}{3}h_{0}^{2}k_{1}^{2})&-\omega\frac{\alpha}{3}h_{0}^{2}k_{1}k_{2}\\ gk_{2}+\frac{\alpha-1}{3}gh_{0}^{2}k_{2}(k_{1}^{2}+k_{2}^{2})&-\omega\frac{\alpha}{3}h_{0}^{2}k_{1}k_{2}&-\omega(1+\frac{\alpha}{3}h_{0}^{2}k_{2}^{2})\end{array}}\right)\left({\begin{array}[]{*{20}{c}}\hat{h}\\ \hat{u}\\ \hat{v}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}0\\ 0\\ 0\end{array}}\right),$$

$$\omega ={\omega }_{\pm }=\pm |k|\sqrt{g{h}_0(1+\frac{\alpha -1}{3}{h}_0^2|k{|}^2){(1+\frac{\alpha }{3}{h}_0^2|k{|}^2)}^{-1}}$$

$$\omega ={\omega }_{\pm }=\pm \sqrt{g|k|\tanh (h_0|k|)},$$

$$u=0,v=0,h+b=\text{mbox}constant.$$

$$hP$$

$$hQ$$

$${\mathbf{U}}_t+\mathbf{F}(\mathbf{U},u,v;b{)}_x+\mathbf{G}(\mathbf{U},u,v;b{)}_y=\mathbf{S}(\mathbf{U},u,v;b),$$

$$\mathbf{F}(\mathbf{U},u,v;b)$$

$$\mathbf{G}(\mathbf{U},u,v;b)$$

$$\mathbf{S}(\mathbf{U},u,v;b)$$

$${\mathcal{V}}^C$$

$${\mathcal{V}}^D$$

$${\mathcal{W}}^C$$

$${\mathcal{W}}^D$$

$${\int }_{C_{ij}}{\mathbf{U}}^{n+1,C}\cdot \mathbf{V}dxdy$$

$${\int }_{\Omega }\left(\frac{\alpha }{3}(h^{n+1,C}{)}^3(u^{n+1,C}{)}_x+\frac{\alpha }{3}(h^{n+1,C}{)}^3(v^{n+1,C}{)}_y-\frac{\alpha }{2}(h^{n+1,C}{)}^2{v}^{n+1,C}(b^C{)}_y\right){\hat{u}}_{x}dxdy$$

$${\int }_{\Omega }\left(\frac{\alpha }{3}(h^{n+1,C}{)}^3(u^{n+1,C}{)}_x+\frac{\alpha }{3}(h^{n+1,C}{)}^3(v^{n+1,C}{)}_y-\frac{\alpha }{2}(h^{n+1,C}{)}^2{u}^{n+1,C}(b^C{)}_x\right){\hat{v}}_{y}dxdy$$

$$f(h^{n+1,C},u^{n+1,C},v^{n+1,C};b^C)$$

$$g(h^{n+1,C},u^{n+1,C},v^{n+1,C};b^C)$$

$$hP$$

$$hQ$$

$$u$$

Full text

††journal: arXiv

1 Introduction

In an incompressible and inviscid fluid, the propagation of surface waves is governed by the Euler equation with nonlinear boundary conditions at the free surface and the bottom. In its full generality, this problem is very complicated to be solved, both mathematically and numerically. Usually, simplified models have been derived to describe the behavior of the solution in some physical specific regimes, such as the nonlinear shallow water equations, the Korteweg-de Vries equations, the Boussinesq type models and the Green-Naghdi models, and to name a few.

The nonlinear shallow water equations (also called Saint-Venant equations) are a set of hyperbolic partial differential equations (or parabolic if viscous shear is considered) that describe the flow below a pressure surface in a fluid. They can model the propagation of strongly nonlinear waves up to breaking and run-up in near-shore zones. However, they fail to properly describe wave propagation in deep water or wave shoaling because they do not incorporate frequency dispersion. The Boussinesq systems carry weak frequency dispersion but are typically restricted to small amplitude waves with relatively weak nonlinearity. Extensions of the shallow water equations that incorporate frequency dispersion can be traced back to Serre ([37]) who derived a one-dimensional (1D) system of equations for fully nonlinear weakly dispersive waves over the flat bottom in 1953. In 1976, Green and Naghdi ([19]) presented the two-dimensional (2D) counterpart of these equations for wave propagation over variable bottom topography. The Green-Naghdi model is a class of fully nonlinear weakly dispersive shallow water wave equations, including dispersive effects and supporting traveling solitary wave solutions. Therefore, it can simulate the long-time propagation of solitary waves with relatively large amplitude.

Due to the significance of the Green-Naghdi model, there has been increasing interest in the numerical simulation of the Green-Naghdi model in the past decade. A fourth-order compact finite volume scheme was proposed for solving the fully nonlinear and weakly dispersive Boussinesq type equations ([7, 8]). A hybrid numerical method using a Godunov type scheme was presented to solve the Green-Naghdi model over the flat bottom ([24]). A pseudo-spectral algorithm was developed for the solution of the rotating Green-Naghdi shallow water equations ([35]). A hybrid finite volume and finite difference splitting approach was presented for numerical simulation of the fully nonlinear and weakly dispersive Green-Naghdi model ([1, 2, 3, 23]). In particular, numerical investigations on a dispersive-effect-improved Green-Naghdi model have been reported in [3, 23]. A well-balanced central discontinuous Galerkin (CDG) method coupling with the finite element (FE) method ([28]) was employed to solve the 1D fully nonlinear weakly dispersive Green-Naghdi model over varying topography. Recent works of other kinds of discontinuous Galerkin (DG) methods for the Green-Naghdi model include [16, 34].

The DG method is a class of high order finite element methods, which was originally introduced in 1973 by Reed and Hill ([36]) for the neutron transport equation. Thereafter, it has achieved great progress in a series of pioneer papers ([9, 10, 11, 12]). The DG method has its own advantages in dealing with numerical solutions for many problems in sciences and engineering, consisting of being very flexible to achieve high order of accuracy and handle complicated geometry and boundary conditions, and to name a few. As a variant of the DG method, the CDG method is also one of popular high order numerical methods defined on overlapping meshes, which was originally introduced for hyperbolic conservation laws ([32]), and then for diffusion equations ([33]). By evolving two sets of numerical solutions defined on overlapping meshes which provide more information on numerical solutions, the CDG method does not rely on any exact or approximate Riemann solver at element interfaces as in the simulation of the DG method. The CDG method has been successfully applied to solve various partial differential equations, such as the Hamilton-Jacobi equations ([25]), the ideal MHD equations ([27, 26]), the Euler equations ([30]) and the shallow water equations ([28, 29]). Recently, a reconstructed CDG method has been developed in [15] for the improvement of computational efficiency, and two kinds of CDG methods defined on unstructured overlapping meshes have also been presented in [43, 31] for the treatment of complex computational domains.

In this paper, we first derive a formulation describing the propagation of 2D fully nonlinear weakly dispersive water waves, using the same approach as in [38] where the authors derived the 1D equation over the flat bottom. Based on this formulation, we derive a Green-Naghdi model with an enhanced dispersive effect, and then carry out a linear dispersion analysis for the corresponding linear equation to show the improvement on dispersive effects. This improvement on computational modeling is of great importance during the procedure of seeking more accurate numerical solutions of shallow water waves. In simulation of the derived 2D Green-Naghdi model, we usually encounter several difficulties: dealing with the mixed spatial and temporal derivative terms in the flux gradient and the source term, particularly in the dispersive-effect-improved model to be considered in the current work; maintaining the nonnegativity of the water depth and preserving the still-water stationary solution; computational efficiency and accuracy in 2D simulation. Concerning the first challenge, we reformulate the model into a hybrid system of the pseudo-conservative form with a nonhomogeneous source term and the pseudo-elliptic equation ([24, 28]). Since the resulting dynamic equations are no longer hyperbolic, we utilize the CDG method ([32]) as the base scheme which has an important feature of being free of Riemann solvers. Meanwhile, we employ the FE method to solve the elliptic part, and point out that there is no need for the CDG-FE hybrid method to make special treatments on the continuity and discontinuity due to the fact of two copies of solutions on overlapping meshes. Secondly, the reason for naming the pseudo-elliptic equation in our numerical model lies in the fact that the computed variable of water depth appears in the coefficients of the second-order derivative term and may be trivial in the dry area. On the other hand, it explains that the nonnegativity of computational water depth is extremely important in our numerical model and numerical scheme, and this issue has not been solved in the previous 1D work ([28]) yet. In this work, we design a CDG method for the Green-Naghdi model not only maintaining the well-balanced property for the stationary solution but also preserving the nonnegativity for the solution at the dry area. In addition to these issues on numerical modeling and schemes, computational efficiency ([42, 23]) is also significantly important for 2D water wave simulations. Since the CDG method is a local scheme similar to the DG method, the computation is implemented element by element at each time step and is thus friendly to the parallel execution. The computation load of solving the FE equation is dominant in our computation, and its fast solver will be considered in the future work.

The remainder of the paper is organized as follows. In section 2, we present a fully nonlinear strongly dispersive water wave model in 2D domain. In section 3, we propose a family of high order schemes, coupling positivity-preserving well-balanced CDG methods and continuous FE methods. In section 4, we perform a series of numerical experiments to demonstrate the well-balanced property, positivity-preserving property, high order accuracy as well as the capability of the Green-Naghdi model to describe the propagation of strongly nonlinear and dispersive waves. Some concluding remarks are given in section 5.

2 Mathematical models

2.1 Green-Naghdi model

In this paper, we consider the following 2D Green-Naghdi model

[TABLE]

where

[TABLE]

The derivation of the model is shown in Appendix.

2.2 Green-Naghdi model with enhanced dispersion effect

We can observe from the last two equations in (1) that ([3])

[TABLE]

and

[TABLE]

Replacing $u_{t}$ and $v_{t}$ with $\alpha u_{t}+(\alpha-1)(uu_{x}+u_{y}v+g(h+b)_{x})$ and $\alpha v_{t}+(\alpha-1)(vv_{y}+uv_{x}+g(h+b)_{y})$ in (2) and (3), respectively, we obtain a modified Green-Naghdi model as follows

[TABLE]

with

[TABLE]

It is apparent that the standard Green-Naghdi model corresponds to a particular case of the modified Green-Naghdi model with $\alpha=1$.

2.3 Linear dispersive analysis

We will carry out a linear dispersion analysis for the modified Green-Naghdi model given in the previous section. We first linearize the model (4) for the flat-bottom case ($b=\text{constant}$) with the trivial solution $h=h_{0}>0$ ($h_{0}$ is a constant), $u=u_{0}=0$ and $v=v_{0}=0$, i.e. looking for solutions of the form

[TABLE]

where $(\tilde{h},\tilde{u},\tilde{v})$ are small perturbations, and retain only linear terms. As a result, we have

[TABLE]

Then we look for perturbations as plane waves taking the form

[TABLE]

where $k_{1}$ and $k_{2}$ are the wave numbers, and $\omega$ is the frequency of the perturbations. Without loss of generality, $k_{1}$ and $k_{2}$ are taken to be real and positive. We end up with the linear system

[TABLE]

and it has nontrivial solutions provided that the determinant of the coefficient matrix is zero, i.e.

[TABLE]

with $|k|=\sqrt{k_{1}^{2}+k_{2}^{2}}$. This gives the linear dispersion relation for the modified Green-Naghdi equations (4) in the flat-bottom case with respect to the trivial solution. In Figure 1, we compare the linear dispersion relation (7) ($\alpha=1.0$ and $\alpha=1.159$) and that of the full water wave problem in finite depth,

[TABLE]

with $g=1$, $h_{0}=1$. It can be observed from this figure that compared with the case $\alpha=1.0$, the linear dispersion relation of the modified model with $\alpha=1.159$ shows a better agreement with that of the full water wave problem.

3 Numerical methods for the Green-Naghdi models

The first difficulty in designing numerical schemes for the Green-Naghdi model comes from the appearance of mixed spatial and temporal derivatives in the equations. To tackle this, we make a reformulation from the original equations through introducing auxiliary variables, and the procedure will be given in Section 3.1. This new model will be adopted for the 2D shallow water wave simulations. In addition, in the case of a variable bottom, these equations admit still-water stationary solutions for the system (4) which are given by

[TABLE]

Numerical errors in discretization will produce the spurious oscillations which lead to a wrong set of solutions. Another major difficulty in the simulation of shallow water waves has to do with the appearance of dry areas. In particular, this case happens frequently in so many important applications involving rapidly moving interfaces between wet and dry areas, for instance, the wave run-up on beaches or over man-made structures, dam break and tsunami. If no special care is taken, non-physical phenomena may arise and it will lead to a breakdown of computation. Following the work in [41, 29], we design in Section 3.2 a high order positivity-preserving well-balanced CDG-FE method for the 2D simulation of the modified Green-Naghdi model in Section 3.1.

3.1 Reformulation of modified Green-Naghdi model

In order to remove the mixed spatial and temporal derivatives in the equations, we introduce two new unknown variables $P$ and $Q$ which satisfy:

[TABLE]

and

[TABLE]

Then the system (4) can be reformulated into a balance law

[TABLE]

where $\mathbf{U}=(h,hP,hQ)^{\top}$ is the unknown vector,

[TABLE]

and

[TABLE]

are the flux terms, and

[TABLE]

is the source term.

Based on the new equations, the solution of the modified Green-Naghdi model (4) amounts to finding the unknowns ${(h,hP,hQ)^{\top}}$ based on (11) and the unknowns ${(u,v)^{\top}}$ from (9)-(10). We remark that the linear dispersive relationship of the reformulated Green-Naghdi model (9)-(11) is the same as the one of (4).

3.2 CDG-FE methods

In this section, we develop the numerical method for the solution of the equations (9)-(11). Let $\mathcal{T}^{C}=\{C_{ij},\forall i,j\}$ and $\mathcal{T}^{D}=\{D_{ij},\forall i,j\}$ define two overlapping meshes for the computational domain $\Omega=[x_{\min},x_{\max}]\times[y_{\min},y_{\max}]$, with $C_{ij}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$, $D_{ij}=[x_{i-1},x_{i}]\times[y_{j-1},y_{j}]$, $x_{i+\frac{1}{2}}=\frac{1}{2}(x_{i}+x_{i+1})$ and $y_{j+\frac{1}{2}}=\frac{1}{2}(y_{j}+y_{j+1})$, where $\{x_{i}\}_{i}$ and $\{y_{j}\}_{j}$ are partitions of $[x_{\min},x_{\max}]$ and $[y_{\min},y_{\max}]$ respectively. Associated with each mesh, we define the following discrete spaces

[TABLE]

To approximate $u$ and $v$, we define two continuous finite element spaces

[TABLE]

We only present the schemes with the forward Euler method for time discretization. High-order time discretizations will be discussed in Section 3.3. The proposed methods evolve two copies of numerical solution, which are assumed to be available at $t=t_{n}$, denoted by ${\bf{U}}^{n,\star}=(h^{n,\star},(hP)^{n,\star},(hQ)^{n,\star})^{\top}\in\mathcal{V}^{\star}$, and we want to find the solutions at $t=t_{n+1}=t_{n}+\Delta t_{n}$. Only the procedure to update ${\bf{U}}^{n+1,C}$ will be described. We project the bottom topography function $b$ into $P^{k}(C_{ij})$ on $C_{ij}$ (resp. into $P^{k}(D_{ij})$ on $D_{ij}$) in the $L^{2}$ sense, and obtain an approximation $b^{C}$ (resp. $b^{D}$) throughout the computational domain.

3.2.1 Standard CDG-FE method

We first apply to (11) the standard CDG methods of Liu et al. ([32]) for space discretization and the forward Euler method for time discretization. That is, we look for ${\bf{U}}^{n+1,C}=(h^{n+1,C},(hP)^{n+1,C},(hQ)^{n+1,C})^{\top}\in\mathcal{V}^{C,k}$ such that for any ${\bf{V}}\in\mathcal{V}^{C,k}|_{C_{ij}}$ with any $i$ and $j$,

[TABLE]

Here $\theta=\Delta t_{n}/\tau\in[0,1]$ with $\tau$ being the maximal time step allowed by the CFL restriction ([32]). In general, this numerical scheme does not necessarily maintain the still-water stationary solution (8) and preserve the non-negativity of water depth.

Once $\mathbf{U}^{n+1,C}$ is available, we can obtain $u^{n+1,C}$ and $v^{n+1,C}$by applying a continuous finite element method to (9) and (10): look for $u^{n+1,C},v^{n+1,C}\in\widetilde{W}^{C,k}$ such that for any $\widehat{u},\widehat{v}\in\widehat{W}^{C,k}$,

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

and $\widetilde{W}^{C,k}$, $\widehat{W}^{C,k}$ are variants of $\mathcal{W}^{C,k}$ with consideration of the boundary conditions ([28]).

Remark 3.1**.**

When the bottom is flat ($b=\mbox{constant}$), the unique solvability of this FE method can be obtained in a straightforward manner if $h\geq h_{0}>0$ and $\alpha\geq\alpha_{0}>0$. When the bottom is not flat, the unique solvability of this FE scheme is more difficult to be determined. Even though $h\geq h_{0}>0$ and $\alpha\geq\alpha_{0}>0$, a constraint condition for the $b$ as in the 1D case in [28] is needed for the unique solvability of the FE method. Due to the complexity of the 2D case, the condition can not be obtained explicitly. However, except for the tests in Sections 4.2(Case B) and 4.3 which involve a dry area or a near dry area, the FE system of the other numerical applications in Section 4 (satisfy $h\geq h_{0}>0$ and $\alpha\geq\alpha_{0}>0$) was found to be uniquely solvable in numerics.

Remark 3.2**.**

For the cases in Sections 4.2(Case B) and 4.3 which involve a dry area or a near dry area, the FE equation is ill-conditioned, so we shall evaluate the velocity $u$ using an approximation technique. We solve the FE equation on the cells with $h^{n+1,C}\geq h_{0}>0$, while for other cells, we omit the high order terms in (9)-(10) and get

[TABLE]

Then the velocity $u$ and $v$ can be evaluated by

[TABLE]

To deal with the singularity in (22)-(23), we employ the regularization technique presented in [22] and thus we have:

[TABLE]

where

[TABLE]

with $\varepsilon=\min((\Delta x)^{4},(\Delta y)^{4})$. In the computation, $h_{0}=\max((\Delta x)^{k+1},(\Delta y)^{k+1})$ for $k=1,2$.

3.2.2 Well-balanced CDG-FE method

In this subsection, we propose a family of high order well-balanced CDG-FE schemes for the model (9)-(11), which exactly preserves the still-water steady state solution (8). The well-balance property can be achieved by adding some terms to the scheme (15),

[TABLE]

where the correction terms are given by

[TABLE]

[TABLE]

[TABLE]

Here, $\gamma^{n,D}_{ij}$ is a special constant which represents the average value of the water surface $\eta_{h}^{n,D}=h_{h}^{n,D}+b_{h}^{D}$ in the element $C_{ij}$. Particularly, we can take

[TABLE]

With the following decomposition of the source term

[TABLE]

the scheme (3.2.2) can be rewritten as

[TABLE]

Proposition 3.1**.**

The numerical scheme, defined in (31), (16) and (17) and their counterparts for ${\bf{U}}_{h}^{n+1,D}$, $u_{h}^{n+1,D}$ and $v_{h}^{n+1,D}$ , to solve the 2D Green-Naghdi model (9)-(11) is well-balanced, in the sense that it preserves the still-water stationary solution (8).

Proof.

The proof can be easily obtained by the mathematical induction with time step $n$ as the Proposition in [28]. ∎

3.2.3 Positivity-preserving CDG-FE method

We now discuss the positivity-preserving CDG-FE method for (9)-(11). Firstly, let $\hat{L}_{i}^{1,x}=\{\hat{x}_{i}^{1,\beta},\beta=1,2,...,{\hat{N}}\}$ and $\hat{L}_{i}^{2,x}=\{\hat{x}_{i}^{2,\beta},\beta=1,2,...,{\hat{N}}\}$ be the Legendre Gauss-Lobatto quadrature points on $[x_{i-\frac{1}{2}},x_{i}]$ and $[x_{i},x_{i+\frac{1}{2}}]$ respectively, while $\hat{L}_{j}^{1,y}=\{\hat{y}_{j}^{1,\beta},\beta=1,2,...,{\hat{N}}\}$ and $\hat{L}_{j}^{2,y}=\{\hat{y}_{j}^{2,\beta},\beta=1,2,...,{\hat{N}}\}$ represent the Legendre Gauss-Lobatto quadrature points on $[y_{j-\frac{1}{2}},y_{j}]$ and $[y_{j},y_{j+\frac{1}{2}}]$ respectively, $\forall i,j$. The corresponding quadrature weights on the reference element $[-\frac{1}{2},\frac{1}{2}]$ are $\hat{\omega}_{\beta},\beta=1,2,...,{\hat{N}}$, and $\hat{N}$ is chosen such that $2\hat{N}-3\geq k$. In addition, let $L_{i}^{1,x}=\{x_{i}^{1,\alpha},\alpha=1,2,...,N\}$ and $L_{i}^{2,x}=\{x_{i}^{2,\alpha},\alpha=1,2,...,N\}$ denote the Gaussian quadrature points on $[x_{i-\frac{1}{2}},x_{i}]$ and $[x_{i},x_{i+\frac{1}{2}}]$ respectively, while $L_{j}^{1,y}=\{y_{j}^{1,\alpha},\alpha=1,2,...,N\}$ and $L_{j}^{2,y}=\{y_{j}^{2,\alpha},\alpha=1,2,...,N\}$ represent the Gaussian quadrature points on $[y_{j-\frac{1}{2}},y_{j}]$ and $[y_{j},y_{j+\frac{1}{2}}]$ respectively. The corresponding quadrature weights $\omega_{\alpha},\alpha=1,2,...,N$ are distributed on the interval $[-\frac{1}{2},\frac{1}{2}]$ and $N$ is chosen such that the Gaussian quadrature is exact for the integration of univariate polynomials of degree $2k+1$. Define $L_{i,j}^{l,m}=(L_{i}^{l,x}\otimes\hat{L}_{j}^{m,y})\cup(\hat{L}_{i}^{l,x}\otimes L_{j}^{m,y})$ with $l,m=1,2$. By taking the test function ${\bf{V}}=(\frac{1}{\Delta x\Delta y},0,0)^{\top}$, we get the equation satisfied by the cell average of $h^{n,C}$,

[TABLE]

where $\bar{h}_{ij}^{n,C}$ denotes the cell average of the CDG solution $h^{C}$ on $C_{ij}$ at time $t_{n}$.

In the numerical implementation, the definite integrals in the intervals $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ and $[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$ are evaluated by applying the Gaussian quadrature rule described above to each half of the interval (also see Section 3.2 in [6]), the scheme (32) becomes

[TABLE]

herein, we used $\hat{x}_{i}^{1,1}=x_{i-\frac{1}{2}},\hat{x}_{i}^{2,\hat{N}}=x_{i+\frac{1}{2}},\hat{y}_{j}^{1,1}=y_{j-\frac{1}{2}},\hat{y}_{j}^{2,\hat{N}}=y_{j+\frac{1}{2}}$. Now we have the following result.

Proposition 3.2**.**

For any given $n\geq 0$, we assume $\bar{h}_{ij}^{n,C}\geq 0$ and $\bar{h}_{ij}^{n,D}\geq 0$, $\forall i,j$. Consider the scheme in (33) and its counterpart for $\bar{h}_{ij}^{n+1,D}$, if $h^{C}(x,y,t_{n})\geq 0$ and $h^{D}(x,y,t_{n})\geq 0$, $\forall(x,y)\in L_{i,j}^{l,m},\forall i,j$ with $l,m=1,2$, then $\bar{h}_{ij}^{n+1,C}\geq 0$ and $\bar{h}_{ij}^{n+1,D}\geq 0$, $\forall i,j$, under the CFL condition

[TABLE]

where $\lambda_{x}=\Delta t_{n}/\Delta x$, $\lambda_{y}=\Delta t_{n}/\Delta y$, $a_{x}=\max(\|u^{n,C}\|_{\infty},\|u^{n,D}\|_{\infty})$, $a_{y}=\max(\|v_{h}^{n,C}\|_{\infty},\|v_{h}^{n,D}\|_{\infty})$.

Proof.

Since the numerical solution $h^{n,D}$ is a piecewise polynomial with degree $k$, the integral of $h^{n,D}$ in cell $C_{ij}$ in (33) is exactly evaluated in our numerical implementation. However, in order to discuss the non-negativity of the cell average $\bar{h}_{ij}^{n+1,C}$. We equivalently evaluate the integral using a combination of the Gauss quadrature rule and the Legendre Gauss-Lobatto quadrature rule as follows:

[TABLE]

Plugging (3.2.3) into (33), one obtains

[TABLE]

Here, we used $\hat{\omega}_{\hat{N}}=\hat{\omega}_{1}$. A few observations can be made. Firstly, $\bar{h}_{ij}^{n+1,C}$ is a linear combination of $\bar{h}_{ij}^{n,C}$, $h^{n,D}\left(x_{i}^{l,\alpha},\hat{y}_{j}^{m,\beta}\right)$ and $h^{n,D}\left(\hat{x}_{i}^{l,\beta},y_{j}^{m,\alpha}\right)$, $l,m=1,2,\alpha=1,2,...,N,\beta==1,2,...,\hat{N}$, which are all non-negative according to the conditions in this Proposition. Secondly, the CFL condition (34) and $\theta\in[0,1]$ imply that all coefficients in the linear combination are non-negative. Therefore $\bar{h}_{ij}^{n+1,C}\geq 0$, $\forall i,j$. Similarly, one can show $\bar{h}_{ij}^{n+1,D}\geq 0$, $\forall i,j$. ∎

Remark 3.3**.**

Although the line to prove Proposition 3.2 is similar to the one in [29], there is an explicit difference in the proof. In [29], the water depth $h$ and the momentum $hu$ may be not consistent due to the numerical error, namely, $hu\neq 0$ if $h\equiv 0$. This issue hampers the proof of the non-negativity of the water depth. To overcome the issue, a numerical technique has been employed to modify the momentum to be consistent with the water depth. While according to the numerical scheme in the present paper, the velocity is always equal to zero as the water depth is equal to zero. Therefore, the proof for Proposition 3.2 is more natural and simpler than the one in [29].

Next, we give the positivity-preserving limiters which modify the CDG solution polynomials $h^{n,C}$ and $h^{n,D}$ into $\tilde{h}^{n,C}$ and $\tilde{h}^{n,D}$ which satisfy the sufficient condition given in Proposition 3.2. In fact, the limiters are the same as in [44, 45, 41], as long as the notation $K$ and $\hat{L}_{K}$ are re-defined as follows: On the primal mesh, $K$ denotes a mesh element $C_{ij}$ and $\hat{L}_{K}$ represents the set of relevant quadrature points in $K$, namely $\hat{L}_{K}=\cup_{l,m=1}^{2}L_{i,j}^{l,m}$. On the dual mesh, $K$ denotes a mesh element $D_{ij}$ and $\hat{L}_{K}$ represents the set of relevant quadrature points in $K$, namely $\hat{L}_{K}=L_{i,j}^{1,1}\cup L_{i,j-1}^{1,2}\cup L_{i-1,j}^{2,1}\cup L_{i-1,j-1}^{2,2}$. Following [41], the positivity-preserving limiter is given as follows: On each mesh element $K$, we modify the water depth $h^{n,\star}$ into $\tilde{h}^{n,\star}=\alpha_{K}(h^{n,\star}-\bar{h}^{n,\star})+\bar{h}^{n,\star}$, with $\alpha_{K}=\min_{x\in\hat{L}_{K}}\left\{1,|\bar{h}^{n,\star}/(\bar{h}^{n,\star}-h^{n,\star}(x))|\right\}$ and $\star=C,D$.

3.2.4 Positivity-preserving well-balanced CDG-FE method

Finally, we study the positivity-preserving and well-balanced CDG-FE method for the 2D Green-Naghdi model (9)-(11). We start with the well-balanced CDG method satisfied by the cell average of the numerical solution $h^{C}$, which is obtained by taking the test function ${\bf{V}}=(\frac{1}{\Delta x\Delta y},0,0)^{\top}$ in the scheme (31),

[TABLE]

where $\bar{b}_{ij}^{C}$ (resp. $\bar{b}_{ij}^{D}$) is the cell average of the bottom topography $b^{C}$ (resp. $b^{D}$) on the element $C_{ij}$ (resp. $D_{ij}$).

Proposition 3.3**.**

For any given $n\geq 0$, we assume $\bar{h}_{ij}^{n,C}\geq 0$ and $\bar{h}_{ij}^{n,D}\geq 0$, $\forall i,j$. Consider the scheme in (37) and its counterpart for $\bar{h}_{ij}^{n+1,D}$, if

[TABLE]

and $h^{C}(x,y,t_{n})\geq 0$, $h^{D}(x,y,t_{n})\geq 0$, $\forall(x,y)\in L_{i,j}^{l,m}$, $\forall i,j$ with $l,m=1,2$, then $\bar{h}_{ij}^{n+1,C}\geq 0$ and $\bar{h}_{ij}^{n+1,D}\geq 0$, $\forall i,j$, under the CFL condition

[TABLE]

Proof.

The proof is a direct result of Proposition 3.2 and condition (38). ∎

Remark 3.4**.**

To enforce the sufficient conditions in Proposition 3.3, we also need to modify the approximations to the bottom and use the positivity-preserving limiter in Section 3.2.3. The new approximations to the bottom topography $b(x,y)$, still denoted by $b^{C}$ and $b^{D}$, are obtained by solving a constrained minimization problem similar to the one in [29] with the Lagrange multiplier method.

3.3 High-order time discretizations and nonlinear limiters

To achieve better accuracy in time, the strong stability preserving (SSP) high-order time discretizations ([18]) will be used in the numerical simulations . Such discretizations can be written as a convex combination of the forward Euler method, and therefore the resulting SSP schemes are also well-balanced and positivity-preserving. In this paper, we use the third order TVD Runge-Kutta method for the time discretization.

When the CDG method is applied to nonlinear problems, nonlinear limiters are often needed to prevent numerical instabilities. In this work, we use the total variation bounded (TVB) minmod slope limiter with parameter $M=10$ ([12]) in a componentwise way as it is needed. This limiter is applied to $(h+b,hP,hQ)^{\top}$ and it is used prior to application of the positivity-preserving limiter.

4 Numerical examples

4.1 Accuracy test

In this example, we test the convergence rate of the proposed CDG-FE method by varying the mesh size. The Green-Naghdi model (1) with $b=0$ has an exact solution given by ([38])

[TABLE]

where $h_{1}$ is the typical water depth, $h_{2}$ corresponds to the solitary wave crest and $D=\sqrt{gh_{2}}$ is the wave speed. In this test, we employ a solitary wave with $h_{1}=1$ and $h_{2}=2.25$ in (40) which is initially located at $x=0$ and propagating in the positive $x$-direction. The computational domain is $[-30,50]\times[-1,1]$ and the final time is $1$. An outgoing boundary condition is used in the $x$-direction and a periodic boundary condition is used in the $y$-direction. We use regular meshes with $\Delta x=\Delta y=1,0.5,0.25,0.125$. The time step is $\Delta t=0.1\Delta x$. We present $L^{2}$ errors and orders of accuracy for $h$ and $u$ in Table 1. The results show that the CDG-FE method is $(k+1)$st order accurate for $P^{k}$ with $k=1,2$ and therefore it is optimal with respect to the approximation properties of the discrete spaces.

4.2 Stationary solution

In this test, we validate the well-balanced feature and the positivity-preserving property of the proposed method as applied to continuous and discontinuous variable bottoms. The initial conditions are

[TABLE]

and the continuous bottom profile (Case A) is defined by

[TABLE]

with $r=\sqrt{x^{2}+y^{2}}$, while the discontinuous bottom profile (Case B) is given by

[TABLE]

We choose $[-1,1]\times[-1,1]$ as the computational domain, divided into $20\times 20$ elements, and use outgoing boundary conditions. We compute the solution up to $t=10$ by the well-balanced CDG methods. Notice that there exists a near dry area for the second case.

For these cases, the standard CDG-FE method usually fails to preserve the still-water stationary solution exactly. Especially for the second case, the standard CDG-FE method will produce negative water depth due to the numerical oscillation, and thus the computation will blow down. To demonstrate that the positivity-preserving well-balanced CDG-FE scheme indeed preserves the still-water stationary solution exactly (i.e., up to machine precision), we perform the computation in both single and double precision. The corresponding $L^{2}$ errors on the water surface $h+b$ and velocity $(u,v)$ are listed in Table 2 for both topographies. We see that the errors have orders of a magnitude consistent with the machine single and double precision, and thus the numerical results verify the well-balanced property and the positivity-preserving property.

4.3 Solitary wave overtopping a seawall

In this test, we consider the simulation of a solitary wave overtopping a seawall, which has been studied experimentally in [21] and numerically in [23]. This example is also used to investigate the validity of our positivity-preserving well-balanced scheme. The initial solitary wave ($h_{1}=0.2,h_{2}=0.27$ in (40)) and the bottom topography including the seawall are shown in Figure 2 (along $y=0$). The computational domain is $[-5,20]\times[-0.2,0.2]$ discretized into $500\times 8$ uniform elements. The outgoing boundary condition is used in the $x$-direction and the periodic boundary condition is used in the $y$-direction. For the standard CDG-FE method, negative water depth was generated during the simulation at the dry or near dry areas, this causes inaccurate velocity which is used as the boundary condition when solving the elliptic equations, and then the numerical solution blows up within a few dozen time steps. For the positivity-preserving well-balanced CDG-FE method, however, the water depth remains non-negative during the entire simulation. The numerical surface profiles (along $y=0$) at time $t=5,7.5,12.5,20$ are plotted in Figure 3. The time series of the wave elevation at several positions of gauges ($x=5.9,7.6,9.644,10.462,10.732$ and $11.12$, $y=0$) are illustrated in Figure 4. The time origin has been shifted in order to compare with results reported in [21]. Both models with $\alpha=1$ and $\alpha=1.159$ give almost the same results. Our numerical results are similar to experimental data and the numerical results in [21].

4.4 Periodic waves propagation over a submerged bar

In this example, we investigate the robustness of the modified Green-Naghdi model. We first consider the propagation of periodic Stokes waves over a submerged bar with plane slopes. The bottom variation is specified by

[TABLE]

and is also exhibited in Figure 5 (along $y=0$) in which we also label the positions of $10$ gauges used in [14].

As shown in experimental work ([14]), regular waves break up into higher-frequency free waves as they propagate past a submerged bar. As the waves travel up the front slope of the bar, higher harmonics are generated due to nonlinear interactions, causing the waves to steepen. These harmonics are then released as free waves on the downslope, producing an irregular pattern behind the bar. This experiment is particularly difficult to simulate because it includes nonlinear interactions and requires accurate propagation of waves in both deep and shallow water over a wide range of depths. Therefore it has often been used as a discriminating test case for nonlinear models of surface wave propagation over variable bottom ([3, 14, 20]).

In the simulation, the computational domain is $[0,25]\times[-0.2,0.2]$, divided into $500\times 8$ uniform cells. At initial time, $h+b=0$ and $u=v=0$ in the computational domain. The incident wave (entering from the left) is a third-order Stokes wave ([17]) given by

[TABLE]

where $T_{0}$, $a_{0}$ and $\lambda$ denote the wave period, amplitude and wavelength, respectively. We choose $(T_{0},a_{0},\lambda)=(2.02,0.01,3.73)$ corresponding to one of the experiments in [14]. An absorbing boundary condition is applied at the right boundary and a periodic boundary condition is used at the upper and bottom boundaries.

Figure 6 depicts the time histories of the water surface at the first $6$ gauges ($x=2$, $4$, $10.5$, $12.5$, $13.5$ and $14.5$, $y=0$) and Figure 7 depicts the time histories of the water surface at the last $4$ gauges ($x=15.7$, $17.3$, $19$ and $21$, $y=0$). The time origin has been shifted in order that the numerical results match the measurements for the first gauge at $x=2,y=0$. We compare three sets of data: experimental data ([14]), numerical solutions from the original Green-Naghdi model ($\alpha=1$), numerical solutions from the improved Green-Naghdi model ($\alpha=1.159$). It can be seen from Figure 6 that before the crest of the bar ($x\leq 14.5$), the numerical results from both models ($\alpha=1$ and $\alpha=1.159$) match well with each other and compare well with experimental data. However, it can be seen from Figure 7 that, compared with the experimental data, discrepancies in amplitude and phase can be observed for gauges beyond the crest of the bar ($x\geq 15.7$) for the numerical solutions with $\alpha=1$, these discrepancies should be attributed to the weakly dispersive character of the original Green-Naghdi model, while for the modified Green-Naghdi model, the numerical solutions match well with the experimental data.

In this test, we further study the dispersive effect of the modified Green-Naghdi model. We consider a periodic wave propagation over a submerged bar with elliptic slope (see Figure 8), which is also given by

[TABLE]

with $r=\frac{(x-12.5)^{2}}{100}+\frac{y^{2}}{16}$.

The computational domain is $[0,25]\times[-1,1]$, which is discretized into $125\times 20$ uniform elements. At initial time, $h+b=0$ and $u=v=0$ in the computational domain. The incident wave (entering from the left) is a third-order Stokes wave given in (45) with $(T_{0},a_{0},\lambda)=(3,0.01,3.73)$. Solid wall boundary conditions are used at the top and bottom boundaries, an absorbing boundary condition is used at right boundary. The numerical free surfaces at $t=30$ for both models ($\alpha=1$ and $\alpha=1.159$) are shown in Figure 9. Overall, both numerical surfaces compare well with each other. To observe the discrepancies more clearly, the time series of numerical surfaces at several positions ($A_{1}(8,0)$, $A_{2}(9,0.5)$, $A_{3}(21,0)$, $A_{4}(18,0)$, $A_{5}(19,0.5)$ and $A_{6}(22,0.5)$) are shown in Figure 10. It can be seen that both numerical surfaces compare well with each other for gauges before the crest of the bar ($x\leq 9$). However, the discrepancies in amplitude and phase can be observed for gauges beyond the crest of the bar ($x\geq 18$). These discrepancies again should be attributed to the dispersive character of the models.

4.5 Solitary wave propagation over a composite beach

To further investigate the robustness of the modified Green-Naghdi model, we simulate the propagation of solitary waves over a composite beach, which consists of three piece-wise linear segments, terminated with a vertical wall on the left. The slopes of the topography are defined as follows ([23]):

[TABLE]

We consider the propagation of a solitary wave: $h_{1}=0.22,h_{2}=1.73h_{1}$. The initial solitary waves, produced by using (40), are located at $x=0$ and propagate to the right. The computational domain is $[-5,23.23]\times[-0.2,0.2]$ discretized into $500\times 8$ uniform elements. The bottom along with $y=0$ and the initial solitary wave is shown in Figure 11. The outgoing boundary condition is used at the left boundary, the reflective boundary condition is employed at the right boundary and the periodic boundary condition is used at the upper and bottom boundaries. We observe the propagation over the beach, reflection on the vertical wall before traveling back to the left boundary. The time series of the wave elevation at several positions of gauges ($x=15.04,19.4,22.33$, $y=0$) are shown in Figure 12. Overall, both sets of the numerical results obtained from the Green-Naghdi models ($\alpha=1$ and $\alpha=1.159$) compare well together before the wave encounters the vertical wall. However, discrepancies in amplitude and phase can be observed after the wave reflects on the vertical wall. Especially, discrepancies are clearer when the positions of gauges are far away from the vertical wall. Our numerical solution with $\alpha=1.159$ matches well with the one reported in [23].

5 Conclusions

In this work, we derive a Green-Naghdi model with enhanced dispersive property and then develop a family of high order positivity-preserving and well-balanced numerical methods for its numerical solutions. These methods are based on the reformulation of the original system into a pseudo-conservation law coupled with an elliptic system. Numerical experiments are presented to demonstrate the accuracy and achievement of expected properties for the proposed schemes, and the capability of the Green-Naghdi equations to model a wide range of shallow water wave phenomena. Application of the proposed schemes to Green-Naghdi models taking into account the irrotational effect and the Coriolis effect due to the Earth’s rotation, and proposal of fast solvers for the elliptic part of the model will be envisioned for our future work.

Acknowledgments

Maojun Li is partially supported by NSFC (Grant Nos. 11501062, 11701055, 11871139). Liwei Xu is partially supported by a Key Project of the Major Research Plan of NSFC (Grant No. 91630205) and a NSFC (Grant No. 11771068).

Appendix

Here, we derive the fully nonlinear and weakly dispersive shallow water equations over the non-flat bottom in 2D space (1) along the same lines in [38] where the authors derived the one-dimensional equations on the flat bottom. Let $\Omega(t)$ be the domain in $R^{3}$ occupied by the water at time $t$. The propagation of the water is described by the fully nonlinear Euler equations

[TABLE]

and the continuity equation

[TABLE]

where $(\tilde{u},\tilde{v},\tilde{w})$ denotes the velocity of the water, $p$ is the pressure and $g$ is the gravitational constant. The subscript $t$ denotes the partial derivative with respect to the time variable, $x$, $y$, and $z$ denote the partial derivatives with respect to the space variables. We assume that the density is taken as one. The boundary conditions are given by

(BC1) the kinematic condition at the free surface

[TABLE]

(BC2) the impermeability of the bottom

[TABLE]

where $h(x,y,t)$ denotes the depth of the water and $b(x,y)$ is the bottom topography. The superscript $(s)$ and $(b)$ denote the quantities evaluated at the free surface and the bottom, respectively.

Let $u(x,y,t)$ and $v(x,y,t)$ denote the vertically averaged horizontal velocity in the x- and y-directions, respectively, and be defined by

[TABLE]

Along the same lines as in [38], we want to find a set of equations which governs the evolution of the water depth $h(x,y,t)$, the average velocity $u(x,y,t)$ and $v(x,y,t)$ under the following assumptions:

(A1) we assume that the vertical movement of a particle is small compared with the horizontal movement, that is, we use the shallow water hypothesis, so that we can write

[TABLE]

(A2) the dynamic condition at the free surface is assumed to be

[TABLE]

Integrating (51) with respect to $z$ gives

[TABLE]

Here, we have used (53). Equation (58) also implies

[TABLE]

Multiplying (54) by $h$ and taking the derivative with respect to $x$, multiplying (55) by $h$ and taking the derivative with respect to $y$, and then adding them together, one arrives at

[TABLE]

Utilizing (52), (53) and (59), one gets from (60)

[TABLE]

This equation gives the evolution of $h$ provided we know the evolution of $hu$ and $hv$.

Integrating (48) with respect to $z$ from $z=b$ to $z=h+b$, performing an integration by parts for the fourth term and utilizing (51), (52) and (53), one obtains

[TABLE]

Similarly, one has from (49)

[TABLE]

Integrating (50) with respect to $z$ from $z$ to $h+b$ yields

[TABLE]

thus one gets

[TABLE]

Here, we have used (58) and (51).

Now, substituting (65) into (62) and (63), respectively, and using the assumptions (56)-(57), one obtains

[TABLE]

and

[TABLE]

where

[TABLE]

Therefore, we get a set of equations by combining (61), (66) and (67)

[TABLE]