Large eddy simulation of turbidity currents in a narrow channel with different obstacle configurations

Turbidity currents are frequently observed in natural and man-made environments, with the potential of adversely impacting the performance and functionality of hydraulic structures through sedimentation and reduction in storage capacity and an increased erosion. Construction of obstacles upstream of hydraulic structures is a common method of tackling adverse effects of turbidity currents. This paper numerically investigates the impacts of obstacle’s height and geometrical shape on the settling of sediments and hydrodynamics of turbidity currents in a narrow channel. A robust numerical model based on LES method was developed and successfully validated against physical modelling measurements. This study modelled the effects of discretization of particles size distribution on sediment deposition and propagation in the channel. Two obstacles geometry including rectangle and triangle were studied with varying heights of 0.06, 0.10 and 0.15 m. The results show that increasing the obstacle height will reduce the magnitude of dense current velocity and sediment transport in narrow channels. It was also observed that the rectangular obstacles have more pronounced effects on obstructing the flow of turbidity current, leading to an increase in the sediment deposition and mitigating the impacts of turbidity currents.

www.nature.com/scientificreports/ particles deposition downstream of the downslope and on the upslope of the humps 30 . Oshaghi et al. (2013) demonstrated that the obstacle height and the upstream velocity of turbidity currents are inversely related 19 . Khavasi et al. (2012 and) studied the effects of particle size, bed slope and inlet Froude number on the stability of turbidity currents. They demonstrated that an increase in the particles size, bed slope and inlet Froude number can diminish the stability of the dense flow regime [31][32][33][34] . Numerical approaches have also been used to study the dynamics of turbidity currents. Toniolo et al. (2007) developed a numerical framework to predict the trapping efficiency of turbidity currents in reservoirs and showed the impacts of topology on the reduction of fine particles settling efficiency 35 . Oehy et al. (2007) compared solid and porous obstacles confronting turbidity currents and showed a slight reduction in the trapping efficiency for porous obstacles 36 .
Several turbulence models and simulation approaches have been used to study turbidity currents including Direct Numerical Simulation (DNS), Large Eddy Simulation (LES) and Reynolds-averaged Navier-Stokes equations (RANS). RANS turbulence models are less computationally expensive in comparison with LES and DNS. RANS models with k-ε turbulence closure have been used in several studies of turbidity currents [37][38][39][40][41] . However, RANS models usually fail to accurately resolve flow zones with intense shear (near walls or obstacles) and flow of low to moderate Reynolds number 42,43 . Additionally, the numerical constants in RANS models need careful tuning procedures for the specific flow conditions in order to improve the accuracy of the solution, given that the Reynolds stresses in the RANS equations depend on the boundary and flow conditions. In the LES models, filtered-out eddies are not influenced by the flow conditions as the governing equations are derived based on the physical properties of the flow 37,44 . This study, for the first time, develops a numerical simulation framework using an LES turbulence model to investigate turbidity currents confronting obstacles of various geometrical configurations in a narrow channel. Also, for the first time, this study investigates the effects of discretization of particles size distribution on sediment deposition and propagation in narrow channels. The flow hydrodynamics and sediment concentration of turbidity currents over two obstacles of varying geometrical configurations in a narrow channel are investigated using the validated method described in this paper. This study highlights the capabilities of the LES numerical approaches for robust and accurate prediction of turbidity currents.

