Wave modelling

Ocean wave modelling efforts and applications can be broadly classified into two large groups: i) phase resolving (or direct) models; and ii) phase average (usually spectral) models. Direct models can explicitly simulate basic equations of fluid mechanics for the water, air, or even two-phase media, and therefore extend the analytical research beyond its traditional range of approximate and asymptotic solutions of such equations. At oceanic scales, however, such models are not practical and not feasible, and therefore spectral models are employed for wind-wave forecasts.

In the next subsections, for both deep and shallow water analysis is given a general description, mathematical model, limitations, and main applications for each type of model.

Spectral models

Evolution of wind-generated waves in water of finite depth d can be described by the wave action N=F/ω balance equation:

where F(_ω,k) is the wave energy density spectrum, _ω is intrinsic (from the frame of reference relative to any local current) radian frequency, k is wavenumber (bold symbols signify vector properties). In the linear case, temporal and spatial scales of the waves are linked through the dispersion relationship (see Eq. 8.2).

The left-hand side of Eq. 8.11 represents time/space evolution of the wave action density because of the energy source terms on the right. On the left, c_g is group velocity, c_k means the spectral advection velocity, U is the current speed, and we note that c =ω/_k is phase speed of the waves. ∇ here is the horizontal divergence operator, and ∇_k is such an operator in spectral space.

On the right, source terms are physically represented by wind energy input from the wind, I; nonlinear interactions of various orders within the wave spectrum, L, whose role is to redistribute the energy within the spectrum; dissipation energy sinks, D; wave-bottom interaction processes, B; and more sources are possible in specific circumstances. Note that all the source terms, as well as the group and advection velocities, and the advection current are spectra themselves. Please refer to Cavaleri et al. (2007) for further details.

Among the source functions, L is a conservative term, i.e. its integral is zero, but the other integrals define energy fluxes in and out the wave system:

is the total flux of energy from the wind to the waves. Note that, depending on the relative speed of wind U₁₀ and wave speeds c(ω,k)=ω / k, contributions to the total flux can be both positive (from the wind to the waves if U10>c) and negative (from the waves to the wind if U₁₀>c). In the tropics, for example, where the wave climate is dominated by swells produced at high latitudes, the local winds are typically light and therefore the wind climate can be actually dominated by wave-induced winds (Hanley et al., 2010).

It should be noted that the energy input to the waves is generally accepted as a purely atmospheric exchange. In principle, however, energy input from the ocean side to the surface waves of scales accommodated in Eq. 8.11 is perceivable. For example, upper-ocean currents, tides, or internal waves can provide such dynamics. Given the amount of energy stored in the ocean movements, this could have large impacts on surface wave fields, even if localised, but it is fair to say that it has not been considered by the wave-ocean modelling community in practical terms. 

Integrating the momentum-input spectrum gives the total momentum flux:

which is an important measure of wind-wave interactions (Tsagareli et al., 2010). Together with the tangential viscous stress τ_v it forms the total wind stress at the ocean surface

and this stress is known independently (usually through empirical parameterisations of the so-called drag coefficient) and thus can be used as a constraint or for validation of the wind input term I. On the other hand, the total stress is often the main, if not the only property which expresses dynamic exchanges in large-scale air-sea models. Apart from situations of light winds, the wave-induced form drag (Eq. 8.12) provides a dominant contribution to this total stress (Kudryavtsev et al., 2001) and thus, if the wave-model physics is well defined and validated, such models can provide explicit rather than empirical estimates of fluxes for general circulation models if those are appropriately coupled with wave models.

The dissipation function D has a similar meaning in the context of wave-ocean dynamic exchanges, but with some essential distinctions. First, the integral

is the total flux of energy out of the wave field. The energy passed to the ocean is largely spent on generating turbulence near the surface and on work against buoyancy forces acting on bubbles injected during the wave breaking.

Unlike the input, however, which only occurs on the air side of the interface, the loss (8.15) can go both to the ocean below and to the atmosphere above the ocean surface. Numerical simulations of Iafrati et al. (2013) showed that up to 80% of wave energy due to breaking can be actually dissipated through the atmospheric turbulence.

The momentum-loss integral of dissipation function gives the so-called radiation stress:

which is presumed to be going to the currents (although some of it may in fact be going back to the wind, or to the bottom in shallow areas). In the present wave models, radiation stress is parameterized in terms of wave-height difference along the propagation direction. Obviously, such parameterization does describe the energy dissipation, and can then be used to estimate the momentum loss, but only in the areas where dissipation (Eq. 8.14) is much larger than the energy input (Eq. 8.14), i.e. usually in shallow waters. In deep water, the mean wave height is not a proxy for the energy loss. In fact, it may grow under wind action or not change if this action is balanced by the whitecapping dissipation, but the integral (Eq. 8.16) and hence the radiation stress is not zero.

Wave-ocean-bottom interactions in infinite depths, depicted by term B in Eq. 8.11, are very rich. Finite depths are characterised by the condition of kd~1 (wavelength is comparable with the water depth d), and shallow non-dispersive environments by kd<<1. Dispersive-wave nonlinear dynamics slowdown in finite and shallow depths, weaken or cease, but other nonlinear behaviours come into existence.

Wave exchanges with the bottom include bottom friction, formation of ripples, sediment suspension and transport if the sea bed is sandy, generation of bottom waves if the bottom is muddy, and percolation. Long-shore, cross-shore, and rip currents result from radiation stresses (Eq. 8.15), infragravity waves are produced by combined action of wave breaking and nonlinear wave groups, which can be subsequently reflected back to the deep ocean or trapped by coastal bays.

An example of this deep-water approach is the WAM (Hasselmann et al.,1988), perhaps the first one proposed as third-generation model, able to explicitly represent all the physics relevant for the development of the sea state in two dimensions, such as wind generation, whitecapping, quadruplet wave-wave interactions, and bottom dissipation. This modes is mainly forced by a two-dimensional ocean wave spectrum that develops freely with no constraints on the spectral shape, so that: a) a transfer source function of the same degree of freedom as the spectrum itself need to be developed; and b) the energy balance had to be closed by defining the dissipation source function. Hasselmann et al., (1985) and Komen et al., (1984) were employed to deal with these aspects, respectively. The dissipation was selected in order to replicate the observed fetch-limited wave growth and the fully developed Pierson-Moskowitz spectrum (WAMDI group, 1988).

Constant improvements and updates have led to a third-generation WAM model. A third-generation wave model explicitly represents all the physics relevant to the development of the sea state in two dimensions. Numerical solutions of the momentum balance of air flow over growing surface gravity waves have been presented in a series of studies by Janssen et al. (1989), and Janssen (1991). The main conclusion was that the growth rate of the waves generated by wind depends on the ratio of friction velocity and phase speed and on several additional factors, such as the atmospheric density stratification, wind gustiness, and wave age. This work has also introduced the surface stress dependency with the sea state, and the feedback of wave-induced stress on the wind profile in the atmospheric boundary layer.

WAM is an Eulerian phase-averaged model. Designed as a deep-water model, it can be used to predict directional spectra and wave properties (significant wave height, mean wave direction and frequency, swell wave height). The model can be used in finite depth as well by introducing bottom dissipation source function and refraction. The model runs on a spherical latitude-longitude grid.

The first WAVEWATCH model was developed at TU Delft (Tolman, 1989; Tolman 2014), followed by the NASA Goddard Space Flight Centre in 1992. Recently, WAVEWATCH III was presented as a worldwide used and full-spectral third-generation wind-wave model. It was developed at NOAA/NCEP and it is based on the first WAM model’s principles. This latest version includes many improvements in the governing equations, model structure, numerical schemes, and physical parameterizations. The model solves the random phase spectral action density balance equation for wavenumber-direction spectra. The medium properties, namely the water depth and current properties, as well as the wave field, vary in time and space in scales much larger than a single wave. WAVEWATCH is an open-source model that is freely available 🔗3, including the whole source code and all documentation.

The discretization of the wave energy spectra in all directions is achieved by using a constant directional increment and a spatially varying wavenumber grid, which corresponds to an invariant logarithmic intrinsic frequency. In order to achieve high accuracy, both first order and third order schemes are available for wave propagation. For the integration of source terms in time, a semi-implicit scheme is used similar to that used in WAM, which includes a dynamically adjusted time stepping algorithm.

