Progress in Earth and Planetary Science
High-resolution simulations of turbidity currents
We employ direct numerical simulations of the three-dimensional Navier-Stokes equations, based on a continuum formulation for the sediment concentration, to investigate the physics of turbidity currents in complex situations, such as when they interact with seafloor topography, submarine engineering infrastructure and stratified ambients. In order to obtain a more accurate representation of the dynamics of erosion and resuspension, we have furthermore developed a grain-resolving simulation approach for representing the flow in the high-concentration region near and within the sediment bed. In these simulations, the Navier-Stokes flow around each particle and within the pore spaces of the sediment bed is resolved by means of an immersed boundary method, with the particle-particle interactions being taken into account via a detailed collision model.
Turbidity currents
Navier-Stokes simulations
Continuum formulation
Grain-resolving simulations
Turbidity currents are particle-laden flows in the ocean that are driven by gravity (Meiburg and Kneller 2010). Particle concentrations are usually sufficiently low far away from the sediment bed so that particle-particle interactions play a small or negligible role throughout most of the body of the current. In this region, the Boussinesq approximation of the Navier-Stokes equations, in conjunction with a continuum formulation for the sediment concentration, is well-suited to capture the dynamics of the flow. However, near the sediment bed particle concentrations can be very high, which can potentially result in complex non-Newtonian behavior, hindered settling, and other effects. Here, we describe the above two different simulation approaches, along with representative results, which open up a path towards multiscale flow simulations via the μ(I) rheology (Cassar et al. 2005; Boyer et al. 2011; Aussillous et al. 2013).
Continuum approach
Physical model and governing equations
In many situations of interest, compositional gravity currents and turbidity currents are driven by small density differences not exceeding O(1%). Under such conditions, the Boussinesq approximation can be employed, which treats the density as constant in the momentum equation with the exception of the body force terms. When dealing with turbidity currents, we account for the dispersed particle phase by means of a Eulerian-Eulerian formulation, which means that we employ a continuum equation for the particle concentration field, rather than tracking particles individually in a Lagrangian fashion.
In the following, it will be important to carefully distinguish between dimensional and dimensionless variables. Towards this end, we will employ the tilde symbol to indicate a dimensional variable, whereas variables without the tilde symbol are dimensionless. Under the Boussinesq approximation, the dimensional governing equations for compositional gravity currents driven by salinity and/or temperature gradients can be written as
$$\begin{array}{@{}rcl@{}} \frac {\partial \widetilde{u}_{j}}{\partial \widetilde{x}_{j}} &=& 0, \end{array} $$
$$\begin{array}{@{}rcl@{}} \frac{\partial \widetilde{u}_{i} }{\partial \widetilde{t}} +\frac{\partial\left(\widetilde{u}_{i}\widetilde{u}_{j}\right)}{\partial \widetilde{x}_{j}} &=& -\frac{1}{\widetilde{\rho}_{1}} \frac{\partial {\widetilde{p}} }{\partial \widetilde{x}_{i}} + \widetilde{\nu}\frac{\partial^{2} \widetilde{u}_{i}}{\partial \widetilde{x}_{j} \partial \widetilde{x}_{j}} +\frac{\widetilde{\rho}\widetilde{g}}{\widetilde{\rho}_{1}}e_{i}^{g}, \end{array} $$
$$\begin{array}{@{}rcl@{}} \frac{\partial {\widetilde{\rho}} }{\partial \widetilde{t}} +\frac{\partial\left({\widetilde{\rho}}{\widetilde{u}}_{j}\right)}{\partial \widetilde{x}_{j}} &=& \widetilde{\alpha}\frac{\partial^{2} {\widetilde{\rho}}}{\partial \widetilde{x}_{j} \partial \widetilde{x}_{j}}\ . \end{array} $$
Here, \(\widetilde {u}_{i}\) denotes the velocity vector, \(\widetilde {p}\) the pressure, \(\widetilde {\rho }\) the density, \(\widetilde {g}\) the gravitational acceleration, \(e_{i}^{g}\) the unit vector pointing in the direction of gravity, \(\widetilde {\nu }\) the kinematic viscosity, and \(\widetilde {\alpha }\) the molecular diffusivity of the density field. We nondimensionalize the above Eqs. (1)–(3) by a reference length scale, such as the domain half height \(\widetilde {H}/2\) of a lock-exchange flow (Meiburg et al. 2015; Nasr-Azadani and Meiburg 2014; Necker et al. 2002, 2005), the current density \(\widetilde {\rho }_{1}\), and the buoyancy velocity \(\widetilde {u}_{b}\)
$$ \widetilde{u}_{b} = \sqrt{\widetilde{g'}\;\widetilde{H}/2} \ . $$
Here, \(\widetilde {g'}\) indicates the reduced gravity
$$ \widetilde{g'}=\widetilde{g} \;\frac{\widetilde{\rho}_{1}\; - \widetilde{\rho}_{2}}{\widetilde{\rho}_{1}} \ . $$
where \(\widetilde {\rho }_{2}\) represents the ambient density. After nondimensionalization, we obtain
$$\begin{array}{@{}rcl@{}} \frac {\partial u_{j}}{\partial x_{j}} &=& 0, \end{array} $$
$$\begin{array}{@{}rcl@{}} \frac{\partial u_{i} }{\partial t} +\frac{\partial\left(u_{i}u_{j}\right)}{\partial x_{j}} &=&-\frac{\partial {p} }{\partial x_{i}} + \frac{1}{Re}\frac{\partial^{2} {u}_{i}}{\partial x_{j} \partial x_{j}} +\rho e_{i}^{g}, \end{array} $$
$$\begin{array}{@{}rcl@{}} \frac{\partial {\rho} }{\partial t} +\frac{\partial\left({\rho}{u}_{j}\right)}{\partial x_{j}} &=&\frac{1}{ReSc}\frac{\partial^{2} {\rho}}{\partial x_{j} \partial x_{j}} \ . \end{array} $$
Here, the nondimensional pressure p and density ρ are given by
$$ p = \frac{\widetilde{p}}{\widetilde{\rho}_{1}\widetilde{u}_{b}^{2}} \, \ \ \ \rho = \frac{\widetilde{\rho}-\widetilde{\rho}_{2}}{\widetilde{\rho}_{1}\; - \widetilde{\rho}_{2}} \ . $$
The nondimensionlization of the governing equations gives rise to two dimensionless parameters in the form of the Reynolds number Re and the Schmidt number Sc
$$ Re = \frac{\widetilde{u}_{b} \widetilde{H}}{2\widetilde{\nu}} \, \ \ \ Sc=\frac{\widetilde{\nu}}{\widetilde{\alpha}} \ . $$
While the Reynolds number indicates the ratio of inertial to viscous forces, the Schmidt number represents the ratio of kinematic fluid viscosity to molecular diffusivity of the density field.
When the driving density difference is due to gradients in particle loading, rather than salinity or temperature gradients, the above set of equations no longer provides a full description of the flow. Particles settle within the fluid, so that the scalar concentration field no longer moves with the fluid velocity. In addition, particle-particle interactions can result in such effects as hindered settling (Ham and Homsy 1988), increased effective viscosity, and non-Newtonian dynamics (Guazzelli and Morris 2011), thereby further complicating the picture. However, away from the sediment bed, turbidity currents are often quite dilute, with the volume fraction of the suspended sediment phase being well below O(1%). Under such conditions, particle-particle interactions can usually be neglected, so that the particle settling velocity remains the key difference (along with erosion) that distinguishes turbidity currents from compositional gravity currents.
Due to the small particle volume fraction of dilute turbidity currents, the volumetric displacement of fluid by the particulate phase can usually be neglected, allowing us to consider the fluid velocity field to be divergence-free. Rather, the particle-fluid interaction occurs primarily through the exchange of momentum, so that it suffices to account for the presence of the particles in the fluid momentum equation. In the following, we assume that the particle diameter \(\widetilde {d}_{p}\) is smaller than the smallest length scale of the flow, such as the Kolmogorov scale in turbulent flow. In addition, we consider only particles whose aerodynamic response time \(\widetilde {t}_{p}\) is significantly smaller than the smallest time scale of the flow \(\widetilde {t}_{f}\), so that the particle Stokes number \(St = \widetilde {t}_{p} / \widetilde {t}_{f} \ll O(1)\) (Raju and Meiburg 1995). Here, the aerodynamic response time is defined as
$$ \widetilde{t}_{p} =\frac{\widetilde{\rho}_{p} \widetilde{d}_{p}^{2}}{18 \widetilde{\mu}} \, $$
with \(\widetilde {\rho }_{p}\) indicating the particle material density and \(\widetilde {\mu }\) denoting the dynamic viscosity of the fluid. Such particles can then be assumed to move with a velocity \(\widetilde {u}_{p,i}\) that is obtained by superimposing the local fluid velocity \(\widetilde {u}_{i}\) and the particle settling velocity \(\widetilde {u}_{s} e_{i}^{g}\)
$$ \widetilde{u}_{p,i} = \widetilde{u}_{i} + \widetilde{u}_{s} e_{i}^{g} \, $$
where \(\widetilde {u}_{s}\) follows from balancing the gravitational force with the Stokes drag force
$$ \widetilde{F}_{i} = 3\pi\widetilde{\mu}\widetilde{d}_{p}(\widetilde{u}_{i}-\widetilde{u}_{p,i}) $$
$$ \widetilde{u}_{s} = \frac{\widetilde{d}_{p}^{2}(\widetilde{\rho_{p}}-\widetilde{\rho})\widetilde{g}}{18\widetilde{\mu}} \ . $$
Note that this implies that the particle velocity field is single-valued and divergence-free, so that monodisperse particles do not, for example, accumulate near stagnation points or get ejected from vortex centers. Hence, we can describe the spatio-temporal evolution of the particle number concentration field \(\widetilde {c}\) in a Eulerian fashion by the transport equation
$$ \frac{\partial {\widetilde{c}} }{\partial \widetilde{t}} +\frac{\partial\left(\widetilde{c}\left(\widetilde{u}_{j}+\widetilde{u}_{s}e_{j}^{g}\right)\right)}{\partial \widetilde{x}_{j}} =\widetilde{\alpha}\frac{\partial^{2} \widetilde{c}}{\partial \widetilde{x}_{j} \partial \widetilde{x}_{j}} \ . $$
The diffusion term in Eq. (15) represents a model for the decay of concentration gradients due to the hydrodynamic diffusion of particles and/or slight variations in particle size and shape (Davis and Hassen 1988; Ham and Homsy 1988).
The motion of the fluid phase is described by the incompressible continuity equation and the Navier-Stokes equation augmented by the force exerted on the fluid by the particles, which is equal and opposite to the Stokes drag force acting on the particles. In a dimensional form, these equations read
$$\begin{array}{@{}rcl@{}} \frac{\partial \widetilde{u}_{i} }{\partial \widetilde{t}} +\frac{\partial\left(\widetilde{u}_{i}\widetilde{u}_{j}\right)}{\partial \widetilde{x}_{j}} &=&-\frac{1}{\widetilde{\rho}} \frac{\partial {\widetilde{p}} }{\partial \widetilde{x}_{i}} + \widetilde{\nu}\frac{\partial^{2} \widetilde{u}_{i}}{\partial \widetilde{x}_{j} \partial \widetilde{x}_{j}} +\frac{\widetilde{c}}{\widetilde{\rho}} \widetilde{F}_{i}, \end{array} $$
As we had done for compositional gravity currents, we use the domain half height \(\widetilde {H}/2\) and buoyancy velocity \(\widetilde {u}_{b}\) for nondimensionalization. The reduced gravity \(\widetilde {g}'\) appearing in the calculation of \(\widetilde {u}_{b}\) can now be calculated as
$$ \widetilde{g}' = \frac{\pi(\widetilde{\rho}_{p}-\widetilde{\rho})\widetilde{c}_{0}\widetilde{d}_{p}^{3}} {6\widetilde{\rho}} \ \widetilde{g} \, $$
where \(\widetilde {c}_{0}\) indicates a reference number concentration of particles in the suspension. After nondimensionalization, we obtain
$$\begin{array}{@{}rcl@{}} \frac {\partial u_{j}}{\partial x_{j}} &\,=\,& 0 \, \end{array} $$
$$\begin{array}{@{}rcl@{}} \frac{\partial u_{i} }{\partial t} +\frac{\partial\left(u_{i}u_{j}\right)}{\partial x_{j}} &\,=\,& \,-\,\frac{\partial {p} }{\partial x_{i}} \,+\, \frac{1}{Re}\frac{\partial^{2} {u}_{i}}{\partial x_{j} \partial x_{j}} \,+\,c e_{i}^{g}, \end{array} $$
$$\begin{array}{@{}rcl@{}} \frac{\partial {c} }{\partial t} +\frac{\partial\left({c}\left({u}_{j}+u_{s}e_{j}^{g}\right)\right)}{\partial x_{j}} &\,=\,& \frac{1}{ReSc}\frac{\partial^{2} {c}}{\partial x_{j} \partial x_{j}} \ . \end{array} $$
For polydisperse suspensions containing particles of different sizes, the above approach can easily be extended by solving one concentration equation for each particle size and corresponding settling velocity (Nasr-Azadani and Meiburg 2014). Note that the set of governing equations for turbidity currents (19)–(21) differs from the corresponding set for compositional gravity currents (6)–(8) only by the additional settling velocity term in the concentration equation. In the following, we employ Eqs. (19)–(21) for both types of currents, with the tacit assumption that the settling velocity vanishes for compositional gravity currents.
Direct numerical simulations (DNS) represent the most accurate computational approach for studying gravity currents. In DNS, all scales of motion, from the integral scales dictated by the boundary conditions down to the dissipative Kolmogorov scale determined by viscosity, are explicitly resolved. However, for the case of turbidity currents, when the particle diameter is smaller than the Kolmogorov scale, the fluid motion around each particle is usually not resolved, due to the prohibitive computational cost. Nevertheless, the drag law accurately captures the exchange of momentum between the two phases at scales smaller than the Kolmogorov scale, so that the approach described above is still referred to as DNS.
Consistent with the above arguments, the grid spacing required for DNS is of the order of the Kolmogorov scale, while the time step needs to be of the same order as the time scales of the smallest eddies. Due to the large disparity between integral and Kolmogorov scales at high Reynolds numbers, the computational cost of DNS scales as Re 3, so that the DNS approach is effectively limited to laboratory scale Reynolds numbers. The first DNS simulations of gravity currents in a lock-exchange configuration were reported by Härtel et al. (2000) for Re = 1225. Necker et al. (2002) extended this work to turbidity currents at Re = 2240. More recent simulations of lock-exchange gravity currents by Cantero et al. (2008) were able to reach Re = 15,000, which corresponds to a laboratory scale current of height 0.5 m with a front velocity of 3 cm/s.
DNS simulations can provide detailed information on the structure and statistics of the flow, on the various components of its energy budget, on the mixing behavior, and many additional aspects. As a case in point, the simulations by Härtel et al. (2000) explored the detailed flow topology near the current front and demonstrated that the stagnation point is located a significant distance behind the nose of the current. DNS results are furthermore very useful for testing the accuracy and identifying any deficiencies in larger-scale LES and RANS models (Yeh et al. 2013). Thus, while they are currently limited to laboratory scale currents, DNS simulations represent an excellent research tool for exploring the detailed physics of moderate Reynolds number gravity currents and for constructing larger-scale models for higher Reynolds number applications.
We illustrate the ability of the continuum approach to reproduce lab-scale experiments by presenting the results of highly resolved simulations of a turbidity current moving down a slope into a stratified saline ambient. The numerical setup directly replicates the experiments conducted by Snow and Sutherland (2014) and is presented in Fig. 1. The density inside the ambient increases linearly from ρ T at the top to ρ B at the bottom such that
Problem setup and configuration. a Particle-laden fluid. b Ambient stratified fluid. c Solid region. H and L denote the height and length of the ramp, L d is the horizontal length of the domain, h 0 and L l denote the height and horizontal length of the lock, and m is the slope. ρ 1 is the bulk density of the lock, ρ T is the initial density at the top of the ambient, ρ B is the initial density at the bottom of the ambient, and ρ A is the initial density in the ambient at half-lock depth
$$ \rho_{2}(y)=\rho_{B}+\left(\rho_{T}-\rho_{B}\right)\cdot\frac{y}{H} $$
The channel has a constant width denoted as W. The lock region is initially at rest with density ρ 1, chosen such that ρ B > ρ 1 > ρ T . At t=0, the lock is released and the particle-laden flow moves down the slope forming a turbidity current interacting with the ambient fluid. Here, both the particle concentration c and salinity s contribute to the Boussinesq term in Eq. (20) such that Eq. (21) has to be solved for each scalar field, with u s = 0 in the case of the salinity field. When a settling velocity is used, the particles (concentration field) are allowed to settle through the lower boundary so that an erodible bed is not formed. We define a nondimensional buoyancy frequency N to quantify stratification such that
$$ N=\sqrt{\frac{\tilde{\rho}_{0}(\tilde{\rho}_{B}-\tilde{\rho}_{T})}{\tilde{\rho}_{T}(\tilde{\rho}_{1}-\tilde{\rho}_{A})}\frac{\tilde{H}}{\tilde{h}_{0}}}, $$
where \(\tilde {\rho }_{0}\) is a reference density.
A finite difference code is used to solve the equations, and the MPI library is used for parallelization. A third-order Runge-Kutta scheme with a three substep method is used to discretize the equations in time. The wall-normal viscous and diffusive terms are solved implicitly while the convective terms and the remaining viscous and diffusive terms are treated explicitly. To impose incompressibility, a projection method is used (Spalart et al. 1991) and a direct solver is used for the resulting Poisson equation. Slip-wall boundary conditions are used at the top and right walls and z-periodic boundary conditions are used at the lateral walls. The domain is assumed to be sufficiently long to neglect boundary effects in x, and the width of the domain is chosen so that the periodic boundary condition does not impact the flow development. Finally, an immersed boundary method is used to impose the no-slip condition on the slope (Nasr-Azadani and Meiburg 2011).
In Fig. 2, we present a time series depicting the evolution of a typical turbidity current moving down a slope and intruding when its density matches that of the stratified ambient. The spanwise-averaged particle concentration is represented on a linear gray scale for various times. Upon release of the lock, the current starts moving down the slope and a trail of large Kelvin-Helmholtz rollers forms in the tail. These large instabilities then break into fully three-dimensional turbulence creating smaller dissipative vortices (t > 10). The absence of large distinct structures indicates the presence of fully developed turbulence in the tail of the current.
Spanwise averaged particle concentration pseudo-color plot for various times (Re = 15000, v s = 0.001, N = 2.11, m = 0.0744). This represents the case of a turbidity current with non-zero settling velocity moving down a uniformly stratified ambient. The density of the lock is initially larger than the ambient fluid outside of the lock such that the current starts moving downslope upon release. It then encounters region of denser and denser saline ambient, until its density reaches that of the ambient and it intrudes horizontally. The flow is 2D at very early stages but quickly grows unstable under the effect of a spanwise instability at the head. This leads to the breakdown of 2D vorticies into much smaller 3D mixing structures at the head and all along the body of the current at the interface with the ambient
The direct impact of stratification is seen at later times (t≈15) when the current intrudes into the ambient, i.e. separates from the surface of the slope. The effects of stratification on intrusion depth are key in understanding the evolution of the suspended mass, deposition profiles and energy budgets of turbidity currents. Intrusion only occurs when the density of the current reaches the density of the ambient.
Using numerical simulations, we are able to investigate the fundamental mechanisms that control the propagation of the turbidity current and monitor all the relevant dynamic variables. For instance, we can investigate the initial perturbation that leads to the spanwise breakdown of the large Kelvin-Helmoltz structures that initially appear in the tail of the current. Figure 3 is a representation of the concentration isosurface c = 0.25 at t 1 = 10 for a typical turbidity current at Re = 6000. The flow at that instant is not yet fully turbulent but displays strong spanwise instabilities characterized at the head by the lobe-and-cleft instability. This instability is responsible for the breakdown of the large Kelvin-Helmholtz rollers that only high-resolution 3D simulations are able to capture.
Concentration isosurface c = 0.25 at t = 10 (Re = 6000, v s = 0.005, N = 2.09, m = 0.0744). The iso-surface reveals the global structure of the envelope of a hyperpycnal current upon destibalization. The lobe-and-cleft instability is clearly visible at the head of the current as hills and crests form in the spanwise direction. This instability propagates in the body of the current and destabilises the large 2D vortical structures initially present. The largest wavelength corresponds to the initial dominant mode of the instability and dictates the width of the domain necessary to observe the instability at this given Reynolds number
Quantitative analysis of the velocity of the current as a function of the buoyancy frequency shows very good agreement between the numerical simulations and the experiments. The front velocity was measured by Snow and Sutherland (2014) in a series of experiments of which we report three. The margin of error for the experimental measurements is typically of the order of ±20% while the relative difference with the numerical simulations was found to be of 14, − 11, and 16% for three widely different settling velocities, Reynolds numbers, and buoyancy frequencies. The simulation parameters and relative error for the front velocity are summarized in Table 1. Numerical results also agree almost perfectly with analytical results in the limit of no stratification, where the Froude number is expected to be Fr = 0.5. Numerical results consistently yielded a Froude number of Fr = 0.496 for N = 0, which corresponds to a relative error of 0.8% when compared to the analytical result.
Comparison of numerically and experimentally measured front velocity
v s
(U e − U s )/U e
− 11%
Data from Snow and Sutherland (2014). The experimentally measured and numerically computed front velocities are denoted as U e and U s respectively
The depth at which the current intrudes into the ambient also reveals a good agreement between the two approaches and validates the ability of numerical simulations to reproduce the dynamical features of turbidity currents moving into a stratified ambient. While it is extremely challenging to experimentally measure the velocity, particle concentration, and salinity fields of such 3D turbulent flows, direct numerical simulations give access to an entirely new set of data and opens the door to more accurate prediction tools and a deeper understanding of the underlying physics of gravity and turbidity currents in realistic environments at the scale of the lab.
Grain-resolving approach
When the concentration of particles grows large, particle-particle interactions become important and the aforementioned continuum approach is no longer applicable. For such cases, we have to account for the rheology of dense suspensions. A key element of progress with regard to the rheology of dense suspensions over the last decade has been the development of the so-called μ(I) approach, cf. (Guazzelli and Morris 2011; Boyer et al. 2011). Grain-resolving simulations of the type to be discussed in the following are expected to provide a tool for further investigating the validity of the assumptions underlying the derivation of the μ(I) rheology. One way to approach the simulation of dense suspensions is to fully resolve the particles interacting with the fluid by tracking each individual particle, evaluating the fluid no-slip condition at the particle surface and accounting for all the forces acting on the particle. Such simulations typically use grid resolutions of 10–25 grid cells per particle diameter to resolve the flow and are thus limited in scope to domains whose dimensions measure only tens to hundreds of diameters in length (Vowinckel et al. 2017; Kidanemariam and Uhlmann 2017). Thus, if simulating sand grains 100 μm in diameter, the domain dimensions would range in length from 1 mm to 1 cm. However, the idea is to use grain-resolved simulations to develop better models of sediment transport to be used in larger-scale simulations.
We employ an immersed boundary method (IBM) developed by Uhlmann (2005) and Kempe and Fröhlich (2012) with a modified collision model adjusted for the present context, as explained below. This method solves the Navier-Stokes equations everywhere in the domain, including nearby and within the particles:
$$\begin{array}{@{}rcl@{}} \frac{\partial{\mathbf{u}}}{\partial{t}}+\nabla\cdot(\mathbf{u}\mathbf{u}) = -\frac{1}{\rho_{f}}\:\nabla p + \nu_{f} \nabla^{2} \mathbf{u} + \mathbf{f}_{{IBM}} \end{array} $$
where ρ f is the fluid density and f IBM is the IBM force, which acts as a source term to enforce the no-slip condition at the particle surfaces. This force effectively couples the particle and fluid momentum equations. Though there are many ways to carry out this coupling, the method we employ for the particles uses regularized Dirac delta functions, which interpolate fluid velocities onto the particle surface and spread f IBM onto the fluid (Roma et al. 1999). Note that this implementation is different from that used to create the sloped lower wall in the turbidity current simulations of the previous section.
The equations of motion for the particles are given by the momentum equations for the translational velocity u p = (u p ,v p ,w p ) T
$$ m_{p} \: \frac{\mathrm{d}\mathbf{u}_{p}}{\mathrm{d} t} = \underbrace{\oint_{\Gamma_{p}} \boldsymbol{\tau} \cdot \mathbf{n}\: {\mathrm{d}A}}_{=\mathbf{F}_{h}} + \underbrace{V_{p}\:(\rho_{p}-\rho_{f})\: \mathbf{g}}_{=\mathbf{F}_{g}} + \mathbf{F}_{l} + \mathbf{F}_{c}, $$
the angular velocity ω p =(ω p,x ,ω p,y ,ω p,z ) T
$$ I_{p} \:\frac{\mathrm{d}\boldsymbol{\omega}_{p}}{\mathrm{d} t} = \underbrace{\oint_{\Gamma_{p}} \mathbf{r}\times(\boldsymbol{\tau}\cdot\mathbf{n})\:{\mathrm{d}A}}_{=\mathbf{T}_{h}} + \mathbf{T}_{l} + \mathbf{T}_{c}, $$
and the position x p =(x p ,y p ,z p ) T
$$ \frac{\mathrm{d}\mathbf{x}_{p}}{\mathrm{d} t} = \mathbf{u}_{p}. $$
Here, m p is the particle mass, Γ p the fluid-particle interface, τ the hydrodynamic stress tensor, ρ p the particle density, V p the particle volume, g the gravitational acceleration, \(I_{p}~=~8\pi \rho _{p} R_{p}^{5}/15\) the moment of inertia, and R p the particle radius. Furthermore, the vector n is the outward-pointing normal on the interface Γ p , r = x − x p is the position vector of the surface point with respect to the center of mass x p of a particle, F l and T l are the force and torque due to lubrication forces, and F c and T c are the force and torque due to particle collisions. We evaluate the IBM force f IBM as well as the hydrodynamic force, F h , and torque, T h , using the approach of Kempe and Fröhlich (2012), fully resolving the hydrodynamic effects of the fluid on the particles as well as the particles on the fluid. The lubrication force, F l , and contact force, F c , model close-range particle-particle interactions. With the exception of the tangential lubrication force, the methods used to evaluate these forces are described and validated in detail by Biegert et al. (2017), but here, we present them briefly.
The lubrication force
$$\begin{array}{@{}rcl@{}} \mathbf{F}_{l} =&-&6 \pi \rho_{f} \nu_{f} R_{\text{eff}} \left(\frac{R_{\text{eff}}}{\max(\zeta_{n},\zeta_{\text{min}})} \mathbf{g}_{n} + F^{*}_{t} \mathbf{g}_{t}\right. \\ &+&F_{r}^{*} \left.(R_{p} \boldsymbol{\omega}_{p}\times \mathbf{n} + R_{q} \boldsymbol{\omega}_{q}\times \mathbf{n}){\vphantom{\frac{R_{\text{eff}}}{\max(\zeta_{n},\zeta_{\text{min}})}}}\right) \end{array} $$
and torque
$$ {}\mathbf{T}_{l} = 8 \pi \rho_{f} \nu_{f} R_{\text{eff}}^{2} \left[\mathbf{g}_{t} T_{t}^{*}+T_{r}^{*} (R_{p} \boldsymbol{\omega_{p}} \times \mathbf{n} + R_{q} \boldsymbol{\omega}_{q}\times \mathbf{n}) \right]\times \mathbf{n} \qquad $$
are added to account for short-range hydrodynamic forces that are unresolved by the fluid grid. Here, R eff = R p R q /(R p + R q ) is an effective radius accounting for size differences between particles p and q, g n and g t are the relative velocities in the normal and tangential directions, respectively, between the two particle surfaces at the point of contact, ζ n is the surface distance between the two particles, and ζ n,min = 3 × 10−3 R p is a limiter preventing the lubrication force from reaching its singularity at ζ n →0. The terms \(F_{t}^{*}\), \(F_{r}^{*}\), \(T_{t}^{*}\), and \(T_{r}^{*}\) were obtained via asymptotic expansions by Goldman et al. (1967):
$$\begin{array}{@{}rcl@{}} F_{t}^{*} &\sim& \frac{8}{15} \ln\left(\frac{\max(\zeta_{n},\zeta_{\text{min}})}{R_{\text{eff}}}\right) - 0.9588 \end{array} $$
$$\begin{array}{@{}rcl@{}} F_{r}^{*} &\sim& -\frac{2}{15} \ln\left(\frac{\max(\zeta_{n},\zeta_{\text{min}})}{R_{\text{eff}}}\right) - 0.2526 \end{array} $$
$$\begin{array}{@{}rcl@{}} T_{t}^{*} &\sim& -\frac{1}{10} \ln\left(\frac{\max(\zeta_{n},\zeta_{\text{min}})}{R_{\text{eff}}}\right) - 0.1895 \end{array} $$
$$\begin{array}{@{}rcl@{}} T_{r}^{*} &\sim& \frac{2}{ 5} \ln\left(\frac{\max(\zeta_{n},\zeta_{\text{min}})}{R_{\text{eff}}}\right) - 0.3817. \end{array} $$
As indicated in (25) and (26), we also account for particle-particle contacts through F c and T c . These contact forces are composed of components normal and tangential to the particle surface, represented by F n and F t , respectively, which act at the point of contact between the two particles so that the resulting force and torque on the particle are given by
$$\begin{array}{@{}rcl@{}} \mathbf{F}_{c} &=& \mathbf{F}_{n} + \mathbf{F}_{t} \end{array} $$
$$\begin{array}{@{}rcl@{}} \mathbf{T}_{c} &=& R_{p} \, \mathbf{n} \times \mathbf{F}_{t}, \end{array} $$
where n is the outward-pointing normal vector from the contact point. A nonlinear spring-dashpot model is used for the normal contact force
$$ \mathbf{F}_{n} = -k_{n} |\zeta_{n}|^{3/2} \mathbf{n} - d_{n} \mathbf{g}_{n}, $$
where the stiffness and damping coefficients, k n and d n , respectively, are adaptively calibrated for every collision. A linear spring-dashpot model is used for the tangential contact force
$$ \mathbf{F}_{t} = \min \left(-k_{t} \boldsymbol{\zeta}_{t} - d_{t} \mathbf{g}_{t}, ||\mu \mathbf{F}_{n}|| \mathbf{t} \right), $$
where ζ t is the tangential displacement vector representing accumulated slip between the two surfaces, μ is the coefficient of friction between the surfaces, and t is the unit normal vector in the tangential direction. The Coulomb friction criterion, represented by ||μ F n ||, allows the two surfaces to slip past one another when large stresses are present. Similar to the normal coefficients, the tangential coefficients of stiffness and damping, k t and d t , are also adaptively calibrated.
Grain-resolving results
Pressure-driven flow over dense sediment
To address the bulk behavior of a dense granular bed sheared by a laminar Poiseuille flow, we carried out numerical simulations to reproduce the experimental results of Aussillous et al. (2013), who studied pressure-driven flows over glass spheres with a mean diameter D p = 1.1 mm and a standard deviation of σ(D p ) = 0.1 mm as a sediment material. This experimental work provides investigations over a range of submergences h f /D p and Reynolds numbers in the laminar regime, where h f is the height of the clear-water layer above the sediment bed illustrated in Fig. 4. We define h f to be the height above which the average particle volume fraction ϕ < 0.05, which is the threshold for negligible impact of particle-particle interaction on the flow (Capart and Fraccarollo 2011).
Pressure-driven flow over a dense sediment bed. The left figure shows an instantaneous slice through the domain, the color scale illustrating streamwise fluid velocity, and the gray scale illustrating particle velocity. The right figure shows the time-and-space-averaged streamwise velocity of the fluid and particles
We executed several simulations in an attempt to match the experimental results of Aussillous et al. (2013) at different flow rates and fluid heights. To this end, we simulated a monodisperse granular sediment bed sheared by a pressure-driven Poiseuille flow with periodic conditions in both the streamwise (x) and spanwise (z) directions, respectively. A no-slip condition was applied at the top and bottom wall as well as at the particle surface. A detailed comparison validating the simulation results has been presented in great detail in Biegert et al. (2017). Here, we extend this work to show data from the same physical setup, but with an increased flow rate, such that the gross of the particles are set into motion and interact in a complex network. The physical and numerical parameters associated with this simulation are listed in Table 2.
Simulation parameters of the pressure-driven flow scenario, where \(u_{\tau }=\sqrt {\tau _{w}/\rho _{f}}\) is the friction velocity at the fluid/particle interface, U b is the bulk (average) velocity of the fluid, τ w is the shear stress at the fluid/particle interface, ν f the kinematic viscosity, L x , L y , and L z are the spatial extents of the computational domain in the Cartesian space, and h is the grid cell size
Re = U b L y /ν f
D + = u τ D p /ν f
Sh = τ w /[(ρ p −ρ f )gD]
ρ p /ρ f
L x × L y × L z
11.26D p × 22.52D p × 11.26D p
h f /D p
D p /h
Some qualitative inferences can already be drawn from Fig. 4, where we have nondimensionalized velocities by the bulk fluid velocity \(U_{b}~=~\frac {1}{L_{y}}\int _{0}^{L_{y}} u \, \mathrm {d}y\). In the clear-water layer, a parabolic profile obeying the analytical solution of the classical Poiseuille flow can be observed. The lower end of this parabolic region, however, is not a no-slip wall, but a moving granular bed, which causes the symmetry axis of the flow profile to shift from h f /2 to a lower position. Inside the granular bed, a linear shear flow profile develops and since all particles are moving, this profile continues all the way to the bottom wall of the domain. This interesting behavior and the wealth of data obtained from the grain-resolving simulations opens up a wide range of analytical tools in terms of statistical description as well as physical modeling, which will be our focus in the future.
To further investigate the rheologic properties of dense suspensions, we simulated a shear flow of two parallel walls with a spacing of H moving in opposite directions. A no-slip condition is applied for the moving walls at the top and bottom of the domain and at the particle surface. The walls move with a relative velocity Δ U w = U t − U b , where U t is the velocity of the top wall and U b is the velocity of the bottom wall. Periodic conditions are applied in streamwise (x) and spanwise (z) directions. The fluid in between the two walls is a dense mixture of spherical particles, where “dense” indicates a volume fraction ϕ v = N p V p /V f > 0.4. Here, \(V_{p}~=~4\, \pi R_{p}^{3} / 3\) is the volume of a single particle, N p is the number of particles in the computational domain V f , and R p is the particle radius. The dimension of the computational domain can be treated as part of the physical problem. We have considered two computational domains: scenario A with dimensions L x × L y × L z = 2H × 1H × 1H and scenario B with L x × L y × L z = 1H × 2H × 1H. Thus, the relative submergence becomes H/D p = 10 and H/D p = 20, respectively, where D p is the particle diameter. For both scenarios A and B, the shear rate \(\dot {I}~=~\Delta U_{w} / H\) was kept constant so that the particle Reynolds number \(Re_{p}~=~\dot {I} D_{p}^{2}/\nu _{f}~=~0.1\) also stays constant but the channel Reynolds number for scenario B is increased by a factor of 4 with respect to scenario A. This yields channel Reynolds numbers of Re = 10 (scenario A) and Re = 40 (scenario B). The particles have a density of ρ p /ρ f = 1.011, which is close to neutrally buoyant conditions. We assume material parameters corresponding to glass or silicate materials as thoroughly validated in Biegert et al. (2017). In particular, we choose e dry = 0.97, μ k = 0.15, μ s = 0.8, and ν = 0.22, where e dry is the wall-normal restitution coefficient for dry collisions, μ k and μ s are the kinetic and static friction coefficients for oblique collisions, and ν is Poisson’s ratio. Every particle is discretized by 25.6 grid cells per diameter. The particle parameters are summarized in Table 3. Three simulations were conducted with varying volume fractions of the mixture and different relative gap sizes H/D p to explore the effects of these two parameters. All simulations were initialized and run for a start-up time t s until a true steady steady state had been established. Subsequently, data was collected for the averaging time t a to reach converged statistics for the profiles presented in the following. The different scenarios are displayed in Table 4.
e dry
μ k
μ s
ζ min
H/D p
3·10−3 R p
Simulation scenarios
ϕ v
\(t_{s} \dot {I}\)
\(t_{a} \dot {I}\)
Re10p42
2H × 1H × 1H
An instantaneous snapshot of the particle distribution colored by the particle velocity is given in Fig. 5. As desired, particles are dragged along the moving walls whenever they collide with them. These particles moving with the wall transfer kinetic energy through collisions towards the channel center. In addition, the moving walls establish a background profile for the fluid velocity, which should be close to the linear shear profile commonly observed in Couette-type flows. The two mechanisms from collision and hydrodynamic interactions establish a shear flow profile within the suspension. Looking at the wall-normal profiles of the porosity, we can see a distinct pattern of oscillations (Fig. 6). This crystal-like layering of the particles reflects the fact that all particles are the same size. While a strongly layered structure is visible for all three simulations close to the wall, less pronounced layers form in the channel center for the two scenarios Re10p42 and Re40p54. For these two scenarios, particles have more space to rearrange due to the lower volume fraction and the larger relative gap size, respectively. Horizontal-and-time-averaged profiles of the streamwise component of the fluid and particle velocity show that particles move with almost the same velocity as the fluid flow (Fig. 7). Slight distortions can be seen, especially for case Re10p54, illustrating local effects of the particle clustering on the global velocity profile.
Shearing of dense suspensions. Color bar indicates streamwise velocity of particles. Shown are instantaneous snapshots for a Re = 10 and ϕ v = 0.42, b Re = 10 and ϕ v = 0.54, and c Re = 40 and ϕ v = 0.54
Porosity profiles averaged in the streamwise and spanwise directions and time for the three simulation scenarios: a Re = 10 and ϕ v = 0.42, b Re = 10 and ϕ v = 0.54, and c Re = 40 and ϕ v = 0.54
Average streamwise particle velocity profiles for the three simulation scenarios: a Re = 10 and ϕ v = 0.42, b Re = 10 and ϕ v = 0.54, and c Re = 40 and ϕ v = 0.54
The present study of a Couette-type flow supplements our simulations of pressure-driven flow described in the previous section to fully understand the rheologic behavior of dense suspensions of particles with different inertia in flows with different momentum supply.
We also studied the hydrodynamic forces acting on a fully resolved sediment bed induced by a gravity current. Here, the key issue is to explore how a jump in the hydrostatic pressure traveling along the surface of the sediment bed propagates within the bed. To this end, we solved the advection-diffusion Eq. (8) and supplemented (24) with the term stemming from the Boussinesq approximation. The Peclet number was chosen to be Pe = ReSc = 104. The initial configuration is similar to a lock-exchange, but the particles are submerged in a layer of high concentration as well. The geometry is L x /H × L y /H × L z /H = 6 × 1 × 1, where H is the channel height from the no-slip wall on the bottom to the free-slip wall on the top, which is equivalent to L y . The Reynolds number for the present scenario was chosen to be Re = u b h/ν f = 164, where \(u_{b}~=~\sqrt {g' h}\) is the buoyancy velocity, h is the lock height, \(g' = \frac {\rho _{p} - \rho _{f}}{\rho _{f}} g\) is the specific gravity, and g is the gravitational acceleration. Free-slip walls were applied to the left, back, and front wall (Fig. 8 a). The right boundary was set to be a convective outflow condition. The relative submergence of a particle is H/D p = 10. Every particle is discretized by D p /h = 16.
Instantaneous snapshots for the internal wave study. Fixed particles are arranged in a hexagonal packing immersed in the heavy (red) fluid. a Initially, a lock of dense fluid is arranged at one end of the domain, which then b generates a bore propagating to the right at a later time. The vertical slice shows the concentration profile of the heavy (red) and light (blue) fluids. Particles are colored by the lift force acting on them; red and blue indicate positive and negative forces, respectively
Immediately after being released, the block of heavy fluid starts to propagate along the rough wall, forming an internal wave at the interface between light and heavy fluid, which are indicated in Fig. 8 as blue and red fluid, respectively. Particles in Fig. 8 b are colored by the lift forces acting in a vertical direction. We can see that particle in front of the wave start to experience a lift force even though the wave front has not yet reached it. The grain-resolving simulation approach now allows us to track the drag and lift on individual particles as a function of time to elucidate this effect in more detail. This has been done for the particles colored in red in Fig. 9 a. The lift force normalized by the buoyant weight of the particles F g = ρ f V p g is shown in Fig. 9 b. Every curve represents a particle shown in Fig. 9 a, and it becomes obvious that the lift experienced by the particles along the transect appears to be similar with decaying intensity along the flow direction.
Lift/drag of the internal propagating wave acting on the centerline particles. a shows a top view of the particle bed, where lift forces are measured at the red-colored centerline particles. b shows the lift force normalized by the gravitational force, F L /F G , versus time for the centerline particles, each particle represented by a separate curve. c shows the position of the internal wave front versus time, where the position is determined either from the concentration profile, the drag force acting on the particle bed, or the lift force acting on the particle bed
We can track the front of the current using the location in the horizontal profile where the interface between light and heavy fluid starts to increase in height. Alternatively, we can track the front using the location where particles start to experience an enhanced lift force. A comparison of these two methods is shown in Fig. 9 c for the first few time units simulated. It can be seen that, indeed, during the initial stage of the simulation, where we can see a steep front of the internal wave, the force signal propagates quicker through the sediment bed than the actual propagation speed of the current would suggest. This effect, however, levels off over time as the wave continues to travel over the rough bed, constantly losing energy due to viscous dissipation.
The modeling of dilute, non-eroding turbidity currents has reached a mature level, as evidenced by the fact that high-resolution simulations have been able to reproduce many of the observations made in laboratory experiments (e.g., Nasr-Azadani et al. 2013). We are now able to account for some topographical complexity via the immersed boundary method. Some of the remaining challenges concern the extension to the very large Reynolds number values of field-scale flows and the frequent interaction with ambient phenomena in the ocean such as internal waves and tides, as well as the accurate modeling of erosion and resuspension in such high Reynolds number flows. However, a similar level of maturity has not yet been achieved with regard to the modeling of highly concentrated turbidity currents with significant erosion, resuspension, and bedload transport. Especially the dynamics of the near-bed region of such high-concentration currents in the form of dense suspensions is still poorly understood, as it is governed by intense particle-fluid and particle-particle interactions that give rise to strongly non-Newtonian dynamics and to mass and momentum exchanges between the current and the sediment bed. As a result, insight into the erosional and depositional behavior of such currents and the coupling between the motion of the current above the sediment bed and the fluid flow inside the bed is just beginning to emerge. Key progress has been accomplished with regard to understanding the rheology of dense suspensions over the last decade, through the development of the so-called μ(I) approach, cf. (Guazzelli and Morris 2011; Boyer et al. 2011). Grain-resolving simulations based on the approach outlined here will provide a tool for further investigating the validity of the assumptions underlying the derivation of the μ(I) rheology. The computational approach outlined and tested in the present paper holds great promise, as it is able to capture the grain-resolved dynamics of thick, mobile sediment beds and their coupled dynamics with the flow above. Simulations on this basis provide the opportunity to understand erosion and dense suspension rheology from a fundamental perspective, which can lead to better models for use at larger scales. This multiscale approach would thus further enrich our understanding of turbidity currents.
Ham, JM, Homsy GM (1988) Hindered settling and hydrodynamic dispersion in quiescent sedimenting suspensions. Int J Multiphase Flow 14(5):533–546.View ArticleGoogle Scholar
Härtel, C, Meiburg E, Necker F (2000) Analysis and direct numerical simulation of the flow at a gravity-current head. Part 1. Flow topology and front speed for slip and no-slip boundaries. J Fluid Mech 418:189–212.View ArticleGoogle Scholar
Kempe, T, Fröhlich J (2012) An improved immersed boundary method with direct forcing for the simulation of particle laden flows. J Comput Phys 231(9):3663–3684.View ArticleGoogle Scholar
Kidanemariam, AG, Uhlmann M (2017) Formation of sediment patterns in channel flow: minimal unstable systems and their temporal evolution. J Fluid Mech 818:716–743.View ArticleGoogle Scholar
Meiburg, E, Kneller B (2010) Turbidity currents and their deposits. Annu Rev Fluid Mech 42:135–156.View ArticleGoogle Scholar
Meiburg, E, Radhakrishnan S, Nasr-Azadani M (2015) Modeling gravity and turbidity currents: computational approaches and challenges. Appl Mech Rev 67(4):040802.View ArticleGoogle Scholar
Nasr-Azadani, MM, Meiburg E (2011) TURBINS: an immersed boundary, Navier-Stokes code for the simulation of gravity and turbidity currents interacting with complex topographies. Comput Fluids 45(1, SI):14–28.View ArticleGoogle Scholar
Nasr-Azadani, MM, Hall B, Meiburg E (2013) Polydisperse turbidity currents propagating over complex topography: comparison of experimental and depth-resolved simulation results. Comput Geosci 53:141–153.View ArticleGoogle Scholar
Nasr-Azadani, MM, Meiburg E (2014) Turbidity currents intracting with three-dimensional seafloor topography. J Fluid Mech 745:409–443.View ArticleGoogle Scholar
Necker, F, Härtel C, Kleiser L, Meiburg E (2002) High-resolution simulations of particle-driven gravity currents. Int J Multiphase Flow 28(2):279–300.View ArticleGoogle Scholar
Necker, F, Härtel C, Kleiser L, Meiburg E (2005) Mixing and dissipation in particle-driven gravity currents. J Fluid Mech 545:339.View ArticleGoogle Scholar
Raju, N, Meiburg E (1995) The accumulation and dispersion of heavy-particles in forced 2-dimensional mixing layers.2. The effect of gravity. Phys Fluids 7(6):1241–1264.View ArticleGoogle Scholar
Roma, AM, Peskin CS, Berger MJ (1999) An adaptive version of the immersed boundary method. J Comput Phys 153(2):509–534.View ArticleGoogle Scholar
Snow, K, Sutherland BR (2014) Particle-laden flow down a slope in uniform stratification. J Fluid Mech 755:251–273.View ArticleGoogle Scholar
Spalart, PR, Moser RD, Rogers MM (1991) Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions. J Comput Phys 96(2):297–324.View ArticleGoogle Scholar
Uhlmann, M (2005) An immersed boundary method with direct forcing for the simulation of particulate flows. J Comput Phys 209(2):448–476.View ArticleGoogle Scholar
Vowinckel, B, Nikora V, Kempe T, Fröhlich J (2017) Spatially-averaged momentum fluxes and stresses in flows over mobile granular beds: a DNS-based study. J Hydraulic Res 55(2):208–223.View ArticleGoogle Scholar
Yeh, T-H, Cantero M, Cantelli A, Pirmez C, Parker G (2013) Turbidity current with a roof: success and failure of RANS modeling for turbidity currents under strongly stratified conditions. J Geophys Res Earth Surf 118(3):1975–1998.View ArticleGoogle Scholar