## Abstract

This paper, which is largely the fruit of an invited talk on the topic at the latest International Conference on Coastal Engineering, describes the state of the art of modelling by means of Boussinesq-type models (BTMs). Motivations for using BTMs as well as their fundamentals are illustrated, with special attention to the interplay between the physics to be described, the chosen model equations and the numerics in use. The perspective of the analysis is that of a physicist/engineer rather than of an applied mathematician. The chronological progress of the currently available BTMs from the pioneering models of the late 1960s is given. The main applications of BTMs are illustrated, with reference to specific models and methods. The evolution in time of the numerical methods used to solve BTMs (e.g. finite differences, finite elements, finite volumes) is described, with specific focus on finite volumes. Finally, an overview of the most important BTMs currently available is presented, as well as some indications on improvements required and fields of applications that call for attention.

## 1. Introduction and motivations

Over the last 40 years, Boussinesq-type models (BTMs) have enjoyed increasing favour within the coastal engineering community. The main reasons for this are the optimal blend of physical adequacy (i.e. their ability to represent all main physical phenomena) and computational ease (i.e. their mathematical well-posedness and numerical cheapness). Hence, Boussinesq-type equations (BTEs), beating the competition of non-linear shallow-water equations (NSWEs), have become the most favoured approximations of Navier–Stokes equations for coastal-type computations.

Such computations are essential tools for any design activity in which water waves play a significant role. Examples of such activities can be found in the field of offshore engineering (design of offshore platforms, pipelines, underwater cables, etc.), near-shore engineering (design of harbours, coastal defence structures, etc.), naval engineering (design of vessels with optimized properties, navigation channels, etc.) and environmental management (evaluation of morphodynamic evolution, assessment of pollutant and nutrients flows, etc.). Applications are now also increasing in the field of riverine engineering.

The interest in BTMs within the scientific community can be simply illustrated by the shear number of journal papers published on the subject (figure 1), beyond the many interesting review papers on the subject [1–4]. Such a success would have surprised my PhD mentor Prof. D. H. Peregrine, who, having provided the fundamental thrust to Boussinesq-type modelling [5], near the end of my PhD (late 1990s) was becoming less interested in BTEs.

However, Peregrine, at that time, was looking at BTEs from a purely scientific point of view, while, as mentioned above and detailed in the following, the fortune of BTEs comes from a mixture of both scientific and applicative interests. This is testified also by the interplay between physics, mathematics and numerics that characterizes the history of BTMs; any progress from the seminal work of Boussinesq [6] has come as the balance of contributions from the mentioned disciplines.

Figure 2 illustrates such an evolution on the basis of some of the most influential researchers in BTMs. Boussinesq, facing the problem of solving the water-wave problem in shallow waters, resorted to perturbation expansions as the most suited techniques of the time. When computers became available for scientific-type calculations, it was possible to couple the power of mathematical and numerical techniques, and hence the pioneering contribution by Peregrine [5]. Under the thrust of coastal engineering needs, in particular, beach protection, BTMs were extended to cover the dynamics of the surf zone. This required some important effort in properly implementing the mechanics of wave breaking within BTMs. The concept of ‘surface roller’ first used by Svendsen [7] has been fundamental. Other fundamental steps towards the BTMs that are currently used as both research and design tools were the extension of BTMs to cover a wide range of the coastal area (i.e. from intermediate to shallow waters); this was largely the fruit of Per Madsen & Sørensen [8] and the implementation of fully nonlinear BTMs, among which the one of Kirby and co-workers is certainly the mostly used [9,10].

The rest of the paper is organized as follows. The fundamentals of Boussinesq-type modelling are introduced in §2, and §3 describes the fundamental improvements that have enhanced the value and validity of such modelling over the years. The fundamental applications of BTMs are illustrated in §4 (e.g. modelling of harbour waves, tsunamis, flow mixing, etc.). The most important numerical issues related with BTM are discussed in §5 and a summary of currently available BTMs is given in §6. §7 summarizes all the analyses and proposes lines for future research.

## 2. Fundamentals

The fundamentals of BTMs are proposed here using a schematic approach, which aims at highlighting the fundamental principles of the modelling and the approach for the choice of the model equations, the characteristics of the chosen equations and the initial applications of the models.

### (a) The modelling

Any BTM to be used for predictive purposes is built by:

— proper account of the physics to be described (e.g. wave propagation, fluid–sediment interactions, etc.);

— choice of the most suited model equations (e.g. best description of nature or idealized dynamics);

— use of appropriate boundary conditions (BCs); and

— use of the most suited numerical methods for the chosen equations and BCs.

### Remark 2.1

