An Efficient Two-Layer Non-hydrostatic Approach for Dispersive Water Waves

Abstract

In this paper, we propose a two-layer depth-integrated non-hydrostatic system with improved dispersion relations. This improvement is obtained through three free parameters: two of them related to the representation of the pressure at the interface and a third one that controls the relative position of the interface concerning the total height. These parameters are then optimized to improve the dispersive properties of the resulting system. The optimized model shows good linear wave characteristics up to $kH\approx 10$, that can be improved for long waves. The system is solved using an efficient formally second-order well-balanced and positive preserving hybrid finite volume/difference numerical scheme. The scheme consists of a two-step algorithm based on a projection-correction type scheme. First, the hyperbolic part of the system is discretized using a Polynomial Viscosity Matrix path-conservative finite-volume method. Second, the dispersive terms are solved using finite differences. The method has been applied to idealized and challenging physical situations that involve nearshore breaking. Agreement with laboratory data is excellent. This technique results in an accurate and efficient method.

Appendices

A Hyperbolicity of the Underlying Hydrostatic System

1.1 A.1 Hyperbolicity

Let us consider the underlying hydrostatic system (20) written in quasi-linear form

$$\begin{aligned} \partial _t \varvec{U}+\mathcal {A}\partial _x \varvec{U} = \varvec{G}(\varvec{U})\partial _x H, \end{aligned}$$

(47)

where $\mathcal {A} = \mathcal {J}_F + \varvec{B}$, being $\mathcal {J}_F$ the Jacobian matrix of the flux $\varvec{F}$. The characteristic polynomial of the matrix $\mathcal {A}$ is given by

$$\begin{aligned} \mathcal {P} (\lambda ,l_1) = \displaystyle \frac{1}{2}\left( u_1 - \lambda \right) \left( u_2 - \lambda \right) \mathcal {Q}(\lambda ,l_1), \end{aligned}$$

being $\mathcal {Q}(\lambda ,l_1)$ a third order polynomial on $\lambda$ given by

$$\begin{aligned} \mathcal {Q}(\lambda ,l_1)= & {} f(\lambda ) - R(l_1),\\ f(\lambda )= & {} \left( 3u_1-u_2-2\lambda \right) \left( (u_2-\lambda )^2-gh\right) ,\\ R(l_1)= & {} l_1K,\\ K= & {} \left( u_1-u_2\right) \left( \left( u_1-u_2\right) ^2-4gh\right) . \end{aligned}$$

For the sake of simplicity on the notation, we do not write explicitly the dependence on $\varvec{U}$. Let us study the hyperbolicity of the hydrostatic system. It is easily to check that

$$\begin{aligned} \lambda _1=u_1,\ \lambda _2 = u_2 \end{aligned}$$

are eigenvalues of the system for every $l_1\in (0,1)$. It remains to check if the cubic polynomial $\mathcal {Q}(\lambda )$ has three distinct roots.

We will prove that the cubic polynomial it has always three different roots for every $l_1\in [0,1],$ and in particular for every $l_1\in (0,1)$ as we just request. Let us give a sketch of the proof:

• 1. Let us remark that $\mathcal {Q}$ is a cubic polynomial on $\lambda$ satisfying

  $$\begin{aligned} \mathcal {Q}(-\infty ,l_1) = +\infty ,\ \mathcal {Q}(\infty ,l_1) = -\infty . \end{aligned}$$

• 2. Note that $f(\lambda )$ is a cubic polynomial that does not depend on $l_1$. Moreover, it has two local extrema given by the roots of $f'(\lambda )$:

  $$\begin{aligned}\lambda ^\pm = \displaystyle \frac{u_1+u_2}{2} \pm \sqrt{\frac{1}{3} + \left( \displaystyle \frac{u_1-u_2}{2} \right) ^2}.\end{aligned}$$

• 3. A sufficient and necessary condition for the existence of three real and distinct roots of the cubic polynomial $\mathcal {Q}(\lambda ,l_1)$ is that:

  $$\begin{aligned} f(\lambda ^-) < R(l_1) \text{ and } f(\lambda ^+) > R(l_1). \end{aligned}$$

