Advances in Numerical Reynolds-Averaged Navier–Stokes Modelling of Wave-Structure-Seabed Interactions and Scour

: This review paper presents the recent advances in the numerical modelling of wave– structure–seabed interactions. The processes that are involved in wave–structure interactions, which leads to sediment transport and scour effects, are summarized. Subsequently, the three most common approaches for modelling sediment transport that is induced by wave–structure interactions are described. The applicability of each numerical approach is also included with a summary of the most recent studies. These approaches are based on the Reynolds-Averaged Navier–Stokes (RANS) equations for the ﬂuid phase, and mostly differ in how they tackle the seabed response. Finally, future prospects of research are discussed.


Introduction
Coastal structures, such as seawalls and breakwaters, provide protection from the effects of the sea and create the conditions for economic growth in coastal environments. These maritime structures also protect ports and coasts against sea dynamics, and their function will be even more important in the upcoming years due to sea level rise and coastal regression as a result of global warming [1]. Regarding maritime structures and nearby environments, the effects of sea level rise will include: increasing levels of sediment transport and erosion, higher impacts of wave energy on structures, the loss of stability, and increasing overtopping events [2]. On the other hand, increased maritime traffic and human activities along coasts make the proper study and modelling of wave-structureseabed interactions indispensable in achieving a balance between human intervention, costs, and environmental impacts.
The design of maritime structures must satisfy the project requirements against windwave actions [3,4] during the structure service life. The uncertainty in maritime structure design increases when its environment is also considered, especially the seabed, which is the foundation of the breakwater. Some structure failure examples due to seabed failure have been reported in the literature [5,6]. Most of these failures were caused by progressive sediment transport and, thus, erosion of the seabed, i.e., by the scouring process of the seabed. In fact, historical studies reveal that scour is one of the main reasons behind the failure of coastal maritime structures, such as seawalls, groynes, breakwaters, revetments, and artificial reefs [7,8].
Sediment transport has been the subject of extensive research over recent decades, particularly wave-structure interactions and, more importantly, their effects on the seabed. It includes mechanical models [9,10], analytical solutions of wave-structure-seabed interactions [11][12][13], and their verification with numerical models and experimental tests.
Experimentally, both field and laboratory measurements focus on the structure foundation. Physical models of laboratory experiments face very serious problems with sediment and wave scale, even when the prototype with the Froude scale in a laboratory is well achieved. Despite the problems and issues that still arise in the experimental testing of the seabed response, several studies have investigated sediment transport and the scour around maritime structures. Sutherland et al. [14], Pearce et al. [15], Tsai et al. [16], Jayaratne et al. [17] studied scour around vertical and sloping seawalls, and Sutherland et al. [18,19], Sumer and Fredsøe [20], Fausset [21], Temel and Dogan [22] tested the seabed response around rubble-mound breakwaters.
Numerical models appear to be a valuable tool in designing protective structures and performing parametric studies because of the prohibitive cost of field and laboratory experiments. Numerical modelling requires an integral method that incorporates the processes and effects (displacement, stresses, pressures, etc.) between the wave, seabed, and maritime structure. Few numerical models address sediment transport at scales of practical engineering interest or consider multiphase effects and coupled interaction responses. In addition, there is no real consensus regarding an ideal modelling approach. This is usually attributed to either the lack of understanding of the underlying physics or to the fact that multi-scale physics become too computationally and time expensive when dealing with problems at a practical scale. The latter not only delays the proper characterization of wave-structure-seabed interactions, but also makes it difficult to verify laboratory tests and extrapolate it in the structure design.
Several review papers on sediment transport and seabed response have been published; the topics include: (i) analytical formulas for calculating the scour depth around piers [23,24], (ii) analytical and experimental studies of scour around coastal structures [25,26], and (iii) the analysis of the seabed response under wave action [27]. However, given the great progress made in numerical modelling in the last decade, there is a need to put the most widely used numerical approaches for modelling the wave-structure-seabed interactions into context and focus. The numerical model selected and the resolution method mainly depend on the coastal study area, the importance and dependence of the processes involved, the characteristic temporal and spatial scales, and the required computational time. For example, for the swash zone, situated at the landward edge of the inundated part of a coastal region, or for numerical simulations of near-shore regions (where the water depth is much smaller than the wavelength, h/L << 1), the use of depth-averaged models that are based on the Shallow Water Equations (SWE) [28] allow for modelling short wave events and studying the morphological evolution of the coastal zone [29], and the scour and deposition around structures [30,31]. The Shallow Water Equations may be used in two different approaches [32]: those that assume that the pressure is hydrostatic and those that do not, which leads to formulations, such as Boussinesq equations or Non-hydrostatic shallow water equations [33]. These depth-averaged models provide practical use for engineering purpose, since they offer acceptable accuracy with very low computational effort. When the depth profiles and turbulence effects are important, more complex and accurate models can be used based on the solution of the statistically Reynolds-Averaged Navier-Stokes equations (RANS), which provide the mean flow field through the domain as well as some turbulent quantities (turbulence kinetic energy and its dissipation rate in most of the models). Large-Eddy Simulation (LES) is another approach, which resolves a filtered version of the Navier-Stokes equations, and it provides significant accuracy as coherent 3D hydrodynamic structures that are larger than the grid scale are resolved. However, its computational cost is still too high for most practical problems. In the seabed response and sediment transport around maritime structures, the models that are based on RANS equations give an adequate description of the wave-structure interactions, as well as the turbulent flow and bottom boundary layer.
Hence, this paper presents a comprehensive review of the numerical approaches for modelling the seabed response, in particular the scour, around maritime structures that are based on RANS equations for modelling the wave-structure interactions. The review focuses, for the most part, on numerical model papers that were published in the last 10 years and the main keywords used for this review are: scour, sediment transport, numerical modelling, breakwaters, piles, RANS equations, Biot equations, multiphase theory, and fluid-structure-seabed-interactions (FSSI).
This paper is structured, as follows: Section 2 gathers definitions and physical processes that are involved in the wave-structure interactions that contribute to scour. Section 3 describes the general methodology of each numerical approach. This Section also discusses the applicability, performances, and comparisons between the numerical approaches. Finally, Section 4 presents the closing comments and future scopes of research in this area.