Any BTM is more than the constituting BTEs, and the interplay between the above four ‘ingredients’ is an essential feature of each specific BTM.

### (b) The approach

The above fundamental ingredients of any BTM can be balanced as functions of the chosen modelling approach. The following two main approaches can be singled out, which largely depend on the scientific community.

— Mathematicians use model equations for inspection of mathematical behaviours (integrability, well-posedness, symmetries, etc.). Examples are works such as those by Bona & Sachs [11], Bona & Tzvetkov [12], Kishimoto [13], etc.

— Physicists/engineers use equations derived from physical arguments and for practical applications. In this case, many fundamental physical phenomena are to be accounted for and mathematical properties are subordinate to physical requirements.

The present contribution focuses on the latter approach; hence, specific references are proposed throughout the entire publication.

Once the approach has been established, which gives a prominent role to the physical mechanisms to be modelled, there is still room for different modelling choices based on specific mechanisms that are believed to be the dominant ones or the most relevant to the specific purposes at hand. These choices essentially lead to two main classes of equations.

#### (i) Equations for irrotational flows

The first class of BTEs is derived upon the assumption of irrotational flow. This means that the effects of vorticity injection at the bottom and free-surface boundary layers are taken as negligible. Hence, this class of equations best reproduces the dynamics of deep-to-intermediate depth flows. For these BTEs, fundamental dissipative mechanisms are introduced *a posteriori* (and somehow disjoinedly with the potential flow assumption), and are the externally added wave-breaking dissipation and parametrized bottom friction (typically of Chezy type). Externally added breaking-induced dissipation has been successfully achieved on use of the ‘surface roller’ concept; extra terms can be added to the momentum equation in the form of either viscous-type terms [14,15] or convective-type terms [16].

The choice of focusing on the dynamics of flow ‘irrotational at large’ is at the basis of the most widely used BTMs. One of these comes from the Danish Hydraulic Institute (DHI) Danish Technical University (DTU) school [8,16] and is known as MIKE21, whereas the other comes from the Delaware school [9,15] and is known as FUNWAVE. These two models have largely evolved in time from their original formulations (see below), and have found great fortune within the scientific (more FUNWAVE than MIKE21) and practitioner (more MIKE21 than FUNWAVE) communities. FUNWAVE is an open-source solver.

The core of the general procedure [17] to obtain BTEs to be used for this class of BTMs can be summarized as follows (dimensionless variables used):

— make an expansion of the velocity potential in powers of the vertical coordinate and in terms of the depth-averaged velocity or the particle velocity components at a chosen level (e.g. the still-water level, the sea bottom, mid-depth, etc.);

— satisfy the Laplace equation and establish recurrence relations;

— satisfy the kinematic bottom condition and link the horizontal and vertical velocity components ( and );

— apply the established power series for the velocity potential in combination with one of the following:

(i) depth-integrated continuity and momentum equation;

(ii) kinematic and dynamic surface BCs; and

— time step the result (i.e. evolve in time the solution).

#### (ii) Equations for rotational and turbulent flows

The second class of BTEs (and of BTMs) is derived by relaxing the assumption of irrotational flow. In this case, dissipative mechanisms induced by wave breaking, rather than by being externally added, naturally stem from the equation derivation procedure [18], and the model equations seem to be more suited to reproduce surf zone flows. Because of the derivation procedure, model equations require closures for both vorticity and turbulence. Solution of such extra equations can be pursued either by semi-analytical means [19] or by numerical calculations [20].

Derivation of this class of BTEs is achieved through the following main actions.

— Depth averaging, using the surface BCs, the continuity and Reynolds equations. Care must be taken because turbulence and space averaging do not commute.

— Decomposing the velocity field into a sum of irrotational and rotational contributions.

— Solving the Poisson equation, for either the stream function

*ψ*or the vertical flow velocity*w*, and expanding in*μ*(=*kh*)≪1.— Relating the bottom velocity

*u*_{b}to the reference velocity (here the depth average of*u*) and writing equations in terms of .— Adding closures for vorticity and turbulence.

The class of ‘naturally rotational’ BTMs, though being used in various scientific studies, has found much less favour than that described above and practically no use is made of these models at the practitioner's level. The reasons for such a partial success are varied, not least the difficulty in physically describing the required closures and numerically implementing them into an efficient code.

Focusing more on the nearbed vorticity effects, Kim *et al.* [21] have devised an approach that seems of interest for the inclusion into BTMs of bottom boundary-layer dynamics.

### (c) The model equations

The above analyses have set the scene for a more detailed description of the model equations. We use a Cartesian coordinate system with the *x** and *z** axes lying on the still-water level and, respectively, pointing seawards and upwards, i.e. *x**=0 gives the still-water shoreline (figure 3). In the following, asterisks are used to label dimensional quantities, whereas plain variables give dimensionless quantities.