• 4. Note that when $l_1=0$, the polynomial $\mathcal {Q}(\lambda ,0)$ has three roots:

  $$\begin{aligned} \lambda _3 = \displaystyle \frac{3u_1-u_2}{2},\ \lambda _{4,5} = u_2\pm \sqrt{gh}, \end{aligned}$$

  and thus $f(\lambda ^ -) < R(0)$ and $f(\lambda ^ +) > R(0)$. Similarly, when $l_1=1$, the polynomial $\mathcal {Q}(\lambda ,1)$ has three roots:

  $$\begin{aligned} \lambda _3 = \displaystyle \frac{-u_1+3u_2}{2},\ \lambda _{4,5} = u_1\pm \sqrt{gh}, \end{aligned}$$

  and therefore $f(\lambda ^ -) < R(1)$ and $f(\lambda ^ +) > R(1)$.

Thus, if we assume that $K\ge 0,$ then $R(1) \ge R(l_1) \ge R(0)$ and therefore

$$\begin{aligned} f(\lambda ^ +)> & {} R(1) \ge R(l_1) \\ f(\lambda ^ -)< & {} R(0) \le R(l_1). \end{aligned}$$

If we assume that $K\le 0,$ then $R(1) \le R(l_1) \le R(0)$ and therefore

$$\begin{aligned} f(\lambda ^ +)> & {} R(0) \ge R(l_1) \\ f(\lambda ^ -)< & {} R(1) \le R(l_1). \end{aligned}$$

This concludes the proof.

1.2 A.2 A First Order Approximation for the Eigenvalues

Let us denote the eigenvalues that depends on $l_1$ as

$$\begin{aligned} \lambda _1(l_1)=u_1,\ \lambda _2(l_1)=u_2, \end{aligned}$$

as the known eigenvalues for any $l_1\in (0,1)$, and

$$\begin{aligned} \lambda _3(l_1),\ \lambda _4(l_1),\ \lambda _5(l_1) \end{aligned}$$

as the eigenvalues that are roots of the cubic polynomial $\mathcal {Q}(\lambda ,l_1)$. As a particular case, we have found an explicit form of the eigenvalues for $l_1=1/2,$

$$\begin{aligned} \lambda _1 = u_1,\ \lambda _2 = u_2,\ \lambda _3 = \displaystyle \frac{u_1+u_2}{2},\ \lambda _{4,5} = \displaystyle \frac{u_1+u_2}{2}\pm \sqrt{gh + \displaystyle \frac{3}{4}(u_1-u_2)^2}. \end{aligned}$$

Let us consider

$$\begin{aligned} \lambda (l_1)\in \{\lambda _3(l_1),\ \lambda _4(l_1),\ \lambda _5(l_1)\}, \end{aligned}$$

an eigenvalue that depends on $l_1$ and is a root of the cubic polynomial $\mathcal {Q}(\lambda ,l_1)$. We propose the following approximation of the eigenvalues, that gives the exact roots of the cubic polynomial $\mathcal {Q}(\lambda , l_1)$ for $l_1\in {\{0,1/2,1\}}$.

$$\begin{aligned} \lambda _3(l_1) \approx \widetilde{\lambda _3}(l_1)= & {} \left( \frac{3}{2} - 2l_1\right) u_1 + \left( 2l_1-\frac{1}{2}\right) u_2, \\ \lambda _{4,5}(l_1) \approx \widetilde{\lambda _{4,5}}(l_1)= & {} l_1u_1+l_2u_2 \pm \sqrt{gh + 3l_1l_2\left( u_1-u_2\right) ^2}. \end{aligned}$$

Another approximation for the eigenvalues is proposed in the following. Since $\lambda (l_1)$ is a root of $\mathcal {Q}(\lambda ,l_1),$ then

$$\begin{aligned} \mathcal {Q}(\lambda (l_1),l_1)=0, \end{aligned}$$

and deriving with respect to $l_1$ it yields

