Open access peer-reviewed chapter

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

Written By

José Simão Antunes do Carmo

Reviewed: October 28th, 2015 Published: December 2nd, 2015

DOI: 10.5772/61866

Chapter metrics overview

1,483 Chapter Downloads

View Full Metrics


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.


  • Extended Boussinesq equations
  • Extended Serre equations
  • Dispersive characteristics
  • Intermediate waters
  • Numerical methods

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 σ2 (σ=h0/λ, where h0 and λ 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 σ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ρ/dt=div v=0), irrotationality (rot v=0) and perfect fluid (dynamic viscosity, μ=0)]:


with p=0 at z=η(x,y,t), w=ηt+uηx+vηy at z=η(x,y,t), and w=ξt+uξx+vξy at z=h0+ξ(x,y,t). In these equations ρ is density, t is time, g is gravitational acceleration, p is pressure, η is free surface elevation, ξ is bottom, and u, v, w are velocity components. Defining the dimensionless quantities ε=a/h0 and σ=h0/λ, in which a is a characteristic wave amplitude, h0 represents water depth, and λ is a characteristic wavelength, we proceed with suitable nondimensional variables:


where c0 represents critical celerity, given by c0=(gh0)1/2, and, as above,η is free surface elevation, ξ 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

  1. Boundary conditions


Integrating the first equation 2(a) between the bed 1+ξ and the free surface εη, taking into account 3(a) and 3(b), yields the continuity equation (4):

[η(1/ε) ξ]t+[(1ξ+εη) u¯]x+[(1ξ+εη) v¯]y=0E4

where the bar over the variables represents the average value along the vertical. Then, accepting the fundamental hypothesis of the shallow water theory, σ=h0/λ<<1, and developing the dependent variables in power series of the small parameter σ2, that is


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


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]:


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

u¯=u0+σ2u1*(σ2/6)(1+εη0ξ0)2A0x          +(σ2/2)(1+εη0ξ0)(ξ0xA0+w0x*)+O(σ4)v¯=v0+σ2v1*(σ2/6)(1+εη0ξ0)2A0y          +(σ2/2)(1+εη0ξ0)(ξ0yA0+w0y*)+O(σ4)E10

On the other hand, taking into account that,


from (5) and (9) we obtain:

u**=u¯(σ2/3)(1+εηξ)2A¯x          +(σ2/2)(1+εηξ)(ξxA¯+wx*)+O(σ4)v**=v¯(σ2/3)(1+εηξ)2A¯y          +(σ2/2)(1+εηξ)(ξyA¯+wy*)+O(σ4)E12

Representing by Γ=wt+εuwx+εvwy+εwwz the vertical acceleration of the particles, we get Γ=w0t+εu0w0x+εv0w0y+εw0w0z+O(σ2), and from (6), (7) and (11) the following approach is obtained:

Γ=(z+1ξ)[A¯t+εu¯A¯x+εv¯A¯yεA¯2]              +[wt*+εu¯wx*+εv¯wy*]+O(σ2)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:


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


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

(ut+ εuux+ εvuy+ εwuz)**+ ηx(1+εσ2Γ**)=0(vt+ εuvx+ εvvy +εwvz)**+ ηy(1+εσ2Γ**)=0E16

or even, given that (fs)**=fs**ε(fz)**ηs, where f=(u,v) and s=(x,y,t):

ut**+ εu**ux**+ εv**uy**+ ηx(1+εσ2Γ**)=0vt**+ εu**vx**+ εv**vy**+ ηy(1+εσ2Γ**)=0E17

By developing expressions (17) in second approach (order 2 in σ2), the following equations of motion (18) are obtained (for details see [34]):

u¯t+εu¯u¯x+εv¯u¯y+ ηx +σ2{ [(2/3)(εηξ)x+(1/2) ξx] P+(1/3)(1+εηξ)Px} +σ2[εηxQ+(1/2)(1+εηξ) Qx]+σ4=0v¯t+εu¯v¯x+εv¯v¯y+ ηy +σ2{ [(2/3)(εηξ)y+(1/2) ξy] P+(1/3)(1+εηξ)Py} +σ2[εηyQ+(1/2)(1+εηξ) Qy]+σ4=0P=(1+εηξ)(εA¯2εu¯A¯xεv¯A¯yA¯t)Q=wt+ε u¯wx+εv¯wyw=(1/ε) ξt+u¯ξx+v¯ξyA¯=u¯x+v¯yE18

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

ht+(hu¯)x+(hv¯)y=0u¯t+u¯u¯x+v¯u¯y+gηx +[(2/3)hx+(1/2) ξx]P+(1/3)hPx+hxQ+(1/2) hQx=0v¯t+u¯v¯x+v¯v¯y+gηy +[(2/3)hy+(1/2) ξy]P+(1/3)hPy+hyQ+(1/2) hQy=0P=h (A¯2u¯A¯xv¯A¯yA¯t)Q=wt+u¯wx+v¯wyw=u¯ξx+v¯ξyA¯=u¯x+v¯yE19

where h=h0ξ+η is total water depth. The one-dimensional form (1HD) of the equation system (19) is written, also with a fixed bottom:

ht+(u¯h)x=0hu¯t+hu¯u¯x+gh ηx+[h2(P/3+Q/2)]x+ξxh(P/2+Q)=0P=h (u¯xt+u¯u¯xxu¯x2)Q=ξx(u¯t+u¯u¯x)+ξxxu¯2E20

Assuming additionally a relative elevation of the surface due to the waves (ε=a/h0) having a value close to the square of the relative depth (σ=h0/λ), i.e. O(ε) = O(σ2), from the system of equations (18), and at the same order of approximation, the following approach is obtained, in dimensional variables:

