Main characteristics of the simulations performed.
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
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  Berkhoff et al., 1982,  Kirby and Dalrymple, 1983,  Booij, 1983,  Kirby, 1984, and  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 , even at that time models based on the Saint-Venant equations (see also  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 ). 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  “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 (, where 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 ) 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 (), irrotationality () and perfect fluid (dynamic viscosity, )]:
with at , at , and at . 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 and , in which is a characteristic wave amplitude, represents water depth, and is a characteristic wavelength, we proceed with suitable nondimensional variables:
where represents critical celerity, given by , and, as above,is free surface elevation, represents bathymetry, , and are velocity components, and is pressure. In dimensionless variables, without the line over the variables, the fundamental equations and the boundary conditions are written :
where the bar over the variables represents the average value along the vertical. Then, accepting the fundamental hypothesis of the shallow water theory, , and developing the dependent variables in power series of the small parameter , that is
where the simple and double asterisk represent the variables values at the bottom and at the surface, respectively. Of 2(e) we obtain, successively :
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,
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:
or even, given that , where and :
where, likewise, the bar over the variables represents the average value along the vertical. In dimensional variables and with a solid/fixed bottom (), the complete set of equations is written, in second approach:
where 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 () having a value close to the square of the relative depth (), i.e. = , from the system of equations (18), and at the same order of approximation, the following approach is obtained, in dimensional variables:
where is the water column height at rest, and are given by and . The momentum equations are written as:
with , 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)  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 , other than shallow waters,  introduced higher-order dispersive terms into the governing equations to improve linear dispersion properties. By redefining the dependent variable,  achieved the same improvement without the need to add such terms to the equations.
Following , the extended Bussinesq equations obtained by  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 by , with and using Leibnitz’ rule, the boundary condition 3(a) at gives:
Considering a Taylor series expansion (30) of , about , where the horizontal velocity denotes the velocity at depth , , this Taylor expansion is integrated through the depth from to , yielding (31) (for details see ):
Differentiating equation (33) with respect to z and noting that both and are :
Repeated differentiation of this expression will produce expressions for the higher derivatives of u with respect to z and show them to be of order or greater. Substituting equations (33) and (34) back in equation (32):
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 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 :
Using Liebnitz’ Rule and the kinematic boundary conditions at the bed 3(a) and at the free surface 3(b) gives:
From the expression (36) for the horizontal velocity u:
where represents the two-dimensional gradient operator with respect to the horizontal coordinates (x, y) (), and the velocity vector represents the two-dimensional velocity field at depth .
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   presented a new approach that allows for applications until values of to the order of 0.25, and still with acceptable errors in amplitude and phase velocity up to values of near 0.50. By an addition and subtraction process, using the approximation 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,  introduced two tuning parameters, and , so that . The nonlinearity in the previous Boussinesq-type models was improved by Liu and Sun adding higher-order terms accurate to the order of . 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:
with the linear dispersion relation , using the approach (60)
allows to obtain the best values for the parameters and : and .
Considering appropriate values for the tuning parameters and , we can identify within the equation system (57):
The standard Boussinesq equations by setting
The Beji and Nadaoka equations considering
The Liu and Sun equations with and
A visual comparison of numerical results of the extended Boussinesq approximation (57), with and ), shown in [1, 2] and , with a similar study performed by , using the extended Boussinesq (53)–(54) model (Nwogu’s approach, with ) shows no relevant differences in the graphs. In this regard, it is worth remembering  “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)  used the Nwogu’s approach to derive a Boussinesq-type model which retains full nonlinearity by including terms not considered in the Nwogu’s (53)–(54) system, and thus improving the nonlinearity to . Wei et al., and later  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:
where, as above, z is the vertical coordinate starting at the still water level and pointing upwards, scaled by a typical depth , and is the water surface displacement scaled by a representative amplitude a. The two dimensionless parameters and are defined as and , with a representative wave number , so . Time t is scaled by , and , the velocity potential, is scaled by . We use the nondimensional water level instead of . Integrating the first equation of (61) over the water column, and using the appropriate boundary conditions, the continuity equation is obtained:
where . Retaining terms to , and denoting as the value of at , an approximate velocity potential is given by:
Similarly, substituting (63) into the third equation of (61), a momentum equation is obtained in terms of the velocity potential. Then, given that at , , a fully nonlinear version of a Boussinesq type model in terms of and 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  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 , which consists of using quadratic polynomials, matched at interfaces that divide the water column into layers. In this regard, it is worth mentioning  "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 and , the model uses the following approach for the continuity equation (to compute values), in nondimensional variables:
where = horizontal velocity vector, , and are coefficients to be defined by the user. The index 1 means one-layer model. To compute the velocity components , 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. ), has been examined and applied to a significant extent. The weakly nonlinear version of (69)–(71) (i.e., assuming , thereby neglecting all nonlinear dispersive terms) was first derived by . Through a linear and first-order nonlinear analysis of the model equations, Nwogu recommended that , and that value has been recommended and adopted by most researchers who use these equations. Refs.  and later  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 . The above one-layer equations (69) and (70) are identical to those derived by .
3.2.2. Mathematical model for two-layers
Details about the two-layer mathematical model are presented in  and . 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 . In dimensional variables, this set of equations is written ():
where , , , and . = breaking-related dissipation term, accounts for bottom friction, where = velocity evaluated at the seafloor, and = bottom friction coefficient, typically in the range of , = constant eddy viscosity, , = evaluation level for the velocity , = layer interface elevation, s = evaluation level for the velocity , and = free surface elevation.
3.3. 1HD Serre’s equations with improved dispersive performance
To allow applications in a greater range of , other than shallow waters, a new set of extended Serre equations, with additional terms of dispersive origin, is developed and tested in  and  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 and considering the parameters , and , with , 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 , using the approach (60), values of and are obtained, so that . It should be noted that the Serre’s equation system (20) is recovered by setting .
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 . This model is named WACUP (a 2HD WAve Plus CUrrent Boussinesq-type model). With dimensional variables, taking mean quantities of the horizontal velocity components and , where index c represents current and and , the final set of these equations are written as follows:
where and represent stresses on surface and at bottom, respectively.
The standard Boussinesq model (24) and the extended system of equations (76)–(78) are solved in  and , 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 (for details see ):
The equation (79) allows us to predict the values of variable h(), considering the known values of h, U and V at time t in the whole domain.
Equation (79) allows us to compute the depth h at time (values of ) considering the values of variables , and 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 to a set of weighting functions , i.e.,
where the general form of the weighting functions applied to these equations is defined as:
and where the coefficients , and are functions of: (i) the local velocities U and V; (ii) the ratio of the wave amplitude to the water depth, and (iii) the element length. To illustrate this procedure, a complete solution of equations (79) (80) and (82) is presented here in detail. Introducing in equation (79), the approximated values are given by:
the following residual is obtained:
The error minimization leads to the following equation:
In the matrix form, this equation may be written as:
The residual is written as:
It should be noted that the term of equation (80) is of order or greater. For this reason, it is not considerd in the numerical developments. Similarly, and for the same reason, considering , all terms involving in the numerical developments are omitted.
Retaining terms up to order , the error minimization leads to:
In the matrix form, this equation may be written as:
The residual is written as:
According to Galerkin’s procedure, after using integration by parts (or Green's theorem) to reduce the second derivatives, the error minimization leads to the following equation:
The last term of (102) can be written as:
where , being the number of nodes in the corresponding element side coincident with the boundary domain, represents the element sides within the domain, with the corresponding integral null because the resulting element contributions are equal, but with opposite signals, and 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:
As recommended in  “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.
Element side lower than the local depth.
Minimum of 20 to 25 elements per wave length.
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).
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.
4.1.3. COULWAVE one-layer model
A numerical scheme similar to that of  and  is utilized by , with the inclusion of extra nonlinear terms. Basically using the same high-order predictor–corrector scheme,  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 , “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” . Also in accordance with , “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 ):
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 and , can be reformulated using the relations:
The explicit predictor equations are:
All first-order spatial derivatives are differenced with fourth-order () 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, and , are:
As noted in , “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 , “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 is the previous iterations value. The value of the error is set to . Linear stability analysis performed by [38, 18] and  show that to ensure stability, where c is the celerity.
4.1.4. Research study using COULWAVE model
For the analysis concerning coastal protection, the mean currents around a submerged structure (artificial reef) are analysed. The output of the COULWAVE model corresponds to the velocity values at a depth 0.531 hr below the water surface. The velocity at this depth is used to determine the velocity cells near the shoreline that could give an indication of the sediment transport. Divergent cells indicate erosion near by the shoreline and convergent cells indicate sedimentation.
The numerical simulations to study the 2HD behavior of the hydrodynamics around the reef have been done for four cases: two reef geometries, varying the reef angle (45o and 66o), and two wave conditions. The characteristics of the simulations for the different cases are described in Table 1.
|Number of grid points per wavelength||Grid size (m)||Time step (s)|
As recorded in , “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 = =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 ( 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 , 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 , a morphological study should be done in order to confirm these results.
4.2. 2HD Serre’s standard model
4.2.1. Numerical formulation
where , , ; and the bottom friction terms, and , are obtained through (128):
In order to apply the MacCormack’s method, equations (125) –(127) are split into two systems of three equations throughout the Ox and Oy directions. The corresponding operators and take the following form :
Considering the generic variable , the solution at time , for the computational point , is obtained from the known solution through the following symmetric application:
where each operator, and , is composed of a predictor–corrector sequence and represents a generic time . 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 , the values of the velocities are updated and the values of the vertical accelerations, and , are recalculated (for details see ).
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.
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  and  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]:
To compute the solution of equation system (75) (values of the variables h and u at time ), 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 and , , in the whole domain at time , the equations (136c) and (136d) are used to obtain the first values of and in the whole domain. Then, we continue with the following steps, in which the index p means predicted values (see also [1, 2] and ):
The first equation (136a) is used to predict the variable values at time (), in the whole domain.
The second equation (136b) makes it possible to predict the variable values at time (), taking into account the values , namely for in the whole domain.
The third equation (136c) makes it possible to compute the mean-averaged velocities at time , taking into account the predicted values and , namely for in the whole domain.
The first operation (step 1) is repeated in order to improve the accuracy of the variable values at time (), using the values in the whole domain.
Finally, the second operation (step 2) is repeated in order to improve the accuracy of the variable values at time (), taking into account the values and 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 and 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:
4.3.2. Boundary conditions
We often prescribe an influx on the left boundary, usually a mono- or bi-chromatic wave flow. The initial condition for this (mono-chromatic) influx is:
where a is the wave amplitude and the tanh is used as a ramp function. It increases smoothly from 0 to 1, to create a smooth start. For t large, the boundary condition reduces to:
If we want to avoid reflections at the right boundary (output), the domain is extended with a damping region of length Ldamp. In this case, terms like and may be added to the continuity (136a) and momentum (136b) equations, respectively. The length of the damped region is chosen such that we do not see any significant reflections.
4.3.3. Solitary wave travelling up a slope and reflection on a vertical wall
Experimental data and numerical results are available for a solitary wave propagating on the bathymetry shown in Figure 8 [1, 2]. It shows a constant depth before x = 55 m and a slope 1:50 between x = 55 m and x = 75 m. An impermeable vertical wall is placed at x = 75 m, corresponding to fully reflecting boundary conditions. A solitary wave of amplitude 0.12 m is initially centered at x = 25 m. The computational domain was uniformly discretized with a spatial step m. A zero friction factor has been considered. Computations were carried out with a time step 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  and ).
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 , contrary to the Serre model, which is . A visual comparison of numerical results of the extended Boussinesq approximation with a similar study performed by , using the extended Boussinesq model developed by , shows no relevant differences in the graphs (see [1, 2] and ).
4.3.4. Periodic wave over an underwater bar
Beji and Battjes (1993)  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 and , and the extended Serre equations (75) (SERIMP model, in  and ), 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 m. A time step s was used. Globally, numerical results of the improved Serre and Boussinesq models agree quite well with the measured data (for details see  and ).
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 . In this experiment, the dispersion parameter () 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 () 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 ). 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.
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.