Physics of Cohesive Sediment Flocculation and Transport: State-of-the-Art Experimental and Numerical Techniques

Abstract

Due to climate change, sea level rise and anthropogenic development, coastal communities have been facing increasing threats from flooding, land loss, and deterioration of water quality, to name just a few. Most of these pressing problems are directly or indirectly associated with the transport of cohesive fine-grained sediments that form porous aggregates of particles, called flocs. Through their complex structures, flocs are vehicles for the transport of organic carbon, nutrients, and contaminants. Most coastal/estuarine models neglect the flocculation process, which poses a considerable limitation of their predictive capability. We describe a set of experimental and numerical tools that represent the state-of-the-art and can, if combined properly, yield answers to many of the aforementioned issues. In particular, we cover floc measurement techniques and strategies for grain-resolving simulations that can be used as an accurate and efficient means to generate highly-resolved data under idealized conditions. These data feed into continuum models in terms of population balance equations to describe the temporal evolution of flocs. The combined approach allows for a comprehensive investigation across the scales of individual particles, turbulence and the bottom boundary layer to gain a better understanding of the fundamental dynamics of flocculation and their impact on fine-grained sediment transport.

1. Introduction

Cohesive sediment transport is a tightly coupled system driven by hydrodynamic forcing, resuspension, deposition, and flocculation. Turbulence intensity in the bottom boundary layers and the water column mixed layers controlled by currents (e.g., river outflows, tidal currents), surface waves, and estuarine circulations can vary by several orders of magnitude in terms of turbulent dissipation rate. Hence, the turbulent shear rate defined as $G=\sqrt{\epsilon /\nu }$ is the most important parameter for flocculation models. In the estuarine tidal boundary layers, the turbulent dissipation rate is around $\mathcal{O}\left({10}^{-6}\sim {10}^{-4}\right)m^2/s^3$ ($G=1\sim 10s^{-1}$; e.g., [1, 2]). Near estuarine fronts with shear instabilities, turbulent dissipation rates can increase up to $\mathcal{O}\left({10}^{-3}\right)m^2/s^3$ ($G=30s^{-1}$; [3]). During storm condition when wave motion become dominant, turbulent dissipation rates in the thin wave bottom boundary layer can exceed $\mathcal{O}\left({10}^{-3}\right)$ ($G=30\sim 100s^{-1}$, [4, 5]). Under breaking waves in the upper ocean or in the surf zone, the turbulent dissipation rate can be as high as $\mathcal{O}\left({10}^{-3}\sim {10}^{-2}\right)$$(G=100\sim 200s^{-1}$, [6, 7]).

Computer simulation models are commonly the chosen tools that coastal managers use to predict sediment transport rates. In these continuum models, the fluid motion is computed via the continuity equation and the Navier-Stokes equations or spatially and temporally averaged variants thereof, whereas the sediment is represented as a concentration field [8]. This approach allows to simulate large spatial scales covering entire estuaries, because the governing equations are solved on the same Eulerian grid. However, these types of models require closures to account for the unresolved physics of the sediment dynamics, in particular vertical sedimentary distributions [9, 10] and mass fluxes. The latter is the product of the concentration and the settling velocity. Manning and Bass [11] found that mass settling fluxes can vary over four or five orders of magnitude during a tidal cycle in mesotidal and macrotidal estuaries; therefore, a realistic representation of flux variations is crucial to an accurate depositional model.

The specification of the flocculation term within numerical models depends on the sophistication of the model structure. Until recently, even the conceptual relationship between floc size, suspended particulate matter (SPM) concentration and turbulent shear stress proposed by Dyer [12] (see Figure 1a) remained largely unproven. Hence, much more work needs to be done in order to arrive at robust simulation tools with predictive capacity. The present contribution, therefore, provides a review on the state-of-the-art of floc measurements in both the field and laboratory. In addition, we review a newly emerging technique of particle-resolved simulations that can provide a promising alternative avenue to generate data of small-scale sediment dynamics to derive the missing constitutive equations for continuum models. Finally, we provide a summary of current techniques that are used in continuum models to account for cohesive sediment dynamics, where we explicitly point to the components that deserve more research in the future.

Figure 1.

An illustration of the dependence of floc size on shear stress and SPM concentration is given. (a) A conceptual diagram showing the relationship between floc modal diameter, suspended sediment concentration and shear stress (from [12]). Subplots (b–e): illustrative examples of real estuarine floc images showing ambient shear stress, SPM concentration and settling velocity; (b) shear stress 0.3 N ${\mathrm{m}}^{-2}$, SPM concentration 3.5 g ${\mathrm{l}}^{-1}$, and settling velocity 5 mm ${\mathrm{s}}^{-1}$; (c) shear stress 0.3 N ${\mathrm{m}}^{-2}$, SPM concentration 3.5 g ${\mathrm{l}}^{-1}$, and settling velocity 8 mm ${\mathrm{s}}^{-1}$; (d) shear stress 0.45 N ${\mathrm{m}}^{-2}$, SPM concentration 0.25 g ${\mathrm{l}}^{-1}$, and settling velocity 1.8 mm ${\mathrm{s}}^{-1}$; (e) shear stress 0.45 N ${\mathrm{m}}^{-2}$, SPM concentration 0.25 g ${\mathrm{l}}^{-1}$, and settling velocity 0.2 mm ${\mathrm{s}}^{-1}$. Illustrations (b–e) are modified from Manning.

2. Assessment of temporal floc evolution

2.1 Floc measurements

2.1.1 Direct floc size measurements

The presence of large estuarine macroflocs was initially observed in situ using underwater photography [13]. However, floc breakage occurs during sampling in response to the additional shear created by the instrumentation [14]. To overcome this problem, less-invasive techniques for measuring floc properties in situ have been developed. Usually, these can be divided into devices that solely measure floc size ($D$) e.g. Lasentec (Figure 2a) [15], LISST [16], LISST-Holo [17], and InSiPid [18]; and those devices that can provide measurements both of floc size and settling velocity ($W_s$) e.g. VIS [19], HR Wallingford video camera system [20], VIL [21], INSSEV: IN-Situ SEttling Velocity instrument [22, 23] (Figure 2b–c), LabSFLOC—Laboratory Spectral Flocculation Characteristics—instrument [25, 26, 27], INSSEV-LF [28], and PICS—Particle Imaging Camera System [29].

Figure 2.

Some examples are given of floc measuring instrumentation. (a) Schematic of the adapted Lasentec par-tec 100 probe unit showing: 1. The light guide, scanning mechanism and focusing lens, 2. The PVC cylinder, 3. Watertight cable termination, and 4. Electronic circuitry and power supplies (from [15]). (b, c) the INSSEV instrument (from [24]); (b) side view of INSSEV mounted on a metal deployment frame; (c) front view of INSSEV (right side of image), together with optical backscatter (OBS) sensors and an acoustic Doppler velocimeter (ADV) positioned on a vertical pole (left side of image). The ADV provides high frequency turbulence data that can be directly related to the floc populations. (d) Views of a LabSFLOC-2 together with a schematic illustration of the instrument after [25]).

The strength of video-based floc measurements is that they minimize the number of assumptions used during the data processing and interpretation stages. Devices that only measure the size component require additional gross and often incorrect assumptions regarding the relationship between settling velocity, floc size, and floc density. The settling velocity of a floc is a function of both its size and effective density, and both of these floc components can display variations spanning three to four orders of magnitude within any one floc population [30, 31, 32].

Of note, the LabSFLOC suite of high-resolution, low intrusive, underwater video camera systems for the past 20+ years have been regarded internationally as a benchmark device for the sampling and dynamical simultaneous measurement of sizes and settling velocities for entire populations of flocs and a range of bio-sedimentary particles, and this enables individual floc effective density values to be determined by applying a modified Stokes Law [33]. This also enables the calculation of additional floc properties including: structural composition (porosity, fractal dimension), shape, sedimentary mass flux, and floc population mass-balancing, all from within a wide range of aquatic environments. In conclusion, selection of the most appropriate instrumentation is paramount when attempting to parameterise flocculated cohesive sediments. Manning et al. [34] provide a detailed review of many of these floc measuring systems.

2.1.2 Parameterizing floc data

To aid the interpretation of floc characteristics and for their inclusion in sediment transport models, each floc population can be segregated into various subgroupings according to floc size. Sample-mean floc values can be computed (i.e., a single value per floc population) to show generalized floc property trends.

Dyer et al. [35] reported that a single mean or median settling velocity did not adequately represent an entire floc spectrum, especially in considerations of a flux to the bed. Dyer et al. [35] recommended that the best approach for accurately representing the settling characteristics of a floc population was to split a floc distribution into two or more components, each with their own mean settling velocity. Both Eisma [13] and Manning [32] concur with this finding by suggesting that a more realistic and accurate generalization of floc behavior can be derived from the macrofloc and microfloc fractions. These two floc fractions form part of Krone’s [36] classic order-of-aggregation theory and produce two floc property values per floc population.

