Two-layer formulation for long-runout turbidity currents: theory and bypass flow case

Abstract

Turbidity currents, which are stratified, sediment-laden bottom flows in the ocean or lakes, can run out for hundreds or thousands of kilometres in submarine channels without losing their stratified structure. Here, we derive a layer-averaged, two-layer model for turbidity currents, specifically designed to capture long-runout. A number of previous models have captured runout of only tens of kilometres, beyond which thickening of the flows becomes excessive, and the models without a lateral overspill mechanism fail. In our framework, a lower layer containing nearly all the sediment is a faster, gravity-driven flow that propels an upper layer, where sediment concentration is nearly zero. The thickness of the lower layer is controlled by competition between interfacial water entrainment due to turbulent mixing and water detrainment due to sediment settling at the interface. The detrainment mechanism, first identified in experiments, is the key feature that prevents excessive thickening of the lower layer and allows long-runout. Under normal flow conditions, we obtain an exact solution to the two-layer formulation, revealing a constant velocity and a constant thickening rate in each of the two layers. Numerical simulations applied

to gradually varied flows on both constant and exponentially declining bed slopes, with boundary conditions mimicking field observations, show that the predicted lower layer thickness after 200 km flow propagation compares with observed submarine channel depths, whereas previous models overestimate this thickness three- to fourfold. This formulation opens new avenues for modelling the fluid mechanics and morphodynamics of long-runout turbidity currents in the submarine setting.

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. 2002; Schwenk et al. 2003). Other very large fans include the Congo Fan at $\sim1000$ km (Picot et al. 2016), 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 1993; Deptuck & Sylvester 2017). 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 1936; Kuenen 1938; Pirmez & Imran 2003). The channels are commonly highly sinuous (Imran et al. 1999; Schwenk et al. 2003), and build the fan itself by avulsing across the fan surface (Jobe et al. 2020).

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. (2011) 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.

Figure 1. Amazon Submarine Fan and channels: (a) the fan itself (after Mikkelsen et al. 1997); (b) a $\sim200$ km long reach of the submarine channel (the red box in figure 1 a) (from IFREMER, France).

Figure 2. (a) Down-channel long profiles of 20 canyon–fan systems (after Covault et al. 2011). (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 2003).

2. Existing models of turbidity currents

Meiburg & Kneller (2010) 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. (1985) and Parker et al. (1986) in the context of submarine canyons. These models are based on an extension of the 1-D model of Ellison & Turner (1959) for the downstream evolution of sediment-free dense underflows, such as those driven by thermohaline effects. Ellison & Turner (1959) 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. (1985) and Parker et al. (1986) adapt concepts from Pantin (1979) and Parker (1982) 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 (1994), Bonnecaze et al. (1993), Fay (2012), Hu et al. (2012, 2015), Cao et al. (2015) and Bolla Pittaluga et al. (2018).

The 3-equation and 4-equation layer-averaged models of Fukushima et al. (1985) and Parker et al. (1986) have been used to study the formation of sediment waves and submarine gullies on the seafloor (Izumi 2004), and cyclic step instability within the flow (Kostic & Parker 2006; Wu & Izumi 2022). 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. 1998; Halsey & Kumar 2019). Wahab et al. (2022) have applied the 4-equation model to the morphodynamics of submarine fans. Traer et al. (2018a , 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. 2017).

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 (1989), Sequeiros et al. (2009), Yeh et al. (2013), Luchi et al. (2018) and Iwasaki & Parker (2020). High-fidelity models have been presented by Cantero et al. (2009a , b ), Biegert et al. (2017), Salinas et al. (2019a , b , 2020, 2021a , 2022) and Xie et al. (2023b ). Wells & Dorrell (2021) 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. (2018) 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. 1998; Halsey & Kumar 2019).

The Reynolds-averaged model of Luchi et al. (2018) and the high-fidelity model of Salinas et al. (2021b ) 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.

Figure 3. 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. (1986), 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

(2.1) \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}

(2.2) \begin{equation} \frac {{\partial \delta }}{{\partial t}} + \frac {{\partial U\delta }}{{\partial x}} = {e_{ws}}U, \end{equation}

(2.3) \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. 1987)

(2.4) \begin{equation} {e_{ws}}={e_{ws}}[Ri_b]= \frac {{0.075}}{{\sqrt {1 + 718\,Ri_b^{2.4}} }}, \end{equation}

(2.5) \begin{equation} R{i_b} = \frac {{RgC\delta }}{{{U^2}}} = \frac {{Rg{q_s}}}{{{U^3}}}, \end{equation}

(2.6) \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

(2.7) \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}

(2.8) \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}

(2.9) \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:

(2.10) \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

(2.11) \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:

(2.12) \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 (1959), and in such a state, Richardson number $Ri_b$ and velocity $U$ attain constant values, and thickness $\delta$ increases linearly with distance downstream:

(2.13) \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}

(2.14) \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