Following the work of Battjes and Janssen (1978), WAM and WAVEWATCH III models have been upgraded to account for the dissipation by wave breaking induced by depth in the surf zone. However, wave models still have difficulties with strong three-wave interactions that occur in finite-depth and shallow waters. That has led to simplified empirical calculations with large errors, especially for complex wave trains with multi-model spectra. In addition, both models lack of diffraction processes, which implies that only open coastal zones could be solved accurately, plus only linear behaviour of wave propagation could be assessed and non-linear corrections to linear wave should be imposed, by triad and quadruplet wave-wave interactions in shallow waters, where the waves break (Booij et al., 1999).

Spectral models

For shallow water domains and wave propagation (Eckart, 1952), the SWAN model could be a good choice. This also is a third-generation wave model developed at the Delft University of Technology, with the purpose of obtaining realistic estimates of wave parameters in coastal areas, lakes, and estuaries from given wind, bottom, and current conditions. The SWAN model can be also used on any scale relevant for wind-generated surface gravity waves. The model equations are based on the wave action balance equation with sources and sinks.

SWAN has been developed to simulate coastal wave conditions (with friction, breaking, whitecapping, triad, and quadruplet wave-wave interaction). SWAN can be also coupled with previous models such as WAM or WAVEWATCH III, and inherit the boundary conditions. SWAN can provide a computational representation of directional and no directional spectrum at one point, and several spectral and time-dependent parameters of waves, such as significant wave height, peak or mean period, direction, and direction of energy transport. SWAN is a freely available 🔗4 open-source software. 

SWAN model is based on the spectral action balance equation, which describes the evolution of the wave spectrum (Booij et al., 1999). 

In Cartesian coordinates the evolution of the action density is governed by the following balance equation:

where σ is the wave frequency, θ is the wave direction component, t is the time, x and y the 2D coordinates in space, N the wave action density spectrum defined as:

where E is the wave energy density spectrum; S_total is the source term and C,S are the wave propagation velocities in space and wavenumber, given by:

where k is the wavenumber, C_g is the group velocity; s is a coordinate in θ direction and m is a coordinate perpendicular to s; h is the mean water depth and K the wavenumber vector.

The left hand side of Eq. 8.16 corresponds to the kinematic terms, as derivatives for the propagation in space; and are the propagation velocities. The term with the derivative with respect to θ is the refraction term. The term with respect to σ causes a change of frequency. The right hand side is the source term and contains the effects of wind generation, whitecapping, dissipation, bottom friction, surf breaking, and nonlinear wave-wave interaction. This equation is implemented with finite difference schemes in all directions: time, geographic space, and spectral space.

The essential input data to run the model is the bathymetry for a sufficiently large area, the incident wave field, and the wind field. Various general and nested grids can be selected, depending on the availability of high-resolution data and the computational efficiency. Nesting is a very important implementation that can save computational time and increase accuracy. The model is validated with analytical solutions, field observations and experimental measurements, and has shown good agreement (Booij et al., 1999). Moreover, SWAN can operate with unstructured grids as well. Zijlema (2009) presented a method of vertex-based, fully implicit, and finite differences that is designed for unstructured meshes with high variability in geographic resolution. It is useful for complex bottom topographies in shallow areas and irregular shorelines.

SWAN is basically designed for applications in open coastal scale, with no-diffraction effects. That means that the model should be used in areas where variations in wave height are large within a horizontal scale of a few wavelengths.

SWAN organises its output in tables, maps (Figures 8.19 and 8.20) and time series, as well as 1D and 2D spectra, significant wave height and periods, average wave direction and directional spreading, one- and two-dimensional spectral source terms, root-mean-square of the orbital near-bottom motion, dissipation, wave induced force (based on the radiation-stress gradients), set-up, diffraction parameter, etc.

Mild slope equations models

MSE originally developed to describe the propagation of the waves over low gradient seabeds. MSE is commonly used in coastal engineering, since it can account well the effects of simultaneous diffraction and refraction of the waves due to coastlines or structures (Berkhoff, 1972). Mild-slope equations are a type of depth-averaged equation, within a x-y domain (2DH), applied in both deep and shallow waters for monochromatic waves (Lin, 2008).