Neglecting boundary-layer dynamics (i.e. modelling the flow above the dotted line of figure 3), a vast majority of BTEs of the above-mentioned classes can be formally written as follows:
2.1
2.2
2.3
2.4*a* is a reference wave amplitude, *k* is the wavenumber and the comma denotes partial differentiation.

The continuity equation (2.1) is exact, whereas the momentum equation (2.2) is given in conservation form. Equations (2.3) and (2.4) are solved for the total water depth (*d*) and for a reference velocity () once closures (2.3) are given. Because these may be of many different forms, they are given here by a generic formal representation.

As for all BTEs, solution of the system (2.1)–(2.4) is a function of the size of the wave nonlinearity (*ϵ*) and frequency dispersion (*μ*) parameters.

### (d) Early applications: harbour waves 1

The first applications of BTMs to real-life problems date back to the 1970s and 1980s, and were related to the description of harbour free-surface oscillations.

Long-wave oscillations in harbours (seiches) arise as the result of nonlinear interactions between narrow-banded waves and swell. Such oscillations may cause excessive ship motion and negatively interfere with loading and unloading operations at port facilities. Hence, predictions of such oscillations (e.g. figure 4*a*), of the related ship motions (e.g. figure 4*b*) and mooring forces were made by means of weakly nonlinear and weakly dispersive BTMs [22–25].

## 3. Improvements

The early BTMs, i.e. those developed up until the 1980s, had so many limitations in their range of applicability that their use was basically restricted to harbour waves and related ship motions. A deep analysis of weakly nonlinear, weakly dispersive BTMs can be found in the book by Dingemans [26].

However, because of the increasingly pressing requirements of modelling tools for coastal engineering design purposes, BTMs flourished to become the most widely used modelling instrument for engineers and practitioners of the coastal zone.

### (a) The 1990s

As recalled earlier, the temporal evolution of BTMs has been largely governed by practical requirements. Among the most pressing needs of the coastal engineering community were those of extending the region of validity of BTMs towards the intermediate-to-deep waters, on one hand, as well as towards the near-shoreline dynamics. The 1990s have seen the most rapid improvements, in both directions, undergone by BTMs in their history, and at the end of the 1990s, largely because of the contributions from the DHI group, BTMs were capable of satisfying most of the modelling requests raised by practitioners.

#### (i) Dispersive properties: frequency dispersion (*O*(*μ*^{2})=1)

The early BTMs, all stemming from the work of Peregrine [5], were weakly dispersive and weakly nonlinear models, i.e. strictly speaking, such that (*O*(*μ*^{2})=*O*(*ϵ*)≪1). This meant that computations were limited to the case of very long waves (e.g. seiches), whereas much of the coastal zone design activities (related to, e.g. beach protection, wave–structure interactions, etc.) are largely influenced by the actions of short waves, i.e. waves (such as wind waves) not much longer than the local water depth.

*Extra terms in the momentum equation.* Improvement of frequency dispersion properties can be, in principle, achieved by simply retaining contributions to the governing equations of order higher than *O*(*μ*^{2}). For example, clear improvements are obtained by retaining *O*(*ϵμ*^{2},*μ*^{4}) terms in the momentum equation. However, this also means introducing in the momentum equations derivatives of order higher than the third (typical of *O*(*μ*^{2})≪1 models), hence largely increasing the model implementation difficulties and the computational costs.

Madsen *et al.* [27], and followers, fully understanding this problem, found a smart approach for improving dispersive properties without increasing the order of the derivatives appearing in the momentum equation. In brief, they applied a linear operator () to the momentum equation, thus obtaining extra dispersive contributions, proportional to the free parameter *B*, which was chosen to match the [2,2] Padé approximation of the dispersion relation. These extra term, however, only contained third-order derivatives such as the original weakly dispersive, weakly nonlinear BTEs.

*Chosen reference velocity.* The method mentioned above was not the only one used to improve dispersive properties without introducing higher-order derivatives. Witting [28] and Nwogu [29] resorted to the use of reference velocities different from the depth-averaged ones. In particular, Nwogu [29], exploiting the fact that dispersive properties strongly depend on the actual velocity profile, introduced a free parameter by choosing the velocity evaluated at an arbitrary *z*-location as the reference velocity. The resulting equations contained third-order derivatives in the mass and momentum equation, and the free coefficient was chosen to match the Padé [2,2] approximation.

Figure 5 illustrates the improvements brought to the linear dispersive relation by use of the two above-mentioned methods.

### Remark 3.1

