Investigating cohesive sediment dynamics in open waters via grain-resolved simulations

Abstract

Cohesive particulate flows play an important role in environmental fluid dynamics, as well as in a wide variety of civil and process engineering applications. However, the scaling laws, constitutive equations and continuum field descriptions governing such flows are currently less well understood than for their non-cohesive counterparts. Grain-resolved simulations represent an essential tool for addressing this shortcoming, along with theoretical investigations, laboratory experiments and field studies. Here we provide a tutorial introduction to simulations of fine-grained sediments in viscous fluids, along with an overview of some representative insights that this approach has yielded to date. After a brief review of the key physical concepts governing van der Waals forces as the main cohesive effect for subaqueous sediment suspensions, we discuss their incorporation into particle-resolved simulations based on the immersed boundary method. We subsequently describe simulations of cohesive particles in several model turbulent flows, which demonstrate the emergence of a statistical equilibrium between the growth and break-up of aggregates. As a next step, we review the influence of cohesive forces on the settling behaviour of dense suspensions, before moving on to submerged granular collapses. Throughout the article, we highlight open research questions in the field of cohesive particulate flows whose investigation may benefit from grain-resolved simulations.

Impact Statement

The erosion, transport and deposition of particles by a fluid feature prominently in a wide range of environmental situations and engineering applications. The dynamics of such flows are frequently governed by physical mechanisms at the particle scale, which need to be properly accounted for in order to develop accurate continuum-scale descriptions for their behaviour. In particular, the dynamics of small, cohesive particles such as clay, silt and dust can be quite complex, as they give rise to aggregation. Here we review recently developed computational approaches that account for the effects of cohesion at the particle level, with a focus on grain-resolved simulations using the immersed boundary method. We subsequently discuss recent simulation results regarding the settling behaviour of cohesive sediment, its flocculation and break-up in turbulent flows, and the way cohesive granular collapses are distinct from their non-cohesive counterparts. These computational approaches provide a path towards gaining further insight into such disparate phenomena as turbidity currents, pyroclastic and debris flows, dust storms and cohesive powders.

1. Introduction

The erosion, transport and deposition of particles by a fluid play an important role in a variety of environmental situations and engineering applications. Generally speaking, sediment is solid material that can be mobilized, transported and deposited by a fluid flow. It can comprise materials from rocks, but also organic matter originating from plants and animals. Typically, sediment is classified by its grain size, where it is common to distinguish between boulders, cobbles, gravel, sand, silt and clay, in order of decreasing grain size. Sediment transport, therefore, constitutes a prime example of multiphase particulate flows (Vowinckel, 2021).

An interesting distinction arises for the two smallest classes of grains, viz. silts with grain diameters less than 63 $\mathrm {\mu }$m, and clays with grain sizes below 2 $\mathrm {\mu }$m. Compared with larger grains such as sand, the ratio of surface to body forces is higher for such small grains, which has important implications. First, such small particles are easier to mobilize and to maintain in suspension. This is the reason, for example, that dry silt can be transported as a wind-blown aeolian suspension. Similar considerations hold for subaqueous processes where fine-grained sediments tend to be transported in suspension or even as wash load, i.e. the grains rarely come into contact with the coarser sediment bed while being convected along with the flow (cf. table 1). Second, grains of silt and clay size are sufficiently small for cohesive forces to become relevant, as they are no longer obscured by the weight of the particles. These attractive forces can have different origins such as van der Waals forces, liquid bridging and electrostatics. The latter two are more relevant for dust particles in gases or industrial flows. Such attractive interaction causes the particles to bond to each other and form complex structures, such as flocs and aggregates (Maggi, Mietta, & Winterwerp, 2007; Mehta, Hayter, Parker, Krone, & Teeter, 1989; Partheniades, 2009; Winterwerp, 2002, 1998). Different terminologies exist for this process, and it has become common to differentiate between aggregation as the formation of more densely packed structures, and flocculation, which refers to much looser structures that are more prone to breaking apart again (Berg, 2010). Figure 1 provides an overview of relevant physical mechanisms of cohesive sediment dynamics in open waters. These include flocculation/aggregation and break-up by turbulent mixing, settling, re-suspension and erosion as well as sedimentation and consolidation.

Table 1. Summary of grain sizes including size range and dominant transport mode in open waters.

Figure 1. Schematic of cohesive sediment dynamics in open waters (figure inspired by Guillou et al., 2011; Maggi, 2005).

Their cohesive properties, and the ease with which they can be resuspended, render fine-grained sediments essential components of many sediment transport processes. Furthermore, their tendency to form flocs and aggregates enables them to act as vehicles for the transport, dispersion and sequestration of contaminants and pollutants. As porous objects, flocs can absorb microscopic materials and substances that would otherwise be in suspension, emulsion or fully dissolved in the fluid. Understanding the dynamics of cohesive sediments is hence crucial for the development of accurate local (Van De Velde, Van Lancker, Hidalgo-Martinez, Berelson, & Meysman, 2018) and global carbon cycle models (Shen, Rosenheim, Törnqvist, & Lang, 2021). At the same time, flocs represent delicate structures whose properties, shape and composition are difficult to investigate experimentally in a non-intrusive way. Consequently, computational research can provide alternative avenues for obtaining insight into the relevant physical processes at play. It is for those reasons that multiphase flows of cohesive sediments have recently been identified as one of the grand challenges in environmental fluid mechanics (Dauxois et al., 2021).

Flocculation has a long tradition as an important tool in process engineering. It can be used in fluidized bed reactors to optimize mixing and chemical processes (e.g. Sundaresan, 2003), or in pneumatic conveying through pipelines (Kuang, Zhou, & Yu, 2020; Molerus, 1996). In wastewater treatment various types of salts or polymers with certain charge densities are frequently used as flocculants in order to remove unwanted substances. They accomplish this by neutralizing the charges of the suspended particles, so that the particles can approach each other, flocculate and settle out (Lee, Robinson, & Chong, 2014). Such process engineering applications are relevant for situations where either air, water or oil (e.g. in pipelines) can act as the carrier fluid. Here the ratio of particle to fluid density, $\rho _p/\rho _f$, plays a crucial role in determining the relative importance of particle inertia. It enters into the definition of the Stokes number $St = T_p / T_0$, which denotes the ratio of the time scales governing particle and fluid motion, $T_p=\rho _p D_p^2/(18 \rho _f \nu _f)$ and $T_0=L_0/U_0$, respectively. Here, $U_0$ and $L_0$ are the characteristic velocity and length scale of the flow, $D_p$ denotes the particle diameter and $\nu _f$ is the kinematic viscosity of the fluid. Very small Stokes numbers, $St\rightarrow 0$, indicate that the particle approximately follows the fluid motion, and thus acts as an ideal tracer. For larger Stokes numbers, particle inertia becomes more influential and the particle trajectory will deviate from the fluid motion. Owing to the lower density ratio of silica in water, subaqueous environmental fluid dynamics tends to deal with lower $St$ flows, while atmospheric research on the motion of particles in air, such as dust and sand storms, or droplets and aerosols in clouds, frequently involves larger $St$ values. On the other hand, cohesive forces are often characterized by the adhesion parameter that reflects the ratio of cohesiveness to kinetic energy of the particles (Li & Marshall, 2007; Marshall & Li, 2014). Since the present review focuses mostly on cohesive sediments suspended in water, the focus will be on applications with Stokes numbers that are smaller and adhesion parameters that are larger than those for particles suspended in gas.

Apart from suspended mineral sediments in marine and riverine open water bodies, an illustrative example of the importance of cohesive sediments in environmental fluid dynamics is the occurrence of marine snow, i.e. organic material falling from the upper ocean layers to the deep ocean (Alldredge & Silver, 1988). This topic attracted considerable attention in the aftermath of the Deepwater Horizon oil spill in 2010, since a substantial fraction of the oil released into the ocean was bound by marine snow and transported from the water surface to the seafloor (Daly, Passow, Chanton, & Hollander, 2016). A further problem of great current interest concerns the sediment plumes caused by autonomous underwater vehicles (AUVs) during deep sea mining operations that gather nodules containing mineral deposits from the seafloor (Ouillon, Muñoz-Royo, Alford, & Peacock, 2022a, 2022b; Peacock, Alford, & Stevens, 2018). It is important to assess the potential for such mining operations to damage sea floor ecosystems by blanketing them with a thin film of sediments. While the sediment suspended by the AUVs is known to propagate as a turbidity current generated by a moving source (Meiburg & Kneller, 2010; Ouillon, Kakoutas, Meiburg, & Peacock, 2021), it is not yet well understood how flocculation influences the dynamics of these flows, and how cohesion affects the permanence of the sediment blanket. Existing turbidity current models are commonly based on continuum model approaches that neglect flocculation and assume a Stokesian settling velocity for the particle concentration field (Biegert, Vowinckel, Ouillon, & Meiburg, 2017b; Necker, Härtel, Kleiser, & Meiburg, 2002, 2005).

For flocculation to occur, particles must come into near contact. This is due to the short range of several nanometres over which cohesive forces act (Berg, 2010; Israelachvili, 1992). There are three main mechanisms that cause particles to approach each other (Partheniades, 2009; Winterwerp, 2002): (i) Brownian motion of nanoscale particles, (ii) differential settling of particles under gravity, and (iii) fluid shear in a viscous flow. The first mechanism requires the particles to be of colloidal size (${<}2\,\mathrm {\mu }$m), so that Brownian motion can compensate for the weight of the particles in an otherwise quiescent fluid. This can be realized in clay suspensions with added dispersants that prevent the primary particles from flocculating. Such suspensions exist in some chemical applications (Lagaly, Schulz, & Zimehl, 2013), but they are not of relevance to most environmental fluid mechanics scenarios. Without dispersant particles will settle, and since not all particles settle with the same speed, differential settling will cause particles to come into contact. Note that polydispersity is not required for this behaviour, as particles interact by hiding and shading mechanisms that alter the effective fluid drag during the settling motion. Such a finite-Reynolds-number effect is nicely illustrated by the drafting-kissing-tumbling scenario described in Fortes, Joseph, and Lundgren (1987). Differential settling is therefore a useful scenario for investigating flocculation due to different settling speeds as a function of the water chemistry that affects the van der Waals forces (Krahl, Vowinckel, Ye, Hsu, & Manning, 2022; Rommelfanger et al., 2022; Sutherland, Barrett, & Gingras, 2015). Nevertheless it represents a somewhat idealized setting because fluid shear is present in nearly all open water bodies. For moderate Stokes number flows, flocculation tends to occur near stagnation points, while the break-up of flocs is typically observed in regions of high shear. In general, even for moderately vigorous flows we can assume that particle contact due to fluid shear will contribute more strongly to flocculation than differential settling (Partheniades, 2009). Under such conditions, we may expect to see a competition between flocculation in areas of preferential concentration, and break-up in regions of high turbulent stresses (Maxey, 1987).

Nevertheless, the analysis of differential settling is important for certain applications, e.g. in waste water treatment. From a fundamental point of view, it provides important information regarding the effective settling velocity of particles in concentrated suspensions. In turn, this affects the runout distance of turbidity currents, and hence, represents an important input parameter for large-scale computational sediment transport models (Biegert et al., 2017b; Cantero, Balachandar, Cantelli, Pirmez, & Parker, 2009; Francisco, Espath, Laizet, & Silvestrini, 2018; Necker et al., 2002, 2005; Olsen, 2021). For dilute suspensions of fine-grained sediments, the Stokes settling velocity and related empirical relationships for grains of different sizes and shapes (Ferguson & Church, 2004; Strom & Keyvani, 2011) serve as adequate approximations. For higher concentrations, the settling velocity decreases as a result of the change of the relative particle buoyancy in the suspension, the upward counterflow of fluid generated by the settling grains and particle–particle interactions including collisions and contact (Dankers & Winterwerp, 2007; Ham & Homsy, 1988; Winterwerp, 2001).

Based on these considerations, empirical correlations for the hindered settling velocity as a function of the local volume fraction have been derived for non-cohesive sediments (Dey, Ali, & Padhi, 2019; Richardson & Zaki, 1954; Shajahan & Breugem, 2020). However, the influence of cohesion on the effective settling rate of flocculated particles is less clear (Liu & Hrenya, 2018; te Slaa, van Maren, He, & Winterwerp, 2015; Winterwerp, 2002). This represents an important knowledge gap, as it has recently been shown that in open water bodies fine-grained sediments exist mostly in a flocculated state (Krahl et al., 2022; Lamb et al., 2020; Rommelfanger et al., 2022). When compared with their non-cohesive counterparts, the settling behaviour of cohesive sediments is quite distinct even for low turbidity levels. Their sedimentation process can be subdivided into three successive stages, viz. flocculation, settling and consolidation, depending on the prevailing processes of the particle dynamics (Adachi, Kawashima, & Ghazali, 2020; Dankers & Winterwerp, 2007; Garcia, 2008).

The dynamics of cohesive particles in turbulent shear flows, and the associated formation and break-up of larger aggregates, plays an important role in the transport of sediments in rivers and oceans, in the erosion of soil by wind and in a wide range of engineering applications including medical devices, as outlined in Zhao et al. (2021). In all of these applications the flocculation process is strongly affected by turbulence, so that the dynamic equilibrium between floc growth and break-up is governed by a complex and delicate balance of hydrodynamic and inter-particle forces.

In order to explore the mechanisms governing this balance, past experimental studies have addressed such aspects as the floc growth rate (Biggs & Lant, 2000; Kuprenas, Tran, & Strom, 2018; Xiao, Yi, Pan, Zhang, & Lee, 2010; Yu, Wang, Ge, Yan, & Yang, 2006), equilibrium size distribution (Bouyer, Line, & Do-Quang, 2004; Chaignon, Lartiges, El Samrani, & Mustin, 2002; Lee, Hyeong, & Cho, 2020; Rahmani, Dabros, & Masliyah, 2004) and transient floc shape (Guérin, Coufort-Saudejaud, Liné, & Frances, 2017; He, Nan, Li, & Li, 2012; Maggi et al., 2007). Following the early work by Levich (1962), population balance models have been developed for describing the floc evolution (Maggi et al., 2007; Shin, Son, & Lee, 2015; Son & Hsu, 2008, 2009; Winterwerp, 1998). Alternative approaches propose statistical collision equations (Ives & Bhole, 1973; Klassen, 2017; Smoluchowski, 1936; Yang et al., 2013).