The equations can be found in various forms, including the effects of wave breaking, nonlinearity of waves, wave-current interactions, and seabed friction. They calculate the wave amplitude or wave height but, if there is a constant water depth, the mild-slope equation reduces to the Helmholtz equation for wave diffraction. First introduced by Berkhoff (1972), the MSE assumed that the wave is linear and the slope is mild, obtaining the following main equation, improved by including the effects of friction dissipation and wave breaking:

where C is the wave celerity and C_g the group velocity; η̂ is the complex wave surface function; k is the wavenumber; σ is the wave frequency; w is a friction factor and γ is a wave breaking parameter. Friction is then obtained with:

where a is the wave amplitude, and f_r is a Reynolds dependent friction coefficient related to the bottom roughness; n is the Manning dissipation coefficient.

For weave breaking parameter, the following formulation is commonly used:

The original MSE has limitations because it is only applicable to linear waves and on mild bottom geometry. In addition, the equation does not contain energy dissipation, but in recent years there have been numerical advances to include energy dissipation and weakly non-linear waves with steeper bottom slopes. Mild-slope equation has been developed with different formulations that can be described by hyperbolic (Dingemans, 1997), elliptic (Berkhoff, 1972), and parabolic (Lin, 2008) formulation of the mild-slope equation respectively.

The practical application of wave transformation usually requires the simulation of directional random waves; thus, the principle of superposition of different wave frequency components can be applied. In general, MSE models for spectral wave conditions require inputs of the incoming directional random sea at the offshore boundary. The two-dimensional input spectra are discretized into a finite number of frequency and direction wave components. 

For the parabolic approach, the evolution of the amplitudes of all the wave components is computed simultaneously. Based on the calculations for all components and assuming a Rayleigh distribution, statistical quantities such as the significant wave height Hs can be calculated at every grid point. Figure 8.20 shows an example for a near-coast wave propagation obtained with a parabolic approximation of the mild slope equation for spectral wave conditions.

When wave reflection becomes relevant for wave propagation and transformation (i.e. within bays, harbours, sheltered areas, etc.), models should be based on the elliptical approximation of the mild-slope equation (Berkhoff, 1972; Madsen and Larsen, 1987; Tsay et al., 1989). This approach allows engineers to obtain the energetic response of reflected (totally or partially) waves, under the penetrating wave action.

Elliptic mild slope models solve the extended mild-slope equation to reproduce the main processes that control dynamics of waves when approaching coastal areas and entering into harbours (Figure 8.21): geometric refraction, shoaling, diffraction by obstacles, and full or partial reflection. Radiation conditions and free infinite outflow conditions are also available in the model. It also considers the complete spectral frequency distribution and the directional spreading of the wave energy spectrum. 

In addition to the above mechanisms, nonlinear waves may be simulated by incorporating amplitude-dependent wave dispersion, which has been demonstrated to be important in certain situations (Kirby and Dalrymple, 1983). 

This practical approach for harbour agitation and wave propagation can be assessed with the following, among others, commercial and non-commercial models: CGWAVE; ARTEMIS MIKE21; PHAROS, and MSP.

Phase resolving models (SWE, NSWE, and Boussinesq)

The Shallow Water Equations (SWE) are applied when water waves enter very shallow domains. Particles move basically horizontally and the vertical accelerations are negligible.

In this case, the propagation of the wave can be described by the SWE (Holthuijsen, 2007). These equations are derived from averaging the depth of the Navier-Stokes equations (NSE) assuming that the horizontal length scale is much greater than the vertical. The profile is uniform in depth and the vertical components very small. Using the conservation of mass, it can be shown that the vertical velocity is small, while using the momentum equation the vertical pressure gradients are hydrostatic. Therefore, the velocity profile is uniform in depth and the vertical components very small, and this is the reason for which SWE are also known as “long-wave equations”, given that they can be applied only to waves which are much larger to the bottom depth.

In the case of ignoring the Coriolis force, the frictional and viscous forces, the formulas of SWE are:

Equation 8.26 is derived from mass conservation and Eqq. 8.27 and 8.28 from momentum conservation, where η is the total fluid column height, (u,v) - a 2D vector - is the fluid’s horizontal velocity in the xy 2D domain.