The above improvements are important examples of how physical/mathematical insight influences the subsequent numerics; dispersive properties of BTMs can be improved with no need to use derivatives higher than third order.

#### (ii) Wave-breaking dissipation (externally added)

This has been the most fundamental improvement that made BTMs become the standard for near-shore flow circulation modelling; many coastal engineering activities occur within the surf zone, and the ability to reproduce the surf zone hydrodynamics was essential for any modelling tool.

As illustrated in §2, virtually all the BTMs of the 1990s assumed the flow to be irrotational. Hence, dissipative dynamics associated with wave breaking had to be introduced by means of externally added contributions.

Two main approaches were used to suitably prescribe such extra contributions, both of which were based on the evidence that shallow-water breaking waves soon become ‘bores’, i.e. waves of quasi-permanent form whose shore-facing front is characterized by a recirculating mass of water known as the ‘surface roller’.

*The roller model.* The concept of ‘surface roller’ by Svendsen (see [7]), is the basis for this approach. The roller is taken as a volume of water, of thickness *δ*, passively riding, at the wave celerity *C*, on the front face of the wave (figure 6). The effect of the roller is included through the extra convective term *ϵ*(Δ*M*),*x* in the momentum equation, which leads to an energy dissipation (Δ*E*) comparable to that of a hydraulic jump of height *H*_{b} equal to that of the breaking wave (Δ*E*≈*H*^{3}_{b}) [16].

*The eddy-viscosity model.* In this case, the dissipative mechanisms related to wave-breaking global turbulence are introduced by means of a diffusive-type extra term () in the momentum equation. The eddy viscosity is computed from integration of the transport equation of the turbulent kinetic energy and the mixing length hypothesis [14,30].

### Remark 3.2

This improvement is another important example of the fundamental role of physical modelling on the subsequent numerics; the dissipation model can include either first-order derivatives *ϵ*(Δ*M*)_{,x} or second-order ones .

### (b) The 2000s

Obviously, the improvements of the 1990s were pushed forwards also into the 2000s. For example, with reference to dispersive properties, Schäffer [31] and Karambas & Memos [32] derived BTMs with exact linear dispersion characteristics. However, despite their improved dispersion relation, the 1990 extended Boussinesq equations were still restricted to situations with weakly nonlinear interactions. In many practical cases, however, the effects of nonlinearity are too large to be treated as a weak perturbation to a primarily linear problem. As waves approach the shore, wave height increases owing to the effect of shoaling, and wave breaking occurs on most gentle natural slopes. The wave height to water depth ratios accompanying this physical process are inappropriate for weakly nonlinear Boussinesq models, and thus extensions to the model are required in order to obtain a computational tool that is locally valid in the vicinity of a steep, almost breaking or breaking wave crest.

This was the focus of research in BTMs over the first decade of the 2000s. The main improvements made for the analysis of the propagation of nonlinear waves and of breaking waves are summarized in the following.

#### (i) Nonlinear properties: large bed slopes: shoaling

For practical calculation reasons, i.e. to establish finite-order equations that can be solved numerically, the infinite-series expansions used in irrotational-type BTEs (see §2) 3.1are approximated by either truncated Taylor's series or Pade's expansion around an arbitrary level .

Significant improvements are obtained for both linear properties (i.e. linear frequency dispersion, linear shoaling, velocity profiles) and nonlinear ones (i.e. amplitude dispersion). An example of the improvements in amplitude dispersion as functions of the chosen truncation approach is illustrated in figure 7, taken from Madsen *et al.* [33].

If the expansions used for the irrotational-type BTEs (see §2 and above) are made without using the mild-slope approximation (∇_{H}*h*(*x*,*y*)≪1), i.e. if the spatial derivatives of are also retained in the expansions, BTEs for flows interacting with rapidly varying bathymetries can be obtained [34,35].

The model obtained enables good shoaling performances, as well as a good description of the interaction between waves and very steep slopes. The example of figure 8, taken from Bingham *et al.* [35], illustrates the performances of various models in the case of waves reflected at a steep shelf (plane slope of horizontal length *b*_{0} connecting two constant-depth regions).

#### (ii) Wave-breaking dissipation (intrinsic)