Over the last decade, large-scale numerical simulations have provided new insight into the interplay of hydrodynamic, inertial and inter-particle forces during the growth, deformation and break-up of aggregates in turbulent flows. While it has not yet been possible to conduct direct numerical simulations over long times that fully resolve a large number of particles much smaller in size than the Kolmogorov scale, various approximations have proved useful, e.g. Marshall and Li (2014) and Dizaji, Marshall, and Grant (2019). Dizaji and Marshall (2016) introduce a stochastic vortex structure method, while Dizaji and Marshall (2017) show that the aggregation process influences the background turbulence only weakly. The simulations of Chen, Li, and Marshall (2019) account for the effects of Stokes drag, lubrication and adhesive contact forces in analysing the early stages of cohesive particle aggregation in homogeneous isotropic turbulence. The follow-up study by Chen and Li (2020) investigates the collision-induced break-up of agglomerates in homogeneous isotropic turbulence. However, as the simulations employ particles of size comparable to the Kolmogorov scale, they do not clarify the role of the Kolmogorov scale in limiting the floc size, which had been observed experimentally (Braithwaite, Bowers, Nimmo Smith, & Graham, 2012; Coufort, Dumas, Bouyer, & Liné, 2008; Fettweis, Francken, Pison, & Van den Eynde, 2006; Kuprenas et al., 2018). In addition, the authors model the cohesive van der Waals force as a ‘sticky force’ that acts only on contact, while previous studies indicated that this attractive force extends over a finite range even before the particles come into contact, so that it can affect the probability that two close-by particles will collide (Israelachvili, 1992; Visser, 1989; Vowinckel, Withers, Luzzatto-Fegiz, & Meiburg, 2019b; Wu, Ortiz, & Jerolmack, 2017).

Certain geophysical applications can give rise to higher sediment volume fractions in excess of 40 %. An important example in this regard are mudslides and debris flows, which pose a major threat to communities, civil infrastructure and human lives. As a result of climate change, hydro-meteorological extremes will likely occur more often in the future, which will result in more frequent successions of drought, fire and heavy rain. As much of the vegetation is destroyed by wildfires, heavy rains can more easily mobilize the soil, which increases the danger of mudslides and debris flows (Kostynick et al., 2022). Such mudslides contain large amounts of suspended fine-grained sediments, and they act as a non-Newtonian carrier fluid that is able to lift and transport large boulders (Ancey, 2007). These flows represent prime examples of how micromechanical interactions at the particle level affect the macroscopic behaviour (Guazzelli & Pouliquen, 2018). Hence, the formulation of quantitative models for the rheological properties of concentrated suspensions of cohesive materials poses an important challenge for the multiphase flow community. It is difficult to extract such rheological information from experiments alone, so that careful simulations can be expected to make a substantial contribution in this regard (Rettinger & Rüde, 2022; Vowinckel, Biegert, Meiburg, Aussillous, & Guazzelli, 2021). At present, the effect of cohesion on the rheology of suspension flows and sediment beds remains largely unexplored.

Further important examples of cohesive sediment flows are found in the context of geotechnical applications, including the scour of consolidated sediments around marine infrastructure such as pipelines, wellheads and offshore wind farms (Whitehouse, 1998), the dredging of river channels and harbour basins (Mehta, 1973; Ravens & Gschwend, 1999), and the removal of sediment from reservoirs, which is frequently accomplished by flushing out sediments (Schleiss, Franca, Juez, & De Cesare, 2016). All of these activities can result in the formation of sediment deposits with unstable slopes where gravitational mass wasting such as landslides and earthfall can occur. As pointed out by Winterwerp (2002), the description of erosion processes is still based primarily on stochastic methods that lack predictive capabilities. Hence, additional research is needed in order to better understand the stability of granular packings and their runout dynamics after failure, especially in the presence of cohesive sediment.

The collapse of granular columns has long served as a canonical test case for studying the mechanisms that govern geophysical granular flows, and for identifying the different rheological flow regimes and scaling laws to which they give rise (Balmforth & Kerswell, 2005; Lajeunesse, Mangeney-Castelnau, & Vilotte, 2004; Lajeunesse, Monnier, & Homsy, 2005; Lube, Huppert, Sparks, & Freundt, 2005, 2007; Lube, Huppert, Sparks, & Hallworth, 2004; Siavoshi & Kudrolli, 2005; Staron & Hinch, 2005). In such flows, the pore pressure feedback mechanisms proposed by Iverson et al. (2000) plays an important role in determining the role of the initial particle volume fraction $\phi$. In submerged cohesionless granular collapses, dense packings result in slow dynamics and short runout distances, while loose packings are associated with more rapid dynamics and longer runout distances (Rondon, Pouliquen, & Aussillous, 2011; Topin, Dubois, Monerie, Perales, & Wachs, 2011; Yang, Jing, Kwok, & Sobral, 2020). The extent to which these mechanisms are modified by cohesive forces is not yet well understood. A few studies have observed a significant influence of cohesion on dry granular collapses, e.g. Rognon, Roux, Wolf, Naaïm, and Chevoir (2006), Mériaux and Triantafillou (2008), Berger, Azéma, Douce, and Radjai (2016) and Mandal, Nicolas, and Pouliquen (2020). In contrast, submerged cohesive collapses have received considerably less attention to date, in spite of their relevance for industrial and environmental applications including sediment transport problems (Baas, Best, & Peakall, 2011; Hampton, 1972; Kuenen, 1951; Marr, Harff, Shanmugam, & Parker, 2001). Very recently, Zhu, He, Zhao, Vowinckel, and Meiburg (2022) took a first step in this direction, by employing particle-resolved simulations, as will be discussed below.

All of the above examples demonstrate how mechanisms acting at the grain scale can dominate the macroscopic behaviour of concentrated suspensions. Since such flows are usually opaque, and since their large-scale dynamics can be quite destructive, some of their aspects are difficult to investigate experimentally in laboratory studies or field measurements. Hence, there exists a strong motivation to develop accurate computational tools for the numerical investigation of cohesive sediment dynamics (e.g. Le Hir, Cayocca, & Waeles, 2011). In this regard, particle-resolved direct numerical simulations (pr-DNS) have emerged as a relatively new tool that holds great promise owing to the rapidly increasing available computational power. This technique allows for the spatial and temporal resolution of the fluid motion, while tracking the motion of the individual grains. Such simulations, which were initially developed for non-cohesive grains (Kempe, Vowinckel, & Fröhlich, 2014; Kidanemariam & Uhlmann, 2017; Uhlmann, 2005; Vowinckel, Kempe, & Fröhlich, 2014), can currently deal with up to $O(10^6)$ grains, so that they can capture a significant range of scales. In recent years they have been extended to cohesive grains as well (Vowinckel, Biegert, Luzzatto-Fegiz, & Meiburg, 2019a; Vowinckel et al., 2019b; Zhu et al., 2022). Tracking the cohesive primary particles in space and time enables us to study, among other things, the microstructure of individual flocs as a function of local shear conditions, or the packing density and stability of sediment deposits, in order to gain a quantitative understanding of the interplay between hydrodynamic and particle contact stresses. Such information is essential for the further development of population balance equations for cohesive sediment flocs (Sherwood et al., 2018; Verney, Lafite, Brun-Cottan, & Le Hir, 2011a), and for modelling debris flows (Garres-Díaz, Bouchut, Fernández-Nieto, Mangeney, & Narbona-Reina, 2020; Iverson, 1997) and the consolidation and stability of cohesive soils (Grasso, Le Hir, & Bassoullet, 2015; Toorman, 1999).

It is evident from the above discussion that cohesive sediment flows pose a number of challenging research problems. A better understanding is needed of the particle–particle interactions, and of the role they play in the flocculation process, for different types of flow fields. Both of these microscale processes strongly affect the macroscale dynamics of cohesive suspension flows, such as their rheology and effective settling rates, along with the resulting sediment deposits. Hence, capturing the dynamical behaviour at the particle level is crucial for the proper macroscopic modelling of such flows. In this regard, grain-resolving simulations can make unique contributions in terms of providing a data base for the formulation of scaling relations from which constitutive equations can be derived that have the ability to predict the large-scale suspension dynamics.

The remainder of this paper is structured as follows. After reviewing a few fundamental properties of cohesive sediments in § 2, in § 3 we provide an overview of state-of-the-art discrete element models (DEM) for cohesive sediments that are coupled to computational fluid dynamics (CFD). Section 4 gives a summary of previous simulation results of CFD-DEM campaigns on flocculation in turbulence, hindered settling and granular collapse. We conclude by giving a brief summary and an outlook in § 5.

2. Physics of cohesive forces

2.1 Short-range interactions of small grains

The interaction of micron-sized grains can be best understood using the classical theory by Derjaguin–Landau–Verwey–Overbeek (DLVO; Derjaguin & Landau, 1941; Verwey & Overbeek, 1948) for colloidal particles, i.e. particle diameters smaller than $2\,\mathrm {\mu }$m. According to the DLVO theory, grains of this size exhibit surface charges that can be expressed as different, additive interaction potentials of either a repulsive or attractive nature. A sketch of the spatial distribution of the relevant potentials as a function of surface distance is given in figure 2.

Figure 2. Sketch of the contributions of different interaction potentials to the net potential of the DLVO theory, where $\varPhi _r$, $\varPhi _a$, $\varPhi _{sr}$ and $\varPhi _{net}$ are the repulsive, attractive, strongly repulsive and net potentials, respectively. While the interaction potentials can have different magnitudes depending on fluid and particle properties, the length scales of $\zeta _n$ plotted logarithmically are based on choices of $\zeta _0=1$ nm (Israelachvili, 1992).

For two approaching surfaces with the same charge, a repulsive potential $\varPhi _r$ is created. This repulsive potential is commonly described by the Gouy–Chapman model as reviewed by Berg (2010). The Gouy–Chapman model is based on the assumption that ions dissolved in the fluid are point charges and surface charges are uniformly distributed across the particle surface. For such conditions, the charged particle surfaces create a so-called double layer, where the first layer (the Stern layer) is a monolayer of counterions binding to the charged surfaces and the second layer (the Gouy layer) prescribes an exponential decay with increasing distance $\zeta _n$ from the surface that is caused by ion diffusion in the free fluid. This yields ${\varPhi _{r} \propto \exp ({-\kappa \zeta _n})}$, where the Debye length $\kappa ^{-1}$ prescribes the length scale for which $\varPhi _r(\zeta _n=\kappa ^{-1})=0.37\varPhi _r( \zeta _n=0)$. The Debye length is, therefore, a function of the fluid properties temperature, salinity and the valency of the dissolved salt ions.

On the other hand, it has been shown by Hamaker (1937) that two approaching surfaces will generate an attractive potential that is commonly known as the van der Waals interaction potential $\varPhi _a$. This potential is caused by the deflection of electrons spinning in orbits around the molecules formed by covalent bonds (Russel, Russel, Saville, & Schowalter, 1991). The deflection can cause a local charge reversal that may persist if the particle surfaces are close enough. The van der Waals potential scales inversely with $\zeta _n$ and is a function of the particle size and shape as well as the Hamaker constant $A_h$, which in itself is dependent on both particle and fluid properties (Lifshitz, 1955). A summary of typical values is given in Vowinckel et al. (2019b) and a common choice for silica materials in water is $A_h\approx 1\times 10^{-20}$ J.

The two potentials $\varPhi _r$ and $\varPhi _a$ are relevant for distances larger than the limiting length scale $\zeta _0$, which represents the microstructure of the surface roughness of a particle and is typically of the order of a nanometre in water. At smaller surface distances, molecular considerations are important, as it becomes increasingly hard to push the last layer of fluid molecules out of the gap. This induces a very steep and strong repulsive potential $\varPhi _{sr}$. It has been shown by Parsons, Walsh, and Craig (2014), however, that the microstructure of the surface roughness usually exceeds this length scale.

For $\zeta _n>\zeta _0$, a competition arises in the net potential $\varPhi _{net}=\varPhi _r+\varPhi _a+\varPsi _{sr}$. Depending on the fluid and particle properties, $\varPhi _{net}$ exhibits three characteristic regions. Very close to the particle, there is a primary minimum that allows for very stable agglomeration. Farther away from the particle surface, there exists a secondary minimum that is less distinct. Particles interacting via the secondary minimum therefore experience weaker attractive forces and we follow the nomenclature by Berg (2010) and refer to such bonding as flocculation. The primary and secondary maximum are separated by an energy barrier induced by the different scalings of $\varPhi _r$ and $\varPhi _a$. It is important to note, however, that this energy barrier does not always exist, but essentially becomes a function of the surface charge and the Debye length. An increase in salinity lowers or even eradicates the energy barrier. The salt concentration for which the energy barrier vanishes is commonly referred to as the critical coagulation concentration (CCC). For distances larger than 100 nm, both potentials decay to zero so that the DLVO theory is relevant for short-range interactions only.

2.2 Assessing the spherical approximation of cohesive grains

As mentioned in the introduction, from a hydromechanical and geotechnical perspective, sediments are typically characterized by their grain size. Following the classification of ISO (2017), these include fractions of clay, silt, sand and gravel (in order of increasing grain size). The first two fractions in this list, silt and clay, are considered to be cohesive. Silt is a granular material consisting of silica that can, to leading order, be described as approximately spherical. The grains are smaller than 63$\,\mathrm {\mu }$m but larger than 2 $\mathrm {\mu }$m. On one hand, we can assume for these types of grains that their weight is sufficiently small for the electrostatic forces described in § 2.1 to remain relevant, albeit of decreasing influence for increasing grain sizes. On the other hand, the grains are significantly larger than colloidal particles, for which Brownian motion would play a role. Consequently, the simulation set-up that will be described below in §§ 3.2 and 3.3 is well suited for silt grains.