To represent the ocean waves frequencies and physical behaviour, an improvement within the original SWE is needed, including the non-linearity terms and dispersive functions. The solution for this is the NSWE, as a non-hydrostatic waveflow solution model. It can be used for predicting transformation of dispersive surface waves from offshore to the beach, solving the surf zone and swash zone dynamics, wave propagation and agitation in bays and harbours, and rapidly varied shallow water flows typically found in coastal flooding (e.g. dike breaks, tsunamis and flood waves, density driven flows in coastal waters), as well as large-scale ocean circulation, tides and storm surges (typically solved by the original SWE models).

Main governing equation considers a 2DH wave motion over a domain represented in a Cartesian coordinate system (x,y). The depth-averaged, non-hydrostatic, free-surface flow can be described by the NSWE and comprise the conservation of mass and momentum. These equations are given by:

where t is the time; ζ is the free surface elevation, d is the water depth and h=d+ζ, u and v are depth-averaged flow velocities, q is the non-hydrostatic pressure, g the gravitational acceleration, Cf the bottom friction coefficient, and the group of τ are the horizontal turbulent stress terms.

The SWASH (Zijlema et al., 2011) is one of the latest worldwide available 🔗5 NSWE models. It is a numerical tool for simulating unsteady, non-hydrostatic, free-surface, rotational flow, and transport phenomena in coastal waters as driven by waves, tides, buoyancy, and wind forces. It provides a general basis for describing wave transformations from deep water to the beach, port or harbour, as well as complex changes to rapidly varied flows, and density driven flows in coastal seas, estuaries, lakes, and rivers. SWASH is an efficient and robust model that allows the application of a wide range of time and space scales of surface waves and shallow water flows in complex environments (Figure 8.22). The model can be also employed to resolve the dynamics of wave transformation, buoyancy flow, and turbulent exchange of momentum, salinity, heat, and suspended sediment in shallow seas, coastal waters, estuaries, reefs, rivers, and lakes.

SWASH may be run in depth-averaged mode or multi-layered mode in which the computational domain is divided into a fixed number of vertical terrain-following layers. SWASH improves its frequency dispersion by increasing the number of layers rather than increasing the order of derivatives of the dependent variables like Boussinesq-type wave models do.

BE can be applied as an alternative to NSWEs as the region between deep and shallow waters can be also well described by the Boussinesq model. In BE models, the horizontal component of the velocity is assumed to be constant in the water column and the vertical component of the velocity varies almost linearly over depth (2DH hypothesis). Essentially, these equations are the shallow-water equations with corrections for the vertical acceleration, and third order derivatives are the result of the Laplace equation forcing the vertical velocity of the velocity potential function to be expressed in terms of the horizontal velocity distribution. These equations can be readily expanded into two horizontal dimensions.

Researchers have introduced many different implementations of the Boussinesq equations, creating Boussinesq-type models to be applied for propagation in deep water and the process of wave-breaking (Brocchini, 2013). A vast majority of Boussinesq equations models (for fully non-linear approach) can be presented as follows:

with

where index of t denotes time; h is the equilibrium depth; η is the free-surface elevation, V is the horizontal velocity, and ∇ is the 2DH gradient operator. N and E respectively represent bottom drag and diffusion (artificial).

On the other hand, a similar family of equations exist and are applied in the region between deep and shallow waters; the Boussinesq equation-based model. For this approach, the main hypothesis is that the horizontal component of the velocity is assumed to be constant in the water column, and the vertical component of the velocity varies almost linearly over depth.

One of the most complete Boussinesq models, the fully nonlinear Boussinesq wave model (FUNWAVE) in its TVD version known as FUNWAVE-TVD model (Fengyan, et al., 2012), was developed at the Centre for Applied Coastal Research at the University of Delaware (USA). It includes several enhancements: i) a more complete set of fully nonlinear Boussinesq equations; ii) a MUSCLE-TVD finite volume scheme together with adaptive Runge Kutta time stepping; iii) shock-capturing wave breaking scheme, iv) wetting-drying moving boundary condition with HLL construction method for the scheme; and v) code parallelization using MPI method. The development of the FUNWAVE-TVD was prompted by the need to model tsunami waves in regional and coastal scale, coastal inundation, and wave propagation at basin scale (Figure 8.23). FUNWAVE is an open-source model available to the public 🔗6.

