Consolidation of freshly deposited cohesive and non-cohesive sediment: particle-resolved simulations

TL;DR

This study uses particle-resolved Navier-Stokes simulations to analyze the consolidation process of freshly deposited cohesive and non-cohesive sediments, providing detailed insights into the stress balance and settling behavior.

Contribution

It introduces a parameterized, first-principles simulation approach that links fluid-particle stress balance to classical effective stress concepts without arbitrary calibration.

Findings

Complete parameterization of the Gibson equation for sediment consolidation

Detailed evaluation of stress terms in fluid-particle mixtures

Comparison of cohesive and non-cohesive sediment settling behaviors

Abstract

We analyze the consolidation of freshly deposited cohesive and non-cohesive sediment by means of particle-resolved direct Navier-Stokes simulations based on the Immersed Boundary Method. The computational model is parameterized by material properties and does not involve any arbitrary calibrations. We obtain the stress balance of the fluid-particle mixture from first principles and link it to the classical effective stress concept. The detailed datasets obtained from our simulations allow us to evaluate all terms of the derived stress balance. We compare the settling of cohesive sediment to its non-cohesive counterpart, which corresponds to the settling of the individual primary particles. The simulation results yield a complete parameterization of the Gibson equation, which has been the method of choice to analyze self-weight consolidation.

Figures

Equations

$$\frac{\partial {\varphi }_v}{\partial t}=\frac{\partial }{\partial y}\left(\frac{k({\rho }_p-{\rho }_f){\varphi }_v^2}{{\rho }_f}+\frac{k{\varphi }_v}{\rho g}\frac{\partial {\sigma }_{\text{eff}}}{\partial y}\right),$$

$$\frac{\partial \text{u}}{\partial t}+\nabla \cdot (\text{u}\text{u})=-\frac{1}{{\rho }_f}\nabla p+{\nu }_f{\nabla }^2\text{u}+{\text{f}}_{\text{IBM}},$$

$$\nabla \cdot \text{u}=0$$

$$m_p\frac{\text{d}{\text{u}}_p}{\text{d}t}=\underset{={\text{F}}_{h,p}}{\underbrace{{\oint }_{{\Gamma }_p}\mathbf{\tau }\cdot \text{n}\text{d}A}}+\underset{={\text{F}}_{g,p}}{\underbrace{V_p({\rho }_p-{\rho }_f)\text{g}}}+{\text{F}}_{c,p},$$

$$I_p\frac{\text{d}{\mathbf{\omega }}_p}{\text{d}t}=\underset{={\text{T}}_{h,p}}{\underbrace{{\oint }_{{\Gamma }_p}\text{r}\times (\mathbf{\tau }\cdot \text{n})\text{d}A}}+{\text{T}}_{c,p}.$$

$${\text{F}}_{c,p}=\sum _{q,q\ne p}^{N_p}\left({\text{F}}_{l,pq}+{\text{F}}_{d,pq}+{\text{F}}_{\text{coh},pq}\right)+$$

$${\text{F}}_{l,pw}+{\text{F}}_{d,pw}+{\text{F}}_{\text{coh},pw},$$

$$F_{\text{DLVO}}=\underset{F_{\text{rep}}}{\underbrace{A_R{R}_{\text{eff}}\exp \left(-\frac{{\zeta }_n}{\kappa }\right)}}-\underset{F_{\text{att}}}{\underbrace{\frac{A_H{R}_{\text{eff}}}{12({\zeta }_n^2+{\zeta }_0^2)}}},$$

