*2.1. Numerical Model: REEF3D*

The numerical modelling developed in this research was conducted using the numerical model known as REEF3D [32,33], which is a computational fluid dynamic (CFD) tool used to solve Reynolds-averaged Navier–Stokes equations (RANS). This tool is able to simulate hydraulic, coastal, and estuarine phenomena for both compressible and incompressible fluids.

The REEF3D model provides a three-dimensional solution for governing equations composed of Equation (2), which corresponds to the continuity of an incompressible flow, and Equation (3), which is the momentum conservation.

$$\frac{\partial \mu\_i}{\partial \mathbf{x}\_i} = 0 \tag{2}$$

$$\frac{\partial \mathbf{u}\_i}{\partial t} + \mathbf{u}\_j \frac{\partial \mathbf{u}\_i}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ (\boldsymbol{\nu} + \boldsymbol{\nu}\_T) \frac{\partial \boldsymbol{u}\_i}{\partial \mathbf{x}\_j} \right] + \mathbf{g}\_{i\boldsymbol{\nu}} \tag{3}$$

where *i*, *j* = 1,2,3, *ui* is the mean velocity vectoral component, *xi* is the spatial vectoral component, *t* is time, ρ is water density, *p* is pressure, υ is kinematic viscosity, υ*T* is kinematic eddy viscosity, and *g* is gravity.

The turbulence closure model used was the well-known *k* − ω model [34], in which υ*T* is defined according to Equation (4), and both the turbulent kinetic energy (*k*) and the kinetic energy for the specific turbulent dissipation per unit of turbulence (ω) are determined with a transport equation, according to Equations (5) and (6), respectively. The constant values in *k* − ω were taken according to β*k* = 9/100, α = 5/9, β = 3/40, σω = 1/2, and σ*k* = 1; and *Pk* is the turbulent production rate defined by Equation (7).

$$\nu\_T = \frac{k}{\alpha} \tag{4}$$

$$\frac{\partial \mathbf{k}}{\partial t} + u\_j \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j} = \frac{\partial}{\partial \mathbf{x}\_j} \left[ \left( \nu + \frac{\nu\_T}{\sigma\_k} \right) \frac{\partial \mathbf{k}}{\partial \mathbf{x}\_j} \right] + P\_k - \beta\_k k \omega \tag{5}$$

$$\frac{\partial \boldsymbol{\omega}}{\partial t} + \boldsymbol{u}\_{j} \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \bigg[ (\boldsymbol{\nu} + \frac{\boldsymbol{\nu}\_{T}}{\sigma\_{\boldsymbol{\omega}}}) \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{x}\_{j}} \bigg] + \frac{\boldsymbol{\omega}}{k} \boldsymbol{a} \boldsymbol{P}\_{k} - \beta \boldsymbol{\omega}^{2} \tag{6}$$

$$P\_k = \nu\_T \frac{\partial u\_i}{\partial \mathbf{x}\_j} \left[ \frac{\partial u\_i}{\partial \mathbf{x}\_j} + \frac{\partial u\_j}{\partial \mathbf{x}\_i} \right]. \tag{7}$$

To solve the flow around complex structures, the ghost-cell method [35,36] was applied to impose the boundary conditions when a pile is included in the numerical domain. This approach uses fictional cells that are incorporated in the domain (specifically on the obstacles) and corresponds to a particular case of the immersed boundary method [37].

The numerical treatment of the governing equations was based on the second order of the Runge–Kutta method for temporal discretization. Meanwhile, the convective terms were solved by the applied weight essentially non-oscillatory (WENO) scheme [38]. The velocities and pressures were determined under a staggered grid, using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) [39].

The model configuration to represent waves in the numerical domain was adopted similar to the method presented by Ahmad et al. [31], using a Dirichlet boundary condition and the second-order Stokes waves theory [40]. We defined the surface elevation (η) and the incident velocity (*<sup>u</sup>*, *w*) according to Equations (8) and (9), respectively.

$$\eta(\mathbf{x}, z, t) = \frac{H}{2}\cos\theta + \frac{H^2K}{16}\frac{\cos h(Kh)}{\sin h^3(Kh)}(2 + \cos h(2Kh))\cos(2\theta) \tag{8}$$

$$u(\mathbf{x}, z, t) = \frac{\partial \Phi}{\partial \mathbf{x}} \; , \; w(\mathbf{x}, z, t) = -\frac{\partial \Phi}{\partial \mathbf{z}} \tag{9}$$

$$\Phi(\mathbf{x}, z, t) = \frac{H\mathbf{g}}{2\alpha} \frac{\cos h[\mathbf{K}(h+z)]}{\cos h(\mathbf{K}h)} \sin \theta + \frac{3}{32} H^2 \alpha \frac{\cos h[2\mathbf{K}(h+z)]}{\sin h^2(\mathbf{K}h)} \sin(2\theta),\tag{10}$$

where *H* is wave height, *K* is the wave number, *h* is the water depth, ϑ is the wave phase, is the angular frequency, and Φ is the velocity potential.

In order to avoid the wave reflection, active wave absorption (AWA), was used in the outlet according to the method described by Schäffer and Klopman [41]. In this methodology, the waves that reach the outlet cancel out the reflected waves, prescribing the velocity as in [31]:

$$
\mu\_0 = -\sqrt{\frac{g}{h}} \eta\_r \tag{11}
$$

$$
\eta\_r = \eta\_{\rm int} - h\_\prime \tag{12}
$$

