Turbidity currents propagating down an inclined slope: particle auto-suspension

Abstract

The turbidity current (TC), a ubiquitous fluid–particle coupled phenomenon in the natural environment and engineering, can transport over long distances on an inclined terrain due to the suspension mechanism. A large-eddy simulation and discrete element method coupled model is employed to simulate the particle-laden gravity currents over the inclined slope in order to investigate the auto-suspension mechanism from a Lagrangian perspective. The particle Reynolds number in our TC simulation is $0.01\sim 0.1$ and the slope angle is $1/20 \sim 1/5$. The influences of initial particle concentration and terrain slope on the particle flow regimes, particle movement patterns, fluid–particle interactions, energy budget and auto-suspension index are explored. The results indicate that the auto-suspension particles predominantly appear near the current head and their number increases and then decreases during the current evolution, which is positively correlated with the coherent structures around the head. When the turbidity current propagates downstream, the average particle Reynolds number of the auto-suspension particles remains basically unchanged, and is higher than that of other transported particles. The average particle Reynolds number of the transported particles exhibits a negative correlation with the Reynolds number of the current. Furthermore, the increase in particle concentration will enhance the particle velocity, which allows the turbidity current to advance faster and improves the perpendicular support, thereby increasing the turbidity current auto-suspension capacity. Increasing slope angle will result in a slightly larger front velocity, while the effect of that on the total force is insignificant.

1. Introduction

Turbidity currents (TCs) have been recognized as a widespread geophysical phenomenon in, for example, oceans, reservoirs, estuaries and lakes, which can transport large numbers of deposited particles long distances downstream (Middleton 1966; Meiburg & Kneller 2010; Wells & Dorrell 2021). With large-scale and high-intensity sediment transport, TC is inextricably linked with the rapid changes in riverbed morphology (Parker et al. 1987), the formation of submarine oil and gas (Meiburg, Radhakrishnan & Nasr-Azadani 2015) and the stability of underwater structures (Fine et al. 2005), and has attracted the attention of many researchers.

In recent years, researchers (Simpson 1982; Meiburg et al. 2015; Ouillon et al. 2021) stated that it is the horizontal pressure difference between the muddy and ambient fluid area, caused by the particle suspension, which drives the horizontal movement of the flow. The suspension regime in TC can be divided into two types: re-suspension and auto-suspension. Both of them are essentially the upward particle suspension in the bed-normal direction. Re-suspension is equivalent to erosion and refers to particles being lifted from the deposited substrate by turbulence (Meiburg & Kneller 2010; Strauss & Glinsky 2012), while auto-suspension corresponds to the event that the particles could remain suspended without additional fluid net energy expenditure while being transported (Bagnold 1962). This paper focuses on auto-suspension processes of the TCs rather than re-suspension (erosion), that is, we explore how the transported particles remain in the suspension state.

