Modeling of Breaching-Generated Turbidity Currents Using Large Eddy Simulation

: Breaching ﬂow slides result in a turbidity current running over and directly interacting with the eroding, submarine slope surface, thereby promoting further sediment erosion. The investigation and understanding of this current are crucial, as it is the main parameter inﬂuencing the failure evolution and fate of sediment during the breaching phenomenon. In contrast to previous numerical studies dealing with this speciﬁc type of turbidity currents, we present a 3D numerical model that simulates the ﬂow structure and hydrodynamics of breaching-generated turbidity currents. The turbulent behavior in the model is captured by large eddy simulation (LES). We present a set of numerical simulations that reproduce particular, previously published experimental results. Through these simulations, we show the validity, applicability, and advantage of the proposed numerical model for the investigation of the ﬂow characteristics. The principal characteristics of the turbidity current are reproduced well, apart from the layer thickness. We also propose a breaching erosion model and validate it using the same series of experimental data. Quite good agreement is observed between the experimental data and the computed erosion rates. The numerical results conﬁrm that breaching-generated turbidity currents are self-accelerating and indicate that they evolve in a self-similar manner.


Introduction
Turbidity currents are buoyancy-driven underflows generated by the action of gravity on the density difference between a fluid-sediment mixture and the ambient fluid. The excess hydrostatic pressure within the turbidity current drives the current downstream while complicated interactions with the surrounding environment take place; it interacts with the ambient fluid at the upper boundary and with the bed at the lower boundary, producing turbulence at both boundaries [1]. Turbidity currents are vital agents of sediment transport that deliver sediment from the river mouths to deeper waters [2]. Moreover, they pose a serious threat to flood defense structures, such as dikes and embankments [3], and submarine structures placed at the seafloor, such as oil pipelines, and communication cables [4].
Turbidity currents can be triggered through several mechanisms, such as hypopycnal river plumes [5], internal waves or tides [6], and submarine slope failures [7]. One of the complex failure mechanisms of submarine slopes is the flow slide [8], which takes place when a large amount of sediments in an underwater slope is destabilized and consequently runs down the slope as a dense fluid. Two categories are distinguished: liquefaction flow slides, which occur in loosely packed sand, and breaching flow slides, which mostly occur in densely packed sand [9,10]. The former results in slurry-like or hyper-concentrated density flows, while the latter results in turbidity currents [11]. Here, the focus is on the turbidity current accompanying breaching flow slides, referred to as a breaching-generated turbidity current. These currents have proved very difficult to measure in the field, as their occurrence is unpredictable while they can also destroy the measuring instruments [12,13]. Alternatively, laboratory experiments are widely used (e.g., [14][15][16][17]) and can be scaled and conducted under more controlled conditions to develop a better understanding of the behavior of turbidity currents, and to provide measurements for the validation of numerical models.
To this end, Alhaddad et al. [18] have recently conducted large-scale experiments, obtaining direct measurements of breaching-generated turbidity currents and the associated sediment transport. Their analysis showed that breaching-generated turbidity currents are self-accelerating; the current strengthens itself by the accumulated erosion of sediment from the breach face. The results also suggested that the velocity profiles of these currents are self-similar. Analysis of particle concentration profiles showed that the concentration decays exponentially from the breach face until the upper boundary of the current. Near-bed concentrations were found to be high, reaching 0.4 by volume or even higher in some cases. Moreover, investigation of the slope failure indicated that its evolution is largely three-dimensional. Sand erosion rates in the middle of the tank width, where turbidity currents were measured, were found to be higher than at the tank glass wall, where the erosion rates were measured. A key finding was that the turbidity current is the main parameter controlling the evolution of the breaching failure and the fate of eroded sediment. This implies that a thorough understanding of the behavior of this current is needed to enhance the knowledge about breaching. Owing to the several difficulties associated with breaching experiments, measurements of turbulence quantities of the flow were not possible. The absence of such measurements hinders the estimate of the flow-induced bed shear stress and hence the predication of erosion during breaching. This highlights the need for using advanced numerical models as a complementary tool to the experimental work, to gain new insights into the behavior and structure of breaching-generated turbidity currents.
Many numerical models for density or gravity currents have been proposed in the literature (e.g., [19][20][21][22][23]). To date, however, there are very few numerical investigations dealing with breaching-generated turbidity currents [3,11,24]. Moreover, these numerical efforts were mostly restricted to layer-averaged, one-dimensional models. The model of [24] was applied to simulate a flushing event in Scripps Submarine Canyon, and showed that breaching-generated turbidity currents can excavate a submarine canyon. Similarly, the model of [11] was applied to a flushing event in Monterey submarine canyon. These depth-averaged models require several empirical closure relations (e.g., the near-bed concentration, bed shear stress, and water entrainment at the upper boundary), which reduces the accuracy of the simulation results [25]. Alhaddad et al. [3] applied the one-dimensional model of [26] to a typical case of a breaching slope, demonstrating that the results are highly sensitive to the type of breaching closure relation used.
To reduce these uncertainties, this study presents large eddy simulations of breaching-generated turbidity currents. Large eddy simulations have the advantage that the larger turbulent scales-containing the bulk of the turbulent kinetic energy-are resolved. In this manner, the influence of density gradients on turbulence production is captured adequately, while including the non-isotropic character of turbulence. Therefore, large eddy simulations can provide a wealth of insights into the structure and physical behavior of turbidity currents, in particular into the turbulence structure including turbulent kinetic energy and Reynolds stresses. The numerical model that we use has been applied earlier to gain insights into the complicated dredge plume behavior close to a dredging vessel where density differences, turbulent mixing and sediment settling play an important role [27][28][29]. The model has been validated for a wide range of flow cases relevant to the present study, such as the front speed of density currents radially spreading and density currents running over both flat and sloping beds, the deposition from high-concentration suspended sediment flow at a flat bed [30], and high sediment concentration conditions encountered in hopper sedimentation [31]. In this paper, the implementation is designed to capture the turbidity current running down the slope surface (or 'breach face'), considering various steep slope inclinations, which were tested in the laboratory experiments. The triggering mechanism of turbidity currents in this work is different from the standard lock-exchange configuration mostly used in the numerical models; sediment particles are released from the bed surface, generating the flow. An adequate prediction of this process has been always difficult, as it involves both geotechnical and hydraulic features [3]. To date, furthermore, no validation of existing breaching erosion models has been presented.
The paper proceeds as follows. We first provide some background of the numerical model in Section 2, and propose a breaching erosion closure model in Section 3. Section 4 revisits the experimental measurements of [18], and validates the performance and limitations of the currently proposed numerical model in characterizing the breaching-generated turbidity currents based on the experimental findings. Following that, Section 5 discusses the flow and turbulence structure and analyses the sensitivity of the numerical results to some initial conditions.

