1. 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.
3. 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.
3.1. 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
[
] is the erosion rate of pure breaching,
[
] is the internal friction angle,
[-] is the in situ porosity of the sand,
[
] is the sand permeability at the loose state,
[
] is the density of water, and
is the relative change in porosity, in which
[-] is the maximum porosity of the sand.
Pure breaching is particularly sensitive to the magnitude of sand permeability
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
, which leads to the expression of the corrected pure breaching
:
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.
3.2. 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.
3.3. 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
is the total erosion velocity including pure breaching and flow-induced erosion,
is Shields velocity for sand grains, and
is an empirical non-dimensional pick-up function, which does not account for the bed dilative behavior nor the sloping bed:
where
[
)] is the sediment pick-up rate perpendicular to the bed surface in which
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
and
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
when another value of
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
to account for the effect of the near-bed concentration
. 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
, which is modified from the existing function of [
52]:
where
is a dimensionless particle diameter defined by
, in which
(m
/s) is the kinematic viscosity of water,
is the Shields parameter,
is the critical Shields parameter for sediment motion and
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
.
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
(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
. This means that the erosion velocities resulting from the simulations are the net magnitudes.
4. 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 (
= 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.
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 mm to account for wall resistance of the current.
4.1. 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.
4.2. 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 for the first velocity point located at ). The average computational time of the runs presented in this study was about 4 days.
5. 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.
5.1. 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
0.1. In the remainder of this section, time-averaged model results are used in the comparison with the experimental results.
5.2. 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
was not uniform all over the breach face (see
Section 6.6 for the effect of
). 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
and the corresponding reduction factor
R;
is the concentration of the first cell above the bed, while
is the average concentration in the inner region. The value of
tends to become constant downslope and consequently the reduction factor
also remains constant. In contrast,
rapidly increases in the downstream direction and hence the reduction factor
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.
5.3. 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
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
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.
5.4. 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
. 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].
5.5. 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 that the time-averaged concentration profiles predicted by the model fall within the scatter range of the corresponding profiles resulted from the laboratory experiments. Also, the very high near-bed concentrations in the order of
c = 0.4–0.52 are captured in the numerical model. This indicates that the numerical model can adequately predict the vertical density distribution of the considered turbidity current.
5.6. 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.