Modeling of Wave Propagation from Arbitrary Depths to Shallow Waters – A Review

Abstract

Numerical models are a useful instrument for studying complex superposition of wave–wave and wave–current interactions in coastal and estuarine regions and to investigate the interaction of waves with complex bathymetries or structures built in nearshore areas. Moreover, since their applications are significantly less expensive and more flexible than the construction of physical models, they are a convenient tool to support design. The ability of the standard Boussinesq and Serre or Green and Naghdi equations to reproduce these nonlinear processes is well known. However, these models are restricted to shallow water conditions, and addition of other terms of dispersive origin has been considered since the 1990s, particularly for approximations of the Boussinesq-type. To allow applications in a greater range of h0 / λ, other than shallow waters, where h0 is the water depth at rest and λ is the wavelength, a new set of Serre-type equations, with additional terms of dispersive origin, is developed and tested with the available data and with a numerical solution of a Boussinesq-type model, also improved with dispersive characteristics. Explicit and implicit methods of finite–difference are implemented to solve both approximations of Boussinesq and Serre types with improved dispersive characteristics. A finite element method is also implemented to solve an extension of a Boussinesq-type model that takes into account wave–current interactions. Application examples to solve real-world problems are shown and discussed. The performances of both 1HD models are compared with experimental data of very demanding applications, namely: (i) a highly nonlinear solitary wave propagating up a slope and reflecting from a vertical wall and (ii) a periodic wave propagating in intermediate-depth waters upstream a trapezoidal bar, followed by very shallow waters over the bar, and again in intermediate-depth water conditions downstream.

1. Introduction

In recent decades, significant advances have been made in developing mathematical and numerical models to describe the entire phenomena observed in shallow water conditions. Indeed, not only our understanding of the phenomena has significantly improved but also the computational capabilities that are available have also increased considerably. In this context, we can now use more powerful and more reliable tools in the design of structures commonly used in coastal environments.

By the end of the 1970s, due to the lack of sufficiently deep knowledge, but above all for lack of computing power, the use of the linear wave theory for the simulation of phenomena, such as refraction and diffraction of waves, was common practice. In the 1980s, other models that take into account not only the refraction but also the diffraction process have been proposed and commonly used by [12] Berkhoff et al., 1982, [21] Kirby and Dalrymple, 1983, [13] Booij, 1983, [20] Kirby, 1984, and [15] Dalrymple, 1988, among many others. However, as they are based on the linear theory, those models should not be utilized in shallow water conditions.

As noted in [6], even at that time models based on the Saint-Venant equations (see also [32] Saint-Venant, 1871) “were frequently used in practical applications. However, as has been widely demonstrated, in shallow water conditions and for some types of waves, models based on a non-dispersive theory, of which the Saint-Venant model is an example, are limited and are not usually able to compute satisfactory results over long periods of analysis” (see also [33]). Nowadays, it is generally accepted that for practical applications the combined gravity wave effects in shallow water conditions must be taken into account. In addition, the refraction and diffraction processes, the swelling, reflection and breaking waves, all have to be considered.

Also according to [6] “a number of factors have made it possible to employ increasingly complex mathematical models”. Indeed, not only has there been great improvement in our theoretical knowledge of the phenomena involved but also the numerical methods have been used more efficiently. The great advances made in computer technology, especially since the 1980s, improving information processing and enabling large amounts of data to be stored, have made possible the use of more mathematical models, of greater complexity and with fewer restrictions. Indeed, only models of order $${\sigma }^2$$ ($$\sigma =h_0/\lambda$$, where $$h_0$$ and $$\lambda$$ represent, respectively, depth and wavelength characteristics) or greater, of the Boussinesq or Serre types [14, 35], are able to describe the entire dispersive and nonlinear interaction process of generation, propagation and run-up of waves resulting from wave–wave and wave–current interactions. It is also worth to point out that in more complex problems, such as wave generation by seafloor movement, propagation over uneven bottoms, and added breaking effects, high-frequency waves can arise as a consequence of nonlinear interaction.

In the past few years, the possibility of using more powerful computational facilities, and the technological evolution and sophistication of control systems have required thorough theoretical and experimental research designed to improve the knowledge of coastal hydrodynamics. Numerical methods aimed for the applications in engineering fields that are more sophisticated and with a higher degree of complexity have also been developed.

In Section 2, the general shallow water wave theory is used to develop different mathematical approaches, which are nowadays the basis of the most sophisticated models in hydrodynamics and sedimentary dynamics. Extensions of the equations for intermediate water of the more general approaches (order $${\sigma }^2$$) are presented in Section 3. Numerical formulations of the models and appropriate application examples are presented in Section 4. Firstly, a submerged dissipation platform to protect the Bugio lighthouse, situated at the Tagus estuary mouth, Lisbon, is designed and tested numerically. The second application is a real-world problem concerning coastal protection, using a submerged structure to force the breaking waves offshore. The third example shows the agitations established in the port of Figueira da Foz, Portugal, for different sea states, with the objective of designing new protective jetties, or extend the existing ones. Finally, the fourth example shows the performance of the extended Serre equations for the propagation of a wave in very demanding conditions. Conclusions and future developments are given in Section 5.

2. Mathematical formulations

We start from the fundamental equations of the Fluid Mechanics, written in Euler’s variables, relating to a three-dimensional and quasi-irrotational flow of a perfect fluid [Euler equations, or Navier-Stokes equations with the assumptions of non-compressibility ($$d\rho /dt=div\overrightarrow{v}=0$$), irrotationality ($$rot\overrightarrow{v}=0$$) and perfect fluid (dynamic viscosity, $$\mu =0$$)]:

$$\begin{array}{l}{u}_x+v_y+w_z=0\\ u_t+u{u}_x+v{u}_y+w{u}_z=-p_x/\rho \\ v_t+u{v}_x+v{v}_y+w{v}_z=-p_y/\rho \\ w_t+u{w}_x+v{w}_y+w{w}_z=-p_z/\rho -g\\ u_z=w_x;v_z=w_y;v_x=u_y\end{array}$$E1

with $$p=0$$ at $$z=\eta \left(x,y,t\right)$$, $$w={\eta }_t+u{\eta }_x+v{\eta }_y$$ at $$z=\eta \left(x,y,t\right)$$, and $$w={\xi }_t+u{\xi }_x+v{\xi }_y$$ at $$z=-h_0+\xi \left(x,y,t\right)$$. In these equations $$\rho$$ is density, t is time, g is gravitational acceleration, p is pressure, $$\eta$$ is free surface elevation, $$\xi$$ is bottom, and u, v, w are velocity components. Defining the dimensionless quantities $$\epsilon =a/h_0$$ and $$\sigma =h_0/\lambda$$, in which $$a$$ is a characteristic wave amplitude, $$h_0$$ represents water depth, and $$\lambda$$ is a characteristic wavelength, we proceed with suitable nondimensional variables:

$$\begin{array}{l}x\text{'}=x/\lambda ,y\text{'}=y/\lambda ,z\text{'}=z/h_0,\eta \text{'}=\eta /a,\xi \text{'}=\xi /h_0,t\text{'}=t\sqrt{g{h}_0}/\lambda =t{c}_0/\lambda ,p\text{'}=p/\left(\rho g{h}_0\right),\\ u\text{'}=u/\left(a\sqrt{g/h_0}\right)=u{h}_0/\left(a{c}_0\right),v\text{'}=v/\left(a\sqrt{g/h_0}\right)=v{h}_0/\left(a{c}_0\right),w\text{'}=w\lambda /\left(a{h}_0\sqrt{g/h_0}\right)=w\lambda /\left(a{c}_0\right),\end{array}$$BB1

where $$c_0$$ represents critical celerity, given by $$c_0={\left(g{h}_0\right)}^{1/2}$$, and, as above,$$\eta$$ is free surface elevation, $$\xi$$ represents bathymetry, $$u$$, $$v$$ and $$w$$ are velocity components, and $$p$$ is pressure. In dimensionless variables, without the line over the variables, the fundamental equations and the boundary conditions are written [5]:

1. Fundamental equations

$\begin{array}{ll}{u}_x+v_y+w_z=0\hfill & \text{a}\hfill \\ \epsilon u_t+{\epsilon }^{2}u{u}_x+{\epsilon }^{2}v{u}_y+{\epsilon }^{2}w{u}_z=-p_x\hfill & \text{b}\hfill \\ \epsilon v_t+{\epsilon }^{2}u{v}_x+{\epsilon }^{2}v{v}_y+{\epsilon }^{2}w{v}_z=-p_y\hfill & \text{c}\hfill \\ \epsilon {\sigma }^2{w}_t+{\epsilon }^2{\sigma }^{2}u{w}_x+{\epsilon }^2{\sigma }^{2}v{w}_y+{\epsilon }^2{\sigma }^{2}w{w}_z=-p_z-1\hfill & \text{d}\hfill \\ u_z={\sigma }^2{w}_x;v_z={\sigma }^2{w}_y;v_x=u_y\hfill & \text{e}\hfill \end{array}$E2

1. Boundary conditions

$$\begin{array}{ll}w=\left(1/\epsilon \right){\xi }_t+u{\xi }_x+v{\xi }_y,z=-1+\xi \hfill & \text{a}\hfill \\ w={\eta }_t+\epsilon u{\eta }_x+\epsilon v{\eta }_y,z=\epsilon \eta \hfill & \text{b}\hfill \\ p=0,z=\epsilon \eta \hfill & \text{c}\hfill \end{array}$$E3

Integrating the first equation 2(a) between the bed $$-1+\xi$$ and the free surface $$\epsilon \eta$$, taking into account 3(a) and 3(b), yields the continuity equation (4):

$${\left[\eta -\left(1/\epsilon \right)\xi \right]}_t+{\left[\left(1-\xi +\epsilon \eta \right)\overline{u}\right]}_x+{\left[\left(1-\xi +\epsilon \eta \right)\overline{v}\right]}_y=0$$E4

where the bar over the variables represents the average value along the vertical. Then, accepting the fundamental hypothesis of the shallow water theory, $$\sigma =h_0/\lambda <<1$$, and developing the dependent variables in power series of the small parameter $${\sigma }^2$$, that is

$$f={\displaystyle \sum _{i=0}^{\infty }{\left({\sigma }^2\right)}^i{f}_i},\text{for}f=\left(u,v,w,p,\eta ,\xi ,A\right)$$E5

where $$A=u_x+v_y$$, from continuity 2(a), and with 3(a) and 3(b), the following equations are obtained:

$$w_0=-\left(z+1-{\xi }_0\right)A_0+w_0^{*}$$E6

$$w_0^{**}=-\left(1+\epsilon {\eta }_0-{\xi }_0\right)A_0+w_0^{*}$$E7

where the simple and double asterisk represent the variables values at the bottom and at the surface, respectively. Of 2(e) we obtain, successively [34]:

$$\begin{array}{l}{u}_0=u_0\left(x,y,t\right)\\ v_0=v_0\left(x,y,t\right)\end{array}$$E8

$$\begin{array}{l}{u}_1=-\left(1/2\right){\left(z+1-{\xi }_0\right)}^2{A}_{0x}+\left(z+1-{\xi }_0\right)\left({\xi }_{0x}{A}_0+w_0^{*}{}_x\right)+u_1^{*}\\ v_1=-\left(1/2\right){\left(z+1-{\xi }_0\right)}^2{A}_{0y}+\left(z+1-{\xi }_0\right)\left({\xi }_{0y}{A}_0+w_0^{*}{}_y\right)+v_1^{*}\end{array}$$E9

so that the average values of the horizontal components of the velocity, on the vertical, are given by:

$$\begin{array}{l}\overline{u}=u_0+{\sigma }^2{u}_1^{*}-\left({\sigma }^2/6\right){\left(1+\epsilon {\eta }_0-{\xi }_0\right)}^2{A}_{0x}\\ +\left({\sigma }^2/2\right)\left(1+\epsilon {\eta }_0-{\xi }_0\right)\left({\xi }_{0x}{A}_0+w_{0x}^{*}\right)+O\left({\sigma }^4\right)\\ \overline{v}=v_0+{\sigma }^2{v}_1^{*}-\left({\sigma }^2/6\right){\left(1+\epsilon {\eta }_0-{\xi }_0\right)}^2{A}_{0y}\\ +\left({\sigma }^2/2\right)\left(1+\epsilon {\eta }_0-{\xi }_0\right)\left({\xi }_{0y}{A}_0+w_{0y}^{*}\right)+O\left({\sigma }^4\right)\end{array}$$E10

On the other hand, taking into account that,

$$f=f_0+O\left({\sigma }^2\right)\text{for}f=\left(u,v,\eta ,\xi ,w^{*}\right)$$E11

from (5) and (9) we obtain:

$$\begin{array}{l}{u}^{**}=\overline{u}-\left({\sigma }^2/3\right){\left(1+\epsilon \eta -\xi \right)}^2{\overline{A}}_x\\ +\left({\sigma }^2/2\right)\left(1+\epsilon \eta -\xi \right)\left({\xi }_x\overline{A}+w_x^{*}\right)+O\left({\sigma }^4\right)\\ v^{**}=\overline{v}-\left({\sigma }^2/3\right){\left(1+\epsilon \eta -\xi \right)}^2{\overline{A}}_y\\ +\left({\sigma }^2/2\right)\left(1+\epsilon \eta -\xi \right)\left({\xi }_y\overline{A}+w_y^{*}\right)+O\left({\sigma }^4\right)\end{array}$$E12

Representing by $$\mathrm{\Gamma }=w_t+\epsilon u{w}_x+\epsilon v{w}_y+\epsilon w{w}_z$$ the vertical acceleration of the particles, we get $$\mathrm{\Gamma }=w_{0t}+\epsilon u_0{w}_{0x}+\epsilon v_0{w}_{0y}+\epsilon w_0{w}_{0z}+O\left({\sigma }^2\right)$$, and from (6), (7) and (11) the following approach is obtained:

$$\begin{array}{l}\mathrm{\Gamma }=-\left(z+1-\xi \right)\left[{\overline{A}}_t+\epsilon \overline{u}{\overline{A}}_x+\epsilon \overline{v}{\overline{A}}_y-\epsilon {\overline{A}}^2\right]\\ +\left[w_t^{*}+\epsilon \overline{u}{w}_x^{*}+\epsilon \overline{v}{w}_y^{*}\right]+O\left({\sigma }^2\right)\end{array}$$E13

in which the terms within the first parenthesis (straight parenthesis) represent the vertical acceleration when the bottom is horizontal, and the terms within the second parenthesis (straight parenthesis) represent the vertical acceleration along the real bottom. It should be noted that equation 2(d) can be written as:

$$\epsilon {\sigma }^2\mathrm{\Gamma }=-p_z-1$$E14

where, by vertical integration between the bottom and the surface, the pressure p on the surface is obtained as:

$$\begin{array}{l}{p}_x^{**}=\epsilon {\eta }_x\left(\epsilon {\sigma }^2{\mathrm{\Gamma }}^{**}+1\right)\\ p_y^{**}=\epsilon {\eta }_y\left(\epsilon {\sigma }^2{\mathrm{\Gamma }}^{**}+1\right)\end{array}$$E15

which, along with 2(b) and 2(c), allow us to obtain [34]:

$$\begin{array}{l}{\left(u_t+\epsilon u{u}_x+\epsilon v{u}_y+\epsilon w{u}_z\right)}^{**}+{\eta }_x\left(1+\epsilon {\sigma }^2{\mathrm{\Gamma }}^{**}\right)=0\\ {\left(v_t+\epsilon u{v}_x+\epsilon v{v}_y+\epsilon w{v}_z\right)}^{**}+{\eta }_y\left(1+\epsilon {\sigma }^2{\mathrm{\Gamma }}^{**}\right)=0\end{array}$$E16

or even, given that $${\left(f_s\right)}^{**}=f_s^{**}-\epsilon {\left(f_z\right)}^{**}{\eta }_s$$, where $$f=\left(u,v\right)$$ and $$s=\left(x,y,t\right)$$:

$$\begin{array}{l}{u}_t^{**}+\epsilon u^{**}{u}_x^{**}+\epsilon v^{**}{u}_y^{**}+{\eta }_x\left(1+\epsilon {\sigma }^2{\mathrm{\Gamma }}^{**}\right)=0\\ v_t^{**}+\epsilon u^{**}{v}_x^{**}+\epsilon v^{**}{v}_y^{**}+{\eta }_y\left(1+\epsilon {\sigma }^2{\mathrm{\Gamma }}^{**}\right)=0\end{array}$$E17

By developing expressions (17) in second approach (order 2 in $${\sigma }^2$$), the following equations of motion (18) are obtained (for details see [34]):

$$\begin{array}{l}{\overline{u}}_t+\epsilon \overline{u}{\overline{u}}_x+\epsilon \overline{v}{\overline{u}}_y+{\eta }_x\\ +{\sigma }^2\left\{\left[\left(2/3\right){\left(\epsilon \eta -\xi \right)}_x+\left(1/2\right){\xi }_x\right]P+\left(1/3\right)\left(1+\epsilon \eta -\xi \right)P_x\right\}\\ +{\sigma }^2\left[\epsilon {\eta }_{x}Q+\left(1/2\right)\left(1+\epsilon \eta -\xi \right)Q_x\right]+{\sigma }^4=0\\ {\overline{v}}_t+\epsilon \overline{u}{\overline{v}}_x+\epsilon \overline{v}{\overline{v}}_y+{\eta }_y\\ +{\sigma }^2\left\{\left[\left(2/3\right){\left(\epsilon \eta -\xi \right)}_y+\left(1/2\right){\xi }_y\right]P+\left(1/3\right)\left(1+\epsilon \eta -\xi \right)P_y\right\}\\ +{\sigma }^2\left[\epsilon {\eta }_{y}Q+\left(1/2\right)\left(1+\epsilon \eta -\xi \right)Q_y\right]+{\sigma }^4=0\\ P=\left(1+\epsilon \eta -\xi \right)\left(\epsilon {\overline{A}}^2-\epsilon \overline{u}{\overline{A}}_x-\epsilon \overline{v}{\overline{A}}_y-{\overline{A}}_t\right)\\ Q=w_t+\epsilon \overline{u}{w}_x+\epsilon \overline{v}{w}_y\\ w=\left(1/\epsilon \right){\xi }_t+\overline{u}{\xi }_x+\overline{v}{\xi }_y\\ \overline{A}={\overline{u}}_x+{\overline{v}}_y\end{array}$$E18

where, likewise, the bar over the variables represents the average value along the vertical. In dimensional variables and with a solid/fixed bottom ($${\xi }_t=0$$), the complete set of equations is written, in second approach:

$$\begin{array}{l}{h}_t+{\left(h\overline{u}\right)}_x+{\left(h\overline{v}\right)}_y=0\\ {\overline{u}}_t+\overline{u}{\overline{u}}_x+\overline{v}{\overline{u}}_y+g{\eta }_x\\ +\left[\left(2/3\right)h_x+\left(1/2\right){\xi }_x\right]P+\left(1/3\right)h{P}_x+h_{x}Q+\left(1/2\right)h{Q}_x=0\\ {\overline{v}}_t+\overline{u}{\overline{v}}_x+\overline{v}{\overline{v}}_y+g{\eta }_y\\ +\left[\left(2/3\right)h_y+\left(1/2\right){\xi }_y\right]P+\left(1/3\right)h{P}_y+h_{y}Q+\left(1/2\right)h{Q}_y=0\\ P=h\left({\overline{A}}^2-\overline{u}{\overline{A}}_x-\overline{v}{\overline{A}}_y-{\overline{A}}_t\right)\\ Q=w_t+\overline{u}{w}_x+\overline{v}{w}_y\\ w=\overline{u}{\xi }_x+\overline{v}{\xi }_y\\ \overline{A}={\overline{u}}_x+{\overline{v}}_y\end{array}$$E19

where $$h=h_0-\xi +\eta$$ is total water depth. The one-dimensional form (1HD) of the equation system (19) is written, also with a fixed bottom:

$$\begin{array}{l}{h}_t+{\left(\overline{u}h\right)}_x=0\\ h{\overline{u}}_t+h\overline{u}{\overline{u}}_x+gh{\eta }_x+{\left[h^2\left(P/3+Q/2\right)\right]}_x+{\xi }_{x}h\left(P/2+Q\right)=0\\ P=-h\left({\overline{u}}_{xt}+\overline{u}{\overline{u}}_{xx}-{\overline{u}}_x^2\right)\\ Q={\xi }_x\left({\overline{u}}_t+\overline{u}{\overline{u}}_x\right)+{\xi }_{xx}{\overline{u}}^2\end{array}$$E20

Assuming additionally a relative elevation of the surface due to the waves ($$\epsilon =a/h_0$$) having a value close to the square of the relative depth ($$\sigma =h_0/\lambda$$), i.e. $$O\left(\epsilon \right)$$ = $$O\left({\sigma }^2\right)$$, from the system of equations (18), and at the same order of approximation, the following approach is obtained, in dimensional variables:

$$\begin{array}{l}{h}_t+{\left(h\overline{u}\right)}_x+{\left(h\overline{v}\right)}_y=0\\ {\overline{u}}_t+\overline{u}{\overline{u}}_x+\overline{v}{\overline{u}}_y+g{\eta }_x-\left[\left(1/6\right){\xi }_x\right]P+\left(1/3\right)h_r{P}_x+\left(1/2\right)h_r{Q}_x=0\\ {\overline{v}}_t+\overline{u}{\overline{v}}_x+\overline{v}{\overline{v}}_y+g{\eta }_y-\left[\left(1/6\right){\xi }_y\right]P+\left(1/3\right)h_r{P}_y+\left(1/2\right)h_r{Q}_y=0\end{array}$$E21

where $$h_r=h_0-\xi$$ is the water column height at rest, $$P$$ and $$Q$$ are given by $$P=-\left(h_0-\xi \right){\left({\overline{u}}_x+{\overline{v}}_y\right)}_t$$ and $$Q={\left(\overline{u}{\xi }_x+\overline{v}{\xi }_y\right)}_t$$. The momentum equations are written as:

$$\begin{array}{l}{\overline{u}}_t+\overline{u}{\overline{u}}_x+\overline{v}{\overline{u}}_y+g{\eta }_x+\left(1/6\right)h_r{\xi }_x{\left({\overline{u}}_x+{\overline{v}}_y\right)}_t-\left(1/3\right)h_r^2{\left({\overline{u}}_x+{\overline{v}}_y\right)}_{xt}\\ +\left(1/3\right)h_r{\xi }_x{\left({\overline{u}}_x+{\overline{v}}_y\right)}_t+\left(1/2\right)h_r{\left(\overline{u}{\xi }_x+\overline{v}{\xi }_y\right)}_{xt}=0\end{array}$$E22

$$\begin{array}{l}{\overline{v}}_t+\overline{u}{\overline{v}}_x+\overline{v}{\overline{v}}_y+g{\eta }_y+\left(1/6\right)h_r{\xi }_y{\left({\overline{u}}_x+{\overline{v}}_y\right)}_t-\left(1/3\right)h_r^2{\left({\overline{u}}_x+{\overline{v}}_y\right)}_{yt}\\ +\left(1/3\right)h_r{\xi }_y{\left({\overline{u}}_x+{\overline{v}}_y\right)}_t+\left(1/2\right)h_r{\left(\overline{u}{\xi }_x+\overline{v}{\xi }_y\right)}_{yt}=0\end{array}$$E23

with $${\xi }_t=0$$, the complete system of equations (24) is obtained:

$$\begin{array}{l}{h}_t+{\left(h\overline{u}\right)}_x+{\left(h\overline{v}\right)}_y=0\\ {\overline{u}}_t+\overline{u}{\overline{u}}_x+\overline{v}{\overline{u}}_y+g{\eta }_x-\left(1/3\right)h_r^2\left({\overline{u}}_{xxt}+{\overline{v}}_{xyt}\right)+h_r{\xi }_x{\overline{u}}_{xt}\\ +\left(1/2\right)h_r\left({\xi }_{xx}{\overline{u}}_t+{\xi }_y{\overline{v}}_{xt}+{\xi }_x{\overline{v}}_{yt}+{\xi }_{xy}{\overline{v}}_t\right)=0\\ {\overline{v}}_t+\overline{u}{\overline{v}}_x+\overline{v}{\overline{v}}_y+g{\eta }_y-\left(1/3\right)h_r^2\left({\overline{u}}_{xyt}+{\overline{v}}_{yyt}\right)+h_r{\xi }_y{\overline{v}}_{yt}\\ +\left(1/2\right)h_r\left({\xi }_{yy}{\overline{v}}_t+{\xi }_y{\overline{u}}_{xt}+{\xi }_x{\overline{u}}_{yt}+{\xi }_{xy}{\overline{u}}_t\right)=0\end{array}$$E24

Further simplifying the equations of motion (18), retaining only terms up to order 1 in $$\sigma$$, i.e., neglecting all terms of dispersive origin, this system of equations is written in dimensional variables:

$$\begin{array}{l}{h}_t+{\left(h\overline{u}\right)}_x+{\left(h\overline{v}\right)}_y=0\\ {\overline{u}}_t+\overline{u}{\overline{u}}_x+\overline{v}{\overline{u}}_y+g{\eta }_x=0\\ {\overline{v}}_t+\overline{u}{\overline{v}}_x+\overline{v}{\overline{v}}_y+g{\eta }_y=0\end{array}$$E25

Approaches (19), (24) and (25) are known as Serre equations, or Green & Naghdi, Boussinesq and Saint-Venant, respectively, in two horizontal dimensions (2HD models). The classical Serre equations (19) [17] are fully nonlinear and weakly dispersive. Boussinesq equations (24) only incorporate weak dispersion and weak nonlinearity and are valid only for long waves in shallow waters. As for the Boussinesq-type models, also Serre’s equations are valid only for shallow water conditions.

3. Derivation of higher-order equations

3.1. Weakly nonlinear approaches with improved dispersive performance

3.1.1. Nwogu’s approach

To allow applications in a greater range of $$h_0/\lambda$$, other than shallow waters, [27] introduced higher-order dispersive terms into the governing equations to improve linear dispersion properties. By redefining the dependent variable, [30] achieved the same improvement without the need to add such terms to the equations.

Following [36], the extended Bussinesq equations obtained by [30] are derived in this section, using the nondimensionalised scaled equation system (2) – (3) as the starting point, rather than the procedure presented by Nwogu. For simplicity reasons, and without loss of generality, the one-dimensional extended Boussinesq equation system is derived. For consistency with the previous work (Section 2), the continuity equation 2(a) in two vertical dimensions (x, z) is integrated through the depth:

$${\displaystyle \underset{-1+\xi }{\overset{z}{\int }}{w}_{z}d}z=-{\displaystyle \underset{-1+\xi }{\overset{z}{\int }}{u}_{x}dz}$$E26

$${w|}_z-{w|}_{-1+\xi }=-{\left({\displaystyle \underset{-1+\xi }{\overset{z}{\int }}udz}\right)}_x+{u|}_{-1+\xi }{\xi }_x-{u|}_z{z}_x$$E27

Denoting $${w|}_z$$ by $$w$$, with $${\xi }_t=0$$ and using Leibnitz’ rule, the boundary condition 3(a) at $$z=-1+\xi$$ gives:

$$w=-{\left({\displaystyle \underset{-1+\xi }{\overset{z}{\int }}udz}\right)}_x$$E28

Substituting (28) in the irrotationality equation 2(e), $$u_z={\sigma }^2{w}_x$$, yields:

$$u_z={\sigma }^2{w}_x=-{\sigma }^2{\left({\displaystyle \underset{-1+\xi }{\overset{z}{\int }}udz}\right)}_{xx}$$E29

Considering a Taylor series expansion (30) of $$u\left(x,z,t\right)$$, about $$z=z_{\alpha }$$, where the horizontal velocity $$u_{\alpha }\left(x,t\right)$$ denotes the velocity at depth $$z_{\alpha }$$, $$u\left(x,z_{\alpha },t\right)$$, this Taylor expansion is integrated through the depth from $$-1+\xi$$ to $$z$$, yielding (31) (for details see [36]):

$$u=u_{\alpha }+\left(z-z_{\alpha }\right){u_z|}_{z=z_{\alpha }}+{\frac{{\left(z-z_{\alpha }\right)}^2}{2}{u}_{zz}|}_{z=z_{\alpha }}+\cdots$$E30

$$\begin{array}{l}{\displaystyle \underset{-1+\xi }{\overset{z}{\int }}udz}=\left(z+1-\xi \right)u_{\alpha }+{\left[\frac{{\left(z-z_{\alpha }\right)}^2}{2}-\frac{{\left(1-\xi +z_{\alpha }\right)}^2}{2}\right]u_z|}_{z=z_{\alpha }}\\ +{\left[\frac{{\left(z-z_{\alpha }\right)}^3}{6}+\frac{{\left(1-\xi +z_{\alpha }\right)}^3}{6}\right]u_{zz}|}_{z=z_{\alpha }}+\cdots \end{array}$$E31

Substituting (31) in equation (29) gives:

$$u_z=-{\sigma }^2{\left\{\left(z+1-\xi \right)u_{\alpha }+{\left[\frac{{\left(z-z_{\alpha }\right)}^2}{2}-\frac{{\left(1-\xi +z_{\alpha }\right)}^2}{2}\right]u_z|}_{z=z_{\alpha }}+\cdots \right\}}_{xx}$$E32

Differentiating equation (32) with respect to z, noting from (29) that $$u_z=O\left({\sigma }^2\right)$$:

$$\begin{array}{l}{u}_{zz}=-{\sigma }^2{\left\{u_{\alpha }+\left(z-z_{\alpha }\right){u_z|}_{z=z_{\alpha }}+\frac{{\left(z-z_{\alpha }\right)}^2}{2}{u_{zz}|}_{z=z_{\alpha }}+\cdots \right\}}_{xx}\\ =-{\sigma }^2{u}_{\alpha }{}_{xx}+O\left({\sigma }^4\right)\end{array}$$E33

Differentiating equation (33) with respect to z and noting that both $$u_z$$ and $$u_{zz}$$ are $$O\left({\sigma }^2\right)$$:

$$u_{zzz}=-{\sigma }^2{\left\{{u_z|}_{z=z_{\alpha }}+\left(z-z_{\alpha }\right){u_{zz}|}_{z=z_{\alpha }}+\cdots \right\}}_{xx}=O\left({\sigma }^4\right)$$E34