The auto-suspension mechanism of tiny particles was firstly proposed by Knapp (1938) and Bagnold (1962). Bagnold (1962) pointed out that the particles could remain suspended when the vertical component of bed slope velocity $u_p^{//}$ is greater than the vertical sinking velocity of the particles $w_p$, i.e. $u_p^{//}\sin \theta >w_p$, in which $\theta$ is the slope angle. With the auto-suspension of the TC and horizontal pressure difference, the flow will become more inclined to enter the self-acceleration mode (Pantin & Franklin 2011). Then, the TC will erode the bottom bed and entrain the deposited sediments continuously. Because of the decrease of the horizontal pressure difference, the sedimentary bed can no longer be eroded or the terrain change. At the same time, the self-acceleration behaviour or ignition behaviour (Parker 1982; Parker, Fukushima & Pantin 1986) may be restricted or return to auto-suspension mode. In essence, auto-suspension is a prerequisite and necessary condition for self-acceleration (Parker 1982).

Since the TC in the field is quite strong (Meiburg et al. 2015), it is very difficult to investigate auto-suspension by measuring the currents directly with instruments (Xu et al. 2002; Xu, Noble & Rosenfeld 2004; Liu et al. 2012). Fortunately, after TCs, visible on-site traces, such as the source of sediments and the changes in riverbed elevation, provide convincing indirect evidence for the long-distance transport (Andrieux, Cooper & Wood 2013; Li & Gong 2018), which also proves the existence and importance of the auto-suspension of the TC.

Previous laboratory and numerical studies had focused more on the hydrodynamic characteristics (such as Froude number), fluid velocity profile, particle concentration and sediment sequence to explain auto-suspension in the TC (Southard & Mackintosh 1981; Parker et al. 1987; Pantin 2001; Sequeiros et al. 2009; Sequeiros, Mosquera & Pedocchi 2018; Pantin & Franklin 2011). Pantin (2001) used the method of continuous inflow to conduct gravity flow experiments in a tubular channel on a steep slope ($\tan \theta =0.36$). He obtained that the evidence of auto-suspension is the sharp erosional upper bed surface, the increasing sediment proportion and the accelerated flow velocity. Pantin & Franklin (2011) improved the experimental technique of auto-suspension in the laboratory on the basis of Pantin (2001). In both of the above studies, the interaction between the TC and the ambient fluid was avoided by a closed circular tube, which was different from the actual TC process. In contrast, the experiments of Sequeiros et al. (2009) and Sequeiros et al. (2018) allowed TC to entrain the ambient water during the propagation. Sequeiros et al. (2009) used lightweight plastic particles to achieve auto-suspension on a gentle slope ($\tan \theta =0.05$) and pointed out that the finer bed material was easier to entrain, aiding the self-acceleration in the TC process. Sequeiros et al. (2018) further revealed the positive feedback mechanism between the bed sediment entrainment and the flow velocity. And they proposed that the supercritical flow regime was a necessary and insufficient condition for self-acceleration and the thinning and fining rates along the turbidite could be adopted as an indicator to identify the acceleration of turbidity current.

The Eulerian–Eulerian model is widely used to simulate the evolution of TCs, including the Reynolds-averaged Navier–Stokes (RANS) model, large-eddy simulation (LES) model and direct numerical simulation (DNS) model. They are all capable of giving general features of the TCs, however, RANS lacks an adequate description of the turbulent structures. The advantage of DNS is that it can provide more accurate results prior to the other two methods, but it requires a lot of computational resources. By comparison, LES gives relatively satisfactory results with less computational cost, which makes it widely adopted (Kyrousi et al. 2018; Goodarzi et al. 2020; Koohandaz et al. 2020; Pelmard, Norris & Friedrich 2018). In recent years, the coupled model of solving individual particle motion with the discrete element method (DEM) and fluid motion has been widely used to simulate two-phase flows, such as Schmeeckle (2014), Sun & Xiao (2016), Jing et al. (2016), Maurin, Chauchat & Frey (2018), Pähtz & Durán (2018a), Pähtz & Durán (2018b), Pähtz & Durán (2020) and Zhu et al. (2020). This approach does not evolve the Boussinesq approximation, different from the Eulerian–Eulerian model (He et al. 2018). For TCs, Biegert et al. (2017) coupled a two-fluid model with a particle resolved model and explored how a TC travelling along the surface of the sediment bed propagates within the bed. Wildt et al. (2021) employed an Eulerian–Lagrangian two-way coupled LES to simulate sediment plume development. Xie et al. (2022) employed the analysis of the dispersed TC particle movement and particle acceleration/deceleration to understand the evolution of the TC. However, such particle-based research is currently rare, and the effects of the complex interaction between dispersed particles and fluids in TC are still far from being understood, which is beneficial for the understanding of auto-suspension.

In order to investigate auto-suspension from the perspective of fluid–particle interactions and the movement of individual particles, this paper employs the LES-DEM model to simulate lock-exchange TCs on inclined slopes. In comparison with our earlier study on TC evolution over a flat slope (Xie et al. 2022), this present work focuses on the particle auto-suspension regimes and investigates the effects of different particle concentrations and slope angles on the auto-suspension regimes. We now additionally analyse the trajectories and spatial distribution of auto-suspension particles, the spatial characteristics of the forces acting on the auto-suspension particles and the auto-suspension index. The remainder of the paper is as follows. Section 2 introduces the governing equations of the LES-DEM model, the fluid–particle interaction forces and the numerical settings. In § 3, the fluid and particle regimes of TC, the auto-suspension mechanism by particle statistics and fluid–particle interaction regimes, the energy budget and the auto-suspension criterion are demonstrated. Summary and conclusions are drawn in § 4.

2. Methodology

In the LES-DEM approach, an Eulerian LES based on the open source code OpenFOAM is adopted for a fluid phase coarser than the particle size to predict the dynamic flow process, and the DEM based on the LIGGGHTS is adopted for the dispersed particle phase to accurately track the individual particles. These two models proceed independently, while the coupling of them at certain time intervals is implemented through the momentum exchange via CFDEMcoupling$^{\circledR}$, developed by Kloss et al. (2012). This coupled LES-DEM model has been widely validated and applied (Blais et al. 2016; Yang et al. 2018; Wang & Shen 2022), especially for TC simulations in our previous work (Xie et al. 2022). In the following, the equations for the fluid phase are described in § 2.1, and the equations for the particle phase and the fluid–particle interaction forces are given in § 2.2. Section 2.3 elaborates on the numerical conditions of TC cases performed in the present work.

2.1. Governing equations for fluid phase

The behaviour of the fluid phase is modelled by using LES, which is governed by the Navier–Stokes equations as follows (Chu et al. 2009):

(2.1)\begin{gather} \frac{\partial \alpha_f}{\partial t}+\frac{\partial \alpha _f \widetilde{u_{i,f}}}{\partial x_i} = 0, \end{gather}

(2.2)\begin{gather} \frac{\partial \alpha _f \rho _f \widetilde{u_{i,f}}}{\partial t}+ \frac{\partial \alpha _f \rho _f \widetilde{u_{i,f}} \widetilde{u_{j,f}}}{\partial x_j} ={-}\frac{\partial \tilde{p}}{\partial x_i} + \alpha_f \rho_f g_i + \frac{\partial \alpha _f \widetilde{\tau_{i j}}}{\partial x_j} - R_{i,p f} + \frac{\partial \alpha _f \widetilde{\varGamma_{i j}}}{\partial x_j}, \end{gather}

where $\alpha _f$ is the fluid volume fraction in each computational cell ($\alpha _f = 1 - \alpha _p$), $\alpha _p$ is the particle volume fraction and $\widetilde {u_{i,f}}$, $\tilde {p}$, $\widetilde {\tau _{i j}}$, $\widetilde {\varGamma _{i j}}$ are the filtered variables of fluid velocity, fluid pressure, fluid stress tensor and sub-grid stress tensor of the fluid phase, respectively; $\rho _f$ is the fluid density, and $g_i$ is the gravitational acceleration component; $R_{i, p f}$ is the momentum exchange with the particle phase, which is calculated by $R_{pf}=\sum _{\xi =1}^{k_c}\boldsymbol {F}^f_\xi /V_{cell}$ with $\boldsymbol {F}^f$ being the fluid–particle interaction forces, $k_c$ the quantity of particles contained in the corresponding fluid cell and $V_{cell}$ the volume of a computational fluid cell. The Smagorinsky model, which has been widely utilized to model two-phase Eulerian–Lagrangian flows (Schmeeckle 2014; Elghannay & Tafti 2018; Gui et al. 2018), is adopted to resolve the sub-grid-scale (SGS) stress tensor as follows:

(2.3)\begin{equation} \widetilde{\varGamma_{i j}} = \rho_f (\widetilde{u_{i,f}} \widetilde{u_{j,f}} - \widetilde{u_{i,f} u_{j,f}}) = 2 \mu_t \widetilde{S_{i j}} - \frac{1}{3} \widetilde{\varGamma _{k k}} \delta_{i j}, \end{equation}

where $\widetilde {S_{i j}}= (\partial \widetilde {u_{i,f}} / \partial x_j + \partial \widetilde {u_{j,f}} / \partial x_i)/2$, $\delta _{i j}$ is the Kronecker delta, and the SGS viscosity is obtained by

(2.4)\begin{equation} \mu_t = \rho _f \left(C_s \Delta \right) ^2 \sqrt{2 \widetilde{S_{i j}} \widetilde{S_{i j}}}, \end{equation}

where $\Delta = (V_{cell})^{1/3}$ is the characteristic length scale and $C_s=0.1$ is a constant. The effect of the constant $C_s$ on the simulation results is investigated in § 2.3.

2.2. Discrete element method

The DEM is adopted to predict the track of each particle from the Lagrangian perspective. The governing equations for the translational and rotational motions of particle $i$ are based on Newton's second law and the conservation law of angular momentum, which are calculated by (Schmeeckle 2014; Jing et al. 2016)

(2.5)\begin{gather} m_i \frac{{\rm d}\boldsymbol{u}_{p,i}}{{\rm d} t} = \boldsymbol{G}_i + \boldsymbol{F}^f_i + \sum_{j = 1}^{n^c_i}\boldsymbol{F}^c_{i j}, \end{gather}

(2.6)\begin{gather} I_i \frac{{\rm d}\boldsymbol{\omega}_{p,i}}{{\rm d} t} = \sum_{j = 1}^{n^c_i}\boldsymbol{M}^c_{i j}, \end{gather}

where $m_i$ and $I_i$ are the mass and the moment of inertia of particle $i$, respectively, $\boldsymbol {u}_{p,i}$ and $\boldsymbol {\omega }_{p,i}$ are the translational velocity and angular velocity of particle $i$, respectively, $n^c_i$ represents the number of contacting particles around particle $i$, $\boldsymbol {F}^c_{i j}$ and $\boldsymbol {M}^c_{i j}$ are, respectively, the contact force and the torque acting on the particle $i$ by particle $j$ or the boundary wall, $\boldsymbol {F}^f_i$ denotes the fluid–particle interaction force acting on the particle $i$ and $\boldsymbol {G}_i = m_i \boldsymbol {g}$ is the gravity of particle $i$. Note that fluid-induced torque is not involved in (2.6). This is because the fluid cell is always several times larger than the diameter of the particle, and the fluid variables are locally averaged over the fluid cell (Jing et al. 2016).

In the research, the fluid–particle interaction force $\boldsymbol {F}^f$ comprises the buoyancy $\boldsymbol {F}^b$, drag force $\boldsymbol {F}^d$, lift force $\boldsymbol {F}^l$ and added mass force $\boldsymbol {F}^{add}$, i.e. $\boldsymbol {F}^f=\boldsymbol {F}^b+\boldsymbol {F}^d+\boldsymbol {F}^l+\boldsymbol {F}^{add}$. The buoyancy $\boldsymbol {F}^b$ acting on a single particle with diameter $d_p$ is calculated by

(2.7)\begin{equation} \boldsymbol{F}^b={-}\tfrac{1}{6}{\rm \pi}\rho_f d_p ^3 \boldsymbol{g}, \end{equation}

with gravitational acceleration $\boldsymbol {g} = [0,0,-9.81] \,\mathrm {m}\,\mathrm {s}^{-2}$.

The drag force $\boldsymbol {F}^d$ is expressed as follows (Di Felice 1994):

(2.8)\begin{equation} \boldsymbol{F}^d=\tfrac{1}{8} C_D \rho _f {\rm \pi}d_p ^2 \left( \boldsymbol{u}_f - \boldsymbol{u}_p \right)|\boldsymbol{u}_f - \boldsymbol{u}_p| \alpha _f ^{-\chi}, \end{equation}

where $C_D$ is the drag coefficient and $\chi$ is the correction factor, which are respectively given by

(2.9)\begin{gather} C_D = \left(0.63+\frac{4.8}{\sqrt{\alpha _f {\textit{Re}}_p}} \right)^2, \end{gather}

(2.10)\begin{gather} \chi = 3.7 -0.65 \exp\left[ -\frac{\left( 1.5-\log_{10}\left(\alpha_f {\textit{Re}}_p \right) \right)^2}{2} \right], \end{gather}

where ${\textit {Re}}_p$ is the particle Reynolds number that can be expressed as

(2.11)\begin{equation} {\textit{Re}}_p = \frac{\rho_f d_p |\boldsymbol{u}_f - \boldsymbol{u}_p|}{\mu _f}. \end{equation}

The lift force $\boldsymbol {F}^l$ here follows Loth & Dorgan (2009)

(2.12)\begin{equation} \boldsymbol{F}^l=\frac{1}{8} C_L \rho _f {\rm \pi}d_p ^2 |\boldsymbol{u}_f - \boldsymbol{u}_p| \left[ \left( \boldsymbol{u}_f - \boldsymbol{u}_p \right) \times \frac{\boldsymbol{\omega}_f}{|\boldsymbol{\omega}_f|} \right], \end{equation}

where $\boldsymbol {\omega }_f$ is the fluid vorticity, and $C_L$ is the lift coefficient given by (McLaughlin 1991; Loth & Dorgan 2009)

(2.13)\begin{equation} C_L=J^{{\ast}}\frac{12.92}{\rm \pi}\sqrt{\frac{\omega^\ast}{{\textit{Re}}_p}} + \varOmega^\ast_{p,e q} C^\ast_{L,\varOmega}, \end{equation}

with $\omega ^\ast =|\boldsymbol {\omega }_f|d_p / |\boldsymbol {u}_f - \boldsymbol {u}_p|$. The function $J^\ast$ in (2.13) reads (Mei 1992)

(2.14)\begin{equation} J^\ast{=} 0.3 \left\{ 1 + \tanh \left[ \frac{5}{2} \left( \log_{10}\sqrt{\frac{\omega^\ast}{{\textit{Re}}_p}} + 0.191 \right) \right] \right\} \left\{ \frac{2}{3} + \tanh\left[ 6 \sqrt{\frac{\omega^\ast}{{\textit{Re}}_p}}\right] \right\}. \end{equation}

Furthermore, the empirical correction $C^\ast _{L,\varOmega }$ and empirical model for $\varOmega ^\ast _{p,eq}$ in (2.13) are given by (Loth & Dorgan 2009)

(2.15) \begin{gather} \varOmega^\ast_{p,e q}=\frac{\omega^\ast}{2}( 1-0.0075{\textit{Re}}_\omega )( 1-0.062\sqrt{{\textit{Re}}_p} - 0.001{\textit{Re}}_p ), \end{gather}

(2.16) \begin{gather} C^\ast_{L,\varOmega} = 1 - \{ 0.675 + 0.15( 1+\tanh[0.28 (\varOmega^\ast_p-2 ) ] ) \} \tanh[0.18 \sqrt{{\textit{Re}}_p}], \end{gather}

where ${\textit {Re}}_\omega =\rho _f |\boldsymbol {\omega }_f|d_p ^2 / \mu _f$ and $\varOmega ^\ast _p = |\boldsymbol {\omega }_p|d_p/|\boldsymbol {u}_f - \boldsymbol {u}_p|$.

The added mass force $\boldsymbol {F}^{a d d}$ is formulated by

(2.17)\begin{equation} \boldsymbol{F}^{a d d} = \frac{1}{6}C_{a d d}\rho _f {\rm \pi}d_p ^2 \left( \frac{{\rm D}\boldsymbol{u}_f}{{\rm D} t} - \frac{{\rm D}\boldsymbol{u}_p}{{\rm D} t} \right), \end{equation}

where $C_{a d d}=0.5$ is the added mass coefficient. Although the added mass force on particles in open channel sediment transport studies is generally considered to be unimportant (Schmeeckle 2014; Pähtz et al. 2021), we found in our previous study (Xie et al. 2022) that it cannot be neglected in TCs. Equation (2.17) assumes that particles are not in close proximity to one another. Due to the low particle volume fraction, this is approximately the case in our TC simulations. Other forces, such as the Basset force, are always secondary (several orders of magnitude smaller) as compared with these predominant forces, and are difficult to quantify analytically (Durán, Claudin & Andreotti 2011; Schmeeckle 2014).

In coupling two phases, the fluid–particle interactions are transmitted only by the momentum exchange term ($R_{i,pf}$ in (2.2)). The effect of particles on SGS stresses of the fluid is ignored due to low particle Reynolds number (${\textit {Re}_{p}}=0.01\sim 0.1$) and small particle–turbulence length-scale ratio, $d_p / \eta \ll 1$, where $\eta \sim O (10^{-3})$ is the Kolmogorov scale (Elghobashi & Truesdell 1992; Bagchi & Balachandar 2004). This means small vortical structures caused by individual particles are negligible. In addition, the turbulence modulates the flow by the SGS model, which then impacts the particle movement, and we do not model the effect of the sub-grid fluid field on the particles, since this effect is relatively weak when the LES velocity field is well resolved (Armenio, Piomelli & Fiorotto 1999).

Assuming that the particles are soft spheres, the particle–particle contact force $\boldsymbol {F}^c_{ij}$ can be evaluated by using an elastic spring and a viscous dashpot (Cundall & Strack 1979). The contact force caused by the collision, which includes the normal component $\boldsymbol {F}^n_{ij}$ and tangential component $\boldsymbol {F}^t_{ij}$, is written as

(2.18) \begin{align} {\boldsymbol{F}}_{ij}^c &= {\boldsymbol{F}}_{ij}^n + {\boldsymbol{F}}_{ij}^t\nonumber\\ &= \left\{ {\begin{array}{@{}cc} {( {{k^n}\delta _{ij}^n - {\gamma ^n}\nu _{ij}^n} ) + ({{k^t}\delta _{ij}^t - {\gamma ^t}\nu _{ij}^t} ),} & {{\boldsymbol{F}}_{ij}^t < {\mu _c}{\boldsymbol{F}}_{ij}^n},\\ {( {{k^n}\delta _{ij}^n - {\gamma ^n}\nu _{ij}^n} ) + {\mu _c}{\boldsymbol{F}}_{ij}^n,} & {{\boldsymbol{F}}_{ij}^t \ge {\mu _c}{\boldsymbol{F}}_{ij}^n} \end{array}}, \right. \end{align}

where $\delta ^n_{ij}$, $\delta ^t_{ij}$, $\nu ^n_{ij}$ and $\nu ^t_{ij}$ are the normal overlap distance, the tangential displacement, the normal relative velocity of particles $i$ and $j$ and the tangential relative velocity of particles $i$ and $j$, respectively. Here, $\mu _c$ ($=0.5$ adopted in this study) is the coefficient of friction, $k^n$, $k^t$, $\gamma ^n$ and $\gamma ^t$ are the elastic constants for normal contact and tangential contact and the viscoelastic damping constants for normal contact and tangential contact, respectively. These parameters are determined by a Young's modulus $Y=5\times 10^6$ Pa, Poisson ratio $\nu =0.45$ and coefficient of restitution $e=0.3$, which are the inherent properties of the particle material.

2.3. Numerical set-up

The set-up of the numerical experiment is shown in figure 1. The numerical simulation domain is a hexahedral tank filled with water with a slope. The angle between the slope and the horizontal plane is defined as $\theta$ and $\tan \theta =L_B/L_x$, in which $L_B$ is the maximum lift height of the slope and $L_x$ represents the length of the hexahedral tank. Different boundary conditions are employed in our simulation. The bottom wall and the top wall in the vertical ($z$) direction are set as a no-slip wall boundary condition and a free-slip boundary condition, respectively. The longitudinal ($x$) walls (the left and right surfaces) and transverse ($y$) walls (the front and back surfaces) have no-slip wall boundary conditions and periodic boundary conditions, respectively. In figure 1, the distance between the gate and the left wall is denoted by $L_g$, and the turbidity current on the left of the gate is filled with dispersed particles. At the beginning of the simulation, the gate is pulled away and the particles move forward during the settling process, thus forming a TC. The TC head refers to the foremost part of the current, the horizontal length of which is defined as 0.15 times the length of the TC. The definition of different head lengths ($0.10\sim 0.20$ times the length of the current) does not substantially affect the quality of the results of the head average forces (Xie et al. 2022), however, its effect on the head auto-suspension index is unknown, which is explored in § 3.4.

Figure 1. Schematic view of the initial configuration of the lock-exchange TC.

In the following, all parameters are dimensionless (except some of the parameters of the particles). We take half of the left water depth ($L_{zl}=L_z - L_B$) in the domain as the characteristic length $L_{zl}/2$ and the buoyancy velocity $u_b$ as the characteristic velocity, which is calculated as follows:

(2.19)\begin{equation} u_b=\sqrt{C_0|\boldsymbol{g}|\frac{\left(\rho_p - \rho_f \right)L_{z_l}}{2\rho_f}}, \end{equation}

where $C_0$ denotes the initial particle concentration. The dimensionless time is given by $2tu_b/L_{zl}$. The forces acting on particles are made dimensionless by the effective gravity $\boldsymbol {G}'$, which is the absolute value of the resultant force of buoyancy $\boldsymbol {F}^b$ and gravity $\boldsymbol {G}$ ($=\rho _p {\rm \pi}d_p ^3 \boldsymbol {g} / 6$), that is $\boldsymbol {G}'=\boldsymbol {G}+\boldsymbol {F}^b = (\rho _p-\rho _f){\rm \pi} d_p^3 \boldsymbol {g}/6$.

As we have a slope in the TC case shown in figure 1, we employ the variables $\psi$ in the bed-parallel and bed-normal directions ($\psi ^{//}$ and $\psi ^\bot$) for the convenience of analysis, which are calculated as follows:

(2.20)\begin{gather} \psi^{/{/}} = \psi_x\cos\theta-\psi_z\sin\theta, \end{gather}

(2.21)\begin{gather} \psi^{\bot}= \psi_x\sin\theta+\psi_z\cos\theta, \end{gather}

with a longitudinal ($x$) variable $\psi _x$ and a vertical ($z$) variable $\psi _z$.

Numerical settings of different cases are shown in table 1. In this study, the effects of different particle concentrations and different inclined slopes are of concern in cases $1,2,3$ and cases $1,4,5,6$, respectively. Fluid grid size for all coordinates is $2\sim 4$ times the particle diameter $d_p$. The dimensional particle diameter is 0.00005 m (50 $\mathrm {\mu }$m) and the particle density $\rho _p$ is $1200\, \mathrm {kg}\,\mathrm {m}^{-3}$. The density of the fluid is $1000\,\mathrm {kg}\,\mathrm {m}^{-3}$, thus the ratio of the densities $s=\rho _p/\rho _f=1.2$. The particle Reynolds number ${\textit {Re}}_p$ is relatively small (approximately $0.01\sim 0.1$), and the particle Stokes number $St=s \rho _f | \boldsymbol {u}_f - \boldsymbol {u}_p | d_p/(9\mu _t)$ is $O(10^{-3}\sim 10^{-1})$, indicating that the particle flow ought to follow the fluid streamlines of the TCs and not be determined by the inertia of the particles. The front Reynolds number ${\textit {Re}}_f$ is given by ${\textit {Re}}_f = \rho _f u_{front} L_h / \mu _t$, where $u_{front}$ is the front velocity, $L_h$ is the height of the TC head and $\mu _t$ is the fluid dynamic viscosity. The Reynolds number ${\textit {Re}}$ is defined by ${\textit {Re}}=\rho _f u_b L_{zl}/(2\mu _t)$. The initial particle concentration $C_0$ varies from 1 % to 2 %.

Table 1. Parameters of numerical simulations for each case.

The DEM time step is set to $10^{-7}$ s, which is lower than the Rayleigh critical time step proposed by Li, Xu & Thornton (2005). To improve computational efficiency, we set the fluid–particle coupling interval $N_t$ to be 100. In other words, the momentum exchange between the particle and fluid phases is integrated every 100 DEM time steps. At this time, the computational fluid dynamics (CFD) time step is equal to $10^{-5}$ s, which corresponds to a Courant number less than 0.01. Here, we briefly discuss the effect of the coupling interval $N_t$ on the simulation results. We keep the DEM time step unchanged and adjust the coupling intervals to 20, 50, 100 and 200. The comparisons of the transverse-averaged fluid velocity at four different positions ($x$ = 3.6, 4.0, 4.6 and 4.8) at $t=6$ for case 1 with four coupling intervals $N_t$ (20, 50, 100, and 200) are shown in figure 2(a). The profiles for various coupling intervals $N_t$ are very similar, indicating that the model results are not sensitive to $N_t$. It is worth mentioning that an excessively large coupling interval may cause particles to move in multiple grids within one CFD time step, which will lead to uncertainty in interphase momentum exchange.

Figure 2. Fluid velocity profiles at the four selected positions ($x$ = 3.6, 4.0, 4.6 and 4.8) at $t=6$ for case 1 (a) with four different coupling intervals ($N_t=$ 20, 50, 100 and 200), (b) with three different computational grid resolutions ($300(N_x)\times 30(N_y)\times 90(N_z)$, $250\times 25\times 80$ and $200\times 20\times 60$) and (c) with three different Smagorinsky constants ($C_s=$ 0.05, 0.10 and 0.15).

The accuracy of the LES model is affected by the quality of the grid resolution. Pelmard et al. (2018) proposed that a good resolution can be obtained if the ratio of the SGS viscosity to the molecular $N_v$ in the low turbulence region above the upper boundary of the TC is lower than 0.3 and the ratio of the SGS shear stress to the resolved Reynolds stress $N_s$ inside the turbulent mixing layer satisfies $N_s<0.05$. When using the grid resolutions in the table 1, $N_v$ is of the order of $O(10^{-4}\sim 10^{-6})$ and $N_s$ of the order of $O(10^{-4})$ in TC simulations, which meet the requirements of Pelmard et al. (2018). This demonstrates that the grid resolutions are fine for LES. Moreover, the effect of different grid resolutions on the results is investigated. The comparisons of the transverse-averaged fluid velocity at four different positions with three meshes ($300 \times 30 \times 90$, $250\times 25\times 80$ and $200\times 20\times 60$) are shown in figure 2(b). The velocity profiles essentially coincide with each other, which means the three grids give consistent results for the current and the numerical model set-up here is reliable. As shown in table 1, the grid of $250\times 25\times 80$ cells for all cases except case 4 is adopted in our simulations. For case 4, the slope angle is large and $L_z$ is 4 in table 1 and the grid is $250\times 25\times 100$.

In general, the constant $C_s$ in the Smagorinsky model needs to be adjusted for different flow events (Katopodes 2018). We discuss the sensitivity of $C_s$ here. Figure 2(c) shows the transverse-averaged fluid velocity at four positions with three constants $C_s$ (0.05, 0.10 and 0.15). As can be observed, the simulation of the TCs in this work is essentially unaffected by the changes in the constant $C_s$. A possible reason for this could be the fact that the particle Reynolds number is very low (${\textit {Re}}_p = 0.01\sim 0.1$), suppressing vortices at the particle scale (note that also the flow Reynolds number is not particularly large, $O(10))$. This indicates that the details of the Smagorinsky model are not very relevant for our cases.

3. Results and discussion

3.1. Development of the TC

The TC front $x_{front}$, an important parameter in TC, is defined as the position of the particle moving furthest in the longitudinal direction. Figure 3 plots the comparison of numerical and experimental results for the time evolution of the front position. The simulation results are consistent with the experimental data by Gladstone, Phillips & Sparks (1998), which shows the feasibility of the model for simulating TCs. It can be seen from the gradient of the curves (the front velocity $u_{front}$) in figure 3(a,b) that the front velocity reaches its maximum rapidly and then remains roughly unchanged. The dimensionless velocities increase slightly with the increase of the particle concentration in figure 3(a). The characteristic velocities $u_b$ are 0.0099 ${\rm m}\,{\rm s}^{-1}$, 0.0121 ${\rm m}\,{\rm s}^{-1}$ and 0.0140 ${\rm m}\,{\rm s}^{-1}$ for cases 1, 2 and 3, respectively. This means that the actual front velocity $u_{front}$ increases positively in relation to the increase in initial particle concentration, which is consistent with the understanding of previous studies (Bonnecaze, Huppert & Lister 1993; Hu et al. 2020). In figure 3(b), as the slope angle increases ($\theta <20^\circ$), the gradient becomes larger, indicating a larger front velocity of the TC (He et al. 2018; Steenhauer, Tokyay & Constantinescu 2017). The steepening of the slope raises the current barycentre, resulting in an increasing available particle potential energy, which will be converted into the kinetic energy of the system.

Figure 3. Temporal evolution of simulated TC fronts under different cases, and the comparison of simulated fronts with the measured data (Gladstone et al. 1998).

In previous studies, the fluid velocity profile of the TC is generally similar (Altinakar, Graf & Hopfinger 1996). The TC area can be divided into two regions, the wall region and the jet region, which is based on the location of the maximum fluid velocity parallel to the slope $u_f^{//,max}$. The main control factors of the two regions are different: the wall region (lower layer) is dominated by the bottom wall friction (Nourmohammadi, Afshin & Firoozabadi 2011), whereas the jet region (upper layer) is dominated by the ambient fluid shear. The fluid velocity distribution of the jet region and the wall region can be described as fitted functions (Altinakar et al. 1996; Nourmohammadi et al. 2011)

(3.1) \begin{equation} \frac{u_f ^{/{/}} \left( Z \right)}{u_f^{/{/}, max}} =\left\{ \begin{array}{@{}ll} \exp \left[ - \beta\left( \dfrac{Z-H_m}{H-H_m} \right)^\gamma \right], & \text{the jet region} ,\\ \left( \dfrac{Z}{H_m} \right)^{1/n}, & \text{the wall region} \end{array}\right., \end{equation}

where $u_f^{//} (Z)$ is the fluid bed-parallel velocity at a distance of $Z$ from the slope, $H_m$ is the distance between the location of $u_f^{//,max}$ and the slope and $H$ denotes the depth-averaged current height which is expressed as follows:

(3.2) \begin{equation} H=\frac{\displaystyle\left(\int_0^{Z_{top}} u_f^{/{/}} \,{\rm d} Z \right)^2}{\displaystyle\int_0^{Z_{top}} ( u_f^{/{/}})^2\,{\rm d} Z}, \end{equation}

where $Z_{top}$ is the farthest distance from the slope for which the fluid bed-parallel velocity is positive. Here, $\beta$, $\gamma$ and $n$ are empirical coefficients. We now proceed to employ (3.1) to fitted fluid velocity profiles at four different locations ($x$ = 3.0, 3.6, 4.2 and 4.8) at $t=10$. Table 2 shows the tuned parameters and the fitting results for all cases. The simulated fluid velocity profiles can be well fitted with (3.1) (the determination coefficients ${R}^2$ are all higher than 0.95). It proves the reasonable performance of the model for reproducing TC processes.

Table 2. Comparison of fitting results for present six cases.

The local average spatial variations of the bed-parallel ($u^{//}_{p,i}$) and bed-normal ($u^\bot _{p,i}$) velocities at $t=2$ and $t=6$ in case 1 are shown in figure 4. In figure 4(a,c), the current is mainly transported downstream along the slope in the TC except for the particles in the top layers at the beginning stage (Nourmohammadi et al. 2011; Wildt et al. 2020; Hitomi et al. 2021). The reason why top particles are transported upstream (blue coloured) is the intrusion flow. Then all particles move downstream at $t=6$. For the bed-normal velocities shown in figure 4(b,d), most particles are gradually settling to the bottom wall due to the dominance of gravity, while $u^\bot _{p,i}$ of the particles in the TC head border region is positive or zero on average. This exhibits that they cannot only stay suspended, but can even rise beyond the suspension state (defining them as auto-suspension particles or highly auto-suspension particles).

Figure 4. Spatial variations of (a,c) local average particle bed-parallel velocity and (b,d) local average particle bed-normal velocity at two selected times. Black line indicates the TC profile.

Given that particle auto-suspension occurs mostly in the head, we display the profiles of particle bed-parallel and bed-normal velocities at 0.05 current lengths upstream from the front at $t=6$ under different cases, as shown in figure 5. For different particle concentrations at the same slope angle, the front positions shown in figure 3 are quite similar to cases 1–3. As shown in figure 5(a), with increasing particle concentration, the particle bed-parallel velocity increases, which agrees well with the previous study (Hu et al. 2020). The different directions of bed-normal velocity of particles in the upper and lower layers of the head exhibit the separation of particles in the TC head. In figure 5(b), it is observed that most particles in the lower layer move toward the slope while those in the upper layer ($z'>0.7$) move away from the slope, which can be directly viewed in the spatial variation of case 1 in figure 4(b,d). In the lower layers, the fluids are subjected to wall viscous shear, resulting in a low fluid velocity, which is insufficient to achieve particle suspension, while the velocity of the fluids in the upper layers is strong enough to help achieve suspension. With the increase of the slope angle shown in figure 5(c), the maximum particle velocity along the slope increases, and the height of the TC grows, which is really different from figure 5(a). It manifests itself as an increase in the ability of auto-suspension of the TC head. Note that the increasing slope angle thickens the lower region (figure 5d) where the particle bed-normal velocity is negative, and it simultaneously inhibits the upward movement of the upper particles and the settling of the lower particles, implying the decrease in the separation rate of the upper and lower particles.

Figure 5. (a,c) Particle bed-parallel velocity profiles and (b,d) particle bed-normal velocity profiles at 0.05 current lengths upstream from the front at $t=6$ under different cases. Here, $z'=z-z_b$ with $z_b$ denoting the bed elevation; ${\textit {Re}}_f$ (case 1) = 27, ${\textit {Re}}_f$ (case 2) = 35, ${\textit {Re}}_f$ (case 3) = 46, ${\textit {Re}}_f$ (case 4) = 37, ${\textit {Re}}_f$ (case 5) = 24, and ${\textit {Re}}_f$ (case 6) = 23.

During the evolution of TCs, typical streamwise and transverse coherent vortical structures are exhibited (Dai & Huang 2022), the qualitative cognition of which is independent of the Reynolds number (Nasr-Azadani & Meiburg 2014; Nasr-Azadani, Meiburg & Kneller 2018). The coherent structure can also reflect the lobe-and-cleft structures of the TC (Espath et al. 2015; Francisco, Espath & Silvestrini 2017), which are related to the mixing dynamics. Here, we attempt to utilize the coherent vortical structures with the $Q$-criterion to survey the suspension regime of the transported particles in TC. The quantity $Q$ is given by (Hunt, Wray & Moin 1988)

(3.3)\begin{equation} Q = \frac{1}{2}\left( {{\varOmega _{ij}}{\varOmega _{ij}} - {S_{ij}}{S_{ij}}} \right), \end{equation}

where the rotation rate tensor $\varOmega _{ij}$ and strain rate tensor $S_{ij}$ are respectively expressed as follows:

(3.4)\begin{gather} {\varOmega _{ij}} = \frac{({u_{f})_{i,j}} - ({u_{f})_{j,i}}}{2}, \end{gather}

(3.5)\begin{gather} {S_{ij}} = \frac{({u_{f})_{i,j}} + ({u_{f})_{j,i}}}{2}. \end{gather}

Figure 6 depicts the coherent vortex structures with $Q=0.125$ at five selected moments in case 1, with the dispersed particles coloured with particle bed-normal velocity. Blue particles represent particles that settle toward the wall and red particles represent suspension particles away from the wall. The settling particles move towards the wall and along the slope due to the advantage of gravity, and the suspension particles are almost in the upper layers of the current head in figure 6.

Figure 6. Coherent vortex structures captured by $Q=0.125$ of case 1 at five selected moments: (a) $t=2$, (b) $t=4$, (c) $t=6$, (d) $t=8$ and (e) $t=10$. The dispersed particles are also drawn here and rendered with the bed-normal particle velocity.

A structurally complete and large-sized structure emerges near the upper interface of the TC in the early stages in figure 6(a,b), then gradually divides into two parts following the time evolution of the current in figure 6(c,d). The front vortex near the head moves forward with the front current during the whole simulation. The vortex at the back stays in place, separates from the mixing layer, shrinks and then vanishes eventually. The time evolution of the vortex structures is similar to the TC over a flat bed in our previous study (Xie et al. 2022). Notably, most of the auto-suspension particles are within the large-scale coherent vortical structure of the head traversing the simulation domain, where $Q>0$ suggests the dominance of the fluid rotation. Thus, it is reasonable to speculate that the uplift of these particles is closely linked to the strong counterclockwise flow around the head, which is consistent with the previous studies (Kyrousi et al. 2018; Lee 2019). The fragmented vortical structures are found in the near-wall region and the mixing layer inside the front vortex. These vortical structures are a manifestation of the combined effects of a complex fluid dynamics with the no-slip condition on the bed (Koohandaz et al. 2020) and the interaction of particle groups with fluids. In the near-wall region, the height of this vortex structure is generally limited to approximately 0.5 dimensionless lengths due to the particle transport along the way in the middle layers, and is limited by the upper interface of the TC. To facilitate the observation of the coherent structure, top-view coherent vortex structures without particles in the near-wall region below the slice $z'=0.5$ are plotted in figure 7. Relatively large-scale vortical structures can be seen in the near-wall region near the head, which are induced by the intense dynamics of the fluid and particles. During the early stage of the downstream propagation ($t=0\sim 6$), small vortex structures of the current (mainly at $x=0\sim 2$) gradually increase. The main reason is that the settling particle groups release their gravitational potential energy to enhance the local fluid kinetic energy, which facilitates the particles advancing. At $t>6$ (figure 7c–f), the small vortex structures in the near-wall region behind the current gradually decline and disappear as a result of dissipation. In addition, the streamwise development of the coherent vortical structures here is limited due to the small ${\textit {Re}}$, and the flow is mainly engaged by the transverse vortical structures, agreeing with the understandings of Nasr-Azadani & Meiburg (2014) and Koohandaz et al. (2020).

Figure 7. Top-view visualization of the coherent vortex structures below the slice $z'=0.5$ illustrated by the $Q$-criterion with an isovalue $Q=0.25$; (a) $t=2$, (b) $t=4$, (c) $t=6$, (d) $t=8$, (e) $t=10$ and ( f) $t=12$.

The inflectional instabilities akin to the inviscid Kelvin–Helmholtz mechanism (Ooi, Constantinescu & Weber 2007; Pelmard, Norris & Friedrich 2020) are not observed in figure 6. This is because the Richardson number in the present work is $Ri=g'H/U^2\sim O(10)$, where $g'=(\rho _m-\rho _f)|\boldsymbol {g}|/\rho _f$ is the effective gravitational acceleration with $\rho _m$ denoting the particle–water mixture density, and $U$ the fluid layer-averaged bed-parallel velocity of the TC. The value is much larger than the critical $Ri_c$, which is of the order of unity (Turner 1986) even for changes in slope angle and the presence of particles (Khavasi & Firoozabadi 2019; Darabian et al. 2021). Furthermore, due to the low Reynolds number of the TCs in this paper, the viscosity of the fluid is not negligible, which in turn suppresses the occurrence of the inflectional instabilities (Khavasi & Firoozabadi 2019).

3.2. Auto-suspension particle statistics and fluid–particle interactions

To understand the kinematics of auto-suspension particles in the evolution of TCs, discrete particle motion is discussed in the following section. Figure 8(a,b) directly visualizes the trajectories of 160 particles of case 1 in the tank reference frame and in the reference frame with the moving head, respectively. In these two panels, the six auto-suspension particles selected are highlighted by coloured lines. The whole process of auto-suspension mainly occurs in the head ($x-x_{front}>-2$) as shown in figure 8(b). The trajectories of the six auto-suspension particles can be described as below:

1. (i) Particles settle and transport downstream with a constant distance from the slope, as shown in figure 8(a), and they go to the front in figure 8(b).

2. (ii) Particles start moving away from the slope (figure 8a) and the front (figure 8b) during the transport downstream and reach peak points away from the slope.

3. (iii) Particles come to the end of their auto-suspension stage, settle and finally enter the lower layers of the current.

Figure 8. Trajectory of 160 particles (a) in the tank reference system and (b) in the reference frame moving with the TC head. The trajectories of the six representative particles which have entered auto-suspension state are distinguished by coloured lines.

Figure 9 plots the time evolution of the auto-suspension particle quantity for all cases. The auto-suspension particle number in all cases increases at $t<4$ and reaches the maximum around $t=4$. Then the value decreases rapidly before $t\approx 6$ and diminishes slowly until the end of the simulation, except for case 4 with a relatively large slope angle $\tan \theta =1/5$, where it increases slightly in the range $t=6\sim 8$. The possible reason is tightly related to the front vortex, which can carry the particles away from the slope at current head. In case 1, the front vortex shown in figure 6 is enhanced before $t=4$, then divides into two vortices at $t=6$ and diminishes gradually after that.

Figure 9. Time evolution of the auto-suspension particle quantity. Here, $N_p$ represents the number of all particles, and $N_{p,Auto}$ denotes the number of auto-suspension particles.

The statistics of particle Reynolds numbers ${\textit {Re}}_p$ are employed to understand the motion patterns of auto-suspension particles in TCs. Substituting the definition of the particle Reynolds number (2.11) into the definition of the particle Stokes number gives $St=\rho _p{\textit {Re}}_p/(9\rho _f)$, which means that ${\textit {Re}}_p$ and $St$ are proportional. We compute the average particle Reynolds number $\overline {{\textit {Re}}_p}$ for the auto-suspension particles and other non-auto-suspension particles in transport, the time evolutions of which are illustrated in figure 10. It is obviously shown that the average auto-suspension particle Reynolds number is always larger than the non-auto-suspension one, which implies that the auto-suspension particles have a larger $|\boldsymbol {u}_f-\boldsymbol {u}_p|$ and suffer a larger drag force. The average particle Reynolds number of auto-suspension particles decreases during the beginning stage at $t=0\sim 2$, and then remains approximately constant (with a slight increase). This indicates that the auto-suspension particles, on average, can maintain a fairly stable motion state during the downstream propagation of the TC. For the non-auto-suspension particles in transport, $\overline {{\textit {Re}}_p}$ remains nearly constant at $t=0\sim 2$, and subsequently decreases. In figure 10(a), as particle concentration increases (${\textit {Re}}$ increases from case 1 to case 3), $\overline {{\textit {Re}}_p}$ becomes smaller, which means the slip velocity between the fluid and particles decreases. This phenomenon is also observed in Sun & Liu (2022). Figure 10(b) shows that the differences of $\overline {{\textit {Re}}_p}$ are quite small, which means the changes in slope (cases 4,1,5,6, with same ${\textit {Re}}$) do not affect $\overline {{\textit {Re}}_p}$.

Figure 10. Time evolution of the average particle Reynolds number $\overline {{\textit {Re}}_p}$ for different cases. The abbreviation ‘AS’ represents ‘auto-suspension’.

The total force $\boldsymbol {F}^T$ acting on the particle in the TC comprises the fluid–particle interaction force $\boldsymbol {F}^f$, the gravity force $\boldsymbol {G}$ and particle contact force $\boldsymbol {F}^c$ (the contact force here is the sum of all contact forces experienced by the particle), and is calculated by

(3.6)\begin{equation} \boldsymbol{F}^T = \boldsymbol{F}^f+\boldsymbol{G}+\boldsymbol{F}^c = ( \boldsymbol{F}^b + \boldsymbol{F}^d + \boldsymbol{F}^l + \boldsymbol{F}^{add})+\boldsymbol{G} + \boldsymbol{F}^c. \end{equation}

Bed-normal total force $\boldsymbol {F}^T_\bot$ is employed to explain the mechanism of the auto-suspension particles. The spatial distributions of $\boldsymbol {F}^T_\bot$ of the auto-suspension particles are presented in figure 11. In the early stage at $t=2.0$, the spatial distribution of particles with a positive total force highly coincides with that with a negative total force. As the current evolves, the auto-suspension region expands rapidly and the positive force region separates from the negative one. The particles with positive total force are more widely dispersed in the lower layers, whereas those with negative total force tend to be distributed in the upper layers and closer to the TC profile. It is not difficult to conclude that, when the auto-suspension particles rise and enter the upper zone, they are more likely to reach a turning point and then exit the auto-suspension stage, as shown in figure 8. As the annular flow weakens at $t>4.0$, the spatial distribution of particles gradually narrows.

Figure 11. Spatial distribution of auto-suspension particles at six selected times. Green dot indicates the negative $\boldsymbol {F}^T_\bot$ particle, while red dot indicates the positive $\boldsymbol {F}^T_\bot$ particle; (a) $t = 2.0$, (b) $t = 3.2$, (c) $t = 4.0$, (d) $t = 6.0$, (e) $t = 8.0$, ( f) $t = 10.0$.

Analysing the TC auto-suspension process by describing the spatial variations of each predominant force is both intriguing and practical. The forces in this work are all non-dimensionalized by dividing with $|\boldsymbol {G}'|$.

Figure 12 shows four different dimensionless bed-normal force components of auto-suspension particles at the dimensionless time 4.0. The effective drag force $\boldsymbol {F}^{Ed}_\bot$, lift force $\boldsymbol {F}^l_\bot$, added mass force $\boldsymbol {F}^{add}_\bot$ and contact force $\boldsymbol {F}^c_\bot$ are taken into consideration. The bed-normal effective gravity $\boldsymbol {G}'$ is the main cause of particle settling, and the drag force $\boldsymbol {F}^d_\bot$ is the most important resistance to settling ($\boldsymbol {F}^d_\bot /|\boldsymbol {G}'|\backsim \textit {O}(1)$ and $|\boldsymbol {G}'\cos \theta |/|\boldsymbol {G}'|\backsim \textit {O}(1)$). Compared with these two, however, the magnitude of the other forces is relatively small. Thus, we consider the resultant force as the effective drag force $\boldsymbol {F}^{Ed}_\bot =\boldsymbol {F}^d_\bot +\boldsymbol {G}'\cos \theta$ and the magnitude is comparable to the other force components.

Figure 12. Four dimensionless bed-normal force components ($\boldsymbol {F}^{Ed}_\bot$, $\boldsymbol {F}^l_\bot$, $\boldsymbol {F}^{add}_\bot$ and $\boldsymbol {F}^c_\bot$) of auto-suspension particles at $t=4.0$. (a,c,e,g) The negative $\boldsymbol {F}^T_\bot$ particles and (b,d, f,h) the positive $\boldsymbol {F}^T_\bot$ particles.

Comparing negative and positive total force particles in figure 12, whether it is the lift force or added mass force, their behaviours on negative and positive $\boldsymbol {F}^T_\bot$ particles are quite similar. The positive lift force of auto-suspension particles near the TC profile is induced by the positive vorticity near the profile and the positive slip velocity perpendicular to slope. Compared with such a lift force, the added mass force is mostly considerably smaller, except when the particle is very close to the front where the added mass force exhibits an overwhelming positive value. This is due to the strong upward flow at the front. However, the effective drag force and contact force show significant differences. The effective drag forces of negative $\boldsymbol {F}^T_\bot$ particles are mostly negative, whereas those of positive $\boldsymbol {F}^T_\bot$ particles are negative near the TC profile area of the head and positive near the central area of the head. The contact force in figure 12(h) is almost larger than that in figure 12(g) in the perpendicular direction away from the slope. Collision processes can be approximately regarded as the reason for the positive contact force on positive $\boldsymbol {F}^T_\bot$ particles and the negative contact force on negative $\boldsymbol {F}^T_\bot$ particles, on average. As a consequence, the bed-normal total force is largely dominated by the effective drag force and contact force.

Combined with the movement trajectory of particles in figure 8 and figure 6, the fluid vortex at the current head, indicating the upward flow, drives the particles upward away from the slope and the particles enter the auto-suspension state. As the particles continue to rise, the fluid–particle interaction force provided by the flow and the particle–particle contact force cannot resist the particle gravity, leading to a negative $\boldsymbol {F}^T_\bot$ in figure 12. It prevents the particles from rising further, and the particles eventually reach their peak points. Henceforth, the particles withdraw from this auto-suspension event. Subsequently, particles gradually settle and enter the lower layers of the current. In essence, the long-term and long-distance cyclic work of the motion system, in which particles continuously enter and exit the auto-suspension state, allows TC to maintain the auto-suspension state from a macro-viewpoint for long distances.

Understanding the average forces on the transported head particles helps us fully explore the auto-suspension mechanism. The transported head particles can be divided into two types: deposited particles and transported particles. The former are very close (within $1.5d_p$) to the bottom wall or deposited particles and their velocity is very small ($|\boldsymbol {u}_p|<10^{-4}\,{\rm m}\,{\rm s}^{-1}$), and the latter are the remaining particles. The head average force is obtained by summing the forces acting on transported (the latter type) particles in the head and dividing by the number of corresponding particles.

Figure 13 shows the temporal evolution of average bed-parallel and bed-normal forces on transported head particles for case 1, including the total force $\boldsymbol {F}^T$, effective drag force $\boldsymbol {F}^{Ed}$, lift force $\boldsymbol {F}^l$, added mass force $\boldsymbol {F}^{add}$ and contact force $\boldsymbol {F}^{c}$. In the two directions, given that the drag force components $\boldsymbol {F}^{d}_{//,\bot }$ and effective gravity components $\boldsymbol {G}'_{//,\bot }$ are overwhelming and far greater than the other forces, we employ the effective drag force $\boldsymbol {F}^{Ed}_{//,\bot }$. For the accelerating particles in figure 13(a), the positive average $\boldsymbol {F}^{T}_{//}$ first increases and then decreases at $t=0\sim 4$, due to the sum of the lift force and contact force exceeding the sum of the negative effective drag force and added mass force. After $t=4$, $\boldsymbol {F}^{T}_{//}$ is basically near zero, suggesting that the head particles are roughly advancing at a constant speed on average, consistent with the evolution of the front position (Blanchette et al. 2005; He et al. 2018). Due to the fundamental driving action of $\boldsymbol {G}'_{//}$, the particle velocity is always greater than the fluid velocity, resulting in a negative $\boldsymbol {F}^{d}_{//}$. The negative $\boldsymbol {F}^{Ed}_{//}$ in figure 13(a) means that the negative drag force is greater than the effective gravity along the slope. The average $\boldsymbol {F}^{l}_{//}$ of the head particles is always positive, which is caused by the positive vorticity around the head (figure 6) and the positive slip velocity perpendicular to the slope. Moreover, the head particles are impacted by the collision of particles behind the head, whereby they experience a positive contact force on average. Slightly negative $\boldsymbol {F}^{add}_{//}$ means the downward particle acceleration is slightly larger than the fluid acceleration along the slope.

Figure 13. Temporal evolution of (a) average bed-parallel forces and (b) average bed-normal forces on transported head particles.

In the bed-normal direction in figure 13(b), the effective drag force $\boldsymbol {F}^{Ed}_{\bot }$ is negative when the stationary particles start to settle down at the beginning, which results in the negative total force $\boldsymbol {F}^{T}_{\bot }$. In the mid-term ($t=1.1\sim 5$), the head particles receive a positive total force $\boldsymbol {F}^{T}_{\bot }$ on average. This positive $\boldsymbol {F}^{T}_{\bot }$ at $t=1.1$ is generated due to the sum of the lift force and added mass force being larger than the sum of the effective drag force and contact force. Subsequently, the positive $\boldsymbol {F}^{T}_{\bot }$ experiences an increase and then a decrease to zero, mainly prompted by the change in the positive contact force $\boldsymbol {F}^{c}_{\bot }$. It is worth mentioning that the positive $\boldsymbol {F}^{add}_{\bot }$ induced by the upward flow around the head plays a comparable role in the bed-normal direction, as can be seen in figure 12 as well.

Figure 14 shows the temporal average variation of the average total force $\boldsymbol {F}^{T}_{//}$, average effective drag force $\boldsymbol {F}^{Ed}_{//}$, average lift force $\boldsymbol {F}^{l}_{//}$ and average contact force $\boldsymbol {F}^{c}_{//}$ acting on the transported head particles in our cases. For the same slope, one can easily observe that all the absolute dimensionless force values ($\boldsymbol {F}^{T}_{//}$, $\boldsymbol {F}^{Ed}_{//}$, $\boldsymbol {F}^{l}_{//}$ and $\boldsymbol {F}^{c}_{//}$) increase with increasing particle concentration in figure 14(a–d). The increasing $\boldsymbol {F}^{T}_{//}$ explains that the higher the concentration, the faster the head of the current (Hu et al. 2020). In figure 14(e–h), increasing the slope has a low impact on these forces parallel to the slope. Only in the later stage of the duration do the negative effective drag force and the positive lift force increase.

Figure 14. Temporal variation of bed-parallel dimensionless components of (a,e) average total force, (b, f) average effective drag force, (c,g) average lift force and (d,h) average contact force under different cases.

For the three increasing force components in all cases, reasons can be explained as follows:

1. (i) The absolute $\boldsymbol {F}^{Ed}_{//}$ increases due to the increasing particle velocity around the head, which is shown in figure 3.

2. (ii) The positive $\boldsymbol {F}^{l}_{//}$ increases due to the increase in positive vorticity around the head.

3. (iii) The positive $\boldsymbol {F}^{c}_{//}$ increases due to the increasing collision probability of the particles in the current with higher particle concentration.

Figure 15 plots average total force $\boldsymbol {F}^{T}_{\bot }$, average effective drag force $\boldsymbol {F}^{Ed}_{\bot }$, average lift force $\boldsymbol {F}^{l}_{\bot }$, average added mass force $\boldsymbol {F}^{add}_{\bot }$ and average contact force $\boldsymbol {F}^{c}_{\bot }$ on transported head particles over time for the different cases. In overall terms, an increase in particle concentration allows an increase in the absolute value of each force perpendicular to the slope, where an increase in the positive total force $\boldsymbol {F}^{T}_{\bot }$ implies an enhancement in the auto-suspension capacity. This leads to an increase in the collision probability of the particles, which in turn drives a corresponding growth in the contact force $\boldsymbol {F}^{c}_{\bot }$. Here, $\boldsymbol {F}^{c}_{\bot }$ and $\boldsymbol {F}^{add}_{\bot }$ are increased due to the enhanced flow strength. The increase of negative $\boldsymbol {F}^{Ed}_{\bot }$ reflects a decrease of the positive drag force $\boldsymbol {F}^{d}_{\bot }$, which is of interest. The greater $\boldsymbol {F}^{l}_{\bot }$ (figure 15c), $\boldsymbol {F}^{add}_{\bot }$ (figure 15d) and $\boldsymbol {F}^{c}_{\bot }$ (figure 15e) provide a stronger effect for particles to suspend and rise up, prompting an increase in positive $u_p^\bot$, which thus results in a decrease in $\boldsymbol {F}^{d}_{\bot }$. The change in slope has a negligible effect on the total force $\boldsymbol {F}^{T}_{\bot }$ (figure 15f), added mass force $\boldsymbol {F}^{add}_{\bot }$ (figure 15i) and contact force $\boldsymbol {F}^{c}_{\bot }$ (figure 15j). However, it allows the negative effective drag force and the positive lift force to be enhanced (figure 15g,h), and the magnitudes of these two changes appear to balance out.

Figure 15. Temporal variation of bed-normal dimensionless components of (a, f) average total force, (b,g) average effective drag force, (c,h) average lift force, (d,i) average added mass force and (e,j) average contact force under different cases.

3.3. Energy budget

The energy of the TC includes the particle gravitational potential energy $E_p^p$, the particle kinetic energy $E_k^p$, the fluid potential energy $E_p^f$, the fluid kinetic energy $E_k^f$ and the dissipated energy $E_{Diss}$ (Xie et al. 2022; Zhu et al. 2022). We herein concentrate on how the changes in potential and kinetic energies ($\Delta E_p^p$, $\Delta E_k^p$, $\Delta E_p^f$ and $\Delta E_k^f$) evolve throughout the history of the TC, where $\Delta$ represents the change between the energy at the current moment and the energy at the initial moment (e.g. $\Delta E_k^p(t) = E_k^p(t) - E_k^p(0)$, where $0$ is the initial moment). The four kinetic and potential energy components can be calculated as follows:

(3.7)\begin{gather} E_p^p\left(t\right) = \sum_{i=1}^{N_p} m_i \left|\boldsymbol{g}\right| z_{p,i}, \end{gather}

(3.8)\begin{gather} E_k^p\left(t\right) = \sum_{i=1}^{N_p} \left( \frac{1}{2} m_i \left|\boldsymbol{u}_{p,i}\right|^2 + \frac{1}{2} I \left|\boldsymbol{\omega}_{p,i}\right|^2 \right), \end{gather}

(3.9)\begin{gather} E_p^f\left(t\right) = \int _\varOmega \alpha_f \rho_f \left|\boldsymbol{g}\right| z \,{\rm d} V, \end{gather}

(3.10)\begin{gather} E_k^f\left(t\right) = \int _\varOmega \frac{1}{2} \alpha_f \rho_f \left|\boldsymbol{u}_f\right|^2 \,{\rm d} V, \end{gather}

where $z_{p,i}$ is the elevation of particle $i$ and $\varOmega$ represents the whole simulation domain.

Figure 16 plots the temporal variations of the energy components, non-dimensionalized by the relative initial particle gravitational potential energy depending on the reference plane ($z=L_B$), which is defined as $E_p^p(0) = \sum _{i=1}^{N_p} m_i |\boldsymbol {g}| (z_{p,i}-L_B)$. All energy components with superscript ‘$\ast$’ in the figure are non-dimensionalized. From the perspective of the energy budget, the energy of the TC on the inclined slope is converted from the particle gravitational potential energy into fluid potential energy, fluid kinetic energy, particle kinetic energy and the dissipated energy due to viscosity (Dai 2015; He et al. 2018). Defining the available potential energy as $E_{p,avail}=-\Delta E_p^p-\Delta E_p^f$, then at $t=10$, approximately $10\sim 20\,\%$ is converted into fluid and particle kinetic energy, and the rest (approximately $80\sim 90\,\%$) is dissipated due to the fluid viscosity and the particle–particle collisions. Figure 16(a–c) shows that, as the particle concentration increases, the non-dimensionalized energy components are similar. The increase in particle concentration boosts the conversion of particle kinetic energy, which means a faster particle velocity. This reflects the fact that, the larger the initial particle concentration is, the faster the TC advances (Farizan et al. 2019; Hu et al. 2020). The possible reason is that, with the increasing particle concentration, all the absolute values of the force components increase (shown in figures 14 and 15) due to the higher collision frequency.

Figure 16. Temporal evolution of non-dimensionalized energy components: (a,d) $\Delta E_p^{p,\ast }$, (b,e) $\Delta E_p^{f,\ast }$, and (c, f) $\Delta E_k^{f,\ast }$, $\Delta E_k^{p,\ast }$.

For different slope angles in cases 1,4,5 and 6, the reference plane is $z=L_B$, which means that the initial particle gravitational potential energies $E_p^p(0)$ are nearly equal. As shown in figure 16(d–f), the increase in slope leads to a significant consumption of the particle gravitational potential energy (Steenhauer et al. 2017), which leads to a considerable increase in the fluid potential energy, fluid kinetic energy and particle kinetic energy (Francisco et al. 2017; He et al. 2018).

3.4. Criterion for auto-suspension

The auto-suspension criterion is generally used for predicting the tendency of auto-suspension in TC on a slope. Based on the balance of force and energy, Bagnold (1962), Pantin (1979) and Parker (1982) proposed different criteria for the auto-suspension index of the TC head as follows:

(3.11) \begin{equation} k_{Auto,h}={-}\frac{w_{p,h}\cos\theta}{u_{p,h}^{/{/}}\sin\theta}\leq k_T, k_T=\left\{ \begin{array}{@{}ll} \cos\theta & \text{(Bagnold, 1962)},\\ 0.01 & \text{(Pantin, 1979)},\\ 1 & \text{(Parker, 1982)}, \end{array}\right. \end{equation}

where $k_{Auto,h}$ is the head auto-suspension index, $w_{p,h}$ is the average vertical velocity of head particles, $u_{p,h}^{//}$ is the average velocity along the slope for the head particles and $k_T$ is the threshold for auto-suspension criterion.

The temporal evolutions of the head auto-suspension index under different initial particle concentrations are plotted in figure 17(a), and those for different bottom slopes are plotted in figure 17(b). The criteria proposed by Bagnold (1962) and Parker (1982) differ little for the cases in this work (for case 4, $\cos \theta =0.98$), and thus the lines in figure 17 indicate criteria by Parker (1982) and Pantin (1979). Auto-suspension regions are located below the lines while non-auto-suspension regions are above the lines. The auto-suspension index is negative (smaller than the Pantin (1979) criterion) only at approximately $1.8\sim 4$ dimensionless time except for case 4 in figure 17. Actually, the negative auto-suspension index lasts only for a short period of time and auto-suspension particles do exist in case 4. The comparison between the present data and the criteria suggests that the Pantin's standard might be excessively strict for auto-suspension (Sequeiros et al. 2009). Thus we use the Parker (1982) criterion hereafter.

Figure 17. Temporal evolutions of head auto-suspension index for (a) different initial particle concentrations (cases 1, 2 and 3) and for (b) different slope angles (cases 4, 1, 5 and 6).

It can be seen in figure 17 that the exact moment when the current starts to fulfil Parker's auto-suspension criterion is roughly the same (approximately $t=1.4$). The time head particles exiting the auto-suspension will be delayed with increasing particle concentration $C_0$ or increasing slope angle $\theta$. This demonstrates that both can extend the period of auto-suspension, exhibiting the behaviour of enhancing the auto-suspension maintenance ability by increasing $C_0$ or $\theta$.

The impact of different head length definitions on the head auto-suspension index are discussed. Figure 18 depicts the time evolution of the head auto-suspension index for six cases under three head length definitions (0.10, 0.15 and 0.20 times the total length of the current). As can be observed, different head length definitions in the range of $0.10\sim 0.20$ times the TC length do not qualitatively affect the evolution of the head auto-suspension index. However, the head auto-suspension index $k_{Auto,h}$ gets smaller as the head length decreases, extending the duration during which the auto-suspension standard is satisfied. This is because particle auto-suspension mostly occurs near the front.

Figure 18. Temporal evolutions of head auto-suspension index for various cases with three head length definition methods. The dotted lines with circles represent 0.10 times the total length of the current, the solid lines 0.15 times the total length and the dotted lines with crosses 0.20 times the total length: (a) $C_0 = 0.01, \tan \theta = 1/10$; (b) $C_0 = 0.01, \tan \theta = 1/5$; (c) $C_0 = 0.015, \tan \theta = 1/10$; (d) $C_0 = 0.01, \tan \theta = 1/15$; (e) $C_0 = 0.02, \tan \theta = 1/10$; ( f) $C_0 = 0.01, \tan \theta = 1/20$.

The depth-averaged auto-suspension index along the slope $k_{Auto,d}$, which also can be said to be the local auto-suspension index, can be expressed, with reference to (3.11), as follows:

(3.12)\begin{equation} k_{Auto,d}={-}\frac{w_{p,d}\cos\theta}{u_{p,d}^{/{/}}\sin\theta}, \end{equation}

where $w_{p,d}$ is the depth-averaged particle vertical velocity and $u_{p,d}^{//}$ denotes the depth-averaged bed-parallel particle velocity.

Figure 19 shows the temporal and spatial distributions of the local auto-suspension index, with the dashed line representing the TC front. The auto-suspension index, which is distinguished by colour in the figure, is judged by using Parker's standard (red means particles in auto-suspension and blue means particles tending to auto-suspension). It can be seen that the red areas meeting Parker's criterion are mainly near the front of the TC, which demonstrates that it is reasonable to focus on the head dynamics when studying the auto-suspension of the TC (Sequeiros et al. 2009). The red auto-suspension areas are not sensitive to the particle concentration and slope angle. The blue areas expand to some extent when the initial particle concentration or slope rises, especially that area in case 4 shown in figure 19(d) that expands a lot. At approximately $t=4$, the spatial expansion of the auto-suspension area reaches its maximum, which shows a good agreement with figure 11. Interestingly, there exists a special blue region away from the dashed line in the range of $x=0\sim 2$ at $t=6\sim 8$, in which the particles tend to auto-suspension. The main reason might be the influence of the reflected wave generated by the backward flow (Bonnecaze et al. 1993; Blanchette et al. 2005; Nasr-Azadani & Meiburg 2014).

Figure 19. Temporal and spatial variations of local auto-suspension index for various cases. Black dotted line indicates the front position of the TC. (a) Case 1 $C_0 = 0.01, \tan \theta = 1/10$. (b) Case 2 $C_0 = 0.015, \tan \theta = 1/10$. (c) Case 3 $C_0 = 0.02, \tan \theta = 1/10$. (d) Case 4 $C_0 = 0.01, \tan \theta = 1/5$. (e) Case 5 $C_0 = 0.01, \tan \theta = 1/15$. ( f) Case 6 $C_0 = 0.01, \tan \theta = 1/20$.

4. Summary and conclusions

A TC can transport for long distances on very gentle slopes, which is inseparable from the auto-suspension mechanism behind it. This paper employs the LES-DEM model to simulate lock-exchange TCs on inclined slopes. The feasibility of the LES-DEM model for simulating TCs on inclined slopes is examined though the quantitative comparison of the temporal variation of the front position with the experimental result (Gladstone et al. 1998), and the quantitative comparison of the fluid velocity profile with an empirical formula. The auto-suspension mechanism of the current is explained from the perspective of the particle flow process and fluid–particle interactions, and the dynamic response of the auto-suspension process with respect to the two key parameters, initial particle concentration $C_0$ and terrain slope angle $\theta$, is discussed.

The value of $\overline {{\textit {Re}}_p}$ of all transported particles is negatively correlated with the ${\textit {Re}}$ of the current. The auto-suspension particles mainly appear near the current head during the current evolution. The auto-suspension particle quantity and area exhibit a tendency to increase first and then decrease, and have a high positive correlation with the coherent vortical structure near the head. The average particle Reynolds number $\overline {{\textit {Re}}_p}$ of auto-suspension particles can remain approximately unchanged as the TC evolves downstream, and is larger than the non-auto-suspension particle Reynolds number that gradually decreases.

The movements of auto-suspension particles during transport are induced by the complicated physical regimes. The particle trajectories and forces are analysed. The rising auto-suspension particles mainly occur near the current front, reach the peak away from the slope, finally enter the lower layers of the current and exit the auto-suspension state. For the auto-suspension particles, the bed-normal lift force and added mass force are mostly positive. The effective drag force (the combination of effective gravity and drag force) and contact force are much larger in positive total force particles than those in negative ones. Throughout the process, particles are constantly entering and leaving the auto-suspension state. It is precisely the long-term existence of this mechanism that allows the TC to transport over long distances.

An increase in particle concentration can significantly increase the positive bed-parallel total force and positive bed-normal total force. The former improves the conversion rate of energy to particle kinetic energy in the TC system, whereby the current achieves an increase in advance velocity. The latter is generated by the increase of positive lift force, added mass force and contact force, which are essential to help to achieve and maintain particle auto-suspension. As a result, the increase in particle concentration drives the current to meet the auto-suspension criterion more easily. For different small slope conditions, the change in slope has an insignificant effect on both the depth-averaged auto-suspension index and the bed-parallel total force and bed-normal total force. Both negative effective drag force and positive lift force can be increased to some extent, but such enhancement is a mutually constrained balance. Additionally, the conversion rates of particle gravitational potential energy and fluid potential energy grow appreciably as the slope increases.

Particle entrainment (re-suspension) on the erodible terrain, provided with more prominent non-Newtonian dynamics, can provide more potential energy for TCs. The effect of this on particle auto-suspension investigated by fluid–particle coupled models, including particle transport intensity, system energy input and transformation, etc., is a work that is worth looking forward to and meaningful in the future.