Damage to Structures and Coastline: The Role of Scour
This section includes a description of the seabed response, scour, and the consequences on maritime structures and coastlines. In addition, the wave-structure interactions processes that induce sediment transport around structures are discussed.

Concept of Scour
The presence of the protection structure changes the flow patterns in its immediate vicinity, inducing wave reflection and breaking, turbulence, and liquefaction [34]. These processes increase the fluid shear stress at the bottom (τ f ,b ) and cause transport when a critical value (τ cr,s ) is reached. Increased local sediment transport leads to one of the greatest problems for coastal regions and marine structures: scour [35,36]. Scour around coastal structures is the additional erosion that takes place due to hydrodynamic forces acting on the seabed (Figure 1). Scour has significant impacts on the economic and service life of many coastal structures [21]. Figure 2 shows an example of scour-induced damage at sloping structures. When scour undermines the structure, it cannot support the armour layer of a breakwater, which then slides and loses the fill material. Scour also impacts vertical-front caissons of breakwaters and other gravity structures ( Figure 2). Seawalls can also settle and collapse as a result of scour. More details on scour around such structures can be found in the textbooks of Whitehouse [36] and Sumer and Fredsoe [34].

Mound breakwater
Vertical caisson Damage to structures is not the only consequence of seabed scour. In fact, the presence of coastal structures can accelerate the erosion (and deposition), thus degrading coastlines [37][38][39], as shown in Figure 1. These examples of structural collapse or a loss of material along the coastline result from both transformation and interaction processes between marine dynamics, especially waves, and the coastal structure.

