*Article* **Hydraulic Characteristics of Lateral Deflectors with Different Geometries in Gentle-Slope Free-Surface Tunnels**

**Jinrong Da 1, Junxing Wang 1,\*, Zongshi Dong <sup>1</sup> and Shuaiqun Du <sup>2</sup>**


**\*** Correspondence: jxwang@whu.edu.cn; Tel.: +86-137-0718-2138

**Abstract:** The gentle-slope tunnel has been adopted in many high dams, and aerators are usually required for high operating heads. For such tunnels, the lateral deflector is superior to the traditional bottom aerator, which loses its efficiency due to cavity blockage and fails to aerate the sidewalls. However, unfavorable flow patterns such as water-wings and shock waves are induced by the lateral deflectors. To address this problem, two novel lateral deflectors are proposed, and their hydraulic characteristics are comparatively investigated together with the triangular deflector by means of model test and numerical simulation. The triangular deflector was revealed to form a wide cavity that allows for the free rise up of the water-wings inside the cavity, leading to the development of a buddle-type shock wave, whereas the two-arc deflector yields a jet with a fluctuating surface, which induces water-wings and further develops into diamond-type shock waves. In contrast, the cavity formed behind the two-arc deflector with a straight downstream guiding line is stabler and shorter, thereby restricting the development of the rising flow and preventing the formation of water-wings and shock waves. Moreover, the two-arc deflector with a straight guiding line exhibits higher energy dissipation capacities and thus is recommended in practical engineering design.

**Keywords:** lateral deflectors; gentle-slope tunnel; water-wing; shock wave; energy dissipation

## **1. Introduction**

Due to complex topographic and geological conditions [1] as well as material transport difficulties, high dams usually adopt arch concrete or local material types [2,3], the dam-body flood-discharge capacity of which is significantly limited [4]. Consequently, complementary flood release structures are of great importance for a sound design of the entire project [4]. In the past century, the tunnel spillway has stood out from other alternatives [4,5], such as chute spillways and shaft spillways owing to its high floodrelease capacity, good terrain compatibility, and decent construction cost [6]. To avoid potential cavitation and denudation damage [1,4,5,7], the flow pattern and velocity in the spillway tunnel have to be strictly controlled. The ideal flow pattern should at least meet the following two requirements: (i) the flow is free-surface flow without shock waves and with adequate space between the aerated surface and the tunnel soffit [8]; (ii) the flow velocity is within 25 m/s [8] and the structure surface is carefully smoothed [8,9]. Under such conditions, technically no aerator is required. In the past decades, the sagging dragon tail tunnel has been proven to effectively fulfill the above requirement and thus has gained popularity in many high dams (e.g., Jinping, Baihetan, Xiaowan [5], Wudongde [10], Shuangjiangkou [11]) in China. The underlying reason for this kind of tunnel's satisfactory hydraulic behavior lies in its longitudinal layout, which is reflected in its name 'sagging dragon tail'. It comprises an overflow ogee weir (dragon head) close to the reservoir, a long gentle-slope free-surface tunnel in the middle (dragon body) and a short open and steep chute connected to a ski-jump bucket (the sagging dragon tail) at the end, as shown

**Citation:** Da, J.; Wang, J.; Dong, Z.; Du, S. Hydraulic Characteristics of Lateral Deflectors with Different Geometries in Gentle-Slope Free-Surface Tunnels. *Water* **2022**, *14*, 2689. https://doi.org/10.3390/ w14172689

Academic Editor: Mouldi Ben Meftah

Received: 16 July 2022 Accepted: 26 August 2022 Published: 30 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

in Figure 1. With such a design, the flow acceleration is negligible inside the gentle-slope tunnel [4,5] and the high-velocity flow only occurs at the short steep chute, where the flow is usually highly aerated and the potential local structure destructions can be easily found and repaired once it occurs.

**Figure 1.** (**a**) Schematic diagram of the experimental spillway tunnel and details of (**b**) the simulated region and (**c**) the three lateral deflectors (unit is cm).

However, as the desired maximum operation head further increased, the flow velocity in the gentle-slope tunnel could exceed the threshold value of some 25 m/s [8], which consequently increases the cavitation risk of the tunnel [5,7,12] and the energy dissipation burden of the downstream chute. Under such conditions, aerators have to be installed to reduce cavitation risk [9,13,14], and restrictions resulting from flow pattern control have to be adequately considered. The main challenges of designing bottom aerators for gentle-slope tunnels are the low aeration efficiency caused by cavity blockage [6,12,15,16] and the cavitation protection of tunnel sidewalls [17,18]. Moreover, Hager discovered that shock waves could also be introduced by local changes in the bottom profile [14], such as offset aerators and bottom deflectors [16,19]. Therefore, to improve the aeration efficiency, lateral deflectors are usually installed together with the bottom aerators [15,17,18,20–23], of which a typical application is the fall-expansion aerators that constructed at the inlet of the free-surface tunnel behind the gate chamber [1]. It has been discovered that the lateral deflectors could not only serve as effective aerators [17,18,24] independent of the tunnel slope, but also help dissipate a certain amount of kinetic energy [24–26] and reduce the energy dissipation pressure of downstream structures. From this point, the hydraulic characteristics of the standalone lateral deflector (SLD, i.e., without bottom aerators) deserve further investigation and it is expected that the SLD could achieve decent hydraulic performance with regard to flow pattern and aeration.

To date, the commonest lateral deflector design is the triangular deflector whose up- and downstream endpoints are connected to the tunnel sidewall abruptly, and without vertical variation. The most noteworthy two of the few exceptions are probably the right-angled tetrahedrons and vertical plate deflectors investigated by Hager [26]. With such designs, water-wings and shock waves are usually induced by the lateral deflectors and thus introducing additional denudation risks. For example, Nie et al. (2006) discovered the jet formed behind the deflector impacts the sidewall and then is deflected upstream, resulting in partial blockage of the lateral cavity and the unstable swirling flow inside the lateral cavity plays a vital role in causing surface erosion [27]. The water-wings [18,20,23,24,28–30] formation was observed downstream of the impact region, which further rises up to the tunnel soffit [31] and could possibly develop into air pocket

flow [16]. Shock waves [14,16,20,24,26,29,32] were also found to be induced by lateral deflection, which could propagate far downstream, mainly in a diamond shape with water crowns colliding alternatively in the middle and on the sidewalls; thereby, causing sidewall erosion with repeating impact and pressure fluctuation. Under highly aerated situations, shockwaves could cause overtopping or trigger conduits choke [14], which dramatically impede the air transportation and thereby harm for the aeration protection [23]. Existing studies mostly focus on the optimization of the size [21,24,26,28,30] and layout [1,26,29–31] of the traditional triangular deflector, and few attempts can be found to improve the flow pattern and energy dissipation by means of modifying the deflector geometries.

In this paper, two new deflector geometries are proposed for the sake of both flow control and energy dissipation. The hydraulic characteristics of the newly proposed deflectors and the traditional triangular deflector are comparatively investigated using the hybrid approach of model test and numerical simulation. Special focus is put on the flow pattern improvement. The energy dissipation characteristics of these three deflectors are also compared based on the simulation results.

## **2. Experiment Setup**