the model
A dense current occurs when a dense fluid propagates into a lighter fluid. The dense current is propagated under the combined influence of its initial momentum and the gravity body force 1 (Fig. 1). A lab-scale narrow channel containing freshwater is considered to investigate the behavior of turbidity currents. The choice of the narrow channel in this study is to characterize the augmented effects of the shear stress caused by side walls. A dense current is released into the channel and the interactions of the dense current with the fresh water (ambient) flow is simulated using LES model described in below.
Fluid flow is governed by the Navier-Stokes equations including continuity and momentum conservation equations. The concentration of particles in the dense current is modelled with a transport equation. The density difference between dense and light (ambient) fluids is assumed to be sufficiently low so that the Boussinesq approximation remains valid for modelling buoyancy forces. The gravity-buoyancy term in the momentum equation is defined as 45 : where C and C 0 represent the normalized particle concentrations [unit less] at density ρ and ρ 0 [kg/m 3 ] for the dense and ambient fluid, respectively. The volumetric coefficient of expansion for the particles is β = 1 [dimensionless] and g represents the gravitational acceleration [m/s 2 ]. In this study,C varies from C 0 = 0 to C max = 1 . To conduct a Large Eddy Simulation (LES), the continuity and momentum conservation equations are defined as:   44,49,[51][52][53][54] . In this study, to maintain reasonable computational costs, the Schmidt number is assumed to be Sc = 1.
The particles inertia forces and particle-particle interactions are not computed as the concentration of suspended solids is relatively low 52,55 . Therefore, particle's transport is simultaneously governed by the flow hydrodynamic and Stokes' settling velocity 56 : where µ is the dynamic viscosity [kg/m.s], d p denotes the diameter of particles and ρ p is the density of particles [kg/m 3 ]. Ten different particle size intervals ranging from 0.5-100 μm are considered to represent the particles with an identical density ρ p . For each particle size, the Eulerian continuum transport equation is implemented according to Eq. (9) 57 , with a constant value for the Stokes' settling velocity U s : Effects of the filtered fluctuations of the flow hydrodynamic are considered by the momentum and concentration residual-stress tensors τ ij and τ C i : Turbulence effects are modelled using Smagorinsky closure model 58 , including a SGS eddy-viscosity υ SGS model to compute the residual tensors: where S ij is the strain tensor, Smagorinsky coefficient is taken as C s = 0.2 , denotes the filtered width and turbulent Schmidt number (which is from the order of unity 59 61 . The numerical flume is then utilized to investigate the hydrodynamic behaviour of turbidity current on a sloping bed subjected to different obstacles configurations. The channel bed is assumed to be smooth with a slope of 1% (Fig. 2). A continuous dense current introduced into the channel, with a constant particle density ρ = 2649 kg/m 3 and a mean diameter D 50 = 11µm. where U 0 is the mean velocity of the turbidity current at the inlet, H inlet is the height of the inlet and the bottom slope is θ . The reduced gravity acceleration is determined as 61 : where ρ and ρ 0 are density of the turbidity current and the ambient fluid, respectively.
Six simulation scenarios are designed to determine the effects of obstacles height and geometrical shape on the behavior of turbidity currents. A summary of simulation cases is shown in Table 1.
To guarantee an appropriate determination of the turbulent boundary layer, the first grid cell adjacent to the solid boundary is resided in the viscous sub-layer (y + < 5) for a robust resolve of the boundary layer 50 . The location of this grid cell in plus units and associated parameters are described by Eqs. (17)(18)(19)(20) 62,63 :  www.nature.com/scientificreports/ where y + is the dimensionless distance from the wall, U 1 represents the velocity at the first cell, y is the distance of the first cell from the solid boundary, τ w and U τ are wall shear-stress [kg/m⋅s 2 ] and associated friction velocity, respectively 64 . The height of the first cell was determined using the well-established empirical correlations described by Eqs. (21)(22)(23) 50,62,63 : where c f is the wall skin friction coefficient, Re x denotes the Reynolds number based on the boundary layer thickness, U x is the inlet velocity [m/s] and c f is an empirical constant computed based on Reynolds number described by Eq. (23).