Governing Equations
The flow of water-sediment mixtures-as in breaching-generated turbidity currents-is mathematically described by the incompressible, variable density Navier-Stokes equations. Using the mixture approach, the concentrations of the individual sediment fractions are solved separately while one set of momentum equations is solved for the water-sediment mixture [32]. Each sediment fraction has its own drift velocity u d , defined as the velocity of the sediment relative to the water-sediment mixture, which involves a correction of the mixture velocity with the settling velocity of the sediment including the influence of hindered settling. Furthermore, the sediment concentration influences the mixture density. The mixture approach combined with the drift velocity is thus two-way coupling and has proved valid to simulate high-concentration suspended sediment flows (e.g., [31,[33][34][35]).
To this end, the corresponding balance equations for the total mass and momentum of the mixture are stated as, respectively, where t denotes time, ρ is the mixture density, ρ w is the water density, u is the mixture velocity vector, P is the modified pressure, τ is the shear stress tensor, g is the gravity vector, and ∇ denotes the spatial gradient operator. The mixtures that we consider are primarily Newtonian in rheology [36] for which the shear stress tensor is proportional to the deviatoric strain rate of the mixture velocity, as follows, where µ is the dynamic viscosity of the mixture and ∇ s (·) = 1 2 ∇(·) + 1 2 ∇(·) T is the symmetric gradient operator. It is to be noted that Equation (3) is valid only if the layer-averaged volumetric sediment concentration is below the Bagnold limit of ca. 9% beyond which the particle-particle interactions will render the mixture non-Newtonian [37]. In our case, observed layer-averaged concentrations are generally below 9% (see Section 6.1), while in the center and outer region of the turbidity current, concentrations are even smaller due to the entrainment of ambient fluid and the subsequent turbulence mixing [38]. This motivates the use of Equation (3), although in the inner region near the bed the relatively high concentration may require the reconsideration of additional effects on the rheology of the fluid in the boundary layer, leaving room for improvement of our model (see also the discussion in Section 5.3). Although the mixture model admits multiple sediment fractions [30], a single sediment fraction suffices here, as uniform sediment will be considered. In this case the mixture density is given by, where ρ s is the density of sediment, and c is the volumetric sediment concentration. The latter satisfies the following transport equation, in which u s = u + u d is the velocity of the sediment fraction, and Γ is the diffusivity which is related to the viscosity µ by, where Sc is the Schmidt number. Following [30,33], Sc = 0.7 is adopted in this study and, based on sensitivity analysis, LES results are insensitive to the value of Schmidt number when Sc is close to 1. The drift velocity u d accounts for the effects of hindered settling [39] and the return flow created by the settling particles [30]. For a well-posed problem, Equations (1), (2) and (5) must be supplemented with boundary conditions, while for time dependent problems also initial conditions for the sediment concentration and mixture velocity must be prescribed. Importantly, the interaction of the mixture flow with a sediment bed involves the bed shear stress τ b , and the densimetric sedimentation flux S and erosion flux E given by, where c b is a representative near-bed sediment concentration, w s is the particle settling velocity, α is the bed slope, ∆ = (ρ s − ρ w )/ρ w is the relative submerged sediment density, d 50 is the median sediment grain size, and φ p,t is a non-dimensional, so-called, pick-up function involving an empirical closure relation. Given the strong feedback between bed erosion and the hydrodynamics of the turbidity current, the formulation of the pick-up function requires special care to accurately predict the temporal evolution of the slope failure. Section 3.3 will further elaborate on the available theoretical background of sediment erosion during breaching, resulting in the proposed erosion closure model.

Turbulence Modeling
In our implementation, Equations (1), (2) and (5) are discretized using a regular rectangular grid with, possibly non-uniform, grid sizes (∆x, ∆y, ∆z in the respective Cartesian directions). Depending on the grid size, the finite resolution of the computed flow field can only partly include the relevant turbulence length scales. To account for the effect of the unresolved scales on the resolved flow field a turbulence closure model must be used.
In this study, Large Eddy Simulation (LES) is adopted to capture the influence of turbulence. LES applies a spatial filter in which all fluctuations smaller than the grid size are replaced by a sub-grid-scale contribution. In this way, turbulence length scales larger than the grid size are resolved, and the smaller isotropic turbulence scales are accounted for by the sub-grid-scale model. By choosing the grid size sufficiently small, the major part of the turbulent kinetic energy is resolved on the grid and only a small part is modeled by the sub-grid-scale model. An advantage of LES over Reynolds-averaged Navier-Stokes modeling (RANS) is that turbulence damping or destruction functions at sharp density gradients are not needed when sufficient grid resolution is used. This is because the influence of density gradients on the resolved turbulent eddies is automatically taken into account in LES. Furthermore, the non-isotropic behavior of the larger turbulence length scales is resolved in LES. Noteworthy to mention in this respect is that the behavior of the inner region of the turbidity currents, as observed in the experiments, is highly affected by both damping and anisotropy of the turbulence motion in this region.
Following the LES approach, the molecular viscosity µ mol is enhanced with an extra sub-grid-scale viscosity µ sgs , as follows, where µ sgs is a function of the strain rate tensor and the grid size, in which C S is the dimensionless Smagorinsky constant, and S is the norm of the (filtered) velocity gradient tensor. For the latter, the WALE (wall adapting eddy viscosity) sub-grid-scale turbulence model is adopted together with a constant C S of 0.325. For details of this implementation see [40]. Turbulence is partly generated at the sediment bed, requiring a formulation of the corresponding boundary conditions consistent with the LES approach. First, the bed shear stress is formulated as a partial slip boundary condition by calculating the local shear velocity u * . This is accomplished by assuming a standard logarithmic velocity profile over the grid cell adjacent to the bed, which gives, where u t is the (filtered) tangential velocity vector in the grid cell adjacent to the bed, κ is the Von Kármán constant, ∆Z is the cell size normal to the bed, and Z 0 = 0.11ν/|u * | + k n /30, in which ν is the kinematic molecular viscosity of water, and k n is the Nikuradse roughness height which is set to a value of 2 d 50 . Next, the bed shear stress, τ b , is computed from, The corresponding magnitude of the shear velocity, u * , is also used in the formulation of the sediment erosion flux E at the bed, which is treated extensively in Section 3.

Numerical Solution Procedure
The spatial discretization of the model equations is based on the finite volume method, using a rectangular staggered grid in which the discrete velocity and pressure variables are located at alternating points [41]. The discretization in time conforms to a pressure-correction algorithm which involves a predictor step, in which an intermediate velocity field is computed using the pressure from the previous time step, followed by a corrector step, where the velocity and pressure are updated in a coupled fashion in order to satisfy the continuity constraint, Equation (1).
In the predictor step an explicit third order Adams-Bashfort time integration scheme is employed, adjusted to be able to apply variable time step sizes. Small time steps are to be used with the Courant-Friedrichs-Lewy number (CFL-number) staying below 0.6. The spatial discretization of the convection terms in the momentum equations is crucial for LES and is performed by a stable scheme with very low diffusion [42]. Likewise, the advection term in the transport equation is discretized with an accurate and robust TVD (Total Variation Diminishing) scheme employing the Van Leer limiter.
In the corrector step, the enforcement of the continuity constraint results in a pressure Poisson equation which is solved by a fast direct solver using the methodology of [43]. Although for incompressible single-phase flows the continuity constraint implies zero divergence of the velocity field (by setting Dρ/Dt = 0 in Equation (1)), this is generally not so for the flow of incompressible mixtures [44]. This owes to the definition of the mixture velocity u in a densimetric manner, while zero divergence only holds for the volumetric mixture velocity. The latter is obtained by correcting the densimetric mixture velocity u with the respective equilibrium settling velocities from all sediment fractions-a necessary step preceding the enforcement of the zero divergence constraint.
Overall, the numerical method is second-order accurate in time and space. For more details on the numerics, the reader is referred to [30].

Breaching Erosion Modeling
Breaching can be defined as a gradual, slow, retrogressive erosion of a steep immersed slope, which is steeper than the internal friction angle of the granular material forming that slope. As noted earlier, breaching is largely encountered in densely packed sand, as it exhibits a dilative behavior, when it is subjected to shear forces [9,45]. Dilatancy refers to the expansion of pore volume of sand under shear deformation, which results in the build-up of a negative pore pressure, with reference to the hydrostatic pressure. This negative pressure holds sand particles together and increases the effective stress [46]. An inward hydraulic gradient is developed, as a result of the pressure difference, compelling the neighboring water to flow into the sand pores, and thus releasing the negative pressure. Consequently, sand particles located at the top of the slope surface (breach face) are exclusively undermined and slowly (∼mm/s) peel off, predominantly one by one [3].
The falling sand particles mix with the ambient water, producing a turbidity current running along and interacting with the breach face, and then flowing down the slope toe. This causes an additional shear stress along the breach face, and hence higher erosion. In the conventional sediment pick-up functions, it is assumed that it is impossible for a grain-by-grain sediment erosion to take place in a submerged slope steeper than the internal friction angle. Rather, the erosion may occur as a sudden collapse of the sand body. Nevertheless, breaching refutes this hypothesis [18,47], showing the need for an erosion model compatible with breaching conditions. It is to be noted that beside the grain-by-grain erosion, an intermittent collapse of coherent sand wedges, termed a surficial slide, was observed in breaching experiments (e.g., [18,47]). The current understanding of these slides remains very limited. In this paper, therefore, we will consider measurements where surficial slides did not take place. This implies that the total erosion will be a combination of particle-by-particle erosion induced by gravity (pure breaching) and sediment erosion induced by the flow motion. In the following, breaching erosion is decomposed into these two main components.

Pure Breaching
The Dutch industry was the first to explore breaching in the 1970s, and used it as an efficient production mechanism for stationary suction dredgers. In that period, breaching was not known as a failure mechanism of underwater slopes outside the dredging field. To estimate the dredging production, Breusers [48] developed a formula for pure breaching: particle-by-particle erosion induced by gravity. The original formula was derived for a vertical breach face; however, it can be adapted to a general form representing the erosion velocity perpendicular to the breach face for variable inclinations [3]: where v e,g [m/s] is the erosion rate of pure breaching, φ [ • ] is the internal friction angle, n 0 [-] is the in situ porosity of the sand, k l [m/s] is the sand permeability at the loose state, ρ w [kg/m 3 ] is the density of water, and δn = (n l − n 0 )/(1 − n l ) is the relative change in porosity, in which n l [-] is the maximum porosity of the sand. Pure breaching is particularly sensitive to the magnitude of sand permeability k l with a linearly proportional relation. This implies that the existence of finer particles within the sand considerably decreases pure breaching, since they fill pore spaces and reduce permeability. Furthermore, the permeability plays a role in the fluid-induced erosion, as it will be shown in the next subsection.
Fortunately, the value of the permeability reported in [18] was measured in the lab, partly leading to a proper validation of the current erosion model.
Alhaddad et al. [18] showed that Equation (13) somewhat overestimates pure breaching. Therefore, we propose here an empirical correction factor of sin(α − φ) 0.55 , which leads to the expression of the corrected pure breaching v e,g,c : Direct measurements of different grain sizes are needed to test the general applicability of this correction factor. Figure 1 depicts the performance of the original and corrected expression of pure breaching. Equation (14) will be used in numerical runs for which no measured pure breaching is available.

Flow-Induced Erosion
Many parameters play a role in sand erosion induced by turbidity currents, such as near-bed velocity gradient, flow turbulent energy, the properties of sand grains and the bulk properties of sand. In breaching, this part of erosion is more complicated than pure breaching, owing to the special conditions of the breaching process including dilatancy-retarded erosion and a very steep bed [3]. On the one hand, some pick-up functions are proposed in the literature to account for the dilative behavior (e.g., [49]). However, these functions do not account for a sloping bed. On the other hand, some pick-up functions account for the sloping bed (e.g., [50]), but not for a slope steeper than the internal friction angle.
The literature shows, to the best of our knowledge, that only two erosion models were suggested to suit the conditions of the breaching problem [24,45]. These two erosion models are an extension of the formula of [48], meaning that they combine both the pure breaching and sediment erosion by the turbidity current. However, the lack of quantitative data on breaching flow slides has resulted in there being no validation of these erosion models under breaching conditions. Moreover, [3] showed that the erosion rate predicted under the same conditions varies considerably between these erosion models.

Total Erosion
In this work, we adopt the erosion model of [24] as a basal point and develop it further to improve its predictive ability of breaching erosion. Their erosion model reads where v e,t is the total erosion velocity including pure breaching and flow-induced erosion, u ss = ∆gD 50 is Shields velocity for sand grains, and φ p, f is an empirical non-dimensional pick-up function, which does not account for the bed dilative behavior nor the sloping bed: where E = ρ s v e, f (1 − n 0 ) [kg/(sm 2 )] is the sediment pick-up rate perpendicular to the bed surface in which v e, f is the velocity of fluid-induced erosion. The general solution of Equation (15) was not reported in [24]. Instead, two solutions for the two extreme cases v e,t /v e,g << 1 and v e,t /v e,g >> 1 were provided. The first extreme case is never encountered in breaching, while the second extreme case does not hold under lab conditions. Alternatively, Equation (15) can be rearranged into a quadratic equation, and after substituting Equation (13) into the resulting quadratic equation, the general solution will read: where f = 1 if Equation (13) is used for pure breaching, whereas f = v e,g,used /v e,g when another value of v e,g,used is used for pure breaching. An important feature of breaching-generated turbidity currents is their high particle concentration, the effect of which should be accounted for in the breaching erosion model. High near-bed concentrations reduce the flow-induced sediment erosion [49,51]. The effect of near-bed concentration can be explained by the continuity principle. The sediments are entrained into the flow by the turbulent eddies, and when a turbulent eddy picks up a volume of sediment-water mixture from the bed, the same volume of near-bed mixture has to be conveyed back to the bed. If the near-bed concentration is low, the backflow will carry few sediment particles back to the bed. However, if the near-bed concentration is high, the backflow will carry more sediment particles back to the bed. When the near-bed concentration is almost equal to the bed concentration, nearly no net sediment erosion will take place. Following this argument, [51] put forward the reduction factor R = (1 − n 0 − c b )/(1 − n 0 ) to account for the effect of the near-bed concentration c b . Nevertheless, there is no clear definition of the reference line for the near-bed concentration.
To close Equation (17), we propose a new pick-up function φ p, f , which is modified from the existing function of [52]: where D * is a dimensionless particle diameter defined by D * = D 50 3 ∆g/ν, in which ν (m 2 /s) is the kinematic viscosity of water, θ is the Shields parameter, θ cr is the critical Shields parameter for sediment motion and f cr = 1 + sin(α − φ) 2 is an amplification factor for the critical Shields parameter, which can be used when α > φ. Lamb et al. [53] demonstrated that there is a clear dependency between the critical Shields stress for sediment motion and the bed slope; the critical Shields stress increases with bed slope. Therefore, we account for this increase by the empirical factor f cr .
Here, we assume that the near-bed concentration is the average value of the particle concentrations within the inner region, where the velocity gradient is positive. The reason behind this choice is to reduce the dependency of the value of the near-bed concentration on the mesh resolution, as higher resolution results in higher concentration of the first cell above the bed. We also show the influence of using the concentration of the first cell above the bed as c b (instead of the average of the inner region) on the erosional characteristics of the flow (see Section 5.2).
In the current numerical tool, the computations are done using a non-dimensional pick-up function rather than an erosion velocity. Therefore, Equation (17) is recast into the total non-dimensional pick-up function, which reads It is worth noting that [24] did not account for sediment deposition in their formula, as it was assumed that sedimentation is negligible compared to erosion. However, this assumption may be valid for low near-bed concentrations. In breaching, the near-bed concentrations can be very large. In the model, therefore, we include the sedimentation flux (Equation (7)), leading to a reduction of the erosion velocity equal to the sedimentation velocity v s = (c b w s cos α)/(1 − n 0 ). This means that the erosion velocities resulting from the simulations are the net magnitudes.

Model Application
To evaluate the applicability and reliability of the present numerical model, the laboratory experiments carried out by [18] are reproduced numerically and the results are compared. A recapitulation of the experimental set-up of [18] is provided as follows. The series of laboratory experiments were conducted in a breaching tank: 4 m long, 0.22 m wide and 2 m high. The front wall of the tank was made of glass, whereas the back wall was made of steel. A sand deposit of a steep slope (50 • -80 • ) was constructed under water. Owing to the over-steepening of the subaqueous slope, it was essentially unstable and, therefore, it was supported by a confining wall, which should be removed to kick-start the breaching failure and subsequently the turbidity current. The length of the breach face in all experiments did not exceed 1.6 m.
In the breaching tank, a false floor of a mild slope, compared to the breach face, was placed next to the slope toe, where the turbidity current made a turn (see Figure 2). To avoid the reflection of the turbidity current at the end of the downstream region, sand-water mixture was drained constantly during the experiment, while clear water of the same rate was supplied into the tank, so as to guarantee a constant water level. The work of [18] did not include measurements for the turbidity current flowing down the toe of the breach face. Therefore, the present simulations consider the breach face and the current running over it, and do not include a slope transition (see Figure 3). The sand used in the experiments (d 50 = 0.135 mm) was uniformly graded. Therefore, the simulations were run using a uniform particle size fraction of 0.135 mm. Table 1 summarizes the sand properties needed for the numerical computations. Table 1. The properties of the sand used in the experiments [18]. A rectangular numerical flow domain is used, which follows the sloping bed. See Section 4.2 for details on the computational domain and grid. The gravity vector is rotated to account for a correct gravity pull on the density current, and sediment settling takes place along the rotated gravity vector. In the lab experiments, the bed erodes and moves backward with a rate equal to the erosion velocity. In the numerical simulations, there is no bed update and the bed does not move backward, but the erosion velocity (∼mm/s) is prescribed as an inflowing boundary condition at the bed. At the free surface, a rigid lid free slip condition is prescribed.
The flow is internally generated in the computational domain and no inflow or outflow is prescribed at the lateral, left or right end. This will result in a flow reflection at the right wall after some time, but a sufficiently large domain is chosen, and the simulations are stopped before that happens. The width of the domain is equal to the experimental width and closed lateral boundary conditions are applied with a partial slip boundary condition employing a wall roughness k n = 0.2 mm to account for wall resistance of the current.

Model Inputs
Some inputs are needed to run the simulations, such as slope angle, slope length and sand properties. The initial conditions of the numerical runs are summarized in Table 2. Upon the start of the numerical simulations, a discharge of sediment-water mixture equivalent to the corresponding pure breaching is introduced to the numerical domain at the first computational grid cell above the bed. Thereupon, the turbidity current starts developing along the breach face.

Computational Grid
The computational geometry used in the simulations is demonstrated in Figure 4. The domain height is 25 cm, deep enough to avoid effects of the overflow above the turbidity current, while the domain width is 22 cm, equal to the width of the experimental setup. As the purpose of the numerical simulations is to reproduce the current running along the breach face (1.5 m long), it was decided to have a total domain length of 3.5 m. The domain is divided into two zones. The first zone (0 to 1.8 m) corresponds to the breach face over which the turbidity current propagates. The second zone (1.8 m to 3.5 m) functions as a sediment sink, where the sand particles settle out, decelerating the flow and preventing the reflection of the turbidity current upstream. The numerical data after x = 1.5 m are not considered, since they are influenced by the sediment sink.
The computational mesh consists of a total number of about 46 million cells. To reduce the computational time, grid clustering was used in x-direction; a width of 2 mm was used for the cells in the first 1.5 m, while the width of the remaining cells was gradually increased with a growth rate of 1.04 with an upper limit of 5 cm. The cell dimensions in the y and z directions were kept constant with a value of 2 mm and 0.5 mm, respectively (leading to maximum ∆z + = 15 for the first velocity point located at 1 2 ∆z). The average computational time of the runs presented in this study was about 4 days.

Model Validation
As noted earlier, the WALE sub-grid-scale model is used in this study. To ensure that the numerical results are independent of the chosen sub-grid-scale model, an extra simulation has been run using the dynamic Smagorinsky sub-grid-scale model (which is more computationally demanding) and the simulation results were compared with those obtained with the WALE sub-grid-scale model. Indeed, no differences have been observed between the results.
In this section, specific quantitative time-averaged numerical results are compared with the corresponding experimental results to test the validity and reliability of the proposed numerical model. However, some instantaneous flow results are first presented to illustrate the type of flow we are dealing with.

Instantaneous Flow Results
In general, turbidity currents are known to be highly turbulent and breaching-generated turbidity current is not an exception. Figure 5 clearly shows the transition of the flow from laminar to turbulent at about x = 20 cm for the 64 • slope, which is in line with visual observations during the experiment. The top plot demonstrates the 3D high turbulent behavior of the turbidity current while the middle plot illustrates the highly turbulent instantaneous concentration and velocity distribution over the vertical. The latter also shows that the inner region is very thin compared to the total layer thickness and that significant velocities can be found in zones with relative low sediment concentrations between 0.01< c <0.1. In the remainder of this section, time-averaged model results are used in the comparison with the experimental results.

Sediment Erosion
To evaluate the proposed breaching erosion model and the morphodynamic response associated with the turbidity currents, simulation results are compared against the corresponding experimental data (see Figure 6). For each of the numerical runs, the constructed erosion velocities are the average magnitudes over a 3-second time span and a lateral distance of the first 10% of the domain width (2.2 cm). This is because the erosion velocities reported in [18] were computed from the video recordings filmed through the glass wall. As it can be seen, the numerically predicted erosion velocities coincide very well with the experimental data. The prediction accuracy of the erosion model is considered high (mean absolute error = 11%) compared to the acceptable accuracy in the field of sediment transport. The erosion lines in Figure 6 begin with a horizontal segment, where the turbidity current is not yet sufficiently energetic to erode sediment from the breach face. Shortly after that, the turbidity current ignites and increasingly erodes sediments from the breach face.
The experimental data suggests a transition in the erosion rate after a certain point. For example, in the case of the 81 • breach face, the observed erosion velocity is found to be almost constant in the streamwise distance 55-110 cm. It could be that the in situ porosity n 0 was not uniform all over the breach face (see Section 6.6 for the effect of n 0 ). Another hypothesis is that the current was in the bypass mode (the current maintains the sediment load), but that is not captured in the numerical model. This may be attributed to the reference line of the near-bed concentration, which should be defined based on the dimensions of the turbulent eddies transporting the sand grains from the turbidity current back to the breach face. The size of those eddies was not considered in [18] as that requires experimental data of higher resolution.
To elaborate on this hypothesis, we show in Figure 7 two different definitions of the near-bed concentration c b and the corresponding reduction factor R; c b1 is the concentration of the first cell above the bed, while c b2 is the average concentration in the inner region. The value of c b2 tends to become constant downslope and consequently the reduction factor R 2 also remains constant. In contrast, c b1 rapidly increases in the downstream direction and hence the reduction factor R 1 rapidly decreases. This implies that at some point along the breach face, the increase in bed shear stress will be balanced by the increase of the near-bed concentration, transforming the flow from the self-acceleration mode to the bypass mode. In conclusion, the definition of the near-bed concentration influences the computed erosion rates and, consequently, might influence the flow mode.

Flow Spatial Evolution
The characterizing layer thicknesses h and depth-averaged velocities U for Run 1, Run 3 and Run 5 are constructed and compared with corresponding magnitudes derived from the laboratory experiments. These flow characteristics, h and U, are calculated using the following relationships [14,26,54]: where u is locally averaged streamwise flow velocity, z is upward-normal distance from the bed and z ∞ is the height at which the local velocity u is zero. In these runs, the densimetric Froude number is greater than one, indicating that the flow is supercritical, which agrees with the experimental results. The densimetric Froude number is calculated using the following relationship [55]: where C denotes the layer-averaged, volume concentration of sediment defined as: in which c is local concentration of suspended sediment. In addition to the characteristics h and U, the layer-peak velocities u max are also constructed in Figure 8. It can be seen that the model predictions of the spatial evolution of the flow agree qualitatively with the experimental data. However, the model fails to accurately predict the layer thickness (mean absolute error = 46%); it is underestimated. We speculate that this underestimation partly relates to the missing feedback from the sand particles to the flow, leading to less momentum exchange and mixing. This underestimation could be one of the reasons for the deviation between the numerically predicted layer-averaged velocities and the experimental ones (mean absolute error = 18%); the predicted flow velocities for the 64 • and 70 • breach faces are somewhat lower than the measured ones. Another important reason for this difference is that the coupled erosion model was calibrated based on the erosion rates measured at the glass wall of the experimental tank. However, Alhaddad et al. [18] found the erosion to be the highest in the middle of the tank width, where the velocity measurements were obtained, declining towards the side walls. This implies that somewhat more sediment should be entrained from the breach face to the turbidity current to gain higher velocities.

Vertical Structure
The non-dimensional vertical profiles of the streamwise velocities for Run 3 and Run 5 are constructed and compared with the corresponding dimensionless profiles derived from the laboratory experiments. In these numerical runs, the profiles are taken at the same distances from the breach crest as in the experiments. The streamwise local velocities are normalized with the depth-averaged velocity U, while the local vertical distances are normalized with the characterizing layer thickness h.
Overall agreement is found between simulation and experimental results in the vertical structure (see Figure 9), but simulation predictions deviate from experiments in the location of the velocity maximum u max . The latter is numerically predicted to be closer to the bed than measured in the experiments, leading to an over-predicted velocity gradient in the inner region. This could be partially attributed to the underestimation of the layer thicknesses by the model. Another possible reason for the differences is the difficulties and uncertainties in pinpointing the bed position in the laboratory experiments, as stated in the work of [18]. The numerical results demonstrate that the velocity profiles are self-similar, as can be inferred from the experimental results [18].

Vertical Density Distribution
We examine here the capability of the model to capture the internal density distribution of the flow through comparing concentration profiles measured with Conductivity-type Concentration Meter (CCM) (single-point device) along different inclinations versus numerical results. It can be seen from Figure 10

Conclusion on Comparison of Numerical Simulations And Experiments
In view of the presented systematic comparison between the numerical and experimental results, it can be concluded that the numerical model gives fairly reasonable predictions of the flow characteristics and the associated morphodynamic response. In the next section, therefore, we investigate further flow characteristics, which were not possible to analyze through the experimental data.

Further Analysis of Numerical Results
In Section 5 we showed that the model was good at simulating the experimental observations. This gives us confidence to closely investigate the flow and turbulence structure, to determine some characterizing parameters and to analyze the sensitivity of the numerical results to some initial conditions. This will be the scope of this section. Figure 11 depicts the spatial development of the layer-averaged concentration C (Equation (23)) for three slope angles. Clearly, steeper slopes result in a higher C owing to the higher sediment erosion. The results show that C increases in the downstream direction, in the same manner as the layer-averaged velocity U (see Figure 8). This supports the conclusion drawn by [18] that breaching-generated turbidity currents are self-accelerating. According to [56], self-accelerating flow is characterized by the downstream increase of flow velocity, which is caused by downstream increase in suspended sediment concentration.

Spatial Evolution of Vertical Density Distribution
The vertical profiles of the sediment concentrations for Run 3 and Run 6 are depicted in Figure 12. As the flow further travels downslope, the near-bed concentration increases. Moreover, steeper slopes result in higher near-bed concentrations, which can be attributed to the higher erosion rates caused by a larger gravity force and more erosive turbidity currents.

Reynolds Stresses
The Reynolds stress distribution corresponds to the velocity gradient within the flow body as maximum stresses occur where the gradient is largest. The normal Reynolds stresses can be obtained from the turbulent fluctuations of downstream u and bed-normal velocity w as follows: The Reynolds stresses are calculated using the turbulent velocity components averaged over a 3-second time span. For three different slopes at x = 1 m, Reynolds stresses are normalized by ρU 2 and their distribution is shown in Figure 13a. The Reynolds stress increases significantly near the breach face, reaching the largest positive Reynolds stress at z/h ∼ 0.025, where the bottom boundary layer ends. Around and below the velocity maximum, z/h = 0.045-0.085, Reynolds stresses are close to zero. Further upwards, in the outer region of the flow, where the velocity gradient is negative, Reynolds stresses become negative, reaching the largest negative Reynolds stress at z/h ∼ 0.45. This elevation has the largest negative velocity gradient. Above this elevation, Reynolds stresses decreases towards nearly zero at the upper boundary of the flow. It is found that Reynolds stresses take a self-similar shape.     Owing to the acceleration of the flow downslope, the flow becomes more turbulent, leading to higher Reynolds stresses in the downstream direction. These stresses, as can be seen in Figure 13c, are higher for steeper breach faces, as these result in higher flow velocities.

Turbulent Kinetic Energy
Profiles of turbulent kinetic energy (TKE) are constructed to analyze the vertical distribution of the TKE within the flow. The total TKE is calculated using the turbulent velocity components averaged over a 3-second time span as follows where u is the streamwise component, v is the across-stream component and w is bed-normal component. Numerical TKE normalized by ρU 2 is plotted in Figure 14a for three different slope angles at x = 1 m. The TKE profiles show a significant increase of TKE in the bottom boundary layer z/h ∼ 0-0.025 and then a decrease until z/h ∼ 0.045, a little below the velocity maximum. Above the elevation of the velocity maximum, the TKE increases again, reaching the largest TKE at z/h ∼ 0.45, coinciding with the largest negative velocity gradient and largest Reynolds stress. Above this elevation, TKE tends to decrease towards the upper boundary of the flow. The results suggest that TKE takes a self-similar shape. Figure 14b illustrates the spatial development of TKE along 64 • and 77 • breach faces. In a similar way to the Reynolds stresses, the TKE increases downslope and is higher for steeper slopes.

Bed Shear Stress and Bed Friction Coefficient
The bed shear stress τ b determines the erosive power of the flow and can be expressed in terms of a bed friction coefficient C f , which relates the bed shear velocity u * to the layer-averaged velocity U as follows This relation is usually used in depth-averaged models. These are more computationally efficient and can be used to make some preliminary computations. In this study, the value of u * is calculated using Equation (11). The average value of the C f is calculated here from the numerical results.
The calculated values of C f show that the bed friction coefficient is not a constant parameter ( Figure 15); it decreases downslope. This is because the thickness of the current increases along the streamwise direction, resulting in an increased bulk Reynolds number and, consequently, a decreased drag coefficient. The results also show that steeper slopes lead to a lower bed friction coefficient. For the considered slope angles and traveling distance, C f ranges from 0.028 to 0.006 with an average value of 0.011. This is consistent with the range of values reported in the literature (e.g. [54]

Influence of In Situ Porosity
The experiments of [18] were all conducted using a constant sand in situ porosity n 0 = 0.4. To investigate the effect of n 0 on the flow and erosion velocity, an additional numerical simulation for 67 • breach face was run using n 0 = 0.44, which corresponds to a sand relative density of 68%. A comparison between the numerical results for n 0 = 0.4 (Run 4) and n 0 = 0.44 (Run 8) is shown in Figure 16. It can be seen from Figure 16a that higher n 0 results in higher pure breaching and hence higher erosion velocities downslope; the average increase in erosion, in the considered case, is significant (47%). Higher n 0 leads to a lower hydraulic gradient (which acts as a stabilizing force) during shearing, increasing the erosion velocity. As a result, the flow becomes denser and runs faster downslope, as depicted in Figure 16b. The difference of the layer-averaged velocities between the two runs is magnified downstream.

Conclusions
We use a 3D numerical model to conduct large eddy simulations of turbidity currents generated during breaching flow slides. The qualitative and quantitative comparison between experimental and simulation results indicates that the proposed numerical tool can reasonably reproduce several substantial aspects of the flow, such as vertical density distribution, and spatial development down the breach face. A limitation of the model is that it underestimates the thickness of the current. Considering the special conditions of breaching, a new breaching erosion model is proposed and verified through published experimental data. Good agreement is observed between the experimental and numerically predicted erosion rates. The numerical results confirm the self-accelerating behavior of breaching-generated turbidity currents. Both Reynolds stresses and TKE decrease sharply below their maximum in the bottom boundary layer. In addition, the numerical results reveal that breaching-generated turbidity currents exhibit a self-similar behavior; velocity, concentration, Reynolds stress, and TKE profiles take a self-similar shape. Based on a sensitivity analysis, sand erosion during breaching is found to be susceptible to the in situ porosity; the lower the in situ porosity, the higher the sand resistance to erosion. If the present numerical model is extended to account for the morphological changes of the breach face, it would serve as a good tool to predict the breaching evolution and stability.