(2.15) \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. (2020), 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. 2024). 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. (2018) 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. (2015) 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. (2015) 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. (2015) 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. 1998; Halsey & Kumar 2019).

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 (1959) 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. 2017), and is built into, for example, the high-fidelity calculations reported in Balachandar et al. (2024). Bonnecaze et al. (1993) 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. 2004). In the course of experiments on turbidity currents flowing into minibasins, Lamb et al. (2006) and Toniolo et al. (2006a , 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. 2024), 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. (2006b ) proposed a formulation according to which (2.2) is amended to

(2.16) \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. (2018) incorporated this formulation into the 3-equation model.

Luchi et al. (2018) 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. (2018) 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. (2018) 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. (2018) 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 (1987 a , b ), originally designed for the treatment of saline wedges. (Several misprints were corrected in Arita 1998.) 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 (1987 b ) for saline wedges to (2.4) (Parker et al. 1987), 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

(3.1) \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}

(3.2) \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}

(3.3) \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. (1986), 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

(3.4) \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}

(3.5) \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 (1987b ) 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

(3.6a) \begin{align} {e_{ws}} & = \frac {{0.075}}{{\sqrt {1 + 718\,Ri_I^{2.4}} }},\\[-10pt]\nonumber \end{align}

(3.6b) \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

(3.7a) \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}

(3.7b) \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}

(3.7c) \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:

(3.8) \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

(3.9a) \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}

(3.9b) \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

(3.10a) \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}

(3.10b) \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 (1959), 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

(3.11a) \begin{align} \frac {{{\rm d}{\delta _L}}}{{{\rm d}x}} = {A_L}, \end{align}

(3.11b) \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 (1959) 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

(3.12a–d) \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

(3.13a) \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}

(3.13b) \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. 2013; Ma et al. 2017, 2020; Simmons et al. 2020). A wide range of values of $S$ and $\tilde v_s$ , compatible with submarine channels (Covault et al. 2011), 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).

Figure 4. 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. 2004; Dorrell et al. 2014; Hughes Clarke 2016; Paull et al. 2018; Pope et al. 2022; Talling et al. 2022), 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 1952; Hsu et al. 2008). 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. 2022). 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. 2017). 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. 2020). 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. (2013).

Konsoer et al. (2013) 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 (2003). 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. (2013). 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.

Figure 5. (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. (2011) 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. 2006), 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. 2011). 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. 1985; Parker et al. 1986), and the 3-equation model modified to include detrainment in (2.1), (2.16) and (2.3) (Toniolo et al. 2006a ; Bolla Pittaluga et al. 2018). 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:

(3.14) \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.

Figure 6. 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.

Figure 7. 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 2003; Fildani et al. 2006).

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. 2006).

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 2003). 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. (2011) and later by Bolla Pittaluga et al. (2018).

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. 2017) and for up to several weeks in the distal part of the system (Baker et al. 2024). 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 (2006) and Cao et al. (2015) can model unsteady flow, whether Froude-supercritical or Froude-subcritical. It was used by Kostic & Parker (2006) to reproduce the migrating turbidity current head and hydraulic jump of one of the experiments of Garcia & Parker (1989). The same formulation could presumably be used to model the hydraulic jumps observed in the field by Sumner et al. (2013) and Clarke & Jhon (2016).

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

(4.1) \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. (1987) and Garcia & Parker (1991) 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 (1998). 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. 2022), flows that grow rapidly in sediment volume by a factor 100–1000 (Pope et al. 2022; Böttner et al. 2024), and the formative conditions for large trains of knickpoints (Heijnen et al. 2020) or smaller trains of upstream-migrating crescentic bedforms (Clarke & Jhon 2016). A morphodynamic version of the two-layer model can also be adapted to model the incision necessary to excavate submarine canyons (Zhang et al. 2017) with the aid of, for example, the sandblasting model of Lamb et al. (2008).

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. 1998; Halsey & Kumar 2019) and the emplacement of channelized submarine fans (Wahab et al. 2022). 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. 2011), 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. (1999) on meandering submarine channels with self-formed levees, in which case sinuosity might abet mixing between the two layers (Straub & Mohrig 2008).

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. 2018; Wang et al. 2020). A first step in overcoming this issue is suggested by the model of Spinewine & Capart (2013) 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. (1985) and Parker et al. (1986). Cao et al. (2015) 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. (1986) 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 (2012).

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. 1987). 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 1982; Parker et al. 1987; Johnson & Hogg 2013). 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. 2019; Marshall et al. 2021; Salinas et al. 2021a ; Lloyd et al. 2022). 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 (1987a ) 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 1976) and causing direct loss of water, e.g. flow stripping (Fildani et al. 2006; Spinewine et al. 2011). 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. (2015) 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. 2019b ). 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. 1987) 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. (2006a ). 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 2017). 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. 2000). More complex grain size distributions and patterns of dispersion can be considered in the future (Xie et al. 2023a ).

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. (2018), 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. 2006a ). 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 (1959) 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.