The 2000s have seen a ground-breaking innovation for BTMs: a conceptually novel model enables us to account for wave-breaking-induced dissipation by means of an intrinsic, rather than external, approach. In other words, the model by Veeramony & Svendsen [18] is built, so that rotational and breaking wave contributions naturally arise in the derivation of the model equations. The structure of such BTEs valid for rotational and turbulent flows (see also §2) can be summarized as follows:
3.2
3.3
and
3.4solved in the form *ω*=*ω*(*ω*_{s},*z*,*h*,*ζ*_{e}), so that *D*_{turb}, Δ*M*, Δ*M*_{1} and Δ*P*=*f*(*h*,*ζ*_{e},*ω*_{s}),
3.5These equations show that the momentum equation is characterized by terms from both irrotational (*u*_{I}) and rotational (*u*_{R}) flow contributions. Hence, at each order [*O*(*ϵ*),*O*(*μ*^{2}),…], beyond the classical Boussinesq terms, dissipative terms (Δ*M*,Δ*M*_{1},Δ*P*,*D*_{turb}) appear, which can be computed once the vorticity transport equation is solved and the vorticity *ω* itself substituted into the expression of the rotational velocity (*u*_{R}=*u*_{R}(*ω*,*μ*^{2})).

In other words, a coupled model is built in which the vorticity transport equation can be solved either analytically [18,19] or numerically [20]. In any case, the solution is found to strongly depend on the surface vorticity *ω*_{s} injected in the fluid body. This is derived from values found at the ‘dividing streamline’ of experimental hydraulic jumps (see the sketch in figure 9).

The model enables a good representation of the internal kinematics of breaking waves both outside and inside the surf zone. This is testified by figure 10, which illustrates the comparison between the results of the fully nonlinear model of Musumeci *et al.* [19] and the experimental results collected by Cox *et al.* [36] at the shoaling zone, at the breaking region and inside the surf zone of a laboratory planar beach. Not only is the water surface elevation reproduced fairly well, but also the computed vertical profiles of the horizontal velocity well match the experimental ones (this is different from most other breaking models, which, being characterized by self-similar vertical profiles, do not allow for a correct description of the ‘undertow’).

Notwithstanding the fundamental theoretical innovation of describing vortical and dissipative dynamics intrinsically to the model equations and the good modelling performances, this approach has not found, as yet, a large favour, neither among scientists nor among practitioners. Among the possible reasons for this lack of success are (i) the strong dependence of the model results on the specific values of surface vorticity *ω*_{s} to be used as BCs in the model, (ii) the need for solving in a coupled way the BTEs and vorticity equations, and (iii) the model being derived only for one-dimensional flow propagation.

### Remark 3.3

This is the third example of the influence of the physics on the numerics: the breaking-induced dissipative terms Δ*M*,Δ*P*,… introduce higher-order derivatives in the model.

### (c) Towards the 2010s

#### (i) Multi-layer, non-hydrostatic models

The most recent improvements in BTMs are being made, in view of the increasing computational power available, to make BTMs become closer and closer to three-dimensional solvers, i.e. capable of also properly describing the vertical structure of the flow, leading to distinctive improvements in the nonlinear-dispersive properties of the model.

The guiding idea is to improve such nonlinear-dispersive properties by using:

— a better flow resolution in the vertical (e.g. multi-layer models of Lynett & Liu [37] and Lynett [38]) and

rather than by improved power series expansion or higher-order derivatives in the momentum equation. In this perspective, the model equations are discretized over the vertical by means of a few layers (typically two to five).

Both approaches lead to excellent results for

— dispersion;

— shoaling; and

— breaking and run-up.

An illustration of such good performance is provided by figure 11, which collects a number of findings (i.e. wave celerity, linear shoaling, two dimensional in the horizontal (2DH) run-up and inundation) obtained by both approaches.

## 4. Applications

The advances and progress in BTMs can be described not only through the analysis of the physical, mathematical and numerical ingredients that are used to build a specific model, but also on the basis of the actual applications of the models. The most important are described here through some examples whose illustration have been made available by the contributions of many colleagues and friends. Other applications, not discussed here for the sake of brevity, are calculations of wave overtopping, of wave–structure interactions and interactions of waves with a porous seabed.

### (a) Harbour waves 2

Harbour waves and ship motions have been the object of BTMs because their initial application (see §2) is still being studied by means of BTMs. Figure 12 illustrates the harbour wave fields induced by wind waves.

### (b) Complex bathymetries

Much more testing are applications associated with the propagation of waves over complex bathymetries. For these computations, most of the characteristics of the models are put under serious scrutiny, such as, for example dispersive properties, shoaling behaviour and interaction with possibly steep slopes. Other properties of importance for this sort of computations are the wave–wave interaction (e.g. at the lee side of an obstacle) and the existence of multiple shorelines. The latter mechanisms are visible in figure 13*b*.

### (c) Tsunamis