where η*r* is the reflected wave amplitude and η*m* is the actual elevation of the free surface.

Complementary information on the numerical model´s hydrodynamic configuration is shown in Table 1.


**Table 1.** Complementary information for the REEF3D hydrodynamic simulation.

In order to ascertain the sediment transport and changes in the bed, both the bed load and the suspended load were considered, for the morphodynamics evolution, the equation for the conservation of sediment initially proposed by Exner [42] and later generalized by Paola and Voller [43], was used.

The bed load sediment transport was computed by using the Meyer–Peter formula [44] according to Equation (13), where *q*∗*b* is the dimensionless bed load, τ<sup>∗</sup> is the dimensionless shear stress on the bed computed by Equation (14), τ∗*c* is the dimensionless critical shear stress defined according to the Equation (15), τ is the bed shear stress, τ*c* is the critical bed shear stress, ρ*s* the sediment density, and *d* is the sediment diameter. The real magnitude of the bed load (*qb*) can be obtained from Equation (16).

$$q\_b^\* = 8(\pi^\* - \tau\_c^\*)^{\frac{\nu}{2}} \tag{13}$$

$$
\pi^\* = \frac{\pi}{(\rho\_s - \rho)gd} \tag{14}
$$

$$
\tau\_c^\* = \frac{\tau\_c}{(\rho\_s - \rho)gd} \tag{15}
$$

$$q\_b^\* = \frac{q\_b}{\sqrt{\frac{(\rho\_s - \rho)gd^3}{\rho}}}.\tag{16}$$

The suspended load was computed using the advection diffusion equation shown in Equation (17), where *c* is the suspended load concentration, and *ws* is the fall velocity of the sediment. To solve Equation (17) two boundary conditions were applied. The first is the zero vertical flux on the surface and the second is the bed load concentration (*cb*) according to Van Rijn [45], defined by Equation (18), where Υ is the relative bed shear stress defined by Equation (19), and *D*∗ is the particle parameter computed according to Equation (20).

$$\frac{\partial \mathcal{C}}{\partial t} + u\_j \frac{\partial \mathcal{C}}{\partial \mathbf{x}\_j} + w\_\delta \frac{\partial \mathcal{C}}{\partial z} = \frac{\partial}{\partial \mathbf{x}\_j} \Big(\nu \tau \frac{\partial \mathcal{C}}{\partial \mathbf{x}\_j}\Big) \tag{17}$$

$$c\_b = 0.015 \frac{d}{a} \Big| \frac{\chi^{1.5}}{D\_\*^{0.3}} \Big| \tag{18}$$

$$\Upsilon = \frac{(\tau - \tau\_c)}{\tau\_c} \tag{19}$$

$$D\_\* = d \left[ \frac{(s-1)g}{\upsilon^2} \right]^{\frac{1}{5}} \tag{20}$$

In Equation (18), *a* is the reference level for computing the suspended load that has been determined according to the methodology proposed by Rouse [46].

The critical bed shear stress was defined according to Shields [47] and parameterized according to Yalin [48], including the slope correction proposed by Dey [49] and, in addition, applying the incipient transport relaxation factor extensively described by Quezada et al. [27].

The morpho-dynamic evolution was determined by the sediment volume conservation equation described in Equation (21), where *zb* is the bed elevation, and *qb*,*<sup>x</sup>* and *qb*,*<sup>y</sup>* are the sediment transport for the bed load in the *x* and *y* directions, respectively. *E* is the sediment entrainment rate from the bed load to the suspended-load, and *D* is the sediment deposition rate from the suspended load onto the bed.

$$
\frac{\partial z\_b}{\partial t} + \frac{1}{(1 - n)} \left[ \frac{\partial q\_{b,x}}{\partial x} + \frac{\partial q\_{b,y}}{\partial y} \right] + E - D = 0. \tag{21}
$$

The difference between *E* and *D* is defined by REEF3D according to Wu et al. [50]. Meanwhile, a more extensive description of the estimation of *qb*,*<sup>x</sup>* and *qb*,*<sup>y</sup>* can be found in Quezada et al. [27].

In order to develop a numerical study of the present investigation, the model was developed as indicated in Table 2, where *L* is the flume length, *W* is the width, *ht* is the height, and Δ*xi* is the cell dimension used in the model; N◦ Cells is the total number of elements comprising the numerical domain, and *ttest* is the simulation time of the numerical model. The W (W1 to W3) and WC series (WC1 to WC3) correspond to a hydrodynamic calibration stage for waves acting alone and waves plus current, respectively, without a pile placed on the flume. Detailed information can be found in Table 3.


**Table 2.** Summary of the REEF3D model configuration applied to each simulated cases.

Δ*xi* is the dimension for the *x*, *y*, and *z* axis, due to the model using regular element definitions.

**Table 3.** Hydrodynamic calibration of waves and currents coexisting without a pile.


Series C (C01 to C04) corresponds to a hydrodynamic stage for waves plus current with a pile. Detailed information for each test can be found in Table 4. Finally, series E (E01 to E08) corresponds to the simulated cases for hydrodynamic analysis due to the combined action of waves and currents; their detailed information can be found in Table 5.


**Table 4.** Hydrodynamic calibration of waves and currents coexisting with pile.

**Table 5.** Simulated cases for the analysis of hydrodynamics due to the combined action of waves and currents.


The following sections provide more exstensive information on the simulated scenarios and the data processing used.