The situation is somewhat more complicated for grain sizes smaller than 2 $\mathrm {\mu }$m, i.e. the clay fraction. Generally speaking, clay particles consist of two elementary building blocks. The first are silicate tetrahedra, where four oxygen anions build a structure that fits one silicon cation in its centre. Due to covalent bonding of many tetrahedral units, the structural formula is SiO$_2$. These covalent bonds also allow for unbounded lateral growth, while limiting the thickness of the structural building block. The same consideration holds for the second type of building block, viz. octahedral structures of hydroxy groups that fit magnesium or aluminum cations in their centres. Similar to the tetrahedra, these structures arrange by covalent bonds in the lateral direction, so that the structural formula becomes either Al$_2$(OH)$_6$ (Gibbsit) or Mg$_3$(OH)$_6$ (Brucit). For these reasons, clay particles are platelets with a width of a couple of micrometres, but only a few nanometres thick (Partheniades, 2009). Two different groups of clay minerals can be distinguished based on the composition of their building blocks. A clay platelet that is built by alternating tetrahedral and octahedral plates is considered two layered (TO structure), whereas a recurring elementary cell of tetrahedral-octahedral-tetrahedral plates is considered a three-layered clay mineral (TOT structure). Prominent reprensentatives of these two types are kaolinite and montmorillonite clays, respectively, both of which are very common types of clay minerals in the environment (Partheniades, 2009; Winterwerp, 2002).

Grains of this size are usually considered colloidal, and they experience randomly fluctuating displacements, i.e. Brownian motion (Berg, 2010), that compensate for the particle weight so that the grains stay in suspension indefinitely. Such small particle length scales pose the challenge of simulating representative volume elements of particle suspensions at relevant scales where continuum assumptions remain valid. In fact, the method of choice for these types of clay platelets has been based on molecular dynamics, where every water molecule and all dissolved ions are modelled that surround the crystalline structure of the clay platelets (e.g. Bourg & Sposito, 2011; Chen, Grabowski, & Goel, 2022).

Such simulations are very expensive, and useful primarily for studying interactions of individual platelets. They confirmed observations first described by Weiss (1959): owing to a chemical process known as isomorphic substitution, the centre cations of the tetrahedra and octahedra can be exchanged from Si$^{+4}$ to Al$^{+3}$, or from Al$^{+3}$ to Mg$^{+2}$, respectively. This exchange yields a negative surface charge and holds potential for different types of short-range electrostatic interactions. Consequently, theories from colloidal chemistry can be applied and particles can interact via hydrogen bonding, van der Waals forces and cation bonds. These interactions introduce the propensity of clay particles to flocculate under the right conditions of the ambient fluid. In fact, dry powder of kaolin clay is already arranged in a booklet-like structure of many platelets stacked onto one another (figure 3a). These structures are much larger than $2\,\mathrm {\mu }$m in size, so that they can be considered non-Brownian.

Figure 3. Kaolinite clay visualized by electron microscopy images: (a) dry kaolinite powder (Dohrmann, 2019), (b) kaolinite swollen in an aqueous sodium chloride solution with the salinity of sea water (35 PSU).

According to the DLVO theory from colloidal chemistry outlined in § 2.1, the controlling factor for the electrostatic forces to become either repulsive or attractive is the salt concentration. The critical concentration for this transition to take place is called the CCC. To understand its importance in environmental fluid mechanics, experiments on the differential settling of clay suspensions with different salt concentrations have been carried out (Krahl et al., 2022; Rommelfanger et al., 2022; Sutherland et al., 2015). These experiments were meant to supplement previous efforts in field studies to understand whether there is a correlation between the floc size and the salt concentration (e.g. Droppo & Ongley, 1994; Gibbs, Tshudy, Konwar, & Martin, 1989; Gibbs, 1983; Mikeš & Manning, 2010; Thill et al., 2001). Owing to the multitude of competing physical mechanisms, it had been challenging in these studies to identify the governing parameter that controls the critical coagulation conditions. Rommelfanger et al. (2022) and Krahl et al. (2022) show that for kaolinite and montmorillonite clay, respectively, the CCC is below the salinity of many open water bodies. In these studies, differential settling of clay suspensions was investigated for varying salinities, and a sudden transition to more rapid settling was observed as soon as the salt concentration exceeded the CCC.

For kaolin clay, this observation holds for nearly all situations outside the laboratory, even when tap water was used instead of de-ionized water. Rommelfanger et al. (2022) reported that tap water for UC Santa Barbara supplied by the Goleta Water District (Goleta, 2019) has a salinity higher than the CCC. Expressing salinity in parts per thousand mass fraction of dissolved salt as the practical salinity unit (PSU), Rommelfanger et al. (2022) determined the CCC as 0.04 PSU. An image of a floc observed in salt water of 35 PSU is shown in figure 3(b). The floc was taken from the same batch of kaolin clay shown as a dry powder in figure 3(a). Clearly, the structure of the flocculated clay has changed from plates stacked in a booklet to needles that arrange in an approximately spherical shape. Those aggregates are about $20\,\mathrm {\mu }$m in size but retain their ability to flocculate with other larger aggregates to form even larger flocs. This behaviour corresponds to the idea of flocs as hierarchical structures that was put forward by Krone (1986). In these structures, primary particles or flocculi are structures of zeroth order, flocs that are made of primary particles are of first order, flocs made of first-order aggregates are second order and so on. This hierarchical model has also led to the description of flocs in terms of their fractal dimension (Kranenburg, 1994).

Krahl et al. (2022) conducted an experimental campaign very similar to Rommelfanger et al. (2022), in order to clarify whether the CCC changes for the highly swellable bentonite clay minerals, which consists mostly of montmorillonite and can be considered as a prototype clay with a high cation exchange capacity. As expected, the CCC increased by a factor of $O(10)$ due to the higher surface area and charge induced by the TOT structure of montmorillonite. This salt concentration is at the threshold for the salinity of fresh water and is lower than most of the salinities reported for natural surface waters, especially in estuarine environments of rivers flowing into the ocean. It can therefore be expected that fine-grained suspended sediments are in the flocculated state even in freshwater.

The results of Rommelfanger et al. (2022) and Krahl et al. (2022) have important consequences for the hydrodynamic properties of suspended clays, such as the settling velocities of flocculated particles. As suggested by figure 3(b), the average floc size appears to increase to values in excess of $20\,\mathrm {\mu }$m in an ambient fluid with some salinity. Consequently, such particles in suspension can no longer be described as individual clay platelets (smaller than $2\,\mathrm {\mu }$m). Instead, the primary particles now have the form of several hundred coagulated platelets (Partheniades, 2009; Zhu, Xiong, Liang, & Zhao, 2018). These observations confirm the findings of Lamb et al. (2020) who presented evidence suggesting that all fine-grained sediments in open water are likely to be flocculated. This challenges the existence of the wash-load concept, where all sediments are transported in suspension and never make contact with the sediment bed (Partheniades, 1977). In the context of the research presented in this review, the experimental evidence discussed above justifies the common simulation approach based on a simplified, compact particle shape. In fact, for most purposes, it suffices to treat the primary particles as spheres. It is, however, important to keep in mind that this simplification refers to flocculi that are non-Brownian and cannot be broken up into smaller pieces by the local shear conditions.

3. Computational approaches

For physically realistic, grain-resolved simulations of cohesive sediment dynamics, we need to achieve several goals: first, the fluid motion has to be accurately coupled to the particle dynamics; second, we require proper models and algorithms to account for particle–particle collisions; and finally, we need to capture the attractive cohesive force between particles in close proximity. Here we review strategies for accomplishing these goals.

3.1 Fluid–particle coupling schemes

Computational approaches that solve the Eulerian equations of motion for the fluid, while tracking the particles in a Lagrangian fashion, are commonly referred to as Euler–Lagrange frameworks. These can be further categorized according to the way in which they account for the coupling of fluid and particle motion (Balachandar & Eaton, 2010; Vowinckel, 2021). A one-way coupled system is obtained if the fluid motion remains unaffected by the particles, which are driven solely by their inertia, fluid drag and gravity, without interacting with each other. Two-way coupled systems account for the feedback of the particle drag on the fluid motion. In a three-way coupled system particles are affected by the fluid and interact with each other through collisions and potentially via cohesive forces, while not modifying the underlying fluid velocity field. This approach can be especially useful for investigating the formation of clusters and the agglomeration/flocculation of small suspended particles in relatively dilute flows (Gimbun, Liew, Nagy, & Rielly, 2016; Zhao et al., 2020), but it is less appropriate for studying the hydro-morphological effects of sediment beds on the flow. Finally, a four-way coupled system accounts for the mutual momentum exchange between fluid and particles, as well as for particle–particle interactions via contact and collision. Note that there exists some ambiguity regarding the definition of three-way coupled systems, as they are sometimes defined as suspensions where short-range effects of the particles on the fluid affect the motion of nearby particles (Loth, 2010). Since such situations will usually involve frequent particle collisions (Sommerfeld, 2017), we consider such situations as another form of four-way coupled systems. Particle–particle interactions for three- and four-way coupling are usually implemented via one of several existing discrete element methods. Particles much smaller than the smallest fluid scales resolved by the Eulerian grid are often modelled as point masses, so that they are represented by their instantaneous location, mass and translational velocity (Loth, 2000). On the other hand, particles larger than an Eulerian grid cell are considered to have a finite spatial extent. Since this review focuses on cohesive effects in particle-laden flows, we will primarily discuss three- and four-way coupled systems, where particle–particle interactions by collision and contact are important. Consequently, we have to solve the Navier–Stokes equation

(3.1)\begin{equation} \rho_f \left(\frac{\partial{\boldsymbol{u}}}{\partial{t}} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}\boldsymbol{u})\right) = \boldsymbol{\nabla} \boldsymbol{\tau} + \boldsymbol{f}_{ext}, \end{equation}

along with the continuity equation

(3.2)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0 \end{equation}

for incompressible fluids, where $\boldsymbol {u}$ denotes the fluid velocity, $t$ is time and $\rho _f$ indicates the fluid density. The fluid stress tensor is given by $\boldsymbol {\tau } = -p \boldsymbol {I} + \eta _f (\boldsymbol {\nabla } \boldsymbol {u} + (\boldsymbol {\nabla } \boldsymbol {u})^{\rm T})$, where $p$ represents the fluid pressure with the hydrostatic component subtracted out and $\eta _f$ denotes the dynamic viscosity of the fluid. The right-hand side includes external volume forces $\boldsymbol {f}_{ext}$ that can comprise source terms that can be employed to generate and maintain statistically stationary turbulence in the fluid $\boldsymbol {f}_t$ or the forces that emerge from the fluid–particle interaction $\boldsymbol {f}_{d}$. The latter represents the opposite of the fluid drag force acting on the surface of the particles within a volume element of the computational grid. Hence, $\boldsymbol {f}_{d}$ accounts for the fluid–particle coupling in the momentum balance of the fluid motion. Note that one-way, but also three-way, coupled schemes such as those used by Zhao et al. (2020) and Zhao et al. (2021) as well as Yu, Yu, and Balachandar (2022), do not account for fluid–particle interactions, which yields $f_{d}=0$.

As will be detailed in § 3.2 below, different types of fluid–particle coupling schemes exists. It is therefore important to emphasize that (3.1) and (3.2) represent single-phase flows or flows for which either the particle geometry is fully resolved by the fluid grid, or the particle concentration is very dilute. For sediment volume fractions larger than 0.1 % and particles that are not resolved by the fluid grid, it is common practice to write those equations in their volume-filtered, averaged form (Capecelatro & Desjardins, 2013; Patankar & Joseph, 2001). As an important consequence, the fluid volume fraction, i.e. the porosity, appears in the fluid phase equations (3.1) and (3.2), which also paves the road to larger-scale two-phase flow models for sediment transport (e.g. Chauchat, Cheng, Nagel, Bonamy, & Hsu, 2017).

3.2 Euler–Lagrange frameworks

To track individual sediment grains, we employ the approximation of § 2.2 and compute the particle motion by the Newton–Euler equations for a spherical particle $p$ with uniform density

(3.3)\begin{equation} m_p \frac{\text{d}\boldsymbol{u}_p}{\text{d} t} = \boldsymbol{F}_{h,p} + \boldsymbol{F}_{g,p} + \boldsymbol{F}_{c,p} \end{equation}

for the translational particle velocity $\boldsymbol {u}_p=(u_p,v_p,w_p)^{\rm T}$ and

(3.4)\begin{equation} I_p \frac{ \text{d}\boldsymbol{\omega}_p}{\text{d} t} = \boldsymbol{T}_{h,p} + \boldsymbol{T}_{c,p} \end{equation}

for the angular particle velocity $\boldsymbol {\omega }_p=(\omega _{p,x},\omega _{p,y},\omega _{p,z})^{\rm T}$. Here, $m_p$ is the mass of particle $p$, $\boldsymbol {F}_{h,p}$ the hydrodynamic force acting on particle $p$ and $\boldsymbol {F}_{g,p}$ the force it experiences due to gravity. Here $I_p$ represents the moment of inertia of a spherical particle, $\boldsymbol {T}_{h,p}$ is the torque due to hydrodynamic forces, and $\boldsymbol {F}_{c,p}$ and $\boldsymbol {T}_{c,p}$ denote the force and torque due to particle collisions, respectively (cf. § 3.3 below).

3.2.1 Point-particle approaches

For particles smaller than the Eulerian grid cells, the particle surface is no longer geometrically resolved. Under these conditions one frequently assumes that it suffices to treat the particle as a mass point with a virtual spatial extent that is smaller than the grid cell size (Balachandar, 2009; Balachandar & Eaton, 2010). Since the particle surface is not resolved by the computational mesh, a simplified approach has to be employed to obtain the relevant contributions to the hydrodynamic force acting on the particles

(3.5)\begin{equation} \boldsymbol{F}_{h,p} = \boldsymbol{F}_{d,p} + \boldsymbol{F}_{l,p} + \boldsymbol{F}_{f,p} + \boldsymbol{F}_{a,p} + \boldsymbol{F}_{b,p} + \cdots, \end{equation}

viz. the drag $\boldsymbol {F}_{d,p}$ and lift forces $\boldsymbol {F}_{l,p}$, the pressure gradient force from the undisturbed flow $\boldsymbol {F}_{f,p}$, the added mass force $\boldsymbol {F}_{a,p}$ and the Basset history force $\boldsymbol {F}_{b,p}$, as derived by Maxey and Riley (1983). Depending on the particular flow field under consideration, some of these contributions may be negligible in size, so that (3.5) can be further simplified, or additional terms may have to be accounted for, as indicated by the $(\cdots )$ notation. Since the flow around the particles is not spatially resolved, closure models are needed for the hydrodynamic force contributions. Semi-empirical expressions for the individual terms are typically derived from simplified situations, such as a sphere settling under its own weight in a quiescent, viscous fluid (Basset, 1890; Boussinesq, 1903; Oseen, 1927). For denser systems of volume fractions larger than $10^{-3}$, corrections exist that account for higher particle concentrations (e.g. Bogner, Mohanty, & Rüde, 2015; Tenneti, Garg, & Subramaniam, 2011). In case of two-way or four-way coupling, different schemes have been proposed to compute the particle forcing in the Navier–Stokes equation (3.1) (e.g. Ferrante & Elghobashi, 2003). A robust method that converges under grid refinement has been put forward by Capecelatro and Desjardins (2013),