Repeated differentiation of this expression will produce expressions for the higher derivatives of u with respect to z and show them to be of $$O\left({\sigma }^4\right)$$ order or greater. Substituting equations (33) and (34) back in equation (32):

$$\begin{array}{l}{u}_z=-{\sigma }^2{\left\{\left(1-\xi +z_{\alpha }\right)u_{\alpha }-{\frac{{\left(1-\xi +z_{\alpha }\right)}^2}{2}{u}_z|}_{z=z_{\alpha }}+O\left({\sigma }^2\right)\right\}}_{xx}\\ =-{\sigma }^2{\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}+O\left({\sigma }^4\right)\end{array}$$E35

Substituting equations (33), (34) and (35) in the Taylor series expansion (30) produces an expression for the horizontal velocity component u:

$$u=u_{\alpha }-{\sigma }^2\left\{\left(z-z_{\alpha }\right){\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}+\frac{{\left(z-z_{\alpha }\right)}^2}{2}{\left(u_{\alpha }\right)}_{xx}\right\}+O\left({\sigma }^4\right)$$E36

Substituting the horizontal velocity (36) in equation (28) leads to an expression for the vertical velocity component w:

$$\begin{array}{l}w=-\{\left(1-\xi +z_{\alpha }\right)u_{\alpha }\\ {+{\sigma }^2\left[\frac{{\left(1-\xi +z_{\alpha }\right)}^2}{2}{\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}-\frac{{\left(1-\xi +z_{\alpha }\right)}^3}{6}{\left(u_{\alpha }\right)}_{xx}\right]\}}_x\\ +O\left({\sigma }^{4\text{}}\right)\end{array}$$E37

Using the velocities u (36) and w (37) in the vertical momentum, equation 2(d) yields:

$$\epsilon {\sigma }^2\left\{-{\left[\left(z+1-\xi \right)u_{\alpha }\right]}_{xt}+O\left({\sigma }^2\right)\right\}+p_z+1+O\left({\epsilon }^2{\sigma }^2\right)=0$$E38

This can be rearranged to give an expression for the pressure p:

$$-p_z=1-\epsilon {\sigma }^2{\left[\left(z+1-\xi \right)u_{\alpha }\right]}_{xt}+O\left(\epsilon {\sigma }^4,{\epsilon }^2{\sigma }^2\right)$$E39

Integrating through the depth from z to the free surface $$\epsilon \eta$$:

$$-{p|}_{\epsilon \eta }+{p|}_z=\epsilon \eta -z+\epsilon {\sigma }^2{\left\{\left[\frac{z^2}{2}+\left(1-\xi \right)z\right]u\right\}}_{xt}+O\left(\epsilon {\sigma }^4,{\epsilon }^2{\sigma }^2\right)$$E40

Using the free surface boundary condition 3(c) for the pressure, and denoting $${p|}_z$$ simply as p gives:

$$p=\epsilon \eta -z+\epsilon {\sigma }^2{\left\{\left[\frac{z^2}{2}+\left(1-\xi \right)z\right]u\right\}}_{xt}+O\left(\epsilon {\sigma }^4,{\epsilon }^2{\sigma }^2\right)$$E41

Substituting the equations for the horizontal velocity component (36), the vertical velocity component (37) and the pressure (41) in the horizontal momentum equation 2(b) gives the Boussinesq momentum equation:

$$\begin{array}{l}{u}_{\alpha }{}_t-{\sigma }^2{\left\{\left(z-z_{\alpha }\right){\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}+\frac{{\left(z-z_{\alpha }\right)}^2}{2}{\left(u_{\alpha }\right)}_{xx}\right\}}_t\\ +\epsilon u_{\alpha }{u}_{\alpha }{}_x+{\eta }_x+{\sigma }^2{\left\{\left[\frac{z^2}{2}+\left(1-\xi \right)z\right]u_{\alpha }\right\}}_{xxt}=O\left(\epsilon {\sigma }^2,{\sigma }^4\right)\end{array}$$E42

or

$$u_{\alpha }{}_t+\epsilon u_{\alpha }{u}_{\alpha }{}_x+{\eta }_x+{\sigma }^2{\left\{z_{\alpha }{\left[\left(1-\xi \right)u_{\alpha }\right]}_{xx}+\frac{z_{\alpha }^2}{2}{\left(u_{\alpha }\right)}_{xx}\right\}}_t=O\left(\epsilon {\sigma }^2,{\sigma }^4\right)$$E43

The second equation of the Boussinesq system is developed by first integrating the continuity equation 2(a) through the depth $$-1+\xi$$:

$${\left({\displaystyle \underset{-1+\xi }{\overset{\epsilon \eta }{\int }}udz}\right)}_x-\epsilon {u|}_{\epsilon \eta }{\eta }_x+{u|}_{-1+\xi }{\xi }_x+{w|}_{\epsilon \eta }-{w|}_{-1+\xi }=0$$E44

Using Liebnitz’ Rule and the kinematic boundary conditions at the bed $$z=-1+\xi$$ 3(a) and at the free surface $$z=\epsilon \eta$$ 3(b) gives:

$${\eta }_t+{\left({\displaystyle \underset{-1+\xi }{\overset{\epsilon \eta }{\int }}udz}\right)}_x=0$$E45

From the expression (36) for the horizontal velocity u:

$$\begin{array}{l}{\displaystyle \underset{-1+\xi }{\overset{\epsilon \eta }{\int }}udz}=\left(1-\xi +\epsilon \eta \right)u_{\alpha }-{\sigma }^2\{{\left[\frac{{\left(z-z_{\alpha }\right)}^2}{2}\right]}_{-1+\xi }^{\epsilon \eta }{\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}\\ +{\left[\frac{{\left(z-z_{\alpha }\right)}^3}{6}\right]}_{-1+\xi }^{\epsilon \eta }{u}_{\alpha }{}_{xx}\}+O\left({\sigma }^4\right)\end{array}$$E46

or

$$\begin{array}{l}{\displaystyle \underset{-1+\xi }{\overset{\epsilon \eta }{\int }}udz}=\left(1-\xi +\epsilon \eta \right)u_{\alpha }-{\sigma }^2\{\left[\frac{z_{\alpha }^2}{2}-\frac{{\left(1-\xi +z_{\alpha }\right)}^2}{2}\right]{\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}\\ +\left[-\frac{z_{\alpha }^3}{6}+\frac{{\left(1-\xi +z_{\alpha }\right)}^3}{6}\right]u_{\alpha }{}_{xx}\}+O\left(\epsilon {\sigma }^2,{\sigma }^4\right)\end{array}$$E47

yielding:

$$\begin{array}{l}{\displaystyle \underset{-1+\xi }{\overset{\epsilon \eta }{\int }}udz}=\left(1-\xi +\epsilon \eta \right)u_{\alpha }-{\sigma }^2\{\left[-\frac{{\left(1-\xi \right)}^2}{2}+\left(1-\xi \right)z_{\alpha }\right]{\left[\left(1-\xi +z_{\alpha }\right)u_{\alpha }\right]}_{xx}\\ +\left[\frac{{\left(1-\xi \right)}^3}{6}+\frac{{\left(1-\xi \right)}^2{z}_{\alpha }}{2}+\frac{\left(1-\xi \right)z_{\alpha }^2}{2}\right]{\left(u_{\alpha }\right)}_{xx}\}+O\left(\epsilon {\sigma }^2,{\sigma }^4\right)\end{array}$$E48

and thus

$$\begin{array}{l}{\displaystyle \underset{-1+\xi }{\overset{\epsilon \eta }{\int }}udz}=\left(1-\xi +\epsilon \eta \right)u_{\alpha }+{\sigma }^2\{\left[z_{\alpha }+\frac{1-\xi }{2}\right]\left(1-\xi \right){\left[\left(1-\xi \right)u_{\alpha }\right]}_{xx}\\ +\left[\frac{z_{\alpha }^2}{2}-\frac{{\left(1-\xi \right)}^2}{6}\right]\left(1-\xi \right){\left(u_{\alpha }\right)}_{xx}\}+O\left(\epsilon {\sigma }^2,{\sigma }^4\right)\end{array}$$E49

Substituting equation (49) in the free surface equation (45) gives the Boussinesq continuity equation:

$$\begin{array}{l}{\eta }_t+{\left[\left(1-\xi +\epsilon \eta \right)u_{\alpha }\right]}_x+{\sigma }^2\{\left(z_{\alpha }+\frac{1-\xi }{2}\right)\left(1-\xi \right){\left[\left(1-\xi \right)u_{\alpha }\right]}_{xx}\\ +{\left(\frac{z_{\alpha }^2}{2}-\frac{{\left(1-\xi \right)}^2}{6}\right)\left(1-\xi \right){\left(u_{\alpha }\right)}_{xx}\}}_x=O\left(\epsilon {\sigma }^2,{\sigma }^4\right)\end{array}$$E50

Returning to dimensional variables, unscaled form, equations (43) and (50) are written:

$$u_{\alpha }{}_t+u_{\alpha }{u}_{\alpha }{}_x+g{\eta }_x+z_{\alpha }{\left(h_r{u}_{\alpha }{}_t\right)}_{xx}+\frac{z_{\alpha }^2}{2}{\left(u_{\alpha }\right)}_{xxt}=0$$E51

$${\eta }_t+{\left[\left(h_0+\eta \right)u_{\alpha }\right]}_x+{\left[\left(z_{\alpha }+\frac{h_r}{2}\right)h_r{\left(h_r{u}_{\alpha }\right)}_{xx}+\left(\frac{z_{\alpha }^2}{2}-\frac{h_r^2}{6}\right)h_r{\left(u_{\alpha }\right)}_{xx}\right]}_x=0$$E52

Setting the arbitrary depth $$z_{\alpha }=\alpha h_0$$, where $$-1\le \alpha \le 0$$, the system (51)–(52) is rewritten:

$$u_{\alpha }{}_t+u_{\alpha }{u}_{\alpha }{}_x+g{\eta }_x+\alpha h_r{\left[h_r{\left(u_{\alpha }\right)}_t\right]}_{xx}+\frac{{\alpha }^2}{2}{h}_r^2{\left(u_{\alpha }\right)}_{xxt}=0$$E53

$${\eta }_t+{\left[\left(h_r+\eta \right)u_{\alpha }\right]}_x+{\left[\left(\alpha +\frac{1}{2}\right)h_r^2{\left(h_r{u}_{\alpha }\right)}_{xx}+\left(\frac{{\alpha }^2}{2}-\frac{1}{6}\right)h_r^3{\left(u_{\alpha }\right)}_{xx}\right]}_x=0$$E54

In two dimensions, the equation system (53)–(54) is written:

$$\begin{array}{l}{\eta }_t+\nabla .\left[\left(h_r+\eta \right)u_{\alpha }\right]\\ +\nabla \left\{\left(\alpha +\frac{1}{2}\right)h_r^2\nabla \left[\nabla .\left(h_r{u}_{\alpha }\right)\right]+\left(\frac{{\alpha }^2}{2}-\frac{1}{6}\right)h_r^3\nabla \left(\nabla .u_{\alpha }\right)\right\}=0\end{array}$$E55

$$u_{\alpha }{}_t+\left(u_{\alpha }.\nabla \right)u_{\alpha }+g\nabla \eta +\alpha h_r\nabla \left[\nabla .\left(h_r{u}_{\alpha }{}_t\right)\right]+\frac{{\alpha }^2}{2}{h}_r^2\nabla \left(\nabla .u_{\alpha }{}_t\right)=0$$E56

where $$\nabla$$ represents the two-dimensional gradient operator with respect to the horizontal coordinates (x, y) ($$\nabla =\partial /\partial x,\partial /\partial y$$), and the velocity vector $$u\left(x,y,t\right)=\left(u,v\right)$$ represents the two-dimensional velocity field at depth $$z=\alpha h$$.

3.1.2. Beji and Nadaoka’s, and Liu and Sun’s approaches

Starting from the standard Boussinesq equations and adopting the methodology introduced by [27] [11] presented a new approach that allows for applications until values of $$h_0/\lambda$$ to the order of 0.25, and still with acceptable errors in amplitude and phase velocity up to values of $$h_0/\lambda$$ near 0.50. By an addition and subtraction process, using the approximation $$u_t=-g{\eta }_x$$ and considering a dispersion parameter $$\beta$$ in the momentum equation of the 1HD system (24), Beji and Nadaoka obtained an improved set of Boussinesq equations for variable depth, with $$\beta$$ value obtained by comparison of the dispersion relation of the linearized form of the resulting equation with a second-order Padé expansion of the linear dispersion relation $${\omega }^2/gk=\tanh \left(k{h}_r\right)$$. In order to improve dispersion and linear shoaling characteristics in the Beji and Nadaoka’s equations, [23] introduced two tuning parameters, $$\alpha$$ and $$\gamma$$, so that $$\beta =1.5\alpha -0.5\gamma$$. The nonlinearity in the previous Boussinesq-type models was improved by Liu and Sun adding higher-order terms accurate to the order of $$O\left(\epsilon {\sigma }^2\right)$$. The 1HD standard Boussinesq equations and the approaches of Beji and Nadaoka and Liu and Sun are identified within the following system of equations (57) for water of variable depth:

$$\begin{array}{l}{h}_t+{\left(hu\right)}_x=0\\ u_t+u{u}_x+g{\eta }_x-\left\{\left(1+\alpha \right)/2-\left(1+\gamma \right)/6\right\}{h}_r^2{u}_{xxt}\\ -\left(\alpha /2-\gamma /6\right)g{h}_r^2{\eta }_{xxx}+\left(1+\alpha \right)h_r{\xi }_x{u}_{xt}+\alpha g{h}_r{\xi }_x{\eta }_{xx}\\ +0.5\left(1+\alpha \right)h_r{\xi }_{xx}{u}_t+0.5\alpha g{h}_r{\xi }_{xx}{\eta }_x=O\left(\epsilon {\sigma }^2\right)\end{array}$$E57

After linearization of the equations system (57), [23] obtained the following dispersion relation (58):

$$\frac{{\omega }^2}{gk}=\frac{k{h}_r\left[1+\left(\alpha /2-\gamma /6\right)k^2{h}_r^2\right]}{1+\left(1+\alpha \right)k^2{h}_r^2/2-\left(1+\gamma \right)k^2{h}_r^2/6}$$E58

Comparing equation (58), written in terms of the phase speed (59)

$$C^2=\frac{{\omega }^2}{k^2}=g{h}_r\left[\frac{1+\left(\alpha /2-\gamma /6\right){\left(k{h}_r\right)}^2}{1+\left(1+\alpha \right){\left(k{h}_r\right)}^2/2-\left(1+\gamma \right){\left(k{h}_r\right)}^2/6}\right]$$E59

with the linear dispersion relation $${\omega }^2/gk=\tanh \left(k{h}_r\right)$$, using the approach (60)

$$C_{Airy}^2=\frac{{\omega }^2}{k^2}=\left(g{h}_r\right)\tanh \left(k{h}_r\right)=g{h}_r\left[\frac{1+{\left(k{h}_r\right)}^2/15}{1+2{\left(k{h}_r\right)}^2/5}\right]+O\left[{\left(kh\right)}^6\right]$$E60

allows to obtain the best values for the parameters $$\alpha$$ and $$\gamma$$: $$\alpha =0.1308$$ and $$\gamma =-0.0076$$.