Macroflocs are large (typically $D>160\mathrm{\mu m}$ [32]) highly porous (typically $>90\%$) and often fragile (and difficult to sample), fast-settling aggregates (see Figure 1b–d), and typically close to the size of the turbulent Kolmogorov microscale [37, 38]. They are recognized as the most important subgroup of flocs because their fast-settling velocities tend to have the strongest influence on the mass settling flux [39]. Macroflocs are progressively broken down as they pass through regions of higher turbulent shear stress and reduced again to their component microfloc substructures [40]; they rapidly attain equilibrium with the local turbulent environment.

The smaller microflocs (typically $D<160\mathrm{\mu m}$ [32]; see Figure 1e) are considered to be the building blocks from which the macroflocs are comprised. Many field studies [22, 41, 42, 43] have shown that the microfloc class of aggregates tend to display a much wider range of effective densities and settling velocities than the macrofloc fraction. Microflocs are much more resistant to breakup by turbulent shear; they tend to have slower settling velocities but exhibit a much wider range of effective densities than the larger macroflocs.

In terms of flocculation kinetics [44], the macroflocs tend to control the fate of purely muddy sediments in an estuary [45]; this is because the smaller microflocs generally settle at less than 1 mm ${\mathrm{s}}^{-1}$, whereas macroflocs settle in the 1–15 mm ${\mathrm{s}}^{-1}$ range, enabling them to deposit to the bed [46]. However, when flocculation of mixed sediments occurs, the microflocs can potentially demonstrate settling velocities comparable to those of the macroflocs [28, 47].

2.1.3 Floc data example

In order to illustrate the spectral variability of floc properties, each floc population can be divided into various size bands. The band divisions can be chosen to best fit the data collected. LabSFLOC data for a mud sample from the Medway Estuary, UK, provide a graphical illustration of size banding and the increasing settling velocities associated with larger flocs (Figure 3). Twelve size bands (SB) have been used to represent the Medway floc data, with SB1 representing microflocs less than $40\mathrm{\mu m}$ in size, whilst SB12 is representative of macroflocs greater than $640\mathrm{\mu m}$ in diameter; SB2 to SB6 range from 40 to $240\mathrm{\mu m}$ in five steps, each of $40\mathrm{\mu m}$, and SB7 to SB11 range from 240 to $640\mathrm{\mu m}$m in five steps of $80\mathrm{\mu m}$.

Figure 3.

An example of various size-banded floc properties for a LabSFLOC SPM sample is given. The complete population of size versus settling velocity data is illustrated in the top-left panel. These data apply to Medway estuary mud slurry (1.6 g ${\mathrm{l}}^{-1}$) that has been sheared in the Southampton oceanography Centre (UK) mini-flume at a shear stress of 0.37 nm^–2 (data from results in [43]).

2.1.4 Floc settling modeling approaches

2.1.4.1 Constant settling velocity

Specification of the flocculation term within numerical models depends on the sophistication of the model. The simplest parameterisation is a single floc settling velocity value that remains constant in both time and space (one coefficient). These fixed settling values are usually in the range of 0.5–5 mm ${\mathrm{s}}^{-1}$ and typically are selected on an arbitrary basis and sometimes used as a tuning parameter to match predicted erosion and deposition patterns to observations for an undisturbed estuary.

2.1.4.2 Power law settling velocity

The next step is to use gravimetric data provided by field settling-tube experiments to relate floc settling velocity to the instantaneous SPM concentration, using a power law with two coefficients (e.g., [48]; see Figure 4a). Empirical results have shown a generally exponential relationship between the mean, or median, floc settling velocity and SPM concentrations for concentrations $<10\mathrm{g}{\mathrm{l}}^{-1}$. This approach sometimes includes hindered settling (see Figure 4b). However, both the constant settling velocity and the power law parameterisation techniques do not include the important and influential effects of turbulence such as floc breakup induced by turbulent shear.

Figure 4.

Some examples of floc settling velocity measurements are shown. (a) Owen tube determinations of median settling velocity as a function of suspended sediment concentration for different estuaries; the bold dashed line represents an exponent of unity (reproduced with minor modifications from [49]). (b) Median settling velocity of Severn estuary mud as a function of SPM concentration; the Owen-tube data are taken from Odd and Roger [50], and the dashed line represents results based on the SandCalc software sediment-transport computational algorithm, which incorporates the hindered settling effect at high concentrations (reproduced with minor modifications from Soulsby [51]).

2.1.4.3 The van Leussen parameterisation

More recently, a number of authors have proposed simple theoretical formulae interrelating a number of floc characteristics that can then be calibrated using empirical studies. Such an approach has been used by van Leussen [19], who utilized a formula that modifies the floc settling velocity in still water by a floc growth factor, due to turbulence, and then reduces it by a turbulent floc disruption factor. The reference settling velocity (taken at low turbulent shear conditions), $W_{s_0}$, is then related to the SPM concentration ($C$) by a power law:

$$W_{s_0}=k\cdot C^m$$E1

where $k$ and $m$ are empirical constants. The van Leussen is a qualitative simplification of a model originally developed for the sanitation industry [52], with only a limited number of interrelated parameters, and hence does not provide a complete description of floc characteristics within a particular sheared environment.

2.1.4.4 The Lick et al. parameterisation

A number of authors have attempted to observe how the floc diameter changes in turbulent environments. In particular, Lick et al. [53] derived an empirical relationship based on laboratory measurements made in a flocculator. They found that the floc diameter varied as a function of the product of the SPM concentration and the turbulence parameter as the turbulent shear rate, $G$:

$$D=c{\left(C\cdot G\right)}^{-d}$$E2

where $c$ and $d$ are empirically determined values. However, this formulation provides no information on the important floc settling or floc dry mass properties.

2.1.4.5 The Manning and Dyer parameterization

The Manning Floc Settling Velocity (MFSV) algorithm for settling velocity [54] is based entirely on empirical observations made in situ using nonintrusive floc and turbulence data acquisition techniques in a wide range of estuarine conditions. The floc population size and settling velocity spectra were sampled using the video-based INSSEV instrument and LabSFLOC data.

The Manning-Dyer algorithms were generated by a parametric multiple regression statistical analysis of key parameters, which were generated from the raw, spectral floc data. Detailed derivations and preliminary testing of the floc-settling algorithms are described by Manning [26, 55]. Although the resulting empirical formulae are not presented in a fully dimensionless form, these formulae have the merit of being based on a large dataset of accurate, in situ settling velocity measurements (157 individually observed floc populations), acquired from different estuaries (Tamar, Gironde and Dollard) and different estuarine locations, such as the turbidity maximum and the intertidal zone.

The algorithms are based on the segregation of flocs into macroflocs ($D>160\mathrm{\mu m}$ [32]) and microflocs ($D<160\mathrm{\mu m}$), which comprise the constituent particles of the macroflocs. This distinction permits the discrete computation of the mass settling flux (MSF) at any point in an estuarine water column.

Equations are given for: the settling velocity of the macrofloc fraction ($W_{s,\mathit{\text{macro}}}$), the settling velocity of microflocs ($W_s$) and the ratio of macrofloc mass to microfloc mass in each floc population—termed the SPM ratio [55]. This type of formulation gives a good compromise between the representation of physicochemical processes and computational simplicity. Eqs (3) and (4) describe $W_{s,\mathit{\text{macro}}}$ ($\mathrm{mm}{\mathrm{s}}^{-1}$) with inputs of SPM in $\mathrm{mg}{\mathrm{l}}^{-1}$ and $\tau$ in Pa:

$$W_{s,\mathit{\text{macro}}}=0.644+0.000471\mathit{SPM}+9.36\tau -13.1{\tau }^2$$E3

for $0.04<\tau <0.6$ Pa

$$W_{s,\mathit{\text{macro}}}=3.96+0.000346\mathit{SPM}-4.38\tau +1.33{\tau }^2$$E4

for $0.6<\tau <1.5$ Pa. These equations require the input of a turbulent shear stress ($\tau$) and an SPM concentration. These regression equations provide a realistic approximation to the field data. Graphical representations of the equations, together with the data, are presented in Manning and Dyer ([54]; see Figure 5). The Manning settling algorithm is valid for SPM concentrations in the range 10–8600 mg ${\mathrm{l}}^{-1}$ and shear stress values of $\tau <2.13$ Pa, with extrapolation extending this range to 5–10 Pa.

Figure 5.