(3.6)\begin{equation} f_d(\boldsymbol{x}) = \sum_{p = 1}^M \boldsymbol{F}_{h,p} G(\boldsymbol{x}-\boldsymbol{x}_p), \end{equation}

where $G(\boldsymbol {x}-\boldsymbol {x}_p)$ is a kernel, often Gaussian, with units inverse volume and $\boldsymbol {x}_p$ is the centre of particle $p$.

3.2.2 Particle-resolved direct numerical simulations

As reviewed by Biegert, Vowinckel, and Meiburg (2017a), there are numerous frameworks to conduct pr-DNS. A particularly successful and popular one has been the immersed boundary method (IBM) with direct forcing (Biegert et al., 2017a; Kempe & Fröhlich, 2012a; Lucci, Ferrante, & Elghobashi, 2010; Uhlmann, 2005). For this method, the underlying Newton–Euler equations for spherical particles (3.3) and (3.4) are solved, where the hydrodynamic stresses on the particle surface are spatially resolved to yield the hydrodynamic stresses and torques as

(3.7)\begin{equation} \boldsymbol{F}_{h,p} = \oint_{\varGamma_p} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n}\,{\text{d}\varGamma} \end{equation}

and

(3.8)\begin{equation} \boldsymbol{T}_{h,p} = \oint_{\varGamma_p} \boldsymbol{r}\times(\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n})\,{\text{d}\varGamma}, \end{equation}

respectively. In contrast to point-particle methods, a key feature of the IBM is that the hydrodynamic forces and torques $\boldsymbol {F}_{h,p}$ and $\boldsymbol {T}_{h,p}$, respectively, become a direct result of the fluid–particle coupling. The fluid acts on the particles through the hydrodynamic stress tensor $\boldsymbol {\tau }$, where $\boldsymbol {r}$ represents the vector from the particle centre to a point on the surface $\varGamma ^p$ and $\boldsymbol {n}$ is the unit normal vector pointing outwards from that point. The net force and torque acting on the particle centre of mass due to collisions and contacts are again given by $\boldsymbol {F}_{c,p}$ and $\boldsymbol {T}_{c,p}$, respectively.

Vice versa, the fluid at an arbitrary location $\boldsymbol {x}$ experiences drag due to particle motion via $\boldsymbol {f}_d = \boldsymbol {f}_{ibm}$,

(3.9)\begin{equation} \boldsymbol{f}_{ibm}(\boldsymbol{x}) = \rho_f\sum_{p}^{N_{tot}} \sum_{l}^{N_{l,p}}\delta_h(\boldsymbol{X}_l-\boldsymbol{x}) \frac{\boldsymbol{U}_{l,p}^d-\boldsymbol{U}_{l,p}}{\Delta t}V_{l,p}, \end{equation}

where $N_{tot}$ is the number of particles in the domain and $N_{l,p}$ is the number of Lagrangian markers employed to discretize the surface of particle $p$. These Lagrangian markers are placed on the particle in equidistant spacing of the order of the grid cell size to discretize its surface. To couple the mesh of these Lagrangian markers to the fluid grid, a regularized Dirac delta function $\delta _h$ is introduced that spreads the forcing from the Lagrangian marker point on the particle surface to the Eulerian grid (Roma, Peskin, & Berger, 1999; Yao et al., 2022). Furthermore, $\boldsymbol {U}^d_{l,p}$ is the desired velocity at the particle surface, $\boldsymbol {U}_{l,p}$ is the interpolated fluid velocity computed prior to the forcing correction and $V_{l,p}$ is the volume element of the thin shell around particle $p$ that is associated with Lagrangian marker $l$, and that is located inside a given fluid grid cell. Hence, $\boldsymbol {f}_{ibm}$ acts as a correction to the fluid velocity to enforce the no-slip condition on the particle surface.

3.3 Particle–particle interaction by lubrication and contact

We evaluate the collision forces and torques according to Biegert et al. (2017a). The resulting collision model involves normal contact forces, $\boldsymbol {F}_{n}$, frictional contact forces, $\boldsymbol {F}_{t}$, and short-range lubrication forces, $\boldsymbol {F}_{l}$, to provide the total collision force acting on particle $p$ as

(3.10)\begin{equation} \boldsymbol{F}_{c,p} = \sum_{q,\ q \neq p}^{N_p} \left( \boldsymbol{F}_{l,pq} + \boldsymbol{F}_{n,pq} + \boldsymbol{F}_{t,pq}\right), \end{equation}

where the subscript $pq$ indicates interactions with particle $q$. Following Cox and Brenner (1967), Biegert et al. (2017b) propose to model the unresolved component of the lubrication force as

(3.11)\begin{equation} \boldsymbol{F}_{l} = {-}6 {\rm \pi}\rho_f \nu_f R_{eff} \left(\frac{R_{eff}}{\max(\zeta_n,\zeta_{min})} \boldsymbol{g}_n + F^ * _t \boldsymbol{g}_t + F_r^{ * } (R_p \boldsymbol{\omega}_{p}\times \boldsymbol{n} + R_q \boldsymbol{\omega}_{q}\times \boldsymbol{n})\right) \end{equation}

and torque

(3.12)\begin{equation} \boldsymbol{T}_{l} = 8 {\rm \pi}\rho_f \nu_f R_{eff}^2 [\boldsymbol{g}_t T_t^{ * } + T_r^{ * } (R_p \boldsymbol{\omega_{p}} \times \boldsymbol{n} + R_q \boldsymbol{\omega}_{q}\times \boldsymbol{n} ) ]\times \boldsymbol{n} \end{equation}

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$, $\boldsymbol {g}_n$ and $\boldsymbol {g}_t$ are the relative velocities in the normal and tangential directions, respectively, between the two-particle surfaces at the point of contact, $\zeta _n$ is the surface distance between the two particles and $\zeta _{min} = 3\times 10^{-3}R_{eff}$ is a limiter as calibrated by Biegert et al. (2017a) preventing the lubrication force from reaching its singularity at $\zeta _n \to 0$. The terms $F_t^*$, $F_r^*$, $T_t^*$ and $T_r^*$ were obtained via asymptotic expansions by Goldman, Cox, and Brenner (1967),

(3.13)\begin{align} F_t^{ * } &\sim \frac{8}{15} \ln\left(\frac{\max(\zeta_n,\zeta_{min})}{R_{{eff}}}\right) - 0.9588, \end{align}

(3.14)\begin{align} F_r^{ * } &\sim{-}\frac{2}{15} \ln\left(\frac{\max(\zeta_n,\zeta_{min})}{R_{eff}}\right) - 0.2526, \end{align}

(3.15)\begin{align} T_t^{ * } &\sim{-}\frac{1}{10} \ln\left(\frac{\max(\zeta_n,\zeta_{min})}{R_{eff}}\right) - 0.1895, \end{align}

(3.16)\begin{align} T_r^{ * } &\sim \frac{2}{ 5} \ln\left(\frac{\max(\zeta_n,\zeta_{min})}{R_{eff}}\right) - 0.3817. \end{align}

It has been shown by Rettinger and Rüde (2022) that the proper implementation of the lubrication force model is key to model the dynamics of low-Reynolds-number particle-laden flows.

The repulsive normal component is represented by a nonlinear spring-dashpot model for the normal direction

(3.17)\begin{equation} \boldsymbol{F}_{n,pq} = {-}k_n |\zeta_n|^{3/2} \boldsymbol{n} - d_n \boldsymbol{g}_{n,cp}, \end{equation}

where $d_n$ and $k_n$ represent stiffness and damping coefficients that are adaptively calibrated for every collision/contact to reproduce a prescribed restitution coefficient of $e_{p} = -u_{out}/u_{in}$ (Kempe & Fröhlich, 2012b). Here, $\boldsymbol {g}_{n,cp}$ is the normal component of the relative particle velocities and $u_{out}$ and $u_{in}$ indicate the respective normal components of the relative particle speed immediately after and right before the particle impact, i.e. $\zeta _n=0$. The forces in the tangential direction are modelled by a linear spring-dashpot model capped by the Coulomb friction law as

(3.18)\begin{equation} \boldsymbol{F}_{t,pq} = \min ({-}k_t \boldsymbol{\zeta}_t - d_t \boldsymbol{g}_{t,cp}, \|\mu_f \boldsymbol{F}_n\| \boldsymbol{t} ), \end{equation}

where $\mu _f$ represents the coefficient of friction between the two surfaces and $\boldsymbol {\zeta }_t$ is the tangential displacement integrated over the time interval for which the two particles are in contact (Thornton, Cummins, & Cleary, 2013). This contact-modelling framework has been calibrated and validated in detail by Biegert et al. (2017a) and Vowinckel et al. (2021) against seminal experimental benchmark data involving glass spheres (Aussillous, Chauchat, Pailha, Médale, & Guazzelli, 2013; Foerster, Louge, Chang, & Allia, 1994; Gondret, Lance, & Petit, 2002; Ten Cate, Derksen, Portela, & Van Den Akker, 2004).

3.4 Cohesive force model for Euler–Lagrange simulations

The modelling framework described in § 3.2 represents the classical approach of combining a CFD technique with a DEM in order to model the dynamics of spherical particles in viscous flows. In the following we review the three most popular approaches for implementing cohesive forces in this framework. One potential implementation is based on the Tabor parameter as the ratio of the elastic and adhesive interaction range (Johnson & Greenwood, 1997)

(3.19)\begin{equation} \lambda_T = \left(\frac{4 R_r\gamma^2}{E_r^2\zeta_0^3}\right)^{1/3}, \end{equation}

where $R_r=[1/R_p + 1/R_q]^{-1}$ denotes the reduced radius, the subscripts $p$ and $q$ denote the particle indices and $r$ stands for reduced. Furthermore, $E_r=[(1-\nu _p^2)/E_p+(1-\nu _q^2)/E_q]^{-1}$ is the elastic modulus of the particles, $\nu _{p/q}$ is the Poisson ratio of particle $p$ and $q$, respectively, $\gamma =A_h/(24{\rm \pi} \zeta _0^2)$ is the potential energy associated with van der Waals force, $A_h$ is the Hamaker constant and $\zeta _0$ is the minimal separation distance that is typically of the order of a couple of nanometres (Israelachvili, 1992; Vowinckel et al., 2019b).

Based on this non-dimensional parameter, two regimes can be differentiated (Marshall & Li, 2014): (i) $\lambda _T\gg 1$, i.e. large, ‘soft’ particles that have a high potential adhesive energy, and (ii) $\lambda _T\ll 1$, i.e. small, ‘stiff’ particles with a low potential adhesive energy. For large $\lambda _T$, the particle surface deforms so strongly that it considerably enlarges its contact surface, which yields a nonlinear coupling of collision and cohesive forces (Yao & Capecelatro, 2021). On the other hand, small $\lambda _T$ indicate minor surface deformations and the cohesive forces remain unaffected. Consequently, the collision and cohesive forces can be treated independently in an additive manner. Sketches of these two principles to model the interaction of cohesive particles with a wall are given in figure 4.

Figure 4. Sketch of the underlying principles for a cohesive particle interacting with a wall: (a) the model by Johnson et al. (1971) where cohesion and contact are nonlinearly coupled through the radius of the contact area $a$ and (b) the additive model of Derjaguin et al. (1975).

In the following we review a selection of different techniques that have been used to model cohesive particles in CFD-DEM applications. The techniques are categorized based on the conditions introduced by $\lambda _T$. These categories can be subsumed under the model approaches of Johnson, Kendall and Roberts (JKR, Johnson, Kendall, & Roberts, 1971) and Derjaguin, Muller and Toporov (DMT, Derjaguin, Muller, & Toporov, 1975), respectively.

3.4.1 The JKR-type models

As mentioned earlier, the JKR theory of Johnson et al. (1971) is appropriate for soft materials with moderate to high surface energy. In fact, the JKR model was first developed for rubber materials (Johnson & Greenwood, 1997). Nevertheless, it has been found useful for DEM simulations to investigate, e.g. the mixing, and (de-)agglomeration of fine powders as reviewed by Chen and Elliott (2020). One model compatible with DEM was put forward by Chokshi, Tielens, and Hollenbach (1993) and implemented in the framework of CFD-DEM simulations by Marshall (2009). Subsequently, this approach was successfully applied to simulate the break-up of particle agglomerates in linear shear flow (Dizaji & Marshall, 2017) and isotropic turbulence (Chen & Elliott, 2020).

According to Marshall (2009), the normal contact forces between particles $p$ and $q$ are expressed as

(3.20)\begin{equation} F_{n,pq} = F_{n,pq}^{E} + F_{n,pq}^{D}, \end{equation}

where $F_{n,pq}^{D}$ is a damping force to set the restitution coefficient to zero (Dizaji & Marshall, 2017) and $F_{n,pq}^{E}$ represents the forces from van der Waals attraction and elastic deformation. These forces are nonlinearly coupled to the radius of the contact region as a function of time, $a(t)$, via

(3.21)\begin{equation} F_{n,pq}^{E} = {-}4F_c\biggl[ \left(\frac{a}{a_0}\right)^3-\left(\frac{a}{a_0}\right)^{3/2}\biggr], \end{equation}

where $F_c$ it the critical normal pull-off force and $a_0$ is the equilibrium contact region radius. The area of the contact region is computed by solving

(3.22)\begin{equation} \delta_n = 6 \delta_c\biggl[ 2\left(\frac{a}{a_0}\right)^2-\frac{4}{3}\left(\frac{a}{a_0}\right)^{1/2}\biggr], \end{equation}

where $\delta _n=-\zeta _n$ is the particle overlap of the two interacting spheres. The pull-off force is the force needed to pull particles apart if they are bonded by cohesive forces only, i.e. $\delta _n = \delta _c$. Owing to the nonlinear coupling of contact forces and adhesive forces through $a(t)$, this yields