Considering appropriate values for the tuning parameters $$\alpha$$ and $$\gamma$$, we can identify within the equation system (57):

• The standard Boussinesq equations by setting $$\alpha =\gamma =0$$

• The Beji and Nadaoka equations considering $$\alpha =\gamma =0.20$$

• The Liu and Sun equations with $$\alpha =0.1308$$ and $$\gamma =-0.0076$$

A visual comparison of numerical results of the extended Boussinesq approximation (57), with $$\alpha =0.1308$$ and $$\gamma =-0.0076$$), shown in [1, 2] and [3], with a similar study performed by [37], using the extended Boussinesq (53)–(54) model (Nwogu’s approach, with $$\alpha =-0.531$$) shows no relevant differences in the graphs. In this regard, it is worth remembering [4] “although the methods of derivation are different, the resulting dispersion relations of these extended Boussinesq equations are similar, and may be thought of as a slight modification of the second-order Padé approximant of the full dispersion relation”.

3.2. Fully nonlinear approaches with improved dispersive performance

Instead of using the horizontal velocity at a certain depth, other extensions of the Boussinesq equations have been made by using the velocity potential on an arbitrary depth, also with a tuning parameter. Wei et al. (1995) [39] used the Nwogu’s approach to derive a Boussinesq-type model which retains full nonlinearity by including $$O\left(\epsilon {\sigma }^2\right)$$ terms not considered in the Nwogu’s (53)–(54) system, and thus improving the nonlinearity to $$O\left(\epsilon \right)=1$$. Wei et al., and later [16] Gobbi et al., derived a fourth-order Boussinesq model in which the velocity potential is approximated by a fourth-order polynomial in z. In terms of non-dimensional variables, both consider the boundary value problem for potential flow, given by:

$$\begin{array}{l}{\phi }_{zz}+{\mu }^2{\nabla }^2\phi =0,-1+\xi \le z\le \epsilon \eta \\ {\phi }_z+{\mu }^2\nabla h_{1-\xi }.\nabla \phi =0,z=-1+\xi \\ \eta +{\phi }_t+\frac{1}{2}\epsilon \left[{\left(\nabla \phi \right)}^2+\frac{1}{{\mu }^2}{\left({\phi }_z\right)}^2\right]=0,z=\epsilon \eta \\ {\eta }_t+\epsilon \nabla \phi .\nabla \eta -\frac{1}{{\mu }^2}{\phi }_z=0,z=\epsilon \eta \end{array}$$E61

where, as above, z is the vertical coordinate starting at the still water level $$h_0\left(x,y\right)$$ and pointing upwards, scaled by a typical depth $$h_0$$, and $$\eta$$ is the water surface displacement scaled by a representative amplitude a. The two dimensionless parameters $$\epsilon$$ and $${\mu }^2$$ are defined as $$\epsilon =a/h_0$$ and $${\mu }^2={\left(k_0{h}_0\right)}^2$$, with a representative wave number $$k_0=2\pi /\lambda$$, so $${\mu }^2={\left(2\pi \sigma \right)}^2$$. Time t is scaled by $${\left[k_0{\left(g{h}_0\right)}^{1/2}\right]}^{-1}$$, and $$\phi$$, the velocity potential, is scaled by $$\epsilon h_0{\left(g{h}_0\right)}^{1/2}$$. We use the nondimensional water level $$h_{1-\xi }$$ instead of $$1-\xi$$. Integrating the first equation of (61) over the water column, and using the appropriate boundary conditions, the continuity equation is obtained:

$${\eta }_t+\nabla .M=0$$E62

where $$M={\displaystyle {\int }_{-1+\xi }^{\epsilon \eta }\nabla \phi dz}$$. Retaining terms to $$O\left({\mu }^2\right)$$, and denoting $${\phi }_{\alpha }$$ as the value of $$\phi$$ at $$z=z_{\alpha }\left(x,y\right)$$, an approximate velocity potential is given by:

$$\phi ={\phi }_{\alpha }\left(x,t\right)+{\mu }^2\left(z_{\alpha }-z\right)\nabla .\left(h_{1-\xi }\nabla {\phi }_{\alpha }\right)+\frac{{\mu }^2}{2}\left(z_{\alpha }^2-z^2\right){\nabla }^2{\phi }_{\alpha }+O\left({\mu }^4\right)$$E63

Substituting equation (63) into (62), a mass flux conservation equation is obtained [39]:

$$\begin{array}{l}{\eta }_t+\nabla .\{\left(h_{1-\xi }+\epsilon \eta \right)\{\nabla {\phi }_{\alpha }+{\mu }^2\{\nabla \left[z_{\alpha }\nabla .\left(h_{1-\xi }\nabla {\phi }_{\alpha }\right)+\frac{z_{\alpha }^2}{2}{\nabla }^2{\phi }_{\alpha }\right]\\ +\frac{\left(h_{1-\xi }-\epsilon \eta \right)}{2}\nabla \left[\nabla .\left(h_{1-\xi }\nabla {\phi }_{\alpha }\right)\right]-\frac{\left[h_{1-\xi }^2-\epsilon \eta h_{1-\xi }+{\left(\epsilon \eta \right)}^2\right]}{6}{\nabla }^2\nabla {\phi }_{\alpha }\}\}\}=0\end{array}$$E64

Similarly, substituting (63) into the third equation of (61), a momentum equation is obtained in terms of the velocity potential. Then, given that at $$z=z_{\alpha }$$, $$u_{\alpha }=\nabla {\phi }_{\alpha }$$, a fully nonlinear version of a Boussinesq type model in terms of $$\eta$$ and $$u_{\alpha }$$ is written:

$$\begin{array}{l}{\eta }_t+\nabla .\{\left(h_{1-\xi }+\epsilon \eta \right)\{u_{\alpha }+{\mu }^2\{\left[\frac{z_{\alpha }^2}{2}-\frac{1}{6}\left(h_{1-\xi }^2-\epsilon \eta h_{1-\xi }+{\left(\epsilon \eta \right)}^2\right)\right]\nabla \left(\nabla .u_{\alpha }\right)\\ +\left[z_{\alpha }+\frac{1}{2}\left(h_{1-\xi }-\epsilon \eta \right)\right]\nabla \left[\nabla .\left(h_{1-\xi }{u}_{\alpha }\right)\right]\}\}\}=O\left({\mu }^4\right)\end{array}$$E65

$$u_{\alpha }{}_t+\epsilon \left(u_{\alpha }.\nabla \right)u_{\alpha }+\nabla \eta +{\mu }^{2}R+\epsilon {\mu }^{2}S=O\left({\mu }^4\right)$$E66

where

$$R=\frac{1}{2}{z}_{\alpha }^2\nabla \left(\nabla .u_{\alpha }{}_t\right)+z_{\alpha }\nabla \left[\nabla .\left(h_{1-\xi }{u}_{\alpha }{}_t\right)\right]-\nabla \left[\frac{1}{2}{\left(\epsilon \eta \right)}^2\nabla .u_{\alpha }{}_t+\epsilon \eta \nabla .\left(h_{1-\xi }{u}_{\alpha }{}_t\right)\right]$$E67

$$\begin{array}{l}S=\nabla \left\{\left(z_{\alpha }-\epsilon \eta \right)\left(u_{\alpha }.\nabla \right)\left[\nabla .\left(h_{1-\xi }{u}_{\alpha }\right)\right]+\frac{1}{2}\left[z_{\alpha }^2-{\left(\epsilon \eta \right)}^2\right]\left(u_{\alpha }.\nabla \right)\left(\nabla .u_{\alpha }\right)\right\}\\ +\frac{1}{2}\nabla \left\{{\left[\nabla .\left(h_{1-\xi }{u}_{\alpha }\right)+\epsilon \eta \nabla .u_{\alpha }\right]}^2\right\}\end{array}$$E68

It should be noted that Nwogu’s approximation is recovered by neglecting higher-order terms. Numerical computations show that this model agrees well with solutions of the full potential problem over the range of relevant water depths.

A high-order, predictor–corrector, finite–difference numerical algorithm to solve this model was developed and is presented in [25] for the one- and two-layer models (see Sections 3.2.1 and 3.2.2). Considering one layer only, the corresponding numerical model that is the basis of the COULWAVE model is presented later, in Sections 4.1.3 and 4.1.4.

To improve the accuracy of numerical models, it has been common practice the use of high-order polynomials to approximate the vertical velocity dependence. However, this requires very elaborate and expensive numerical calculation procedures. A different approach is suggested by [26], which consists of using quadratic polynomials, matched at interfaces that divide the water column into layers. In this regard, it is worth mentioning [26] "this approach leads to a set of model equations without the high-order spatial derivatives associated with high-order polynomial approximations".

3.2.1. Mathematical model for one-layer

Defining the parameters $$S_1=\nabla .u_1$$ and $$T_1=\nabla .\left[\left(1-\xi \right)u_1\right]+\left(1/\epsilon \right)\left(\partial h_{1-\xi }/\partial t\right)$$, the model uses the following approach for the continuity equation (to compute $$\eta$$ values), in nondimensional variables:

$$\begin{array}{l}\frac{1}{\epsilon }\frac{\partial h_{1-\xi }}{\partial t}+\frac{\partial \eta }{\partial t}+\nabla .\left[\left(\epsilon \eta +h_{1-\xi }\right)u_1\right]\\ -{\mu }^2\nabla .\{\left[\frac{{\epsilon }^3{\eta }^3+h_{1-\xi }^3}{6}-\frac{\left(\epsilon \eta +h_{1-\xi }\right)k_1^2}{2}\right]\nabla S_1\\ +\left[\frac{{\epsilon }^2{\eta }^2-h_{1-\xi }^2}{2}-\left(\epsilon \eta +h_{1-\xi }\right)k_1\right]\nabla T_1\}=O\left({\mu }^4\right)\end{array}$$E69

where $$u_1$$ = horizontal velocity vector, $$k_1={\alpha }_{1}h+{\beta }_1\eta$$, $${\alpha }_1$$ and $${\beta }_1$$ are coefficients to be defined by the user. The index 1 means one-layer model. To compute the velocity components $$\left(u,v\right)$$, the following approach of the momentum equation is solved, in nondimensional variables:

$$\begin{array}{l}\frac{\partial u_1}{\partial t}+\epsilon u_1.\nabla u_1+\nabla \eta +{\mu }^2\frac{\partial }{\partial t}\left\{\frac{k_1^2}{2}\nabla S_1+k_1\nabla T_1\right\}\\ +\epsilon {\mu }^2[\left(u_1.\nabla k_1\right)\nabla T_1+k_1\nabla \left(u_1.\nabla T_1\right)+k_1\left(u_1.\nabla k_1\right)\nabla S_1\\ +\frac{k_1^2}{2}\nabla \left(u_1.\nabla S_1\right)]+\epsilon {\mu }^2\left[T_1\nabla T_1-\nabla \left(\eta \frac{\partial T_1}{\partial t}\right)\right]\\ +{\epsilon }^2{\mu }^2\nabla \left(\eta S_1{T}_1-\frac{{\eta }^2}{2}\frac{\partial S_1}{\partial t}-\eta u_1.\nabla T_1\right)\\ +{\epsilon }^2{\mu }^2\nabla \left[\frac{{\eta }^2}{2}\left(S_1^2-u_1.\nabla S_1\right)\right]=O\left({\mu }^4\right)\end{array}$$E70

The horizontal velocity vector is given as:

$$U_1=u_1-{\mu }^2\left\{\frac{z_1^2-k_1^2}{2}\nabla S_1+\left(z_1-k_1\right)\nabla T_1\right\}+O\left({\mu }^4\right)$$E71

This one-layer model, often referred to as the “fully nonlinear, extended Boussinesq equations” in the literature (e.g. [38]), has been examined and applied to a significant extent. The weakly nonlinear version of (69)–(71) (i.e., assuming $$O\left(\epsilon \right)=O\left({\mu }^2\right)$$, thereby neglecting all nonlinear dispersive terms) was first derived by [30]. Through a linear and first-order nonlinear analysis of the model equations, Nwogu recommended that $$z_1=-0.531h$$, and that value has been recommended and adopted by most researchers who use these equations. Refs. [22] and later [38] used the Nwogu’s approach to derive a high-order model which retains "full nonlinearity" by including terms not considered in the Nwogu’s approach. Comparing Liu's and Wei and Kirby's equations, there are some differences that can be attributed to the absence of some nonlinear dispersive terms in Wei and Kirby [19]. The above one-layer equations (69) and (70) are identical to those derived by [22].

3.2.2. Mathematical model for two-layers

Details about the two-layer mathematical model are presented in [25] and [24]. The model consists of a continuity equation, a momentum equation for the upper layer and a matching equation for the velocity in the lower layer. These equations are solved using an appropriate numerical method to compute values for the free surface elevation $$\eta$$ and the velocity components $$\left(u,v\right)$$. In dimensional variables, this set of equations is written ([24]):

$$\begin{array}{l}\frac{\partial \eta }{\partial t}+\nabla .\left[\left(\eta -\zeta \right)u_1+\left(\zeta +h_r\right)u_2\right]\\ -\nabla .\left\{\left[\frac{{\zeta }^3+h_r^3}{6}-\frac{\left(\zeta +h_r\right)k_2^2}{2}\right]\nabla S_2+\left[\frac{{\zeta }^2-h_r^2}{2}-\left(\zeta +h_r\right)k_2\right]\nabla T_2\right\}\\ -\nabla .\left\{\left[\frac{{\eta }^3-{\zeta }^3}{6}-\frac{\left(\eta -\zeta \right)k_1^2}{2}\right]\nabla S_1+\left[\frac{{\eta }^2-{\zeta }^2}{2}-\left(\eta -\zeta \right)k_1\right]\nabla T_1\right\}=0\end{array}$$E72

$$\begin{array}{l}\frac{\partial u_1}{\partial t}+\frac{1}{2}\nabla \left(u_1.u_1\right)+g\nabla \eta \\ +\frac{\partial }{\partial t}\left[\frac{k_1^2}{2}\nabla S_1+k_1\nabla T_1-\nabla \left(\frac{{\eta }^2}{2}{S}_1\right)-\nabla \left(\eta T_1\right)\right]\\ +\nabla [\frac{\partial \eta }{\partial t}\left(T_1+\eta S_1\right)+\left(k_1-\eta \right)\left(u_1.\nabla \right)T_1+\frac{1}{2}\left(k_1^2-{\eta }^2\right)\left(u_1.\nabla \right)S_1\\ +\frac{1}{2}{\left(T_1+\eta S_1\right)}^2]-R_b+R_f+{\nu }_T[\nabla S_1-{\nabla }^2{u}_1\\ -{\nabla }^2\left(\frac{k_1^2}{2}\nabla S_1+k_1\nabla T_1\right)+\nabla \left(\frac{{\eta }^2}{2}{\nabla }^2{S}_1+\eta {\nabla }^2{T}_1\right)]=0\end{array}$$E73

$$u_2+\frac{k_2^2-{\zeta }^2}{2}\nabla S_2+\left(k_2-\zeta \right)\nabla T_2=u_1+\frac{k_1^2-{\zeta }^2}{2}\nabla S_1+\left(k_1-\zeta \right)\nabla T_1$$E74