Because of the growing societal impact of tsunami waves and in virtue of the modelling improvements made over the last 30 years, applications of BTMs to the description of the evolution of tsunamis are growing. Although, theoretically appropriate only for the far-field modelling of such long waves (i.e. far from the generation point of the tsunami wave), BTMs are practically being used to describe tsunamis over their entire lifespan (e.g. [42]). Figure 14 illustrates the computations of laboratory solitary waves (Figure 14*a*) and of field tsunamis (Figure 14*b*).

### (d) Flow mixing

Another recent application of BTMs is in the field of study of shallow-water turbulence and water quality management. BTMs are being used to run numerical experiments whose outputs are used both to analyse Lagrangian statistics of a large number of water particles and to assess the main characteristics of the horizontal flow mixing.

As an example of this sort of application, figure 15 shows the maps of water particles released in a longshore bar–trough system. The initial distribution of a large number (10^{4}) of particles (figure 15*a*) evolves because of the flow field computed by FUNWAVE. The later stage maps (figure 15*b*) enable the visualization of the flow topology, the comprehension of the differential transport (in time–space) and the illustration of the structure of vortices, jets and cross-shore transport.

### (e) Morphodynamics

One fundamental field of application of BTMs is near-shore morphodynamics. Used in conjunction with simplified morphological evolution equations (i.e. Exner equations) and simplified sediment transport closures, BTMs have been used to model many specific features of the near-shore sediment transport, morphodynamics and morphology. Erosion and accretion of a sandy beach (figure 16*a*,*b*) as well as water infiltration/exfiltration at a porous beach (figure 16*c*,*d*) are among the most typical morphodynamic computations made with BTMs.

## 5. Numerics

Once the physics of interest and the mathematics in use have been described, the final aspect to be analysed is the numerics used to carry out suitable BTM computations. This aspect, which was a specific target of the original International Conference on Coastal Engineering talk, is discussed in some detail.

### (a) Early approaches

A brief overview is here given of the numerical approaches used before the 2000s, when finite volumes (FVs) became the dominating method.

#### (i) Finite differences (FDs)

This approach has found much favour because of the wealth of specific literature available and because of the ease of application, as it is a rather intuitive approach for writing difference equations. The best-known FD scheme used in BTMs is that initially proposed by Abbott *et al.* [22], and subsequently used by many authors [8]. In brief, the differential equations were discretized using a time-centred implicit scheme with variables defined on a space-staggered rectangular grid. The method is based on the alternating direction implicit algorithm.

#### (ii) Finite elements (FEs)

This method has found less favour than the FD method because it is relatively more difficult to implement. However, its excellent performances in dealing with complex geometries suggested its use in the case of confined water bodies (e.g. harbours). The early FE BTMs [25,45,46] relied on the Galerkin-weighted residual method and used linear shape functions for interpolation.

### (b) Finite volumes

A more detailed analysis is dedicated to the description of the use of FVs in BTMs because such an approach, which makes the most of the wealth of knowledge in the integration of hyperbolic equations, has become the most widely used in the last 10 years. Following the approach used throughout this paper, rather than illustrating the details of each specific model under investigation, all models are discussed on the basis of some important common features.

Synthetically, all BTEs can be cast in the following matrix form (conservative):
5.1
or
5.2where **U** is the vector of variables (including dispersive terms), **A**_{x},**A**_{y} are components of the advective fluxes, **S**_{B} is the source term owing to the bathymetry, **S**_{D} is the source term owing to dispersive and bed frictional terms, and **D** is the dispersive contribution appearing in the evolution term.

### Remark 5.1

Solution of equations (5.1) and (5.2) entails two main issues.

— Solution for

**U**requires the inversion of a tridiagonal matrix, from which and can be extracted (see underbraced terms on the left-hand side of equation (5.1)).— Schemes for discretizing

**A**_{x},**A**_{y}are, typically, different from those used for**S**_{B}(see overbraced terms on the right-hand side of equation (5.1)). Hence, a numerical imbalance arises, generating numerical noise. From this imbalance comes the need to resort to methods such as the ‘surface gradient method’ (SGM; see [47,48]).

The SGM consists of rewriting the original pressure gradient term as a function of the free surface variable *ζ* [48],
5.3so that *ζ* appears in variables and fluxes (see underbraced terms in equation (5.4)) and **S**_{B} depends on *ζ* and not on *d* (see overbraced terms in equation (5.4)),
5.4
or
5.5

### Remark 5.2

The above is an example of a mathematical preconditioning that positively influences the numerics.

Then, the FV discretization is obtained from the weak form of the original PDE (i.e. equation (5.4)). Hence, integration over a cell-centred (*i*,*j*), at time *n*, space–time volume *Ω*=[*t*^{n},*t*^{n+1}]×[*x*_{i−1/2},*x*_{i+1/2}]×[*y*_{j−1/2},*y*_{j+1/2}] gives
5.6where *k* gives the component and the overbar a cell average.