Processes Involved in Wave-Structure-Seabed Interactions
It is necessary to identify the physical processes contributing to the scour in order to properly design maritime protection structures and predict the final response of the seabed. A state-of-the-art review on scour processes has been presented in Sutherland et al. [18], Whitehouse [35], Sumer and Fredsøe [40]. The reviewed works complement and considerably extend the understanding of the governing processes for certain structures. According to these references, the processes that lead to scour are (acting singularly or in combination): (A) Wave energy transformation by its interactions with the structure [3,41,42]: Wave reflection: sediment transport due to wave reflection is perhaps the most commonly cited process in structure-seabed interactions [43]. The presence of a maritime structure transforms the incident wave energy into reflected energy, which is returned to the sea ( Figure 3). This process involves an increase in wave height and wave energy, in front of the structure, i.e. greater shear stresses on the seabed and, thus, sediment transport. • Wave dissipation: when waves impact a coastal structure, part of its incident wave energy is dissipated by (Figure 3): (a) wave-breaking on the slope or wall, (b) turbulent interactions (circulation and friction) with the armour layer, and (c) wave propagation through the porous media of the structure [44]. Among these, wave-breaking on the wall or slope mainly affects sediment transport, which generates turbulence in front of the structure and then mobilizes sediment.
(B) Wave diffraction and refraction: when waves interact with an obstacle or coastline, part of the incident energy is concentrated or dispersed in a certain direction, which induces variations in the wave height and flow velocity, causing non-uniform erosion or sediment deposition [45,46]. Furthermore, waves that interact with submerged breakwaters are reduced by a number of energy dissipation mechanisms, including wave breaking and frictional dissipation [47,48], which can cause sediment transport around these structures. (C) Turbulent flow around the structure: the flow contraction and velocity field generate turbulence and vortexes around and in front of the structure (Figure 4), which increase the bottom shear stresses (τ f ,b ) and sediment mobilization [34,43]. Some figures of Higuera et al. [49] show the increase of turbulence kinetic energy (k = u 2 +v 2 +w 2 2 ) in front of and around a vertical wall and a porous breakwater. (D) Liquefaction: the pressure differential in the soil may produce the liquefaction of the seabed [13,50,51], which leads to the loss of the seabed load-bearing capacity.
A liquefied soil will be more susceptible to scour. The liquefaction phenomenon happens when the soil effective stresses are cancelled due to an increase of the pore water pressure. Figure 5 shows the stress field and pore pressure in front of a maritime structure. It can be seen how the effective normal stress (σ s ) decreases as the waveinduced pore pressure (p s ) increases.  ) increases at the toe of the seawall. Both of the figures are adapted from the numerical simulations to study the interactions between ocean waves and maritime structures performed by Ye et al. [53] and Higuera [54], respectively.

Numerical Modelling of Seabed Response
This section presents the governing equations and general methodology of three main Eulerian approaches, in this paper called (i) simple approach (Section 3.1), (ii) Biot approach (Section 3.2), and (iii) full multiphase approach (Section 3.3), for modelling the sediment transport and, thus, the seabed response around maritime structures. The numerical approaches attempt to capture the processes that are involved in the interactions of waves with the structure and their consequences on the seabed foundation. Figure 6 shows a general outline of the numerical approaches described hereafter, which are based on the Reynolds-Averaged Navier-Stokes (RANS) equations (the mass and momentum conservation equations) for the wave-structure interactions and mainly differ in the way that they model the seabed behaviour. This section also shows the applicability and summarizes the performance and parameters used for each numerical methods.

Simple Approach
The simple approach decomposes the wave-structure-seabed problem into three modules: (1) wave propagation and transformation with the structure, (2) bed load and suspended sediment transport, and (3) bed morphological changes. Figure 7 shows a summary diagram of the modules, parameters, and the interactions between modules. A detailed description of the modules is given below.

Wave-Structure Module: Rans Equations
The incompressible Reynolds-Averaged Navier-Stokes (RANS) equations are applied as the governing equations of fluid flow. These are time-averaged conservation equations of mass and momentum, which describe the mean flow field, respectively: where the sub-indexes i, j = 1, 2, and 3 represent the x-, y-, and z-directions, respectively; f represents the fluid-phase; u f = (u f ,x , u f ,y , u f ,z ) is the fluid velocity; p f is the pressure; ρ f is the fluid density; ν = µ/ρ f is the fluid kinematic viscosity (being µ the dynamic viscosity); ν T is the turbulent/eddy viscosity; and g i = (0, 0, −g) is the gravitational acceleration. Equations (1) and (2) have the following considerations: • The resulting flow variables (e.g., u f , p f ) describe the mean flow field. Turbulence effects are introduced via the Reynolds's stresses in the momentum equation, which must be modelled using a turbulence closure. • The turbulence models are usually based on the Boussinesq's eddy viscosity hy- . The eddy viscosity ν T is estimated using semi-empirical models. • The free-surface of the waves and flow inside the porous media of the structure should be considered in this module. • Density variations are considered in some instances. In this case, an extra densityweighted averaging step is carried out, which results in the Favre-Averaged Navier-Stokes equations (FANS).