Model verification
The numerical model was developed using finite-volume technique and computational codes were written in C++ with OpenFOAM (V6) open-source license. Second-order limited linear scheme was implemented for discretizing governing equations, except for the transport equation where a second-order QUICK scheme was adopted. The numerical robustness and accuracy of the QUICK scheme have been demonstrated by previous related studies [46][47][48]51 . The Pressure Implicit with Splitting of Operators (PISO) algorithm was implemented to solve the filtered LES equations in a transient mode. A two-steps corrector was considered in the PISO algorithm to guarantee computational robustness and a better convergence (i.e. pressure equation is corrected two times per time-step to satisfy the continuity equation) [65][66][67] . A series of numerical simulations were conducted to verify the numerical method and the developed model. In turbidity currents, where complex phase-coupling (momentum exchange) between particles and fluid exists, velocity profiles determines the sediment deposition and the characteristics of flow hydrodynamic. A common method to evaluate the performance of numerical methods and computational codes is comparing the numerical results with physical modelling measurements 36,47,[68][69][70][71][72][73][74] . Comparison of vertical variations of the numerical velocity profiles with the experimental measurements of Farizan et al. (2018) was conducted to validate the numerical model described in 61 .
The temporally-averaged velocity profiles at 0.5 m before the obstacle were obtained once the steady-state condition is reached, and then are compared with the experimental measurements of Farizan et al. 61 . For the validation case, channel (L:12 m, W: 0.2 m and H: 0.6 m) with a triangle obstacle located at 4.5 m away from the inlet was considered. The inlet geometry has the same width as the channel (= 0.2 m) with the height of 0.07 m. A fully developed flow condition is applied at the inlet for the turbidity current entering into the channel to include turbulent flow fluctuations. Three grid resolutions with 3.900, 4.350 and 5.590 million structured hexahedron cells (namely Mesh 1, Mesh 2 and Mesh 3, respectively) were considered along with a mean particle diameter of D 50 = 11µm , to conduct sensitivity analysis and mesh dependency study. Time-averaged velocity profiles at a distance of 4 m from the inlet were determined for the validation test cases (Mesh 1, Mesh 2 and Mesh 3) with one concentration transport equation, and compared with the laboratory measurements of Farizan et al. 61 (Fig. 3). The comparison of the results presented in Fig. 3 highlights considerable deviations between experimental and numerical velocity profiles. The numerical model with one concentration transport equation overestimates the shear effects in the dense current and the velocity profile near the bed (the lower part of the velocity profile). It was shown that by ignoring the particles larger than D 50 , the velocity of turbidity current near the bed was increased.
To improve the discrepancy between numerical results and the measurements, additional concentration transport equations were considered. Additional simulation cases with two, five and ten concentration transport equations were conducted to produce a more robust estimation of the particles size and distribution (0.5-100 μm). Figure 4 shows the distribution of particles experimentally measured by Farizan et al. (2018) 61 . Figure 5 compares the velocity profiles from numerical simulations with the experimental measurements of Farizan et al. (2018) 61 . The simulation sets presented in Fig. 5 include extra particle size intervals and more transport equations to enhance the computational accuracy of the velocity profiles. The simulations show that for all the cases, the velocity in the upper region of the turbidity current (see Fig. 6) is increased as finer particles are introduced to the flow, while the lower region of the turbidity current (near the channel bed) is slowed down due to the effects of larger particles. The accuracy of the numerical results is improved by increasing the number of concentration transport equations and the best performance was achieved for the case of 10 concentration transport equations. Following grid dependency analysis, to achieve numerical stability, computational Despite the shear stress effects of walls on the fluid might initially seem to be significant, in this narrow channel, the results showed that at the Reynolds number considered in this study the walls effect is negligible, which is in good agreement with the study by Khavasi et al. and Oehy et al. 31,32,36,68 .