ht+(hu¯)x+(hv¯)y=0u¯t+u¯u¯x+v¯u¯y+gηx[(1/6) ξx]P+(1/3)hrPx+(1/2) hrQx=0v¯t+u¯v¯x+v¯v¯y+gηy[(1/6)ξy]P+(1/3)hrPy+(1/2) hrQy=0E21

where hr=h0ξ is the water column height at rest, P and Q are given by P=(h0ξ) (u¯x+v¯y)t and Q=(u¯ξx+v¯ξy)t. The momentum equations are written as:

u¯t+u¯u¯x+v¯u¯y+gηx+(1/6) hrξx (u¯x+v¯y)t(1/3)hr2(u¯x+v¯y)xt +(1/3)hrξx(u¯x+v¯y)t+(1/2)hr(u¯ξx+v¯ξy)xt=0E22
v¯t+u¯v¯x+v¯v¯y+gηy+(1/6) hrξy (u¯x+v¯y)t(1/3)hr2(u¯x+v¯y)yt +(1/3)hrξy(u¯x+v¯y)t +(1/2)hr(u¯ξx+v¯ξy)yt=0E23

with ξt=0, the complete system of equations (24) is obtained:

ht+(hu¯)x+(hv¯)y=0u¯t+u¯u¯x+v¯u¯y+gηx(1/3) hr2(u¯xxt+v¯xyt)+hrξxu¯xt +(1/2) hr(ξxxu¯t+ξyv¯xt+ξxv¯yt+ξxyv¯t)=0v¯t+u¯v¯x+v¯v¯y+gηy(1/3)hr2(u¯xyt+v¯yyt)+hrξyv¯yt +(1/2) hr(ξyyv¯t+ξyu¯xt+ξxu¯yt+ξxyu¯t)=0E24

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


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 h0/λ, 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:

w|zw|1+ξ=(1+ξzu dz)x+u|1+ξξxu|zzxE27

Denoting w|z by w, with ξt=0 and using Leibnitz’ rule, the boundary condition 3(a) at z=1+ξ gives:

w=(1+ξzu dz)xE28

Substituting (28) in the irrotationality equation 2(e), uz=σ2wx, yields:

uz=σ2wx=σ2(1+ξzu dz)xxE29

Considering a Taylor series expansion (30) of u(x,z,t), about z=zα, where the horizontal velocity uα (x,t) denotes the velocity at depth zα, u (x,zα,t), this Taylor expansion is integrated through the depth from 1+ξ to z, yielding (31) (for details see [36]):

u=uα+(zzα) uz|z=zα+(zzα)22uzz|z=zα+E30
1+ξzu dz=(z+1ξ) uα+[(zzα)22(1ξ+zα)22] uz|z=zα +[(zzα)36+(1ξ+zα)36] uzz|z=zα+E31

Substituting (31) in equation (29) gives:

uz=σ2{(z+1ξ) uα+[(zzα)22(1ξ+zα)22] uz|z=zα+}xxE32

Differentiating equation (32) with respect to z, noting from (29) that uz=O(σ2):

uzz=σ2{uα+(zzα) uz|z=zα+(zzα)22uzz|z=zα+}xx =σ2uαxx+O(σ4)E33

Differentiating equation (33) with respect to z and noting that both uz and uzz are O(σ2):

uzzz=σ2{ uz|z=zα+(zzα) uzz|z=zα+}xx=O(σ4)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(σ4) order or greater. Substituting equations (33) and (34) back in equation (32):

uz=σ2{(1ξ+zα) uα(1ξ+zα)22uz|z=zα+O(σ2)}xx =σ2[(1ξ+zα) uα]xx+O(σ4)E35

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

u=uασ2{(zzα) [(1ξ+zα) uα]xx+(zzα)22(uα)xx}+O(σ4)E36

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

w={(1ξ+zα) uα       +σ2[(1ξ+zα)22[(1ξ+zα) uα]xx(1ξ+zα)36(uα)xx]}x       +O(σ4)E37

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

εσ2{[(z+1ξ) uα]xt+O(σ2)} +pz+1+O(ε2σ2)=0E38

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

pz=1εσ2[(z+1ξ) uα]xt+O(εσ4,ε2σ2)E39

Integrating through the depth from z to the free surface εη:

p|εη+p|z=εηz+εσ2{[z22+(1ξ) z] u}xt+O(εσ4,ε2σ2)E40

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

p=εηz+εσ2{[z22+(1ξ) z] u}xt+O(εσ4,ε2σ2)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:

uαtσ2{(zzα)[(1ξ+zα) uα]xx+(zzα)22(uα)xx}t         +εuαuαx+ηx+σ2{[z22+(1ξ)z]uα}xxt=O(εσ2,σ4)E42


uαt+εuαuαx+ηx+σ2{zα[(1ξ) uα]xx+zα22(uα)xx}t=O(εσ2,σ4)E43

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

(1+ξεηu dz)xεu|εηηx+u|1+ξξx+w|εηw|1+ξ=0E44

Using Liebnitz’ Rule and the kinematic boundary conditions at the bed z=1+ξ 3(a) and at the free surface z=εη 3(b) gives:

ηt+(1+ξεηu dz)x=0E45

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

1+ξεηu dz=(1ξ+εη) uασ2{[(zzα)22]1+ξεη[(1ξ+zα) uα]xx +[(zzα)36]1+ξεηuαxx}+O(σ4)E46