Turbulence Model
This review paper does not focus on the different turbulence models and Blondeaux et al. [55] provides a detailed description of the flow-structure by comparing various turbulence models. However, it should be noted that, in the study of wave-structureseabed interactions, the most appropriate turbulence model should provide an accurate description of (i) the flow close to the bottom (seabed surface) to obtain a reliable evaluation of the sediments transported by the flow [56] and (ii) the turbulent flow generated by the interactions with the structure. The turbulence models used in the papers cited in this work are summarized in Sections 3.1.4, 3.2.2, and 3.3.2, respectively.

Free-Surface Capturing
Capturing the free surface is essential for modelling wave motion and wave-breaking. The Volume Of Fluid (VOF) technique is one of the most common methods to track the free surface location [54,57]. This method determines the phase within each cell in the domain by introducing an indicator function, φ (φ = 1 full of water, φ = 0 full of air, interface 0 < φ < 1), and solves the momentum equation (Equation (2)) for a pseudo fluid whose properties are weighted sums based on φ. The distribution of φ is determined by solving a simple VOF advection equation: An alternative interface capturing method is the Level Set Method (LSM) that was proposed by Osher and Sethian [58], in which the interface is implicitly represented by the zero level set of a smooth signed distance function that adopts negative values in a phase, zero in the interphase, and positive values in the other phase. The evolution of this function is captured using a conservation equation that is similar to Equation (3), but an interface thickness is predefined. Across this region, a Heaviside function is used to smooth sharp flow variable changes and ensure numerical stability. In comparison with the VOF method, the LSM is more accurate in the capture of high curvature interfaces. However, it is not capable of ensuring mass conservation, which is a major drawback in comparison with the VOF method [59].

Flow Inside the Permeable Structure
If the coastal structure is permeable, i.e., formed by a porous material, then the numerical approach should calculate the flow inside the porous media. In this case, the RANS equations need an additional volume-based averaging step to account for low porosity materials, n ∈ (0.35-0.65) (being n, the porosity of the permeable structure), such as those that are normally found in coastal engineering structures. Hence, RANS equations are integrated into a porous control volume, obtaining the Volume-Averaged Reynolds-Averaged Navier-Stokes equations (VARANS) [49,60]: Equations (4)-(6) are the mass, momentum, and phase volume fraction conservation equations inside the porous media, respectively. After volume averaging, three new terms appear on the right hand side of Equation (5), which are closure terms to account for the frictional and pressure forces, as well as the added mass of the individual components of the porous media. These three terms are drag forces that characterize the flow through the porous media [61]: the linear term (laminar flow), au f ,i ; the quadratic term (turbulent flow), bu f ,i |u f ,i |; and, the inertial acceleration (unsteady flow) term c ∂u f ,i ∂t . The coefficients a, b, and c depend on the physical properties of the material and control the balance between each of the friction terms [62,63].
Note that, for impermeable structures, a non-slip (u f = 0) boundary condition is imposed on the structure surface and the governing equations in the wave-structure module in the case of RANS equations.

Sediment Transport Module
This module includes the calculation of the bed load and suspended sediment transport. Numerous bed load equations have been developed over the past century, q B . The formulas are often empirical and they are commonly based on excess shear stress, where the shear stress is greater than the critical value for incipient motion (τ b, f > τ cr,s ). Bedload equations can include parameters for the characteristics of the fluid and the seabed, such as bed-forms (flat bed or ripples), slope, and sediment composition [64][65][66]. Common variables include: with τ f ,b the bottom shear stress of the fluid (from the wave-structure module); τ cr,s ∼ d 50 (ρ s − ρ f )g, the critical bed shear stress; and, d 50 and ρ s , the characteristic diameter and density of the sediment, respectively.
For the suspended sediment transport, the suspended concentration (c s ) is often calculated by the advection-diffusion equation (Equation (8)) [67,68], but other papers use empirical formulas for c s as Ahmad et al. [69], who used the empirical formula of Rouse [70]; and, others used empirical formulas for the total sediment transport load (bed and suspended), such as Tofany et al. [71,72].
where w s is the settling velocity of the sediment, whose expression is empirical and depends on the hindering settling effects. D T is the diffusion coefficient, normally related with the eddy viscosity, ν T (from the wave-structure module) [56,73], and u f is the flow velocity (from the wave-structure module). Concerning the suspended concentration, the suspended load transport is calculated by:

Morphological Bed Module
The evolution of the seabed profile is obtained with the Exner formula. It is based on the sediment local mass balance and it involves a non-linear propagation of the bed-level deformation in the direction of the sediment transport and the spacial variation of the sediment fluxes between the bed load and suspended load [69]. The equation is defined as: where q t = q B + q S is the total load transport; n s is the seabed porosity; and, the term (E − D) defines the net sediment movement of the suspended load [74]. Here, E and D are the erosion and deposition rates, whose expressions depend on the settling velocity, w s , and the volumetric reference concentration at the bed (c b ).

Applicability of the Simple Approach
This numerical method is referred to as the simple approach in this paper, because the seabed behaviour is modelled using empirical equations, which are relatively easier to implement in comparison with the other approaches. The simple approach is used, in particular, for estimating the scour around structures with low weight-structure load on the seabed, such as submarine pipelines. Some of the most recent and outstanding studies that used the simple approach for estimating the scour around structures are: Ahmad et al. [69], who investigated the local scour around a composition of pile structures (Figure 8), and, Ahmad et al. [75], who presented the numerical approach for studying scour due to wave impacts on a vertical wall. They showed that the wave-breaking and impact on the structure govern the different stages of the scour formation. Baykal et al. [76] also investigated the flow and scour around a vertical pile. Seabed deformation was modelled while using the mesh deformation technique, which allowed for showing the development of turbulent structures in the wake of the monopile, such as horseshoe vortices. Furthermore, Chen [77] developed a Lagrangian-Eulerian model based on the two-dimensional (2D) Navier-Stokes equations for the wave field, and a bed load sediment transport and morphological model to simulate the scouring process in front of a verticalwall. Their results tend to underestimate the scouring rate. Gislason et al. [78] calculated the 2D flow and sediment transport in front of a breakwater, coupling a morphological model with Navier-Stokes equations. Their results do not show the seabed movement in detail, a drawback of the simple modelling approach. Li et al. [79] investigated the 2D local scour beneath two submarine pipelines in tandem based on the calculation of the bed/suspended load transport and the morphological changes of the seabed. This work investigated the scour around the pipelines under wave-current conditions. Omara et al. [80] simulated the scour processes, both hydrodynamically and morphologically, around vertical and inclined piers in a new version of FLOW-3D that estimates the motion of sediment transport (suspended, settling, entrainment, and bed load). Tofany et al. [71,72] studied the effects of breakwater steepness and wave overtopping on the scouring process in front of impermeable and vertical breakwaters. They used a 2D RANS-VOF model with additional bottom shear stresses in the momentum equations, coupled with turbulence closure, sediment transport load, and morphological model. Klonaris et al. [81] presented the numerical and experimental results of the evolution of a sandy beach in a surf-zone around a permeable submerged breakwater. They used the Boussinesq-type equations instead of Navier-Stokes equations to simulate wave propagation over a porous bed, coupled with a sediment transport module, including both bed and suspended loads.  Table 1 provides a summary of the parameters and performance of the simple model for some of the relevant studies cited. In general, these focus on water-seabed interactions. Thus, they often disregard the air flow. The RANS equations for the liquid phase are solved using the k − ω model of Wilcox et al. [83] as turbulence closure model, which is considered to be a better choice for solving the boundary layer over the other two-equation turbulence models. Wave action is often imposed using source terms on the corresponding wave generating boundary. In some studies, the seabed deformation is accounted for using an interface tracking method [69,69]. The Level Set Method (LSM) [58] is the most common method in this regard, which is considered to be more accurate and stable when dealing with high curvature interfaces than the VOF method. Note that Baykal et al. [76] used a mesh deformation technique, which might be more accurate, but it is noticeably more expensive from a computational perspective. Seabed deformation (either using an interface tracking method or mesh deformation) allows for two-way coupling between the fluid region and seabed. This is often not the case in other more complex approaches, such as the RANS-Biot combination. In general, although the results of these studies are considered to be acceptable, they tend to underestimate the scouring rate and do not show the seabed movement in detail. Moreover, there is no consensus regarding the choice of empirical formulas for the seabed module.