Results and discussion
The velocity and concentration profiles inside the dense layer of the turbidity current are categorized into three distinct regions (Fig. 6): (1) the upper part is known as a shear layer region where the density of the turbidity current decreases and asymptotes to the counterpart value for the ambient fluid; (2) the middle part is known as the suspension zone where the majority of particles are suspended in the fluid; and (3) the lower zone which is a depositional area where particles are settled 31 . In this paper, the upper zone of the channel (part 1 in Fig. 6) is not described as the flow velocity asymptotes to zero 61 .
Velocity and concentration values are measured after reaching a quasi steady-state to avoid temporal fluctuations in the flow parameters.  where z upper is the upper boundary over which the concentration becomes negligible, u(z) is the sum of settling and horizontal velocity and c(z) is the concentration for the dense current.
The turbidity bore is defined as a moving hydraulic jump over the bed of the channel. Obstacles alter the flow regime and can move the internal bores towards the inlet of the channel. Also, the flow hydrodynamic characteristics significantly impact the position and structure of turbidity currents on the channel bed 36,75 . Figures 7  and 8 show temporal evolution of the turbidity current in the channel, indicating multiple reflected bores of the turbidity current as it travels inside the channel and over the obstacle. The horizontal velocity of turbidity current slightly decreases when the flow climbs up the obstacle which is due to the flow-obstacle interactions and the consequent dissipation of turbulent kinetic energy of the current (Figs. 7, 8). The bore of the dense flow becomes thinner immediately after it passes the obstacle and as it moves towards the outlet. For the case of both obstacles (Figs. 7, 8) sediment deposition and formation of an interface between turbidity current and the fluid of lighter density is observed 36 .
Simulations were continued until the sediment deposition behind the obstacle reached a steady height. Temporally-averaged flow characteristics were determined once the quasi steady-state condition is met, to avoid temporal fluctuations of the LES.
The comparison between Figs. 7 and 8 at t = 125 s demonstrates an intensified turbulent hydraulic jump when the turbidity current passes over the rectangular obstacle. However, the flow over the triangular obstacle can be characterize as a quasi-uniform jump with a less disturbed flow.
Turbidity currents can significantly be influenced by the height of the obstacle. Previous studies show the effects of the obstacle on blockage and reflection of turbidity currents 36,61 . The obstacle's height of equal to twice the height of the current, was reported to cause a considerable reflection in the turbidity current flow 76 .
Considering different heights, the stream-wise time-averaged velocity profiles of turbidity current in the quasi steady-state were plotted at 0.5 m upstream of both rectangular and triangular obstacles (Fig. 9). The increase in the height of the obstacles considerably changed the vertical structure of the dense current's velocity profiles. For the case of obstacle height of 0.06 m and for both geometries, the maximum velocity was observed in the settling zone behind the obstacle. For the cases with the obstacle height of 0.10 m, the velocity profile shows a sharp increase from the bottom up to the depth of 0.10 m, then a sharp reduction is seen up to the interface between the dense and ambient fluid. For the cases with the obstacle height of 0.15 m, the velocity profiles and vertical distribution of shear effects have smaller values with a less distributed pattern in comparison to the cases with smaller obstacle height. The patterns of the velocity profiles for both obstacles geometries are very similar, with the maximum velocities for the rectangular obstacle occurring at slightly higher depths from the bed. Rapid changes in the velocity profile of the turbidity current for the rectangular obstacle begin in higher elevation in contrast with triangular obstacle, which is due to the higher sediment decomposition for rectangular obstacle case. Figure 10 shows www.nature.com/scientificreports/ turbidity current behind the rectangular obstacle in comparison with the triangular obstacle, highlighting the effects of obstacle's geometrical shape on the control and mitigation of turbidity currents. Figure 11 shows the suspended sediments flux per unit width of the channel for all the simulation cases. Implementation of the obstacles improved the settling efficiency for the turbidity current upstream of the obstacles. Increasing the height of the obstacles slowed down the vertical variation of horizontal velocity profile which led to a reduction in the sediment flux. The presence of obstacle, regardless of its geometrical shape, had no considerable impact on the settling rate at the downstream of the channel. The settling deposition of particles for the cases with 0.10 m and 0.15 m obstacles in the upstream of the obstacles are almost equal.
The effects of the obstacle's geometrical shape on the sediment deposition was more pronounced for the rectangular obstacles, mainly due to the higher reduction in the dense flow velocity behind the rectangular obstacle. In order to compare the obstacles shape impact, the difference in sedimentation flux-rate (q s ) for triangular and rectangular obstacles are determined ( Table 2). Positive values indicate a dominancy in the settling of suspended sediments for the rectangular obstacle. Accordingly, rectangular obstacles are suggested to be implemented in channels leading to hydraulic control structures. More deposition of the sediments increases the efficiency of hydraulic structures in water systems 44 .
(25) q s_variation = ( q s(triangularobstacle) − q s(rectangularobstacle) q s(triangularobstacle) ) × 100   www.nature.com/scientificreports/ www.nature.com/scientificreports/ conclusion Appropriate understanding and analysis of turbidity currents are vital for sustainable and efficient management and operation of natural and man-made hydraulic structures. This study develops and successfully validates a high-resolution numerical simulation model using mathematical capabilities of Large Eddy Simulation (LES) technique. The effects of the number of concentration transport equations on the robustness and numerical accuracy were studied in detail. The results highlight that discretization of the particles size distribution improves the accuracy of the model in predicting turbidity current hydrodynamics and spatiotemporal structure of turbulence. The effects of obstacle's geometrical shape and height on the turbidity currents characteristics in a narrow three-dimensional channel were modelled. Two obstacle prototypes of rectangular and triangular shape with varying height were investigated. The numerical velocity and concentration profiles were determined for all simulation cases described in Table1. The findings indicate that, for both rectangular and triangular obstacles, by increasing the height of the obstacle, the maximum velocity of the turbidity current was reduced and the shape of the vertical distribution of flow hydrodynamic in the dense layer was changed. The results show that the increased height of obstacle directly impacted the vertical structure of shear and turbulent velocity. Across all test cases, comparison between the two obstacle geometries shows that for both obstacles the overall shapes of flow hydrodynamics are similar.
Furthermore, a direct relationship between the obstacle's height and the settling capability of the obstacles was observed. The numerical results highlight that the shape and height of obstacles significantly change the hydrodynamics, sediment particle distributions and structure of turbidity currents over a smooth bed channel. Installation of a rectangular obstacle is recommended to enhance the deposition and efficiency of hydraulic structures (dams, reservoirs and weirs) in water systems.
The computational framework developed in this study demonstrates that LES modelling can be implemented as computationally robust and reliable numerical technique to investigate the dynamics of turbidity currents in turbulent flow conditions. www.nature.com/scientificreports/