Illustration of the settling velocities of macroflocs and microflocs, plotted against shear stress, for a mixed sediment suspension comprising a ratio of 25 per cent mud to 75 per cent sand and a pure mud suspension, all for a total SPM concentration of 5 gl^–1 (modified from Figure 14 of Manning et al. [56]).

An example of this is the implementation of the algorithm in a TELEMAC-3D numerical model of the Thames Estuary, UK [57], in which it was shown that the use of the Manning algorithm greatly improved the reproduction of observed distributions of SPM concentrations compared with the other formulations, both in the vertical and horizontal dimensions.

The Manning settling algorithms have been extended to cater for mixed sediment flocculation settling, including different ratios of mud:sand ([24, 47, 56]; see Figure 5). These algorithms are a major step forward in establishing a reliable estimate of the settling velocity. It has been developed based on a large and reliable dataset, it caters for the spectrum of hydrodynamic conditions that occur during a typical tidal cycle [58] (a feature often lacking in the settling terms of many estuarine sediment transport models) and has been shown to more accurately reproduce the distribution of suspended sediment compared with simpler settling models.

Soulsby et al. [59] has developed a more ‘physics-based’ version of the empirical model based on the Manning-Dyer formulation, called Soulsby-Manning 2013. It should be noted that for flocculation algorithms and models that include turbulence as a contributing variable, it is vital to ensure that the turbulence data are accurate, otherwise it has significant implications for the accuracy of the calculated floc settling characteristics.

2.1.4.6 Complex population approaches

Lee et al. [60] applied a time-evolving two-class population balance equation (PBE) to determine the spatially and temporally changing distribution of fixed-size microflocs and size-varying macroflocs for bimodal floc distributions, with a fractal relationship between floc size and mass to derive the distribution of settling velocities. However, the authors felt that further intensive investigation of the aggregation and breakage kinetics would be required before their model was generally applicable when compared with the simpler approach of Manning and Dyer [54] and, presumably, Soulsby et al. [59].

Verney et al. [61] applied a time-evolving, multi-fraction model to determine the spatially and temporally changing distribution of the numbers of flocs in each size fraction, with a fractal relationship between floc size and mass to derive the distribution of settling velocities.

A relationship between the floc settling velocity and floc properties and fractal dimensions is given by Winterwerp [62]. A fractal approach has been used by Winterwerp [63] to solve a differential equation that simulates the time-varying representative floc diameter, from which floc density is derived from fractal considerations, and settling velocity obtained from a Stokes-like formula. Winterwerp et al. [64] also used a simplified fractal model to relate settling velocity to a turbulent shear parameter, the instantaneous concentration, and water depth.

The state-of-the-art model for floc structure is to assume a fractal structure. Many studies [65, 66, 67] indicate that the most sensitive parameters in a fractal model controlling the resulting settling velocity are the primary particle properties (primary particle diameter, density and their distributions) and the fractal dimension. For organic-rich particles, evidence suggests that the fractal dimension highly depends on stickiness [68].

Floc breakup in the existing size-class-based PBE formulation is modeled by assuming an invariant floc structure (e.g., fractal dimension $\sim 2$) and key properties, such as floc yield strength, are assumed constant over a wide range of floc size. The importance of modeling floc yield strength as a function of floc size has been demonstrated by Son and Hsu [67] via a fractal concept [69, 70]. Further details are provided in Section 3.1 below.

A fractal model is widely used to parameterize floc structure. Fractals are a mathematical simple approach, where scientists feel computationally ‘comfortable’ and therefore are happy to ‘shoehorn’ all flocs into this framework, and a fractal dimension of around 2 is often used. However, its applicability for heterogeneous sediments remains to be proved. Moreover, for flocs of high organic content, stickiness can be significantly enhanced due to the presence of extracellular polymeric substances (EPS) and transparent exopolymer particles (TEP). Field observations suggest that the fractal dimension for inorganic particles is larger than 2.0 while for organic rich flocs, it can be smaller than 2.0 [68]. There are also empirical formulations suggesting that the fractal dimension depends on floc size [65, 66].

Unlike Verney et al. [61], who use floc diameter, Maggi et al. [71] describe the floc population based on the number of primary particles in the flocs, which appears to make the incorporation of a variable fractal dimension straightforward. Moreover, Maggi et al. [71] adopt a sophisticated collisional efficiency closure that considers the effects of floc size and permeability.

2.1.5 Future floc modeling directions

Real flocs are multi component of different densities; even measuring real fractal dimensions is highly problematic. The emergence theory has been utilized in many disciplines (e.g. [72, 73, 74]) and provides a valid alternative and potentially more realistic approach for representing real multi-component mud floc structures. Both Cranford et al. [75], and Rietkerk and van de Koppel [76] have successfully adapted an emergence approach for application to natural biomaterials and ecosystems (respectively).

By utilizing an emergence framework for flocculation, at one end a simple fractal representation would still operate for basic, geometrically repeating simple floc structures (e.g. flocs composed from a single clay primary particle). Whilst as the flocs before more complex in structure, composition and geometry, and fractal theory become less representative, the emergence would adapt to a more suitable floc representation. It is envisaged that this new emergence approach [77] could cope better and more efficiently and realistically for real flocs at a wide range in resolution scales all using real image data at each scale. This will provide a level of error checks that are not supported by regular fractal approaches. Nonetheless, we are some way off from implementing this approach in a numerical floc model and more fundamental research on floc dynamics, properties and characteristics is required, in particular 2D and 3D floc imaging techniques (e.g. [27, 78, 79, 80, 81, 82]).

2.2 Grain-resolving simulations

Grain-resolving simulations are a powerful tool to obtain detailed, high-fidelity data of complex fluid-particle systems. Despite its rather large computational costs, recent advances in computational power have made it possible to perform grain-resolving simulations on scales that become relevant for sediment transport phenomena. The idea is to compute the trajectories of all the individual grains in a flow. Typically, the flow is computed on a Eulerian grid that is fixed in space and time. It can either be approximated by assuming a prescribed background flow or by solving the Navier–Stokes equations. The movement of the particles is then computed by their equations of motion in a Lagrangian sense, i.e. the particle is free to move within the entire computational domain. Hence, this representation is commonly referred to as the Euler–Lagrange approach.

Depending on their fluid-particle coupling, several schemes can be employed. If the particles are driven by a fluid flow but do not modify the flow field, the scheme is considered one-way coupled, whereas the fluid-particle mixture is two-way coupled if the flow is modified by the particle motion as well. In addition, momentum exchange of particles may be accounted for by means of collision and contact or any other particle related force. Using a two-way coupled simulation that accounts for particle-particle interactions is considered a fully-coupled scheme. In the context of the present study, attractive cohesive forces can contribute to the particle-particle interaction by binding grains into aggregates that are much larger than the individual primary particle. This process can be understood as flocculation. In the following, we will first review the governing equations to compute the particle motion in a Lagrangian sense and then proceed to present two examples to model flocculation of cohesive sediment to investigate aggregation and settling processes.

2.2.1 Computing the particle motion

Regardless of the background flow, we prescribe the motion of each primary cohesive particle $i$ as a sphere moving with translational velocity ${\mathbf{u}}_{p,i}$ and angular velocity ${\omega }_{p,i}$. These are obtained from the Newton-Euler equations

$$m_{p,i}\frac{\mathrm{d}{\mathit{u}}_{p,i}}{\mathrm{d}t}=\underset{{\mathit{F}}_{h,i}}{\underbrace{{\oint }_{{\mathrm{\Gamma }}_{p,i}}\tau \cdot \mathit{n}\mathrm{d}A}}+\underset{{\mathit{F}}_{g,i}}{\underbrace{\frac{\pi D_{p,i}^3\left({\rho }_{p,i}-{\rho }_f\right)}{6}\mathit{g}}}+\underset{{\mathit{F}}_{c,i}}{\underbrace{\sum _{j=1,j\ne i}^{N_p}\left({\mathit{F}}_{\mathit{con},\mathit{ij}}+{\mathit{F}}_{\mathit{lub},\mathit{ij}}+{\mathit{F}}_{\mathit{coh},\mathit{ij}}\right)}},$$E5

$$I_{p,i}\frac{\mathrm{d}{\omega }_{\mathbf{p}\mathbf{,}\mathbf{i}}}{\mathrm{d}t}=\underset{{\mathit{T}}_{h,i}}{\underbrace{{\oint }_{{\mathrm{\Gamma }}_{p,i}}\mathit{r}\times \left(\tau \cdot \mathit{n}\right)\mathrm{d}A}}+\underset{{\mathbf{T}}_{c,i}}{\underbrace{\sum _{j=1,j\ne i}^{N_p}\left({\mathit{T}}_{\mathit{con},\mathit{ij}}+{\mathit{T}}_{\mathit{lub},\mathit{ij}}\right)}},$$E6