Biot Approach
An alternative to the simple method for describing wave--seabed interactions around marine structures is to combine Biot theory [85] (see Equations (11)-(13)) with the Navier-Stokes (NS) equations for the wave field. This numerical approach studies the wavestructure-seabed interactions with two modules: (1) the wave-structure module, based on RANS or VARANS equations as with the simple approach (see Section 3.1.1), and (2) the seabed module (described below) in which the Biot's equations govern the dynamic response of the porous seabed under wave loading. Figure 9 presents a summary of the modules that are presented in this numerical approach.

Seabed Module: Biot's Equations
The Biot's equations are adopted as the governing equations for describing the dynamic response of the isotropic porous seabed under wave loading [85]. This paper provides a brief outline of the Biot's equations, but further details can be found in Jeng and Ou [86].
(1) The momentum conservation equation for the mixture (fluid and sediment) can be written as: (2) When considering the flow to be governed by Darcy's law, the momentum equation of fluid phase is written as: (3) The mass conservation equation for fluid phase can be stated as; In these equations, σ s and τ s are the effective normal and shear stresses of the soil, respectively; p s is the wave-induced oscillatory pore pressure; x s = (x s,x , x s,y , x s,z ) is the seabed displacement; w f is the averaged (Darcy) fluid velocity; k s and n s are the permeability and porosity of the seabed, respectively; ρ = ρ f n s + ρ s (1 − n s ) is the total density; and, K f is the bulk modulus of pore water that is related with the degree of saturation of the seabed, S r . For this module, there are two considerations to be taken into account: (i) the inertial terms and (ii) the rheological model that relates the stresses and deformation of the seabed under wave action.

Inertial Terms
Different expressions for the Biot's equations are implemented, depending on the rate of wave loading (wave period, T w ) and the characteristics of the porous media (permeability, k s ) [87][88][89]: • Fully dynamic: the coupled equations of flow and deformation are formulated to include both acceleration terms: Partly dynamic, which is also known as u − P dynamic: the coupled equations of flow and deformation are formulated when only considering the acceleration of the seabed, Quasi-static: both inertial terms are ignored, resulting in quasi-static coupled flow and deformation formulation.

