1. Introduction: submarine fans and long-runout turbidity currents
Submarine, or deep-sea fans represent the major ultimate sinks for terrestrial sediment. Sediment is transported across these fans through submarine channels that may extend for hundreds or thousands of kilometres into water that is up to thousands of metres deep. The longest of these fans is the Bengal Fan at
$\sim3000$
km (Curray et al. Reference Curray, Emmel and Moore2002; Schwenk et al. Reference Schwenk, Spieß, Hübscher and Breitzke2003). Other very large fans include the Congo Fan at
$\sim1000$
km (Picot et al. Reference Picot, Droz, Marsset, Dennielou and Bez2016), the Indus Fan at
$\sim1500$
km, the Amazon Fan at
$\sim700$
km, the Mississippi Fan at
$\sim500$
km, and the Rhone Fan at
$\sim400$
km (Wetzel Reference Wetzel1993; Deptuck & Sylvester Reference Deptuck and Sylvester2017). Submarine fan systems are emplaced largely by sediment transported through submarine channels, which both convey and are constructed by turbidity currents, i.e. dense bottom flows that obtain their driving power from suspended sediment (Daly Reference Daly1936; Kuenen Reference Kuenen1938; Pirmez & Imran Reference Pirmez and Imran2003). The channels are commonly highly sinuous (Imran et al. Reference Imran, Parker and Pirmez1999; Schwenk et al. Reference Schwenk, Spieß, Hübscher and Breitzke2003), and build the fan itself by avulsing across the fan surface (Jobe et al. Reference Jobe, Howes, Straub, Cai, Deng, Laugier, Pettinga and Shumaker2020).
Views of the Amazon Submarine Fan and the meandering Amazon Channel itself are given in figure 1. The sinuosity of this channel, and many other channels on submarine fans, implies that the down-channel lengths may notably exceed the length of the fan itself.
Covault et al. (Reference Covault, Fildani, Romans and McHargue2011) have documented 20 long profiles of submarine channels as they traverse an upstream canyon and emanate onto the fan below (figure 2 a). Eight of these are 200 km or more in length. It can be seen from the figure that channels on canyon–fan systems have long profiles that vary from approximately constant slopes to strongly upward concave.
Figure 2(b) shows the along-channel profile of the Amazon Channel of figure 1, including channel thalweg, levee crest and canyon top. Distances are measured down-channel and thus include the effect of sinuosity. The reach in the submarine canyon is 120 km long, and the downstream reach on the submarine fan is 760 km long. Bed slope ranges from approximately 0.014 upstream to approximately 0.002 downstream.
The turbidity currents that excavate submarine canyons and emplace submarine fans thus must also run out as much as hundreds or thousands of kilometres without dissipating or becoming so thick and dilute that they cannot coherently channelize themselves. We refer to such currents as long-runout turbidity currents. Although numerous models of turbidity currents have been presented to date, none has had the capability of satisfying this constraint over such lengths. Here, we provide a resolution to this problem.
Amazon Submarine Fan and channels: (a) the fan itself (after Mikkelsen et al. Reference Mikkelsen, Maslin, Giraudeau and Showers1997); (b) a
$\sim200$
km long reach of the submarine channel (the red box in figure 1
a) (from IFREMER, France).
(a) Down-channel long profiles of 20 canyon–fan systems (after Covault et al. Reference Covault, Fildani, Romans and McHargue2011). (b) Long profiles of channel thalweg, levee crest and top of canyon for the Amazon Channel of figure 1. The channel is confined within the Amazon Canyon for the first 120 km, then extends out 760 km on the fan. The distances are measured along the channel thalweg (based on Pirmez & Imran Reference Pirmez and Imran2003).
2. Existing models of turbidity currents
Meiburg & Kneller (Reference Meiburg and Kneller2010) presented an overview of both models of turbidity current dynamics and their objectives. Models to date that predict either spatial or spatiotemporal evolution of such currents fall into three classes. The first of these includes layer-averaged models; the second encompasses Reynolds-averaged models that can resolve the structure of the flow in the upward normal direction and averaged flow fields; and the third encompasses high-fidelity – large eddy simulations (LES) or direct numerical simulations (DNS) – models that resolve all or part of the turbulent structure.
One-dimensional (1-D) layer-averaged models for turbidity currents were presented by Fukushima et al. (Reference Fukushima, Parker and Pantin1985) and Parker et al. (Reference Parker, Fukushima and Pantin1986) in the context of submarine canyons. These models are based on an extension of the 1-D model of Ellison & Turner (Reference Ellison and Turner1959) for the downstream evolution of sediment-free dense underflows, such as those driven by thermohaline effects. Ellison & Turner (Reference Ellison and Turner1959) assumed their ambient fluid to be infinitely deep. In turn, the 1-D assumption implicitly assumes infinitely high vertical frictionless sidewalls. These assumptions were also used in the formulation of the 3-equation and 4-equation models, and is retained in the analysis below. The models of Fukushima et al. (Reference Fukushima, Parker and Pantin1985) and Parker et al. (Reference Parker, Fukushima and Pantin1986) adapt concepts from Pantin (Reference Pantin1979) and Parker (Reference Parker1982) to explain how turbidity currents could ‘ignite’ or self-accelerate via the entrainment of bed sediment. Further developments in layer-averaged modelling have been presented by Garcia (Reference Garcia1994), Bonnecaze et al. (Reference Bonnecaze, Huppert and Lister1993), Fay (Reference Fay2012), Hu et al. (Reference Hu, Cao, Pender and Tan2012, Reference Hu, Pähtz and He2015), Cao et al. (Reference Cao, Li, Pender and Liu2015) and Bolla Pittaluga et al. (Reference Bolla Pittaluga, Frascati and Falivene2018).
The 3-equation and 4-equation layer-averaged models of Fukushima et al. (Reference Fukushima, Parker and Pantin1985) and Parker et al. (Reference Parker, Fukushima and Pantin1986) have been used to study the formation of sediment waves and submarine gullies on the seafloor (Izumi Reference Izumi2004), and cyclic step instability within the flow (Kostic & Parker Reference Kostic and Parker2006; Wu & Izumi Reference Wu and Izumi2022). A two-dimensional (2-D) extension of variants of these models has been used to explain incipient self-channelization of turbidity currents via levee emplacement (Imran et al. Reference Imran, Parker and Katopodes1998; Halsey & Kumar Reference Halsey and Kumar2019). Wahab et al. (Reference Wahab, Hoyal, Shringarpure and Straub2022) have applied the 4-equation model to the morphodynamics of submarine fans. Traer et al. (Reference Traer, Fildani, Fringer, McHargue and Hilley2018a ,Reference Traer, Fildani, Fringer, McHargue and Hilley b ) have used a version of the 4-equation model to study flow stripping over levees. A version of the model was further developed to simulate the excavation of submarine canyons (Zhang et al. Reference Zhang, Parker, Izumi, Cartigny, Li and Wang2017).
Reynolds-averaged models of turbidity current dynamics that can resolve the vertical structure of the flow, and in particular
$k{-}\epsilon$
models, have been presented by Eidsvik & Brørs (Reference Eidsvik and Brørs1989), Sequeiros et al. (Reference Sequeiros, Cantelli, Viparelli, White, Garcia and Parker2009), Yeh et al. (Reference Yeh, Cantero, Cantelli, Pirmez and Parker2013), Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018) and Iwasaki & Parker (Reference Iwasaki and Parker2020). High-fidelity models have been presented by Cantero et al. (Reference Cantero, Balachandar, Cantelli, Pirmez and Parker2009a
,Reference Cantero, Balachandar and Parker
b
), Biegert et al. (Reference Biegert, Vowinckel, Ouillon and Meiburg2017), Salinas et al. (Reference Salinas, Cantero, Shringarpure and Balachandar2019a
,Reference Salinas, Cantero, Shringarpure and Balachandar
b
, Reference Salinas, Balachandar, Shringarpure, Fedele, Hoyal and Cantero2020, Reference Salinas, Balachandar, Shringarpure, Fedele, Hoyal, Zuñiga and Cantero2021a
, Reference Salinas, Zúñiga, Cantero, Shringarpure, Fedele, Hoyal and Balachandar2022) and Xie et al. (Reference Xie, Zhu, Hu, Yu and Pan2023b
). Wells & Dorrell (Reference Wells and Dorrell2021) provide a comprehensive overview of such techniques, and the features that they capture in addition to those captured by layer-averaged approaches.
These different modelling approaches, each of which has its intrinsic value, cannot be used generally to directly predict long-runout turbidity currents over hundreds or thousands of kilometres. This is due to either computational limitations (e.g. DNS) or the configuration studied (e.g. lock exchange). The layer-averaged model of Bolla Pittaluga et al. (Reference Bolla Pittaluga, Frascati and Falivene2018) adopts the mechanism of the flow stripping over pre-existing levees as one element to overcome the overthickening problem, and so capture the long-runout. Such an assumption, however, precludes the possibility that the same model predicts the levee formation process (such as Imran et al. Reference Imran, Parker and Katopodes1998; Halsey & Kumar Reference Halsey and Kumar2019).
The Reynolds-averaged model of Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018) and the high-fidelity model of Salinas et al. (Reference Salinas, Balachandar, Shringarpure, Fedele, Hoyal, Zuñiga and Cantero2021b ) applied to the Froude-subcritical regime have demonstrated that flow conditions exist that would facilitate long-runout turbidity currents, the former due to flow detrainment, and the latter through suppression of turbulence. However, none of these Reynolds-averaged, DNS and LES models can predict the evolution of the current over hundreds or thousands of kilometres, as such simulations require large computational domains that are computationally prohibitive. Traditional layer-averaged models based on 3-equation and 4-equation models do not present such computational difficulty. They nevertheless suffer from a deficiency in the formulation itself, as illustrated below.
Definition diagrams for layer-averaged turbidity currents: (a) single-layer formulation such as used in the 3-equation model; (b) two-layer formulation proposed here.
2.1. Deficiency of layer-averaged approaches to turbidity currents
The deficiency in question is common to both the 3-equation and 4-equation models of Parker et al. (Reference Parker, Fukushima and Pantin1986), so only the 3-equation model is outlined here. The configuration is shown in figure 3(a). The turbidity current is contained within a single layer. It runs over a bed with slope
$S$
, and has thickness
$\delta$
and layer-averaged stream velocity
$U$
. It carries a dilute suspension of sediment with layer-averaged volumetric concentration
$C$
(
$ \ll 1$
). The sediment has submerged specific gravity
$R$
(where
$R = 1.65$
for quartz in water). With
$t$
and
$x$
denoting time and the down-channel coordinate, and
$g$
denoting gravitational acceleration, the governing equations for momentum, fluid mass and suspended sediment conservation are
\begin{equation} \frac {{\partial \delta U}}{{\partial t}} + \frac {{\partial {\delta }U_{}^2}}{{\partial x}} = - \frac {1}{2}Rg\frac {{\partial C\delta ^2 }}{{\partial x}} + RgC\delta S - {C_{fb}}U_{}^2, \end{equation}
\begin{equation} \frac {{\partial \delta }}{{\partial t}} + \frac {{\partial U\delta }}{{\partial x}} = {e_{ws}}U, \end{equation}
\begin{equation} \frac {{\partial C\delta }}{{\partial t}} + \frac {{\partial UC\delta }}{{\partial x}} = {v_s}\left ({E_s} - rC\right ), \end{equation}
where
$C_{fb}$
is a dimensionless coefficient of bed friction, here taken as constant for simplicity,
$e_{ws}$
is a coefficient of water entrainment across the interface between the turbidity current and the ambient fluid,
$v_s$
is the fall velocity of sediment,
$E_s$
is a dimensionless coefficient of sediment entrainment from the bed into suspension, which is in turn a function of near-bed flow, and
$r$
is the ratio of near-bed concentration to layer-averaged concentration. The closure relation for
$e_{ws}$
has been obtained empirically as (Parker et al. Reference Parker, Garcia, Fukushima and Yu1987)
\begin{equation} {e_{ws}}={e_{ws}}[Ri_b]= \frac {{0.075}}{{\sqrt {1 + 718\,Ri_b^{2.4}} }}, \end{equation}
\begin{equation} R{i_b} = \frac {{RgC\delta }}{{{U^2}}} = \frac {{Rg{q_s}}}{{{U^3}}}, \end{equation}
\begin{equation} {q_s} = U\delta C. \end{equation}
Here,
$R_{ib}$
is a bulk Richardson number, and
$q_s$
is the volume transport rate of suspended sediment per unit width
$[L^2T^{-1}]$
. In the case of steady flows that develop spatially downstream, (2.1), (2.2) and (2.3) can be cast in the forms
\begin{equation} \frac {\delta }{U}\frac {{{\rm d}U}}{{{\rm d}x}} = \frac {{ - \left ( {1 + \dfrac {1}{2}R{i_b}} \right ){e_{ws}} + R{i_b}\,S - {C_{fb}} - \dfrac {1}{2}\dfrac {{{v_s}}}{U}r\,R{i_b}\left ( {\dfrac {{{q_{se}}}}{{{q_s}}} - 1} \right )}}{{ {1 - R{i_b}} }}, \end{equation}
\begin{equation} \frac {{{\rm d}\delta }}{{{\rm d}x}} = \frac {{\left ( {2 - \dfrac {1}{2}R{i_b}} \right ){e_{ws}} - R{i_b}\,S + {C_{fb}} + \dfrac {1}{2}\dfrac {{{v_s}}}{U}r\,R{i_b}\left ( {\dfrac {{{q_{se}}}}{{{q_s}}} - 1} \right )}}{{ {1 - R{i_b}} }}, \end{equation}
\begin{equation} \frac {\delta }{{{q_s}}}\frac {{{\rm d}{q_s}}}{{{\rm d}x}} = \frac {{{v_s}}}{U}r\,\left ( {\frac {{{q_{se}}}}{{{q_s}}} - 1} \right ), \end{equation}
where
$q_{se}$
is the value of
$q_s$
that would be in equilibrium with the local flow:
\begin{equation} {q_{se}} = U\delta \frac {{{E_s}}}{r}. \end{equation}
The densimetric Froude number of the flow
$Fr_{db}$
can be defined as
\begin{equation} F{r_{db}} = \frac {U}{{\sqrt {RgC\delta } }} = \frac {{{U^{3/2}}}}{{\sqrt {Rg{q_s}} }} = Ri_b^{ - 1/2}. \end{equation}
In the case of Froude-supercritical flow (
$Fr_{db}\gt 1,\ Ri_b\lt 1$
), (2.7), (2.8) and (2.9) can be integrated downstream upon specification of upstream values
$U$
,
$\delta$
and
$q_s$
. As noted above, these relations can be used to predict self-acceleration upon the assumption of appropriate functional forms for
$E_s$
and
$r$
. The major deficiency of this model is, however, best seen in the case of bypass flow, according to which there is no net exchange of sediment with the bed:
\begin{equation} {q_s} = \text{constant} \quad (\text{e.g. } {q_{se}}). \end{equation}
Such flows can be realized, for example, by running the currents over a sediment-starved bed.
The above equations possess a normal flow solution over a constant bed slope
$S$
for bypass flow. Here, the normal flow state is defined in the sense of Ellison & Turner (Reference Ellison and Turner1959), and in such a state, Richardson number
$Ri_b$
and velocity
$U$
attain constant values, and thickness
$\delta$
increases linearly with distance downstream:
\begin{equation} - \left ( {1 + \tfrac {1}{2}R{i_b}} \right ){e_{ws}}[ {R{i_b}}] + R{i_b}\,S - {C_{fb}} = 0, \end{equation}
\begin{equation} \frac {{{\rm d}\delta }}{{{\rm d}x}} = {e_{ws}} [ {R{i_b}} ]. \end{equation}
For given values of
$S$
and
$C_{fb}$
, (2.13) can be solved in conjunction with (2.4) to obtain the normal flow Richardson number
$Ri_{bn}$
, from which normal velocity
$U_n$
is found to be
\begin{equation} {U_n} = {\left ( {\frac {{Rg{q_s}}}{{R{i_{bn}}}}} \right )^{1/3}}. \end{equation}
The defect associated with these models becomes apparent upon consideration of (2.14). For example, we consider a value
$Ri_{bn}=0.7$
. This corresponds to a Froude-supercritical flow in the sense that the densimetric Froude number
$Fr_{db}$
takes the value
$1.20 \gt 1$
. According to (2.4),
$e_{ws}$
takes the constant value 0.0043. A current that is 5 m thick upstream (
$x = 0$
) and has the normal velocity at that point would attain a thickness of at least 1720 m at
$x = 400$
km. The suspended sediment concentration at
$x = 400$
km would be an order of
$10^{-3}$
times smaller than its upstream value. According to Jobe et al. (Reference Jobe, Howes, Straub, Cai, Deng, Laugier, Pettinga and Shumaker2020), channels on the Amazon Submarine Fan have bankfull depths ranging from 147 m to 10 m down-fan. As noted above in the context of figure 2, the channel is at least 760 km long. Similarly, the channel of the Congo Fan has depth approximately 100–150 m for the first 900 km of the channel (Hasenhündl et al. Reference Hasenhündl2024). Referring to the example with
$Ri_b=0.7$
, there is no obvious way for 1720 m thick turbidity current, with a suspended sediment concentration that is of the order of one-thousandth of its upstream value, to follow a channel that is 10–150 m deep, much less construct it.
The models like Bolla Pittaluga et al. (Reference Bolla Pittaluga, Frascati and Falivene2018) suffer from the same defect of overthickening. They achieve long-runout only by limiting current thickness by means of overflow across pre-existing levees. Were the levees not already confining the flow, the flow would overthicken, and the suspended sediment concentration would become dilute to the point where the flow would be incapable of constructing them.
Cao et al. (Reference Cao, Li, Pender and Liu2015) have succeeded in running a layer-averaged model of turbidity currents in a reservoir over 60 km. Their innovative model is able to capture the plunge point where the river dives into the reservoir to form a bottom turbidity current. During turbidity current events they studied, the deepest part of the reservoir is approximately 60 m, or approximately three times the thickness of the turbidity currents. Due to the shallow environment, Cao et al. (Reference Cao, Li, Pender and Liu2015) added a dynamic formulation of the flow in the ambient water above the current as well as the current itself, calling their formulation a ‘double layer-averaged model’. The volume sediment concentration in the current is approximately 0.085, corresponding to a hyperconcentrated flow. This and the relatively slow-moving flow dictate the value of
$Ri_b$
of the order of hundreds, in which case values of
$e_{ws}$
are so small that thickening over 60 km is negligible. This result, however, cannot be used to formulate a general model of long-runout turbidity currents because such large bulk Richardson numbers with near-vanishing entrainment of ambient water constitute a special case.
The Cao et al. (Reference Cao, Li, Pender and Liu2015) model includes a turbidity current layer and an ambient water layer, the dynamics of which must be considered in the shallow setting of the reservoir that they model. In that sense, our model is a ‘three-layer model’: driving layer, driven layer and ambient layer. In our model, we take the ambient water to be infinitely deep, so can treat it as stagnant, which enables the description of the deep sea environment.
Here, we seek a model that allows long-runout turbidity currents over hundreds or thousands of kilometres for arbitrary values of the Richardson number, which neither overthicken nor become overdilute, and as such would be competent to emplace their own levees far downstream (see Imran et al. Reference Imran, Parker and Katopodes1998; Halsey & Kumar Reference Halsey and Kumar2019).
2.2. Water detrainment from turbidity currents
An examination of the above calculations reveals that in the case of bypass conditions, when sediment fall velocity
$v_s$
is neglected in the problem, the formulation becomes identical to that of Ellison & Turner (Reference Ellison and Turner1959) for a conservative contaminant such as dissolved salt. Yet fall velocity should play a role even in bypass suspensions. An easy way to see this is in terms of the Rouse solution for the equilibrium (bypass) vertical distribution of suspended sediment in an open channel. Even though there is no net bed erosion or deposition, the higher the fall velocity, the more the suspended sediment profile is biased towards the bed. Such grain size bias is also observed in direct measurements of turbidity currents and their deposits in the Monterey Canyon (Symons et al. Reference Symons, Sumner, Paull, Cartigny, Xu, Maier, Lorenson and Talling2017), and is built into, for example, the high-fidelity calculations reported in Balachandar et al. (Reference Balachandar, Salinas, Zuniga, Cantero and Kneller2024). Bonnecaze et al. (Reference Bonnecaze, Huppert and Lister1993) made a first attempt to consider the role of fall velocity in the case of a 1-D formulation of non-turbulent lock-exchange flow driven by suspended sediment.
Further insight into the role of sediment fall velocity can be gained by the study of turbidity currents entering into bowl-like basins with horizontal scales of the order of tens of kilometres, called minibasins. These basins can fill over time due to the delivery of sediment from turbidity currents (Lamb et al. Reference Lamb, Hickson, Marr, Sheets, Paola and Parker2004). In the course of experiments on turbidity currents flowing into minibasins, Lamb et al. (Reference Lamb, Toniolo and Parker2006) and Toniolo et al. (Reference Toniolo, Lamb and Parker2006a
,Reference Toniolo, Parker, Voller and Beaubouef
b
) recognized a phenomenon that they called detrainment. Under the right circumstances, a fully turbulent turbidity current can flow continuously into a minibasin, yet no sediment escapes over the downstream lip of the minibasin. In other words, the turbidity current can form a relatively stagnant pond with a settling interface that equilibrates below the downstream lip of the minibasin. If the pond is sufficiently stagnant, so that there is negligible flow circulation within it (Reece et al. Reference Reece, Dorrell and Straub2024), then water detrains across the interface and then escapes the minibasin at the rate
$v_sA$
, where
$A$
is the surface area of the settling interface within the minibasin.
With the above in mind, Toniolo et al. (Reference Toniolo, Parker, Voller and Beaubouef2006b ) proposed a formulation according to which (2.2) is amended to
\begin{equation} \frac {{\partial \delta }}{{\partial t}} + \frac {{\partial U\delta }}{{\partial x}} = {e_{ws}}U - {v_s}. \end{equation}
In simple terms the fall velocity in (2.16) indicates that the sediment ‘fights back’ against turbulent entrainment into the ambient fluid above. Bolla Pittaluga et al. (Reference Bolla Pittaluga, Frascati and Falivene2018) incorporated this formulation into the 3-equation model.
Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018) further developed this idea using a
$k{-}\epsilon$
model of turbidity currents that naturally accounts for the effect of water detrainment mediated by sediment fall velocity through its presence in the equation of conservation of suspended sediment. Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018) modelled the evolution of the flow down a slope that is uniform in space but developing in time. After a sufficient amount of time, the flow segregates into two layers. The turbulent flow in the bottom layer contains nearly all the suspended sediment, and eventually achieves a near steady-state thickness and streamwise velocity profile. The flow in the top layer is also turbulent but nearly sediment-free, and thickens monotonically in time. They referred to the bottom layer as the ‘driving’ layer, in that the suspended sediment sequestered there provides the impelling force for the flow. They referred to the top layer as the ‘driven’ layer, in that the nearly sediment-free water there is more or less simply dragged along by the driving layer, as in the case of a flow above a plate moving at constant velocity.
The interface between the driving layer and the driven layer in the model of Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018) corresponds to a settling interface. This interface is not necessarily as sharp as that seen in a fully ponded minibasin, because the downward tendency of the settling interface works against turbulent mixing of the flow itself. This notwithstanding, as the driven layer thickens, the mean concentration of suspended sediment in it becomes negligible compared to the driving layer.
3. A two-layer formulation
The
$k{-}\epsilon$
model of Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018) is not easily implemented on a scale of hundreds or thousands of kilometres. The essence of the model results can, however, be cast in terms of a much simpler two-layer, layer-averaged model that does have that capability. This configuration is summarized in figure 3(b). Our model is again 1-D, with the implicit assumtion of infinitely high, vertical frictionless sidewalls. As such, it does not directly pertain to cases including meandering or levee overflow. Guidance for future work in this regard is provided in § 4. We apply the model to bed slopes ranging from a high of 0.03 to a low of 0.003.
The fluid mechanical basis for the two-layer model presented here is the two-layer formulation of Arita & Jirka (Reference Arita and Jirka1987 Reference Arita and Jirka a ,Reference Arita and Jirka b ), originally designed for the treatment of saline wedges. (Several misprints were corrected in Arita Reference Arita1998.) That framework is adapted here, but the characterization of the boundary between the two layers is amended in terms of a settling interface. We also replace the relation for water entrainment used by Arita & Jirka (Reference Arita and Jirka1987 Reference Arita and Jirka b ) for saline wedges to (2.4) (Parker et al. Reference Parker, Garcia, Fukushima and Yu1987), which is more appropriate for density underflows.
Let
$\delta _L$
and
$\delta _U$
denote the thicknesses of the lower (driving) layer and upper (driven) layer in figure 3(b). The corresponding layer-averaged velocities are
$U_L$
and
$U_U$
. The layer-averaged volume suspended sediment concentration in the lower layer is
$C$
; the upper layer is approximated as sediment-free. The friction coefficient at the interface between the two layers is denoted as
$C_{fi}$
. The coefficient of water entrainment across the interface between the lower and upper layer is denoted as
$e_{ws}$
, whereas the corresponding coefficient between the upper layer and the ambient water is denoted as
$e_{wo}$
. Insofar as the upper layer is (to a first approximation) sediment-free, the value of
$e_{wo}$
can be computed from (2.4) as the limiting value in the absence of stratification (
$Ri_b\to 0$
), so that
$e_{wo}\to 0.075$
. The coefficient of bed friction
$C_{fb}$
corresponds to the turbulent rough flow considered here.
The governing equations for the lower layer are
\begin{equation} \frac {{\partial {U_L}{\delta _L}}}{{\partial t}} + \frac {{\partial U_L^2{\delta _L}}}{{\partial x}} = - \frac {1}{2}Rg\frac {{\partial C\delta _L^2}}{{\partial x}} + RgC{\delta _L}S - {C_{fb}}U_L^2 - {C_{fi}}\,{\lvert {{U_L} - {U_U}} \rvert \left ( {{U_L} - {U_U}} \right )}, \end{equation}
\begin{equation} \frac {{\partial {\delta _L}}}{{\partial t}} + \frac {{\partial {U_L}{\delta _L}}}{{\partial x}} = {e_{ws}}\left ( {{U_L} - {U_U}} \right ) - {v_s}, \end{equation}
\begin{equation} \frac {{\partial C{\delta _L}}}{{\partial t}} + \frac {{\partial {U_L}C{\delta _L}}}{{\partial x}} = {v_s}\left ({E_s} - rC\right ). \end{equation}
A derivation of (3.2) from the 2-D (streamwise – upward normal) continuity equation is provided in Appendix A. A corresponding derivation of (3.3) is given in Appendix B. In the case of bypass flows, (3.3) is replaced by (2.12). As opposed to the 3-equation model of Parker et al. (Reference Parker, Fukushima and Pantin1986), however, the effect of sediment does not vanish in the bypass case; it enters through the right-hand side of (3.2). The corresponding forms for the upper layer are
\begin{equation} \frac {{\partial {U_U}{\delta _U}}}{{\partial t}} + \frac {{\partial U_U^2{\delta _U}}}{{\partial x}} = {C_{fi}}\,\lvert {{U_L} - {U_U}}\rvert\, {\left ( {U_L^{} - U_U^{}} \right )}, \end{equation}
\begin{equation} \frac {{\partial {\delta _U}}}{{\partial t}} + \frac {{\partial {U_U}{\delta _U}}}{{\partial x}} = {e_{wo}}{U_U} - {e_{ws}}\left ( {{U_L} - {U_U}} \right ) + {v_s}. \end{equation}
Arita & Jirka (Reference Arita and Jirka1987b
) evaluate the interfacial friction coefficient as
${C_{fi}} = 2{e_{ws}}$
. The Richardson number used in (2.4) must be modified to account for the two-layer structure of the flow. Insofar as it represents a ratio between bulk gravitational (buoyancy) gradient and bulk shear strength, the appropriate generalization is
$Rg\delta\, \varDelta C/ (\varDelta U )^2$
. Insofar as we approximate the upper layer as sediment-free, this reduces to
$RgC\delta _L/ (U_L - U_U )^2$
. We thus amend the relation for the turbulent entrainment coefficient
$e_{ws}$
to
\begin{align} {e_{ws}} & = \frac {{0.075}}{{\sqrt {1 + 718\,Ri_I^{2.4}} }},\\[-10pt]\nonumber \end{align}
\begin{align} R{i_I} & = \frac {{RgC\delta _L}}{{{{\left ( {{U_L} - {U_U}} \right )}^2}}} = \frac {{Rg{q_s}}}{{{U_L}{{\left ( {{U_L} - {U_U}} \right )}^2}}}, \end{align}
which is used to obtain numerical results of the two-layer model.
The assumption of a sediment-free upper layer merits some elaboration. The present model specifically tracks the boundary between the two layers in terms of a settling interface. Thus as sediment is mixed upwards against fall velocity, the interface itself moves upwards, rather than suspended sediment being transferred into the upper layer. In a turbulent flow, the interface can always be expected to be diffuse, and in this sense the upper layer is not entirely sediment-free. But as elaborated below, the thickness of the upper layer tends to grow much more rapidly in the downstream direction than the lower layer, so that the layer-averaged concentration in the upper layer tends to vanish.
To illustrate how the two-layer formulation overcomes the shortcomings of the 3-equation model, it is useful to cast (3.1)–(3.5) into the form for steady, gradually varied flow corresponding to (2.7)–(2.9). The relations for the lower layer are
\begin{align} \frac {{{\delta _L}}}{{{U_L}}}\frac {{{\rm d}{U_L}}}{{{\rm d}x}} & = \frac {1}{{{ {1 - Ri} }}} \left [-\left ( {1 + \frac {1}{2}Ri} \right ){e_{ws}}\frac {{ {{U_L} - {U_U}} }}{{{U_L}}} + \left ( {1 + \frac {1}{2}Ri} \right )\frac {{{v_s}}}{{{U_L}}} \right .\nonumber\\ & \quad \left .{}+ Ri\,S - {C_{fb}} - {C_{fi}}\frac {{{\lvert {U_L} - {U_U}\rvert\, {\left ( {{U_L} - {U_U}} \right )}}}}{{U_L^2}} - \frac {1}{2}\frac {{{v_s}}}{{{U_L}}}r\,R{i_{}}\left ( {\frac {{{q_{se}}}}{{{q_s}}} - 1} \right ) \right ], \end{align}
\begin{align} \frac {{{\rm d}{\delta _L}}}{{{\rm d}x}} & = \frac {1}{{{ {1 - Ri} }}}\left [ \left ( {2 - \frac {1}{2}Ri} \right ){e_{ws}}\frac {{ {{U_L} - {U_U}} }}{{{U_L}}} - \left ( {2 - \frac {1}{2}Ri} \right )\frac {{{v_s}}}{{{U_L}}} \right . \nonumber\\ & \quad \left . {}- Ri\,S + {C_{fb}} + {C_{fi}}\frac {{{\lvert {{U_L} - {U_U}}\rvert\, {\left ( {{U_L} - {U_U}} \right )}}}}{{U_L^2}} + \frac {1}{2}\frac {{{v_s}}}{{{U_L}}}r\,R{i_{}}\left ( {\frac {{{q_{se}}}}{{{q_s}}} - 1} \right ) \right ], \end{align}
\begin{align} \frac {\delta _L }{{{q_s}}}\frac {{{\rm d}{q_s}}}{{{\rm d}x}} = \frac {{{v_s}}}{{{U_L}}}r\left ( {\frac {{{q_{se}}}}{{{q_s}}} - 1} \right ). \end{align}
In (3.7a
) and (3.7b
),
$Ri$
is a bulk Richardson number based on the lower layer:
\begin{equation} Ri = \frac {{RgC{\delta _L}}}{{U_L^2}} = \frac {{Rg{q_s}}}{{U_L^3}}. \end{equation}
The corresponding equations for the upper layer are
\begin{align} \frac {{{\delta _U}}}{{{U_U}}}\frac {{{\rm d}{U_U}}}{{{\rm d}x}} = - {e_{wo}} + {e_{ws}}\frac {{ {{U_L} - {U_U}} }}{{{U_U}}} - \frac {{{v_s}}}{{{U_U}}} + \frac {{{C_{fi}}\,\lvert {U_L^{} - U_U^{}}\rvert\, {{\left ( {U_L^{} - U_U^{}} \right )}}}}{{U_U^2}},\\[-10pt]\nonumber \end{align}
\begin{align} \frac {{{\rm d}{\delta _U}}}{{{\rm d}x}} = 2{e_{wo}} - 2{e_{ws}}\frac {{ {{U_L} - {U_U}} }}{{{U_U}}} + 2\frac {{{v_s}}}{{{U_U}}} - \frac {{{C_{fi}}\,\lvert {U_L^{} - U_U^{}}\rvert\, {{\left ( {U_L^{} - U_U^{}} \right )}}}}{{U_U^2}}. \end{align}
As in the case of the single-layer model, when the flow in the lower layer is Froude supercritical, i.e.
$Ri \lt 1$
, the above equations can be integrated downstream from specified upstream values of
$U_L$
,
$\delta _L$
,
$q_s$
,
$U_U$
and
$\delta _U$
.
3.1. Normal flow for bypass conditions
To show how the introduction of settling detrainment affects the behaviour of a turbidity current, we consider the case of bypass flow, so that (3.7c ) is replaced with (2.12). Accordingly, (3.7a ) and (3.7b ) reduce to
\begin{align} &&\begin{split} \frac {{{\delta _L}}}{{{U_L}}}\frac {{{\rm d}{U_L}}}{{{\rm d}x}} & = \frac {1}{{{ {1 - Ri} }}} \left [-\left ( {1 + \frac {1}{2}Ri} \right ){e_{ws}}\frac {{ {{U_L} - {U_U}} }}{{{U_L}}} + \left ( {1 + \frac {1}{2}Ri} \right )\frac {{{v_s}}}{{{U_L}}} \right .\\[7pt] & \quad \left .{}+ Ri\,S - {C_{fb}} - {C_{fi}}\frac {{{\lvert {U_L} - {U_U}\rvert\, {\left ( {{U_L} - {U_U}} \right )}}}}{{U_L^2}} \right ], \end{split}\\[-10pt]\nonumber \end{align}
\begin{align} &&\begin{split} \frac {{{\rm d}{\delta _L}}}{{{\rm d}x}} & = \frac {1}{{{ {1 - Ri} }}}\left [ \left ( {2 - \frac {1}{2}Ri} \right ){e_{ws}}\frac {{ {{U_L} - {U_U}} }}{{{U_L}}} - \left ( {2 - \frac {1}{2}Ri} \right )\frac {{{v_s}}}{{{U_L}}} \right . \\[7pt] & \quad \left .{} - Ri\,S + {C_{fb}} + {C_{fi}}\frac {{{\lvert {{U_L} - {U_U}}\rvert\, {\left ( {{U_L} - {U_U}} \right )}}}}{{U_L^2}} \right ]. \end{split} \end{align}
Just as in the case of the 3-equation bypass model and Ellison & Turner (Reference Ellison and Turner1959), these equations have a normal flow solution. Setting
${\rm d}U_L/{\rm d}x = 0$
in (3.10a
) and
${\rm d}U_U/{\rm d}x = 0$
in (3.9a
), it is possible to solve for the constant normal values
$U_{Ln}$
and
$U_{Un}$
. From these values, it can be found that (3.10b
) and (3.9b
) reduce to the forms
\begin{align} \frac {{{\rm d}{\delta _L}}}{{{\rm d}x}} = {A_L}, \end{align}
\begin{align} \frac {{{\rm d}{\delta _U}}}{{{\rm d}x}} = {A_U}, \end{align}
where
$A_L$
and
$A_U$
are constants obtained from the right-hand sides of (3.10b
) and (3.9b
). In summary, at normal flow, the velocities of both the lower and upper layers are constant, and the layer thicknesses of the lower and upper layers increase linearly downstream. We find below that the effect of settling renders the downstream growth rate of the lower layer much less than the upper layer, indeed so much less that a turbidity current is likely able to track (and thus potentially make) its own channel.
The two-layer model is compared with the original one-layer model of Ellison & Turner (Reference Ellison and Turner1959) for the case of vanishing fall velocity in Appendix C using sample laboratory-scale input parameters. It is shown there that the results of the two models are in general agreement over a downstream scale of 10 m, but show divergent behaviour at a scale of 200 m. We suggest here that this divergence can be attributed to the more complete physical description of the two-layer model.
The relations for
$U_{Ln}$
and
$U_{Un}$
can be cast in dimensionless form using
\begin{align} R{i_n} = \frac {{Rg{q_s}}}{{U_{Ln}^3}} ,\quad R{i_{In}} = \frac {{Rg{q_s}}}{{U_{Ln}^3\left ( {1 - {\Gamma }} \right )^2}},\quad \Gamma = \frac {{{U_{Un}}}}{{{U_{Ln}}}},\quad {\tilde v_s} = \frac {{{v_s}}}{{{{\left ( {Rg{q_s}} \right )}^{1/3}}}} \end{align}
to give
\begin{align} &\left ( {1 + \tfrac {1}{2}R{i_n}} \right )\left [ {{{\tilde v}_s}Ri_n^{1/3} - {e_{ws}}\left [ {R{i_{In}}} \right ]\left ( {1 - \Gamma } \right )} \right ] + R{i_n}\,S \nonumber\\ &\quad{}- {C_{fb}} - 2{e_{ws}}[R{i_{In}}]{\left ( {1 - \Gamma } \right )^2} = 0,\\[-10pt]\nonumber \end{align}
\begin{align} & -{e_{wo}}{\Gamma ^2} + {e_{ws}}\left [ {R{i_{In}}} \right ]\Gamma \left ( {1 - \Gamma } \right ) - {\tilde v_s}\,R{i^{1/3}}\,\Gamma + 2{e_{ws}}\left [ {R{i_{In}}} \right ]\left ( {1 - {\Gamma }} \right )^2 = 0. \end{align}
From (3.13a
,
b
), it can be found that once the three parameters
$S$
,
$\tilde v_s$
and
$C_{fb}$
are specified, the two variables
$Ri_n$
and
$\Gamma$
can be determined. These values in turn set the dimensional values
$U_U$
and
$U_L$
. Thus
$U_U$
and
$U_L$
at the normal flow condition do not depend on boundary conditions. Here, without losing generality, we set
$C_{fb}= 0.002$
, which is typical for fine-grained channels (Konsoer et al. Reference Konsoer, Zinger and Parker2013; Ma et al. Reference Ma, Nittrouer, Naito, Fu, Zhang, Moodie, Wang, Wu and Parker2017, Reference Ma2020; Simmons et al. Reference Simmons, Azpiroz-Zabala, Cartigny, Clare, Cooper, Parsons, Pope, Sumner and Talling2020). A wide range of values of
$S$
and
$\tilde v_s$
, compatible with submarine channels (Covault et al. Reference Covault, Fildani, Romans and McHargue2011), are chosen to illustrate solutions to the normal flow condition obtained from solving (3.13a
,
b
). The results are shown in figures 4(a–d). The pattern of normal flow solutions divides into a Froude-supercritical regime defined in terms of the lower layer (
$Ri \lt 1$
) for sufficiently large values of
$S$
and
$\tilde v_s$
, and a Froude-subcritical regime (
$Ri \gt 1$
) as
$S$
and
$\tilde v_s$
become small. A clear threshold behaviour can be identified: when
$S \gt 0.0063$
or
${\tilde v_s} \gt 0.0042$
, the flow is always Froude-supercritical regardless of the value of the other parameter. Assuming
$q_s = 0.6$
m
$^2$
s–1 (a value justified below) as an example in the computation of
$\tilde v_s$
, three lines corresponding to grain sizes
$D=31.25$
, 62.5, 125
$\unicode{x03BC}$
m (specific gravity of quartz, so
$R = 1.65$
and water temperature 20
$^\circ$
C) are plotted in all of figures 4(a–d).
Bulk Richardson number of the lower layer
$Ri$
and the velocity ratio
$\Gamma$
of
$U_{U}$
to
$U_{L}$
solved from (3.13a
,
b
) for various combinations of channel slope
$S$
and dimensionless settling velocity
$\tilde{v}_s$
at bypass normal flow. (a) Plot of
$Ri$
as a function of
$S$
and
$\tilde{v}_s$
. Since the lower-layer Froude number is
$Fr_d=1/\sqrt {Ri}$
, the isoline
$Ri=1$
separates the Froude-supercritical and -subcritical flow regimes. A clear threshold behaviour can be identified: when
$S \gt 0.0063$
or
$\tilde{v}_s\gt 0.0042$
, the flow is always Froude-supercritical regardless of the value of the other parameter. (b) Plot of
$\Gamma =U_{Un}/U_{Ln}$
as a function of
$S$
and
$\tilde{v}_s$
. Note that the upper layer is always slower than the lower layer; this effect strengthens as
$S$
and
$\tilde{v}_s$
become small. (c) Plot of
$A_L$
as a function of
$S$
and
$\tilde{v}_s$
. There is a neutral line where water entrainment due to turbulent mixing and water detrainment due to sediment settling zeros out. Below the line where
$S$
is small and
$\tilde{v}_s$
is large, the turbidity currents may subside (negative thickening rate) due to sediment-settling induced drop in the level of the interface. (d) Plot of
$A_U$
as a function of
$S$
and
$\tilde{v}_s$
;
$A_U$
is at least one order of magnitude larger than
$A_L$
because the ambient water entrainment coefficient
$e_{w0}=0.075$
sets the top interface boundary condition for the upper layer, which corresponds to the upper bound for this coefficient. Assuming
$q_s=0.6$
m
$^2$
s–1, three lines corresponding to
$D=31.25$
, 62.5, 125
$\unicode{x03BC}$
m are shown in all plots.
3.2. Calculations at field scale under bypass conditions
Advances in field measurements have shown that turbidity currents come in a variety of shapes and sizes (Xu et al. Reference Xu, Noble and Rosenfeld2004; Dorrell et al. Reference Dorrell, Darby, Peakall, Sumner, Parsons and Wynn2014; Hughes Clarke Reference Clarke and John2016; Paull et al. Reference Paull2018; Pope et al. Reference Pope2022; Talling et al. Reference Talling2022), depending mainly on the sizes of the systems and the transported grain sizes. Long-runout flows are known to have emanated from the Gaoping Canyon and the Grand Banks Canyon, where breakages of submarine telecommunication cables have provided indications of peak velocities at approximately 15–20 m s–1 (Heezen & Ewing Reference Heezen and Ewing1952; Hsu et al. Reference Hsu, Kuo, Lo, Tsai, Doo, Ku and Sibuet2008). More recently, measurements of long-runout flows have been made in the Congo Canyon, where flows can accelerate for over 1200 km to reach peak velocities 8 m s–1 upon reaching the abyssal plain at water depth over 5 km (Talling et al. Reference Talling2022). Unfortunately, this flow was so powerful that it destroyed all the instrumentation, consequently there are no flow discharge measurements from this event. More detailed velocity measurements of smaller flows in the Congo Canyon show that at
$\sim150$
km offshore and water depth almost 2 km, turbidity current events have peak discharges of up to
$\sim16\,000$
m
$^3$
s–1, and typically last for approximately a week (Azpiroz-Zabala et al. Reference Azpiroz-Zabala, Cartigny, Talling, Parsons, Sumner, Clare, Simmons, Cooper and Pope2017). Following the passage of the faster head of the flow, a flow speed of approximately 0.75 m s–1 is typically maintained for 5 days, but occasionally for up to 8 days (Simmons et al. Reference Simmons, Azpiroz-Zabala, Cartigny, Clare, Cooper, Parsons, Pope, Sumner and Talling2020). However, these detailed measurements were taken opportunistically, so are unlikely to be channel-forming turbidity currents. For the rarer and more powerful channel-forming flows, we refer to the reconstruction of channel-forming turbidity currents by Konsoer et al. (Reference Konsoer, Zinger and Parker2013).
Konsoer et al. (Reference Konsoer, Zinger and Parker2013) reconstructed channel-forming flows in turbidity currents by means of an approximate matching of current driving force with rivers. They offer two estimates each for mean sediment concentration
$C$
and the bed resistance
$C_{bf}$
. Of these, we choose
$C_u=0.006$
and
$C_{bf}=0.002$
, where
$C_u$
is the upstream boundary condition for
$C$
. The levee-to-levee channel width of the Amazon Submarine Channel at
$x=270$
km is found to be approximately 2 km in figure 3(d) of Pirmez & Imran (Reference Pirmez and Imran2003). The channel-forming discharge at this width can be estimated to be 200 000 m
$^3$
s–1 according to figure 9 of Konsoer et al. (Reference Konsoer, Zinger and Parker2013). Therefore, we assume a water discharge per unit width
$q_{wu}$
at the upstream end of our calculation of 100 m
$^2$
s–1. Insofar as we consider bypass currents, we hold the volume sediment discharge per unit width
$q_s$
at the constant value
$100\,\text{m}^2\,\text{s}^{-1}\times 0.006 = 0.6\,\text{m}^2\,\text{s}^{-1}$
. This is the justification for using
$q_s=0.6$
m
$^2$
s–1 in figure 4.
We now show numerical results for the spatial, down-canyon development of a bypass turbidity current. All the calculations shown below assume a size 62.5
$\unicode{x03BC}$
m for the suspended sediment, a bed friction coefficient
$C_{fb} = 0.002$
, a suspended sediment transport rate per unit width 0.6 m
$^2$
s–1, an upstream flow discharge
$q_{wu}=100$
m
$^2$
s–1, and an upstream volume suspended sediment concentration
$C_u=0.006$
. We consider two cases for slope: one with constant slope 0.03, and one with a slope that exponentially declines downstream in approximate concordance with the long profile of figure 2(b). The numerical results were obtained by solving (3.10a
,
b
) and (3.9a
,
b
). The cases below are for the Froude-supercritical condition, which allows a stepwise downstream numerical solution to (3.9) and (3.10) using the Euler step method. Sample numerical results under a Froude-subcritical condition can be found in figure 9 of Appendix D.
(a) Spatial evolution of the lower layer and upper layer velocities
$U_L$
and
$U_U$
, over a 5 km reach, starting from two sets of upstream conditions. In all cases, the velocities evolve towards normal flow. (b) Spatial evolution of thicknesses of the lower and upper layers
$\delta _L$
and
$\delta _U$
over a 200 km reach. (c) Spatial evolution of lower and upper layer thicknesses
$\delta _L$
and
$\delta _U$
over a 10 km reach, using two different sets of upstream conditions.
Figure 5 shows the downstream development of a bypass flow over a constant slope
$S=0.03$
. This slope has been chosen insofar as it is representative of the constant-slope channel profiles of Covault et al. (Reference Covault, Fildani, Romans and McHargue2011) shown in figure 2(a). In figure 5, two Froude-supercritical upstream conditions have been chosen for the numerical solution of (3.9a
,
b
) and (3.10a
,
b
):
$(U_L, U_U, \delta _L, \delta _U)$
= (2.5 m s–1, 1.25 m s–1, 40 m, 1 m) and (4.5 m s–1, 2.25 m s–1, 22.22 m, 1 m).
In figure 5(a), velocities converge to the normal values
$(U_L, U_U)$
= (3.083 m s–1, 0.853 m
s
–
1) within approximately 5 km. It is not necessary to confirm that these correspond to long-runout values; on a constant slope, they would not change even hundreds of kilometres downslope. The result
$U_L \gt U_U$
confirms that the lower layer is the driving layer, and the upper layer is the driven layer.
Figure 5(b) shows the development of the layer thicknesses
$\delta _L$
and
$\delta _U$
out to 200 km. By 200 km,
$\delta _U$
has thickened to over 13 000 m, an unreasonable value in line with the overthickening of the original 3-equation model. This issue is considered in more detail in § 4. The lower layer, on the other hand, has thickened to only 520 m. This is comparable with channel depth 165 m, and a levee crest to back-levee elevation difference of at least 270 m at the Shepard Bend of Monterey Channel (Fildani et al. Reference Fildani, Normark, Kostic and Parker2006), which is approximately 140 km down-channel of the canyon head. This system has a mean down-channel slope close to the value 0.03 assumed here (Covault et al. Reference Covault, Fildani, Romans and McHargue2011). Figure 5(c) is identical to figure 5(b) except that the spatial domain has been reduced to 10 km. The results clearly show that in the two-layer model, the lower layer thickens downstream at a much slower rate (factor 0.038) than the upper layer, whereas the upper layer thickens at a rate higher than that predicted by the single-layer 3-equation model.
Figure 6 provides a comparison of the spatial evolution predicted by three models: the two-layer model described here in (3.9a
,
b
) and (3.10a
,
b
), the original 3-equation model in (2.1), (2.2) and (2.12) (Fukushima et al. Reference Fukushima, Parker and Pantin1985; Parker et al. Reference Parker, Fukushima and Pantin1986), and the 3-equation model modified to include detrainment in (2.1), (2.16) and (2.3) (Toniolo et al. Reference Toniolo, Lamb and Parker2006a
; Bolla Pittaluga et al. Reference Bolla Pittaluga, Frascati and Falivene2018). Again, the bed slope
$S$
is held constant at 0.03. All calculations use one of the sets of upstream conditions of figure 5. Figure 6(a) shows spatial evolution of velocity over a 5 km reach. The results for
$U_L$
of the two-layer model,
$U$
of the 3-equation model, and
$U$
of the 3-equation model modified to include detrainment, all show similar spatial evolution, and approach nearly the same normal velocity. Figure 6(b) shows spatial evolution over a 200 km reach of
$\delta _L$
of the two-layer model,
$\delta$
of the 3-equation model, and
$\delta$
of the 3-equation model modified to include detrainment. The predictions for
$\delta$
from both versions of the 3-equation model at 200 km are greatly in excess of that predicted for
$\delta _L$
by the two-layer model. In figure 7, a current running down a long, concave upward profile is considered. Slope declines downstream in accordance with an exponential law that is an approximate fit to figure 2(b) (Amazon Submarine Channel) over 800 km:
\begin{equation} S = {S_u}\,{{\rm e}^{ - (x/{x_e})}}, \end{equation}
where
$x$
and
$x_e$
are in km, with
$x_e$
= 265.8 km and
$S_u$
= 0.0166.
Comparison of spatial development on the slope
$S = 0.03$
. (a) Spatial evolution over a 5 km reach of lower-layer velocity
$U_L$
and upper-layer velocity
$U_U$
of the two-layer model, velocity
$U$
of the 3-equation model, and velocity
$U$
of the 3-equation model modified to include detrainment. (b) Spatial evolution over a 200 km reach of lower-layer thickness
$\delta _L$
of the two-layer model, thickness
$\delta$
of the 3-equation model, and thickness
$\delta$
of the 3-equation model modified to include detrainment.
Bypass calculations based on a simplified profile of the Amazon Canyon–Fan system, using 62.5
$\unicode{x03BC}$
m suspended sediment over a 400 km reach. (a) Spatial evolution of velocities
$U_L$
and
$U_U$
for the two-layer model,
$U$
for the 3-equation model, and
$U$
for the 3-equation model modified to include detrainment. The slope profile is also shown. (b) Spatial evolution of thicknesses
$\delta _L$
and
$\delta _U$
for the two-layer model,
$\delta$
for the 3-equation model, and
$\delta$
for the 3-equation model modified to include detrainment. (c) Densimetric Froude number
$Fr_d$
for the lower layer of the two-layer model, the 3-equation model, and the 3-equation model modified to include detrainment. (d) Spatial evolution of water discharge per unit width
$q_w$
for the lower layer of the two-layer model, the 3-equation model, and the 3-equation model modified to include detrainment. (e) Spatial evolution of the suspended sediment concentration
$C$
in the lower layer of the two-layer model, the 3-equation model, and the 3-equation model modified to include detrainment.
Figure 7(a) shows the downstream evolution over a 400 km reach of
$U_L$
and
$U_U$
of the two-layer model, and
$U$
of the original 3-equation model and the version modified for detrainment. We terminate the calculation where the Froude number declines to the Froude-critical value (
$Fr_{d}=1$
); a hydraulic jump may occur upstream of this point, depending on downslope conditions. It can be seen that both versions of the 3-equation model reach the condition
$Fr_{d}=1$
at distances shorter than 400 km. The two-layer model reaches
$Fr_{d}=1$
at 402 km. The results for
$U_L$
compare well with those for
$U$
of the two versions of the 3-equation model, with values declining to 2.14 m s–1 at
$x$
= 400 km. The predicted values of
$U_U$
are uniformly lower than
$U_L$
, again indicating that the upper layer is driven by the lower layer. Also shown in the diagram is the slope profile
$S = 0.0037$
at
$x$
= 400 km. Figure 7(b) shows the corresponding results for
$\delta _L$
and
$\delta _U$
, and also
$\delta$
predicted by the two versions of the 3-equation model. The predicted values of
$\delta$
of the 3-equation models are far too high to follow any channel so far down the system. The predicted value of
$\delta _L$
, on the other hand, is 250 m at
$x=200$
km, a value that compares reasonably with estimates of channel bankfull depth (see below) (Pirmez & Imran Reference Pirmez and Imran2003; Fildani et al. Reference Fildani, Normark, Kostic and Parker2006).
Figure 7(c) shows the long profiles of the Froude number
$Fr_d\ (= Ri^{-1/2)})$
predicted for the lower layer of the two-layer model and the two versions of the 3-equation model. The Froude number declines downstream towards unity in all three cases. In the case of the 3-equation model, critical flow is attained at
$x=260$
km and 380 km, respectively. In the two-layer model, it is attained at
$x=402$
km. As noted above, the implication is that a hydraulic jump may occur somewhat upstream of this point. The spatially varying model encompassed in (3.9a
,
b
) and (3.10a
,
b
) cannot capture hydraulic jumps. A shock-fitting solution to the primitive equations (3.1), (3.2), (3.4) and (3.5) would, however, capture them (Fildani et al. Reference Fildani, Normark, Kostic and Parker2006).
Figures 7(d) and 7(e) respectively show the down-channel evolution of water discharge per unit width
$q_w$
and suspended sediment concentration
$C$
. Note that
$q_w$
reaches a maximum and then declines after
$x=328$
km, and
$C$
reaches a minimum and then increases after this point. These extreme points are because
$U_L$
can quickly adjust to the local equilibrium value, which decreases with exponentially declining channel slope. This flow slowdown effect dominates farther downstream, where
$\delta _L$
increases more slowly than linear. As a result,
$q_w= U_L\delta _L$
declines downstream of its peak value. Concentration
$C$
has a minimum at the same location because
$C=q_s/q_w$
and
$q_s$
is constant for the bypass condition.
A lower-layer thickness
$\delta _L=250$
m at
$x=200$
km compares well with the observed channel depth of approximately 110–170 m in the Amazon system (Pirmez & Imran Reference Pirmez and Imran2003). It should be expected that flow thickness exceeds channel depth, but is of the same order of magnitude in order to construct the channel and its levees. The present model is 1-D, which means that it corresponds to flow between frictionless vertical walls. In actual submarine channels, flow stripping (i.e. the overflow that builds the levees and confines the channel) should cause streamwise flow discharge to decline downstream. A possible way to incorporate this process into a model of long-runout turbidity currents is presented by Spinewine et al. (Reference Spinewine, Sun, Babonneau and Parker2011) and later by Bolla Pittaluga et al. (Reference Bolla Pittaluga, Frascati and Falivene2018).
4. Discussion
The two-layer model of bypass turbidity currents presented here offers many avenues for future development, allowing us to extend our understanding of the fluid dynamics and morphodynamics of long-runout turbidity currents and the morphologies they create. We enumerate a few of these below.
The example flows that are modelled here are limited to 1-D steady, Froude-supercritical flows that develop in the downstream direction. By definition, such a steady flow that has run out 1000 km must be continuously occupying the channel for 1000 km. Measurements in the Congo Submarine Channel indicate that turbidity current events can last for a week or more in the proximal part of the system (Azpiroz-Zabala et al. Reference Azpiroz-Zabala, Cartigny, Talling, Parsons, Sumner, Clare, Simmons, Cooper and Pope2017) and for up to several weeks in the distal part of the system (Baker et al. Reference Baker2024). As the head of the flows outruns the rest of the flow, these flows are likely to stretch to hundreds of kilometres and last for weeks in the distal part of the system. Such stretching flows can be modelled, at least in part, by abandoning the steady, gradually varied flow assumption, and instead solving the full time-varying equations (3.1), (3.2), (3.4) and (3.5), and for bypass flows, (2.12). An appropriate shock-capturing numerical technique such as the one used by Kostic & Parker (Reference Kostic and Parker2006) and Cao et al. (Reference Cao, Li, Pender and Liu2015) can model unsteady flow, whether Froude-supercritical or Froude-subcritical. It was used by Kostic & Parker (Reference Kostic and Parker2006) to reproduce the migrating turbidity current head and hydraulic jump of one of the experiments of Garcia & Parker (Reference Garcia and Parker1989). The same formulation could presumably be used to model the hydraulic jumps observed in the field by Sumner et al. (Reference Sumner, Peakall, Parsons, Wynn, Darby, Dorrell, McPhail, Perrett, Webb and White2013) and Clarke & Jhon (Reference Clarke and John2016).
A bypass current cannot be used directly to model the morphodynamics of bed evolution. Morphodynamics can, however, be modelled by implementing the full forms of (3.1)–(3.5), along with the Exner equation of bed sediment conservation. That is, where
$\eta$
is bed elevation,
$\lambda$
is bed porosity and
$q_b$
is the volume rate of bedload transport per unit width, we have
\begin{equation} (1 - \lambda )\frac {{\partial \eta }}{{\partial t}} = - \frac {{\partial {q_b}}}{{\partial x}} + {v_s}\left ( {rC - {E_s}} \right ) . \end{equation}
Appropriate closure assumptions are necessary for
$E_s$
,
$r$
and
$q_b$
. For example, Parker et al. (Reference Parker, Garcia, Fukushima and Yu1987) and Garcia & Parker (Reference Garcia and Parker1991) present closures for
$E_s$
and
$r$
, and a closure for
$q_b$
that is valid up to and including the regime of unidirectional bedload sheet flow is given in Ribberink (Reference Ribberink1998). Among the various submarine phenomena that can be revisited with a morphodynamic two-layer formulation are field observations of accelerating flow on constant and decreasing slopes (Talling et al. Reference Talling2022), flows that grow rapidly in sediment volume by a factor 100–1000 (Pope et al. Reference Pope2022; Böttner et al. Reference Böttner, Stevenson, Englert, Schönke, Pandolpho, Geersen, Feldens and Krastel2024), and the formative conditions for large trains of knickpoints (Heijnen et al. Reference Heijnen2020) or smaller trains of upstream-migrating crescentic bedforms (Clarke & Jhon Reference Clarke and John2016). A morphodynamic version of the two-layer model can also be adapted to model the incision necessary to excavate submarine canyons (Zhang et al. Reference Zhang, Parker, Izumi, Cartigny, Li and Wang2017) with the aid of, for example, the sandblasting model of Lamb et al. (Reference Lamb, Dietrich and Sklar2008).
The 3- and 4- equation models have, after extension to a 2-D streamwise-lateral form, been used to explain self-channelization of turbidity currents via levee emplacement (Imran et al. Reference Imran, Parker and Katopodes1998; Halsey & Kumar Reference Halsey and Kumar2019) and the emplacement of channelized submarine fans (Wahab et al. Reference Wahab, Hoyal, Shringarpure and Straub2022). Due to the limitations of these models, such features have been successfully modelled out to only approximately 5–25 km. The new two-layer model, upon extension to two dimensions, offers the possibility of modelling levee emplacement over most of the length of a long-runout turbidity current path. Built into such a model would be a characterization of flow stripping (channel overflow; Spinewine et al. Reference Spinewine, Sun, Babonneau and Parker2011), which would cause flow discharge to decline downstream and bring the flow thickness more in line with the observed levee elevations. A 2-D formulation could also be used to extend the work of Imran et al. (Reference Imran, Parker and Pirmez1999) on meandering submarine channels with self-formed levees, in which case sinuosity might abet mixing between the two layers (Straub & Mohrig Reference Straub and Mohrig2008).
Field measurements have shown that the dynamics of the migrating front or nose of a turbidity current may be too complex to be modelled even using a version of the two-layer model that captures front behaviour, as sediment concentrations can be above 10 % (Paull et al. Reference Paull2018; Wang et al. Reference Wang, Xu, Talling, Cartigny, Simmons, Gwiazda, Paull, Maier and Parsons2020). A first step in overcoming this issue is suggested by the model of Spinewine & Capart (Reference Spinewine and Capart2013) for intense flow and sediment transport at the nose of a dam-break flow. It may be possible to adapt this formulation to the front of a turbidity current otherwise modelled by the present two-layer formulation.
The two-layer model presented here represents an extension of the 3-equation model of Fukushima et al. (Reference Fukushima, Parker and Pantin1985) and Parker et al. (Reference Parker, Fukushima and Pantin1986). Cao et al. (Reference Cao, Li, Pender and Liu2015) referred to their model as a ‘double-layer averaged model’. Their model includes a turbidity current layer and an ambient water layer, the dynamics of which must be considered in the shallow setting of the reservoir that they model. In that sense, our model is a ‘three-layer model’: driving layer, driven layer and ambient layer. In our model, we take the ambient water to be infinitely deep, so can treat it as stagnant, which enables the description of the deep sea environment. Parker et al. (Reference Parker, Fukushima and Pantin1986) also include a 4-equation model, where the extra equation accounts for the balance of turbulent kinetic energy. In principle, the extension of kinetic energy balance to the two-layer model is straightforward, adding one extra equation each to the formulation for the lower and upper layers. It may be useful to revisit the formulation in light of the results of Fay (Reference Fay2012).
Numerical methods that resolve the upward normal structure of the flow, such as
$k{-}\epsilon$
, LES or DNS, may not be feasible to implement for long-runout turbidity currents. They nevertheless could be used to develop refined closures for the two-layer model that enhance its performance and accuracy.
While the entrainment coefficient
$e_{wo}=0.075$
for unstratified turbulent flow is well justified by data, it may not apply to upper layers that are predicted by the present model to become thousands of metres thick (e.g. as shown in figure 5
b). There are a number of potential reasons why
$e_{wo}$
might not attain such a high value. First, a semi-empirical, fully turbulent entrainment coefficient
$e_{wo}= 0.075$
corresponds to the limit of a plane free jet (
$Ri\to 0$
; Parker et al. Reference Parker, Garcia, Fukushima and Yu1987). The rate of production of turbulent energy for such a flow should scale as
${\rm d}u/{\rm d}z$
, where
$z$
is the upward normal coordinate, and
$u$
is local streamwise velocity averaged over turbulence. As the upper layer becomes thicker and thicker, the term
${\rm d}u/{\rm d}z$
may drop to the point where full turbulence can no longer be maintained. Were an entirely sediment-free upper layer flow to become fully laminar,
${\rm d}U$
would scale as
$x^{1/2}$
at normal flow, in accordance with the Prandtl result for laminar flow over a flat plate, rather than the turbulent scaling
${\rm d}U \sim x^1$
used here (corresponding to constant bed resistance coefficient
$C_{fb}$
). Second, we have adopted a relation between entrainment rate
$e_{ws}$
and Richardson number
$Ri_I$
at the interface between the lower and upper layers, again assuming fully turbulent flow (Parker Reference Parker1982; Parker et al. Reference Parker, Garcia, Fukushima and Yu1987; Johnson & Hogg Reference Johnson and Hogg2013). Especially for flow in the Froude-subcritical range, however, there are conditions under which stratification is so strong at the lower/upper interface that turbulence is extinguished there (Dorrell et al. Reference Dorrell, Peakall, Darby, Parsons, Johnson, Sumner, Wynn, Özsoy and Tezcan2019; Marshall et al. Reference Marshall, Dorrell, Keevil, Peakall and Tobias2021; Salinas et al. Reference Salinas, Balachandar and Cantero2021a
; Lloyd et al. Reference Lloyd, Dorrell and Caulfield2022). Under such conditions, mixing at both lower/upper and upper/ambient interfaces is likely to be governed by laminar processes, and can thus be expected to be much weaker than that predicted by (3.6). Information in Arita & Jirka (Reference Arita and Jirka1987a
) suggests that the condition for the domination of entrainment by laminar effects becomes more stringent with increasing Reynolds number of the lower layer. Third, the value
$e_{wo}=0.075$
is for a 2-D (streamwise-upward normal) plane jet where the flow is not allowed to spread in the third, transverse dimension. In a field case, however, the upper layer loses confinement and turns three-dimensional when it overspills the submarine levee, reducing vertical entrainment (Rajaratnam Reference Rajaratnam1976) and causing direct loss of water, e.g. flow stripping (Fildani et al. Reference Fildani, Normark, Kostic and Parker2006; Spinewine et al. Reference Spinewine, Sun, Babonneau and Parker2011). It is this transverse spreading and overspilling, rather than the failure of (2.4) and (3.6) themselves, that is the likely reason why Traer et al. (Reference Traer, Fildani, McHargue and Hilley2015) implied that this relation overpredicts entrainment. This feature represents a limitation of the present 1-D formulation that can be overcome in a 2-D formulation. Moreover, although it is very dilute, there should still be some sediment in the upper layer entrained from the lower layer, which can also reduce the entrainment coefficient (Salinas et al. Reference Salinas, Cantero, Shringarpure and Balachandar2019b
). The present model is 1-D in nature and cannot capture lateral expansion of both the lower and upper layers, and corresponding effects that would help to limit the streamwise increase in flow thickness. The two-layer model may thus merit modifications in light of the above comments.
The formula for water entrainment coefficient
$e_{ws}$
(2.4) was developed empirically (Parker et al. Reference Parker, Garcia, Fukushima and Yu1987) based on laboratory-scale flume experiments (
$\sim10$
m) of both density currents with saline or clay-silt particles (negligible settling velocity) and particulate turbidity currents with non-negligible settling velocity. No adjustment for the effect of detrainment on the data used to develop (2.4) was made in the analysis, because the relevance of detrainment was not generally recognized until the work of Toniolo et al. (Reference Toniolo, Lamb and Parker2006a
). Such an adjustment is warranted in the future; it would result in somewhat higher estimates of the entrainment coefficient
$e_{ws}$
for experiments of silty and sandy turbidity currents. The general form of the relation, however, is anchored towards the upper limit, lower limit and part of the middle range by experiments of particulate currents with negligible fall velocity, and as such, we use it for the present analysis. In the future, experiments on turbidity currents in long flumes (
$\sim100$
m) will be of great advantage to study the water detrainment effect, and thus two-layer structures of turbidity currents in a strong turbulent flow, in light of indiscrimination between density and particulate turbidity currents in small scale flumes (Appendix C).
The two-layer model with a single grain size can be adapted in a straightforward way to a three-layer model with two grain sizes. Choosing one of these grain sizes to be in the sand range and the other to be in the mud range can help to clarify the nature of morphodynamic channel–levee interaction (Deptuck & Sylvester Reference Deptuck and Sylvester2017). It can also further quantify the role of mud that has a relatively small settling velocity in maintaining sand that has a relatively high fall velocity, in suspension, so as to transport the sand farther downstream (Salaheldin et al. Reference Salaheldin, Imran, Chaudhry and Reed2000). More complex grain size distributions and patterns of dispersion can be considered in the future (Xie et al. Reference Xie, Hu, Zhu, Yu and Pähtz2023a ).
5. Conclusions
Turbidity currents are bottom density flows driven by the excess weight of suspended sediment. Turbidity currents in the deep sea are known to sculpt leveed channels that are hundreds or thousands of kilometres long. To construct such channels, a current must run out at least that far, consistently following the channel that is sculpted by the current itself. No existing model of turbidity current dynamics is capable of accomplishing this. Here, informed by the
$k{-}\epsilon$
model of Luchi et al. (Reference Luchi, Balachandar, Seminara and Parker2018), we quantify the problem in terms of a simpler two-layer model. The lower (driving) layer is where nearly all suspended sediment is sequestered. This sediment is assumed to have a single, constant settling velocity and to be transported as a dilute suspension. The upper layer is nearly sediment-free, and is dragged along by the lower (driving) layer. It is essential to introduce the concept of detrainment across the two-layer interface to quantify how sediment resists upward turbulent mixing via its fall velocity, an issue that was first studied systematically in the context of turbidity currents entering a minibasin (Toniolo et al. Reference Toniolo, Lamb and Parker2006a
). We apply the model to sediment bypass conditions, such that there is no net flux of sediment at the bed. The model presents a normal flow solution analogous to that of Ellison & Turner (Reference Ellison and Turner1959) for thermohaline bottom density flows. At normal flow, both lower and upper velocities
$U_L$
and
$U_U$
attain constant values, and both lower and upper layer thicknesses
$\delta _L$
and
$\delta _U$
increase linearly downstream. Under normal flow conditions, two thresholds are identified: one related to channel slope, and the other to the non-dimensional settling velocity. Exceeding either threshold causes the turbidity current to become Froude-supercritical. Although both layers thicken linearly downstream, we show that the effect of detrainment mediated by fall velocity dramatically slows the thickening rate of the lower layer. Our calculations using gradually varied flow show that the lower layer of a turbidity current can run out 400 km without overthickening to the point that it would lose track of its own channel. The calculations terminate there only because the flow reaches the Froude-critical condition, beyond which the computational method is no longer applicable. This issue can be overcome in the future by solving the parent unsteady, non-uniform version of the model (3.1)–(3.5). The model opens up further future avenues in the study of the fluid dynamics and morphodynamics of long-runout turbidity currents, including non-bypass flows, levee construction and the effect of multiple sediment sizes on grain-size-specific sediment runout.
Funding.
H.M. is supported by China’s MST Key R&D grant 2023YFC3206204; M.C. is supported by a Royal Society Research Fellowship (DHF∖R∖231008).
Declaration of interests.
The authors report no conflict of interest.
Author contributions.
G.P. conceived the study, developed the theoretical formulations and led the writing of the first version of the manuscript. H.M. obtained the exact solution and conducted the numerical simulations. E.V. and H.M. verified the formulations and analyses. G.P., M.C., S.B., E.V. and H.M. interpreted the results. H.M., G.P., M.C., S.B., E.V., X.F. and R.L. contributed to the discussion of the research and to the writing and revision of the manuscript.
Data availability statement.
The data/code that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix A. Derivation of (3.2)
The flow is incompressible turbulent and uniform in the transverse direction, and contains a dilute suspension of sediment with fall velocity
$v_s$
. It is assumed that the flow velocity averaged over turbulence is
$ (u, w )$
, where
$u$
is the velocity in the
$x$
direction, and
$w$
is the velocity in the
$z$
direction. The equation of continuity is
\begin{equation} \frac {{\partial u}}{{\partial x}} + \frac {{\partial w}}{{\partial z}} = 0. \end{equation}
The velocity
$ (u_s, w_s )$
of a sediment particle is taken to be
\begin{equation} ({u_s},{w_s}) = (u,w - {v_s}). \end{equation}
Integrating (A1) from
$z = 0$
to
$z = \delta _L$
yields
\begin{equation} \frac {\partial }{{\partial x}}\int _0^{{\delta _L}} {u\,{\rm d}z} - {\left . u \right |_{\delta _L}}\frac {{\partial {\delta _L}}}{{\partial x}} + {\left . w \right |_{{\delta _L}}} = 0 . \end{equation}
The settling interface is a material interface following the sediment, not the fluid. An appropriate version of the kinematic boundary condition for this case is
\begin{equation} \frac {{\partial {\delta _L}}}{{\partial t}} + {\left . u \right |_{{\delta _L}}}\frac {{\partial {\delta _L}}}{{\partial x}} - \left ( {{{\left . w \right |}_{{\delta _L}}} - {v_s}} \right ) = {u_e}, \end{equation}
where
$u_e$
is a turbulent entrainment velocity of fluid across the interface. Set
\begin{equation} {u_e} = {e_{ws}}({U_L} - {U_U}), \end{equation}
where
$e_{ws}$
is a coefficient of turbulent entrainment. Define
\begin{equation} {U_L}{\delta _L} = \int _0^{{\delta _L}} {u\,{\rm d}z}. \end{equation}
\begin{equation} \frac {{\partial {\delta _L}}}{{\partial t}} + \frac {{\partial {U_L}{\delta _L}}}{{\partial x}} = {e_{ws}}({U_L} - {U_U}) - {v_s} . \end{equation}
Appendix B. Derivation of (3.3)
The analysis follows from Appendix A. With
$c (t,x,z )$
denoting volume suspended sediment concentration averaged over turbulence, the relevant 2-D conservation equation is
\begin{equation} \frac {{\partial c}}{{\partial t}} + \frac {{\partial uc}}{{\partial x}} + \frac {{\partial wc}}{{\partial z}} - {v_s}\frac {{\partial c}}{{\partial z}} = - \frac {{\partial {F_R}}}{{\partial z}},\quad {F_R} = \overline {c'w'}, \end{equation}
where the primes denote fluctuating quantities, and
$F_R$
denotes the upward normal flux of suspended sediment due to turbulence. Integrating from
$z = 0$
to
$z = \delta _L$
and using (A4), it is found upon some reduction that
\begin{equation} \frac {{\partial {\delta _L}C}}{{\partial t}} + \frac {{\partial {U_L}{\delta _L}C}}{{\partial x}} - {\left . c \right |_{{\delta _L}}}\left ( {\frac {{\partial {\delta _L}}}{{\partial t}} + {{\left . u \right |}_{{\delta _L}}}\frac {{\partial {\delta _L}}}{{\partial x}} - {{\left . w \right |}_{{\delta _L}}} + {v_s}} \right ) = - {\left . {{F_R}} \right |_{{\delta _L}}} + {v_s}\left ( {E - rC} \right ), \end{equation}
where
\begin{equation} {C}{\delta _L} = \int _0^{{\delta _L}} {c\,{\rm d}z},\quad {U_L}{\delta _L}C = \int _0^{{\delta _L}} {uc\,{\rm d}z},\quad {F_{Rb}} = {v_s}E,\quad r = \frac {{{c_b}}}{C}, \end{equation}
and the subscript
$b$
denotes a near-bed value. Substituting (A4) and the closure hypothesis (B4) into (B2) reducing with the closure hypothesis, (B4), (B2) yields (3.3)
\begin{equation} {\left . {{F_R}} \right |_{{\delta _L}}} = {\left . c \right |_{{\delta _L}}}{w_e}, \end{equation}
Appendix C. Comparison of the two-layer model with the one-layer model of Ellison & Turner (Reference Ellison and Turner1959) for the case of a density flow with vanishing fall velocity
Here, we consider a comparison of the results of the two-layer model applied to the case of vanishing fall velocity with the corresponding version of the 3-equation one-layer model. In this case, the 3-equation model corresponds to that of Ellison & Turner (Reference Ellison and Turner1959), who performed their experiments in a laboratory channel with length 2 m. We run the models with the same values
$S = 0.03$
and
$C_{fb} = 0.002$
as the field-scale runs of figure 6, but with laboratory-scale values as follows:
$q_{wu} = 0.3$
m
$^2$
s–1,
$q_s = 0.01$
m
$^2$
s–1,
$U_{3u} = 0.8$
m s–1 and
$\delta _{3u} = 0.375$
m for the 3-equation, single-layer model; and
$q_{wu} = 0.3$
m
$^2$
s–1,
$q_s = 0.01$
m
$^2$
s–1,
$U_{Lu} = 0.8$
m s–1,
$U_{Uu} = 0.05$
m s–1,
$\delta _{Lu} = 0.375$
m and
$\delta _{Uu} = 0.1$
m for the two-layer model. The results are compared in figure 8.
Figures 8(a–c) show results for velocity, layer thickness and sediment concentration profiles for
$x \leq 10$
m. The results for the lower layer of the two-layer model and the 3-equation, single-layer model show relatively little difference at this scale. Figures 8(d–f) show the corresponding results for
$x \leq 200$
m, where differences are more apparent.
Numerical results for vanishing fall velocity under laboratory-scale conditions. (a–c) Flow velocity, layer thickness and sediment volumetric concentration at the laboratory scale (
$x\lt 10$
m). The difference between single-layer and two-layer models is small. (d–f) Flow velocity, layer thickness and sediment volumetric concentration at the far field (
$x\lt 200$
m). The difference between single-layer and two-layer models is more visible. The upstream boundary conditions are
$(U_L, U_U, \delta _L, \delta _U)$
= (0.8 m s–1, 0.05 m s–1, 0.375 m, 0.1 m).
Numerical results of long profiles of gradually varied (a) flow velocity and (b) layer thickness under Froude-subcritical conditions. The blue line represents the lower layer of the two-layer model. The black dashed line represents the original 3-equation model, and the black solid line represents the 3-equation model with the water detrainment term. The red dashed line represents the upper layer velocity. The solutions were obtained by integrating upstream from the downstream boundary
$(U_L, U_U, \delta _L, \delta _U)$
= (0.5 m s–1, 0.275 m s–1, 200 m, an arbitrary large value) at
$x=70\,000$
m.
Appendix D. Gradually varied flow under a Froude-subcritical condition
The numerical results under the Froude-subcritical condition are shown in figure 9. The channel slope is set one order of magnitude smaller than the Froude-supercritical case (figures 5and 6), and grain size is halved. Under the Froude-supercritical condition, the boundary condition is given downstream, and (3.10 a , b ) and (3.9a , b ) are integrated upstream to obtain the numerical results.
Open access