$$\begin{aligned} \lambda '(l_1) = \displaystyle \frac{\left( \left( u_1 - u_2\right) ^2-4 g h\right) \left( u_1 - u_2\right) }{2 \left( g h - 3 u_1 u_2 + 3 \left( u_1 + u_2 - \lambda (l_1)\right) \lambda (l_1)\right) }. \end{aligned}$$

Thus, we propose to approximate the eigenvalues that are roots of $\mathcal {Q}(\lambda ,l_1)$ with the first order approximation

$$\begin{aligned} \lambda _i(l_1) \approx \widetilde{\lambda _{i}} = \lambda _i(1/2) + \lambda _i'(1/2)(l_1-1/2), \ i\in \{3,4,5\}, \end{aligned}$$

that can be explicitly computed, since $\lambda _i(1/2)$ are known:

$$\begin{aligned} \lambda _3(l_1) \approx \widetilde{\lambda _{3}}= & {} \displaystyle \frac{u_1+u_2}{2} + \left( u_1 - u_2\right) \left( 1-2 l_1 \right) \frac{g h-\frac{1}{4}\left( u_1 - u_2\right) ^2}{gh + \frac{3}{4}\left( u_1-u_2 \right) ^2} + \mathcal {O}(l_1^2),\\ \lambda _{4,5}(l_1) \approx \widetilde{\lambda _{4,5}}= & {} \displaystyle \frac{u_1+u_2}{2} \pm \sqrt{gh + \displaystyle \frac{3}{4}(u_1-u_2)^2}\nonumber \\&+ \left( u_1 - u_2\right) \left( l_1 -\frac{1}{2} \right) \frac{g h-\frac{1}{4}\left( u_1 - u_2\right) ^2}{gh + \frac{3}{4}\left( u_1-u_2 \right) ^2} + \mathcal {O}(l_1^2). \end{aligned}$$

This procedure is more rigorous and lead to a more sophisticated expressions of the approximated eigenvalues.

B Linear Dispersion Properties

1.1 B.1 Linear Dispersion Relation

Substituting (23) into (24) yields the linear dispersion relation:

$$\begin{aligned} \displaystyle \frac{C^2}{g H} = \frac{N_0 + N_1 \left( kH \right) ^2 }{D_0 + D_1 \left( kH \right) ^2 + D_2 \left( kH \right) ^4}, \end{aligned}$$

(48)

where