Rheological Models
The response of the wave-seabed interactions can be evaluated by rheological models, which relate the normal and shear stresses, (σ s , τ s ), and the strain, s = ∂x s,i ∂x i , of the soil under the cyclic wave action. Most of the works analyse the wave-structure interactions on non-cohesive seabeds, i.e. composed by sands and/or gravels. For cohesive seabeds (silts and clays), the lack of laboratory or field tests makes it difficult to calibrate and validate the numerical models. A detailed description of rheological models for cohesive and non-cohesive soils is given in Díaz-Carrasco [90].
Most of the papers study the seabed response that is composed of non-cohesive sediments that behave primarily as elastic materials that, in contact with water, have a pore-elastic (dense seabed) [91,92] or pore-elastoplastic (loose seabed) behaviour [86,[93][94][95].

Interactions between the modules
The steps to interact the results of the wave-structure module with the seabed module are the following: 1.
The wave-induced pressure and bottom shear stresses determined in the wave module are the boundary conditions for the seabed module: p s = p f and τ f ,b = τ s at the seabed surface.

2.
The initial consolidation state of the seabed due to static wave and structure loading has to be determined before the dynamic wave loading is applied in the numerical model. This initial consolidation state of the seabed will be the initial stress state for the dynamic seabed response under wave loading.

3.
The seabed's feedback effect on the wave and structure motions (two-way coupling) is neglected.

Applicability of the Biot Approach
There are many studies that use this approach for estimating the seabed movement around maritime structures, in particular breakwaters. Their results focus more on the mechanical behaviour (stress field and pressure) of the bearing seabed under wave loading and its interactions with the structure. Jeng et al. [91], Zhang et al. [96,97], Zhao et al. [98] studied the response of the porous seabed around a permeable breakwater while using the VARANS equations that were coupled with the dynamic Biot equations for a porous elastic seabed, and Zhang et al. [99] around a submarine pipeline. He et al. [100], Ye et al. [101] used the integrated numerical model FSSI-CAS 2D with VARANS equations and the dynamic Biot equations that were developed by Ye et al. [102] for studying the harbour zone of Yantai port in China as an engineering case. Zhao and Jeng [92] developed an integrated model (VARANS equations coupled with Biot equations) to investigate the wave-induced sloping seabed response in the vicinity of breakwaters. Their numerical results showed that the breakwater permeability, density, and slope of the seabed significantly affect the potential liquefaction, which leads to scour effects ( Figure 10). Ye et al. [93], Zhao et al. [95] used the same approach to study the wave-structure-seabed interactions. They modelled the relation between stresses and deformation of the seabed with a poro-elastoplastic behaviour. Lin et al. [103,104], Zhao et al. [105] investigated the seabed response around a mono-pile foundation simulating the wave field with two-phase incompressible RANS equations and a wave-induced dynamic seabed response with a quasi-static version of the Biot equations. More recently, Jeng et al. [106] developed a combined approach that solved the Biot equations using a Radial Basis Function method to investigate the wave-induced soil response around a submerged breakwater. This method is meshless, which simplifies the solution of the seabed region and allows its application to larger configurations, such as around offshore pipelines Wang et al. [107]. Li et al. [84] developed an open-source numerical toolbox for modelling the porous seabed interactions with waves and structures in the finite-volume-method (FVM)-based OpenFOAM framework. This toolbox incorporates the incompressible Navier-Stokes equations for the wave module and the Biot equations for the seabed module considering the anisotropic seabed characteristics. Li et al. [108] studied the wave-induced seabed response and liquefaction risk around a hexagonal gravity-based offshore foundation. As can be seen, the abundance of works regarding the application of the RANS-Biot approach demonstrates its viability for estimating the scour around relatively large structures. Although most of them limit impermeable structures and its results focus on the stress state and pressure on the seabed. In fact, these works generally study the liquefaction phenomenon that is induced by wave action and how it affects the structure.  Table 2 summarizes the performance and parameters for developing the RANS-Biot approach of the most relevant studies found. These are noticeably more common in the literature than the other two approaches, which shows its good compromise with regards to accuracy versus complexity. Most of these studies consider the effects of wave breaking and wave-structure interactions, for which they use the VARANS equations in combination with the VOF method. Moreover, the standard k − model is often chosen for turbulence closure. This approach does not consider seabed effects on the fluid phase (oneway coupling), as mentioned above. Most of these studies use the u − p approximation for the inertial terms in the Biot's equations, which offers improved results over the quasi-static formulation [88]. Elsafti and Oumeraci [109] observed that this approximation significantly reduces the computational time, and it should be considered whenever pore fluid acceleration is negligible. Otherwise, the fully-dynamic formulation is recommended.

Full Multiphase Approach
The movement of sediments due to the fluid flow is a two-phase phenomenon [110]. The starting point for developing any multiphase model is the local instantaneous formulation, which consists of field equations that are applied to the distinct phases complemented with jump conditions to match the solution at the interface and source terms to account for interactions between phases.

Governing Equations
The governing equations for multiphase models are the conservation equations of mass and momentum for each phase. Generally, the system is simplified using a pseudofluid approach for the water-air region in combination with a free surface capture method (see Section 3.1.1). A specific momentum equation is used for the sediment (solid phase). However, specific continuity equations are necessary for each phase since the system now includes three phases. The final mass and momentum conservation equations are written as, (1) Mass conservation equations: (2) Momentum conservation equations: where the sub-indexes f and s refer to the fluid and sediment phases, respectively; ρ is the density; c is the volume concentration of sediment; u is the velocity; ν is the kinematic viscosity; S c is the Schmidt number; p is the pressure; T represents the viscous and turbulent stress tensor; and, τ p = ρ s ρ f −ρ s (1 − w s ) 2 g is the response time of the particles, with w s being the settling velocity of the sediment. The last two terms in the momentum equations are related to momentum interchanges between the sediment and fluid phases: the first term accounts for the drag force and the second term for turbulent dispersion.

Applicability of the Full Multiphase Approach
The application of the full multiphase model allows for greater accuracy regarding sediment movement and seabed-wave-structure response. In this case, the seabed region is modelled as a continuous medium, which provides greater detail in the description of pressure, stresses, and interchange terms between phases. This numerical approach is the most recent, but it has a high numerical complexity. For the high sediment concentration region (near the bed), the model requires a very refined mesh, which leads to a higher computational cost. The studies that are based on this approach mainly focus on smallscale problems that are related with the estimation of the sediment transport around small structures, such as monopiles.
Many models have been developed for the two phases, i.e. water-sediment [111][112][113][114][115][116], and few of them include air as a third phase to handle the free surface tracking. These studies mostly differ in the way they treat closure terms for particle stresses, turbulence, and momentum exchanges. In this regard, Hsu et al. [111] show that adding a turbulence modulation term due to particle presence to the k − model improves the accuracy in the boundary layer. Given the grid requirements and complexity of the models, these are usually limited to 1D or 2D domains. Lee et al. [115] propose the inclusion of extra rheological terms that account for enduring-contact, as well as viscous and interstitial effects, which extend the range of application of their model to higher Reynolds numbers. Lee et al. [110] and Ouda and Toorman [117] used a three-phase model to study the local scour caused by a submerged wall jet and around a submarine pipeline. In both studies, the numerical results exhibited good agreement with the experimental data and the model was able to capture the key processes that govern the scour effect ( Figure 11).  Table 3 proposes a summary of the performance and parameters of the most relevant studies using this approach. The multiphase approach requires turbulence modulation due to the presence of sediment particles. Hence, these studies generally use a modified k − model. Moreover, momentum exchanges between phases are accounted for by using source terms in the momentum equations. Thus, this approach offers two-way coupling between the sediments and the water phase. Only one study considered a domain big enough to also include the air-water interface [117]. However, in this study, the authors simplified the system by using a mixture model, wherein a single system of continuity and momentum equations is used and the phase locations are tracked using a concentration marker. As can be seen, the use of a full multiphase approach requires the consideration of several closure terms for the particle stresses, turbulence stresses, and phase interaction mechanisms, most of which are based on empirical equations. Table 3. Summary of the most recent and leading works that use the full multiphase approach to model the wave-structureseabed interactions.

Author
Year

Conclusion and Future Prospects
The main objective of this paper is to present a comprehensive review of the numerical approaches for modelling the wave-structure-seabed interactions, while focusing on the seabed response, i.e., scour. This review is intended to complement the review papers on sediment transport and seabed response around maritime structures, which are mainly based on analytical and experimental studies. For this purpose, a literature review of the main and most recent works (over the last 10 years) on numerical modelling was carried out.
There are three main Eulerian approaches for modelling the sediment transport and the seabed response around maritime structures. In this paper, they are referred to as: (i) simple approach, (ii) Biot approach, and (ii) full multiphase approach. These three approaches are based on Reynolds-Averaged Navier-Stokes equations for the wave-structure module and they differ mostly in how they tackle the seabed module.

1.
The simple approach considers: (a) the bed load transport with empirical formulas; (b) the suspended concentration with the advection-diffusion equation and suspended load transport; and, (c) the seabed level with the Exner formula. The wave module is mainly based on RANS or VARANS equations, depending on whether the structure is permeable or impermeable. The studies that use the simple modelling approach mainly focus on small and low weight-structures to simulate the movement and deformation of the seabed (two-coupled approach). For large structures, such as breakwaters, these studies tend to underestimate the scouring results.

2.
The Biot approach involves the RANS or VARANS equations for the wave module and the Biot's equations for the seabed module. Biot's equations are used to obtain the seabed displacement and stresses. The wave model is responsible for the wave generation and propagation, and it determines the pressure, p f , and stresses, τ f ,b , acting on th seabed and marine structures. The outputs of the wave module are used as boundary conditions for the seabed module. The studies that used Biot's equations for the seabed module focus more on the stress field and pressures, which is, the mechanical behaviour of the seabed.

3.
The full multiphase approach is the more complex approach and is difficult to implement numerically. It solves each phase (sediment, water, and air) with the Navier-Stokes equations (mass and momentum) at the same time and space. Most of the studies are based on small-scale problems given the very high computational costs that are associated with its use.
This paper describes the general concept and methodology of each numerical approach for modelling wave-structure-seabed interactions. However, there are certain aspects when developing a numerical model that need to be studied in more detail and that were not the focus of this review: (i) the turbulence model that best characterizes the near-bed turbulent flow-structure and the turbulent flow interactions with the structure; (ii) the particle stresses, which is, the rheological models that provides a functional relation between the stresses and strains of the seabed under wave action; and (iii) the optimal boundary and initial conditions for each numerical approach.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.