Numerical solutions of Boussinesq equations can be significantly corrupted if truncation errors, arising from the differencing of the leading order wave equation terms, are allowed to grow in size and become comparable to the terms describing the weak dispersion effects. All errors involved in solving the underlying nonlinear SWE are reduced to 4th order in grid spacing and time step size. Due to non-linear interaction in the model, higher harmonic waves will be generated as the program runs. These super harmonic waves could have very short wavelengths and the classic Boussinesq model is not valid. For this reason, a numerical filter suggested by Shapiro (1970) can be used.

In summary, both Boussinesq and NSWEs modelling approaches are the preferred solutions in their respective physical regions: Boussinesq where nonlinearity and dispersion are both significant, typically prior to breaking; and NSWE where nonlinearity predominates, from the mid-surf to inner surf zone shoreward, although it should be noted that there can be a significant overlap of these regions. Therefore, the NSWE models, which work well from the surf zone shoreward and naturally model wave breaking and the moving shoreline, find their main weakness in the absence of frequency dispersion, so that in deeper water waves will propagate incorrectly at the shallow water wave speed and, sooner or later, break again, which is not usual and correct in this region.

Two- and three-dimensions wave structure interaction model

CFD utilises numerical approaches to examine fluid flows, heat transfer and chemical reactions. Therefore, within wave propagation and structure interaction problems, the CFD term mainly refers to computer codes that solve the fully nonlinear Navier-Stokes equations in all three dimensions (3D).

CFD is then a state-of-the-art techniques for industrial and research applications, although its often high computation cost demands the use of high-performance computers. Within the wave propagation and wave structure interaction field, two of the most used CFD programs are: IH2VOF (two-dimension approach, derived from COBRAS original model) and OpenFOAM (three-dimension approach); these two codes are also well validated for many marine and ocean engineering applications.

As a classic Eulerian approach, both models are based on the RANS equations. These equations represent the continuum properties of the flow. By averaging the Navier-Stokes equations, more recent VARANS equations are obtained. The VARANS equations can have different terms, depending on the assumptions applied; for example, they include a k-ω turbulence model closure within the porous media, which make them the most suitable formulation for coastal engineering as the advantages of VARANS equations are numerous. The solving process yields very detailed solutions, both in time and space. Pressure and velocity fields are obtained cell-wise, even inside the porous zones, so that the whole three-dimensional flow structure is solved. Furthermore, non-linearity is inherent to the equations, and therefore all the complex interactions among the different processes are also taken into consideration. Finally, the effects of turbulence within the porous zones can be also easily incorporated with closure models.

IH2VOF model (Lara et al., 2006) solves 2D RANS equations for the oscillatory fluid and VARANS equations for the porous media. This 2D model can simulate the most relevant hydrodynamic near-field processes that take place in the interaction between waves and low-crested breakwaters. It considers wave reflection, transmission, overtopping, and breaking due to transient nonlinear waves, including turbulence in the fluid domain and in the permeable regions for any kind of geometry and number of layers. This model is highly validated, with different wave conditions and breakwater configurations, achieving a high degree of agreement with all the studied magnitudes, free surface displacement, pressure inside the porous structure, and velocity field.

IH2VOF is based on the decomposition of the instantaneous velocity and pressure fields into mean and turbulent components, the κ-ε equations for the turbulent kinetic energy κ, and its dissipation rate ε. This permits the simulation of any kind of coastal structure (e.g. rubble mound, vertical or mixed breakwaters). The free surface movement is tracked by the volume of fluid (VOF) method for one phase only, water and void. In order to replicate solid bodies immersed in the mesh instead of treating them as sawtooth shape, the model uses a cutting cell method. The main purpose of this technique is to use an orthogonal structured mesh in the simulations to save computational cost.

IH2VOF includes a complete set of wave generation boundary conditions, which cover most water depth ranges. These include a Dirichlet boundary condition and a moving boundary method, which are linked with an active wave absorption system to avoid an increase of the mean water level and the agitation. An internal source function can be also used to generate waves, but it has to be linked with a dissipation zone.