where $$S_1=\nabla .u_1$$, $$T_1=\zeta \left(S_2-S_1\right)+T_2$$, $$S_2=\nabla .u_2$$, and $$T_2=\nabla .\left(h_r{u}_2\right)$$. $$R_b$$ = breaking-related dissipation term, $$R_f=\left[f/\left(h_r+\eta \right)\right]u_b\left|u_b\right|$$ accounts for bottom friction, where $$u_b$$ = velocity evaluated at the seafloor, and $$f$$ = bottom friction coefficient, typically in the range of $${10}^{-3}-{10}^{-2}$$, $${\nu }_T$$ = constant eddy viscosity, $${\nabla }^2=\left({\partial }^2/\partial x^2,{\partial }^2/\partial y^2\right)$$, $$k_1=-0.127h_r$$ = evaluation level for the velocity $$u_1$$, $$\zeta =-0.266h_r$$ = layer interface elevation, s = evaluation level for the velocity $$u_2$$, and $$\eta$$ = free surface elevation.

3.3. 1HD Serre’s equations with improved dispersive performance

To allow applications in a greater range of $$h_0/\lambda$$, other than shallow waters, a new set of extended Serre equations, with additional terms of dispersive origin, is developed and tested in [1] and [2] by comparisons with the available test data. The equations are solved using an efficient finite–difference method, whose consistency and stability are tested in that work by comparison with a closed-form solitary wave solution of the Serre equations.

From the equation system (20), by adding and subtracting terms of dispersive origin, using the approximation $$u_t=-g{\eta }_x$$ and considering the parameters $$\alpha$$, $$\beta$$ and $$\gamma$$, with $$\beta =1.5\alpha -0.5\gamma$$, allows to obtain a new system of equations with improved linear dispersion characteristics:

$$\begin{array}{l}{h}_t+{\left(uh\right)}_x=0\\ u_t+u{u}_x+g{\left(h+\xi \right)}_x+\left(1+\alpha \right)\left(\mathrm{\Omega }{u}_t-h{h}_x{u}_{xt}\right)-\left(1+\beta \right)\frac{h^2}{3}{u}_{xxt}\\ +\alpha g\mathrm{\Omega }{\left(h+\xi \right)}_x-\alpha gh{h}_x{\left(h+\xi \right)}_{xx}-\beta g\frac{h^2}{3}{\left(h+\xi \right)}_{xxx}-h{h}_{x}u{u}_{xx}\\ +\frac{h^2}{3}\left(u_x{u}_{xx}-u{u}_{xxx}\right)+h{\left(u_x\right)}^2{\left(h+\xi \right)}_x+{\xi }_{xx}{u}^2{\left(h+\xi \right)}_x\\ +\left(\mathrm{\Omega }+h{\xi }_{xx}\right)u{u}_x+\frac{h}{2}{\xi }_{xxx}{u}^2=0\end{array}$$E75

After linearization of the equation system (75), the dispersion relation (58) is obtained. As for the Boussinesq approach obtained by Liu and Sun, equating equation (58) with the linear dispersion relation $${\omega }^2/gk=\tanh \left(k{h}_r\right)$$, using the approach (60), values of $$\alpha =0.1308$$ and $$\gamma =-0.0076$$ are obtained, so that $$\beta =0.20$$. It should be noted that the Serre’s equation system (20) is recovered by setting $$\alpha =\beta =0$$.

4. Numerical formulations and applications

4.1. 2HD Boussinesq-type approaches

4.1.1. WACUP numerical model

An extension of the Boussinesq model (24) to take into account wave–current interactions has been derived and presented in [5]. This model is named WACUP (a 2HD WAve Plus CUrrent Boussinesq-type model). With dimensional variables, taking mean quantities of the horizontal velocity components $$U=\left(u+u_c\right)$$ and $$V=\left(v+v_c\right)$$, where index c represents current and $$u=U-u_c$$ and $$v=V-v_c$$, the final set of these equations are written as follows:

$$h_t+{\left(hU\right)}_x+{\left(hV\right)}_y=0$$E76

$$\begin{array}{l}{U}_t+U{U}_x+V{U}_y+g{\left(h+\xi \right)}_x-\frac{h_r^2}{3}\left(U_{xxt}+V_{xyt}\right)\\ -\frac{h_r^2}{3}{\left[u_c\left(U_{xx}+V_{xy}\right)+v_c\left(U_{xy}+V_{yy}\right)\right]}_x+h_r{\left(\frac{1}{2}{\xi }_{tt}+u_c{\xi }_{xt}+v_c{\xi }_{yt}\right)}_x\\ -{\nu }_l\left(U_{xx}+U_{yy}\right)-\frac{{\tau }_{s_x}}{\rho h}+\frac{{\tau }_{b_x}}{\rho h}=0\end{array}$$E77

$$\begin{array}{l}{V}_t+U{V}_x+V{V}_y+g{\left(h+\xi \right)}_y-\frac{h_r^2}{3}\left(U_{xyt}+V_{yyt}\right)\\ -\frac{h_r^2}{3}{\left[u_c\left(U_{xx}+V_{xy}\right)+v_c\left(U_{xy}+V_{yy}\right)\right]}_y+h_r{\left(\frac{1}{2}{\xi }_{tt}+u_c{\xi }_{xt}+v_c{\xi }_{yt}\right)}_y\\ -{\nu }_l\left(V_{xx}+V_{yy}\right)-\frac{{\tau }_{s_y}}{\rho h}+\frac{{\tau }_{b_y}}{\rho h}=0\end{array}$$E78

where $${\tau }_s$$ and $${\tau }_b$$ represent stresses on surface and at bottom, respectively.

The standard Boussinesq model (24) and the extended system of equations (76)–(78) are solved in [8] and [5], respectively, using an efficient finite element method for spatial discretization of the partial differential equations. Firstly, the (U, V) derivatives in time and third spatial derivatives are grouped in two equations. This means that an equivalent system of five equations is solved instead of the original equation system (76)–(78). The final equation system takes the following form:

$$h_t+h{U}_x+U{h}_x+h{V}_y+V{h}_y=0$$E79

$$\begin{array}{l}{r}_t+u_c{r}_x+v_c{r}_y=-u{U}_x-v{U}_y-g{\left(h+\xi \right)}_x\\ +\frac{h_r^2}{3}\left[{\left(u_c\right)}_x\left(U_{xx}+U_{yy}\right)+{\left(v_c\right)}_x\left(V_{xx}+V_{yy}\right)\right]-h_r{\left(\frac{\text{1}}{\text{2}}{\xi }_{tt}+u_c{\xi }_{xt}+v_c{\xi }_{yt}\right)}_{\text{x}}\\ +\left[{\nu }_l+\frac{2}{3}{h}_r\left({\xi }_t+u_c{\xi }_x+v_c{\xi }_y\right)\right]\left(U_{xx}+U_{yy}\right)+\frac{{\tau }_{s_x}}{\rho h}-\frac{{\tau }_{b_x}}{\rho h}=0\end{array}$$E80

$$\begin{array}{l}{s}_t+u_c{s}_x+v_c{s}_y=-u{V}_x-v{V}_y-g{\left(h+\xi \right)}_y\\ +\frac{h_r^2}{3}\left[{\left(u_c\right)}_y\left(U_{xx}+U_{yy}\right)+{\left(v_c\right)}_y\left(V_{xx}+V_{yy}\right)\right]-h_r{\left(\frac{\text{1}}{\text{2}}{\xi }_{tt}+u_c{\xi }_{xt}+v_c{\xi }_{yt}\right)}_{\text{y}}\\ +\left[{\nu }_l+\frac{2}{3}{h}_r\left({\xi }_t+u_c{\xi }_x+v_c{\xi }_y\right)\right]\left(V_{xx}+V_{yy}\right)+\frac{{\tau }_{s_y}}{\rho h}-\frac{{\tau }_{b_y}}{\rho h}=0\end{array}$$E81

$$U-\frac{h_r^2}{3}\left(U_{xx}+U_{yy}\right)=r$$E82

$$V-\frac{h_r^2}{3}\left(V_{xx}+V_{yy}\right)=s$$E83

It should be noted that weakly vertical rotational flows were assumed, which strictly correspond to a limitation of the numerical method.

As the values of variables h, U, V, r and s are known at time t, we can use a numerical procedure based on the following steps to compute the corresponding values at time $$t+\mathrm{\Delta }t$$ (for details see [5]):

1. The equation (79) allows us to predict the values of variable h($$h_p^{t+\mathrm{\Delta }t}$$), considering the known values of h, U and V at time t in the whole domain.

2. Equations (80) and (81) make it possible to predict the values of variables r ($$r_p^{t+\mathrm{\Delta }t}$$) and s ($$s_p^{t+\mathrm{\Delta }t}$$), taking into account the values of $$U^t$$, $$V^t$$, $$r^t$$, $$s^t$$ and $${\tilde{h}}^{t+\mathrm{\Delta }t}=0.5h^t+0.5h_p^{t+\mathrm{\Delta }t}$$ known for the whole domain.

3. Solutions of equations (82) and (83) give us the values of the mean-averaged velocity components U and V ($$U^{t+\mathrm{\Delta }t}$$ and $$V^{t+\mathrm{\Delta }t}$$), taking into account the predicted values of r and s ($$r_p^{t+\mathrm{\Delta }t}$$ and $$s_p^{t+\mathrm{\Delta }t}$$, respectively).

4. Equation (79) allows us to compute the depth h at time $$t+\mathrm{\Delta }t$$ (values of $$h^{t+\mathrm{\Delta }t}$$) considering the values of variables $$h^t$$, $$U^{t+05\mathrm{\Delta }t}=0.5U^t+0.5U^{t+\mathrm{\Delta }t}$$ and $$V^{t+05\mathrm{\Delta }t}=0.5V^t+0.5V^{t+\mathrm{\Delta }t}$$ known for the whole domain.

5. Equations (80) and (81) allow us to compute the values of variables r and s at time $$t+\mathrm{\Delta }t$$ (values of $$r^{t+\mathrm{\Delta }t}$$ and $$s^{t+\mathrm{\Delta }t}$$), taking into account the values of $$r^t$$, $$s^t$$ $$h^{t+05\mathrm{\Delta }t}=0.5h^t+0.5h^{t+\mathrm{\Delta }t}$$, $$U^{t+05\mathrm{\Delta }t}=0.5U^t+0.5U^{t+\mathrm{\Delta }t}$$ and $$V^{t+05\mathrm{\Delta }t}=0.5V^t+0.5V^{t+\mathrm{\Delta }t}$$ known for the whole domain.

The Petrov-Galerkin procedure is utilized to achieve solutions for the unknowns h, r and s. According to the weighted residual technique, minimization requires the “orthogonality” of the residual $$R_J$$ to a set of weighting functions $$W_{i,J}$$, i.e.,

$${\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,J}{R}_J}d{\mathrm{\Delta }}^e=0,i=1,\cdots ,n$$E84

where the general form of the weighting functions applied to these equations is defined as:

$$W_{i,J}=N_i+{\delta }_u{}_i{\left(N_x\right)}_i+{\delta }_v{}_i{\left(N_y\right)}_i,i=1,\cdots ,n$$E85

and where the coefficients $${\delta }_u{}_i$$, and $${\delta }_v{}_i$$ are functions of: (i) the local velocities U and V; (ii) the ratio of the wave amplitude to the water depth, and (iii) the element length. To illustrate this procedure, a complete solution of equations (79) (80) and (82) is presented here in detail. Introducing in equation (79), the approximated values are given by:

$$p\approx \widehat{p}={\displaystyle \sum _{1=1}^n{N}_i{p}_i}$$E86

the following residual $$R_{79}$$ is obtained:

$$R_{79}={\widehat{h}}_t+\widehat{h}{\widehat{U}}_x+\widehat{U}{\widehat{h}}_x+\widehat{h}{\widehat{V}}_y+\widehat{V}{\widehat{h}}_y$$E87

The $$R_{79}$$ error minimization leads to the following equation:

$$\begin{array}{l}{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}{R}_{79}}d{\mathrm{\Delta }}^e={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}\left({\widehat{h}}_t+\widehat{h}{\widehat{U}}_x+\widehat{U}{\widehat{h}}_x+\widehat{h}{\widehat{V}}_y+\widehat{V}{\widehat{h}}_y\right)}d{\mathrm{\Delta }}^e\\ ={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}[{\displaystyle \sum _{j=1}^n{N}_j{\left(h_t\right)}_j}+{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k{U}_k}{\displaystyle \sum _{j=1}^n{N}_j}{h}_j+{\displaystyle \sum _{k=1}^n{N}_k{U}_k}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{h}_j}\\ +{\displaystyle \sum _{k=1}^n{\left(N_y\right)}_k{V}_k}{\displaystyle \sum _{j=1}^n{N}_j}{h}_j+{\displaystyle \sum _{k=1}^n{N}_k{V}_k}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{h}_j]d{\mathrm{\Delta }}^e=0\end{array}$$E88

In the matrix form, this equation may be written as:

$$A{h}_t+Bh=0$$E89

where

$$\begin{array}{l}{a}_{i,j}={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}{N}_j}d{\mathrm{\Delta }}^e;i,j=1,\cdots ,n\\ b_{i,j}={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}{\displaystyle \sum _{k=1}^n{N}_k{U}_k{\left(N_x\right)}_j}}d{\mathrm{\Delta }}^e+{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k{U}_k{N}_j}}d{\mathrm{\Delta }}^e\\ +{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}{\displaystyle \sum _{k=1}^n{\left(N_y\right)}_k{V}_k{N}_j}}d{\mathrm{\Delta }}^e+{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,79}{\displaystyle \sum _{k=1}^n{N}_k{V}_k{\left(N_y\right)}_j}}d{\mathrm{\Delta }}^e\text{; }i,j=1,\cdots ,n\end{array}$$E90

The residual $$R_{80}$$ is written as:

$$\begin{array}{l}{R}_{80}={\widehat{r}}_t+{\widehat{u}}_c{\widehat{r}}_x+{\widehat{v}}_c{\widehat{r}}_y+\widehat{u}{\widehat{U}}_x+\widehat{v}{\widehat{U}}_y\\ +g{\left(\widehat{h}+\widehat{\xi }\right)}_x-\frac{{\widehat{h}}_r^2}{3}\left[{\left({\widehat{u}}_c\right)}_x\left({\widehat{U}}_{xx}+{\widehat{U}}_{yy}\right)+{\left({\widehat{v}}_c\right)}_x\left({\widehat{V}}_{xx}+{\widehat{V}}_{yy}\right)\right]\\ +{\widehat{h}}_r\left[\frac{1}{2}{\left({\widehat{\xi }}_{tt}\right)}_x+{\left({\widehat{u}}_c\right)}_x{\left({\widehat{\xi }}_t\right)}_x+{\widehat{u}}_c{\left({\widehat{\xi }}_t\right)}_{xx}+{\left({\widehat{v}}_c\right)}_x{\left({\widehat{\xi }}_t\right)}_y+{\widehat{v}}_c{\left({\widehat{\xi }}_t\right)}_{xy}\right]\\ -{\nu }_l\left({\widehat{U}}_{xx}+{\widehat{U}}_{yy}\right)-\frac{{\tau }_{s_x}}{\rho h}+\frac{{\tau }_{b_x}}{\rho h}\end{array}$$E91