(3.23)\begin{equation} F_c = 3{\rm \pi} \gamma R_r \end{equation}

and $\delta _c$ is the critical distance for which necking of the two particles remains possible in the contact region. Note that (3.21) replaces the terms of (3.10) due to the nonlinear coupling of cohesive force and surface deformation during the collision of the soft spheres. Similarly, the torque induced by cohesion upon contact is rewritten as

(3.24)\begin{equation} T_{c,p} = {-}k_R \zeta_t, \end{equation}

where $k_R=4F_c(a/a_0)^{3/2}$ is the tangential stiffness for cohesive materials and $\zeta _t$ is the same tangential displacement that enters (3.18). With this type of contact modelling, particles remain in contact until the tensile normal forces exceed the critical pull-off force $F_c$ and the particles move apart beyond $\delta _c$. Even though the JKR model has been derived for soft materials with $\lambda _T\gg 1$, Johnson and Greenwood (1997) claim that similarly good results can be obtained for materials with a smaller Tabor parameter.

A somewhat similar approach to the JKR model is the so-called hit-and-stick assumption. Under this assumption, particles will always stick together and form larger aggregates as they come into contact. It was pointed out by Chen et al. (2019) that such an assumption is only valid if there are exclusively low energy collisions in the system. The hit-and-stick assumption has been a popular tool for dust aggregation (Okuzumi, Tanaka, & Sakagami, 2009) and even for planet formation in protoplanetary disks (Güttler, Blum, Zsom, Ormel, & Dullemond, 2010), but these types of studies are usually done for DEM simulations in dry conditions only.

To model the potential energy suddenly released to the particles upon detachment in pr-DNS, Derksen (2014) proposed a square-well potential for agglomerated particles. This potential applies to two particles bonded by cohesive forces until the tensile force pulls them apart and the potential energy is converted into kinetic energy. Delenne, El Youssoufi, Cherblanc, and Bénet (2004) proposed a contact force model that introduces cohesive bonds that establish a rigid bridging of the particles. These bonds fail as soon as the yield load is reached and the usual contact and friction laws are used.

3.4.2 The DMT-type models

Even though the JKR model has been applied with success to the deagglomeration of cohesive particles, many studies deal with materials that are quite stiff, such as silica grains and clay flocculi described in § 2.2. In addition, the nonlinear coupling of forces due to collision and cohesion poses challenges to the numerical stability of classical CFD-DEM solvers. For these reasons, many studies have adopted the additive DMT model that is more in line with the underlying idea of the Newton–Euler equations (3.3) and (3.4). For example, van Wachem, Thalberg, Remmelgas, and Niklasson-Björn (2017) model collisions by a spring-dashpot system and introduce an additional cohesive contact force $F_{vdW}$ that acts normal to the surface and is proportional to the particle overlap $\delta$ at a finite range $\delta _0$ representing the equilibrium distance for cohesive contact. For larger overlaps, the cohesive forces are equivalent to the critical pull-off force $F_c/2$.

For the DMT model, $F_c$ is computed via the Hertzian contact theory, which yields $F_c=4{\rm \pi} R_r \gamma$, where again $\gamma =A_h/(24{\rm \pi} \delta _0^2)$ is the surface energy due to van der Waals forces. Compared with the JKR model (3.23), this slightly larger value is the direct consequence of the additive and nonlinear coupling of the DMT- and JKR-framework approaches, respectively. A detailed derivation of $F_c$ for both frameworks can be found in Marshall and Li (2014). Owing to the additive nature of the DMT model, this cohesive force $F_{vdW}$ can subsequently be added to the resulting force from particle interaction and contact (3.10).

Similarly, Yao and Capecelatro (2021) followed Gu, Ozel, and Sundaresan (2016) and defined cohesive forces as $F_{vdW}\propto \zeta _2^{-2}$ that act in a finite gap in between particles and becomes constant for particle overlap. These forces are added to the classic spring-dashpot model for contact, so that they can directly enter (3.10). Full details are provided in Gu et al. (2016) and Yao and Capecelatro (2021). Yao and Capecelatro (2021) investigated the deagglomeration of cohesive aggregates by isotropic turbulence and determined $0.19\leq \lambda _T\leq 0.98$. They found that their results are not sensitive to the choice of using either a DMT or a JKR model for the cohesive forces.

A recent study by Yu et al. (2022) followed the example of Yao and Capecelatro (2021) but, in addition, these authors claim that adhesive forces according to the JKR model need to be included. They follow the reasoning by Parteli et al. (2014) and propose to add an additional force during contact that is proportional to the radius of the contact region $a$ (Kendall, 1971). Therefore, this approach represents a mixture of DMT and JKR modelling and it contradicts the rationale put forward by Marshall and Li (2014) and Yao and Capecelatro (2021), who argued that the JKR model implies that collision and cohesion do not represent additive effects. Nevertheless, Yu et al. (2022) have successfully used this model to generate data for particles flocculating in isotropic turbulence to quantify rates of aggregation and disaggregation, respectively.

3.4.3 Cohesive force model compatible with point-particle and particle-resolved simulations

For geometrically resolved CFD-DEM applications, i.e. pr-DNS, the modelling approach poses some additional difficulties that go beyond the DEM approaches summarized in §§ 3.4.1 and 3.4.2. Owing to the high resolution of pr-DNS, the computational model needs to properly reflect the following three phases of particle–particle interaction: (a) particle approach/flocculation, (b) capture/steady state contact and (c) separation/deagglomeration in the presence of external forces. Such a model provides detailed information on the work performed by the inter-particle forces upon approach and separation. To this end, Vowinckel et al. (2019b) propose a model that follows the DMT theory, where collision and cohesion are treated as additive effects. Hence, the computational model recovers the original DEM scheme by Biegert et al. (2017a) for cohesionless grains in the framework of IBM simulations of particle-laden flows.

To guarantee the compatibility with the soft-sphere DEM of Biegert et al. (2017a), Vowinckel et al. (2019b) developed a model that does not alter the original algorithm for contact and collision and is consistent with the DLVO theory as sketched in figure 2. In particular, Vowinckel et al. (2019b) aim to provide an algebraic expression that mimics the secondary maximum of the DLVO curve, where flocculation is possible. According to the DLVO theory, the secondary maximum yields an attractive inter-particle force within the interval $2\,\text {nm} \leq \zeta _n \leq 10\,\text {nm}$ that has a local maximum at $\zeta _n \approx 4\,\text {nm}$. Length scales of $\zeta _n < 2\,\text {nm}$ are not considered as they are part of the surface roughness. To this end, the model is designed such that cohesive forces decay to zero for $\zeta _n = 0$, while the repulsive forces act during contact through (3.17). Such a model is well represented by a parabolic spring force

(3.25) \begin{equation} \boldsymbol{F}_{coh} = \begin{cases} -k_{coh}(\zeta_n^2 - \zeta_n \lambda) \boldsymbol{n} & 0, < \zeta_n \leq \lambda, \\ 0 & \text{otherwise}, \end{cases} \end{equation}

where $k_{coh}$ denotes the stiffness constant. As desired, (3.25) has 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. As an important parameter, $\lambda$ represents the range over which the cohesive force is smeared. This measure becomes necessary because length scales of a few nanometres are typically not resolved in CFD.

As was shown by Vowinckel et al. (2019b), the simulation results are insensitive to the exact value of $\lambda$. For pr-DNS, $\lambda$ has to be large enough so that cohesive forces can be numerically resolved, but much smaller than the particle size, e.g. $D_{50} / \lambda = 20$, where $D_{50}$ is the median grain size of a polydisperse grain size distribution. Owing to the substepping routine proposed by Biegert et al. (2017a), $\lambda$ can be as small as the grid cell size $h$, where a typical resolution of pr-DNS is 20 grid cells per diameter. Such a choice still guarantees a proper resolution of cohesive effects in space and time. This modelling approach is also consistent with the experimental observations of Delenne et al. (2004), who used the same data to derive their cohesion model outlined in § 3.4.1. In these experiments, rods of $D = 8$ mm in diameter were coated with epoxy resin to glue them together. These rods where then put under tension to determine the cohesive forces. It was found in this study that the cohesive force increases with gap size $\zeta _n$ to a maximum at $D / \zeta _n \approx 80$. For larger gap sizes, the force decreases and eventually, the rods detach.

Vowinckel et al. (2019b) determine the stiffness $k_{coh}$ of the model by preserving the energy contained in the van der Waals forces (Israelachvili, 1992)

(3.26)\begin{equation} F_{vdW} = \frac{A_h R_{eff}}{6 \zeta_0^2}, \end{equation}

which yields

(3.27)\begin{equation} k_{coh} = \frac{A_h R_{eff}}{\zeta_0 \lambda^3 } \end{equation}

to obtain the dimensional cohesive force model

(3.28) \begin{equation} \boldsymbol{F}_{coh} = \begin{cases} -\dfrac{A_h R_{eff}}{\zeta_0 \lambda^3 }(\zeta_n^2 - \zeta_n \lambda) \boldsymbol{n}, & 0 < \zeta_n \leq \lambda, \\ 0 & \text{otherwise}. \end{cases} \end{equation}

This form provides an expression for cohesive forces in pr-DNS via an additional force term in the collision model (3.10)

(3.29)\begin{equation} \boldsymbol{F}_{c,p} = \sum_{q,\ q \neq p}^{N_p} \left( \boldsymbol{F}_{l,pq} + \boldsymbol{F}_{n,pq} + \boldsymbol{F}_{t,pq} + \boldsymbol{F}_{{coh},pq}\right) + \boldsymbol{F}_{l,pw} + \boldsymbol{F}_{n,pw} + \boldsymbol{F}_{t,pw} + \boldsymbol{F}_{{coh},pw}, \end{equation}

where $pq$ and $pw$ indicate collisions of particle $p$ with another particle of index $q$ or the wall, respectively. Such a framework can be used in conjunction with the DEM of Biegert et al. (2017a) without any further modification and the hydrodynamic forces and the particle weight, $\boldsymbol {F}_{h,p}$ and $\boldsymbol {F}_{g,p}$, respectively, remain a result of the original IBM. Such a model resolves the particle bonding process in space and time and, at the same time, it uses the DMT rationale to retain the distinction between the individual inter-particle force components via (3.29).

For pr-DNS, it is convenient to write the dimensional form (3.28) in terms of non-dimensional quantities. This especially applies to the proper parameterization of the Hamaker constant $A_h$ with respect to particle inertia. To this end, suitable non-dimensional numbers are needed to define cohesive effects for arbitrary systems. As usual, one can render the Navier–Stokes equation (3.1) dimensionless by choosing a proper reference length scale, velocity and density ($L_0$, $U_0$ and $\rho _0$, respectively), which yields the Reynolds number $Re=U_0 L_0 / \nu _f$ as a dimensionless parameter (Biegert et al., 2017b) where again $\nu _f$ represents the kinematic viscosity. Applying the same logic to the particle equation of motion (3.3), non-dimensional values emerge for the lubrication force and the cohesive force scales. By writing the algebraic expression for cohesive forces (3.28) in dimensionless form, Vowinckel et al. (2019b) obtained the cohesive number

(3.30)\begin{equation} Co = \frac{\text{max}(\| \boldsymbol{F}_{{coh},50} \|)}{U_0^2L_0^2\rho_0} \end{equation}