1+ξεηu dz=(1ξ+εη) uασ2{[zα22(1ξ+zα)22][(1ξ+zα) uα]xx             +[zα36+(1ξ+zα)36] uαxx} +O(εσ2,σ4)E47


1+ξεηu dz=(1ξ+εη) uασ2{[(1ξ)22+(1ξ) zα][(1ξ+zα) uα]xx +[(1ξ)36+(1ξ)2zα2+(1ξ) zα22](uα)xx}+O(εσ2,σ4)E48

and thus

1+ξεηu dz=(1ξ+εη) uα+σ2{[zα+1ξ2](1ξ)[(1ξ) uα]xx                      +[zα22(1ξ)26] (1ξ)(uα)xx }+O(εσ2,σ4)E49

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

ηt+[(1ξ+εη) uα]x+σ2{(zα+1ξ2)(1ξ)[(1ξ) uα]xx +(zα22(1ξ)26)(1ξ)(uα)xx}x=O(εσ2,σ4)E50

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

uαt+uαuαx+gηx+zα(hr uαt)xx+zα22(uα)xxt=0E51
ηt+[(h0+η) uα]x+[(zα+hr2)hr(hr uα)xx+(zα22hr26) hr(uα)xx]x=0E52

Setting the arbitrary depth zα=αh0, where 1α0, the system (51)–(52) is rewritten:

uαt+uαuαx+gηx+αhr [hr(uα)t]xx+α22hr2(uα)xxt=0E53
ηt+[(hr+η) uα]x+[(α+12)hr2(hr uα)xx+(α2216) hr3(uα)xx]x=0E54

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

ηt+.[(hr+η) uα]     +{(α+12)hr2[.(hr uα)]+(α2216) hr3(.uα)}=0E55
uαt+(uα.)uα+gη+αhr[.(hr uαt)]+α22hr2 (.uαt)=0E56

where represents the two-dimensional gradient operator with respect to the horizontal coordinates (x, y) (=/x,/y), and the velocity vector u(x,y,t)=(u,v) represents the two-dimensional velocity field at depth z=α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 h0/λ to the order of 0.25, and still with acceptable errors in amplitude and phase velocity up to values of h0/λ near 0.50. By an addition and subtraction process, using the approximation ut=g ηx and considering a dispersion parameter β in the momentum equation of the 1HD system (24), Beji and Nadaoka obtained an improved set of Boussinesq equations for variable depth, with β 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 ω2/gk=tanh(khr). In order to improve dispersion and linear shoaling characteristics in the Beji and Nadaoka’s equations, [23] introduced two tuning parameters, α and γ, so that β=1.5α0.5γ. The nonlinearity in the previous Boussinesq-type models was improved by Liu and Sun adding higher-order terms accurate to the order of O(εσ2). 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:

ht+(hu)x=0ut+ uux+gηx {(1+α)/2(1+γ)/6}hr2 uxxt (α/2γ/6) g hr2 ηxxx +(1+α) hr ξxuxt +α g hr ξxηxx +0.5 (1+α) hr ξxxut+0.5 α g hr ξxxηx=O(εσ2)E57

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

ω2gk=khr[1+(α/2γ/6 ) k2hr2]1+(1+α) k2hr2/2(1+γ) k2hr2/6E58

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

C2=ω2k2=ghr[1+(α/2γ/6 ) (khr)21+(1+α) (khr)2/2(1+γ) (khr)2/6]E59

with the linear dispersion relation ω2/gk=tanh(khr), using the approach (60)


allows to obtain the best values for the parameters α and γ: α=0.1308 and γ=0.0076.

Considering appropriate values for the tuning parameters α and γ, we can identify within the equation system (57):

  • The standard Boussinesq equations by setting α=γ=0

  • The Beji and Nadaoka equations considering α=γ=0.20

  • The Liu and Sun equations with α=0.1308 and γ=0.0076

A visual comparison of numerical results of the extended Boussinesq approximation (57), with α=0.1308 and γ=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 α=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(εσ2) terms not considered in the Nwogu’s (53)–(54) system, and thus improving the nonlinearity to O (ε)=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:

φzz+μ22φ=0 ,                            1+ξzεηφz+μ2h1ξ.φ=0 ,                        z=1+ξη+φt+12ε[(φ)2+1μ2(φz)2]=0 ,    z=εηηt+εφ.η1μ2φz=0 ,                  z=εηE61

where, as above, z is the vertical coordinate starting at the still water level h0(x,y) and pointing upwards, scaled by a typical depth h0, and η is the water surface displacement scaled by a representative amplitude a. The two dimensionless parameters ε and μ2 are defined as ε=a/h0 and μ2=(k0h0)2, with a representative wave number k0=2π/λ, so μ2=(2πσ)2. Time t is scaled by [k0(gh0)1/2]1, and φ, the velocity potential, is scaled by εh0(gh0)1/2. We use the nondimensional water level h1ξ instead of 1ξ. Integrating the first equation of (61) over the water column, and using the appropriate boundary conditions, the continuity equation is obtained:


where M=1+ξεηφdz. Retaining terms to O(μ2), and denoting φα as the value of φ at z=zα(x,y), an approximate velocity potential is given by:


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

ηt+.{(h1ξ+εη){φα+μ2{[zα.(h1ξφα)+zα222φα]   +(h1ξεη)2[.(h1ξφα)][h1ξ2εηh1ξ+(εη)2]62φα}}}=0E64

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α, uα=φα, a fully nonlinear version of a Boussinesq type model in terms of η and uα is written:

ηt+.{(h1ξ+εη){uα+μ2{[zα2216(h1ξ2εηh1ξ+(εη)2)](.uα)                              +[zα+12(h1ξεη)][.(h1ξuα)]}}}=O(μ4)E65


S={(zαεη)(uα.)[.(h1ξuα)]+12[zα2(εη)2](uα.)(.uα)}     +12{[.(h1ξuα)+εη.uα]2}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 S1=.u1 and T1=.[(1ξ) u1]+(1/ε)(h1ξ/t), the model uses the following approach for the continuity equation (to compute η values), in nondimensional variables:

1ε h1ξt+ηt+.[(εη+h1ξ) u1]                 μ2.{[ε3η3+h1ξ36(εη+h1ξ) k122]S1                 +[ε2η2h1ξ22(εη+h1ξ) k1]T1}=O(μ4)E69

where u1 = horizontal velocity vector, k1=α1h+β1η, α1 and β1 are coefficients to be defined by the user. The index 1 means one-layer model. To compute the velocity components (u,v), the following approach of the momentum equation is solved, in nondimensional variables:

u1t+ε u1.u1+η+μ2t{k122S1+k1T1}       +εμ2[(u1.k1)T1+k1(u1.T1)+k1(u1.k1)S1       +k122(u1.S1)]+εμ2[T1T1(ηT1t)]       +ε2μ2(η S1T1η22S1tη u1.T1)       +ε2μ2[η22(S12u1.S1)]=O(μ4)E70

The horizontal velocity vector is given as:

U1=u1μ2{z12k122S1+(z1k1) T1}+O(μ4)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(ε)=O(μ2), 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 z1=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 η and the velocity components (u,v). In dimensional variables, this set of equations is written ([24]):

ηt+.[(ηζ) u1+(ζ+hr) u2]      .{[ζ3+hr36(ζ+hr)k222]S2+[ζ2hr22(ζ+hr)k2]T2}           .{[η3ζ36(ηζ)k122]S1+[η2ζ22(ηζ)k1]T1}=0E72
u1t+12(u1.u1)+gη       +t[k122S1+k1T1(η22S1)(ηT1)] +[ηt(T1+ηS1)+(k1η)(u1.) T1+12(k12η2)(u1.)S1 +12(T1+ηS1)2]Rb+Rf+νT[S12u1      2(k122S1+k1T1)+(η222S1+η2T1)]=0E73
u2+k22ζ22S2+(k2ζ) T2=u1+k12ζ22S1+(k1ζ) T1E74

where S1=.u1, T1=ζ(S2S1)+T2, S2=.u2, and T2=.(hru2). Rb = breaking-related dissipation term, Rf=[f/(hr+η)]ub|ub| accounts for bottom friction, where ub = velocity evaluated at the seafloor, and f = bottom friction coefficient, typically in the range of 103102, νT = constant eddy viscosity, 2=(2/x2,2/y2), k1=0.127hr = evaluation level for the velocity u1, ζ=0.266hr = layer interface elevation, s = evaluation level for the velocity u2, and η = free surface elevation.

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

To allow applications in a greater range of h0/λ, 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 ut=g ηx and considering the parameters α, β and γ, with β=1.5α0.5γ, allows to obtain a new system of equations with improved linear dispersion characteristics:

ht+(uh)x=0ut+uux+g (h+ξ)x+(1+α)(Ω uthhxuxt)(1+β)h23uxxt +α g Ω(h+ξ)xα ghhx(h+ξ)xxβ gh23(h+ξ)xxxhhxuuxx +h23(uxuxxuuxxx)+h (ux)2(h+ξ)x+ξxxu2(h+ξ)x + (Ω +h ξxx) uux+h2ξxxxu2=0E75

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 ω2/gk=tanh(khr), using the approach (60), values of α=0.1308 and γ=0.0076 are obtained, so that β=0.20. It should be noted that the Serre’s equation system (20) is recovered by setting α=β=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=(u+uc) and V=(v+vc), where index c represents current and u=U-uc and v=V-vc, the final set of these equations are written as follows:

Ut+UUx+VUy+g(h+ξ)xhr23(Uxxt+Vxyt)      hr23[uc(Uxx+Vxy)+vc(Uxy+Vyy)]x+hr(12ξtt+ucξxt+vcξyt)x       νl(Uxx+Uyy) τsxρh+τbxρh=0E77
Vt+UVx+VVy+g(h+ξ)yhr23(Uxyt+Vyyt)      hr23[uc(Uxx+Vxy)+vc(Uxy+Vyy)]y+hr(12ξtt+ucξxt+vcξyt)y       νl(Vxx+Vyy)τsyρh+τbyρh=0E78

where τs and τ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:

rt+ucrx+vcry=uUxvUyg(h+ξ)x      +hr23[(uc)x(Uxx+Uyy)+(vc)x(Vxx+Vyy) ] hr(12ξtt+ucξxt+vcξyt)x       +[νl+23hr(ξt+ucξx+vcξy)](Uxx+Uyy)+τsxρhτbxρh=0E80
st+ucsx+vcsy=uVxvVyg(h+ξ)y      +hr23[(uc)y(Uxx+Uyy)+(vc)y(Vxx+Vyy) ]hr(12ξtt+ucξxt+vcξyt)y       +[νl+23hr(ξt+ucξx+vcξy)](Vxx+Vyy) +τsyρhτbyρh=0E81

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+Δt (for details see [5]):

  1. The equation (79) allows us to predict the values of variable h(hpt+Δ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 (rpt+Δt) and s (spt+Δt), taking into account the values of Ut, Vt, rt, st and h˜t+Δt=0.5ht+0.5hpt+Δ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 (Ut+Δt and Vt+Δt), taking into account the predicted values of r and s (rpt+Δt and spt+Δt, respectively).

  4. Equation (79) allows us to compute the depth h at time t+Δt (values of ht+Δt) considering the values of variables ht, Ut+05Δt=0.5Ut+0.5Ut+Δt and Vt+05Δt=0.5Vt+0.5Vt+Δ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+Δt (values of rt+Δt and st+Δt), taking into account the values of rt, st ht+05Δt=0.5ht+0.5ht+Δt, Ut+05Δt=0.5Ut+0.5Ut+Δt and Vt+05Δt=0.5Vt+0.5Vt+Δ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 RJ to a set of weighting functions Wi,J, i.e.,


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


and where the coefficients δui, and δvi 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:


the following residual R79 is obtained:


The R79 error minimization leads to the following equation:

ΔeWi,79R79dΔe=ΔeWi,79(h^t+h^U^x+U^h^x+h^V^y+V^h^y) dΔe          =ΔeWi,79[j=1nNj(ht)j+k=1n(Nx)kUkj=1nNjhj+k=1nNkUkj=1n(Nx)jhj                            +k=1n(Ny)kVkj=1nNjhj +k=1nNkVkj=1n(Ny)jhj] dΔe=0E88

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



ai,j=ΔeWi,79NjdΔe;i,j=1,,nbi,j=ΔeWi,79k=1nNkUk(Nx)jdΔe+ΔeWi,79k=1n(Nx)kUkNjdΔe        +ΔeWi,79k=1n(Ny)kVkNjdΔe+ΔeWi,79k=1nNkVk(Ny)jdΔe;  i,j=1,,nE90

The residual R80 is written as:

R80=r^t+u^cr^x+v^cr^y+u^U^x+v^U^y       +g(h^+ξ^)x h^r23[(u^c)x(U^xx+U^yy)+(v^c)x(V^xx+V^yy)]       +h^r[12(ξ^tt)x+(u^c)x(ξ^t)x+u^c(ξ^t)xx+(v^c)x(ξ^t)y+v^c(ξ^t)xy]       νl(U^xx+U^yy)τsxρh+τbxρhE91

It should be noted that the term 23hr(ξt+ucξx+vcξy)(Vxx+Vyy) of equation (80) is of order σ4 or greater. For this reason, it is not considerd in the numerical developments. Similarly, and for the same reason, considering f=hr2/3, all terms involving f in the numerical developments are omitted.

The Green’s theorem is used to solve the second derivatives present in equation (80) (residual R80), and in equations (81) to (83), i.e., considering p^=(U^,V^), from p^xx+p^yy we obtain:

ΔeWi,l(p^xx+p^yy) dΔe=ΔeWi,lj=1n[(Nxx)j+(Nyy)j] pj dΔeE92

and so

ΔeWi,lj=1n[(Nxx)j+(Nyy)j] pj dΔe=Δe(Wx)i,lj=1n(Nx)jpj dΔe             Δe(Wy)i,lj=1n(Ny)jpj dΔe+ΓeWi,lj=1n(Nn)jpj dΓe ;  p=(U,V)E93

Retaining terms up to order σ3, the R80 error minimization leads to:

ΔeWi,80R80dΔe =Δe{Wi,80[j=1nNj(rt)j+k=1nNk(uc)kj=1n(Nx)jrj         +k=1nNk(vc)kxj=1n(Ny)jrj +k=1nNkukj=1n(Nx)jUj+k=1nNkvkj=1n(Ny)jUj        +gj=1n(Nx)j(h+ξ)j]l=1nNl(hr2)l3k=1n(Nx)k[(uc)kUx2y2+(vc)kVx2y2]         +Wi,80[k=1nNk(hr)k2j=1n(Nx)j(ξtt)j+l=1nNl(hr)lk=1n(Nx)k(uc)kj=1n(Nx)j(ξt)j                        +l=1nNl(hr)lk=1n(Nx)k(vc)kj=1n(Ny)j(ξt)j]         +l=1nNl(hr)l[k=1nNk(uc)k(ξ˜t)xx+k=1nNk(vc)k(ξ˜t)xy]νlUx2y2          Wi,80k=1nNk[(τsx)k(τbx)k]/j=1hnNjhj}dΔe=0E94





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



ci=Δe{Wi,80[k=1nNkukj=1n(Nx)jUj+k=1nNkvkj=1n(Ny)jUj        +gj=1n(Nx)j(h+ξ)j] l=1nNl(hr2)l3k=1n(Nx)k[(uc)kUx2y2+(vc)kVx2y2]        +Wi,80[k=1nNk(hr)k2j=1n(Nx)j(ξtt)j+l=1nNl(hr)lk=1n(Nx)k(uc)kj=1n(Nx)j(ξt)j                        +l=1nNl(hr)lk=1n(Nx)k(vc)kj=1n(Ny)j(ξt)j]        +l=1nNl(hr)l[k=1nNk(uc)k(ξ˜t)xx+k=1nNk(vc)k(ξ˜t)xy] νlUx2y2        Wi,80k=1nNk[(τsx)k(τbx)k]/j=1hnNjhj}dΔe=0 ,  i=1,.nE100

The residual R82 is written as:


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

ΔeWi,82R82dΔe=Δe{Wi,82j=1nNjUj                +k=1nNk(hr2)k3[(Wx)i,82j=1n(Nx)j+(Wy)i,82j=1n(Ny)j]Uj                Wi,82j=1nNjrj}dΔe ΓeNk(hr2)k3Wi,82(Nn)jUj dΓeE102

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

ΓeNk(hr2)k3Wi,82(Nn)jUj dΓe=ΓieNp(hr2)p3Wq,82(Nn)rUr dΓie                                                      +ΓeeNp(hr2)p3Wq,82Un dΓeeE103

where p,q=1,,ne, ne being the number of nodes in the corresponding element side coincident with the boundary domain, Γie represents the element sides within the domain, with the corresponding integral null because the resulting element contributions are equal, but with opposite signals, and Γee represents the element sides coincident with the boundary domain. Accordingly, an equivalent form of equation (102) may be written as follows:

ΔeWi,82R82dΔe=Δe{Wi,82j=1nNj            +k=1nNkh(r2)k3[(Wx)i,82j=1n(Nx)j+(Wy)i,82j=1n(Ny)j]}Uj dΔe            =ΔeWi,82j=1nNjrjdΔe+ΓeNp(hr2)p3Wq,82Un dΓeep,q=1,,neE104

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



ai,j=Δe{Wi,82j=1nNj+k=1nNk(hr2)k3[(Wx)i,82j=1n(Nx)j+(Wy)i,82j=1n(Ny)j]} dΔeE106
bi=ΔeWi,82j=1nNjrjdΔe+ΓeNp(hr2)p3Wq,82Un dΓep,q=1,,ne.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+k2η22uxx+(kη)(hru)xxηx[η ux+(hru)x]E108
V=v+k2η22vyy+(kη)(hrv)yyηy[η vy+(hrv)y]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 [η(.(hruα)t+hrtt/ε)] and [(η2/2).uαt], can be reformulated using the relations:

[η(.(hruα)t+hr ttε)]=[η(. (hruα)+hr tε)]t[ηt(. (hruα)+hr tε)]E110
(η22.uαt)=(η22. uα)t(ηηt. uα)E111

The explicit predictor equations are:

Ui,jn+1=Ui,jn+Δt12(23Fi,jn16Fi,jn1+5Fi,jn2)+2 (F1)i,jn3 (F1)i,jn1+ (F1)i,jn2E113
Vi,jn+1=Vi,jn+Δt12(23Gi,jn16Gi,jn1+5Gi,jn2)+2 (G1)i,jn3 (G1)i,jn1+ (G1)i,jn2E114


E=hrt[(η+hr) u]x[(η+hr) v]y +{(η+hr) [(16(η2ηhr+hr2)12k2)Sx+(12(ηhr)k) Tx]}x +{(η+hr) [(16(η2ηhr+hr2)12k2)Sy+(12(ηhr)k) Ty]}yE115
F=12[(u2)x+(v2)x]gηxkhrxttkthrxt+(Ehrt+η hrtt)x    [E (ηS+T)]x[12(k2η2) (uSx+vSy)]x    [(kη) (u Tx+v Ty)]x12[(T+η S)2]xE116
F1=η2k22vxy(kη) (hrv)xy+ηx[η vy+(hrv)y]E117
G=12[(u2)y+(v2)y]gηykhryttkthryt+(Ehrt+η hrtt)y    [E (η S+T)]y[12(k2η2) (uSx+vSy)]y    [(kη) (u Tx+v Ty)]y12[(T+η S)2]yE118
G1=η2k22uxy(kη) (hru)xy+ηy[η ux+(hru)x]E119



All first-order spatial derivatives are differenced with fourth-order (Δx4=Δy4) 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, η, and the horizontal velocities, u and v, are:

Ui,jn+1=Ui,jn+Δt24(9Fi,jn+1+19Fi,jn5Fi,jn1+Fi,jn2)+ (F1)i,jn+1 (F1)i,jnE122
Vi,jn+1=Vi,jn+Δt24(9Gi,jn+1+19Gi,jn5Gi,jn1+Gi,jn2)+ (G1)i,jn+1 (G1)i,jnE123

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 η, 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 |wn+1w*n+1wn+1|<ε100or|wn+1w*n+1||wn+1|<εE124

In these expressions, w represents any of the variables η, u and v, and w is the previous iterations value. The value of the error is set to 106. Linear stability analysis performed by [38, 18] and [40] show that Δt<Δx/(2 c) 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 hr 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 (45o and 66o), and two wave conditions. The characteristics of the simulations for the different cases are described in Table 1.

Reef angle
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 Δx = Δ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:

Pt+(uP)x+(vP)y+{ [(g+β)/2+α/3]h2}x =(g+β+α/2) hξxτbx+R div (h grad u)E126
Qt+(uQ)x+(vQ)y+{ [(g+β)/2+α/3]h2}y =(g+β+α/2) hξyτby+R div (h grad v)E127

where P=hu, Q=hv, α=d2h/dt2; β=d2ξ/dt2 and the bottom friction terms, τbx and τby, are obtained through (128):


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

Pt+(uP)x+{ [(g+β)/2+α/3]h2}x=(g+β+α/2) hξxτbx+R PxxE130
Qt+(uQ)x=R QxxE131

Operator Ly

Qt+(vQ)y+{ [(g+β)/2+α/3]h2}y=(g+β+α/2) hξyτby+R QyyE134

Considering the generic variable F, the solution at time (n+1)Δt, for the computational point (i,j), is obtained from the known solution Fi,jn through the following symmetric application:

Fi,jn+1=Lx(Δt4)Ly(Δt4)Lx(Δt4)Ly(Δt4)Ly(Δt4)Lx(Δt4)Ly(Δt4)Lx(Δt4) Fi,jnE135

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 (u,v) are updated and the values of the vertical accelerations, α and β, 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 λ = 147 m, and direction Φ = 260° W. The second case corresponds to a mean wave height H = 4.80 m, period T = 17.5 s, wavelength λ = 173 m, and direction Φ = 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 λ = 147 m and direction Φ = 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 λ = 173 m and direction Φ = 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]:

ht+(uh)x=0aqt+{uq12[u2+(1+2α) h2(ux)2+(1+2α)(ξx)2u2h ξx(u2)x]}x +[g (1+αΩ)+αh u uxx] ηxαghhxηxxβgh23ηxxxb α2[(h ξxxu2)x+hx ξxxu2h ξxxu ux]+(αβ3) h2uxuxx       +βh23u uxxx+τb/(ρh)=0[1+(1+α) Ω] u(1+α) h hxux(1+β)h23uxx=qcΩ=ξxηx+12h  ξxx+(ξx)2dE136

To compute the solution of equation system (75) (values of the variables h and u at time t+Δ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 hi and ui, i=1,N, in the whole domain at time nΔt, the equations (136c) and (136d) are used to obtain the first values of qi and Ω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 hpi at time t+Δt (hpit+Δt), in the whole domain.

  2. The second equation (136b) makes it possible to predict the variable values qpi at time t+Δt (qpit+Δt), taking into account the values h˜it+Δt=0.5(hit+ hpit+Δt), namely for Ωi in the whole domain.

  3. The third equation (136c) makes it possible to compute the mean-averaged velocities uit+Δt at time t+Δt, taking into account the predicted values hpit+Δt and qpit+Δt, namely for Ωi in the whole domain.

  4. The first operation (step 1) is repeated in order to improve the accuracy of the variable values hi at time t+Δt (hit+Δt), using the values u¯it+Δt=0.5 (uit+ uit+Δt) in the whole domain.

  5. Finally, the second operation (step 2) is repeated in order to improve the accuracy of the variable values qi at time t+Δt (qit+Δt), taking into account the values h¯it+Δt=0.5(hit+ hit+Δt) and u¯it+Δt=0.5 (uit+ uit+Δt) 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 (uh)x and (uq)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:

(uh)x=uit(hi+1thi1t+hi+1t+Δthi1t+Δt4Δx)+12(hit+hit+Δt) (ui+1tui1t2Δx)E137
(uq)x=uit(qi+1tqi1t+qi+1t+Δtqi1t+Δt4Δx)+12(qit+qit+Δt) (ui+1tui1t2Δx)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:


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:


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:


If we want to avoid reflections at the right boundary (output), the domain is extended with a damping region of length Ldamp. In this case, terms like m(x) η and m(x) 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 Δx=0.05 m. A zero friction factor has been considered. Computations were carried out with a time step Δ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 (); 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(ε)<<1, contrary to the Serre model, which is O(ε)=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 α=0.1308 and γ=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 Δx=0.025 m. A time step Δ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 h0/λ=0.05. In this experiment, the dispersion parameter (σ=h0/λ) 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 (σ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) () 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.


  1. 1. Antunes do Carmo JS. Boussinesq and Serre type models with improved linear dispersion characteristics: applications. J Hydraulic Res 2013a;51(6):719–27:814090. DOI: 10.1080/00221686.814090.
  2. 2. Antunes do Carmo JS. Extended Serre equations for applications in intermediate water depths. Open Ocean Eng J 2013b;6:16–25.
  3. 3. Antunes do Carmo JS. Applications of Serre and Boussinesq type models with improved linear dispersion characteristics. In: Jesús MB, Irene A, Alberto P, José MG, Nuno S, Miguel S. (Eds.) Congreso de Métodos Numéricos en Ingeniería; 25–28 June; Bilbao, Spain. Barcelona, Spain: International Center for Numerical Methods in Engineering (CIMNE); 2013c. DOI: ISBN 978-84-941531-4-3.
  4. 4. Antunes do Carmo JS. Closure to “Boussinesq- and Serre-type models with improved linear dispersion characteristics: applications”. J Hydraulic Res 2015;53:284–5.
  5. 5. Antunes do Carmo JS, Seabra-Santos FJ. On breaking waves and wave-current interaction on shallow water: a 2DH finite element model. Int J Num Method Fluids 1996;22:429–44.
  6. 6. Antunes do Carmo JS, Seabra-Santos FJ. Wave-current interactions over bottom with appreciable variations in both space and time. Adv Eng Software 2010;41:295–305.
  7. 7. Antunes do Carmo JS, Seabra-Santos FJ, Almeida AB. Numerical solution of the generalized Serre equations with the MacCormack finite-difference scheme. Int J Num Method Fluids 1993b;16:725–38. DOI: 10.1002/fld.1650160805.
  8. 8. Antunes do Carmo JS, Seabra-Santos FJ, Barthélemy E. Surface waves propagation in shallow-water: a finite element model. Int J Num Method Fluids 1993a;16:447–59. DOI: 10.1002/fld.1650160602.
  9. 9. Antunes do Carmo JS, Voorde MT, Neves MG. Enhancing submerged coastal constructions by incorporating multifunctional purposes. J Coastal Conserv 2011;15(4):531–46.
  10. 10. Beji S, Battjes JA. Experimental investigations of wave propagation over a bar. Coastal Eng 1993;19(1,2):151–62.
  11. 11. Beji S, Nadaoka K. A formal derivation and numerical modelling of the improved Boussinesq equations for varying depth. Ocean Eng 1996;23(8):691–704.
  12. 12. Berkhoff JCW, Booij N, Radder AC. Verification of numerical wave propagation models for simple harmonic linear water waves. Coastal Eng 1982;6:255–79.
  13. 13. Booij N. A note on the accuracy of the mild-slope equation. Coastal Eng 1983;7:191–203.
  14. 14. Boussinesq J. Théorie des ondes et des remous qui se propagent le long d’un canal rectangulaire horizontal. J Math Pures Appliq 1872;2(17):55–108.
  15. 15. Dalrymple RA. Model for refraction of water waves. J Waterway Port CoastalOcean Eng 1988;114(4):423–35.
  16. 16. Gobbi MFG, Kirby JT, Wei G. A fully nonlinear Boussinesq model for surface waves. Part 2. Extension to O(kh)4. J Fluid Mechanics 2000;405:181–210.
  17. 17. Green AE, Naghdi PM. A derivation of equations for wave propagation in water of variable depth. J Fluid Mechanics 1976;78(2):237–46.
  18. 18. Hsiao SC. Permeable effects on nonlinear water waves [thesis]. Cornell University:2000.
  19. 19. Hsiao SC, Liu PLF. Permeable effects on nonlinear water waves. In: Royal Society of London, editor. London. London:2002.
  20. 20. Kirby JT. A note on linear surface wave-current interaction over slowly varying topography. J Geophys Res 1984;89(C):745–47.
  21. 21. Kirby JT, Dalrymple RA. A parabolic equation for the combined refraction-diffraction of stokes waves by mildly varying topography. J Fluid Mechanics 1983;136:435–66.
  22. 22. Liu PLF. Model equations for wave propagation from deep to shallow water. In: Liu PLF, (Ed.) Advances in Coastal Engineering. World Scientific. 1994, pp. 125–157.
  23. 23. Liu ZB, Sun ZC. Two sets of higher-order Boussinesq-type equations for water waves. Ocean Eng 2005;32:1296–310.
  24. 24. Lynett P. Nearshore wave modeling with high-order Boussinesq-type equations. J Waterway Port Coastal Ocean Eng 2006;132(5):348–57.
  25. 25. Lynett P, Liu PLF. (Eds.) Modeling Wave Generation, Evolution, and Interaction with Depth Integrated, Dispersive Wave Equations COULWAVE Code Manual. 1st ed. Ithaca, NY: Cornell University Long and Intermediate Wave Modeling Package; 2002.
  26. 26. Lynett P, Liu PLF. Linear analysis of the multi-layer model. Coastal Eng 2004;51:439–54. DOI: 10.1016/coastaleng.2004.05.004.
  27. 27. Madsen PA, Murray R, Sørensen OR. A new form of the Boussinesq equations with improved linear dispersion characteristics. Coastal Eng 1991;15(4):371–88.
  28. 28. Mendonça AM, Fortes CJ, Capitão R, Neves MG, Moura T, Antunes do Carmo JS. Wave hydrodynamics around a multi-functional artificial reef at Leirosa. J Coastal Conserv 2012a;16(4):543–53. DOI: 10.1007/s11852-012-0196-1.
  29. 29. Mendonça AM, Fortes CJ, Capitão R, Neves MG, Antunes do Carmo JS, Moura T. Hydrodynamics around an Artificial Surfing Reef at Leirosa, Portugal. J Waterway Port Coastal Ocean Eng 2012b;138(3):226–35. DOI: 10.1061/(ASCE)WW.1943-5460.0000128.
  30. 30. Nwogu O. Alternative form of Boussinesq equations for nearshore wave propagation. J Waterway Port Coastal Ocean Eng 1993;119:618–38.
  31. 31. Press WH, Flannery BP, Teukolsky SA. Numerical Recipes. Cambridge University Press. 1989, pp. 569–72.
  32. 32. Saint-Venant B. Theory of unsteady water flow, with application to river floods and to propagation of tides in river channels. Computes Rendus Acad. Sci., Paris. 1871;73:148–54, 237–40.
  33. 33. Seabra Santos FJ. Contribution a l’ètude des ondes de gravité bidimensionnelles en eau peu profonde [thesis]. Université Scientifique et Médicale et Institut National Polutechnique de Grenoble, France:1985.
  34. 34. Seabra Santos FJ. As aproximações de Wu e de Green & Naghdi no quadro geral da teoria da água pouco profunda (in Portuguese). In: APRH, editor. Simpósio Luso-Brasileiro de Hidráulica e Recursos Hídricos (4º SILUSBA); 14–16 June; Lisbon, Portugal, LNEC, Lisbon:1989. pp. 209–219.
  35. 35. Serre F. Contribution à l’étude des écoulements permanents et variables dans les canaux. La Houille Blanche 1953;8:374–88.
  36. 36. Walkley MA. A Numerical Method for Extended Boussinesq Shallow-Water Wave Equations [thesis]. Leeds, UK: University of Leeds, School of Computer Studies:1999.
  37. 37. Walkley MA, Berzins M. A finite element method for the one-dimensional extended Boussinesq equations. Int J Num Method Fluids 1999;29:143–57.
  38. 38. Wei G, Kirby JT. Time-dependent numerical code for extended Boussinesq equations. J Waterway Port Coastal Ocean Eng 1995;121(5):251–61.
  39. 39. Wei G, Kirby JT, Grilli ST, Subramanya R. A fully nonlinear Boussinesq model for surface waves. I. Highly nonlinear, unsteady waves. J Fluid Mechanics 1995;294:71–92.
  40. 40. Woo S-B. Finite element modeling of the fully-nonlinear extended Boussinesq equations [thesis]. Ithaca, NY: Cornell University:2002.

Written By

José Simão Antunes do Carmo

Reviewed: October 28th, 2015 Published: December 2nd, 2015