It should be noted that the term $$\frac{2}{3}{h}_r\left({\xi }_t+u_c{\xi }_x+v_c{\xi }_y\right)\left(V_{xx}+V_{yy}\right)$$ of equation (80) is of order $${\sigma }^4$$ or greater. For this reason, it is not considerd in the numerical developments. Similarly, and for the same reason, considering $$f=h_r^2/3$$, all terms involving $$\nabla f$$ in the numerical developments are omitted.

The Green’s theorem is used to solve the second derivatives present in equation (80) (residual $$R_{80}$$), and in equations (81) to (83), i.e., considering $$\widehat{p}=\left(\widehat{U},\widehat{V}\right)$$, from $${\widehat{p}}_{xx}+{\widehat{p}}_{yy}$$ we obtain:

$${\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,l}\left({\widehat{p}}_{xx}+{\widehat{p}}_{yy}\right)}d{\mathrm{\Delta }}^e={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,l}{\displaystyle \sum _{\text{j}=1}^{\text{n}}\left[{\left(N_{xx}\right)}_j+{\left(N_{yy}\right)}_j\right]}}{p}_{j}d{\mathrm{\Delta }}^e$$E92

and so

$$\begin{array}{l}{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,l}{\displaystyle \sum _{\text{j}=1}^{\text{n}}\left[{\left(N_{xx}\right)}_j+{\left(N_{yy}\right)}_j\right]}}{p}_{j}d{\mathrm{\Delta }}^e=-{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{\left(W_x\right)}_{i,l}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j{p}_j}d{\mathrm{\Delta }}^e}\\ -{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{\left(W_y\right)}_{i,l}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{p}_{j}d{\mathrm{\Delta }}^e}+{\displaystyle {\oint }_{{\mathrm{\Gamma }}^e}{W}_{i,l}{\displaystyle \sum _{j=1}^n{\left(N_n\right)}_j}{p}_{j}d{\mathrm{\Gamma }}^e}\text{ ; }p=\left(U,V\right)\end{array}$$E93

Retaining terms up to order $${\sigma }^3$$, the $$R_{80}$$ error minimization leads to:

$$\begin{array}{l}{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,80}{R}_{80}}d{\mathrm{\Delta }}^e={\displaystyle {\int }_{{\mathrm{\Delta }}^e}\{W_{i,80}[{\displaystyle \sum _{j=1}^n{N}_j{\left(r_t\right)}_j}+{\displaystyle \sum _{k=1}^n{N}_k{\left(u_c\right)}_k}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{r}_j}\\ +{\displaystyle \sum _{k=1}^n{N}_k{\left(v_c\right)}_{k_x}}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{r}_j+{\displaystyle \sum _{k=1}^n{N}_k{u}_k}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{U}_j+{\displaystyle \sum _{k=1}^n{N}_k{v}_k}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{U}_j\\ +g{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{\left(h+\xi \right)}_j]-{\displaystyle \sum _{l=1}^n{N}_l}\frac{{\left(h_r^2\right)}_l}{3}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k}\left[{\left(u_c\right)}_k{U}_{x^2{y}^2}+{\left(v_c\right)}_k{V}_{x^2{y}^2}\right]\\ +W_{i,80}[{\displaystyle \sum _{k=1}^n{N}_k\frac{{\left(h_r\right)}_k}{2}}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{\left({\xi }_{tt}\right)}_j+{\displaystyle \sum _{l=1}^n{N}_l{\left(h_r\right)}_l}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k}{\left(u_c\right)}_k{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{\left({\xi }_t\right)}_j\\ +{\displaystyle \sum _{l=1}^n{N}_l{\left(h_r\right)}_l}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k}{\left(v_c\right)}_k{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{\left({\xi }_t\right)}_j]\\ +{\displaystyle \sum _{l=1}^n{N}_l{\left(h_r\right)}_l}\left[{\displaystyle \sum _{k=1}^n{N}_k}{\left(u_c\right)}_k{\left({\tilde{\xi }}_t\right)}_{xx}+{\displaystyle \sum _{k=1}^n{N}_k}{\left(v_c\right)}_k{\left({\tilde{\xi }}_t\right)}_{xy}\right]-{\nu }_l{U}_{x^2{y}^2}\\ -W_{i,80}{\displaystyle \sum _{k=1}^n{N}_k\left[{\left({\tau }_{s_x}\right)}_k-{\left({\tau }_{b_x}\right)}_k\right]}/{\displaystyle \sum _{j=1h}^n{N}_j{h}_j}\}d{\mathrm{\Delta }}^e=0\end{array}$$E94

where

$${\left({\tilde{\xi }}_t\right)}_{xx}=-{\left(W_x\right)}_{i,80}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j{\left({\xi }_t\right)}_j},{\left({\tilde{\xi }}_t\right)}_{xy}=-{\left(W_x\right)}_{i,80}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j{\left({\xi }_t\right)}_j}$$E95

and

$$p_{x^2{y}^2}=-\left[{\left(W_x\right)}_{i,J}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}+{\left(W_y\right)}_{i,J}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}\right]p_j,p=\left(U,V\right)$$E96

In the matrix form, this equation may be written as:

$$A{r}_t+Br=C$$E97

where

$$a_{i,j}={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,80}{N}_j}d{\mathrm{\Delta }}^e,i,j=1,\cdots ,n$$E98

$$b_{i,j}={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,80}\left[{\displaystyle \sum _{k=1}^n{N}_k{\left(u_c\right)}_k{\left(N_x\right)}_j}+{\displaystyle \sum _{k=1}^n{N}_k{\left(v_c\right)}_k{\left(N_y\right)}_j}\right]}d{\mathrm{\Delta }}^e,i,j=1,\cdots ,n$$E99

$$\begin{array}{l}{c}_i=-{\displaystyle {\int }_{{\mathrm{\Delta }}^e}\{W_{i,80}[{\displaystyle \sum _{k=1}^n{N}_k{u}_k}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{U}_j}+{\displaystyle \sum _{k=1}^n{N}_k{v}_k}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{U}_j\\ +g{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{\left(h+\xi \right)}_j]-{\displaystyle \sum _{l=1}^n{N}_l}\frac{{\left(h_r^2\right)}_l}{3}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k}\left[{\left(u_c\right)}_k{U}_{x^2{y}^2}+{\left(v_c\right)}_k{V}_{x^2{y}^2}\right]\\ +W_{i,80}[{\displaystyle \sum _{k=1}^n{N}_k\frac{{\left(h_r\right)}_k}{2}}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{\left({\xi }_{tt}\right)}_j+{\displaystyle \sum _{l=1}^n{N}_l{\left(h_r\right)}_l}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k}{\left(u_c\right)}_k{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}{\left({\xi }_t\right)}_j\\ +{\displaystyle \sum _{l=1}^n{N}_l{\left(h_r\right)}_l}{\displaystyle \sum _{k=1}^n{\left(N_x\right)}_k}{\left(v_c\right)}_k{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}{\left({\xi }_t\right)}_j]\\ +{\displaystyle \sum _{l=1}^n{N}_l{\left(h_r\right)}_l}\left[{\displaystyle \sum _{k=1}^n{N}_k}{\left(u_c\right)}_k{\left({\tilde{\xi }}_t\right)}_{xx}+{\displaystyle \sum _{k=1}^n{N}_k}{\left(v_c\right)}_k{\left({\tilde{\xi }}_t\right)}_{xy}\right]-{\nu }_l{U}_{x^2{y}^2}\\ -W_{i,80}{\displaystyle \sum _{k=1}^n{N}_k\left[{\left({\tau }_{s_x}\right)}_k-{\left({\tau }_{b_x}\right)}_k\right]}/{\displaystyle \sum _{j=1h}^n{N}_j{h}_j}\}d{\mathrm{\Delta }}^e=0\text{ , }i=1,\dots .n\end{array}$$E100

The residual $$R_{82}$$ is written as:

$$R_{82}=\widehat{U}-\frac{{\widehat{h}}_r^2}{3}\left(\frac{{\partial }^2\widehat{U}}{\partial x^2}+\frac{{\partial }^2\widehat{U}}{\partial y^2}\right)-\widehat{r}$$E101

According to Galerkin’s procedure, after using integration by parts (or Green's theorem) to reduce the second derivatives, the $$R_{82}$$ error minimization leads to the following equation:

$$\begin{array}{l}{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,82}{R}_{82}}d{\mathrm{\Delta }}^e={\displaystyle {\int }_{{\mathrm{\Delta }}^e}\{W_{i,82}{\displaystyle \sum _{j=1}^n{N}_j}{U}_j}\\ +{\displaystyle \sum _{k=1}^n{N}_k\frac{{\left(h_r^2\right)}_k}{3}\left[{\left(W_x\right)}_{i,82}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}+{\left(W_y\right)}_{i,82}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}\right]U_j}\\ -W_{i,82}{\displaystyle \sum _{j=1}^n{N}_j}{r}_j\}d{\mathrm{\Delta }}^e-{\displaystyle {\oint }_{{\mathrm{\Gamma }}^e}{N}_k}\frac{{\left(h_r^2\right)}_k}{3}{W}_{i,82}{\left(N_n\right)}_j{U}_{j}d{\mathrm{\Gamma }}^e\end{array}$$E102

The last term of (102) can be written as:

$$\begin{array}{l}{\displaystyle {\oint }_{{\mathrm{\Gamma }}^e}{N}_k}\frac{{\left(h_r^2\right)}_k}{3}{W}_{i,82}{\left(N_n\right)}_j{U}_{j}d{\mathrm{\Gamma }}^e={\displaystyle {\oint }_{{\mathrm{\Gamma }}_i^e}{N}_p}\frac{{\left(h_r^2\right)}_p}{3}{W}_{q,82}{\left(N_n\right)}_r{U}_{r}d{\mathrm{\Gamma }}_i^e\\ +{\displaystyle {\oint }_{{\mathrm{\Gamma }}_e^e}{N}_p}\frac{{\left(h_r^2\right)}_p}{3}{W}_{q,82}{U}_{n}d{\mathrm{\Gamma }}_e^e\end{array}$$E103

where $$p,q=1,\cdots ,n^e$$, $$n^e$$ being the number of nodes in the corresponding element side coincident with the boundary domain, $${\mathrm{\Gamma }}_i^e$$ represents the element sides within the domain, with the corresponding integral null because the resulting element contributions are equal, but with opposite signals, and $${\mathrm{\Gamma }}_e^e$$ represents the element sides coincident with the boundary domain. Accordingly, an equivalent form of equation (102) may be written as follows:

$$\begin{array}{l}{\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,82}{R}_{82}}d{\mathrm{\Delta }}^e={\displaystyle {\int }_{{\mathrm{\Delta }}^e}\{W_{i,82}{\displaystyle \sum _{j=1}^n{N}_j}}\\ +{\displaystyle \sum _{k=1}^n{N}_k\frac{h{\left({}_r^2\right)}_k}{3}\left[{\left(W_x\right)}_{i,82}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}+{\left(W_y\right)}_{i,82}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}\right]}\}{U}_{j}d{\mathrm{\Delta }}^e\\ ={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,82}{\displaystyle \sum _{j=1}^n{N}_j}{r}_j}d{\mathrm{\Delta }}^e+{\displaystyle {\oint }_{{\mathrm{\Gamma }}^e}{N}_p}\frac{{\left(h_r^2\right)}_p}{3}{W}_{q,82}{U}_{n}d{\mathrm{\Gamma }}_e^e\\ p,q=1,\cdots ,n^e\end{array}$$E104

In the matrix form, this equation may be written as:

$$AU=B$$E105

where

$$a_{i,j}={\displaystyle {\int }_{{\mathrm{\Delta }}^e}\{W_{i,82}{\displaystyle \sum _{j=1}^n{N}_j}}+{\displaystyle \sum _{k=1}^n{N}_k\frac{{\left(h_r^2\right)}_k}{3}\left[{\left(W_x\right)}_{i,82}{\displaystyle \sum _{j=1}^n{\left(N_x\right)}_j}+{\left(W_y\right)}_{i,82}{\displaystyle \sum _{j=1}^n{\left(N_y\right)}_j}\right]}\}d{\mathrm{\Delta }}^e$$E106

$$\begin{array}{l}{b}_i={\displaystyle {\int }_{{\mathrm{\Delta }}^e}{W}_{i,82}{\displaystyle \sum _{j=1}^n{N}_j}{r}_j}d{\mathrm{\Delta }}^e+{\displaystyle {\oint }_{{\mathrm{\Gamma }}^e}{N}_p}\frac{{\left(h_r^2\right)}_p}{3}{W}_{q,82}{U}_{n}d{\mathrm{\Gamma }}^e\\ p,q=1,\cdots ,n^e.\end{array}$$E107

As recommended in [5] “a suitable grid is normally crucial to the success of a finite element model. In our case, the following rules must be fulfilled for its generation.

1. Element side lower than the local depth.

2. Minimum of 20 to 25 elements per wave length.

3. Courant number always lower than one in the whole domain”.

4.1.2. Real case study using the WACUP model

The fortification of S. Lourenço da Cabeça Seca (lighthouse of Bugio) – Tagus estuary (Portugal) – has endured over the course of four centuries the continuous action of waves and currents, as well as bathymetric modifications resulting from the movement of significant quantities of sand in the area where it is located.

With the intention of preventing the destruction of this fortification, several studies were conducted to evaluate the best protection structure. The studies have led to a protective structure which consists of a circular dissipation platform, with a level of 2 m (HZ) and about 80 m radius (Figure 1).

Figure 1.

Submerged dissipation platform to protect the Bugio lighthouse, situated at the mouth of the Tagus estuary.

The wave–current Boussinesq-type model WACUP was used to obtain the time-dependent hydrodynamic characteristics of the joint action of a relatively common wave over the high flow tide current. The wave characteristics are: wave height, H = 3.5 m, period T = 12 s and direction = 180°. Figure 2 shows the free surface level when the wave approaches the platform and its action on the Bugio lighthouse. Figure 3 shows a free surface water level, in a quasi-stationary state.

As can be seen, the wave action on the fortification was drastically reduced as a consequence of the wave breaking, refraction and diffraction on the dissipation platform constructed around the fortification.

Figure 2.

Perspective of the free surface obtained by simulation around the Bugio lighthouse, situated in the Tagus estuary, Portugal, in a transient condition.

Figure 3.

Protection of the Bugio lighthouse situated at the mouth of the Tagus estuary. Perspective of the free surface obtained with a numerical simulation, in a quasi-steady state.

4.1.3. COULWAVE one-layer model

A numerical scheme similar to that of [38] and [39] is utilized by [25], with the inclusion of extra nonlinear terms. Basically using the same high-order predictor–corrector scheme, [25] developed a numerical code (COULWAVE) based on Nwogu’s equations for one and two layers. Parameterizations of bottom friction and wave breaking have been included in the code, as well as a moving boundary scheme to simulate wave runup and rundown. A finite–difference algorithm is used for the general one- and two-layer model equations.

According to [9], “the equations are solved using a high-order predictor-corrector scheme, employing a third order in time explicit Adams-Bashforth predictor step, and a fourth order in time Adams-Moulton implicit corrector step” [31]. Also in accordance with [9], “the implicit corrector step must be iterated until a convergence criterion is satisfied”.

In order to solve numerically the nondimensional equations (69)–(71), these are previously rewritten in dimensional variables. Then, to simplify the predictor–corrector equations, the velocity time derivatives in the momentum equations are grouped into the dimensional form (for details see [25]):