Integration of the above evolution equation requires

— an average reconstruction (represented by the overlines in equation (5.6));

— flux/source computation (underbraced terms in equation (5.6)); and

— time integration (superscripts

*n*and*n*+1 in equation (5.6)).

In the majority of currently available FV BTMs:

—

*Weakly nonlinear BTEs are integrated*[8,29]. Examples of exceptions are the models of Kim*et al.*[21] and Shi*et al.*[49], which integrate fully nonlinear equations.— FV methods are used for the conservative (left-hand side) part of the equations, while FD schemes are used for the source terms

*S*_{B}and*S*_{D}:*hybrid FV–FD models*. Examples of exceptions are the models of Roeber*et al.*[50], Roeber & Cheung [51] and Dutykh*et al.*[52], which use fully FV methods.— The data reconstruction (averages) is performed by means of the

*fourth-order monotone upstream-centred schemes for conservation laws (MUSCL) method*. Examples of exceptions are the models of Cienfuegos*et al.*[53] (Padè interpolation) and Roeber*et al.*[50] (multi-dimensional limiting process).— The flux computation is performed by using the

*approximate Harten, Lax and Van Leer finite volume scheme (HLL) and the Harten, Lax and Van Leer finite volume scheme for contact waves*(*HLLC*). Examples of exceptions are the models of Erduran*et al.*[54] (Roe scheme), Roeber*et al.*[50] and Dutykh*et al.*[52] (various schemes).— The time stepping is performed by either the

*fourth-order predictor-corrector Adams–Bashforth–Moulton*(ABM) scheme [50,54–56] or the*Runge–Kutta scheme*[49,52,53]. Examples of exceptions are the models of Ning*et al.*[57] and Shiach & Mingham [58] (second-order MUSCL–Hancock, claiming second order is enough for time stepping).

The above analysis reveals that only some of the inspected FV models are able to deal with wave-breaking dissipation. Furthermore, it is interesting to note that two main approaches are used to such a purpose, i.e.

— The BTE+NSWE matching. This is based on the assumption that in very shallow waters, dispersive terms become negligible and BTEs collapse onto NSWEs. Hence, wave-breaking-induced dissipation is handled by Rankine–Hugoniot conditions associated with the NSWE [49,55,56].

### Remark 5.3

This is a further, important example of how physical arguments (i.e. BTEs collapsing onto NSWEs in shallow waters) influence the numerics (i.e. modelling of wave breaking through Rankine–Hugoniot conditions).

— The eddy-viscosity approach. If a fully FV approach is used instead of a hybrid approach, typical models for externally added breaking can be used [50]. The most favoured one is the eddy-viscosity approach.

### Remark 5.4

Another example of the interplay between physics and numerics comes from the relationships between the intrinsic dissipation owing to Rankine–Hugoniot conditions and eventual externally added breaking models. No analysis for such an interplay is still available.

### Remark 5.5

Surprisingly enough, notwithstanding the physical importance of the roller model, this seems to be abandoned by modellers who apply FV methods.

### (c) Finite elements

Although still lagging behind FV methods, FE methods are increasingly being applied to modern BTMs in virtue of the power and flexibility of the discontinuous Galerkin (DG) approach.

For the sake of simplicity, the derivation of the scheme is briefly given starting from a simple, generic, scalar conservation law, 5.7

The main steps can summarized as follows.

— Define the FE subdivision of the domain 5.8

— Use of local basis functions. The DG approach makes use of elemental basis functions (figure 17

*a*) instead of nodal basis functions (figure 17*b*), 5.9 5.10— Write the discontinuous (local) Galerkin residual (

*R*^{DG}_{i}) and minimization 5.11and, using integration by parts, write the local weak form 5.12where the intercell fluxes (IFs) are computed such as in the FV method.

Among the important models and applications, the following can be mentioned.

— The model of Eskilsson & Sherwin [59]. This solves the equations by Madsen & Sørensen [8] using the HLLC scheme for the IF. The model performances are assessed here on the basis of the propagation of waves over a semicircular shoal (figure 18

*a*). In particular, the amplitudes of the first three harmonics obtained from the model are compared with those from the laboratory experiments of Whalin [60] to show a fairly good match (figure 18*b*).— The model of Engsig-Karup

*et al.*[61]. This model solves the equations of Madsen*et al.*[62] through a Lax–Friedrich scheme for the IF. Application of the model to the scattering of linear waves about a vertical cylinder in open water (figure 19*a*) and comparison with the analytical solution owing to McCamy & Fuchs [63] reveal the good performance of the model.

## 6. Summary of models