$${\text{F}}_{\text{coh}}=\{\begin{array}{ll}-k_{\text{coh}}({\zeta }_n^2-{\zeta }_n\lambda )\text{n}& 0<{\zeta }_n\le \lambda \\ 0& \text{otherwise},\end{array}$$

$${\text{F}}_{\text{coh}}=\{\begin{array}{ll}-\text{Co}\frac{8R_{\text{eff}}}{{\lambda }^2}({\zeta }_n^2-{\zeta }_n\lambda )\text{n}& 0<{\zeta }_n\le \lambda \\ 0& \text{otherwise}.\end{array}$$

$$\text{Co}=\frac{\text{max}(||{\text{F}}_{\text{coh},50}||)}{m_{50}{g}^{\prime }}.$$

$$\underset{{\Omega }_{CV}^{+}}{\int }{\rho }_f\frac{\partial \text{u}}{\partial t}\text{d}V+\underset{{\Gamma }_{CV}^{+}}{\int }{\rho }_f(\text{u}\text{u})\cdot {\text{n}}^{+}\text{d}A=\underset{{\Gamma }_{CV}^{+}}{\int }{\mathbf{\tau }}^{+}\cdot {\text{n}}^{+}\text{d}A,$$

$$\underset{L_{CV}}{\int }{\text{f}}_{\text{IBM}}\text{d}V=\underset{{\Gamma }_{CV}^p}{\int }{\mathbf{\tau }}^{+}\cdot {\text{n}}^{+}\text{d}A+\underset{{\Gamma }_{CV}^p}{\int }{\mathbf{\tau }}^{-}\cdot {\text{n}}^{-}\text{d}A,$$

$$\underset{{\Gamma }_y^{-}}{\int }{\rho }_f(\text{u}\text{u})\cdot {\text{n}}^{-}\text{d}A=\underset{{\Gamma }_{CV}^p}{\int }{\mathbf{\tau }}^{-}\cdot {\text{n}}^{-}\text{d}A+\underset{{\Gamma }_y^{-}}{\int }{\mathbf{\tau }}^{-}\cdot {\text{n}}^{-}\text{d}A.$$

$$\underset{\text{Acceleration term}}{\underbrace{\underset{{\Omega }_{CV}^{+}}{\int }{\rho }_f\frac{\partial \text{u}}{\partial t}\text{d}V}}\underset{\text{External force}}{\underbrace{\underset{{\Gamma }_w}{\int }{\mathbf{\tau }}^{+}\cdot {\text{n}}^{+}\text{d}A+\underset{{\Omega }_{CV}}{\int }{\text{f}}_b\text{d}V}}=$$

$$\underset{\text{Fluid force}}{\underbrace{-\underset{{\Gamma }_y}{\int }\mathbf{\tau }\cdot \text{n}\text{d}A+\underset{{\Gamma }_y}{\int }{\rho }_f(\text{u}\text{u})\cdot \text{n}\text{d}A}}\underset{\text{Particle force}}{\underbrace{-\underset{L_{CV}}{\int }{\text{f}}_{\text{IBM}}\text{d}V}},$$

$$⟨\theta {|}_y⟩=\frac{1}{L_x{L}_z}{\int }_0^{L_z}{\int }_0^{L_x}\theta (x,y,z,t)\text{d}x\text{d}z,$$

$$\underset{\text{Local acceleration}}{\underbrace{\frac{1}{A_w}\int {\text{displaylimits}}_{{\Omega }_{CV}}{\rho }_f\frac{\partial v}{\partial t}\text{d}V}}+\underset{\text{Stress on the top wall}}{\underbrace{⟨p{|}_{L_y}⟩+2{\nu }_f{\rho }_f⟨{\frac{\partial v}{\partial y}|}_{L_y}⟩}}=$$

$$\underset{\text{Fluid stress}}{\underbrace{-⟨p{|}_y⟩+2{\nu }_f{\rho }_f⟨{\frac{\partial u}{\partial y}|}_y⟩-{\rho }_f⟨vv{|}_y⟩}}$$

$$\underset{\text{Particle stress}}{\underbrace{-{\int }_y^{L_y}⟨f_{IBM,x}⟩\text{d}y}},$$

$${\sigma }_{t,f}={\sigma }_{h,f}-{\sigma }_{w,f}+{\sigma }_{i,f}.$$

$${\sigma }_{t,f}(y,t)=\frac{1}{A_w}\int {\text{displaylimits}}_{{\Omega }_{CV}}{\rho }_f\frac{\partial v}{\partial t}\text{d}V$$

$${\sigma }_{w,f}(t)=-⟨p_w{|}_{L_y}⟩+2{\rho }_f{\nu }_f⟨{\frac{\partial v}{\partial y}|}_{L_y}⟩$$

$${\sigma }_{h,f}(y,t)=-⟨p{|}_y⟩+2{\nu }_f{\rho }_f⟨{\frac{\partial v}{\partial y}|}_y⟩-{\rho }_f⟨vv{|}_y⟩$$

$${\sigma }_{i,f}(y,t)=-\frac{1}{A_w}\int {\text{displaylimits}}_{S_{int}}\mathbf{\tau }\cdot \text{n}\text{d}A$$

$$\underset{=⟨{\sigma }_{t,p}{⟩}_y}{\underbrace{\frac{1}{A_w}\sum _{N_{CV}}{m}_p\frac{\text{d}{v}_p}{\text{d}t}}}=\underset{=⟨{\sigma }_{i,p}{⟩}_y=⟨p{|}_y⟩}{\underbrace{\frac{1}{A_w}\sum _{N_{CV}}{\oint }_{{\Gamma }_p}{\tau }_{yy}\cdot n_y\text{d}A}}+$$

$$\underset{=⟨{\sigma }_{\text{tot}}{⟩}_y}{\underbrace{\frac{1}{A_w}\sum _{N_{CV}}{V}_p({\rho }_p-{\rho }_f)g_y}}+\underset{=⟨{\sigma }_{\text{eff}}{⟩}_y}{\underbrace{\frac{1}{A_w}\sum _{N_{CV}}{F}_{c,p,y}}},$$

$$-⟨{\sigma }_{\text{tot}}{⟩}_y=⟨p{|}_y⟩+⟨{\sigma }_{\text{eff}}{⟩}_y,$$

$$⟨{\sigma }_{\text{tot}}{⟩}_y=({\rho }_p-{\rho }_f)g{\int }_y^{L_y}{\varphi }_v(y,t)\text{d}y,$$

$$⟨p{|}_y⟩=\frac{1}{A_w}{\int }_0^{L_z}{\int }_0^{L_x}p(x,y,z,t)\text{d}x\text{d}z.$$

$$⟨{\sigma }_{\text{eff}}{⟩}_y=\frac{1}{A_w}{\int }_0^{L_z}{\int }_0^{L_x}\frac{\Delta V_p}{V_p}{F}_{c,p,y}(x,y,z,t)\text{d}x\text{d}z,$$

$$(1-{\varphi }_v)(v_f-v_p)=-k\frac{1}{g{\rho }_f}\frac{\partial ⟨p{|}_y⟩}{\partial y},$$

Full text

††preprint: APS/123-QED

I Introduction

Fine grained sediments interact via attractive electric forces, commonly referred to as van der Waals (vdW) forces [1, 2], and adhesive forces due to extracellular polymeric substances (EPS) such as biofilms [3, 4]. Cohesive sediment thus behaves very differently from its cohesionless granular counterpart. This is especially true for sediment mixtures containing large quantities of cohesive sediment and organic matter, which are also known as ‘mud’ [5, 6, 7]. Mud deposition can lead to siltation of marine and riverine infrastructure or it can bind contaminants, which has important implications for ecology, sedimentology, and civil engineering as it is ubiquitous in various aquatic environments such as lakes, estuaries, and benthic habitats [8, 9]. The consolidation of mud is also important in the context of deep-sea hydrocarbon exploration [10].

The deposition of mud can be subdivided into two processes: hindered settling [11, 12, 13, 14, 15, 16] and consolidation [17, 18, 19, 20, 21, 22]. Both processes can happen simultaneously in the water column as sediment is still in the process of settling until all suspended grains have made contact with the sediment bed that is supported by a bottom wall [23, 24, 14]. Consolidation is characterized by a contracting sediment bed due to the weight of the overlying deposits. This yields an excess pore pressure, which imposes an upward counterflow through the porous bed that is governed by the bed’s permeability and the effective stress of the sediment. These processes have been the basis of the Gibson equation [25]

[TABLE]

which is derived from one-dimensional mass and momentum conservation principles for the fluid and the sediment to predict the change of the horizontally averaged volume fraction $\phi_{v}$ of consolidating soils over time [24]. Here, $t$ is time, $y$ is the vertical coordinate, $\rho_{p}$ and $\rho_{f}$ are the particle and fluid densities, $\rho=\phi_{p}\rho_{p}+(1-\phi_{p})\rho_{f}$ is the averaged density of the suspension in a given control volume, $g$ is the gravitational acceleration, $k$ is the permeability of the porous medium in units of length per time, and $\sigma_{\text{eff}}$ is the effective normal stress component in the $y$-direction.

The effective stress has classically been assessed by means of the ‘effective stress concept’ [26], which states that the total stress is the sum of the pore pressure and the effective stress. While it is a common definition that the effective stress reflects the part of the sediment weight that is supported by inter-particle contact, the experimental assessment of this physical quantity has been subject to debate. Winterwerp and Van Kesteren [27] argued that this term is used as a mathematical concept to close the stress balance and Sills [19] even claimed that there is no physical meaning to effective stress. Hence, while simulations using the Gibson equation [28, 29] show excellent agreement with experimental data, the model lacks general predictive capabilities since the empirical fitting of the parameters entering the Gibson equation has to be performed separately for each new experiment.

There are several reasons for the limitations of the Gibson model. First, the effective stress concept remains one-dimensional, whereas it has been acknowledged that viscous fingering due to spatially heterogeneous particle concentrations can trigger instabilities that lead to complex three-dimensional flow features [30, 31]. Furthermore, it was pointed out that every experimental setup inevitably contains sidewalls [19]. These sidewalls can lead to preferential drainage that introduces artifacts in the analysis. In recent years, experimental efforts have focused on the settling of sand-mud mixtures [32, 33, 34, 35, 36, 37], where particles are too large to exhibit Brownian motion. Since the Gibson theory was derived for pure muds in the first place, applying the model to sand-mud mixtures might be problematic [38]. At the same time, the larger sized particles render the physical configuration attractive for phase-resolved simulations that allow for a full description of the fluid and the particle motion [39, 40, 41, 42, 43, 44].

The present study addresses these issues from a computational perspective. Following our earlier particle-resolved simulations [45, 46, 47], we present simulation data of settling cohesive and noncohesive sediment where we compute the motion of every particle in a fully resolved three-dimensional flow field without side walls. We present a detailed stress balance for the fluid-particle mixture that allows for a direct transfer of the governing equations to the classical ‘effective stress concept’. As a result, we provide a way to parameterize the Gibson equation in a straightforward fashion and illustrate the effects of cohesive forces on macroscopic particles by means of intergranular stresses in the sediment packing.

The paper is structured as follows: First, we briefly review our fully coupled computational approach to simulate particles in a viscous flow in section II. Then, we derive the ’effective stress concept’ from our governing equations in section III, before we present results for the settling behavior of cohesive grains of silt size and its noncohesive counterpart in section IV.

II Computational Approach

II.1 Particle-resolved Simulations

We solve the unsteady Navier-Stokes equations for an incompressible Newtonian fluid, given by

[TABLE]

along with the continuity equation

[TABLE]

on a uniform rectangular grid with grid cell size $\Delta x=\Delta y=\Delta z=h$ following the scheme of [43]. Here, bold and italic symbols represent vectors and scalar quantities, respectively, where $\textbf{u}=(u,v,w)^{T}$ designates the fluid velocity vector in Cartesian components, $p$ denotes pressure with the hydrostatic component subtracted out, $\nu_{f}$ is the kinematic viscosity, and $\textbf{f}_{\textit{IBM}}$ represents a volume force introduced by the Immersed Boundary Method (IBM) [48, 49]. This volume force acts in the vicinity of the inter-phase boundaries and couples the fluid phase to the particle motion.

Within the framework of the IBM, we consider spherical particles as an approximation of primary particles from the grain-size fraction of silt. We calculate the motion of each individual spherical particle by solving an ordinary differential equation for its translational velocity $\textbf{u}_{p}=(u_{p},v_{p},w_{p})^{T}$

[TABLE]

and its angular velocity $\bm{\omega}_{p}=(\omega_{p,x},\omega_{p,y},\omega_{p,z})^{T}$

[TABLE]

Here, $m_{p}$ is the particle mass, $\Gamma_{p}$ the fluid-particle interface, $\bm{\tau}$ the hydrodynamic stress tensor, $\rho_{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 $\Gamma_{p}$, $\textbf{r}=\textbf{x}-\textbf{x}_{p}$ is the position vector of the surface point with respect to the center of mass $\textbf{x}_{p}$ of a particle, and $\textbf{F}_{c,p}$ and $\textbf{T}_{c,p}$ are the force and torque due to particle interactions, respectively, which are computed using the Discrete Element Method (DEM) as described by [43]. Furthermore, note the designation of the hydrodynamic force and torque as $\textbf{F}_{h,p}$ and $\textbf{T}_{h,p}$, respectively, as well as $\textbf{F}_{g,p}$ the force due to gravity.

We employ the approach of [49] for evaluating the IBM forces and solve (4) and (5) according to [43]. The particles are explicitly coupled to the fluid motion through the hydrodynamic stress tensor $\bm{\tau}$ comprising viscous and pressure drag as a direct result of the IBM. The integration scheme subdivides the fluid time step into a total of 15 substeps to integrate (4) and (5) in time. It was shown by Biegert et al. [43] that this is necessary to resolve short-range particle-particle interactions such as lubrication and cohesive forces.

Our simulation approach was validated by [43] for the fluid-particle coupling of the method against experimental data of a sphere settling in an unbounded quiescent fluid [50] as well as towards a wall [51]. The particle-contact model was validated against benchmark data of [52] and [53]. The collective motion of a sediment bed sheared by a viscous flow was compared to the experimental data of [54] and we found satisfactory agreement.

II.2 Particle-particle Interaction

We use the computational approach of Biegert et al. [43] for modeling cohesionless particle-particle interactions. This reference provides validation results for various benchmark experiments. In addition, the present study employs the cohesive force model proposed and validated by Vowinckel et al. [45]. The particle-particle interaction comprises short-range effects due to unresolved hydrodynamic lubrication forces $\textbf{F}_{l}$ and cohesive forces $\textbf{F}_{\text{coh}}$, as well as direct contact forces $\textbf{F}_{d}=\textbf{F}_{n}+\textbf{F}_{t}$ acting in the normal and tangential directions, denoted as $\textbf{F}_{n}$ and $\textbf{F}_{t}$. The resulting force on particle $p$ is the sum of all these effects

[TABLE]

where the subscripts $pq$ and $pw$ indicate interactions with particle $q$ or a wall, respectively. Detailed information about how we model the lubrication and direct contact forces is provided in [43]. For the present study, we have chosen the same parameterization determined for silicate materials in experiments and used in our previous simulations [43, 55, 56, 39, 41, 42, 45].

To account for cohesive forces, we take an approach that is consistent with the theory for colloids developed by Derjaguin-Landau-Verwey-Overbeek (DLVO) [57, 58], which states that there are two dominant short-range forces that can be interpreted as opposing potentials surrounding particles with grain sizes in the micro- to nanometer range. On the one hand, there exists a repulsive force when equally charged surfaces are in close proximity. On the other hand, as one particle causes correlations in the fluctuating polarization of a nearby particle surface an attractive force is generated. The former effect is usually called the repulsive ‘double-layer’ (DL) force, while the latter effect is commonly referred to as van-der-Waals (vdW) force. These forces become important for gap sizes $\zeta_{0}<\zeta_{n}<\zeta_{\infty}$, where $\zeta_{0}$ defines the microscopic size of surface asperities and $\zeta_{\infty}$ is the distance for which these forces decay to zero [59]. The repulsive DL force and the attractive vdW force due to polarization scale as $F_{\text{rep}}\propto e^{-\zeta_{n}}$ and $F_{\text{att}}\propto\zeta_{n}^{-2}$, respectively. The superposition of the two potentials yields a net force as a function of the gap size $\zeta_{n}$. A model incorporating both effects has been proposed by Pednekar et al. [60]

[TABLE]

where $A_{R}$ is a repulsive force scale for $\zeta_{n}=0$ as a measure of the particles’ surface potential, $\kappa$ is the Debye length, $A_{H}$ is the Hamaker constant, and $\zeta_{0}$ is the surface roughness preventing $F_{\text{att}}$ from diverging to infinity for vanishing gap size. These four parameters need to be adjusted according to the physical system.

To account for these effects in particle resolved simulations, Vowinckel et al. [45] proposed a model with the following properties: (i) it decays to zero as the gap size goes to zero, (ii) it has a maximum at a gap width orders of magnitude smaller than the particle diameter, and (iii) it decays to zero for larger gap sizes, without any discontinuous jumps. The physical idea of this model is described in [45]. The above properties are fulfilled with the mathematically simple model of a parabolic spring force

[TABLE]

where $k_{\text{coh}}=\frac{A_{H}R_{\text{eff}}}{\zeta_{0}\lambda^{3}}$ denotes the stiffness constant, $\lambda$ represents the range over which the cohesive force is distributed, $A_{H}$ is the Hamaker constant, $\zeta_{0}$ is the minimal separation distance [61] and $R_{\text{eff}}=R_{p}R_{q}/(R_{p}+R_{q})$ is the effective radius. The cohesive range $\lambda$ can be interpreted as a Debye length.

The dimensional form (8) still requires the proper parameterization of the empirical parameters $A_{H}$ and $\zeta_{0}$. However, we can replace these empirical constants by writing (8) with respect to the maximum cohesive force $\text{max}(||\textbf{F}_{\text{coh},50}||)$. Choosing the median grain diameter $D_{50}$, the buoyancy velocity $u_{s}=\sqrt{g^{\prime}D_{50}}$, the characteristic time scale $\tau_{s}=D_{50}/u_{s}$ and the characteristic mass $m_{50}=\rho_{f}\pi D_{50}^{3}/6$, the characteristic force scale for particles settling under gravity in an otherwise quiescent fluid becomes the specific weight $m_{50}g^{\prime}$, where $g^{\prime}=(\rho_{p}-\rho_{f})g/\rho_{f}$ denotes the reduced gravity. After normalizing (8) with the specific weight we obtain

[TABLE]

Hence, the characteristic parameter to define cohesive forces becomes the cohesive number

[TABLE]

It represents the ratio of the maximum cohesive force for particles of diameter $D_{50}$ to the characteristic gravitational force scale of the problem [62].

We can transfer the DLVO-theory (7) to our simpler model (9) by choosing the following parameters: (i) $A_{H}=1\cdot 10^{-20}$ J, which reflects silica materials in water according to [64], (ii) $R_{\text{eff}}=\frac{R_{p}R_{q}}{R_{p}+R_{q}}=\frac{R_{p}}{2}=5\mu$m for monodisperse silt particles of grain size $D_{p}=20\mu$m, (iii) $\zeta_{0}=4.7$nm; (iv) we determine $\kappa$ using the approximation for the monovalent salt sodium chloride given by [65] as $\kappa=\frac{0.304\cdot 10^{-9}\text{m}^{-1}}{|z|\sqrt{C_{\text{salt}}}}$ in meters, where $z$ is the valency of the salt and $C_{\text{salt}}$ is the salt concentration in mol/liter. Here, we choose the salinity of sea water with 35ppt. These parameters yield $C_{\text{salt}}=0.6$mol/liter and $\kappa=0.393$nm; (v) Since a key feature of our model is to have vanishing forces for particle contact, i.e. $\zeta_{n}=0$, we set $A_{R}=A_{H}/(12\zeta_{0}^{2})$. The DLVO curve for this case is displayed in figure 1. In this scenario, the total force $F_{\text{DLVO}}$ follows the attractive forces with a distinct minimum at $\zeta_{n}\approx 1$nm. The minimum force is $|\min(F_{\text{DLVO}})|=3\cdot 10^{-10}$N, while the weight becomes $F_{g}=\pi\,g(\rho_{p}-\rho_{f})D_{p}^{3}/6=6.78\cdot 10^{-11}$N, where we set gravitational acceleration, particle density, and fluid density to be $g=9.81\text{m/s}^{2}$, $\rho_{p}=2650\text{kg/m}^{3}$, and $\rho_{f}=1000\text{kg/m}^{3}$, respectively. This yields a cohesive number of $\text{Co}=|\min(F_{\text{DLVO}})|/F_{g}=5.0$. We found our parabolic spring model to be a good approximation of the curve shown by the solid line in figure 1.

Consequently, the cohesive number used in our simulations corresponds to the properties of fine to medium sized grains of silt settling in salt water. Note that the present approach can easily be extended to different types of cohesive sediment given that the temporal discretization is able to resolve the parabolic spring model during particle-particle interaction. In such case, the parameter needed to determine the cohesive number is $|\min(F_{\text{DLVO}})|$, i.e. the maximum attractive force. To investigate biofilms, for example, one can no longer use the present analogy of the DLVO-theory. Instead, the maximum attractive force would have to be determined experimentally [63].

II.3 Simulation Scenario

To explore the influence of cohesive forces on the sedimentation process of a large, polydisperse ensemble of particles, we further analyze the simulation data first presented in Vowinckel et al. [45], who investigated the hindered settling of silt particles. In this reference, a polydisperse mixture of particles with a relative density of $\rho_{p}/\rho_{f}=2.6$ was placed in a tank of viscous fluid to obtain an initial volume fraction of $\phi_{v}=V_{s}/V_{0}=0.155$ (Figure 2a), where $V_{s}$ denotes the volume occupied by the particles and $V_{0}=L_{x}\times L_{y}\times L_{z}=13.1D_{50}\times 40.0D_{50}\times 13.1D_{50}$ is the computational domain size.

Consistent with the experiments of [37], we chose a Reynolds number of $\text{Re}=u_{s}D_{50}/\nu_{f}=1.35$. The size distribution of the polydisperse particles was continuous with a ratio of ${\max\{D\}/\min\{D\}=4}$ obeying a log-normal size distribution, which yields $\min\{D\}/D_{50}=0.6$ and $\max\{D\}/D_{50}=2.4$ for the smallest and largest particle, respectively. A total of 1261 particles was placed randomly to obtain an almost uniform initial profile of $\phi_{v}$ (Figure 2b). We conducted preliminary tests of dry settling, i.e. neglecting fluid forces, with varying particle distributions and domain sizes to investigate the dependency of the final deposit on the initial particle distribution. With the present domain being $40D_{50}$ in height, we found no dependence of the final configuration on the initial particle distribution, either. We impose a no-slip condition at the bottom wall ($y=0$) and at the particle surfaces as well as a free-slip condition at the top wall ($y=L_{y}$), along with periodic boundary conditions in the wall-parallel $x$- and $z$-directions. This means that sidewalls are absent in our numerical simulations. Gravity is pointing towards the bottom wall in the negative $y$-direction. The minimum, median and maximum primary particle sizes are discretized by $\min\{D\}/h=11$, $D_{50}/h=18.25$ and $\max\{D\}/h=44$ grid cells, respectively. It was shown by Vowinckel et al. [45] by a comparison to the analytical particle settling velocity that the grid resolution of our computational domain is fine enough to capture the settling behavior of all particle sizes at the particle Reynolds numbers encountered in the simulations presented in section IV. We spread the cohesive forces over a shell of thickness $\lambda=h$, but it was shown by Vowinckel et al. [45] that the simulation results are not sensitive to this choice provided that $\lambda<R_{p}$.

Two simulations were performed for different values of the cohesive number: (i) cohesionless grains with $\text{Co}=0$, and (ii) cohesive sediment with $\text{Co}=5$. For both simulations, the particles with identical initial particle distributions were released from rest in a quiescent fluid to guarantee a straightforward comparison of (i) and (ii). Subsequently, the particles settle under the influence of gravity undergoing different settling behavior due to cohesive forces. Following the scaling argument of [60] presented in section IIB above, Co=5 corresponds to the properties of fine to medium sized silt in a saline ambient, i.e. the settling of macroscopic silica particles in ocean water. For example, choosing a median grain size of $D_{50}=25\mu$m renders the simulation domain 1mm tall. As will be shown in section IV below, this choice allows for a simulation domain that is large enough to capture all relevant processes of particle settling so that the results can be evaluated along the lines of the ‘effective stress concept’.

III Stress Balance for the Fluid-particle Mixture

To understand the settling and the consolidation of the fluid-particle mixture, we analyze the balance of the wall-normal stress components for the two phases separately. According to Biegert et al. [46] and Biegert [47], we can write the momentum balance of the fluid (2) in an integral sense to obtain fluid stresses for a control volume $\Omega_{CV}$ that extends from the top-wall ($y=L_{y}$) to an arbitrary height $y$ in the vertical direction and encompasses the entire domain in the $x$- and $z$-directions (figure 3). We can write the integral form of (2) as

[TABLE]

where $\Gamma_{\mathit{CV}}^{+}=\Gamma_{w}\cup\Gamma_{s}\cup\Gamma_{y}^{+}\cup\Gamma_{\mathit{CV}}^{p}$ comprises all surfaces of the control volume shown for a single particle as a minimal example in Figure 3b and $\textbf{n}^{+}$ is the normal vector pointing outwards from $\Omega_{\mathit{CV}}^{+}$. In addition, we recast the pressure and viscous terms using the fluid stress tensor, $\bm{\tau}^{+}=-p\textbf{I}+\mu_{f}(\nabla\textbf{u}+(\nabla\textbf{u})^{T})$, where I is the identity matrix. Here we neglect the effect of the immersed boundary force, assuming that it is implicitly handled by the fluid stress at the fluid/particle interface.

All of these terms except for the fluid stress at the particle surface are straightforward to calculate. However, we can evaluate the fluid stress indirectly using the IBM force, as was done to obtain the particle equations of motion (4) and (5). That is, the IBM force acts as a jump in stress:

[TABLE]

where we are careful to distinguish between $\textbf{n}^{+}$, the outward surface normal for the volume $\Omega_{\mathit{CV}}^{+}$, and $\textbf{n}^{-}$, the outward surface normal for the volume $\Omega_{\mathit{CV}}^{-}$, which point in opposite directions.

To evaluate the fluid stress $\bm{\tau}^{-}\cdot\textbf{n}^{-}$ in the particle interior, we can perform a stress balance on $\Omega_{\mathit{CV}}^{-}$ (figure 3a). The integral form of the Navier-Stokes equations together with divergence theorem give us

[TABLE]

Using (13) together with (12) and (11), we obtain

[TABLE]

where $\Omega_{\mathit{CV}}=\Omega_{\mathit{CV}}^{+}\cup\Omega_{\mathit{CV}}^{-}$, $\Gamma_{w}$ is the area of the top wall, and $\Gamma_{y}=\Gamma_{y}^{+}\cup\Gamma_{y}^{-}$.

Using the definition of the horizontal, superficial average [66],

[TABLE]

we can rewrite (III) for the $y$-velocity component as

[TABLE]

where $A_{w}=L_{x}L_{z}$ is the area of the bottom wall. For this derivation, we have used the fact that $\mu_{f}$ and $\rho_{f}$ are constant throughout the domain.

This yields the stress balance of the fluid phase for the horizontally averaged, vertical component of the normal stress

[TABLE]

Here,

[TABLE]

is the stress due to local acceleration. Further,

[TABLE]

is the hydrodynamic stress on the top wall. The angular brackets denote horizontal averaging. The hydrodynamic stress due to fluid motion becomes

[TABLE]

which comprises pressure and viscous forces as well as convection. Finally,

[TABLE]

is the interfacial stress exerted on the fluid by the particles, where $S_{int}$ is the total area of the fluid particle interface enclosed in the control volume $\Omega_{CV}$. In our simulations, we found that $\sigma_{w,f}-\sigma_{h,f}=\sigma_{i,f}$ holds at all $y$-locations and all times (e.g. figure 2c). Furthermore, the fluid stress was dominated by the pressure term, so that we conclude $\sigma_{w,f}-\sigma_{h,f}\approx\langle p\lvert_{y}\rangle$, provided that viscous stresses and stresses due to convection are small. Since the hydrostatic component is subtracted out for (2), $p$ is equivalent to the excess pressure in the water column [27].

For the particle phase, we sum over all particles within the control volume for the different terms in the particle equation of motion (4) to get integral quantities:

[TABLE]

where $N_{CV}$ is the number of particles enclosed in $\Omega_{CV}$ and the operator $\langle\cdot\rangle_{y}$ indicates the volume average over $\Omega_{CV}$, which is to be distinguished from the horizontal averaging operator $\langle\cdot\rangle$ defined by equation (15).

The explicit link between the two phases via the IBM yields that $\sigma_{i,f}$ is equal to the static excess pressure (cf. figure 2c). Hence, (17) and (III) are coupled through the equality $-\sigma_{i,f}=\langle p\rvert_{y}\rangle=\langle\sigma_{i,p}\rangle_{y}$. Using this equality and assuming $\langle\sigma_{t,p}\rangle_{y}$ to be small, we can interpret (III) in terms of the ’effective stress concept’ [26]. It was described in this reference that

[TABLE]

where $-\langle\sigma_{\text{tot}}\rangle_{y}$ is the weight of the particles submerged within $\Omega_{CV}$, $\langle p\rvert_{y}\rangle$ is the pore water pressure at location $y$ and $\langle\sigma_{\text{eff}}\rangle_{y}$ is the effective stress due to particle interactions. For this reason, we normalize all stresses by the total submerged weight of all the particles $\langle\sigma_{\text{tot}}\rangle_{0}$ in the following section.

The weight of the sediment can be computed in a straightforward manner:

[TABLE]

where $\phi_{v}$ is particle volume fraction, and $L_{y}$ is the length of the domain in $y$-direction. The pressure is computed as the horizontal average in plane $y$, which is the lower bound of the control volume under consideration:

[TABLE]

To compute depth-resolved distributions of effective stresses, we introduce a horizontal averaging operator that uses a step function to distinguish between particles and fluid:

[TABLE]

where $A_{w}$ is the area of the bottom wall and $\Delta V_{p}$ is the volume of particle $p$ cut by the horizontal slice, and $V_{p}$ is the total volume of this particle. This operator is very similar to the one used by [45] to compute averages of particle velocities.

IV Results

IV.1 Settling Behavior

After the initial release of the particles from rest, the grains start to settle towards the bottom. As shown by Vowinckel et al. [45], the kinetic energy of the particles peaks at $t=17.6\tau_{s}$. The stress balance for this moment is shown in Figure 2c. Since both the cohesionless and the cohesive simulation are initialized identically, their results collapse for this initial stir-up phase. During this phase, particles have not yet settled out and all the static pressure induced by the particle weight is transferred to the fluid pressure. As the particle distribution is still fairly uniform (figures 2a and b), the two profiles shown in figure 2c are approximately linear in the $y$-direction. The excess pore pressure builds up throughout the water column to reach a balance with $\langle\sigma_{tot}\rangle_{0}$ at the bottom wall.

Figure 4 shows the settling process over time as snapshots taken at $t/\tau_{s}=25\%\,T_{\text{sim}}$, $t/\tau_{s}=50\%\,T_{\text{sim}}$ and $t/\tau_{s}=100\%\,T_{\text{sim}}$, where $T_{\text{sim}}$ is the total simulation time. The figures show both, the instantaneous particle distribution and a translucent contour slice cutting through the front of the domain showing the vertical fluid velocity component. Particles are also colored by their wall-normal velocity component following the color scheme of the contour slices. Since particles appear as opaque objects, the grains in the front are blocking the view on the particles in the back of the domain. Note, however, that since we are resolving the motion of each particle individually, the number of particles remains constant throughout the simulations.

As particles start to settle, they replace fluid at the bottom of the tank and generate an upward counterflow (red regions in figure 4 at $t=120\tau_{s}$). For the current simulations, this counterflow is sufficiently strong to sweep smaller particles upward. This effect is more pronounced for cohesionless particles. The upwelling fluid represents one mechanism for hindered settling and for the segregation of the grain sizes for very large water columns [37]. Indeed, it was shown by Vowinckel et al. [45] that the present data yields excellent agreement with the classical hindered settling functions of Richardson and Zaki [12] and Winterwerp [14].

As time progresses the impact of cohesive forces on the settling behavior becomes evident (figure 4 at $t=240\tau_{s}$). The cohesive sediment starts to built aggregates as particles attach to each other once they come into close enough contact. The flocculation process is exemplified in figure 5, which shows a zoom into figure 4b at $t=240\tau_{s}$. For the current physical configuration, cohesive particles are able to form chains of several particles with varying diameter. Hence, small cohesive particles that bond to bigger ones move with the settling velocity of the large particle. As a result, smaller particles settle faster than individual primary particles would. This observation is in line with experimental evidence [67, 68, 69].

This observation is addressed in a statistical sense in figure 6. Figure 6 shows the number of settling particles $n_{s}$ (conditioned by $v_{p}>0.01u_{s}$). Initially, particles start to accelerate and those that were initially close to the bottom are immediately filtered out. Later in time, a steady decrease of $n_{s}$ over time can be observed for the cohesionless sediment, whereas for cohesive sediment $n_{s}$ shows a slight increase reaching a local maximum at $t\approx 100\tau_{s}$. During this time, particles are moving rather fast through the domain, and once they get into contact, they form larger aggregates such as the one depicted in figure 5. Subsequently, the number of settling particles decreases more rapidly for cohesive sediments as the aggregates settle faster than the noncohesive sediments.

The aggregation process is further addressed by counting all particle-particle interactions of all settling particles. Throughout the entire simulation time, the average number of direct contacts $n_{\text{con}}$ and the average number of short-range interactions $n_{coh}$ (within the distance $\lambda$) remains fairly constant for cohesionless grains. Note that for this type of sediment, short-range interactions are influenced by lubrication forces only. Cohesive sediment, on the other hand, forms larger aggregates with up to three different particles on average considering direct contacts and cohesive short-range interaction combined. As larger aggregates settle faster than individual particles, these aggregates make contact with the bottom wall earlier in time, so that the number of particle-particle interactions summed over all settling particles starts to decay at $t\approx 300\tau_{s}$. As the cohesive sediment settles out faster than the noncohesive sediment, the counterflow decays earlier in time (figure 4b $t=480\tau_{s}$).

IV.2 Transition from Hindered Settling to Consolidation

The enhanced settling speed for $\text{Co}=5$ is confirmed by the horizontally averaged concentration profiles plotted over time (Figure 7a and b) following the analysis of Been and Sills [17]. Up to $t=120\tau_{s}$, the volume fraction contours for the two simulations are nearly identical, as cohesive forces have not yet had sufficient time to cause a noticeable change. As time progresses, two fronts become clearly visible for both simulations. The upper front marks the transition between the clear fluid and the suspended sediment, whereas the second front shows the transition between the suspended sediment and the sediment bed. The isoline marking this front is called the gelling concentration [14]. Since cohesive forces result in the formation of flocs with larger settling speeds, particles accumulate at the bottom of the tank more quickly. At $t\approx 350\tau_{s}$, the two fronts merge into one for the simulation data of the cohesive sediment. This point in time is called the point of contraction (PoC) and marks the transition between hindered settling and consolidation [27]. The cohesionless sediment, on the other hand, has not yet reached the PoC by the end of the simulation.

The transition from hindered settling to consolidation can also be analyzed in terms of the effective stress concept (23). To this end, we compute all stress components for $\Omega_{CV}=V_{0}$, i.e. the entire computational domain. To evaluate the effective stress, we compute the normal stresses due to contact forces in the wall-normal direction. Note that we also computed the normal stress components in the periodic $x$- and $z$-directions as will be shown in section IV D below. We found these two wall-parallel components to be identical, which yields transversal isotropic conditions. This proves that the width of the simulation domain is large enough to capture all relevant effects of particle-particle interactions in these directions.

We plot all components entering the effective stress balance (equation (23)) over time in figure 7c. It was shown in section IV.1 that particles are deposited on the bottom more rapidly for cohesive sediment. These dynamics are also reflected in the stress balance of the fluid-particle mixture. Recall that $\langle p\rvert_{0}\rangle$ is the fluid stress acting on the particle phase, which after the initial phase is equivalent to the weight of those particles that are still suspended [17, 24]. On the other hand, $\langle\sigma_{\text{eff}}\rangle_{0}$ captures all particle-particle interactions. Since the forces between two interacting particles in suspension are opposite and equal, these do not cause a net force on the bottom wall. Thus, $\langle\sigma_{\text{eff}}\rangle_{0}$ reflects the weight of those particles that are supported by the external contact forces with the bottom wall. After the initial increase of $\langle p\rvert_{0}\rangle$ up to $\langle\sigma_{\text{tot}}\rangle_{0}$, when all of the particle weight is supported by fluid forces (cf. figure 2c), $\langle p\rvert_{0}\rangle$ decays over time, whereas $\langle\sigma_{\text{eff}}\rangle_{0}$ increases at the same rate. Hence, less and less of the particle weight is supported by fluid forces, and more and more of it is supported by interparticle forces within the sediment bed, which illustrates the transition from hindered settling to consolidation [27]. While $\langle p\rvert_{0}\rangle$ and $\langle\sigma_{\text{eff}}\rangle_{0}$ behave very similarly for cohesive and noncohesive sediment until $t=180\tau_{s}$, cohesive sediment experiences a systematic shift at this time that illustrates the enhanced settling due to flocculation. Particles that assemble in flocs are seen to make contact with the wall at earlier times. Apart from the initial increase, we obtain $\langle p\rvert_{0}\rangle+\langle\sigma_{\text{eff}}\rangle_{0}\approx\langle\sigma_{\text{tot}}\rangle_{0}$, whereas the rate of change for both components is approximately linear. This suggests that the local acceleration does not add to the stress balance of the particle motion even though the processes under investigation are inherently transient. Interestingly, we obtain $\langle\sigma_{\text{eff}}\rangle_{0}=0.71\langle\sigma_{\text{tot}}\rangle_{0}$ for the cohesive sediment at the PoC ($t\approx 350\tau_{s}$). This is the same value observed for cohesionless sediment at the end of the simulation time $t=480\tau_{s}$, where the PoC has not been reached (figure 7a), but the height of the cohesive sediment bed is much thicker compared to the cohesionless sediment bed. Hence, comparing figures 7a and b, we conclude that flocculation promotes larger pore spaces due to different vertical distributions of effective stresses.

While (III) provides a measure for the volume averaged stresses, we can also evaluate the stress terms as vertical profiles for all components of (23). The instants shown in figure 8a correspond directly to the situations in figure 4. As expected, we obtain $\langle p\rvert_{y}\rangle=\langle\sigma_{\text{tot}}\rangle_{y}$ far away from the wall, where all particles are fully supported by the fluid. Initially, $\langle p\rvert_{y}\rangle$ decreases linearly with height ($t=120\tau_{s}$). Whereas $\langle\sigma_{\text{tot}}\rangle_{y}$ decreases linearly in the suspended region, the profile for cohesive sediment becomes convex as flocs start to form, which indicates an accelerated settling process.

The increase of the sediment bed thickness is shown by regions with $\langle\sigma_{\text{eff}}\rangle_{y}>0$. Again, $\langle\sigma_{\text{eff}}\rangle_{y}$ increases linearly with height. In addition, $\langle\sigma_{\text{eff}}\rangle_{y}$ and $\langle p\rvert_{y}\rangle$ add up to $\langle\sigma_{\text{tot}}\rangle_{y}$ by an out-of-balance of less than 1%, illustrating the capability of our method to explicitly compute the stress balance of the entire fluid-particle mixture.

IV.3 Parameterization of the Gibson Model

The Gibson equation (1) was derived by [24] for the simultaneous treatment of sedimentation and self-weight consolidation by using continuity principles and the Darcy-Gersevanov law:

[TABLE]

where $v_{f}$ is the horizontally averaged vertical fluid velocity and $v_{p}$ is the horizontally averaged particle velocity. Solving (1) requires constitutive relationships for the parameters $k$ and $\sigma_{\text{eff}}$ [70]. According to [24], it has been a common assumption that $k$ and $\sigma_{\text{eff}}$ depend on the volume fraction only, even though it was discussed in this reference that the assumption can not hold for $\sigma_{\text{eff}}$, which must be zero in a suspension while $\phi_{v}$ is not.

We can test this hypothesis by computing $k$ and $\sigma_{\text{eff}}$ based on (27) and (III). If $k$ and $\sigma_{\text{eff}}$ were functions of $\phi_{v}$ only, then the data from different instants in time would collapse on a single master curve. This is illustrated in figure 8b for the data plotted in figures 4 and 8a. The curves for the different quantities and different times do not collapse and the lines are highly convoluted in some regions. While this is due to the semi-logarithmic plot style as $k$ and $\sigma_{\text{eff}}$ approach zero, we can clearly identify three characteristic regions for both quantities: (i) a suspension region ($\phi_{v}<0.25$) with a high permeability, which decreases rapidly as $\phi_{v}$ decreases. This behavior is most pronounced at $t/\tau_{s}=120$, when most of the particles are still in suspension. Since there are hardly any particle-particle interactions in this region, $\sigma_{\text{eff}}$ approaches zero. In addition, the lines of $k$ and $\sigma_{\text{eff}}$ are highly convoluted around $\phi_{v}\approx 0.15$. Hence, as expected, the assumption of $k$ and $\sigma_{\text{eff}}$ being a function of the volume fraction only does not hold in this region. (ii) A consolidating region ($0.25<\phi_{v}<0.58$) with an exponential decay of $k$ and increase of $\sigma_{\text{eff}}$, respectively, as $\phi_{v}$ increases. Hence, this region serves very well to derive scaling laws for $k$ and $\sigma_{\text{eff}}$. (iii) A jamming region ($\phi_{v}=0.58$), where large values for $\sigma_{\text{eff}}$ represent the layer formed at the very bottom of the domain that holds the weight of the entire overlying sediment. As the fluid flow ceases and the soil becomes fully consolidated, $k$ approaches zero as well.

As a result, we obtain a mean permeability of $k/u_{s}=5.7\cdot 10^{-4}$ within the sediment bed. We can convert the nondimensional permeability into dimensional quantities choosing $D_{50}=20\mu$m, $\rho_{p}=2600\text{kg/m}^{3}$, $\rho_{f}=1000\text{kg/m}^{3}$, and $g=9.81\text{m/s}^{2}$. This yields $u_{s}=\sqrt{D_{50}g(\rho_{p}-\rho_{f})/\rho_{f}}=1.8\cdot 10^{-2}$m/s and a dimensional value of ${k=1\cdot 10^{-5}}$m/s. This value corresponds very well to the permeability of a semi-pervious medium such as very fine sand or silt [71], which is exactly in line with the chosen values for Co and Re that are meant to represent primary particles with a median grain size of silt.

Le Hir et al. [28] and Grasso et al. [29] have used power laws as constitutive relationships to parametrize $k$ and $\langle\sigma_{\text{eff}}\rangle$. Based on the observation above, we can immediately conclude that fitting a power law function over all three regions will result in a poor fit. Hence, we performed a regression of a power law for regions (ii) and (iii) in figure 8b for all data that are marked by the circles. The fitted functions have different coefficients for cohesive and noncohesive sediments. A fair agreement for region (ii) can be achieved, but region (iii) is not well described for $k$. The coefficient of determination $R^{2}$ ranges between 0.67 and 0.81, which is satisfactory. Hence, the presented data highlights the capability of our simulation approach to improve parameterization strategies for the Gibson equation. However, larger simulations for a wider range of cohesive numbers and larger Reynolds numbers are needed to obtain a better scaling of $k$ and $\sigma_{\text{eff}}$.

IV.4 Interparticle Stress

While the analysis in section IV.3 provides insight into how the support of the granular weight transitions from hydrodynamic forces to collision forces with the bottom wall, it does not characterize internal granular stresses parallel to the wall since sidewalls are absent in our numerical simulations. As a consequence, all interparticle forces are opposite and equal in this direction. On the other hand, the weight of the particles $\textbf{F}_{g}$ does not cancel out for contacts in the wall-normal direction so that contact forces at a given height reflect the weight of the overlying sediment. We applied operator (26) to all contact forces summarized by equation (II.2). We denote $\sigma_{\perp}=\sigma_{yy}$ as the normal stress acting in $y$-direction, whereas $\sigma_{\parallel}=\frac{1}{2}(\sigma_{xx}+\sigma_{zz})$ denotes the normal stresses acting parallel to the wall. Having this analysis in place, we can further subdivide the effective stress into its components comprising stresses due to direct contact $\sigma_{d}$, unresolved lubrication $\sigma_{l}$, and cohesion $\sigma_{\text{coh}}$.

Since we focus on the sediment deposit, figures 9a - c show a close-up of the bottom part of the domain. For cohesionless sediment, the data was taken at the final simulation stage, i.e. $t=480\tau_{s}$. At this stage, $71\%$ of the particle weight is supported by contact forces. For a meaningful comparison, we take the data from the cohesive sediment simulation at $t=346\tau_{s}$, when it has the same value of $\langle\sigma_{\text{eff}}\rangle_{0}$. For both simulations, it becomes immediately obvious that stresses due to direct contact (Figure 8a) are one order of magnitude larger than the other two components (Figure 8b and c). As expected, the nearly linear profiles of $\sigma_{d,\perp}$ reflect effective stress as the weight of the sediment beds. The steeper slope of the curve for cohesionless sediment is consistent with the fact that it is packed much more densely. The horizontal component $\sigma_{d,\parallel}$ shows significant stresses within the sediment deposit that peak at $y/D_{50}\approx 2$. These are more pronounced for cohesive grains. Throughout the sediment column, unresolved lubrication stresses remain small (Figure 8b). This holds for both simulations and is expected as most of the lubrication forces are resolved by the IBM. Figure 8c shows significant cohesive forces throughout the entire sediment deposit. Interestingly, the vertical component $\sigma_{\text{coh},\perp}$ is smaller than the horizontal component $\sigma_{\text{coh},\parallel}$. From a physical point of view, the horizontal intergranular stress due to cohesion $\sigma_{\text{coh}\parallel}$ prevents particles from arranging themselves into the densest packing formation, while the vertical cohesive stress $\sigma_{\text{coh},\perp}$ works towards a more rapid consolidation of the sediment. In the present simulation the horizontal component exceeds the vertical one, and we conclude that even small cohesive forces enable the deposited particles to remain in flocs as they are subjected to the weight of the particles settling from above. This effect accounts for the lower sediment volume fraction at the bottom of the tank for the cohesive sediment simulation (Figure 7).

V Conclusions

The present study successfully applied phase-resolved simulations to the situation of consolidation of freshly deposited cohesive and non-cohesive sediment. To this end, we have derived a stress balance based on the governing equations that makes a direct connection to the effective stress concept. The simulations fully resolve the three-dimensional flow field and provide an efficient way for excluding sidewalls from the analysis by employing periodic boundary conditions. The highly resolved data yields a physical interpretation of the effective stress as the part of the sediment weight that is supported by external contact forces with the bottom wall, as well as a full parametrization of the Gibson equation.

An analysis of the intergranular stresses clarified the respective roles of direct contact, unresolved lubrication, and cohesive forces during consolidation and dewatering of a freshly deposited bed. For the presented data, we find that direct contact forces dominate over the other two components. As a result of the attractive Van der Waals forces, cohesive sediment experiences larger intergranular stresses due to direct contact forces. Within the sediment deposit, cohesive forces yield intergranular stresses that lead to larger pore spaces than obtained for the consolidation of cohesionless sediment.

Hence, the first test case presented here encourages further studies with larger physical domains and a larger number of particles to achieve more realistic settling conditions. For example, larger domains could accommodate three-dimensional instabilities of particle-induced viscous fingering in the flow field [30, 31] and the stress balance presented here could further enhance the development of two-phase flow model closures such as the $\mu(I)$-rheology [73, 74] or kinetic theory [75]. In addition, thicker sediment beds could lead to self-weight consolidation in a creeping motion [19, 72]. To extend the present simulation approach from macroscopic silica grains towards aggregates of mud, a more realistic description of the particles is needed [6]. For example, the physical properties of the spherical particles used in the present study could be extended towards particles that are porous and compressible [76], so that a description of creeping consolidation and dewatering becomes possible.

Acknowledgements.

This research is supported by the National Science Foundation (NSF) through grant CBET-1638156 and through Army Research Office grant W911NF-18-1-0379. BV gratefully acknowledges the Feodor-Lynen scholarship provided by the Alexander von Humboldt Foundation, Germany. The authors thank J. Israelachvili for stimulating discussions on the interaction between colloids. Computational resources for this work used the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by the National Science Foundation, USA, Grant No. TG-CTS150053.