$$U=u+\frac{k^2-{\eta }^2}{2}{u}_{xx}+\left(k-\eta \right){\left(h_{r}u\right)}_{xx}-{\eta }_x\left[\eta u_x+{\left(h_{r}u\right)}_x\right]$$E108

$$V=v+\frac{k^2-{\eta }^2}{2}{v}_{yy}+\left(k-\eta \right){\left(h_{r}v\right)}_{yy}-{\eta }_y\left[\eta v_y+{\left(h_{r}v\right)}_y\right]$$E109

where subscripts denote partial derivatives. For reasons of stability and less iterations required in the process of convergence, the nonlinear time derivatives, arisen from the nonlinear dispersion terms $$\nabla \left[\eta \left(\nabla .{\left(h_r{u}_{\alpha }\right)}_t+h_r{}_{tt}/\epsilon \right)\right]$$ and $$\nabla \left[\left({\eta }^2/2\right)\nabla .u_{\alpha t}\right]$$, can be reformulated using the relations:

$$\nabla \left[\eta \left(\nabla .{\left(h_r{u}_{\alpha }\right)}_t+\frac{h_{rtt}}{\epsilon }\right)\right]=\nabla {\left[\eta \left(\nabla .\left(h_r{u}_{\alpha }\right)+\frac{h_{rt}}{\epsilon }\right)\right]}_t-\nabla \left[{\eta }_t\left(\nabla .\left(h_r{u}_{\alpha }\right)+\frac{h_{rt}}{\epsilon }\right)\right]$$E110

$$\nabla \left(\frac{{\eta }^2}{2}\nabla .u_{\alpha t}\right)=\nabla {\left(\frac{{\eta }^2}{2}\nabla .u_{\alpha }\right)}_t-\nabla \left(\eta {\eta }_t\nabla .u_{\alpha }\right)$$E111

The explicit predictor equations are:

$${\eta }_{i,j}^{n+1}={\eta }_{i,j}^n+\frac{\mathrm{\Delta }t}{12}\left(23E_{i,j}^n-16E_{i,j}^{n-1}+5E_{i,j}^{n-2}\right)$$E112

$$U_{i,j}^{n+1}=U_{i,j}^n+\frac{\mathrm{\Delta }t}{12}\left(23F_{i,j}^n-16F_{i,j}^{n-1}+5F_{i,j}^{n-2}\right)+2{\left(F_1\right)}_{i,j}^n-3{\left(F_1\right)}_{i,j}^{n-1}+{\left(F_1\right)}_{i,j}^{n-2}$$E113

$$V_{i,j}^{n+1}=V_{i,j}^n+\frac{\mathrm{\Delta }t}{12}\left(23G_{i,j}^n-16G_{i,j}^{n-1}+5G_{i,j}^{n-2}\right)+2{\left(G_1\right)}_{i,j}^n-3{\left(G_1\right)}_{i,j}^{n-1}+{\left(G_1\right)}_{i,j}^{n-2}$$E114

where

$$\begin{array}{l}E=-h_r{}_t-{\left[\left(\eta +h_r\right)u\right]}_x-{\left[\left(\eta +h_r\right)v\right]}_y\\ +{\left\{\left(\eta +h_r\right)\left[\left(\frac{1}{6}\left({\eta }^2-\eta h_r+h_r^2\right)-\frac{1}{2}{k}^2\right)S_x+\left(\frac{1}{2}\left(\eta -h_r\right)-k\right)T_x\right]\right\}}_x\\ +{\left\{\left(\eta +h_r\right)\left[\left(\frac{1}{6}\left({\eta }^2-\eta h_r+h_r^2\right)-\frac{1}{2}{k}^2\right)S_y+\left(\frac{1}{2}\left(\eta -h_r\right)-k\right)T_y\right]\right\}}_y\end{array}$$E115

$$\begin{array}{l}F=-\frac{1}{2}\left[{\left(u^2\right)}_x+{\left(v^2\right)}_x\right]-g{\eta }_x-k{h}_r{}_{xtt}-k_t{h}_r{}_{xt}+{\left(E{h}_r{}_t+\eta h_r{}_{tt}\right)}_x\\ -{\left[E\left(\eta S+T\right)\right]}_x-{\left[\frac{1}{2}\left(k^2-{\eta }^2\right)\left(u{S}_x+v{S}_y\right)\right]}_x\\ -{\left[\left(k-\eta \right)\left(u{T}_x+v{T}_y\right)\right]}_x-\frac{1}{2}{\left[{\left(T+\eta S\right)}^2\right]}_x\end{array}$$E116

$$F_1=\frac{{\eta }^2-k^2}{2}{v}_{xy}-\left(k-\eta \right){\left(h_{r}v\right)}_{xy}+{\eta }_x\left[\eta v_y+{\left(h_{r}v\right)}_y\right]$$E117

$$\begin{array}{l}G=-\frac{1}{2}\left[{\left(u^2\right)}_y+{\left(v^2\right)}_y\right]-g{\eta }_y-k{h}_r{}_{ytt}-k_t{h}_r{}_{yt}+{\left(E{h}_r{}_t+\eta h_r{}_{tt}\right)}_y\\ -{\left[E\left(\eta S+T\right)\right]}_y-{\left[\frac{1}{2}\left(k^2-{\eta }^2\right)\left(u{S}_x+v{S}_y\right)\right]}_y\\ -{\left[\left(k-\eta \right)\left(u{T}_x+v{T}_y\right)\right]}_y-\frac{1}{2}{\left[{\left(T+\eta S\right)}^2\right]}_y\end{array}$$E118

$$G_1=\frac{{\eta }^2-k^2}{2}{u}_{xy}-\left(k-\eta \right){\left(h_{r}u\right)}_{xy}+{\eta }_y\left[\eta u_x+{\left(h_{r}u\right)}_x\right]$$E119

and

$$S=u_x+v_y;T=h_r{}_t+{\left(h_{r}u\right)}_x+{\left(h_{r}v\right)}_y$$E120

All first-order spatial derivatives are differenced with fourth-order ($$\mathrm{\Delta }{x}^4=\mathrm{\Delta }{y}^4$$) accurate equations, which are five-point differences. Second-order spatial derivatives are approximated with three-point centered finite-difference equations, which are second-order accurate. The fourth-order implicit corrector expressions for the free surface elevation, $$\eta$$, and the horizontal velocities, $$u$$ and $$v$$, are:

$${\eta }_{i,j}^{n+1}={\eta }_{i,j}^n+\frac{\mathrm{\Delta }t}{24}\left(9E_{i,j}^{n+1}+19E_{i,j}^n-5E_{i,j}^{n-1}+E_{i,j}^{n-2}\right)$$E121

$$U_{i,j}^{n+1}=U_{i,j}^n+\frac{\mathrm{\Delta }t}{24}\left(9F_{i,j}^{n+1}+19F_{i,j}^n-5F_{i,j}^{n-1}+F_{i,j}^{n-2}\right)+{\left(F_1\right)}_{i,j}^{n+1}-{\left(F_1\right)}_{i,j}^n$$E122

$$V_{i,j}^{n+1}=V_{i,j}^n+\frac{\mathrm{\Delta }t}{24}\left(9G_{i,j}^{n+1}+19G_{i,j}^n-5G_{i,j}^{n-1}+G_{i,j}^{n-2}\right)+{\left(G_1\right)}_{i,j}^{n+1}-{\left(G_1\right)}_{i,j}^n$$E123

As noted in [25], “the system is solved by first evaluating the predictor equations, then u and v are solved via (108) and (109), respectively. Both (108) and (109) yield a diagonal matrix after finite differencing. The matrices are diagonal, with a bandwidth of three (due to three-point finite differencing), and the efficient Thomas algorithm can be utilized. At this point in the numerical system, we have predictors for $$\eta$$, $$u$$ and $$v$$ ”. Next, the corrector expressions are evaluated, and again u and v are determined from (108) and (109).

Also in accordance with [25], “the error is calculated, in order to determine if the implicit correctors need to be reiterated. The error criteria employed is a dual calculation, and requires that either“

$$\max \left|\frac{w^{n+1}-w_{*}^{n+1}}{w^{n+1}}\right|<\frac{\epsilon }{100}\text{or}\frac{{\displaystyle \sum \left|w^{n+1}-w_{*}^{n+1}\right|}}{{\displaystyle \sum \left|w^{n+1}\right|}}<\epsilon$$E124

In these expressions, w represents any of the variables $$\eta$$, u and v, and $$w_{\ast }$$ is the previous iterations value. The value of the error is set to $${10}^{-6}$$. Linear stability analysis performed by [38, 18] and [40] show that $$\mathrm{\Delta }t<\mathrm{\Delta }x/\left(2c\right)$$ to ensure stability, where c is the celerity.

4.1.4. Research study using COULWAVE model

For the analysis concerning coastal protection, the mean currents around a submerged structure (artificial reef) are analysed. The output of the COULWAVE model corresponds to the velocity values at a depth 0.531 h_r below the water surface. The velocity at this depth is used to determine the velocity cells near the shoreline that could give an indication of the sediment transport. Divergent cells indicate erosion near by the shoreline and convergent cells indicate sedimentation.

The numerical simulations to study the 2HD behavior of the hydrodynamics around the reef have been done for four cases: two reef geometries, varying the reef angle (45^o and 66^o), and two wave conditions. The characteristics of the simulations for the different cases are described in Table 1.

	Reef angle; (°)	H; (m)	T; (s)	Number of grid points per wavelength	Grid size (m)	Time step (s)
C1	45	4.0	15	60	2.77	0.11999
C2	66	4.0	15	70	2.37	0.10285
C3	45	1.5	9	43	2.14	0.09288
C4	66	1.5	9	43	2.14	0.09288

Table 1.

Main characteristics of the simulations performed.

As recorded in [28], “the computational domain is around 1870 m in the long-shore and 1670 m in the cross-shore direction with a constant node spacing of around $$\mathrm{\Delta }x$$ = $$\mathrm{\Delta }y$$ =2.0 m and a Courant number of 0.5. The total simulation time was 800 s. A flat bottom is placed in front of the slope where waves are generated using the source function method ([39] Wei et al., 1995). The source function is located at x=80 m along the y direction” (Figure 4).

Figure 4.

Schematic representation of the simulated geometry (not at scale) [28].

Two sponge layers are used, one in front of the offshore boundary to absorb the outgoing wave energy, and the other on the beach, both with a width of 0.5 times the wavelength of the incident wave. The numerical results obtained by the model are the time series of the free surface elevation, the two velocity components, u and v, and the wave breaking areas (Figure 5).

Figure 5.

Velocity patterns of cases C1 to C4 [28].

Results of this simulation are described in [29], including the observed phenomena, such as “along the reef, an increase of the wave height is observed, in opposition to the situation without the reef, owing to the decrease of the depth in the reef zone. Moreover, owing to the increase in wave heights, the wave breaking occurs earlier (and in general over the reef) in comparison with the situation without a reef. From Figure 5, it is clear that the presence of the reef significantly alters the wave heights, with the wave height increasing along the reef as a consequence of decreasing water depth”.

Figure 5 also shows for all cases that convergent cells appear, indicating a possible sedimentation near the shoreline, suggesting that the chosen geometries (Figure 4) are advantageous for both coastal protection and improving surf conditions. Anyway, as referred in [28], a morphological study should be done in order to confirm these results.

4.2. 2HD Serre’s standard model

4.2.1. Numerical formulation

The equation system (19) is solved in [7] using an explicit finite–difference method based on the MacCormack time-splitting scheme. For this purpose, the equations are written in the following form:

$$h_t+P_x+Q_y=0$$E125

$$\begin{array}{l}{P}_t+{\left(uP\right)}_x+{\left(vP\right)}_y+{\left\{\left[\left(g+\beta \right)/2+\alpha /3\right]h^2\right\}}_x\\ =-\left(g+\beta +\alpha /2\right)h{\xi }_x-{\tau }_{bx}+Rdiv\left(hgradu\right)\end{array}$$E126

$$\begin{array}{l}{Q}_t+{\left(uQ\right)}_x+{\left(vQ\right)}_y+{\left\{\left[\left(g+\beta \right)/2+\alpha /3\right]h^2\right\}}_y\\ =-\left(g+\beta +\alpha /2\right)h{\xi }_y-{\tau }_{by}+Rdiv\left(hgradv\right)\end{array}$$E127

where $$P=hu$$, $$Q=hv$$, $$\alpha =d^{2}h/d{t}^2$$; $$\text{β}={\text{d}}^{\text{2}}\text{ξ}/{\text{dt}}^{\text{2}}$$ and the bottom friction terms, $${\tau }_{bx}$$ and $${\tau }_{by}$$, are obtained through (128):

$${\tau }_{bx}=g\frac{P\sqrt{P^2+Q^2}}{K^2{h}^{7/3}}\text{and}{\tau }_{by}=g\frac{Q\sqrt{P^2+Q^2}}{K^2{h}^{7/3}}$$E128

In order to apply the MacCormack’s method, equations (125) –(127) are split into two systems of three equations throughout the Ox and Oy directions. The corresponding operators $$Lx$$ and $$Ly$$ take the following form [7]:

Operator $$Lx$$

$$h_t+P_x=0$$E129

$$P_t+{\left(uP\right)}_x+{\left\{\left[\left(g+\beta \right)/2+\alpha /3\right]h^2\right\}}_x=-\left(g+\beta +\alpha /2\right)h{\xi }_x-{\tau }_{bx}+R{P}_{xx}$$E130

$$Q_t+{\left(uQ\right)}_x=R{Q}_{xx}$$E131

Operator $$Ly$$

$$h_t+Q_y=0$$E132

$$P_t+{\left(vP\right)}_y=R{P}_{yy}$$E133

$$Q_t+{\left(vQ\right)}_y+{\left\{\left[\left(g+\beta \right)/2+\alpha /3\right]h^2\right\}}_y=-\left(g+\beta +\alpha /2\right)h{\xi }_y-{\tau }_{by}+R{Q}_{yy}$$E134

Considering the generic variable $$F$$, the solution at time $$\left(n+1\right)\mathrm{\Delta }t$$, for the computational point $$\left(i,j\right)$$, is obtained from the known solution $$F_{i,j}^n$$ through the following symmetric application:

$$F_{i,j}^{n+1}=Lx\left(\frac{\mathrm{\Delta }t}{4}\right)Ly\left(\frac{\mathrm{\Delta }t}{4}\right)Lx\left(\frac{\mathrm{\Delta }t}{4}\right)Ly\left(\frac{\mathrm{\Delta }t}{4}\right)Ly\left(\frac{\mathrm{\Delta }t}{4}\right)Lx\left(\frac{\mathrm{\Delta }t}{4}\right)Ly\left(\frac{\mathrm{\Delta }t}{4}\right)Lx\left(\frac{\mathrm{\Delta }t}{4}\right)F_{i,j}^n$$E135

where each operator, $$Lx$$ and $$Ly$$, is composed of a predictor–corrector sequence and $$n$$ represents a generic time $$t$$. In the application (135) of eight predictor–corrector sequences, alternately backward and forward space differences are used. After each predictor and each corrector of the application $$F$$, the values of the velocities $$\left(u,v\right)$$ are updated and the values of the vertical accelerations, $$\alpha$$ and $$\beta$$, are recalculated (for details see [7]).

4.2.2. Real case study

