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

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 expen‐ sive and more flexible than the construction of physical models, they are a conven‐ ient 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 h 0 / λ , other than shallow waters, where h 0 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 solu‐ tion of a Boussinesq-type model, also improved with dispersive characteristics. Ex‐ plicit and implicit methods of finite–difference are implemented to solve both approximations of Boussinesq and Serre types with improved dispersive character‐ istics. A finite element method is also implemented to solve an extension of a Bous‐ sinesq-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 ap‐ plications, 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.


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 (σ = h 0 / λ, where h 0 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.

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 [ 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 = − h 0 + ξ(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 / h 0 and σ = h 0 / λ, in which a is a characteristic wave amplitude, h 0 represents water depth, and λ is a characteristic wavelength, we proceed with suitable nondimensional variables:

x y' y z' z h ' a ' h t' t gh tc p' p gh u' u a g h uh ac v' v a g h vh ac w' w ah g h w ac
where c 0 represents critical celerity, given by c 0 = (gh 0 ) 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]: a. Fundamental equations ( ) ( ) ( ) where the bar over the variables represents the average value along the vertical.Then, accepting the fundamental hypothesis of the shallow water theory, σ = h 0 / λ < < 1, and devel- oping the dependent variables in power series of the small parameter σ 2 , that is , for where A = u x + v y , 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: On the other hand, taking into account that, ( ) ( ) from ( 5) and ( 9) we obtain: Representing by Γ = w t + εuw x + εvw y + εww z the vertical acceleration of the particles, we get Γ = w 0t + εu 0 w 0x + εv 0 w 0 y + εw 0 w 0z + O ( σ 2 ) , and from ( 6), ( 7) and ( 11) the following approach is obtained: 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]: or even, given that ( f s ) ** = f s ** − ε( f z ) ** η s , where f = (u, v) and s = (x, y, t): h es e e h es By developing expressions (17) in second approach (order 2 in σ 2 ), the following equations of motion (18) are obtained (for details see [34]): New Perspectives in Fluid Dynamics 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 = h 0 − ξ + η 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 / h 0 ) having a value close to the square of the relative depth (σ = h 0 / λ), 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 h r = h 0 − ξ is the water column height at rest, P and Q are given by P = − (h 0 − ξ)(u ¯x + v ¯y) t and 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: ( ) ( )

Derivation of higher-order equations
3.1.Weakly nonlinear approaches with improved dispersive performance

Nwogu's approach
To allow applications in a greater range of h 0 / λ, other than shallow waters, [27] introduced higher-order dispersive terms into the governing equations to improve linear dispersion New Perspectives in Fluid Dynamics 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 | z by w, with ξ t = 0 and using Leibnitz' rule, the boundary condition 3(a) at z = − 1 + ξ gives: Substituting (28) in the irrotationality equation 2(e), u z = σ 2 w x , 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 Differentiating equation (33) with respect to z and noting that both u z and u zz are O ( σ 2 ) : 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): 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: Using the velocities u (36) and w (37) in the vertical momentum, equation 2(d) yields: This can be rearranged to give an expression for the pressure p: Integrating through the depth from z to the free surface εη: Using the free surface boundary condition 3(c) for the pressure, and denoting p | z simply as p gives: ( ) ( ) 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: 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: and thus 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 α = αh 0 , where − 1 ≤ α ≤ 0, the system (51)-( 52) is rewritten: In two dimensions, the equation system ( 53)-( 54) is written: 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 twodimensional velocity field at depth z = αh .