where the primary particle $i$ moves in response to the hydrodynamic force ${\mathit{F}}_{\mathit{h}\mathbf{,}\mathit{i}}$, the gravitational force ${\mathit{F}}_{\mathit{g}\mathbf{,}\mathit{i}}$ and the particle-particle interaction force ${\mathit{F}}_{\mathit{c}\mathbf{,}\mathit{i}}$ which accounts for the direct contact force ${\mathit{F}}_{\mathit{con}\mathbf{,}\mathit{ij}}$ in both the normal and tangential direction, as well as for short-range normal and tangential forces due to lubrication ${\mathit{F}}_{\mathit{lub},\mathit{ij}}$ and cohesion ${\mathit{F}}_{\mathit{coh},\mathit{ij}}$, where the subscript $\mathit{ij}$ indicates the interaction between particles $i$ and $j$. The hydrodynamic torque is denoted by ${\mathit{T}}_{h,i}$, while ${\mathit{T}}_{c,i}$ represents the torque due to particle-particle interactions, where we distinguish between direct contact torque ${\mathit{T}}_{\mathit{con},\mathit{ij}}$ and lubrication torque ${\mathit{T}}_{\mathit{lub},\mathit{ij}}$. Here, $m_{p,i}$ denotes the particle mass, $D_{p,i}$ the particle diameter, ${\rho }_{p,i}$ the particle density, ${\rho }_f$ the fluid density, $\mathit{g}$ the gravitational acceleration, ${\mathrm{\Gamma }}_{p,i}$ the fluid–particle interface, $\tau$ the hydrodynamic stress tensor, $N_p$ the total number of particles in the flow and $I_{p,i}={\mathit{\pi \rho }}_{p,i}{D}_{p,i}^5/60$ the moment of inertia of the particle. Furthermore, the vector $\mathit{n}$ represents the outward-pointing normal on the interface ${\mathrm{\Gamma }}_{p,i}$, $\mathit{r}$ is the position vector of the surface point with respect to the centre of mass of the particle.

Following Biegert et al. [83, 84], Zhao et al. [85, 86] represent the direct contact force ${\mathit{F}}_{\mathit{con},\mathit{ij}}$ by means of spring-dashpot functions, while the lubrication force ${\mathit{F}}_{\mathit{lub},\mathit{ij}}$ is accounted for based on lubrication theory [87] as implemented in Zhao et al. [85, 86]. The model for the cohesive force ${\mathit{F}}_{\mathit{coh},\mathit{ij}}$ is based on the work of [88], which assumes a parabolic force profile, distributed over a thin shell of thickness $h_{\mathit{co}}$ surrounding each particle:

$${\tilde{\mathit{F}}}_{\mathit{coh},\mathit{ij}}=\left\{\begin{array}{ll}-4\mathit{Co}\frac{{\zeta }_{n,\mathit{ij}}^2-h_{\mathit{co}}{\zeta }_{n,\mathit{ij}}}{h_{\mathit{co}}^2}\mathit{n},& {\zeta }_{\mathit{min}}<{\zeta }_{n,\mathit{ij}}⩽h_{\mathit{co}},\\ 0,& \text{otherwise}.\end{array}\right.$$E7

Here, ${\zeta }_n$ is the gap size in between particle surfaces, ${\zeta }_{\mathit{min}}$ is a limiter and $\mathit{Co}$ is the cohesive number that has to be defined according to a given problem as will be detailed below. We remark that, based on Eqs. (5) and (6), the configuration of the primary particles within a floc can change with time in response to fluid forces, since the cohesive bonds are not rigid.

2.2.2 Aggregation

The flocculation process is strongly affected by the turbulent nature of the underlying fluid flow. Small-scale eddies modify the collision dynamics of the primary particles and hence the growth rate of the flocs, while turbulent stresses can result in the deformation and breakup of larger cohesive flocs. Hence, the dynamic equilibrium between floc aggregation and breakage is governed by a complex and delicate balance of hydrodynamic and inter-particle forces.

In the spirit of earlier investigations [89, 90], Zhao et al. [85] apply a simple model flow in order to investigate the effects of turbulence on the dynamics of cohesive particles. These authors consider the one-way coupled motion of small spherical particles in the two-dimensional, steady, spatially periodic cellular flow field commonly employed as initial condition for simulating Taylor-Green vortices (cf. Figure 6a), with fluid velocity field ${\mathit{u}}_f={\left(u_f,v_f\right)}^{\mathrm{T}}$

Figure 6.

(a) Streamlines of the doubly periodic background flow given by Eq. (8); (b) typical floc configuration with individual flocs distinguished by color (figure taken from [85].

$$u_f=\frac{U_0}{\pi }sin\left(\frac{\mathit{\pi x}}{L_0}\right)cos\left(\frac{\mathit{\pi y}}{L_0}\right),v_f=-\frac{U_0}{\pi }cos\left(\frac{\mathit{\pi x}}{L_0}\right)sin\left(\frac{\mathit{\pi y}}{L_0}\right),$$E8

where $L_0$ and $U_0$ represent the characteristic length and velocity scales of the vortex flow. By computing the fluid flow via this idealized flow field, the hydrodynamic torque ${\mathit{T}}_{h,i}$ can be omitted while the hydrodynamic force ${\mathit{F}}_{h,i}$ is generally replaced by a simple stokes drag force ${\mathit{F}}_{d,i}=-3\pi D_{p,i}\mu \left({\mathit{u}}_{p,i}-{\mathit{u}}_{f,i}\right)$ when simulations are one-way coupled in the sense that the particles do not modify the fluid flow. Here, ${\mathit{u}}_{f,i}$ indicates the fluid velocity evaluated at the particle center, $\mu$ the dynamic viscosity of the fluid.

The dynamics of the primary particles are characterized by the Stokes number $\mathit{St}=U_0{\rho }_{p,i}{D}_{p,i}^2/\left(18L_0\mu \right)$, and the settling velocity $W_s=v_s/U_0$ where $v_s=\left({\rho }_{p,i}-{\rho }_f\right)D_{p,i}^2\mathit{g}/\left(18\mu \right)$ is the Stokes settling velocity of an individual primary particle $i$, as well as the cohesive number

$$\mathit{Co}=\frac{max\left(\Vert {\mathit{F}}_{\mathit{coh},\mathit{ij}}\Vert \right)}{U_0^2{L}_0^2{\rho }_f}=\frac{A_H\left(D_{p,i}+D_{p,j}\right)}{32{\mathit{\lambda \zeta }}_0}\frac{1}{U_0^2{L}_0^2{\rho }_f},$$E9

where the Hamaker constant $A_H$ is a function of the particle and fluid properties, $\lambda =\left(D_{p,i}+D_{p,j}\right)/40$ represents the range of the cohesive force and ${\zeta }_0=\left(D_{p,i}+D_{p,j}\right)/8000$ the characteristic distance. Representative values of $A_H$ for common natural systems can be found in [88].

Zhao et al. [85] employ a computational domain with periodic boundaries. All particles have identical diameters $D_p$ and densities ${\rho }_p$. Initially they are at rest and separated, and randomly distributed throughout the domain. When the distance between two particles is less than $\lambda /2$, the particles are considered as part of the same floc and the number of flocs $N_f$ is tracked as a function of time, with an individual particle representing the smallest possible floc (Figure 6b). Based on the simulation results, Zhao et al. [85] propose a new flocculation model to predict the temporal evolution of the floc size. For flocs of fractal dimension $n_f$, the mean floc size ${\overline{D}}_f$ is related to the average number of primary particles per floc ${\overline{N}}_{p,\mathit{\text{local}}}=N_p/N_f$,

$${\overline{D}}_f={\left({\overline{N}}_{p,\mathit{\text{local}}}\right)}^{\frac{1}{n_f}}{D}_p,$$E10

$${\overline{N}}_{p,\mathit{\text{local}}}=\frac{1}{\left(1/{\overline{N}}_{p,\mathit{\text{local}},\mathit{int}}-1/{\overline{N}}_{p,\mathit{\text{local}},\mathit{max}}\right)e^{\mathit{bt}}+1/{\overline{N}}_{p,\mathit{\text{local}},\mathit{max}}},$$E11

where ${\overline{N}}_{p,\mathit{\text{local}},\mathit{int}}$ denotes the initial number of particles per floc and the average number of particles per floc during the equilibrium stage ${\overline{N}}_{p,\mathit{\text{local}},\mathit{max}}$ is defined as,

$${\overline{N}}_{p,\mathit{\text{local}},\mathit{max}}=\left\{\begin{array}{ll}{N}_p,& {\overline{N}}_{p,\mathit{\text{local}},\mathit{max}}⩾N_p,\\ 8.5a_1{\mathit{St}}^{0.65}{\mathit{Co}}^{0.58}{D}_p^{-2.9}{\varphi }^{0.39}{\rho }_s^{-0.49}{\left(W_s+1\right)}^{-0.38},& \text{otherwise},\end{array}\right.$$E12

where ${\rho }_s={\rho }_p/{\rho }_f$ denotes the density ratio, $\varphi$ the volume fraction of particles. The agglomeration rate $\left|b\right|$ with the constraint $b\le 0$ is obtained by

$$b=\left\{\begin{array}{ll}-0.7a_2{\mathit{St}}^{0.36}{\mathit{Co}}^{-0.017}{D}_p^{-0.36}{\varphi }^{0.75}{\rho }_s^{-0.11}{\left(W_s+1\right)}^{-1.4},& \mathit{St}⩽0.7,\\ -0.3a_2{\mathit{St}}^{-0.38}{\mathit{Co}}^{0.0022}{D}_p^{-0.61}{\varphi }^{0.67}{\rho }_s^{0.033}{\left(W_s+1\right)}^{-0.46},& \mathit{St}>0.7.\end{array}\right.$$E13

For the present cellular model flow the values $a_1=a_2=1$ in Eqs. (12) and (13) yield optimal agreement with the simulation data with the fitting deviation of $\pm 30\%$. For real turbulent flows, $a_1$ and $a_2$ need to be determined by calibrating with experimental data. The new model, Eqs. (10)–(13), with a constant fractal dimension $n_f=2$ for predicting the floc size has been employed and successfully validated with experimental data in our earlier work [85] (cf. Figure 7). The recent study by Zhao et al. [86] even goes beyond this assumption by showing that the fractal dimension becomes a function of the floc size when particles undergo flocculation in isotropic turbulence.

Figure 7.

Calibration of the empirical coefficients for the models of Winterwerp [62] ($k_A^{\prime }=1.35$ and $k_B^{\prime }=1.29\times {10}^{-5}$), Kuprenas et al. [91] ($k_A^{\prime }=0.45$ and $k_B^{\prime }=1.16\times {10}^{-6}$), and for our Eqs. (10)–(13) ($a_1=500$ and $a_2=35$); comparison between experimental data and predictions by the models. The experimental parameters are measured by Tran et al. [92], $D_p=5\mathrm{\mu m}$, ${\rho }_p=2650\mathrm{kg}/{\mathrm{m}}^3$, ${\rho }_f=1000\mathrm{kg}/{\mathrm{m}}^3$, $\mu =0.001\mathrm{Ns}/{\mathrm{m}}^2$, the shear rate $G=50{\mathrm{s}}^{-1}$, concentration $C=200\mathrm{mg}/\mathrm{L}$.

2.2.3 Hindered settling

It is well known that the settling behavior of a dense suspension differs substantially from the settling behavior of a single grain. Particles settling in a dense suspension induce a counterflow and experience friction by colliding with other particles. These processes yield the so-called hindered settling, which is substantially slower than the settling of an individual particle and depends on fluid and particle properties. Nevertheless, the Stokes settling velocity of an individual grain is still widely used to quantify the settling speed of sediment in particle-laden turbidity currents (e.g., [84]). Hence, constitutive equations to predict the settling speed as a function of the local flow conditions can enhance existing computational frameworks for the analysis of turbidity currents. To investigate the effects of the settling behavior of flocculating cohesive sediment by means of grain-resolving simulations [88, 93], it is important to not only account for frictional contact between particles in a dense suspension but also for the modifications of the fluid flow that is caused by settling particles displacing the fluid underneath them [94]. In this case, one needs to solve the Navier-Stokes equations for an incompressible Newtonian fluid:

$$\frac{\partial \mathbf{u}}{\partial t}+\nabla \cdot \left(\mathbf{uu}\right)=-\frac{1}{{\rho }_f}\nabla p+{\nu }_f{\nabla }^2\mathbf{u}+{\mathbf{f}}_{\mathit{IBM}},$$E14

along with the continuity equation

$$\nabla \cdot \mathbf{u}=0,$$E15

where $\mathbf{u}={\left(u,v,w\right)}^T$ designates the fluid velocity vector in Cartesian components, $p$ denotes the pressure, ${\nu }_f$ is the kinematic viscosity, $t$ the time, and ${\mathbf{f}}_{\mathit{IBM}}$ represents an artificial volume force introduced by the Immersed Boundary Method (IBM) [95, 96]. This volume force is introduced in the vicinity of the inter-phase boundaries to enforce a no-slip condition on the particle surface and to modify the fluid motion according to the particle motion. This measure also yields the hydrodynamic forces and torques, ${\mathit{F}}_{h,i}$ and ${\mathit{T}}_{h,i}$ in Eqs. (5) and (6), respectively, as a direct result of this coupling scheme.

While this set of equations represents a fully coupled system, it is important to note that the relevant scales that define the system have changed compared to Section 2.2.2. In the scenario that investigates hindered settling of polydisperse cohesive sediment, the fluid flow is driven by moving particles. Hence, the relevant scale becomes $m_{50}{g}^{\prime }$, where $m_{50}$ is the mass of the median grain size and $g^{\prime }=\left({\rho }_p/{\rho }_f-1\right)g$ is the specific gravity of the sediment. This scaling yields a modified cohesive number:

$$\mathit{Co}=\frac{max\left(\Vert {\mathit{F}}_{\mathit{coh},\mathit{ij}}\Vert \right)}{m_{50g^{\prime }}}$$E16

and a characteristic particle Reynolds number$\mathit{Re}=D_{50}{u}_s/{\nu }_f$, where $D_{50}$ is the median grain size and $u_s=\sqrt{g^{\prime }{D}_{50}}$ is the buoyancy velocity of the sediment. These two non-dimensional numbers, $\mathit{Co}$ and $\mathit{Re}$, fully define the physical system under consideration.

Particles that are placed in a tank will settle due to gravity thereby displacing the fluid and accumulate at the bottom of this tank. The induced counter flow as well as the particle-particle interaction yields frictional contacts and flocculation due to cohesive forces, which is the desired situation for hindered settling [88]. For small particle sizes, where cohesive forces remain relevant, Vowinckel et al. [88, 93] obtain a faster settling behavior for cohesive sediment as compared to its non-cohesive counterpart (Figure 8). During the settling process, the sediment will transform from a suspended state, where the weight is fully supported by the fluid pressure, to a deposited state, where the weight is supported through contact chains of the deposited sediment that extent all the way to the bottom of the tank [93]. This process is described by the effective stress concept, which states that the total stress, i.e. the submerged weight of sediment, is supported by either the particle pressure or the effective stress due to particle contact [97].

Figure 8.

Particle configurations during the settling and deposition process. Top row: Non-cohesive ($\mathit{Co}=0$)and bottom row: Cohesive ($\mathit{Co}=5$ according to Eq. (16)). Left column: $t=17.6{\tau }_s$, which corresponds to the time at which the particle phase has its maximum kinetic energy. Here, ${\tau }_s=D/u_s$ is the non-dimensional reference time. From left to right, the columns are separated by time intervals of $72.5{\tau }_s$. The gray shading reflects the vertical particle velocity (figure taken from [88]).

3. Coastal modeling

Coastal modeling typically refers to the modeling of regional scale (10–100 km²) coastal and estuarine processes over a timescale of hours to months. Due to the large spatial and temporal scales that need to be covered, a Reynolds-averaged Navier-Stokes (RANS) model is adopted (e.g., [98, 99, 100]). Turbulence dissipation and mixing are parameterized with two-equation closure models via a diffusion process. Moreover, when surface waves are present, the individual wave-phase is often not resolved. The generation, transformation and dissipation of the random wave field are represented by a wave spectrum and solved by using the spectral wave action balance equation [101]. Wave-period-averaged wave statistics are then coupled with the coastal models. Consequently, the wave bottom boundary layer processes cannot be directly resolved and additional parameterizations are needed, such as the apparent roughness [102], i.e. the effect of the wave bottom boundary layer on the current resolved by the coastal model, wave-driven small-scale seabed morphological features (e.g., ripples), and the near-bed sediment transport processes (often called bedload or near-bed load). The suspended load transport for a range of non-cohesive sediment classes above the wave bottom boundary layer can be resolved in the coastal models via the conservation of mass. One of the main challenges for extending these coastal models for simulating cohesive sediment transport is the parameterization of settling velocity due to flocculation.

Since the recognition of flocculation in controlling the settling velocity of cohesive sediment in coastal and marine environments [12, 41], significant progress has been made, particularly in the physical parameters controlling the floc dynamics, floc size distribution and their relationships with settling velocity statistics. To name a few, the role of turbulent shear and particle concentration in determining the aggregation rate and the resulting floc size has been quantified [62, 103, 104, 105] (see also Figure 1a). Particularly, in many tidal boundary layers and laboratory experiments of homogeneous turbulence, the mean floc size is observed to be limited by the Kolmogorov length scale (e.g., [106, 107]). Moreover, the relationship between floc size population and settling velocity (or floc density) in both idealized and realistic conditions has been revealed (e.g., [108, 109]) and the fractal dimension has been applied to model these relationships [65, 66, 69, 71] (see also Figure 3). Finally, researchers have begun to understand complex floc characteristics in estuaries dominated by organic particles (e.g., [110]), high cohesion due to TEP (e.g., [111]) and the presence of sand (e.g., [112, 113], see also Figure 5). A more complete discussion on these observational-based empirical parameterizations is provided in Section 2.1.4.

Advancements have also been made in modeling flocculation processes of cohesive sediments. Following the summary presented in Section 2.1, this section focuses on a more in-depth discussion of the complex population approach. Pioneering work by Winterwerp [62, 94] established a robust single-size class (averaged floc size) flocculation modeling framework. This framework has been refined by Kuprenas et al. [91] to limit the floc size growth by the Kolmogorov length scale. A more sophisticated flocculation model [114] based on the population balance equations (PBE) has been incorporated into the Princeton Ocean Models (POM) by Xu et al. [115] to study the dynamics of Estuarine Turbidity Maximum (ETM). Recently, Sherwood et al. [116] incorporated the PBE flocculation model FLOCMOD by Verney et al. [2] into the Regional Ocean Modeling system (ROMS), which is part of the Coupled Ocean-Atmosphere-Wave-Sediment Transport Modeling System (COAWST, [98]). The model is used to study cohesive sediment transport in an idealized setting and a realistic application in the York River estuary. Through a direct simplification of the PBE type model, a tri-modal flocculation model was recently developed in the coastal model TELEMAC [117]. Last but not least, the empirical parameterization of the floc settling velocity MFSV algorithm (see Section 2.1.4) suggested by Manning and Dyer [54] has been incorporated into TELEMAC-3D [57], while more recently the nondimensional version of MFSV proposed by Soulsby et al. [59] has been incorporated into the Finite Volume Community Ocean Model (FVCOM, [99]). These closure models cover a wide range of flocculation physics incorporated or neglected, and the resulting computational cost also varies from minimal to very significant. It is worth to point out that besides the turbulence-averaged models discussed so far, a PBE formulation for flocculation dynamics has recently been incorporated into a turbulence-resolving large-eddy simulation model by Liu et al. [118] to study floc dynamics in the upper-ocean mixing layer subject to Langmuir turbulence.

3.1 Structure of continuum model

The Reynolds-averaged suspended sediment mass concentration $C\left(x_i,t\right)$ is solved by the conservation of mass in an advection-diffusion equation:

$$\frac{\partial C}{\partial t}+\frac{\partial u_{i}C}{\partial x_i}=\frac{\partial {\mathit{CW}}_s{\delta }_{i2}}{\partial x_i}+\frac{\partial }{\partial x_i}\left[\left(\nu +{\nu }_t\right)\frac{\partial C}{\partial x_i}\right]$$E17

where $u_i$ is the fluid velocity, $W_s$ is the settling velocity, $\nu$ is the fluid viscosity and ${\nu }_t$ is the turbulent (eddy) viscosity. The Kronecker delta ${\delta }_{\mathit{ij}}$ is used here with $j=2$ representing the direction of gravitational acceleration. A challenge in modeling cohesive sediments is that the settling velocity, $W_s$, depends on the flocculation process. To include flocculation, there are generally two approaches [119]. The first approach is called the distribution based approach (e.g., [116, 120]), which directly models a bulk settling velocity as a variable due to flocculation. The simplest formulations are empirical formulas characterizing settling velocity as a function of mean floc size, sediment concentration, and turbulent shear rate [64]. Most notably, the framework provided by Winterwerp [62, 94] solves flocculation using a mean floc size, or number concentration of floc (having a mean floc size) that explicitly includes aggregation and breakup terms. Floc density and settling velocity can then be estimated by assuming a constant fractal dimension $n_f$. Such framework has been extended for a variable fractal dimension [67] and applied to model sediment resuspension in the Ems/Dollard estuary by Son and Hsu [121]. Flocculation process involved complex interaction between floc aggregation and breakup of different sizes, ranging from primary particles, microflocs, and macroflocs (see Section 2.1.2 for their definition). Using a mean value of floc property to describe these complex flocculation process may be too simplistic. For example, field observations suggest that natural flocs sometimes show a bimodal distribution (e.g., [18]). Maerz et al. [120] developed a flocculation model that solves the first moment of floc size distribution. Shen and Maa [119, 122] further use the quadrature of moments to model the evolution of floc size distribution. As discussed by Sherwood et al. [116], since settling velocity calculated from the distribution-based approach are based on the statistics of floc properties, the resulting bulk settling velocity must be a variable in time and space. Many existing coastal models utilize an efficient advection scheme that is designed for a constant and uniform settling velocity [115, 123]. Since these coastal models have been designed and routinely used to model suspended sediment (non-cohesive) transport of multiple size classes, it is more straightforward to extend the multiple-size class model into a PBE formulation to model flocculation processes.

3.2 Population balance equation

In the population balance formation, the sediment mass concentration $C$ is partitioned into $N$ classes and each class is represented by an index $k$, where $k$ and $N$ are positive integers. $N$ must be a sufficiently large number to resolve the distribution. The sediment mass concentration in each class $c_k\left(x_i,t\right)$ is calculated by its mass conservation equation:

$$\frac{\partial c_k}{\partial t}+\frac{\partial u_i{c}_k}{\partial x_i}=\frac{\partial c_k{w}_{s,k}{\delta }_{i2}}{\partial x_i}+\frac{\partial }{\partial x_i}\left[\left(\nu +{\nu }_t\right)\frac{\partial c_k}{\partial x_i}\right]+G_{1,k}-L_{1,k}+\dots ,k=1,2,\dots N$$E18

where $w_{s,k}$ is the settling velocity of class $k$. To include the flocculation processes, additional gain and loss terms, e.g., $G_{1,k},L_{1,k},\dots$, due to aggregation of flocs (and primary particles) and breakup of flocs need to be modeled, which will be discussed later. It is important to point out that the flocculation processes only re-distribute the floc mass among different size classes and care must be taken to ensure the total sediment mass conservation when modeling these terms (i.e., these gain and loss terms must cancel each other after the summation of all size classes in Eq. (18)), including the numerical treatment (e.g., using a logarithmically distributed size class and a mass-weighted interpolation).

As mentioned before, the sediment mass concentration distribution can be described by floc size class [61] or number of primary particles in the floc [66]. Here, we focus on the more popular one using floc size class with each class having a floc diameter of $D_{f,k}$. The sediment mass concentration of each class can be related to number concentration of flocs in each size class $n_k\left(x_i,t\right)$ as:

$$c_k\left(x_i,t\right)=m_k{n}_k\left(x_i,t\right)$$E19

where $m_k$ is the mass of a floc in size class $k$. Assuming a floc is formed by spherical primary particles of diameter $D_p$ and density ${\rho }^s$, the mass of a primary particle can be calculated as

$$m_p=\frac{\pi }{6}{\rho }^s{D}_p^3.$$E20

By, furthermore, assuming that the floc structure follows a fractal relationship, we can calculate the mass of a floc as [69]:

$$m_k=m_p{\left(\frac{D_{f,k}}{D_p}\right)}^{n_f}.$$E21

Therefore, with the given fractal dimension and primary particle properties, $m_k$ can be explicitly calculated for a given floc size class $D_{f,k}$. Following the fractal theory, the floc density of size class $k$ can be readily calculated as [69].

$${\rho }_k^f={\rho }^w+\left({\rho }^s-{\rho }^w\right){\left(\frac{D_{f,k}}{D_p}\right)}^{n_f}$$E22

where ${\rho }^w$ is the density of water (or seawater). With the known floc diameter and floc density, the settling velocity of the flocs in each size class can be calculated. The simple Stokes law is used here following Sherwood et al. [116]

$$w_{s,k}=\frac{\left({\rho }_k^f-{\rho }^w\right){\mathit{gD}}_{f,k}^2}{18\mu }$$E23

where $\mu$ is the dynamic viscosity of water (seawater).

It is worth noting that in the population balance formulation, the floc settling velocity of a particular size class is treated as a constant determined by the given fractal dimension, primary particle properties, and fluid properties. Therefore, it is more suitable for typical coastal modeling systems due to their numerical treatment of advection. In typical field conditions, a size-class based population balance formulation is reported to require at least 10–20 size classes [61, 116].

After proposing the appropriate gain and loss terms in Eq. (18), the full dynamics of floc transport, settling and re-distribution of sediment mass among all floc size classes due to flocculation can be modeled. In practice, some models solve a system of partial differential equation for the number concentration $n_k\left(x_i,t\right)$ directly by substituting Eq. (19) into Eq. (18) (e.g., Liu et al. [118]). For typical coastal models (e.g., Sherwood et al. [116]), the numerical treatment of flocculation and the advection-diffusion-settling processes are split into two steps. The zero-dimensional (homogenous turbulence condition without the advection, diffusion and settling terms) number concentration equations that only include the gain and loss terms are solved at every grid point over each (baroclinic) time step size of the coastal model to redistribute floc mass between different classes. Then, the newly calculated floc mass due to flocculation in each size class is updated by the advection-diffusion-settling equation (Eq. (18)) without the gain and loss (flocculation) terms. One advantage of such approach is that many zero-dimensional flocculation models developed elsewhere can be easily coupled into the coastal models. For example, Sherwood et al. [116] couple an existing zero-dimensional size-class based flocculation model FLOCMOD developed by Verney et al. [61] into the COAWST coastal modeling framework.

In this paper, the essential model formulations and closures of Verney et al. [61] are reviewed. The purpose here is not to discuss the model details, since they can be found in the cited references. Rather, this section is intended to bridge the discussions in Section 1 and Section 2 by focusing on the key model elements that are sensitive to the model results and hence may require more physical understanding. The governing equation for floc number concentration $n_k$ (S.I. unit: ${\mathrm{m}}^{-3}$) of size class $k$ is written as

$$\frac{\partial n_k}{\partial t}=G_a\left(k\right)-L_a\left(k\right)+G_{\mathit{bs}}\left(k\right)-L_{\mathit{bs}}\left(k\right)$$E24

in which $G_a\left(k\right)$ and $L_a\left(k\right)$ represent the gain and loss of flocs in size class $k$ due to aggregation, while $G_{\mathit{bs}}\left(k\right)$ and $L_{\mathit{bs}}\left(k\right)$ represent the corresponding gain and loss due to shear-driven breakup. Verney et al. [61] also include terms due to collision-induced breakup [70]. However, this mechanism is shown to be of minor importance for pure clay flocculation and it is not discussed here for brevity. The model assumes aggregation that is driven by turbulent shear and binary collision. Other known mechanisms driving aggregation, namely, the Brownian motion (important for particle smaller than $1\mathrm{\mu m}$) and the differential settling are neglected by assuming that coastal and estuarine environments are dominanted by turbulent shear. The $G_a\left(k\right)$ and $L_a\left(k\right)$ terms are modeled as

$$G_a\left(k\right)=\frac{1}{2}\sum _{i+j=k}{\alpha }_{\mathit{ij}}A\left(i,j\right)n_i{n}_j$$E25

and

$$L_a\left(k\right)=\sum _i^N{\alpha }_{\mathit{ik}}A\left(i,k\right)n_i{n}_k$$E26

where ${\alpha }_{\mathit{ij}}$ represents the collisional efficiency (dimensionless) and the shear-driven binary collision probability function is written as

$$A\left(i,j\right)=\frac{1}{6}G{\left(D_i+D_j\right)}^3.$$E27

The quantity $G_a\left(k\right)$ shown in Eq. (25) represents the gain of flocs in size $k$ due to aggregation of smaller flocs at size class $i$ and $j=k-i$ while $L_a\left(k\right)$ shown in Eq. (26) expresses the loss of flocs in size class $k$ due to aggregation with flocs in other size classes. In Verney et al. [61], ${\alpha }_{\mathit{ij}}$ is set to be a constant ${\alpha }_{\mathit{ij}}=\alpha$ for simplicity. This point will be revisited later. Gain and loss of flocs in size class $k$ due to breakup driven by turbulent shear are calculated as

$$G_{\mathit{bs}}\left(k\right)=\sum _{i=k+1}^N{\mathrm{\Pi }}_{\mathit{ki}}{B}_i{n}_i$$E28

and

$$L_{\mathit{bs}}\left(k\right)=B_k{n}_k.$$E29

Essentially, gain of flocs $G_{\mathit{bs}}\left(k\right)$ in size class $k$ due to fragmentation can only occur in floc size classes larger than $k$. The fragmentation probability function is written as

$$B_i={\beta }_i{G}^{3/2}{D}_i{\left(\frac{D_i-D_p}{D_p}\right)}^{3-n_f}$$E30

where ${\beta }_i$ is the (dimensional: ${\mathrm{m}}^{-1}{\mathrm{s}}^{1/2}$) fragmentation rate and it is assumed to be a constant ${\beta }_i=\beta$ in Verney et al. [61]. More discussions on ${\beta }_i$ will be given later. The breakup distribution function is defined as ${\mathrm{\Pi }}_{\mathit{ki}}$ and it characterizes different breakup scenarios. Our understanding on the shear-induced breakup scenarios is currently limited. Nevertheless, the breakup distribution functions ${\mathrm{\Pi }}_{\mathit{ki}}$ provided in Verney et al. [61] include binary (fragmentation of floc mass $m_k$ into two small flocs of equal mass $m_k/2$), ternary (fragmentation of floc mass $m_k$ into three smaller flocs, one floc of $m_k/2$ and two flocs of $m_k/4$), and erosion (fragmentation of floc mass $m_k$ into $r+1$ flocs, one larger floc and $r$ smaller flocs of equal size).

As demonstrated in Verney et al. [61], their PBE-based flocculation model can predict several key features of floc dynamics observed in the field. For instance, the model is capable of reproducing the observed slower aggregation and more rapid fragmentation process (so-called clock-wise hysteresis of aggregation/fragmentation process) during a tidal cycle. As the floc size directly controls the settling velocity, capturing the hysteresis of floc aggregation/fragmentation is essential to further predict the net sediment fluxes during a tidal cycle. Moreover, Verney et al. [61] showed that a bimodal distribution of flocs often observed in the field can be reproduced by the PBE model by including a mix of different breakup distribution functions. Although the PBE-based flocculation models provide a promising modeling framework for cohesive sediment transport, there are limitations that require future investigations.

The most sensitive empirical parameters in FLOCMOD are the collisional efficiency $\alpha$ and the fragmentation rate $\beta$, which are often assumed to be constants. As demonstrated systematically by Verney et al. [61], to match the mean equilibrium floc size, there is no unique set of optimum $\alpha$ and $\beta$ values and it is the ratio of $\alpha /\beta$ for these two empirical parameters that controls the mean equilibrium floc size. A similar finding is reported in the single size-class flocculation model of Winterwerp [62]. Hence, to fully benefit from the capability of PBE-based models that provide the temporal evolution of the full spectrum of floc sizes, the model calibration should go beyond just using equilibrium mean floc size. Sherwood et al. [116] provide a limited validation of Verney et al. [61] model for the temporal evolution of mean floc size. Such validation should be expanded for different data sets and for different floc size statistics in order to better constrain the empirical model parameters.

The sensitivity of the modeled mean floc diameter to the prescribed fractal dimension has also been discussed in Verney et al. [61]. More recently, Penaloza-Giraldo et al. [124] report that the temporal evolution of floc sizes (timescale to reach the equilibrium floc sizes) are sensitive to the fractal dimension $n_f$. As shown in Eq. (30), when $n_f$ is smaller, the fragmentation probability function becomes larger, and hence it takes a shorter time to reach the equilibrium floc size. Even when $\alpha /\beta$ and $n_f$ are calibrated to match the measured mean equilibrium floc size and flocculation timescale, it may not be straightforward to match the model results with measured data for the entire floc size distribution. Specifically, since the choice of breakup distribution function ${\mathrm{\Pi }}_{\mathit{ki}}$ is not constrained, the calibrated $\alpha /\beta$ and $n_f$ may also depend on the ${\mathrm{\Pi }}_{\mathit{ki}}$-function. In summary, the sensitivity study of Verney et al. [61] and preliminary findings reported by Penaloza-Giraldo et al. [124] imply that the following flocculation physics controlling the empirical parameters in the PBE warrant future studies:

3.2.1 Parameterizing the floc structure with a fractal dimension

To estimate floc mass and floc density in a given floc size, further assumptions on floc structures is needed. The state-of-the-art model for floc structure is to assume a fractal structure, which renders Eqs. (21) and (22) useful. In the PBE models for floc dynamics, fractal dimension directly affects the breakup term via the fragmentation probability function (see Eq. (30)). While most of the existing flocculation models assume a constant fractal dimension, examining the field and laboratory data by relating measured floc settling velocity and floc size (see Eqs. (22) and (23)) suggests that the fractal dimension may not be a constant, especially when considering the PBE equations are very sensitive to the prescribed fractal dimension value. Recent grain-resolving simulations (see Section 2.2.2 and the studies of [85, 86] referenced therein) also confirms that the fractal dimension depends on floc size. Researchers have proposed to model fractal dimension as a function of floc size [65, 71], however, whether it significantly improves the modeled flocculation processes remain to be proven.

3.2.2 More complete descriptions of collisional efficiency and fragmentation rate

As discussed by Hill and Nowell [125], the collision efficiency is practically treated as an empirical tuning parameter that parameterizes three main processes: encounter, contact and sticking. Encounter and contact are physical processes while sticking is associated with chemical-biological processes. From a physical perspective, only sticking efficiency is solely an empirical parameter. Maggi et al. [71] used a more complex collisional efficiency formulation that depends on the size and porosity of two colliding flocs. Through detailed laboratory experiments, Soos et al. [126] proposed a collisional efficiency formulation ${\alpha }_{\mathit{ij}}$ that depends on the size of the two colliding flocs and the turbulent shear rate. More systematic studies on the impact of collisional efficiency on the modeled flocculation in the PBE-based formulation are required.

While the existing studies mostly treat the fragmentation rate to be a model constant, physically the fragmentation rate ${\beta }_i$ further depends on the floc breakage force $F_y$ and the dynamic viscosity of the fluid [71]:

$${\beta }_i=E{\left(\frac{\mu }{F_y}\right)}^{1/2}$$E31

with $E$ an empirical constant. It is clear that the only way to justify that the fragmentation rate is a constant, i.e., ${\beta }_i=\beta$, is to assume that the floc breakage force $F_y$ is independent of floc size, which theoretically is consistent with the fractal theory [69]. Indeed, in most cohesive sediment transport literature (e.g., [61, 71, 115]), $F_y$ is assumed to be $F_y=\mathcal{O}\left({10}^{-10}\right)N$ following Winterwerp [62]. Since assuming floc structure follows the fractal theory (independent of the size of the aggregates) is an approximation, $F_y$ may need to be considered as a variable in practice. Jarvis et al. [127] presented a review of many laboratory measurements of floc breakage force per unit area (it was called floc strength in their paper). Due to a wide range of cohesive sediment samples tested and different measurement techniques used, the quantitative results are not conclusive. However, it is clear that $F_y/D_f^2$ decreases as floc size increases. This conclusion does not exclude the possibility that $F_y$ is a variable and it may be insensitive, or has a certain degree of dependence on floc size. More quantitative investigations on floc breakup force for different mineral composition and turbulence intensity are needed.

3.2.3 Improved physical understanding of breakup distribution functions

In the coastal sediment literature, detailed studies on floc breakup, particularly from the observational perspective, are rare. In the water quality literature, Jarvis et al. [127] provide some insights into the breakup distribution function. First, they discern the fragmentation mechanism, similar to the binary breakup, as the most likely scenario to occur when flocs are subjected to tensile stress acting across the floc. Secondly, the erosion mechanism in floc breakup is likely due to shear stress acting tangentially to the floc surface. Based on this argument, researchers hypothesize that the floc breakage types may depend on the ratio of floc size to the Kolmogorov length scale (smallest turbulent eddy size). However, the results are not conclusive and more comprehensive studies on how turbulent eddies interact with flocs and causing floc breakage are warranted.

4. Conclusions

We have presented an overview covering different types of floc analyses based on experimental measurements and grain-resolved simulations. These tools are currently emerging and show a very promising perspective to generate the data needed to account for unresolved cohesive sediment dynamics in continuum models with high fidelity. More work will be needed in the future to cover the different aspects laid out in this chapter. Those are in particular, the effects of biofilms, the settling velocity of different types of flocs, as well as the aggregation and break-up efficiencies governing the exchange between different classes of PBE-type models. The knowledge to be gained can lead to a new generation of continuum models that enable simulations with predictive power for entire estuaries, which will bring inestimable advantages for these attractive settlement areas, both in economic terms and in terms of an increased quality of life.

Acknowledgments

BV gratefully acknowledges the support through the German Research Foundation (DFG) grant VO2413/2-1. KZ is supported by the National Natural Science Foundation of China through the Basic Science Center Program for Ordered Energy Conversion (51888103). LY, AJM, TJH and EM gratefully acknowledge supports by US National Science Foundation through a collaborative research between the University of Delaware (OCE-1924532) and the University of California Santa Barbara (OCE-1924655). AJM’s contribution towards this research was partly supported by the US National Science Foundation under grants OCE-1736668 and OCE-1924532, TKI-MUSA project 11204950-000-ZKS-0002, and HR Wallingford company research FineScale project (ACK3013_62).

Conflict of interest

The authors declare no conflict of interest.

Nomenclature

αij=α	collisional efficiency
βi=β	fragmentation rate
Γp,i	fluid-particle interface of particle i
δi,j	Kronecker delta
ζ0	characteristic distance
ζn	gap size in between particles
ζmin	limiting gap size in between particles
λ	cohesive force range
μ	dynamic fluid viscosity
ν	kinematic fluid viscosity
νt	turbulent viscosity
Πki	breakup distribution function
ρf	fluid density
ρkf	floc density of class k
ρp,i	density of particle i
ρs	density ratio
ρs	primary particle density
ρw	density of water
τ	turbulent shear stress
τ	hydrodynamic stress tensor
τs	nondimensional reference time
ϕ	particle volume fraction
ωp,i	Angular velocity vector of particle with index i
Aij	Shear-driven binary collision probability function
Ah	Hamaker constant
a1	model constant
a2	model constant
b	agglomeration rate
Bi	fragmentation probability function
C	SPM concentration
c	empirical constant
ck	sediment mass concentration of class k
Co	cohesive number
D	floc size
D50	median grain size
D¯f	mean floc size
Df,k	floc diameter of class k
Dp	primary particle diameter
Dp,i	diameter of particle i
E	empirical constant
Fc,i	particle-interaction force vector of particle with index i
Fcoh,ij	cohesive force vector of particle i interacting with particle j
Fcon,ij	normal contact force vector of particle i interacting with particle j
Fh,i	hydrodynamic force vector of particle with index i
Fg,i	gravitational force vector of particle with index i
fIBM	artificial volume forced introduced by the IBM i
Flub,ij	lubrication force vector of particle i interacting with particle j
Fy	floc breakage force
G	turbulence parameter / shear rate
G1,k	mass gain due to aggregation of flocs
Ga	gain of flocs in a size class due to aggregation
Lbs	gain of flocs in a size class due to shear driven break-up
g	gravitational acceleration
g′	specific gravity
hco	thickness of cohesive shell
Ip,i	moment of inertia of particle i
k	empirical constant
k	index of class for sediment mass concentration
kA′	agglomeration constant
kB′	breakup constant
L0	characteristic length
L1,k	mass loss due to breakup of flocs
La	loss of flocs in a size class due to aggregation
Lbs	loss of flocs in a size class due to shear driven break-up
m	empirical constant
mk	mass of an individual floc in class k
N	number of classes for sediment mass concentration
nf	fractal dimension
Nf	total number of flocs in the flow
nk	number concentration of class k
Np	total number of particles in the flow
N¯p,local	average number of primary particles per floc
m50	mass of the median grain size
mp,i	mass of particle i
mk	mass of a floc
r	radial position vector
Re	Reynolds number
St	Stokes number
t	time
Tc,i	Torque vector of particle with index i due to particle-interaction
Tcon,ij	collision torque vector of particle i interacting with particle j
Th,i	hydrodynamic torque vector of particle with index i
Tlub,ij	lubrication torque vector of particle i interacting with particle j
U0	characteristic velocity
uf,i=u	fluid velocity vector
up,i	translational velocity vector of particle with index i
us	buoyancy velocity
vs	Stokes settling velocity
Ws	settling velocity of microflocs
ws,k	settling velocity of of class k
Ws,macro	settling velocity of macroflocs

Abbreviations

ADV	acoustic doppler velocimeter
COAWST	Coupled Ocean-Atmosphere-Wave-Sediment Transport Modeling System
DNS	direct numerical simulation
EPS	extracellular polymeric substance
ETM	estuarine turbidity maximum
FVCOM	finite volume community ocean model
IBM	immersed boundary method
INSSEV	IN-Situ SEttling Velocity instrument
LabSFLOC	laboratory spectral flocculation charactersitics
MFSV	manning floc settling velocity
MSF	mass settling flux
OBS	optical backscatter
PBE	population balance equation
PICS	particle imaging camera system
POM	Princeton ocean models
RANS	Reynolds-averaged Navier-Stokes
ROMS	regional ocean modeling systems
SB	size band
SPM	suspended particulate matter
TEP	transparent exopolymer particles