The experiments were carried out with the physical model of a sagging dragon tail tunnel constructed at the State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University. In the experiment, water was supplied by a circulating system composed of an underground reservoir, pumps, and pipelines. The nonpressured section of the tunnel was a *B* = 18.75 cm wide, *h* = 17.75 cm deep flume inclined with a *i* = 3% bottom slope. The corresponding scale fell into the range of 1/24~1/80 with regard to the typical width being 4.5~15 m of the high unit-width-discharge freesurface tunnels (e.g., Yele, Xiaowan, Xiluodu) in China [8]. Experimental studies adopting similar geometric scale or model tunnel width can be found in [33–36]. The scale effect is considered acceptable as the main focus of this paper are the flow pattern and other macroscopical flow characteristics such as energy dissipation behaviors [36,37], which are much less sensitive to the model scale compared to two-phase flow characteristics such as air concentration, bubble size and air vent discharge [19,23]. The lateral deflectors were installed symmetrically (i.e., with identical geometry and streamwise location) on the sidewalls with our bottom aerators. The flume and deflectors are both made of plexiglass to provide a transparent view of the flow pattern. Three types of deflectors (Figure 1c) were investigated, the width *b* of which were all 0.94 cm and the streamwise length of deflectors A and B was 3.75 cm, while deflector C was further extended 1.25 cm downstream with a straight guiding line parallel to the side walls. The detailed geometries of the three types of deflectors are sketched in Figure 1. Type A is the traditional triangular deflector, deflector B is composed of two arcs tangent to each other: one being negative and has a radius of 5.625 cm and the other being positive and features a radius of 2.35 cm. As for deflector C, it is the same as deflector B except for the aforementioned additional straight line. Thus, the deflectors feature an identical contraction ratio of 10% (2*b*/*B*) and a tangent value of 1/4, and the dimensionless tail extension of deflector C with regard to the streamwise length of the curved section l is 1/3.

In the experiment, the inflow conditions were controlled at the cross-section 33.75 cm upstream of the deflector. All the measurements in this paper were conducted under the situation of inflow depth *h*in = 11.25 cm and volume flow rate *Q*in = 57.65 L/s, featuring a Reynolds number *Re* = *Qin hinR<sup>υ</sup>* <sup>=</sup> 1.2 <sup>×</sup> <sup>10</sup><sup>5</sup> (*<sup>R</sup>* represents the hydraulic diameter calculated as *bhin <sup>b</sup>*+2*hin* and *<sup>υ</sup>* = 1 <sup>×</sup> <sup>10</sup>−<sup>6</sup> <sup>m</sup>2/s is the kinematic viscosity of water). This *Re* value being larger than 1 × 105 indicates the scale effect arising from viscous stress can be neglected according to [38]. The corresponding depth-width ratio of 0.6 in the experiment is a representative value for the real-world high-discharge flood tunnels [15,22,23,33–36] in China. In this study, the flow depth was measured using a fluviograph (accuracy ± 0.18 mm). The water discharge was monitored using an electromagnetic flowmeter (IFM4080K, Jiangsu Runyi Instrument Co., Ltd., Huaian, China), featuring accuracy of 0.1 L/s. The flow velocity was measured with a propeller-type flow meter (LS300-A, Fuzhou Lesida Information Technology Co., Ltd., Fuzhou, China) featuring accuracy of 0.01 m/s.

#### **3. Numerical Models and Simulation Setup**

Numerical simulations were performed in this study to obtain full-field velocities and turbulence properties for the analysis of the underlying mechanisms of flow pattern improvement and energy dissipation characteristics. The simulations were conducted using the commercial software FLOW-3D, which is claimed to have advantages over other opponents for free-surface flows and has been widely used for spillway and tunnel flows [39–41].

## *3.1. Governing Equations*

FLOW-3D utilizes the one-fluid framework for free-surface flow modeling. This was achieved by using the Tru-VOF [42] technique to dynamically track the interface and, in the meantime, impose proper boundary conditions at the free surface. In this way, the computational cost is significantly reduced. The continuity and momentum equations are:

$$\frac{\partial \mu\_l}{\partial x\_l} = 0 \tag{1}$$

$$\frac{\partial u\_i}{\partial t} + u\_j \frac{\partial u\_i}{\partial \mathbf{x}\_j} = -\frac{1}{\rho} \frac{\partial p}{\partial \mathbf{x}\_i} + f\_i + v \frac{\partial^2 u\_i}{\partial \mathbf{x}\_j \partial \mathbf{x}\_j} - \frac{\partial}{\partial \mathbf{x}\_j} (\frac{2}{3} k \delta\_{ij} - 2 v\_l \mathbf{S}\_{lij}) \tag{2}$$

Here, *u* is velocity, *p* is the pressure, *ρ* and *υ* are the fluid density and kinematic viscosity, respectively. *f* stands for the body force, *k* and *μ<sup>t</sup>* stand for the turbulent kinetic energy, and turbulent viscosity. *δij* is the Kronecker delta, and *Sij* is the mean rate of strain tensor calculated as <sup>1</sup> <sup>2</sup> ( *<sup>∂</sup>ui ∂xj* + *∂uj ∂xi* ).

The VOF equation for interface tracking reads:

$$
\frac{\partial \alpha}{\partial t} + \mu\_i \frac{\partial \alpha}{\partial x\_i} = 0 \tag{3}
$$

where *α* is the volume fraction of the simulated fluid.

The RNG *k*-*ε* turbulence model [43] was adopted to account for the turbulence contribution to the time-averaged momentum transport, in which a transport equation was solved for the turbulent kinetic energy *k* and the turbulent dissipation rate *ε*, respectively:

$$\frac{\partial k}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_l}(ku\_l) = \frac{\partial}{\partial \mathbf{x}\_{\dot{\jmath}}} \left( (\upsilon + \sigma\_k \upsilon\_l) \frac{\partial k}{\partial \mathbf{x}\_{\dot{\jmath}}} \right) + P\_k - \varepsilon \tag{4}$$

$$\frac{\partial \varepsilon}{\partial t} + \frac{\partial}{\partial \mathbf{x}\_l} (\varepsilon u\_l) = \frac{\partial}{\partial \mathbf{x}\_j} \left( (\upsilon + \sigma\_\varepsilon \upsilon\_l) \frac{\partial \varepsilon}{\partial \mathbf{x}\_j} \right) + \mathbf{C}\_{1\varepsilon} \mathbf{P}\_k \frac{\varepsilon}{k} - \mathbf{C}\_{2\varepsilon}^\* \rho \frac{\varepsilon^2}{k} \tag{5}$$

where *υ<sup>t</sup>* = *C<sup>μ</sup> k*2 *<sup>ε</sup>* , *Pk* = 2*υtSijSij*, *C*<sup>∗</sup> <sup>2</sup>*<sup>ε</sup>* = *C*2*<sup>ε</sup>* + *Cμη*<sup>3</sup> <sup>1</sup><sup>−</sup> *<sup>η</sup> η*0 <sup>1</sup>+*βη*<sup>3</sup> , *<sup>η</sup>* <sup>=</sup> <sup>2</sup>*Sij*·*Sij*1/2 *<sup>k</sup> <sup>ε</sup>* , and model parameters are: *σ<sup>k</sup>* = *σε* = 1.39, *C<sup>μ</sup>* = 0.085, *C*1*<sup>ε</sup>* = 1.42, *C*2*<sup>ε</sup>* = 1.68, *η*<sup>0</sup> = 4.38, and *β* = 0.012.

## *3.2. Simulation Setup*

The computational domain adopted the same coordinate system as the experiment and extends from *x* = −33.75 to *x* = 135 cm in the streamwise direction. A grid convergence study involving three mesh schemes was carried out in advance to select a proper mesh for the simulation. The details of the grid convergence study were shown in Appendix A. One structured mesh block was used to discretize the domain into 10.25 million cuboid cells, the majority of which feature an average size of 0.34 cm × 0.125 cm × 0.2 cm (*x* × *y* × *z*). To

ensure the spatial resolution of the deflector geometry using the FAVOR® technique [44,45], the computational domain was rotated by the flume slope of 3% to be aligned with the *x*-coordinate and the gravity vector *g* was adjusted accordingly. Moreover, the domain was locally refined to 0.188 cm × 0.094 cm × 0.188 cm in the region of −7.5 < *x* < 30 cm, −10 < *y* < −6.25 cm and 6.25 < *y* < −10 cm, and 0 < *z* < 12.5 cm.

The outlet boundary was configured with the outflow boundary condition (BC), implying a zero-gradient condition since the outflow is supercritical. The bottom and sidewalls were set as non-slip walls with equivalent roughness *k*s = 0.015 mm according to the plexiglass surface properties. The top boundary is specified with fixed relative pressure *p* = 0 (i.e., atmosphere pressure). As for the inlet boundary, a pre-simulation of the entire physical model (i.e., from the reservoir to the downstream cushion pool) was conducted first, and then all the flow parameters at *x* = −33.75 cm were mapped onto the inlet of the short domain using the grid overlay BC. The illustrative diagram of the simulation setup is shown in Figure 2.

**Figure 2.** Schematic diagram of the simulation setup.

## **4. Results and Discussion**

## *4.1. Model Validation*

The simulated and measured water surface profiles at two longitudinal slices are comparatively shown in Figure 3a. It can be seen that the simulation results are in good agreement with the experimental data. Only small deviations can be found near the deflectors (−5 < *x <* 10 cm at *y* = −8.375 cm). Moreover, vertical and horizontal profiles of the measured and calculated streamwise velocity at *x* = 112.5 cm are presented in Figure 3b and c respectively. The simulated vertical velocity profile overlaps well with the measured data in the middle of the flume (i.e., *y* = 0 cm), whereas the calculated velocities near the sidewall (*y* = −8.38 cm) are approximately 10% smaller than the experimental data, which can be attributed to the wall function effects. As for the horizontal velocity distribution, the simulated velocities are also generally consistent with the measured data, with the bottom values (*z* = 1 cm) noticeably lower than those in the middle (*z* = 6 cm) and near the surface (*z* = 11 cm). While the bottom and middle velocities exhibit typical U-shape horizontal profiles, the surface velocities show a remarkable attenuation from the middle to both sides due to influence of shock waves, which are also reflected in the fluctuant measured data. Nevertheless, the overall accuracy of the simulation is considered acceptable to support the analysis of the flow pattern and energy dissipation characteristics.

**Figure 3.** Comparison of the experimental and numerical results: (**a**) water surface profile; (**b**) vertical; and (**c**) horizontal distribution of the streamwise velocity *ux* at *x* = 112.5 cm. The unit of the coordinates is cm.

#### *4.2. Flow Pattern*

Figure 4 shows the side view of the flow pattern downstream of three different deflectors. Lateral jets can be found to form behind the deflectors, and the jets travel further at higher elevations, leading to a non-uniform lateral cavity that is consequently longer and wider at higher elevations. The jet impacting the sidewall is deflected immediately and keeps rising, clinging to the sidewall, associated with noticeable turbulence and aeration. In contrast, the water below the impacting region remains non-aerated. The three deflectors also form different boundary lines between the white and black water regions. With type A, the boundary line is smooth, above which the deflected flow develops into water-wings and continually rises over the sidewalls. With type B, the boundary line is fluctuant, and the deflected flow develops into water-wings that intermittently jump over the sidewall, and a few bubbles could be found in the lower clear water. As for type C, neither distinguished boundary line nor water-wings can be found since the aeration in the lower region is also noticeable. It appears that the deflected flow is submerged in the lateral jet and induces stronger turbulence and intensive aeration since a larger amount of air bubbles can be seen at the impacting region throughout the flow depth. This is probably caused by the rapid close of the lateral cavity and the stronger interaction between the jet and the deflected flow.

**Figure 4.** Side view of the flow pattern downstream of the deflectors.

The jet trajectories behind three different deflectors are marked with white lines and shown in Figure 5. The detaching jet behind deflector A moves in the direction tangent to the hypotenuse, leading to the widest and largest lateral cavity, which further leaves enough space for the lower deflected flow to develop to the upper region without contacting the jet trajectory. However, the jet behind deflector B exhibits a fluctuating phenomenon and forms an unstable cavity smaller and narrower than that behind deflector A. The fluctuating phenomenon can be mainly ascribed to the continuously varying edge slope of the deflector since it could yield non-uniform velocity distribution that exacerbates jet surface fluctuations. A benefit from the additional straight line in comparison to deflector B, the flow at the tail of deflector C has developed to a more uniform condition before detaching away from the deflector. It exhibits a trajectory almost parallel to the sidewall right behind the deflector, forming the narrowest and smallest lateral cavity. The rapidly closed cavity blocks the rising passage of the underneath deflected flow, preventing the formation of water-wings.

**Figure 5.** Top view of the flow pattern near the lateral cavity.

The development of water-wings and the subsequent shock waves downstream of the three deflectors are shown in Figure 6. It can be discovered that deflectors A and B form shock waves in bundle shape and diamond shape, respectively, whereas no shock wave is observed downstream of deflector C. With deflector A, the water-wings from both sides quickly rise up and detach from the main flow, jetting to the middle and colliding with each other in the air, and then fall back to the main flow with a certain vertical velocity and the lateral momentum dissipated. The returned water further produces disturbance on the surface and thereby inducing a bundle-shaped shock wave. With deflector B, the lateral momentum of the water-wings partially counteracts the jet in the cavity, and thus the water-wings fail to contact each other before falling back to the main flow. The returned water still contains certain lateral momentum and moves together with the main flow, producing two repeated reverse oblique developing lines on the flow surface and finally developing into a diamond-shaped shock wave. As for deflector C, the deflected flow is restricted in the relatively small cavity and thus only forms small water crowns, which are restricted in the vicinity of the sidewall and disappear further downstream.

**Figure 6.** Top view of the flow pattern downstream of the deflectors.

Based on the above observations, it can be concluded that a continuous variation of the lateral deflector surface at the tail with an additional flow guiding extension is the key to the elimination of the water-wings and shock waves.

## *4.3. Velocity Distribution*

The near-bottom (*z* = 1 cm) velocity distribution around three deflectors is comparatively shown in Figure 7. The jets behind deflectors A and B are firstly contracted and reach the maximum cavity width at approximately 1/3 of the cavity length. After that, the jets restart to spread to the sidewall. In comparison, the jet behind deflector C continuously moves close to the sidewall. Consequently, deflector A exhibits the largest cavity size, followed by deflectors B and C, which is in agreement with the experimental observation. It is worth emphasizing that the velocity vector distribution around deflector B is apparently more non-uniform compared to those around deflectors A and C, which is the reason for the fluctuating phenomenon of the jet surface. Moreover, the maximum velocity behind deflector C is about 0.2 m/s lower than those behind the other two deflectors, and the near-wall low-velocity region (i.e., the purple dashed rectangle in Figure 7) behind deflector C is obviously larger than those behind deflectors A and B. These phenomena all imply more effective energy dissipation of deflector C compared to A and B.

**Figure 7.** Horizontal velocity distribution around the deflectors at *z* = 1.0 cm. The vectors are generated using the same density and scale for the three types of deflectors.

The near-wall (*y* = −9.3 cm) vertical distribution of velocity behind the deflectors is comparatively shown in Figure 8 to illustrate the kinematic characteristics of the deflected flow. Distinguished by the velocity magnitude and vector direction, the longitudinal flow behind deflectors A and B can be divided into three distinct regions: the main flow region that flows downstream, the water-wing region that rises up while flowing downstream and the impacting region that exhibits the highest velocities and separating the above two

regions. The rising water-wing regions are in agreement with the experimentally observed boundary lines. As for deflector C, the flow direction is almost not affected by the jet since only a small region can be found that has upward velocities and therefore is the underlying reason for the fast recovery of the rising surface and the elimination of shock waves. Moreover, the velocities behind deflector C are lower than those behind deflectors A and B, especially dropping about 50% at the impacting region. This again confirms the more effective energy dissipation of deflector C.

**Figure 8.** Velocity distribution near the sidewall (*y* = −9.3 cm). The vectors are generated using the same density and scale for the three types of deflectors.

Figure 9 shows the streamwise development of the cross-sectional velocity (CSV) with three different and without deflectors, which is computed from 17 flux surfaces with a spacing of 7.5 cm. The no-deflector scheme shows a linearly increased CSV, whereas for the other three schemes, the CSVs increase at the deflector region and decrease at the cavity region because of lateral flow contractions and expansions. With deflectors A and B, the CSVs first increase in the front of the cavity region and then decrease, whereas the CSVs of the flow behind deflector C exhibits a consistent decreasing trend due to continuous flow expansion. Affected by the shock waves, the CSVs of the downstream flow with deflectors A and B exhibit some fluctuations, and the fluctuating amplitude of deflector B is higher than that of deflector A. In contrast, the CSVs of deflector C feature a linear increase behavior similar to that of the no-deflector scheme in this region, although with lower values. Moreover, it can be found that the efficient energy dissipation mainly occurs at the impacting region (i.e., approximately 15 < *x* < 40 cm).

**Figure 9.** Streamwise development of the cross-sectional velocity along the flume with three different deflectors and without deflectors. The brown color indicates regions where the cavity range of C (pink) overlaps those of A and B (green).

## *4.4. Energy Dissipation Characteristics*

Figure 10 shows the streamwise development of the flux-averaged hydraulic head (FAHH) along the flume with three different deflectors and without deflectors. In FLOW-3D, the FAHH is recommended to evaluate the hydraulic head where significant changes occur. It is calculated using the following equation:

$$h = \frac{\int \left(\frac{\mu^2}{2\mathcal{g}} + \frac{p}{\rho\mathcal{g}} + z\right) ds}{\int \mathcal{q} ds} \tag{6}$$

here, *h* is the FAHH, *ϕ* is the flux across a cell surface, and *ds* is the open area of the cell face.

**Figure 10.** Streamwise development of the flux-averaged hydraulic head (FAHH) along the flume with three different and without deflectors. The brown color indicates regions where the cavity range of C (pink) overlaps those of A and B (green).

The FAHH decreases linearly in the flume without lateral deflectors. The FAHH of the flow underwent some certain fluctuations around the no-deflector values when passing deflectors A and B. These fluctuations continue until they are twice the cavity length downstream (i.e., to approximately *x* = 60 cm). Further downstream, the FAHH values of the flow behind deflectors A and B turn out to be lower than the no-deflector values. In contrast, a rapid decrease in FAHH can be found at the deflector and cavity region (0 < *x* < 20) for the flume equipped with deflector C, and further downstream (*x* ≥ 45), the values of FAHH decrease steadily with a similar rate as those of the no-deflector scheme. A quantitative comparison of the energy dissipation effect of the three deflectors can be achieved using the local head loss coefficient *<sup>ζ</sup>* = (FAHH1 − FAHH2) · 2g/*v*<sup>2</sup> between *x* = −17 cm and *x* = 45 cm with *v* being a reference CSV calculated at *x* = −17 cm. The local head loss coefficients for deflectors A, B and C are 2.72%, 2.90% and 4.2%, respectively. This behavior is in agreement with the streamwise development of the CSV and confirms the kinetic energy is mainly dissipated at the deflector and jet region.

To further analyze the energy dissipation mechanism of the lateral deflectors, the distribution of turbulent kinetic energy *k* and turbulence dissipation rate *ε* in the deflector and jet region are comparatively shown in Figure 11. The most noticeable difference in Figure 11 is the larger *k* and *ε* values at the impacting region behind deflector C compared to those behind deflectors A and B. This implies more mean kinetic energies are turned into turbulent energies and are then dissipated for the flow behind deflector C. This is consistent with the intensive aeration observed in the flow photos shown in Figure 6, as higher turbulent levels induce more intensive aeration [44–46].

**Figure 11.** Horizontal distribution of (**a**) turbulent kinetic energy *k*; (**b**) turbulence dissipation rate *ε* in deflector and cavity region at *z* = 1 cm.

## **5. Conclusions**

In this paper, the hydraulic characteristics of lateral deflectors with three different geometries in gentle-slope free-surface tunnels were investigated using the hybrid approach of model test and numerical simulation. Special focus was placed on the flow pattern and energy dissipation features. The main findings are:


caused by the continuously varying curvature of the arcs. Water-wings also develop inside the cavity and eventually produce a diamond-type shock wave downstream. In contrast, the jet behind the two-arc deflector with a straight guiding line at the tail is stabler and travels a shorter distance before impacting the side wall. The jet could thus restrict the development of the rising flow, and thereby eliminate the formation of water-wings and shock waves. Based on these observations, it is concluded that a continuous variation of the lateral deflector surface at the tail with an additional flow guiding extension is the key to the elimination of the water-wings and shock waves.

3. Compared to the triangular deflector and the two-arc deflector, the two-arc deflector with a straight line exhibits more effective energy dissipation, as reflected in the local energy loss coefficient. The underlying reason for its effective energy dissipation is the more intensive turbulence introduced by the stronger interaction between the deflected flow and the jet surface, which also leads to more intensive aeration.

Compared to real-world engineering problems, the present study has some limitations. The physical model is roughly 1/24~1/80 of real-world free-surface tunnels, and thus scale effect is expected, in particular for air–water two-phase flow characteristics [47–49]. However, the scale of the physical model is considered acceptable with regards to the main concern of the current paper, which is the comparative evaluation of the different lateral deflectors in terms of the flow pattern and the energy dissipation behaviors. Particularly, the jet trajectory behind the lateral deflector and its interaction with the deflected flow could be decently reproduced in the physical model. Moreover, considering that the numerical simulation approach remains immature for highly aerated flows but much more reliable for non-aerated flow [13,50,51], the hybrid approaches used in this paper could offer valuable insight into the underlying reason for the flow pattern improvement and higher energy dissipation observed with deflector type C. Aeration is another key concern for lateral deflector designs but is not investigated in detail in this paper due to facility reasons. Nonetheless, the findings from the study provide a novel and preferable lateral deflector design for gentle-slope free-surface tunnels, which could, if not resolve, significantly improve the unwanted flow patterns of water-wings and shock waves. Meanwhile, it could also achieve improved energy dissipation compared to traditional alternatives. The aeration characteristics of the novel deflector will be investigated in future research.

**Author Contributions:** The conceptualization of this research was directed by J.W. and the experimental data were from J.D. and S.D.; the numerical modeling was performed by J.D. under the supervision of J.W. and Z.D.; the validation is carried by J.D. and S.D.; the manuscript was mainly finished by J.D., Z.D. and J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Key Research and Development Program of China (Grant No. 2021YFC3090105), the China Postdoctoral Science Foundation (No. 2021M692754), and Guizhou Science and Technology Joint Support (2019) No. 2890.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Appendix A**

A grid convergence study involving three mesh schemes was carried out to evaluate the grid quality. The cell size (*x* × *y* × *z*) in the refined region of the three mesh schemes was 0.625 cm × 0.313 cm × 0.625 cm (coarse), 0.375 cm × 0.188 cm × 0.375 cm (medium), and 0.188 cm × 0.094 cm × 0.188 cm (fine), respectively.

The grid spacing of the 3D grids is defined as *hk* = <sup>3</sup> *hx*,*k*·*hy*,*k*·*hz*,*k*. The computational solution *f* for a grid *k* is defined as the computed maximum streamwise velocities *u*max in the middle longitudinal section at *x* = 112.5 cm. The *GCI* for each grid *k* is computed as *GCIk* = *Fs*|*Ek*|, where a value of 1.25 is chosen for the factor of safety *Fs*, following the recommendations of Roache [52].

Table A1 shows the grid convergence index (GCI) calculated using the approach of Roache and Salas [52]. The observed order of accuracy is *p* = 1.84, and the *GCI* shows a convergence trend. The fine grid takes a very satisfactory value of 0.224%, which is an order of magnitude higher than the former and therefore is considered fine enough to resolve the flow field.


**Table A1.** Parameter of grid convergence calculation.

The vertical distributions of streamwise velocity *u*x and turbulent kinetic energy *k* in the middle longitudinal section at *x* = 112.5 cm are comparatively shown in Figure A1. From Figure A1, it can be found that the main deviations of the *u*x profiles lie in the region close to the water surface and the bottom wall, whereas in the middle, all three meshes return almost identical velocity values. Moreover, the difference between the medium and fine meshes is less obvious compared to that between the coarse and the medium meshes. As for the *k* profile, the result from the coarse mesh is strikingly different from those from the medium and fine meshes, the difference between which is relatively negligible. Therefore, the fine grid was selected for the simulation.

**Figure A1.** Vertical distributions of streamwise velocity *ux* and turbulent kinetic energy *k* in the middle longitudinal section at *x* = 112.5 cm.

#### **References**


## *Article* **SWAN Modeling of Dredging Effect on the Oued Sebou Estuary**

**Nisrine Iouzzi 1,\*, Laila Mouakkir 1, Mouldi Ben Meftah 2,\* , Mohamed Chagdali <sup>1</sup> and Dalila Loudyi <sup>3</sup>**


**Abstract:** The estuary ecosystem's health and ecological integrity are essential for preserving environmental quality, habitats, and economic activity. The main objective of the present study is to comprehend the wave hydrodynamic impact on the Oued Sebou estuary, which is situated in the Kenitra region on Morocco's north Atlantic coast in North Africa. Specifically, it focused on the dredging effect (caused by sand extraction) on the wave motion and its impact on the estuary environment. Different scenarios of wave-propagation simulations were carried out, varying the significant wave height, in deep water (from 1.5 to 4 m), and considering the bathymetry before and after two dredging cases of 2- and 4-m depths. The change of wave height at the Oued Sebou estuary shoreline was simulated by using the third version of the Simulating Waves Nearshore Model (SWAN). The SWAN model formulates the wave evolution in terms of a spectral energy balance on a structured grid. The effect of dredging on the wave spreading in addition to the flow hydrodynamic structures were extensively analyzed. According to the simulated results, the dredging activities in the Oued Sebou estuary mainly affect the river mouth and the southern breakwater area, increasing the potential erosive action. The areas at the northern coastal strip and near the northern breakwater are subject to possible accumulation of sediments.

**Keywords:** Oued Sebou estuary; bathymetry; dredging; SWAN; wave spreading; flow hydrodynamic structures

## **1. Introduction**

Estuaries are one of the most important interconnections between land and sea. They are frequently important areas for leisure and economic activities. Estuaries are also ecosystems that are very susceptible to pollution and environmental disturbances and must be adequately preserved. In estuaries, rivers and oceans interact through a number of intricate phenomena: exchange of water of different densities, sediments, pollutants, nutrients, organic matter, and biota. Therefore, it is necessary to have a better understanding of the hydrodynamics of estuaries and coasts in order to effectively assess and monitor the environment quality in these regions and to predict coastal evolution [1]. It takes appropriate methodologies combining theoretical analysis and modeling studies to anticipate the hydrodynamics of estuaries and coasts, which are caused by complicated mechanisms linking mass exchange to heat transfer processes [1].

Recently, several research studies were conducted to investigate the hydrodynamics of estuaries [2–12]. Many driving conditions, such as river flow characteristics, tidal amplitude, sediment properties, surges, waves, currents, and bathymetry, affect the dominant physical processes in an estuary. The primary sources of energy input into estuaries are typically tides, surges, and waves. Wave dispersion in shallow water is significantly impacted by the change in bathymetry, leading to refraction, diffraction, reflection, and shoaling. Natural beaches and manmade coastal structures that reflect waves can have a significant impact

**Citation:** Iouzzi, N.; Mouakkir, L.; Ben Meftah, M.; Chagdali, M.; Loudyi, D. SWAN Modeling of Dredging Effect on the Oued Sebou Estuary. *Water* **2022**, *14*, 2633. https://doi.org/10.3390/w14172633

Academic Editor: Giuseppe Oliveto

Received: 14 July 2022 Accepted: 21 August 2022 Published: 26 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

on the hydrodynamic structures and, consequently, the transport of sediment in front of the reflector [13]. According to Chang and Hsu [13], the prediction of wave reflection coefficients is still a challenging task. Before marine structures are created, wave reflection at them is typically evaluated with physical models for engineering practice. Incident waves are also reflected by the sloping bathymetry. The drag coefficient estimation is complicated by the superposition of incident and reflected waves (due to the variation in bathymetry and the existence of any natural or man-made structures). In their study, Chang and Hsu [13] compiled many approaches for predicting the reflection coefficient of the waves suggested in the literature and highlighted their limitations. Most of these methods, according to the authors, are ineffective at predicting the reflection coefficients of waves on a bed with arbitrary or natural bathymetry. Based on a linear wave shoaling theory (without considering wave breaking), determining changes in wave height and phase due to bathymetry variations, Chang and Hsu [13] proposed a simple frequency-domain method for separating incident and reflected waves to account for normally incident linear waves propagating on an inclined bed with arbitrary 2D bathymetry. According to the authors, their method is applicable to both laboratory and field conditions predominantly for normal shores on which deep-water waves are propagated.

In intermediate and shallow waters, waves from deep waters begin to interact with the seabed, undergoing various transformations. These transitions are the results of many phenomena, for example, refraction, diffraction, reflection, shoaling, breaking, and swash. Wave prediction is of crucial importance to human activities as it provides useful information for many coastal engineering applications such as the coastal protection, environmental monitoring, navigation, safe port management, and good control of recreational activities in popular coastal areas [14]. Both observational and measurement data, as well as physical and numerical models, are often used to assess extreme marine events around the world [15]. Chi and Rong [15] confirmed that long-term in situ observations can accurately estimate the sea level variation but are usually spatially limited. The authors also indicated that many numerical models have been developed and validated to reproduce the spatial and temporal variation of wave spreading and uncover the underlying dynamics. The sea wave motion is strongly nonlinear and greatly influenced by many factors, i.e., the seabed, wind velocity, current circulations, and induced radiation stress. The interaction of the waves with different factors leads to complex hydrodynamic behaviors, making their numerical simulation very uncertain, which always requires validations using measured data.

The Simulating Waves Nearshore (SWAN) model is one of the most very popular wind–wave models [16–27], used by many government organizations, research institutes and consulting companies worldwide, to predict the size and forces of waves, allowing for changes in wave propagation from deep water to the surf zone [14]. Based on SWAN manuals and as outlined in Lin [14], the model's primary function is to resolve the spectral action balance equations, which represent the effects of spatial propagation, refraction, shoaling, generation, dissipation and nonlinear wave–wave interactions, without imposing any a priori constraints on the spectral evolution of wind waves. Hoque et al. [16] applied the SWAN model to forecast storm wave conditions in the shallow nearshore region off the Mackenzie Delta in the Canadian Beaufort Sea. The standard setup for the SWAN model was implemented by the authors, comparing different methods for quantifying the wave whitecapping dissipation. Hoque et al. [16] found that, after examining the bottom friction effects and triad interactions in predicting shallow-water waves, the simulated results in terms of significant wave heights and peak period are in good agreement with the observed data. Amunugama et al. [17] analyzed wind waves with SWAN on structured mesh and unstructured mesh during the arrival of a typhoon. After comparing the simulated results with measured data, they confirmed that the wind–wave characteristics obtained by both approaches (structured SWAN and unstructured SWAN) were substantially consistent with some advantages of unstructured SWAN, especially in complex geometries. The authors recommended the combination of both approaches where necessary. Due to

a lack of time-series wave data, Gorman et al. [19] have used SWAN and wind data to hindcast the generation and propagation of deep-water waves incident on the New Zealand coast over a 20-year period (1979–1998). The SWAN model was also used to analyze spatial and temporal variations in cold front events [20], predict waves generated by cyclones [21], simulate wave characteristics at shorelines [22,24], and estimate wave energy potential [23], etc.

This study aims at investigating the effect of extreme dredging depths (due to sand extraction) on the wave characteristics in the Oued Sebou estuary, located in the Kenitra region on the north Atlantic coast of Morocco by applying the SWAN model (Cycle III version 41.20). The wave hydrodynamic characteristics around the river outlet were examined. The dredging impact on wave dispersion, bottom friction velocity field, flow agitation in the river, and stress on the structures housing the estuary and wave energy budget (before breaking) were extensively analyzed.

## **2. Study Area**

This study covers the Oued Sebou estuary area. The Sebou estuary is located on the Atlantic coast in the Kenitra region of Morocco in North Africa. It is a region where numerous coastal developments are located. Two concentrically spaced longitudinal breakwaters channel the Oued Sebou over a distance of about 800 m, extending to a bathymetric line of 7 m. The sandbar that prevents shipping in the canal was reduced by the construction of the breakwaters in 1931 (Figure 1). During low tide, these structures give the tidal current the ability (velocities and forces) to drive the sand along the shoreline in either direction. Additionally, this river is controlled by water agencies that guarantee its stability. The Mehdia beach (of the city side in the southern part of the river mouth) is surrounded by a corniche composed of a wall with an average height of 1.50 m and 2 m for the northern and southern halves, respectively. In the southern and northern beach zones, respectively, the wall's height above sea water level (SWL) is nearly 11 and 13 m. It is also important to note that the Mehdia shoreline dune is occupied by various manmade structures.

Dredging operations were carried out every year with an average volume of almost 460,000 m3/year. This makes it possible to ensure adequate depth, which is necessary for the navigation of the fishing vessels and ships in direction of the Mehdia ports, located inside the channel (Figure 1).

**Figure 1.** *Cont.*

**Figure 1.** Outlet of the Oued Sebou into the Atlantic Ocean (above: Google Maps Platform) and a general overview of the study region (below).

## **3. Model Theoretical Formulation**

The selected SWAN model is a spectral numerical model designed to simulate waves evolving in coastal regions, lakes, and estuaries under defined wind, bathymetry, and current conditions. It is based on the Energy Density Balance equation linking the advection term to the source and sink terms. Therefore, the wave energy evolves in both geographic and spectral space and changes its aspect due to the presence of wind at the surface, friction with the bottom, or during the breaking of the waves. The SWAN model is a stable model based on unconditionally stable numerical schemes (implicit schemes).

SWAN, in its third version, is in stationary mode (optionally non-stationary) and is formulated in Cartesian or spherical coordinates. The unconditional numerical stability of the SWAN model makes its application more effective in shallow water.

In SWAN, the waves are described with the two-dimensional spectrum of the wave action density *N*,

$$N(\mathbf{x}, y, \sigma, \theta) = \frac{E(\mathbf{x}, y, \sigma, \theta)}{\sigma} \tag{1}$$

where *x* and *y* are the horizontal components of geographic space, *σ* is the relative frequency, *θ* is the wave direction, and *E* is the energy density.

The spectrum considered in the SWAN model is that of the wave action density *N* (*σ*, *θ*) rather than the spectrum of the energy density *E* (*σ*, *θ*). This is because, in the presence of currents, the wave action density is conserved while the energy density is not [27].

Because wave action propagates in both geographic and spectral space under the influence of genesis and dissipation terms, wave characteristics are described in terms of two-dimensional wave action density. The action density spectrum balance equation relating the propagation term to the effects of the source and sink terms, in Cartesian coordinates, is [28]

$$\frac{\partial N}{\partial t} + \frac{\partial (\mathcal{C}\_{\mathcal{X}} N)}{\partial \mathbf{x}} + \frac{\partial (\mathcal{C}\_{\mathcal{Y}} N)}{\partial y} + \frac{\partial (\mathcal{C}\_{\mathcal{C}} N)}{\partial \sigma} + \frac{\partial (\mathcal{C}\_{\theta} N)}{\partial \theta} = \frac{S}{\sigma} \,. \tag{2}$$

On the left-hand side of Equation (2), the first term represents the local temporal variation of the wave action density, the second and third terms represent the propagation of wave action in the geographical space of velocities *Cx* and *Cy*, the fourth term represents the shifting of the relative frequency due to variations in bathymetry (with propagation velocity *Cσ*) and currents (with propagation velocity *Cθ*), and the fifth term represents the refraction induced by the combined effects of depth and currents. It is worth mentioning that *Cx*, *Cy*, *Cσ*, *C<sup>θ</sup>* propagation velocities were obtained from linear wave theory. The term in the right-hand side of Equation (2) represents processes that generate, dissipate, or redistribute wave energy, and *S* can be expressed as [24]

$$S = S\_{\rm in} + S\_{\rm uvc} + S\_{\rm brk} + S\_{\rm bot} + S\_{\rm n14} + S\_{\rm n13} \,\,\,\,\,\tag{3}$$

where *Sin* is the wind energy input. The dissipation term of wave energy is represented by the contribution of three terms: the whitecapping *Swc*, bottom friction *Sbot*, and depthinduced breaking *Sbrk. Sn*l4 and *Sn*l3 represent quadruplet interaction and triad interactions, respectively.

Adopting a finite difference scheme for each of the five dimensions: time, geographic space, and spectral space made the numerical implementation in SWAN effective. The following guidelines must be followed in order to obtain the discretization adopted at the SWAN model level for numerical computation: (i) time of a constant and identical time step Δ*t* for the propagation term and the source term, (ii) geographical space of a rectangular grid with constant spatial steps Δ*x* and Δ*y*, (iii) spectral space of a constant directional step Δ*θ* and a constant relative frequency step Δ*σ*/*σ*; (iv) frequencies between a fixed minimum maximum values of 0.04 Hz and 1 Hz, respectively, unlike the WAM and WAVEWATCHIII models, which use dynamic values, and (v) as an option, the direction *θ* can also be delimited by the minimum and maximum values *θ*min and *θ*max.

## **4. Model Forcing Data**

## *4.1. Bathymetry*

The bathymetry of the study area is generated by using the SHOM charts, SHOM7702 SHOM Nautical Charts—Morocco. Available online: https://www.nauticalchartsonline. com/charts/SHOM/Morocco (accessed on 1 February 2021).

Figure 2a shows the bathymetry of the coarse grid domain, while Figure 2b illustrates a more detailed bathymetry. The seabed zone is characterized by comparatively regular isobaths parallel to the coastline. According to the bathymetry of the research region, shallow water extends up to 550 m from the coast with a depth range between 2 and 3 m.

**Figure 2.** *Cont.*

**Figure 2.** Bathymetry of Sebou estuary. (**a**) Bathymetry of the coarse grid domain. (**b**) Detailed isobaths of the study area.

Generally, for wave simulation, the research and commercial models use flexible grid mesh (structured and unstructured grids). In previous studies [17,21], it was observed that simulations with structured or unstructured grids were substantially consistent. In this study, we use a nested regular grid with a resolution of approximately 25 m, as shown in Figure 3. The bathymetry meshing was generated by using the BlueKenue software.

**Figure 3.** Mesh used of the study area: regular grid of approximately 25 m resolution.

## *4.2. Wind and Wave Fields*

Because the accuracy of the wind field has a large impact on the predicted wave fields [16], in this study, the main meteorological parameters were well analyzed. The average annual values of the wind properties indicate that the study area is characterized by: (i) a winter regime for which the dominant wind (or most frequent) mainly comes from the eastern sector (onshore wind). During this period, the strongest wind speed (>9 m/s) comes from the directions ranging between the southwest and west sectors, with an occurrence of almost 3%, (ii) a summer regime (from March to October), for which the dominant wind comes mainly from the western sector (sea wind). In this period, the strongest wind speed (>9 m/s) comes from the west and north sectors, with an occurrence of almost 3%.

The meteorological data (Figure 4) were collected from the nearest weather station to the study area. It is situated at the Kenitra airport, which is 8 km from the study location. The collected measured data correspond to the period from 1990 to 2009. The annual predominant wind direction is from the west and its speed range from 4 to 9 m/s (Figure 4c).

**Figure 4.** *Cont.*

**Figure 4.** Wave rose (**a**) of the significant wave heights, (**b**) corresponding periods of the significant wave, (**c**) and average annual wind rose.

Figure 4 shows the annual wave rose of the significant wave heights (Figure 4a) and their corresponding periods (Figure 4b). The dominant wave (the most frequent) comes from the northwest (300◦N to 320◦N). Storm swells are more westerly and come from the west–northwest sector (280◦N to 300◦N). Five significant wave heights were considered in this study (results of 19-year data analysis), *Hs*<sup>0</sup> = 1.5, 2, 2.5, 3, 3.5, and 4 m. Here, we denote by *Hs*<sup>0</sup> the significant wave height in deep water.

## *4.3. Tidal Current and Water Level*

The average sea level in the region of Kenitra is almost equal to 2.17 m. The tide is of semi-diurnal type with a period of 12 h 25 min. Previous studies confirmed that the tidal current and water level significantly affect wave behaviors [24]. As a result, in this study, the tidal forcing was also considered as SWAN mode input data l. Table 1 provides a summary of the astronomical tide level values for the study area.


**Table 1.** Tidal levels in the study area.

## *4.4. Model Setup*

The site of the Oued Sebou estuary is exposed to the west–northwest direction, which is also the direction of the main dominant wave. Based on the layout of the coastline and the regularity of the hydrodynamic solicitations, the main sediment transport occurs along the channel profile. Monitoring of the estuary area shows that there is no significant littoral transit.

Regarding the primary swell direction (N300), the most exposed areas to wave ac-tion are located downstream of the river. These locations have been given the labels A, B, C, and D, as shown in Figure 5.

**Figure 5.** The different zones subject to analysis: stress on structures (A, C), B-agitation in the river outlet (B), littoral erosion (D).

After considering the wave and tidal conditions, we chose to proceed with a selection of eight (8) modeling cases (see Table 2), combining the considered five significant wave heights, indicated above, and different tide levels ranging between 0 and 3.4 m. The selection of wave characteristics covers the most representative samples. The current bathymetry (without dredging), the bathymetry with a dredged area of uniform extraction depth of 2 m (with respect to the current bathymetry), and the dredged area of uniform extraction depth of 4 m were all taken into consideration. Moreover, the SWAN model was used with the following additional hydrodynamic parameters: energy spectral density distribution (JONSWAP) of a spectrum width parameter at 3.3 and a low angular spread, spectrum discretization of 36 angular sectors and 32 frequency intervals (between 4 and 20 s), the refraction and diffraction phenomena were considered, and the wind turbulence effects were not considered.


**Table 2.** Specification of the hydrodynamic cases considered.

## **5. Results and Discussion**

As an illustration, Figure 6 depicts the simulated outcomes of the local significant wave height (*Hs*) distribution, before (with current bathymetry), as shown in Figure 6a, and

after dredging of 2 m, as shown in Figure 6b. The data illustrated in Figure 6 refer to case 3 (Table 2) of significant wave height in deep water *Hs*<sup>0</sup> = 1.5 m and a period *Tp* = 8 s.

**Figure 6.** Map of local significant wave height distribution for (*Hs*<sup>0</sup> = 1.5 m, *Tp* = 8 s), (**a**) before dredging (with current bathymetry), (**b**) after dredging of 2 m depth. The legend indicates the values of *Hs* in meters.

The results show that the presence of the dredging affects the wave dispersion. There was an increase in *Hs* in and at the edges of the dredging area, which can be explained by the sudden increase in the water depth. In general, the simulations show a decrease in the transfer of wave energy out of the extraction zone, causing a lateral energy flow with respect to the excavation. The waves propagating across the excavation area tend to refract toward the areas of shallow water along the edges, increasing the wave heights (known as wave energy focusing). In terms of wave amplitude, there is a modest decrease downstream (in the wave propagation direction) of the excavation and a slight increase inside it and on its sides. This behavior is also confirmed with extreme waves, as shown in Figure 7. In Figure 7, the data of *Hs* refer to case 8 (Table 2) of *Hs*<sup>0</sup> = 4 m and a period *Tp* = 12 s, simulated with and without the dredging of 2 m depth.

**Figure 7.** Map of local significant wave height distribution for (*Hs*<sup>0</sup> = 4 m, *Tp* = 12 s), (**a**) before dredging, (**b**) after dredging of 2 m depth. The legend indicates the values of *Hs* in meters.

To evaluate the impact of dredging on the estuary environment, two indicators have been considered: (i) impact on sediment transport processes, where the analysis focuses on the evolution of the bottom friction velocities, responsible for the movement of sediment particles and proportional (velocity squared) to the erosion rate, and (ii) impact on structures, where the analysis focuses on the variation of the energy before the wave breaks, at each concerned zone (Figure 5). For this final point, the impact of the significant heights (squared), measured at the bathymetry of 10 m, and on each of the four selected areas was integrated. In practice, the values of *Hs* <sup>2</sup> determined at different locations SA, SB, SC, and SD (Figure 8), along the isobath of 10 m, were performed to indicate the wave energy impact on the selected zones, as shown in Figure 5. A comparison was made between relative configurations with and without dredging.

**Figure 8.** Segmentation for the analysis of energy flows in the concerned zone. The contour lines indicate the bathymetry values.

Figure 9 depicts examples of the relative bottom friction velocity (reflecting the bottom resistance) distribution in the target area. The relative friction velocity is defined as the difference between the simulated friction velocity with bathymetry with the dredged area and that without dredging (current bathymetry). The data shown in Figure 9 refer to the dredging of 2 m depth. Figure 9a shows the results of the relative friction velocity obtained with the most frequent case (*Hs*<sup>0</sup> = 1.5 m, *Tp* = 8 s) (see Table 2), whereas Figure 9b shows that of the highest wave condition (*Hs*<sup>0</sup> = 4 m, *Tp* = 12 s). Figure 9 shows the clear effects of the dredging on the bottom friction velocity distribution.

The bottom friction velocity changes significantly as the wave gets closer to the excavation site. Locally, around the borders of the excavation, a minor increase of bottom resistance is noted, as demonstrated by the positive relative friction velocity values in Figure 9. The bottom friction velocity significantly decreases inside the dredging area, reaching maximum magnitudes (of negative values). The effect of the dredging area on the friction velocities extends downstream of it, reaching zones C and D (Figure 5).

The variation of the bottom friction velocity due to the excavation certainly affects the sediment transport potential [29–31] and the wave energy flux linked to the redirected waves, particularly in zones C and D. Figure 9 indicates a clear increase in the bottom friction velocity at the lateral sides (in the wave direction) of the excavation, which is more pronounced with the high tide condition (*Hs*<sup>0</sup> = 1.5 m, *Tp* = 8 s). The possibility for erosion action increases as bottom shear stress increase. The friction velocity at the sea bottom slightly decreases in zones C and D, which are located downstream of the extraction site. This suggests that zones C and D are likely subject to sediment accumulation following dredging activities. An excessive sediment buildup can cause several environmental problems. It can reduce the seawater depth, preventing the passage of ships. It can also lead to contamination that poses a threat to aquatic plants (*Posidonia oceanica*) and wildlife. Zones A and B of the coastal area in the southern part (from breakwaters) are almost not affected by the excavation site. Additionally, findings indicate that dredging up to 2 m depth has no significant effect along the channel profile.

**Figure 9.** Relative bottom friction velocity, (**a**) (*Hs*<sup>0</sup> = 1.5 m, *Tp* = 8 s), (**b**) (*Hs*<sup>0</sup> = 4 m, *Tp* = 12 s) after dredging of 2 m depth. The legend indicates the values of relative friction velocity in [cm/s]. The contour lines indicate the bathymetry values.

Figure 10 shows the wave height, *Hs* (Figure 10a) and its relative velocity (Figure 10b) distribution maps, with additional excavations up to 4 m deep. Figure 10a demonstrates how waves propagating across the excavation area tend to refract toward the areas of shallow water along the edges, increasing the wave heights (known as wave energy focusing). This causes the amplification of *Hs* at the level of the river mouth, affecting zone B (Figure 5). The increase in wave height at the mouth of the river continues to propagate upstream along the channel. As a kind of obstruction, the excavation causes the incident waves to converge and creates a wake region downstream of it where the wave height decreases.

Similar behavior to the wave height distribution is shown in Figure 10b, where the relative friction velocity oscillates near the bottom. The friction velocities increase along the excavation's lateral edges and decrease downstream of it. Compared to Figure 9a, dredging up to 4 m deep has a greater impact on the estuary environment than dredging up to 2 m deep. With 4 m of dredging, the area is subject to increased erosion. Zones C and D, however, are subject to sediment accumulation.

Estuary ecosystems are significantly impacted by wave energy. The bathymetry in some coastal regions causes a concentration of wave energy, which raises wave height [32]. Figure 11 displays the fraction of wave relative energy percentage that corresponds to zones A, B, C, and D (Figure 5). The energy was estimated at 10 m isobath, as shown in Figure 8, and the relative percentage of energy was determined as the difference between the energy of the simulated wave with the excavation present and without it. The percentage is in relation to the simulated wave energy without excavation. The data shown in Figure 11 relates to a 2 m-deep excavation.

**Figure 10.** Impact of dredging 4 m-deep with (*Hs*<sup>0</sup> = 1.5 m, *Tp* = 8 s) (**a**) significant wave height distribution after dredging, (**b**) relative friction velocity at the seabed after dredging up to 2 m depth. The legends indicate the values of *Hs* in meters and those of relative friction velocity in [cm/s]. The contour lines indicate the bathymetry values.

Figure 11 demonstrates that the cases 1, 5, and 8 (see Table 2), with the largest values of *Hs*<sup>0</sup> (4 m), undergo the greatest values of energy variations. The condition of the lowest sea level (case 1), in the different zones from A to D, consistently exhibits the highest variation in energy. The smallest variations in energy always appear in cases 4 and 7, with the lowest values of *Hs*<sup>0</sup> (1.5 m). Figure 11 further demonstrates that zone A, which has a maximum variation value of order 4%, is the most affected zone by wave energy. The energy impact gradually decreases in the order from zone A to zone D, which is especially prominent in zones C and D. In general, it can be concluded that the dredging activities have a minimal but noticeable impact on the estuary environment of Oued Sebou that is not significant.

The simulated results of wave dispersions and their hydrodynamic structures are useful for estimating the sediment transport rates that will be an extension of this work.

**Figure 11.** Energy percentage change in zones A, B, C and D.

## **6. Conclusions**

The SWAN model was used to simulate a large number of wave motions at the Kenitra site in order to better understand the effects of dredging activities (due to sand extraction) in the Oued Sebou estuary. In various meteorological conditions, comparisons between the simulated results with and without dredging bathymetries were made. Eight configurations that accurately reflected the actual hydrodynamic circumstances that characterize the study area were tested.

The effects of the bathymetric changes, due to the sand extraction, on wave dispersion, bottom friction velocity field, and energy budget (before breaking) were extensively analyzed.

The sharp variation of the bathymetry due to the dredging that increases water depth causes an increase of the local significant wave height in and at the edges of the area susceptible to sand extraction. According to the simulation results, the waves propagating across the excavation area tend to refract toward the areas of shallow water along the edges, increasing the wave heights (known as wave energy focusing). This causes the amplification of *Hs* at the river mouth level and along the southern breakwater structure, which is more pronounced with deeper dredging. The excavation plays the role of a kind of obstacle, causing the incident waves to converge and creates a sort of wake region downstream of it where the wave height decreases.

The results show that zones C (at the northern breakwater) and D (northern coastal area) are subject to possible accumulation of sediments, whereas zones A (at the southern breakwater) and B (river mouth) are subject to an increased potential for erosive action and a risk of scouring processes at the southern breakwater.

Dredging activities in the Oued Sebou estuary mainly affect the river mouth (zone B) and the southern breakwater area (zone A), which is very noticeable with dredging up to 4 m deep.

In general, it can be concluded that the dredging activities show a certain level of impact on the estuary environment of Oued Sebou that is not very significant.

Despite the crucial role dredging plays in the nation's economy and maritime engineering management (i.e., it helps make the water navigable, removes contaminants from seabeds and recreates damaged areas, maintains many marine infrastructures, and many other advantages), dredging could have serious and long-lasting negative impacts on the environment, leading to contamination that poses a threat to aquatic plants (*Posidonia oceanica*) and wildlife

The simulation results, which will be validated by some measured wave characteristic data, are useful for examining sediment variations along the estuary coastal area, a subject we are currently working on.

**Author Contributions:** N.I. and M.C., performed the numerical simulations and methodology; M.B.M. and N.I. formal analysis, study design, writing–original draft preparation; L.M., M.C. and D.L. contributed suggestions, discussions, and reviewed the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are available on request from the first author (Nisrine Iouzzi).

**Acknowledgments:** This study was carried out at Hassan II University of Casablanca, Faculty of Sciences Ben M'Sik, (Morocco).

**Conflicts of Interest:** The authors declare that they have no conflict of interest.

## **References**