This model was applied to compute the agitation established in the port of Figueira da Foz, Portugal, for different sea states. This port is 2,250 m long and 400 m wide, approximately. Its average depth is of the order of 7 m, with approximately 12 m throughout the outer port basin.

Two simulations were performed, considering the state of rest as initial condition in both cases, although for different tidal heights. The first case corresponds to the entrance of a sinusoidal wave in the open sea (input boundary), with a mean wave height H = 3.90 m, period T = 15 s, wavelength $$\lambda$$ = 147 m, and direction $$\mathrm{\Phi }$$ = 260° W. The second case corresponds to a mean wave height H = 4.80 m, period T = 17.5 s, wavelength $$\lambda$$ = 173 m, and direction $$\mathrm{\Phi }$$ = 258.2° W. Figures 6 and 7 show perspective views of the free surface computed in the basin 269.1 s and 360 s after excitation. As can be seen, a zone with stronger agitation is observed in the outer harbor in the second case.

Figure 6.

Port of Figueira da Foz. Perspective view of the free surface computed 269.1 s after excitation, for a mean wave height, H = 3.90 m, period T = 15 s, wavelength $$\lambda$$ = 147 m and direction $$\mathrm{\Phi }$$ = 283°W. Tide height, 3.35 m (HZ).

Figure 7.

Port of Figueira da Foz. Perspective view of the free surface computed 360 s after excitation, for a mean wave height H = 4.80 m, period T = 17.5 s, wavelength $$\lambda$$ = 173 m and direction $$\mathrm{\Phi }$$ = 258.2° W. Tide height, 2.65 m (HZ).

4.3. 1HD Serre’s extension model

4.3.1. Numerical solution

The equation system (75) is solved using an efficient finite–difference method, whose consistency and stability were tested in [1] and [2] by comparison with a closed-form solitary wave solution of the Serre equations. For this purpose, the terms containing derivatives in time of u are grouped. The final system of three equations is rewritten according to the following equivalent form (SERIMP model) [1, 2 3]:

$$\begin{array}{l}{h}_t+{\left(uh\right)}_x=0\text{a}\\ q_t+{\left\{uq-\frac{1}{2}\left[u^2+\left(1+2\alpha \right)h^2{\left(u_x\right)}^2+\left(1+2\alpha \right){\left({\xi }_x\right)}^2{u}^2-h{\xi }_x{\left(u^2\right)}_x\right]\right\}}_x\\ +\left[g\left(1+\alpha \mathrm{\Omega }\right)+\alpha hu{u}_{xx}\right]{\eta }_x-\alpha gh{h}_x{\eta }_{xx}-\beta g\frac{h^2}{3}{\eta }_{xxx}\text{b}\\ -\frac{\alpha }{2}\left[{\left(h{\xi }_{xx}{u}^2\right)}_x+h_x{\xi }_{xx}{u}^2-h{\xi }_{xx}u{u}_x\right]+\left(\alpha -\frac{\beta }{3}\right)h^2{u}_x{u}_{xx}\\ +\beta \frac{h^2}{3}u{u}_{xxx}+{\tau }_b/\left(\rho h\right)=0\\ \left[1+\left(1+\alpha \right)\mathrm{\Omega }\right]u-\left(1+\alpha \right)h{h}_x{u}_x-\left(1+\beta \right)\frac{h^2}{3}{u}_{xx}=q\text{c}\\ \mathrm{\Omega }={\xi }_x{\eta }_x+\frac{1}{2}h{\xi }_{xx}+{\left({\xi }_x\right)}^2\text{d}\end{array}$$E136

To compute the solution of equation system (75) (values of the variables h and u at time $$t+\mathrm{\Delta }t$$), we use a numerical procedure based on the following scheme, itself based on the last equation system (136), for variables h, q and u. Knowing all values of $$h_i$$ and $$u_i$$, $$i=1,N$$, in the whole domain at time $$n\mathrm{\Delta }t$$, the equations (136c) and (136d) are used to obtain the first values of $$q_i$$ and $${\mathrm{\Omega }}_i$$ in the whole domain. Then, we continue with the following steps, in which the index p means predicted values (see also [1, 2] and [3]):

1. The first equation (136a) is used to predict the variable values $$h_{pi}$$ at time $$t+\mathrm{\Delta }t$$ ($$h_{pi}^{t+\mathrm{\Delta }t}$$), in the whole domain.

2. The second equation (136b) makes it possible to predict the variable values $$q_{pi}$$ at time $$t+\mathrm{\Delta }t$$ ($$q_{pi}^{t+\mathrm{\Delta }t}$$), taking into account the values $${\tilde{h}}_i^{t+\mathrm{\Delta }t}=0.5\left(h_i^t+h_{pi}^{t+\mathrm{\Delta }t}\right)$$, namely for $${\mathrm{\Omega }}_i$$ in the whole domain.

3. The third equation (136c) makes it possible to compute the mean-averaged velocities $$u_i^{t+\mathrm{\Delta }t}$$ at time $$t+\mathrm{\Delta }t$$, taking into account the predicted values $$h_{pi}^{t+\mathrm{\Delta }t}$$ and $$q_{pi}^{t+\mathrm{\Delta }t}$$, namely for $${\mathrm{\Omega }}_i$$ in the whole domain.

4. The first operation (step 1) is repeated in order to improve the accuracy of the variable values $$h_i$$ at time $$t+\mathrm{\Delta }t$$ ($$h_i^{t+\mathrm{\Delta }t}$$), using the values $${\overline{u}}_i^{t+\mathrm{\Delta }t}=0.5\left(u_i^t+u_i^{t+\mathrm{\Delta }t}\right)$$ in the whole domain.

5. Finally, the second operation (step 2) is repeated in order to improve the accuracy of the variable values $$q_i$$ at time $$t+\mathrm{\Delta }t$$ ($$q_i^{t+\mathrm{\Delta }t}$$), taking into account the values $${\overline{h}}_i^{t+\mathrm{\Delta }t}=0.5\left(h_i^t+h_i^{t+\mathrm{\Delta }t}\right)$$ and $${\overline{u}}_i^{t+\mathrm{\Delta }t}=0.5\left(u_i^t+u_i^{t+\mathrm{\Delta }t}\right)$$ in the whole domain.

At each interior point i, the first, second and third-order spatial derivatives are approximated through centered differences and the time derivatives are approximated using forward differences. The convective terms $${\left(uh\right)}_x$$ and $${\left(uq\right)}_x$$ in equations (136a) and 136b) are approximated through centered schemes in space and time for variables h and q. At each time t, these terms are written in the following form:

$${\left(uh\right)}_x=u_i^t\left(\frac{h_{i+1}^t-h_{i-1}^t+h_{i+1}^{t+\mathrm{\Delta }t}-h_{i-1}^{t+\mathrm{\Delta }t}}{4\mathrm{\Delta }x}\right)+\frac{1}{2}\left(h_i^t+h_i^{t+\mathrm{\Delta }t}\right)\left(\frac{u_{i+1}^t-u_{i-1}^t}{2\mathrm{\Delta }x}\right)$$E137

$${\left(uq\right)}_x=u_i^t\left(\frac{q_{i+1}^t-q_{i-1}^t+q_{i+1}^{t+\mathrm{\Delta }t}-q_{i-1}^{t+\mathrm{\Delta }t}}{4\mathrm{\Delta }x}\right)+\frac{1}{2}\left(q_i^t+q_i^{t+\mathrm{\Delta }t}\right)\left(\frac{u_{i+1}^t-u_{i-1}^t}{2\mathrm{\Delta }x}\right)$$E138

All finite–difference equations are implicit. Therefore, the solution of system (136) requires, in each time step, the computation of five three-diagonal systems of N-2 equations (steps 1 to 5), which are easily computed using the three-diagonal matrix algorithm (TDMA), also known as the Thomas algorithm. The stability condition to be observed can be expressed in terms of the Courant/CFL number, and is given by:

$$C_R=\sqrt{gh}\frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}<1.0$$E139

4.3.2. Boundary conditions

We often prescribe an influx on the left boundary, usually a mono- or bi-chromatic wave flow. The initial condition for this (mono-chromatic) influx is:

$$\eta \left(0,t\right)=\frac{a}{2}\{\sin \left(\omega t\right)\tanh \left(t-\frac{\pi }{\omega }\right)+\sin \left(\omega t\right)-\frac{\cos \left(\omega t\right)}{\omega }\left[1-{\tanh }^2\left(t-\frac{\pi }{\omega }\right)\right]\}$$E140

$$u\left(0,t\right)=\frac{a}{2}\frac{g}{\omega }\cos \left(\omega t\right)\left[\tanh \left(t-\frac{\pi }{\omega }\right)+1\right]$$E141

where a is the wave amplitude and the tanh is used as a ramp function. It increases smoothly from 0 to 1, to create a smooth start. For t large, the boundary condition reduces to:

$$\eta \left(0,t\right)=\frac{a}{2}\sin \left(\omega t\right)$$E142

$$u\left(0,t\right)=\frac{g}{\omega }\cos \left(\omega t\right)$$E143

If we want to avoid reflections at the right boundary (output), the domain is extended with a damping region of length L_damp. In this case, terms like $$-m\left(x\right)\eta$$ and $$-m\left(x\right)u$$ may be added to the continuity (136a) and momentum (136b) equations, respectively. The length of the damped region is chosen such that we do not see any significant reflections.

4.3.3. Solitary wave travelling up a slope and reflection on a vertical wall

Experimental data and numerical results are available for a solitary wave propagating on the bathymetry shown in Figure 8 [1, 2]. It shows a constant depth before x = 55 m and a slope 1:50 between x = 55 m and x = 75 m. An impermeable vertical wall is placed at x = 75 m, corresponding to fully reflecting boundary conditions. A solitary wave of amplitude 0.12 m is initially centered at x = 25 m. The computational domain was uniformly discretized with a spatial step $$\mathrm{\Delta }x=0.05$$ m. A zero friction factor has been considered. Computations were carried out with a time step $$\mathrm{\Delta }t=0.010$$ s. Figure 9 compares numerical time series of surface elevation and test data at x = 72.75 m.

Figure 8.

Bathymetry for a solitary wave travelling up a slope and its reflection on a vertical wall (not in scale).

Figure 9.

Solitary wave travelling up a slope and its reflection on a vertical wall. Free surface elevation in a depth gauge located at x = 72.75 m. Experimental ($$―$$); Serre extended ($$\cdots$$); Boussinesq extended (^╼╼) [1, 2, 3].

Figure 9 shows two peaks; the first one corresponding to the incident wave, and the second to the reflected wave. The extended Serre model predictions for both peaks agree well with the measurements. RMSE values equal to 0.0090 m and 0.0117 m were found in first and second peaks, respectively, for the wave height. Regarding the phase, there is a loss of approximately 0.05 s and of 0.10 s in those peaks (for details see [1] and [2]).

Predictions of the extended Boussinesq equations for both peaks are less accurate. Particularly for the reflected peak, this is overestimated in about 20%. This result is not surprising, given the lower validity of the Boussinesq model for waves of higher relative amplitude. Indeed, this model assumes $$O\left(\epsilon \right)<<1$$, contrary to the Serre model, which is $$O\left(\epsilon \right)=1$$. A visual comparison of numerical results of the extended Boussinesq approximation with a similar study performed by [37], using the extended Boussinesq model developed by [30], shows no relevant differences in the graphs (see [1, 2] and [3]).

4.3.4. Periodic wave over an underwater bar

Beji and Battjes (1993) [10] conducted experiments in a flume of 0.80 m wide with a submerged trapezoidal bar. The up- and down-wave bottom slopes of the submerged bar are 1:20 and 1:10, respectively. Before and after the bar, the water depth is 0.40 m, with a reduction to 0.10 m above the bar, as shown in Figure 10. Experimental data obtained in this installation are available in the literature, and can be used for comparisons.

Figure 10.

Bathymetry for a periodic wave propagating over a bar (not in scale) [2].

The measured data are compared with numerical results of the extended Boussinesq model (57), with $$\alpha =0.1308$$ and $$\gamma =-0.0076$$, and the extended Serre equations (75) (SERIMP model, in [2] and [3]), both improved with linear dispersive characteristics. Comparisons are made in three wave gauges located at x = 10.5 m, x = 13.5 m and x = 17.3 m. For this purpose, a regular incident wave case with height 0.02 m, period T = 2.02 s and wavelength 3.73 m has been simulated. The computational domain was discretized with a uniform grid interval $$\mathrm{\Delta }x=0.025$$ m. A time step $$\mathrm{\Delta }t=0.0010$$ s was used. Globally, numerical results of the improved Serre and Boussinesq models agree quite well with the measured data (for details see [2] and [3]).

Following is presented a comparison of the standard Serre’s model (20) with the extended Serre equations (75) (SERIMP model). The standard Serre’s model (20) is only valid for shallow waters, thus under conditions up to $$h_0/\lambda =0.05$$. In this experiment, the dispersion parameter ($$\sigma =h_0/\lambda$$) is greater than 0.05 (about 0.11) in front and behind the bar, and therefore affects the validity of the numerical outcomes. Due to the fact that over the bar there are very shallow water conditions ($$\sigma \approx 0.03$$) the standard Serre equations are used considering the input boundary located at section x = 13.5 m, where the input signal is known (measured data). In this way, results of the Serre’s standard model are not influenced, as would happen, by changes arising from the wave propagation before the bar, under intermediate water depths.

Figure 11 shows a comparison of numerical results of the 1HD standard Serre’s model (20) with the extended Serre equations (75), considering, in the first case, the input boundary at x = 13.5 m (gauge signal) (see [2]). The influence of additional terms of dispersive origin included in the extended Serre equations is clearly shown in Figure 11. The standard Serre model results (dashed line) are clearly of lesser quality. It should be noted that this application also demonstrates the good behavior of our numerical model to propagate a complex signal imposed at boundary.

Figure 11.

Periodic wave propagating over a bar. Comparison of test data ($$―$$) with numerical results of the extended Serre model (75) ($$\cdots$$) and the standard Serre equations (20) (^{_ _ _}) [2].

5. Conclusions and recommendations

This work presents some of the most recent methodologies to improve the linear dispersion characteristics of the classical Boussinesq equations for variable depth. A simple procedure to improve the linear dispersive performance of Serre’s equations is presented as well.

Extensions of the classical Boussinesq and Serre equations for intermediate depths and deep-water applications are really important, and have been the subject of major developments in recent decades, since these are the conditions typically found in nearshore areas. Application examples to solve real-world problems clearly demonstrate the capabilities and potential of the Boussinesq- and Serre-type models with improved linear dispersion characteristics.

The influence of the dispersion characteristics is clearly evidenced by the generation and propagation of waves in intermediate water depths, as shown here in cases of very demanding applications over bottoms with considerable slopes. The overall agreement of the extended Serre model with improved dispersion characteristics is very good both in shallow water conditions as in intermediate water depths. The influence of additional terms of dispersive origin included in the extended Serre’s equations is clearly shown.

This work also shows that current models of Boussinesq type with improved dispersive characteristics are accurate enough for applications in extensive real-world conditions.

Although the contribution of this work to improve the dispersive characteristics of the 1HD standard Serre equations (20) is clear, the author is aware of the need to further improve the dispersion properties of the equation system (75).

An extended version of the two-dimensional system of equations (19) with improved linear dispersion characteristics is being developed and will be published soon.