The RANS equations (clear fluid region) are redefined as follows:

Generally, IH2VOF application is within a detailed incident-wave and structure interaction (rubble-mound breakwaters, vertical structures and beaches), taking into account a realistic wave breaking and porous media interaction (see Figure 8.24).

The general VARANS equations include conservation of mass (8.39), conservation of momentum (8.40), and the VOF function advection equation (8.41) as follows:

where u is the extended averaged Darcy velocity; n is the porosity (volume of voids over the total volume); ρ is the density; p is the pressure; g is the acceleration of gravity; ν is the kinematic viscosity, and _α1 is the VOF function indicator (quantity of water per unit of volume at each cell).

The OpenFOAM (Higuera et al., 2014a and 2014b) is an extensive software package that has been widely used in industrial and academic applications. It is freely distributed 🔗7 as an open source CFD Toolbox, and it includes a broad range of features. IHFOAM 2.0 is an extension of the original software for coastal applications, newly developed with a three-dimensional numerical two-phase flow solver, specially designed to simulate coastal, offshore, and hydraulic engineering processes. It contains an advanced multiphysics model, widely used in the industry. A wide collection of boundary conditions, which handle wave generation and active absorption at the boundaries with a high practical application to coastal and harbour engineering (Figure 8.25), makes IHFOAM 2.0 different from the rest of solvers. Maza et al. (2016) have studied and proposed natural-based solutions for coastal protections using IHFOAM.

Various discretization methods are used in water wave solving problems, a brief description for each of them is presented below (for additional references see Sections 5.4.2.4 and 7.2.3.5):

• FDM. Maybe the most used and simplest ways to solve numerically partial differential equations (PDEs). The method establishes the value of the flow variable at a given point based on the number of neighbour points. The numerical domain forms a grid. The governing equations of the fluid are considered in their differential form at each point in the domain, so that the solution is solved by replacing the partial derivatives with approximations by means of the nodal values of the functions. This method is recommended for structured grids and low-order equation schemes.
• FEM for the solution of PDEs employs variational methods to minimise the error of the approximated solution, similarly to the Galerkin method. FEM was used in structural mechanics but this technique developed for computational fluid dynamics applications being introduced to common wave propagation and agitations models. FEM technique, similarly to the FDM, is based on the concept of subdividing a continuum computational domain into elements, forming a grid of triangular or quadrilateral unstructured elements or curved cells (Figure 8.26). Therefore, the method can handle problems with great geometric complexity, such as harbour perimeter definition, concentration of nodes at relevant parts of the domain, etc. 

FEM used variational methods, which in practice means that the solution is assumed to have a prescribed form and to belong to a function space. The function space is built by varying functions, such as linear and quadratic. The varying functions connect the nodal points, which can be the vertices, mid-side points, mid-element points, etc., of the elements. As a result, the geometric representation of the domain plays a crucial role in the outcome of the numerical simulation. The original PDEs are not solved by the FEM. Instead, the solution is approximated locally by an integral form of the PDEs. The integral of the inner product of the residual and the weight functions are constructed. The integral is set to zero and trial functions are used to minimise the residual. The most general integral form is obtained from a weighted residual formulation. The process eliminates all the spatial derivatives from the PDEs and, therefore, differential type boundary conditions for transient problems and algebraic type boundary conditions for steady state problems can be considered, hence the differential equations become algebraic. Only one equation is solved per grid node, which has one variable as unknown. The same variable is also unknown at the neighbouring cells.

• FVM solves PDEs by transforming them to algebraic equations around a control volume (subdivisions of the computational domain). The variables are calculated at the centre of each control volume. General interpolation methods are used to derive the values of the variables at the surfaces of the control volume considering the neighbour control volumes as well. The FVM has two major advantages: i) it is able to accommodate any type of grid, making it applicable for domains of high complexity; and ii) it is conservative by definition, since the control volumes that share a boundary have the same surface integrals, describing the convective and diffusive fluxes. FVMs are very popular in the numerical wave propagation community, succeeding in free surface flow simulations, especially when highly nonlinear processes are involved, such as wave breaking