$$\begin{aligned} {\left\{ \begin{array}{ll} N_0 = 1,\quad N_1 = \displaystyle \frac{l_1 l_2 (-\gamma _1-\gamma _2+2 (\gamma _2-1) l_1+2)}{4 (\gamma _1+\gamma _2)}, \\ \\ D_0 = 1,\qquad D_1 = \displaystyle \frac{{\gamma 1}+ {\gamma 2}+2 ({\gamma 2}-2) l_1^2-2 l_1 ( {\gamma 1}+{\gamma 2}-2)}{4 ({\gamma 1}+ {\gamma 2})},\quad D_2 = \displaystyle \frac{l_1^2l_2^2 ({\gamma 1}- {\gamma 2})}{16 ({\gamma 1}+{\gamma 2})}. \end{array}\right. } \end{aligned}$$

C Breaking Waves Parameters

By taking into account the two vertical velocities equations in (45) and the incompressibility equations, which relates $w_\alpha$ with $u_\alpha$, lead us to write P in terms of the derivatives of U

$$\begin{aligned} A_{(xx)} \partial _{xx} \varvec{U} + A_{(xt)} \partial _{xt} \varvec{U} + A_{(x)} \partial _{x} \varvec{U} + A_{(t)} \partial _{t} \varvec{U} + A \varvec{U} + B = \begin{pmatrix} p_b-p_I \\ \\ \gamma _1 p_b + \gamma _2 p_I \end{pmatrix} +I_{(\varsigma )} \partial _x \varvec{U}, \end{aligned}$$

where

$$\begin{aligned} I_{(\varsigma )} = \begin{pmatrix} -\varsigma _1 l_1h&{} 0\\ \\ 0&{} -\varsigma _2 l_2h \end{pmatrix},\ A_{(\cdot )}\in \mathcal {M}_2(\mathbb {R}). \end{aligned}$$

We propose define $\varsigma _\alpha$ such that

$$\begin{aligned} I_{(\varsigma )}:=Diag(A_{(x)}). \end{aligned}$$

We then proceed to compute $A_{(x)}$. The two continuity equations can be written as

Neglecting mass transfer terms due to $\varGamma _I$, the vertical equations can be written as

$$\begin{aligned} {\left\{ \begin{array}{ll} \partial _t \left( l_1h w_1 \right) + \partial _x \left( l_1hu_1w_1 \right) = p_b-p_I + \varsigma _1 l_1h \partial _x u_1, \\ \\ \partial _t \left( l_2h w_2 \right) + \partial _x \left( l_2hu_2w_2 \right) = \gamma _1 p_b+\gamma _2 p_I+ \varsigma _2 l_2h \partial _x u_2. \end{array}\right. } \end{aligned}$$

(50)

Then, retaining at the left hand side of Eq. (50) only the terms multiplied by $\partial _x u_\alpha$ at the equation concerning to the layer $\alpha$, leads to

$$\begin{aligned} I_{(\varsigma )} = \begin{pmatrix} l_1h\left( -\partial _t l_1h + w_1 - u_1 \partial _x \left( l_1h + H \right) \right) &{} 0 \\ \\ 0 &{} l_2h\left( -\partial _t l_2h + w_2 + u_2 \partial _x \left( z_I- l_2h \right) \right) \end{pmatrix}, \end{aligned}$$

Again, using that $\varGamma _I=0$, $I_{(\varsigma )}$ can be rewritten as

$$\begin{aligned} I_{(\varsigma )} = \begin{pmatrix} l_1h\left( -w_1 + 3 u_1 \partial _x H - 2 \partial _t H \right) &{} 0 \\ \\ 0 &{} l_2h\left( -w_2 - 2 \partial _t H - 3 u_2 \partial _x z_I \right) \end{pmatrix}, \end{aligned}$$

and finally we propose

$$\begin{aligned} {\left\{ \begin{array}{ll} \varsigma _1 = \left( w_1 - 3 u_1 \partial _x H + 2 \partial _t H \right) , \\ \\ \varsigma _2 = \left( w_2 + 3 u_2 \partial _x z_I+ 2 \partial _t H \right) . \end{array}\right. } \end{aligned}$$

(51)

D Coefficients and Matrices of the Linear System

1.1 D.1 Coefficients of the Poisson-Like Equations

The coefficients appearing in (37) and (38) are:

$$\begin{aligned}&\left\{ \begin{array}{l} a_1 = -(l_1h)^2,\\ \\ a_2 = -l_1^2h \partial _x h,\\ \\ a_3 = l_1h\partial _{xx} \left( l_1h+2l_2h-2\eta \right) + \partial _x \left( h_1+2l_2h-2\eta \right) ^2 + 4,\\ \\ a_4 = -(l_1h)^2,\\ \\ a_5 = -l_1h \partial _x \left( 3l_1h+4l_2h-4\eta \right) ,\\ \\ a_6 = -l_1h\partial _{xx} \left( l_1h+2l_2h-2\eta \right) - \partial _x \left( l_1h+2l_2h-2\eta \right) ^2 - 4, \end{array} \right. \end{aligned}$$

(52)

$$\begin{aligned}&\left\{ \begin{array}{l} b_1 = -(l_2h)^2\left( \gamma _1+2\displaystyle \frac{l_1}{l_2} \right) ,\\ \\ b_2 = l_2h \partial _x \left( (4-\gamma _1)l_2h-4\eta \right) ,\\ \\ b_3 = l_2h\partial _{xx}\left( 2l_1h+(2\gamma _1+3)l_2h-2(\gamma _1+2)\eta \right) + \gamma _1\partial _x \left( l_2h-2\eta \right) + 4\gamma _1,\\ \\ b_4 = -(l_2h)^2\left( \gamma _2+2\displaystyle \frac{l_1}{l_2} \right) ,\\ \\ b_5 = l_2h \partial _x \left( 4l_1h + (4-\gamma _2)l_2h - 4\eta \right) ,\\ \\ b_6 = l_2h\partial _{xx}\left( -2l_1h+(2\gamma _2-5)l_2h-2(\gamma _2-2)\eta \right) + \gamma _2\partial _x \left( l_2h-2\eta \right) + 4\gamma _2, \end{array} \right. \nonumber \\&RHS_1 = l_1h^{(\widetilde{k})} \partial _x q_{u,1}^{(\widetilde{k})} - 2 q_{u,1}^{(\widetilde{k})}\partial _x z_1^{(\widetilde{k})} + 2q_{w,1}^{(\widetilde{k})} + 2 h^{(\widetilde{k})} \partial _t H,\nonumber \\&RHS_2 = 2 l_1h^{(\widetilde{k})} \partial _x q_{u,1}^{(\widetilde{k})} + l_2h^{(\widetilde{k})} \partial _x q_{u,2}^{(\widetilde{k})} - 2 q_{u,2}^{(\widetilde{k})}\partial _x z_2^{(\widetilde{k})} + 2q_{w,2}^{(\widetilde{k})} + 2 h^{(\widetilde{k})} \partial _t H. \end{aligned}$$

(53)

1.2 D.2 Matrices of the Linear Systems

After replace (39) and (40) in (37) and (38), one has to solve a linear system

being $T_{(j)}, C_{(j)}$ tridiagonal matrices of dimension $N\times N$ given by:

$$\begin{aligned} T_{(1)}= & {} \displaystyle \frac{A_{(1)}}{\varDelta x^2} T_{(1,-2,1)} +\displaystyle \frac{A_{(2)}}{2\varDelta x} T_{(-1,0,1)} + A_{(3)}I,\\ C_{(1)}= & {} \displaystyle \frac{A_{(4)}}{\varDelta x^2} T_{(1,-2,1)} +\displaystyle \frac{A_{(5)}}{2\varDelta x} T_{(-1,0,1)} + A_{(6)}I,\\ T_{(2)}= & {} \displaystyle \frac{B_{(1)}}{\varDelta x^2} T_{(1,-2,1)} +\displaystyle \frac{B_{(2)}}{2\varDelta x} T_{(-1,0,1)} + B_{(3)}I,\\ C_{(2)}= & {} \displaystyle \frac{B_{(4)}}{\varDelta x^2} T_{(1,-2,1)} +\displaystyle \frac{B_{(5)}}{2\varDelta x} T_{(-1,0,1)} + B_{(6)}I, \end{aligned}$$

where

$$\begin{aligned} T_{(a,b,c)} = \begin{pmatrix} {b} &{} \quad {c} &{} \quad { } &{} \quad { } &{} \quad { 0 } \\ {a} &{} \quad {b} &{} \quad {c} &{} \quad { } &{} \quad { } \\ { } &{} \quad {\ddots } &{} \quad {\ddots } &{} \quad \ddots &{} \quad { } \\ { } &{} \quad { } &{} \quad a &{} \quad b &{} \quad {c}\\ { 0 } &{} \quad { } &{} \quad { } &{} \quad {a} &{} \quad {b}\\ \end{pmatrix}, \end{aligned}$$

gather the centred finite difference matrix of second ($T_{(1,-\,2,1)}$) and first ($T_{(-1,0,1)}$) order, and I the identity matrix of dimension $2N\times 2N$.

The matrices $A_{(j)}$ and $B_{(j)},\ j\in {\{1,\ldots , 6\}}$ are diagonal matrices of dimension $N\times N$

$$\begin{aligned}A_{(j)} = \begin{pmatrix} {a_{j,1}} &{} \quad { } &{} \quad { } &{} \quad { } &{} \quad { 0 } \\ { } &{} \quad {a_{j,2}} &{} \quad { } &{} \quad { } &{} \quad { } \\ { } &{} \quad {} &{} \quad {\ddots } &{} \quad &{} { } \\ { } &{} \quad { } &{} \quad &{} \quad a_{j,N-1} &{} \quad {}\\ { 0 } &{} \quad { } &{} \quad { } &{} \quad {} &{} \quad {a_{j,N}}\\ \end{pmatrix},\ B_{(j)} = \begin{pmatrix} {b_{j,1}} &{} \quad { } &{} \quad { } &{} \quad { } &{} \quad { 0 } \\ { } &{} \quad {b_{j,2}} &{} \quad { } &{} \quad { } &{} \quad { } \\ { } &{} \quad {} &{} \quad {\ddots } &{} \quad &{} \quad { } \\ { } &{} \quad { } &{} \quad &{} \quad b_{j,N-1} &{} \quad {}\\ { 0 } &{} \quad { } &{} \quad { } &{} \quad {} &{} \quad {b_{j,N}}\\ \end{pmatrix}, \end{aligned}$$

where the coefficients $a_{j,i}$ (and $b_{j,i}$ ) are the point value approximations of $a_j$ (and $b_j$) described in Appendix D.2. For example,

$$\begin{aligned} a_{3,i}= & {} h_{1,i}\displaystyle \frac{2H_{i-1}-h_{1,i-1} -2\left( 2H_{i}-h_{1,i} \right) + 2H_{i+1}-h_{1,i+1}}{\varDelta x^2}\\&+ \,\displaystyle \frac{2H_{i+1}-h_{1,i+1} - \left( 2H_{i-1}-h_{1,i-1} \right) }{2\varDelta x} + 4. \end{aligned}$$

1.3 D.3 Analysis of the Linear System for Small Water Heights

If we assume

$$\begin{aligned} h=\epsilon ,\ u_\alpha =w_\alpha =0, \end{aligned}$$

then the coefficients (52) and (53) reduce to

$$\begin{aligned} \left\{ \begin{array}{l} a_1 = -l_1^2\epsilon ^2,\\ \\ a_2 = 0,\\ \\ a_3 = 4(1+(\partial _x H)^2) + 2l_1\epsilon \partial _{xx} H,\\ \\ a_4 = -l_1^2\epsilon ^2,\\ \\ a_5 = -4l_1\epsilon \partial _x H,\\ \\ a_6 = -4(1+(\partial _x H)^2)-2l_1\epsilon \partial _{xx} H, \end{array} \right. \qquad \left\{ \begin{array}{l} b_1 = -l_2((\gamma _1-2)l_2+2)\epsilon ^2,\\ \\ b_2 = 4l_2\epsilon \partial _x H,\\ \\ b_3 = 4\gamma _1(1+(\partial _x H)^2) + 2l_2(\gamma _1+2)\epsilon \partial _{xx}H,\\ \\ b_4 = -l_2((\gamma _2-2)l_2+2)\epsilon ^2,\\ \\ b_5 = -4l_2\epsilon \partial _x H,\\ \\ b_6 = 4\gamma _2(1+(\partial _x H)^2) + 2l_2(\gamma _2-2)\epsilon \partial _{xx} H, \end{array} \right. \end{aligned}$$

and the Right Hand Side vectors reduce to

$$\begin{aligned} \varvec{\mathcal {RHS}} = \begin{pmatrix}\varvec{\mathcal {RHS}}_1 \\ \\ \varvec{\mathcal {RHS}}_2 \end{pmatrix}=\begin{pmatrix}\varvec{0} \\ \\ \varvec{0} \end{pmatrix}. \end{aligned}$$

In the following analysis we will assume:

$$\begin{aligned} \epsilon ^2\approx 0,\ \epsilon \partial _x H \approx 0,\ \epsilon \partial _{xx} H\approx 0, \end{aligned}$$

(54)

and for the sake of simplicity we assume that $\partial _x H = m$. Then the linear system becomes

The matrix $\varvec{A}$ is invertible

since we assume in Remark 4 that $\gamma _1+\gamma _2\ne 0$.

We note that (54) collects the particular case of a slowly varying bathymetry $\partial _x H \approx 0,$ and in particular the case under study in this work, when $\partial _x H = m$ with $\epsilon \cdot m\approx 0$.