Open access peer-reviewed chapter

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

By José Simão Antunes do Carmo

Submitted: October 30th 2014Reviewed: October 28th 2015Published: December 2nd 2015

DOI: 10.5772/61866

Downloaded: 1138


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 h0and λ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=divv=0), irrotationality (rotv=0) and perfect fluid (dynamic viscosity, μ=0)]:


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


where c0represents critical celerity, given by c0=(gh0)1/2, and, as above,ηis free surface elevation, ξrepresents bathymetry, u, vand ware velocity components, and pis 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):


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+εwwzthe 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 pon 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]):


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:


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


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:


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


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


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:


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


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


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


Substituting (31) in equation (29) gives:


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


Differentiating equation (33) with respect to zand noting that both uzand uzzare O(σ2):


Repeated differentiation of this expression will produce expressions for the higher derivatives of uwith respect to z and show them to be of O(σ4)order or greater. Substituting equations (33) and (34) back in equation (32):


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


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:


Integrating through the depth from zto the free surface εη:


Using the free surface boundary condition 3(c) for the pressure, and denoting p|zsimply as pgives:


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



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


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


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



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



and thus

1+ξεηudz=(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:


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


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

ηt+[(hr+η) uα]x+[(α+12)hr2(hruα)xx+(α2216)hr3(uα)xx]x=0E54

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

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

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 ηxand 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:


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


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


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


allows to obtain the best values for the parameters αand γ: α=0.1308and γ=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.1308and γ=0.0076

A visual comparison of numerical results of the extended Boussinesq approximation (57), with α=0.1308and γ=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, zis 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 μ2are defined as ε=a/h0and μ2=(k0h0)2, with a representative wave number k0=2π/λ, so μ2=(2πσ)2. Time tis 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=.u1and 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η, α1and β1are 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:


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

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 ηxand considering the parameters α, βand γ, with β=1.5α0.5γ, allows to obtain a new system of equations with improved linear dispersion characteristics:


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.1308and γ=0.0076are 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 crepresents current and u=U-ucand 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 τsand τbrepresent 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, rand sare 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, Uand Vat time tin 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, stand h˜t+Δt=0.5ht+0.5hpt+Δtknown for the whole domain.

  3. Solutions of equations (82) and (83) give us the values of the mean-averaged velocity components Uand V(Ut+Δtand Vt+Δt), taking into account the predicted values of rand s(rpt+Δtand spt+Δt, respectively).

  4. Equation (79) allows us to compute the depth hat time t+Δt(values of ht+Δt) considering the values of variables ht, Ut+05Δt=0.5Ut+0.5Ut+Δtand Vt+05Δt=0.5Vt+0.5Vt+Δtknown for the whole domain.

  5. Equations (80) and (81) allow us to compute the values of variables rand sat time t+Δt(values of rt+Δtand st+Δt), taking into account the values of rt, stht+05Δt=0.5ht+0.5ht+Δt, Ut+05Δt=0.5Ut+0.5Ut+Δtand Vt+05Δt=0.5Vt+0.5Vt+Δtknown for the whole domain.

The Petrov-Galerkin procedure is utilized to achieve solutions for the unknowns h, rand s. According to the weighted residual technique, minimization requires the “orthogonality” of the residual RJto 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 δviare functions of: (i) the local velocities Uand 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 R79is obtained:


The R79error 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 R80is 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 σ4or 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 fin 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^yywe 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 R80error 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 R82is written as:


According to Galerkin’s procedure, after using integration by parts (or Green's theorem) to reduce the second derivatives, the R82error 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, nebeing the number of nodes in the corresponding element side coincident with the boundary domain, Γierepresents the element sides within the domain, with the corresponding integral null because the resulting element contributions are equal, but with opposite signals, and Γeerepresents 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]):


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

The explicit predictor equations are:



F=12[(u2)x+(v2)x]gηxkhrxttkthrxt+(Ehrt+ηhrtt)x  [E(ηS+T)]x[12(k2η2)(uSx+vSy)]x  [(kη)(uTx+vTy)]x12[(T+ηS)2]xE116
G=12[(u2)y+(v2)y]gηykhryttkthryt+(Ehrt+ηhrtt)y  [E(ηS+T)]y[12(k2η2)(uSx+vSy)]y  [(kη)(uTx+vTy)]y12[(T+ηS)2]yE118



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, uand v, are:


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η,uandv”. Next, the corrector expressions are evaluated, and again uand vare 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


In these expressions, wrepresents any of the variables η, uandv, and wis 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 cis 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 hrbelow 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 wavelengthGrid size (m)Time step (s)

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, uand 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. FromFigure 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:


where P=hu, Q=hv, α=d2h/dt2; β=d2ξ/dt2and the bottom friction terms, τbxand τ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 Oxand Oydirections. The corresponding operators Lxand Lytake the following form [7]:

Operator Lx

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

Operator Ly


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,jnthrough 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, Lxand Ly, is composed of a predictor–corrector sequence and nrepresents 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, periodT= 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 heightH= 4.80 m, periodT= 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 uare 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+αΩ)+αhuuxx]ηxαghhxηxxβgh23ηxxxbα2[(hξxxu2)x+hxξxxu2hξxxuux]+(αβ3)h2uxuxx       +βh23uuxxx+τ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 hand uat time t+Δt), we use a numerical procedure based on the following scheme, itself based on the last equation system (136), for variables h, qand u. Knowing all values of hiand 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 qiand Ωiin the whole domain. Then, we continue with the following steps, in which the index pmeans predicted values (see also [1, 2] and [3]):

  1. The first equation (136a) is used to predict the variable values hpiat time t+Δt(hpit+Δt), in the whole domain.

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

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

  4. The first operation (step 1) is repeated in order to improve the accuracy of the variable values hiat 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 qiat 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)xand (uq)xin 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 ais the wave amplitude and the tanhis used as a ramp function. It increases smoothly from 0 to 1, to create a smooth start. For tlarge, 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) umay 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.05m. A zero friction factor has been considered. Computations were carried out with a time step Δt=0.010s. 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 atx= 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. RMSEvalues 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.1308and γ=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.025m. A time step Δt=0.0010s 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 Serreequations (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.

© 2015 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

José Simão Antunes do Carmo (December 2nd 2015). Modeling of Wave Propagation from Arbitrary Depths to Shallow Waters – A Review, New Perspectives in Fluid Dynamics, Chaoqun Liu, IntechOpen, DOI: 10.5772/61866. Available from:

chapter statistics

1138total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

CFD-Based Investigation of Wind-Strokes over Highway Bridge Section

By Medzid Muhasilovic, Kenan Imsirpasic, Karel Ciahotny and Brano Sirok

Related Book

First chapter

Surface Friction and Boundary Layer Thickening in Transitional Flow

By Ping Lu, Manoj Thapa and Chaoqun Liu

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us