as the ratio of the cohesive force maximum to the characteristic inertial force scale of the problem under consideration. Similar characteristic numbers have been introduced for CFD-DEM applications in the form of an adhesion parameter $Ad=\gamma / (\langle u_p'\rangle ^2 R_p^2 \rho _p)$ as the ratio of surface energy due to the mean kinetic energy based on the granular temperature of the particles (e.g. Dizaji & Marshall, 2017; Li & Marshall, 2007; Yao & Capecelatro, 2021) or the bond number $Bo$ as the ratio of cohesive forces to the particle weight (e.g. Sun, Xiao, & Sun, 2018). In this context, $Bo$ should not be confused with the Bond number, which represents the ratio of gravitational forces to surface tension forces and is unrelated to cohesive bonding. This dimensionless number is named after Wilfrid Noel Bond and is commonly used to characterize the shape of bubbles in a viscous flow (Clift, Grace, & Weber, 2005). For the applications presented in the present review, the choices of the relevant velocity and length scales, $U_0$ and $L_0$ that enter (3.30), very much depend on the problem under considerations as will be detailed in § 3.5 below

To determine the numerator of (3.30), Vowinckel et al. (2019b) compute the maximum cohesive force at $\zeta _n = \frac {1}{2} \lambda$ to obtain, for $R_{eff} = D_{50} / 2$,

(3.31)\begin{equation} \text{max}(\| \boldsymbol{F}_{{coh},50} \|) = {-}\frac{A_h}{\zeta_0} \frac{D_{50}}{2 \lambda^3} \left(\frac{\lambda^2}{4} - \frac{\lambda^2}{2}\right) = \frac{A_h}{\zeta_0} \frac{D_{50}}{8 \lambda}, \end{equation}

which yields the dimensionless cohesive force as

(3.32)\begin{equation} \tilde{\boldsymbol{F}}_{coh} = \begin{cases} - Co \dfrac{8 R_{eff}}{\lambda}^2(\tilde{\zeta}_n^2 - \zeta_n \lambda) \boldsymbol{n}, & 0 < \zeta_n \leq \lambda, \\ 0 & \text{otherwise}. \end{cases} \end{equation}

The stiffness of the cohesive force model is therefore determined as $k_{coh} = 8 {Co} R_{eff} / \lambda ^2$, where the tilde indicates dimensional values for the effective radius and the cohesive force range. As desired, the cohesive forces for a given physical system scale linearly with the cohesive number and the effective radius of the two interacting particles, a scaling that is consistent with the considerations of Visser (1989), Lick, Jin, and Gailani (2004) and Righetti and Lucarelli (2007). According to § 2.2, such a model is meant to represent rough macroscopic (non-Brownian) particles. This model is based on principles of the DLVO theory, and at the same time it is compliant with the DMT theory. Given the advantages outlined above, such a model is especially useful for flocculation processes in pr-DNS (Vowinckel et al., 2019a, 2019b; Zhu et al., 2022), but it can also be applied to point-particle approaches undergoing continuous agglomeration and break-up in a quasisteady state (Zhao et al., 2021, 2020).

3.5 Scaling considerations and governing dimensionless parameters

As mentioned in the previous section and the introduction, a proper parameterization of the problem under consideration requires an appropriate choice of length and velocity scales to obtain a suitable Reynolds number. Depending on the scenario under consideration, these scales can be related to the fluid motion or the particle properties, as will be discussed in more detail below. Hence, the specific choice of scales determines if $Re$ represents the ratio of inertial and viscous forces acting on a fluid element, or on a solid particle suspended in the flow. In particular, the inertial and viscous forces acting on a solid particle govern the ratio of the particle response time $T_p=\rho _s D_p^2/(18 \nu _f)$ to the fluid time scale $T_0=L_0/U_0$ commonly referred to as the Stokes number

(3.33)\begin{equation} St = \frac{T_p}{T_0} = \frac{\rho_s}{18} \frac{D_p^2}{L_0^2} Re, \end{equation}

where $\rho _s$ is the ratio of particle to fluid density $\rho _p/\rho _f$.

For simulating a given physical situation, the relevant characteristic scales should be based on the specific mechanisms that bring particles into close proximity of each other, so that their flocculation can be triggered via the short-range cohesive forces that are active over a range of the order of nanometres (cf. § 2.1). Partheniades (2009) and Winterwerp (2002) list three mechanisms that are able to bring particles into contact: (1) Brownian motion, (2) differential settling and (3) fluid shear. We remark that the term ’fluid shear’ in this context refers to the general existence of fluid velocity gradients, and it does not necessarily imply the presence of vorticity, as it is well known that the irrotational flow near a stagnation point can result in the preferential concentration of particles (e.g. Maxey, 1987; Raju & Meiburg, 1997).

It was discussed in § 2.2 that cohesive sediments in natural open water bodies occur as flocculi, which we consider as the primary particles in the present context. These primary particles have typical sizes of the order of several tens of micrometres, and they cannot be broken down into smaller particles by local hydrodynamic stresses, so that the influence of Brownian motion can be neglected. Hence, our considerations for non-dimensionalization will focus only on differential settling and fluid shear as the mechanisms that bring particles into contact. Differential settling occurs in situations where particles settle at different speeds due to different weights or shapes, or as a result of local fluid velocity gradients caused, for example, by wake effects in a suspension. Vowinckel et al. (2019b) demonstrated that differential settling remains relevant even for monodisperse particles of equal weight.

In the following we illustrate how cohesive forces as introduced in § 3.4 affect the non-dimensionalization for scenarios governed by differential settling or fluid shear. To this end, we follow the reasoning of Vowinckel et al. (2019b) and consider the dimensionless equation of motion for the particles (3.3) with non-dimensional variables $\boldsymbol {u}_p=U_0\tilde {\boldsymbol {u}}_p$, $D_p=L_0\tilde {D}_p$, $\rho _f=\rho _0\tilde {\rho }_f$ and $m_p = \rho _f V_{50} \tilde {m}_p=m_{50}\tilde {m}_p$. Here, $V_{50}$ denotes the volume of the median grain size of a polydisperse particle distribution, and the tilde symbol indicates a dimensionless variable. Using these non-dimensional values yields a characteristic force $m_{50}U_0^2/L_0$ that varies depending on the choice for the characteristic scaling under consideration. Normalizing (3.3) by this scaling yields the non-dimensional cohesive number $Co$ as a pre-factor of the cohesive force term (3.32) that governs its magnitude. For the case of hindered settling of polydisperse particles, we choose the median grain size as the characteristic length scale, i.e. $L_0=D_{50}$, and we define a buoyancy velocity $U_0=u_b=\sqrt {g'D_{50}}$, where $g'=(\rho _p-\rho _f)g/\rho _f$ is the reduced gravity. For these choices, the Reynolds number is identical to the Galileo number $Ga$,

(3.34)\begin{equation} Ga = \frac{D_{50}u_b}{\nu_f}, \end{equation}

and the cohesive number becomes

(3.35)\begin{equation} Co = \frac{\mathrm{max}(\|\boldsymbol F_{coh,50}\|)}{m_{50}g'} = \frac{A_h D_{50}}{8\lambda \zeta_0}\frac{1}{m_{50}g'}. \end{equation}

This definition has a straightforward physical interpretation. Consider a particle stuck under a horizontal solid surface. A critical value of $Co=1$ represents the condition at which this particle begins to detach because the weight exceeds the cohesive force holding the particle against gravity (Vowinckel et al., 2019b).

The cohesive forces need to be scaled differently if it is fluid shear that brings particles into contact, rather than differential settling. For example, in homogeneous isotropic turbulence the intensity of the fluid shear determines the Kolmogorov length scale $\eta$ as the smallest scale of a turbulent eddy for which inertial effects remain significant. Hence, we choose as characteristic scales $L_0=\eta$ and $U_0=u_\eta =\nu _f/\eta$, which yields $Re_\eta =u_\eta \eta /\nu _f=1$ (Pope, 2001). Using this parameterization, the local shear rate becomes $G=\nu _f/\eta ^2$ and the cohesive number takes the form

(3.36)\begin{equation} Co = \frac{{\rm max} (\|{\boldsymbol F_{coh,ij}}\|)}{u_\eta^2 \eta^2 \rho_f} = \frac{A_h D_{50}}{8 \lambda \zeta_{0}}\frac{1}{u_\eta^2 \eta^2 \rho_f}. \end{equation}

Choosing a certain dimensional value for $D_{50}$ and using characteristic values for the Hamaker constant $A_h$, we can link this non-dimensional description to dimensional physical properties. Vowinckel et al. (2019b) provide a review of characteristic values for the Hamaker constant, which is typically of the order of $A_h=1\times 10^{-20}$ J for silica materials in water. Note, however, that a wide range of values for $A_h$ have been reported in the literature. Furthermore, typical values for the minimal separation distance and the cohesive range are determined as $\zeta _0=5\, \times\, 10^{-9}$ m and $\lambda =D_{50}/20$ in this reference. Since $m_{50}\propto D_{50}^3$, (3.30) states that cohesive forces decrease as $D_{50}^{-2}$. For example, for a grain size distribution with a median grain size of $D_{50}=5\,\mathrm {\mu }$m, $Co$ becomes $9.36$, but reaches a value of 0.0047 for $D_{50}=63\,\mathrm {\mu }$m. As desired, $Co\ll 1$ for the latter choice, which marks the desired upper threshold value for which cohesive effects remain important. It is therefore the linear scaling by $R_{eff}$ in (3.32) that determines the cohesiveness of a sediment grain. Note that this argument can be made for various choices of median flocculi size. For example, Vowinckel et al. (2019b) and Vowinckel et al. (2019a) investigated a differential settling scenario with $L_0=D_{50}$ and $U_0=u_b$ and they were able to show that grains with $Co=5$, $Ga=1.35$ and $St=0.09$ roughly correspond to the properties of cohesive silt flocculi with a median grain size of $D_{50}=20\,\mathrm {\mu }$m.

Using the same characteristic quantities $L_0=D_{50}=D_p$ and $U_0=u_b$, Zhu et al. (2022) used this rationale to simulate the granular collapse of monodisperse particles with diameter $D_p$ to obtain $0\leq Co \leq 290$, $Ga=200$ and $St=12.79$. The large values of the non-dimensional numbers correspond to sediment grains on the millimetre scale. Although it was discussed in § 2.2 that grains of this size are generally considered non-cohesive, these choices are useful in the context of experiments employing adhesively coated grains (Brunier-Coulin, Cuellar, & Philippe, 2020; Gans, Pouliquen, & Nicolas, 2020; Jarray, Shi, Scheper, Habibi, & Luding, 2019). Note that even though the concept of introducing the cohesive number is motivated by electrostatic van der Waals forces that scale with the Hamaker constant $A_h$, the concept can be further generalized to additional cohesive and adhesive effects due to biofilms (Krahl et al., 2022; Lamb et al., 2020) or the addition of bentonite to sandy beds (Lick et al., 2004).

4. Euler–Lagrange simulations of cohesive particles

Euler–Lagrange approaches have been employed by a number of authors in recent years for the purpose of investigating the effects of cohesive forces on the collective dynamics of dispersed particulate flows. In the following we will focus primarily on several recent studies that discuss the flocculation and break-up of cohesive sediment in turbulence, as well as particle-resolved simulations of the collective settling dynamics of such sediment, and of cohesive granular collapses.

4.1 Unresolved cohesive particles in turbulence

In the following we will review computational investigations that studied flocculation in turbulent shear flows by means of three-way coupled simulations. The fluid–particle interactions for such an approach were outlined in §§ 3.1 and 3.2.1 and the modelling of the particle–particle interactions were described in §§ 3.3 and 3.4.3. We will begin with a conceptually simple cellular model flow that clarifies some fundamental aspects of the competition between inertial, drag and cohesive forces during the flocculation process (Zhao et al., 2020). Subsequently, we will discuss simulations of flocculation in homogeneous isotropic turbulence (Zhao et al., 2021).

4.1.1 Flocculation in cellular model flows

An instructive setting for investigating the three-way coupled dynamics of flocculation and break-up is given by the Taylor-Green cellular flow shown in figure 5(a), with fluid velocity field ${\boldsymbol u}_f = (u, v)^{\rm T}$,

(4.1a,b)\begin{equation} u = \frac {U_0}{\rm \pi} \sin\left(\frac{{\rm \pi} x}{L}\right) \cos\left( \frac{{\rm \pi} y}{L}\right),\quad v_f = {-}\frac {U_0}{\rm \pi} \cos\left( \frac{{\rm \pi} x}{L}\right) \sin\left( \frac{{\rm \pi} y}{L}\right). \end{equation}

Here $L$ and $U_0$ represent the characteristic length and velocity scales of the vortex flow, which are used to render the problem dimensionless. This steady, doubly periodic flow field can be viewed as a simple analytical model of turbulence at the Kolmogorov scale, and it has been successfully employed for investigating several fundamental aspects of non-cohesive particle–vortex interactions (Bergougnoux, Bouchet, Lopez, & Guazzelli, 2014; Maxey, 1987).

Figure 5. (a) Streamlines of the doubly periodic background flow. (b) Typical floc configuration made up of spherical primary particles, with individual flocs distinguished by colour.

A typical floc configuration for $N_p = 50$ particles in a flow domain of size $L_x \times\, L_y = 2 \,\times\, 2$ is shown in figure 5(b). All primary particles have identical diameter $D_p = 0.1$ and density $\rho _s = 1$, with $St = 0.1$ and $Co = 5 \,\times\, 10^{-4}$. Details regarding the non-dimensionalization are provided in Zhao et al. (2020). Initially the particles are at rest and randomly distributed throughout the domain. When the distance between two particles is less than $\lambda /2$, we consider them as part of the same floc. We then track the number of flocs $N_f$ as a function of time, with an individual particle representing the smallest possible floc. Figure 6(a) demonstrates that the number of flocs decreases rapidly from its initial value $N_p$ due to flocculation, before levelling off around a constant value $N_{f,min}$ that reflects a stable balance between aggregation and breakage. We can fit the transient variation of $N_f(t)$ by an exponential function of the form

(4.2)\begin{equation} N_f(t) = (N_{p} - N_{f,min}) \,{\rm e}^{b t} + N_{f,min}, \end{equation}

where we evaluate $N_{f,min}$ as the average number of flocs during the equilibrium stage $50 \leqslant t \leqslant 200$. A characteristic flocculation time scale $t_{char}$ can then be defined as the time it takes for the number of flocs to decrease from its initial value $N_p$ to $N_{f,char} = N_p /2$. Hence, the corresponding characteristic time can be calculated as $t_{char} = {\rm ln} [(N_p /2 - N_{f,min}) / (N_p - N_{f,min})] / b$.

Figure 6. (a) Typical evolution of the number of flocs $N_f$ as a function of time. (b) Probability density function (PDF) of the floc size distribution during the equilibrium stage $50 \leqslant t \leqslant 200$. Simulation parameters are $N_p = 50$, $D_p = 0.1$, $\rho _s = 1$, $St = 0.1$ and $Co = 5 \times 10^{-4}$.

Figure 6(b) shows the statistical floc size distribution during the equilibrium stage $50 \leqslant t \leqslant 200$, where the floc size $N_{p,local}$ denotes the number of particles in a floc. Here $N_{f,local}$ refers to the number of the flocs of the same size. We find that the floc size distribution closely follows a log-normal distribution, consistent with previous experimental observations (Bouyer et al., 2004; Hill, Boss, Newgard, Law, & Milligan, 2011; Verney, Lafite, Burn-Cottan, & le Hir, 2011b).

Zhao et al. (2020) present a detailed discussion of the flocculation dynamics as a function of the governing dimensionless parameters. In particular, they observe that the number of flocs during the equilibrium stage $N_{f,min}$ decreases for increasing $Co$, along with the characteristic flocculation time. They find that the flocculation time $t_{char}$ has a minimum around $St \sim O(1)$, which reflects the well-known optimal coupling between particle and fluid motion when the particle response time is of the same order as the characteristic time scale of the flow (Wang & Maxey, 1993). This value is larger than the critical value $St_c=1/(8{\rm \pi} )=0.04$ derived by Massot (2007) that marks the transition for which particle inertia allows cohesionless grains to cross between different cellular flow fields (Yao & Capecelatro, 2018). Under these conditions, particles rapidly accumulate near the edges of the vortices, which facilitates the formation of flocs. The authors furthermore discuss the average number of primary particles per floc, the mean floc size and the associated fractal dimension of the floc. These quantities are useful for parameterizing the terms in population balance models, cf. the pioneering work by Levich (1962) and current extensions (Maggi et al., 2007; Shin et al., 2015; Verney et al., 2011b) and simplified versions (Lee, Toorman, Molz, & Wang, 2011; Shen, Lee, Fettweis, & Toorman, 2018; Son & Hsu, 2008, 2009; Winterwerp, 1998).

4.1.2 Flocculation in homogeneous isotropic turbulence

As a next step, we consider three-way coupled simulations with 10 000 primary particles in homogeneous isotropic turbulence (Zhao et al., 2021). As discussed by those authors, the problem is fully characterized by a turbulent Reynolds number $Re$, a characteristic parameter of the random turbulent forcing process $D_s$, the dimensionless particle diameter $D_p$, the density ratio $\rho _s$ and the cohesive number $Co$. Here $Re$ and $D_s$ can equivalently be expressed by the shear rate $G$ of the turbulence, and a suitably defined particle Stokes number $St$. The ratio of the Kolmogorov length scale to the primary particle size takes values up to 3.3 in the simulations of Zhao et al. (2021). Based on the observations by Bosse, Kleiser, and Meiburg (2015) that particle loading can modify the turbulence statistics even for volume fractions as low as $10^{-5}$, we expect two-way coupling effects to have an impact on the flocculation dynamics even in moderately dilute flows. Furthermore, even for globally dilute flows the local volume fraction inside a floc will be $O(1)$, resulting in significant shielding effects, so that the one-way coupled assumption generally will not hold within a floc. Nevertheless, as we will discuss in the following, three-way coupled simulations such as those by Zhao et al. (2021) are able to reproduce several experimentally observed statistical features of the flocculation dynamics.

Figure 7(a) shows the evolution of the number of flocs $N_f(t)$ with time for a representative case with $Co = 1.2 \times 10^{-7}$, $St=0.06$, $G = 0.62$, $\rho _s = 2.65$ and $\eta /D_p = 2.25$. Similar to the above cellular flow simulations of Zhao et al. (2020), $N_f$ decreases rapidly with time from its initial value, before levelling off around a constant value $N_{f,eq}$ that reflects a stable equilibrium between aggregation and breakage. Figure 7(b) shows separately the number of flocs with $N_{p} = 1, 2, 3$ and more than three primary particles. While the number of two- and three-particle flocs initially grows quickly, it soon reaches a peak and subsequently declines, as more flocs of larger sizes form.

Figure 7. (a) Temporal evolution of the number of flocs $N_f$. The vertical dashed line divides the simulation into the flocculation and equilibrium stages. (b) Number of flocs containing $N_p$ primary particles. The number of flocs with a single particle rapidly decreases from its initial value of $N_f = 10\,000$. The numbers of flocs with two or three particles intially grow and subsequently decay, as increasingly many flocs with three or more particles form. All results are for $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $G = 0.62$, $\rho _s = 2.65$ and $\eta /D_p = 2.25$.

To obtain insight into the dynamics of floc growth and breakage, it is useful to keep track of those flocs that maintain their identity over a suitably specified time interval $\Delta T$, those that add additional primary particles while keeping all of their original ones, and all those that have undergone a breakage event during the time interval. Zhao et al. (2021) present detailed results in this regard, which show that all three fractions reach statistically steady states.

The simulations furthermore enable us to track the floc size $D_f$ as the Feret diameter and its fractal dimension $n_f$ with time. Figure 8 shows the evolution of $D_f$ and $n_f$ over time for a representative floc. During the time interval $200 \leqslant t \leqslant 210$, hydrodynamic forces deform the floc so that it becomes more compact, which reduces $D_f$ and increases $n_f$. Later on, near $t \approx 240$, the floc is being stretched, with corresponding changes in $D_f$ and $n_f$.

Figure 8. Temporal evolution of the characteristic diameter $D_f$ and the fractal dimension $n_f$ of a typical floc that maintains its identity over the time interval considered. Three instants are marked by vertical dashed lines, and the corresponding floc shapes are shown. In response to the fluid forces acting on it, the floc first changes from a slightly elongated to a more compact shape, and subsequently to a more strongly elongated one. The floc with seven primary particles is taken from the case with governing parameters $Co = 1.2 \times 10^{-7}$, $St = 0.1$, $G = 0.91$.

Figure 9 shows representative results for different Stokes numbers $St$ of the primary particles. Figure 9(a) indicates that the equilibrium value of $\bar N_{p}$ increases for smaller $St$. The evolution of the average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles, shown in figure 9(b), demonstrates that smaller Stokes numbers result in more compact flocs. By varying the governing parameters, Zhao et al. (2021) find that, as a general trend, weaker turbulence, lower Stokes numbers and higher cohesive numbers result in larger and more compact flocs.

Figure 9. Temporal evolution of (a) the average number of primary particles per floc $\bar N_p$, and (b) the average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles, for different Stokes number values $St$, with $Co = 1.2 \times 10^{-7}$, $G = 0.91$ and $\eta /D_p = 1.85$. During the equilibrium stage, the number of primary particles per floc and the fractal dimension increase for smaller Stokes numbers.

In order to discuss the floc size distribution during the equilibrium stage, we sort the flocs into bins of width $\Delta (D_f/D_p) = 0.7$. Figure 10(a) shows that, for all values of the turbulent shear $G$, the size distribution peaks at the smallest flocs and then decreases exponentially with the floc size. The decay rate is largest for the strongest turbulence, confirming our earlier observation that strong turbulence breaks up large flocs and reduces the average floc size. This finding is consistent with the experimental observations by Braithwaite et al. (2012) in an energetic tidal channel. Corresponding results for different $St$ values display a similar trend.

Figure 10. Floc size distribution during the equilibrium stage, obtained by sorting all flocs into bins of constant width $\Delta (D_f/D_p) = 0.7$. (a) Results for different shear rates $G$, with $Co = 1.2 \times 10^{-7}$ and $St = 0.06$, during the time interval $1000 \leqslant t \leqslant 4000$. (b) Results for different cohesive numbers $Co$, with $St = 0.02$ and $G = 0.29$, for the time interval $15\,000 \leqslant t \leqslant 19\,000$.

We refer the reader to Zhao et al. (2021) for further detailed information regarding the preferential alignment of the flocs with the principal strain directions of the turbulent velocity field, as well as with the local vorticity vector. The authors observe that elongated flocs and the vorticity vector are strongly aligned with the direction corresponding to the largest eigenvalue of the Lagrangian stretching tensor. This alignment is consistent with corresponding previous findings by Parsa et al. (2011) and Ni, Ouellette, and Voth (2014).

Zhao et al. (2021) furthermore discuss how information obtained from the simulations can be employed towards the development of flocculation models. Based on their observation that the diameter and the average fractal dimension of the flocs are typically related by a power law expression, the authors proceed to fit the exponent of this relationship from the simulation results, which avoids the assumption of a constant fractal dimension in earlier investigations (Kuprenas et al., 2018; Winterwerp, 1998; Zhao et al., 2020). Consequently, the authors are able to formulate a new flocculation model, with variable fractal dimension.

In order to relate the cellular flow results to those for homogeneous isotropic turbulence, we note that one can view the cellular fields described in § 4.1.1 (Zhao et al., 2020) such that each cell represents an idealized eddy with the size of the Kolmogorov length scale $L_0$ and a maximum velocity $U_0$. For the case of homogeneous isotropic turbulence, the turbulent eddy size is set by the characteristic time scale $T_0=L_0/U_0$ that defines the magnitude of the force term $\mathbf {f}_t$ entering (3.1). Using this framework, Zhao et al. (2020) and Zhao et al. (2021) consider primary particles with diameter $D_p = 5 \,\mathrm {\mu }{\rm m}$ and demonstrated good agreement with the experiments of Kuprenas et al. (2018) and Maggi et al. (2007), respectively. For the case of homogeneous isotropic turbulence, Zhao et al. (2021) selected $T_0 = 7.81 \times 10^{-5}\,{\rm s}$ for the random turbulent forcing process. By choosing $L_0$, $T_0$ and $\rho _f = 1000\,{\rm kg}\,{\rm m}^{-3}$ as the characteristic length, time and density scales, one can obtain the characteristic velocity scale $U_0 = L_0/T_0 = 8\,{\rm m}\,{\rm s}^{-1}$, which is similar to values employed in previous investigations (Chen & Li, 2020; Chen et al., 2019). The definition of (3.36) yields much smaller values of the cohesive number. For example, Zhao et al. (2021) reported values in the range $10^{-10}\leq Co \leq 10^{-7}$. A conversion between the two definitions of the cohesive number (3.35) and (3.36) requires a choice for the value of the gravitational acceleration. For example, the standard value of $g=9.81$ m s$^{-2}$ yields a conversion factor of $3.7\times 10^8$.

4.2 Particle-resolved simulations of cohesive grain deposits

A key advantage of the method laid out in § 3.4.3 is that it is compatible with the grain-resolved simulation approach described in § 3.2.2. Such simulations provide high-fidelity data of cohesive sediment dynamics resolving the forces from fluid–particle and particle–particle interactions on every individual particle without any closures. Such simulations are costly and, hence, only suitable for small-scale applications, but they provide valuable insights especially for denser systems that are relevant, e.g. scenarios of hindered settling and granular collapse.

4.2.1 Hindered settling

In order to explore the influence of cohesive forces on the sedimentation process, Vowinckel et al. (2019a, 2019b) simulated the settling of an ensemble of 1261 polydisperse particles for the experimental conditions reported by te Slaa et al. (2015). The geometry of the particles was fully resolved, so that the fluid–particle coupling does not require any closures (§ 3.2.2). The particle–particle interactions were modelled the same way as described in § 4.1. A polydisperse mixture with a homogeneous particle volume fraction of $15\,\%$ was placed in a tank of quiescent fluid. To characterize the physical problem under consideration, it is convenient to define a reference velocity based on the buoyancy of a representative particle, which yields $u_{b}=\sqrt {g'D_{50}}$, where $g'=(\rho _p-\rho _f)g/\rho _f$ and $D_{50}$ denotes the median grain size of the particle size distribution. The characteristic time scale based on the buoyancy velocity and the median diameter then becomes $\tau _s = D_{50} / u_{b}$. Using a density ratio of $\rho _p/\rho _f=2.6$, the particle Reynolds number becomes $Re=u_b D_{50}/\nu _f = 1.35$. The grain sizes obey a cumulative log-normal distribution, with a maximum size ratio of ${\max \{D\} / \min \{D\}=4}$. The computational domain is of size $L_x \times L_y \times L_z = 13.1 D_{50} \times 40.0 D_{50}\times 13.1 D_{50}$, with gravity pointing in the negative $y$ direction. We assume periodic boundary conditions in the $x$ and $z$ directions, respectively, along with a no-slip condition at the bottom wall and a free-slip condition at the top wall. The median particle size is discretized by $D_{50}/h = 18.25$ grid cells. Two different values of the Cohesion number ${Co}$ are considered: (i) cohesionless grains with ${Co}=\max (\| \boldsymbol {F}_{{coh},50} \|) / m_{50} g'=0$, and (ii) strongly cohesive sediment with ${Co} = 5$. For both simulations, the particles are released from rest in quiescent fluid, and subsequently settle under the influence of gravity. Collisions are inelastic with $e_{p} = 0.97 < 1$, and subject to friction according to (3.18). Full details regarding the simulation set-up can be found in Vowinckel et al. (2019a).

The impact of cohesive forces on the settling behaviour is illustrated by figure 11. During the early stages the particle configurations are very similar for both simulations. Over the course of the simulations, however, the cohesive sediment is seen to settle faster than its non-cohesive counterpart. This qualitative observation is confirmed by the horizontally averaged concentration profiles shown as contours in figure 12, as a function of time.

Figure 11. Particle configurations during the settling process for $Co=0$ (a) and $Co=5$ (b). Left column: $t=17.6 \tau _s$, which corresponds to the time at which the particle phase has its maximum kinetic energy. From left to right, the columns are separated by time intervals of $72.5 \tau _s$. The particle colouring reflects the vertical particle velocity. The cohesive sediment is seen to settle more rapidly than its non-cohesive counterpart.

Figure 12. Contours of the horizontally averaged particle volume fraction $\phi _v$ over time: (a) cohesionless sediment and (b) cohesive sediment (image taken from Vowinckel et al., 2019a).

In agreement with the experimental observations of Been and Sills (1981), the volume fraction contours for the two simulations are nearly identical until $t = 120\tau _s$, as cohesive particles have not yet had sufficient time to flocculate. Subsequently, two distinct fronts emerge 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. 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 and marks the transition between hindered settling and consolidation (Winterwerp, 2002). The cohesionless sediment, on the other hand, has not yet reached the point of contraction by the end of the simulation.

As the particles settle, they displace fluid at the bottom of the tank and generate an upward counterflow that is sufficiently strong in the current simulations to sweep smaller particles upward. This represents one reason for the reduced settling velocity of the smaller grains, and for the separation of the different grain sizes in very tall water columns (te Slaa et al., 2015). Figure 13 illustrates this effect by showing the final volume fraction profiles for the smallest, intermediate and largest third of the particles. The figure demonstrates that small and intermediate cohesive particles settle much more rapidly than their non-cohesive counterparts, consistent with the observation of Lick et al. (2004), who found intermediate size particles to be most strongly affected by cohesive forces. These types of grains have their peak concentration in the interval $3 D_{50} < y < 10 D_{50}$. In contrast, large cohesive grains have a lower volume fraction near the bottom of the tank than large cohesionless grains (figure 13c). As a consequence, the effect of size segregation is less pronounced for cohesive grains, as flocculation results in particles of different sizes settling with the same velocity (Mehta et al., 1989).

Figure 13. Particle volume fraction profile for different particle radii at $t=480 \tau _s$: (a) small particles with $D_p \leq D_{33}$, (b) medium sized particles in the range $D_{33} < D_p \leq D_{66}$, and (c) large particles with $D > D_{66}$. Note the different horizontal axis scalings for the individual frames. The results in (a,b) were smoothed by a moving average with a filter width of $1.5 D_{50}$ for clarity.

Vowinckel et al. (2019b) compared their simulation results for the effect of cohesion on the settling velocity of polydisperse particles to established empirical relations for silt. They computed the settling behaviour of the particle phase by means of an averaging operator applied to the Eulerian fluid grid that evaluates instantaneous snapshots of the particle velocity distribution $(\phi _v, v_p)$ to obtain horizontally averaged values of $\overline {\langle v_p \rangle }$ and the settling velocity $v_\infty$ (Clift et al., 2005) of a single particle in still fluid at the same particle Reynolds number, where the angular brackets denote horizontal averaging and $\phi _v$ is the part of a cell occupied by solids and $v_p$ is the settling velocity of the particle taking up this space. Note that this quantity still depends on the particle diameter. Vowinckel et al. (2019b) computed its value in an iterative fashion using the force balance of a particle and empirical correlations for the drag coefficients. The ratio $\overline {\langle v_p \rangle }/v_\infty$ is still a function of the volume fraction and can therefore be used to compare with the hindered settling functions available in the literature.

One of the first hindered settling functions was proposed by Richardson and Zaki (1954),

(4.3)\begin{equation} \frac{\overline{\langle v_p \rangle}}{v_{\infty} } = \left(1-\frac{\phi}{\phi_s}\right)^n, \end{equation}

who experimentally investigated the settling behaviour of non-cohesive grains. Here, $\phi _s$ and $n$ are empirical parameters describing the volume fraction of a freshly deposited sediment bed and the particle size and shape, respectively. As argued by Dankers and Winterwerp (2007), cohesive sediment such as mud with a significant amount of clay deposits at the bottom in a gel-like structure with a volume fraction that is lower than the maximum possible volume fraction. Hence, these authors introduced another critical volume fraction $\phi _{max}$ with the property $\phi _s<\phi _{max}$, where the transition from $\phi _s$ to $\phi _{max}$ is governed by gel collapse and self-weight consolidation rather than settling. Due to its simplicity, (4.3) has been very popular in hydraulic engineering. However, as described by te Slaa et al. (2015), this function is known to underestimate the settling velocity for higher concentrations. Instead, these authors have used the hindered settling function of Winterwerp (2002),

(4.4)\begin{equation} \frac{\overline{\langle v_p \rangle}}{v_\infty} = \frac{\left(1-\dfrac{\phi}{\phi_s}\right)^m (1-\phi)}{\left(1-\dfrac{\phi}{\phi_{max}}\right)^{-({5}/{2})\phi_{max}}}, \end{equation}

where the numerator represents the effects of the counterflow and the buoyancy, respectively, while the denominator reflects the increased viscosity of dense suspensions according to Krieger and Dougherty (1959). Here, $m$ is an empirical exponent, which plays a similar role to the parameter $n$ in (4.3). It was shown by te Slaa et al. (2015) that (4.3) and (4.4) both yield good agreement with experimental results for settling coarse silt.

Vowinckel et al. (2019a) compared their simulation results to the hindered settling functions as parameterized by te Slaa et al. (2015) in figure 14(a) to validate their simulations. The results agreed well with the two empirical hindered settling functions, and they demonstrated the enhanced settling velocity due to the cohesive forces for all concentration values. It was shown by Vowinckel et al. (2019a) that choosing $\phi _s = \phi _{max}$ yields good results. This is in contrast to the observations of Dankers and Winterwerp (2007), who found for mud that $\phi _s<\phi _{max}$, but since the simulations of Vowinckel et al. (2019b) deal with coarse silt, consolidation does not play a large role and it is justified to choose $\phi _s=\phi _{max}$. The consolidation of the sediment will continue to squeeze out water from the bed until the particle packing jams. This process will maintain a counterflow over very long time scales (e.g. Houssais, Ortiz, Durian, & Jerolmack, 2016). Based on these observations, Vowinckel et al. (2019b) propose to not calibrate $\phi _s$ but to parameterize critical volumetric concentrations by the maximum value of their concentration data. Using this data, they can immediately parameterize $\phi _s = \phi _{max}=0.7$, which is in line with experimental and computational studies of polydisperse particle packings (Desmond & Weeks, 2014; Sohn & Moreland, 1968). Note that this value is larger than the volume fraction of a randomly closed packing of monodisperse spheres, since we are dealing with polydisperse particles and the employed grid fully resolves the volume fraction over intervals smaller than the particle diameter. Using this physically based parameterization, one obtains a much better fit of (4.4) to the data (dotted line in figure 14b), which illustrates the importance of including the effects of the counterflow, of buoyancy, and of the increased viscosity in the formulation of the hindered settling function, even during the consolidation phase. On the other hand, (4.3) underestimates the settling velocity of cohesive sediment, but is very close to the settling behaviour of the simulated cohesionless sediment. This was expected, as (4.3) was derived for cohesionless sediments in the first place.

Figure 14. Settling velocity normalized with the undisturbed settling velocity. Comparison to empirical relationships (4.3) of Richardson and Zaki (1954) (RZ) and (4.4) of Winterwerp (2002) (W).

4.2.2 Submerged cohesive granular collapse

Submerged non-cohesive granular collapses can be classified into three different categories (free fall, inertial and viscous), depending on the particle and fluid properties (Bougouin & Lacaze, 2018; Jing, Yang, Kwok, & Sobral, 2019). The relevant regime is determined by two dimensionless numbers, the Stokes number $St$ and the density ratio $r$,

(4.5)$$\begin{gather} St = \frac{1}{18\sqrt{2}}\frac{\sqrt{\rho_p(\rho_p-\rho_f)gD_p^3}}{\nu_f\rho_f}, \end{gather}$$

(4.6)$$\begin{gather}r = \sqrt{\frac{\rho_p}{\rho_f}}. \end{gather}$$

For $St\gg 10$ and $r\gg 4$, the collapse is in the free-fall regime, for $St\gg 2.5r$ and $r\ll 4$, it is in the inertial regime, and otherwise, it is in the viscous regime. For cohesive collapses, an additional parameter $Co$ arises in the form of (3.35).

Using the method by Vowinckel et al. (2019b) summarized in § 4.2.1 above, Zhu et al. (2022) discuss particle-resolved cohesive simulation results for both shallow and tall columns, with representative results shown in figure 15. Consistent with earlier observations for non-cohesive collapses (Lee, Huang, & Yu, 2018; Meruane, Tamburrino, & Roche, 2010), cohesive collapses proceed through an acceleration or collapse stage, a constant front velocity stage and a deceleration stage. More cohesive collapses travel over shorter distances than their less cohesive counterparts, and with a smaller front velocity. For $Co=50$ with shallow columns and $Co=250$ with tall columns, the cohesive force is sufficiently large so that the column no longer collapses. This finding is qualitatively similar to previous experimental and numerical observations for dry cohesive granular collapse (Artoni, Santomaso, Gabrieli, Tono, & Cola, 2013; Langlois, Quiquerez, & Allemand, 2015; Santomaso, Volpato, & Gabrieli, 2018).

Figure 15. Evolution of granular collapse for $a=1$ and $Co=0$ and 30. From left to right: $t=5$, 20 and 30. (a–c) Magnitude of the angular velocity $\|\boldsymbol \omega _p\|$ for $Co=0$; (d–f) magnitude of the translational velocity $\|\boldsymbol {u_p}\|$ for $Co=0$; (g–i) $\|\boldsymbol \omega _p\|$ for $Co=30$; (j–l) $\|\boldsymbol {u_p}\|$ for $Co=30$. The red lines indicate the location of the failure surface. The black arrows represent vectors of the average particle velocity.

The angular and translational particle velocity magnitudes are displayed in figure 15. Figure 15(d,j) shows that the velocity magnitude in the $x$ and $y$ directions remains very small in the lower left corner of the columns, as the particles remain approximately at rest there. For both cohesive and non-cohesive (Sun, Zhang, Wang, & Liu, 2020) columns, during the acceleration stage particles near the upper right corner slide down along an inclined failure surface (indicated by a red line in frames 15d,j). The failure surface is defined as the contour where $\|\boldsymbol {u_p}\| = 0.05\|\boldsymbol {u_p}\|_{max}$ (Lacaze & Kerswell, 2009), with $\|\boldsymbol {u_p}\|_{max}$ denoting the maximum translational velocity at the same time. By comparing figure 15(d,j), we note that cohesive forces corresponding to $Co=30$ elevate the location of the failure surface, resulting in the growth of the region of stationary particles in the lower left corner. The angular velocity of the particles remains quite small near the failure surface (cf. frames 15a,b,g,h), which indicates that the particles slide, rather than roll, past each other (Xu, Dong, & Ding, 2019).

Based on the above results, Zhu et al. (2022) are able to establish scaling laws for the quasisteady front velocity and runout length as functions of $a$ and $Co$. They proceed to analyse the deposit morphology as well as the energy budget of the granular collapse events, in terms of the potential and kinetic energy components of the fluid and particulate phases, respectively.

The particle-resolving simulations provide complete information on the evolution of the individual cohesive force bonds, so that these can be analysed in detail. For the initial configurations with $Co=10$ as well as $a=1$ and 8.6, figure 16(a,c) indicates all cohesive bonds between individual particles by straight line segments that connect the particle centres. Figure 16(b,d) shows those of the initial cohesive bonds that have survived until the end of the collapse. We observe a few clear differences between the two aspect ratios. As discussed earlier, for $a=1$, the particles in the lower left corner hardly move at all, so that many of the initial cohesive bonds between them survive the entire collapse process. Near the left wall, and in the entire upper (pink) layer, quite a few of the bonds also survive, whereas this is not true for the lower sections of the deposit near the front. For the particles in that section, the collapse process destroys most of their initial cohesive bonds. For $a=8.6$, on the other hand, cohesive bonds primarily survive in a very small section in the lower left corner, and along the entire top layer of the deposit, including at the very front. For the entire interior section of the final deposit, almost all initial cohesive bonds are destroyed during the course of the collapse.

Figure 16. Cohesive bonds at the initial time for $a=1$, $Co=10$ (a) and $a=8.6$, $Co=25$ (c), respectively. For the same two flows, (b,d) shows those cohesive bonds that have stayed intact during the entire collapse process until the final time $t=60$.

5. Summary and outlook

While the last couple of decades have seen substantial improvements in our understanding of liquid and gaseous flows with non-cohesive particles, comparable advances have not yet been achieved for cohesive particulate flows. Progress in this regard will be essential for enhancing our ability to model the transport of mud in rivers, estuaries and coastal oceans, and for predicting carbon and nutrient fluxes in these aquatic environments. Similar considerations apply to subaerial flows including dust storms, dune migration and snow avalanches, as well as to industrial processes involving powders, or to helicopter downwashes and rocket exhaust plumes in desert environments. Advances in these directions will furthermore enable us to better capture the dynamics of complex media such as soil, in which inter-particle forces vary as a function of environmental conditions such as humidity, salinity and temperature. It is important to note that watery mixtures with a large portion of fine-grained sediments potentially contain some organic material. Such mixtures are typically referred to as mud (Berlamont, Ockenden, Toorman, & Winterwerp, 1993). These organic substances frequently form a biofilm coating of extracellular polymeric substances (Lai, Fang, Huang, He, & Reible, 2018), which may give rise to a range of additional complex phenomena (Gerbersdorf et al., 2020) that have not been discussed in the present review but should be considered in future research.

Grain-resolved simulations, i.e. simulations with meshes fine enough to geometrically resolve the flow around individual particles, represent a powerful tool for advancing our understanding of the dynamics of cohesive materials, and for achieving progress on the above issues. Since they provide information on the forces and stresses at the individual grain scale, along with the resulting motion, they will enable us to quantify the erodibility and to develop rheological models of such materials, with the ultimate goal of deriving reliable scaling laws for their dynamics in terms of the macroscopic friction, a framework that has previously been derived for non-cohesive particulate flows (Boyer, Guazzelli, & Pouliquen, 2011; Guazzelli & Pouliquen, 2018). Furthermore, the information obtained from such simulations can feed into tools developed in related research areas such as network theory, in order to advance our understanding of microstructure of suspended particles at a fundamental level (Papadopoulos, Porter, Daniels, & Bassett, 2018; Papadopoulos, Puckett, Daniels, & Bassett, 2016). An additional attractive pathway towards leveraging the highly resolved data of grain-resolved simulations is to generate high-fidelity inputs that can be used as features for machine learning strategies. Such approaches have shown promise in terms of providing efficient upscaling strategies in complex flows (Moore, Balachandar, & Akiki, 2019; Seyed-Ahmadi & Wachs, 2020).

In this regard, it is very beneficial that grain-resolved simulations allow us to carefully control the initial and boundary conditions of specific configurations, such as their initial volume fraction, the regularity of the initial particle arrangements, the relative strength of cohesive and gravitational forces, and their grain size distributions and shapes, along with the surface roughness and the fluid/particle density ratio. At the same time, additional progress on several fronts is required in order to further increase the potential impact of grain-resolving simulations. Examples in this regard are an improved understanding of the fluid-mediated interaction of individual grains at the microscale and nanoscale under different ambient conditions; the further refinement of numerical multiscale approaches that accurately account for very short particle–particle interactions while also capturing the long time scales over which continuum behaviour emerges; and the formulation of upscaling approaches that yield accurate continuum models for use in such frameworks as direct numerical simulations and large eddy simulations. All of these research directions will require the careful integration of theory, numerical analysis, laboratory experiments and field observations. In this regard, the recently developed capability to tune the cohesion of particles in laboratory experiments offers interesting opportunities (Sharma et al., 2022).

In addition, as reviewed by Vowinckel (2021), recent advances in the field of pr-DNS have introduced techniques to relax the assumption of spherical particles (Ardekani & Brandt, 2019; Eshghinejadfard, Zhao, & Thévenin, 2018; Jain, Tschisgale, & Fröhlich, 2021; Voth & Soldati, 2017). Such considerations could further complicate the picture of cohesive sediment dynamics with important implications that are yet to be investigated.

In this way, it is hoped that such simulations will eventually enable us to address a wide variety of complex scenarios. A prime example in this regard is the dynamics of clay/sand suspensions. Due to the difference in the size of sand and clay particles, and the resulting large range of temporal and spatial scales over which they interact, such problems are truly multiscale in nature. A major challenge will be to unravel the macroscopic behaviour of such materials, including the emergence of a yield stress, the existence of a fluid-to-solid-like transition, and the appearance of microscopic precursors to the catastrophic failure of load-bearing structures. The ability to predict such phenomena will be key to controlling the structure and motion of highly heterogeneous particulate systems. Another important application is the formation of porous flocs as the smallest entities that are transported by the flow. Such flocs behave very differently from rigid, non-porous particles. This is especially true for flocs moving through stratified ambient fluids, a scenario that is highly relevant for the behaviour of marine snow. The porous structure of the flocs, however, will require a coupled-continuum approach where the particle interior is treated as a porous medium, whereas the exterior is still governed by the full Navier–Stokes equations.

Tackling these complex issues may require hybrid computational approaches, in which only the large grains are explicitly resolved, while the smaller ones are parameterized by constitutive equations. Those equations, however, are still largely unknown, and further grain-resolved investigations of cohesive fine particles will be required for their formulation.