A summary of the main characteristics of the many BTMs explored here is given in the form of table 1. Clearly, neither the list of inspected models nor the features chosen to characterize such models are claimed to be exhaustive. However, they provide a rather good synthetic description of the state of the art in BTMs.

## 7. Conclusion and future

The state of the art of BTMs has been illustrated, with special attention on the interplay between the physics to be described, its mathematical modelling and the related numerics. Fundamentals (§2), time evolution (§3), applications (§4) and numerics (§5) of BTMs have been discussed from the point of view of a physicist/engineer rather than of an applied mathematician.

Lines for ongoing and future research close the paper. Following the approach used throughout the paper, these are classed into two main groups: ‘improvements’ and ‘applications’, the former focusing more on basic scientific research issues and the latter more on practical matters.

### (a) Improvements

A number of improvements to BTMs is proposed here as a function of their increasing level of urgency.

#### (i) Low urgency

Offshore and shoreline BCs, issues of some importance, are well-explored matters nowadays, with a wealth of knowledge available. It is now clear that numerical fixes cannot be used to hide physical/mathematical shortcomings in the prescription of BCs [71–74].

Similarly, for the breaking criteria. Although, better modelling of the instability and flow separation at the breaker front must be sought [75], robust and reliable parametrizations of the inception of breaking are available [16,76,77].

#### (ii) Medium urgency

Of medium urgency are issues such as: (i) the implementation of high-order non-irrotational BTMs, (ii) the suitable inclusion of wave breaking within non-hydrostatic models and vertically resolved models and (iii) the proper modelling of rapidly varying bathymetries. These matters are the focus of ongoing research. For example, rather innovative approaches are being made to avoid the artificial BTE + NSWE matching currently used in many BTMs [78] and to build high-order non-hydrostatic BTMs [79].

#### (iii) High urgency

Finally, the most urgent needs in BTMs, i.e. the issues that are believed to bring in the most important progresses, are (i) the joint modelling of intrinsic (i.e. Rankine–Hugoniot conditions) and externally added breaking dissipation, (ii) the adaptive mathematical/numerical description of the flow (i.e. Lagrangian, moving grid, curvilinear-unstructured grids), (iii) boundary-layer (surface and bottom) inclusion (with specific emphasis on vorticity injection mechanisms) and (iv) coupling with oceanographic models.

### (b) Applications

Applications, i.e. dynamics that can be usefully modelled by means of BTMs, are illustrated in their increasing order of societal impacts.

#### (i) Low importance

Sloshing in tanks is a typical application of BTMs [74,80,82], but it does not seem to be a leading-order problem to be tackled by means of BTMs.

#### (ii) Medium importance

Of medium importance are environmental-type applications such as: (i) the description of riverine flows [83], (ii) the flow interaction with vegetation (in rivers, seas and estuaries, e.g. [84]) and (iii) the horizontal mixing of natural streams [85].

#### (iii) High importance

Among the most important applications, I think we should list: (i) the (coupled–decoupled) modelling of the morphology of natural streams [86], (ii) the flow–structure interactions and (iii) the assessment of devices to be used for the extraction of energy from sea waves [87].

I would like to close the paper by responding to a challenging question posed by one reviewer who asked me to speculate on the likely main roles of BTMs going forward, in recognition of the fact that three-dimensional non-hydrostratic models are reaching a degree of maturity and are simpler to implement and understand than higher-order BTMs. In summary: ‘do BTMs have an obviously persistent niche, or will they end up being a developmental stage?’. In line with the perspective motivating the paper, I think that BTMs can still play an important role in the modelling of near-shore dynamics if the interplay between physics, mathematics and numerics is used to design models that are ‘smarter’ (i.e. superior in the physical insight, mathematically simpler and numerically cheaper) than the competing ones.

## Funding statement

I acknowledge the financial support from the Italian RITMARE Flagship Project, a National Research Programme supported by the Italian Ministry of University and Research. Financial support from the EU, through the ENVICOP project (PIRSES-2011-295162) and from the US-ONR, through the NICOP research grant (N62909-13-1-N020) is also gratefully acknowledged.

## Acknowledgements

Many colleagues contributed, in various ways, to the present review work, too many to be recalled here in detail. My most sincere thanks to all of them. Special thanks to A. G. L. Borthwick, K. F. Cheung, J. T. Kirby, M. H. Kobayashi, P. J. Lynett, J. Orszaghova, M. Petti, V. Roeber, F. Shi, P. H. Taylor and M. Tonelli, who provided useful material for this paper. The constructive criticism by Dr P. A. Madsen and Dr R. Briganti is gratefully acknowledged.

- Received July 23, 2013.
- Accepted September 16, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.