Beji and Nadaoka's, and Liu and Sun's approaches
Starting from the standard Boussinesq equations and adopting the methodology introduced by [27] [11] presented a new approach that allows for applications until values of h 0 / λ to the order of 0.25, and still with acceptable errors in amplitude and phase velocity up to values of h 0 / λ near 0.50.By an addition and subtraction process, using the approximation u t = − 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 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: New Perspectives in Fluid Dynamics 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(kh r ), 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".

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.
where, as above, z is the vertical coordinate starting at the still water level h 0 (x, y) and pointing upwards, scaled by a typical depth h 0 , and η is the water surface displacement scaled by a representative amplitude a.The two dimensionless parameters ε and μ 2 are defined as ε = a / h 0 and μ 2 = (k 0 h 0 ) 2 , with a representative wave number k 0 = 2π / λ, so μ 2 = (2πσ) 2 .Time t is scaled by k 0 (gh 0 ) 1/2 −1 , and φ, the velocity potential, is scaled by εh 0 (gh 0 ) 1/2 .We use the nondimensional 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]: New Perspectives in Fluid Dynamics 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: 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.3and 4.1.4.
To improve the accuracy of numerical models, it has been common practice the use of highorder 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".

Mathematical model for one-layer
Defining the parameters S 1 = ∇ .u 1 and , the model uses the following approach for the continuity equation (to compute η values), in nondimensional variables: where u 1 = horizontal velocity vector, k 1 = α 1 h + β 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: 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 z 1 = − 0.531h , and that value has been recommended and adopted by most researchers who use these equations.Refs.[22] and later [38] used the Nwogu's approach to derive a high-order model which retains "full nonlinearity" by including terms not considered in the Nwogu's approach.Comparing Liu's and Wei and Kirby's equations, there are some differences that can be attributed to the absence of some nonlinear dispersive terms in Wei and Kirby [19].The above one-layer equations ( 69) and (70) are identical to those derived by [22].

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

1HD Serre's equations with improved dispersive performance
To allow applications in a greater range of h 0 / λ, other than shallow waters, a new set of extended Serre equations, with additional terms of dispersive origin, is developed and tested in [1] and [2] by comparisons with the available test data.The equations are solved using an efficient finite-difference method, whose consistency and stability are tested in that work by comparison with a closed-form solitary wave solution of the Serre equations.
From the equation system (20), by adding and subtracting terms of dispersive origin, using the approximation u t = − g η x and 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(kh r ), 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.

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 + u c ) and V = (v + v c ), where index c represents current and u = U -u c and v = Vv c , the final set of these equations are written as follows: ( ) ( ) 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: 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(h p t +Δt ), considering the known values of h, U and V at time t in the whole domain.80) and (81) make it possible to predict the values of variables r (r p t +Δt ) and s (s p t +Δt ), taking into account the values of U t , V t , r t , s t and h ˜t+Δt = 0.5h t + 0.5h p t +Δt known for the whole domain.82) and (83) give us the values of the mean-averaged velocity components U and V (U t +Δt and V t +Δt ), taking into account the predicted values of r and s (r p t +Δt and s p t +Δt , respectively).

Equation (79) allows us to compute the depth h at time t + Δt (values of h t +Δt ) considering
the values of variables h t , U t +05Δt = 0.5U t + 0.5U t +Δt and V t +05Δt = 0.5V t + 0.5V t +Δ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 r t +Δt and s t +Δt ), taking into account the values of r t , s t h t +05Δt = 0.5h t + 0.5h t +Δt , U t +05Δt = 0.5U t + 0.5U t +Δt and V t +05Δt = 0.5V t + 0.5V t +Δt known for the whole domain.
The Petrov-Galerkin procedure is utilized to achieve solutions for the unknowns h, r and s.
According to the weighted residual technique, minimization requires the "orthogonality" of the residual R J to a set of weighting functions W i,J , i.e., , 0, 1, , where the general form of the weighting functions applied to these equations is defined as: , and where the coefficients δ u i , and δ v i are functions of: (i) the local velocities U and V; (ii) the ratio of the wave amplitude to the water depth, and (iii) the element length.To illustrate this procedure, a complete solution of equations (79) (80) and ( 82) is presented here in detail.
Introducing in equation ( 79), the approximated values are given by: In the matrix form, this equation may be written as: The residual R 80 is written as:

R r u r v r uU vU
It should be noted that the term 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 = h r 2 / 3, all terms involving ∇ f in the numerical developments are omitted.
The residual R 82 is written as: According to Galerkin's procedure, after using integration by parts (or Green's theorem) to reduce the second derivatives, the R 82 error minimization leads to the following equation: The last term of (102) can be written as: New Perspectives in Fluid Dynamics where p, q = 1, ⋯ , n e , n e being the number of nodes in the corresponding element side coincident with the boundary domain, Γ i e represents the element sides within the domain, with the corresponding integral null because the resulting element contributions are equal, but with opposite signals, and Γ e e represents the element sides coincident with the boundary domain.
Accordingly, an equivalent form of equation ( 102) may be written as follows: In the matrix form, this equation may be written as: where 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.
a. Element side lower than the local depth.

b. Minimum of 20 to 25 elements per wave length.
c. Courant number always lower than one in the whole domain".

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).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 quasistationary 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.

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 finitedifference 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 ∇ η ( ∇ .(hr u α ) t + h r tt / ε ) and ∇ ( η 2 / 2 ) ∇ .uαt , can be reformulated using the relations: ( The explicit predictor equations are: ( ) All first-order spatial derivatives are differenced with fourth-order (Δx 4 = Δy 4 ) accurate equations, which are five-point differences.Second-order spatial derivatives are approximated with three-point centered finite-difference equations, which are second-order accurate.The fourth-order implicit corrector expressions for the free surface elevation, η, and the horizontal velocities, u and 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 η, 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" 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 10 −6 .Linear stability analysis performed by [38,18] and [40] show that Δt < Δx / (2 c) to ensure stability, where c is the celerity.

Research study using COULWAVE model
For the analysis concerning coastal protection, the mean currents around a submerged structure (artificial reef) are analysed.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).
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).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.

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: Considering the generic variable F , the solution at time (n + 1)Δt, for the computational point (i, j), is obtained from the known solution F i, j n through the following symmetric application: 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]).

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]: New Perspectives in Fluid Dynamics 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 h i and u i , i = 1, N , in the whole domain at time nΔt, the equations (136c) and (136d) are used to obtain the first values of q i 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 h pi at time t + Δt (h pi t +Δt ), in the whole domain.

2.
The second equation (136b) makes it possible to predict the variable values q pi at time t + Δt (q pi t +Δt ), taking into account the values h ˜i t +Δt = 0.

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

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 L damp .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.

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

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.
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  Following is presented a comparison of the standard Serre's model (20) with the extended Serre equations (75) (SERIMP model).The standard Serre's model ( 20) is only valid for shallow waters, thus under conditions up to h 0 / λ = 0.05.In this experiment, the dispersion parameter (σ = h 0 / λ) 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.

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 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.

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

Figure 2 .
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 .
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.
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. Figures6 and 7show 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 .
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 .
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).

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

Figure 10 .
Figure 10.Bathymetry for a periodic wave propagating over a bar (not in scale) [2]. ) [16]u'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: [39]et al. (1995)[39]used the Nwogu's approach to derive a Boussinesqtype model which retains full nonlinearity by including O ( εσ 2 ) terms not considered in the .

Table 1 .
The output of the COULWAVE model corresponds to the velocity values at a depth 0.531 h r below the water surface.The velocity at this depth is used to determine the velocity cells near the shoreline that could give an indication of the sediment transport.Divergent cells indicate erosion near by the shoreline and convergent cells indicate sedimentation.The numerical simulations to study the 2HD behavior of the hydrodynamics around the reef have been done for four cases: two reef geometries, varying the reef angle (45 o and 66 o ), and two wave conditions.The characteristics of the simulations for the different cases are described in Table1.Main characteristics of the simulations performed.
[7]v, α = d 2 h / dt 2 ; β = d 2 ξ / dt 2 and the bottom friction terms, τ bx and τ by , are Ox and Oy directions.The corresponding operators Lx and Ly take the following form[7]: of the classical Boussinesq and Serre equations for intermediate depths and deepwater 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.