Floc Size Distributions of Cohesive Sediment in Homogeneous Isotropic Turbulence

Floc size distribution is one of the key parameters to characterize ﬂ occulating cohesive sediment. An Eulerian – Lagrangian framework has been implemented to study the ﬂ occulation dynamics of cohesive sediments in homogeneous isotropic turbulent ﬂ ows. Fine cohesive sediment particles are modeled as the dispersed phase by the discrete element method, which tracks the motion of individual particles. An adhesive contact model with rolling friction is applied to simulate the particle – particle interactions. By varying the physicochemical properties (i.e., stickiness and stiffness) of the primary particles, the dependence of the mathematical form of the ﬂ oc size distribution on sediment properties is investigated. At the equilibrium state, the aggregation and breakup processes reach a dynamic equilibrium, in which construction by aggregation is balanced with destruction by breakup, and construction by breakup is balanced with destruction by aggregation. When the primary particles are less sticky, ﬂ oc size distribution ﬁ ts better with the lognormal distribution. When the primary particles are very sticky, both the aggregation of smaller ﬂ ocs and breakup from larger ﬂ ocs play an equally important role in the construction of the intermediate-sized ﬂ ocs, and the equilibrium ﬂ oc size distribution can be better ﬁ tted by the Weibull distribution. When the Weibull distribution develops, a shape parameter around 2.5 has been observed, suggesting a statistically self-similar ﬂ oc size distribution at the equilibrium state. review, and editing. SB: conceptualization; methodology, and writing — review and editing. AM: writing – review and editing.


INTRODUCTION
The transport of fine-grained cohesive sediment in nearshore and estuarine environments plays a critical role in ecosystem dynamics, water quality, bed morphology, and engineering applications, for example, the rapid siltation in navigation channels and harbors (Hayter and Mehta, 1986;Winterwerp et al., 2000), cohesive sediment transport in salt marsh (Graham and Manning, 2007), depositional rates of contaminated muddy sediments (Ye et al., 2020), and long-term morphology of deltas (Edmonds and Slingerland, 2010). Cohesive sediment can bind together through both physical (van Olphen, 1964;Winterwerp and van Kesteren, 2004) and biological (Tolhurst et al., 2002) cohesion to form large aggregates, namely, flocs. A floc size distribution develops in sediment suspension (Sheremet et al., 2017). Due to the variability in floc's structure and effective density, flocs of different size settle at different velocities (Manning, 2004;Mehta et al., 2014). The larger, low-density macroflocs (van Leussen, 1994) tend to settle faster than smaller microflocs (Eisma, 1986;Manning, 2001), but they are more fragile and more prone to break up by turbulent shear. Macroflocs often dominate the depositional mass flux (Mehta and Lott, 1987;Manning et al., 2006). The floc size distribution is therefore of crucial importance in understanding the spatiotemporal transport patterns of cohesive sediment (e.g., Geyer et al. (2000); Baugh and Manning (2007); Prandle et al. (2005)).
Turbulence plays an important role in the flocculation process of cohesive sediment in natural environments (Dyer, 1989;van Leussen, 1997;Winterwerp, 1998;McAnally and Mehta, 2001;Manning, 2004). On the one hand, turbulence enhances the aggregation through the collision frequency, which scales with the turbulent dissipation rate (Sundaram and Collins, 1997). On the other hand, large flocs break in turbulent flows by turbulent eddies via turbulent shear or hydrodynamic drag (Saha et al., 2016). Several phases exist during the flocculation. Initially, the aggregation dominates with the rapid growth of the floc size. As flocs continue to grow, large flocs with porous structures form. Large flocs are vulnerable to fragmentation by fluid shear (Tambo and Watanabe, 1984). Breakup starts to play an increasingly important role in late stages of flocculation. When the two competing mechanisms, namely, the aggregation and breakup processes, balance, an equilibrium floc size distribution develops Soulsby et al., 2013;Mehta et al., 2014).
Due to the large variability in the floc size, cohesive sediment is often characterized by the floc size distribution. The mathematical properties of floc size distributions have drawn a lot of attention from the cohesive sediment transport research community, and the interest in unifying the properties of floc size distribution has remained strong. Various statistics for floc size distribution have been proposed to serve as indices of the quality of sediment flocs, as well as sludge in waste treatment. However, theoretical studies, field observations, and laboratory experiments yield different statistics. It is important to investigate the physical mechanisms that lead to different floc size distributions, and the potential implication of different mathematical forms of the floc size distribution.
By applying a dimensional analysis, Hunt (1982) showed the steady state floc size distribution follows a power law. Pushkin and Aref (2002) later developed a more rigorous self-similarity theory of stationary coagulation and showed the floc size distribution follows a power law in the coagulating system. In these studies, the system is forced with particle injection at small sizes, and breakup is not considered. The breakup of large flocs can lead to a skewed floc size distribution with a peak (Hunt, 1982). Spicer and Pratsinis (1996) conducted laboratory experiments to study the evolution of floc size distribution induced by shear and showed the steady state floc size distribution normalized by the average floc size to be selfpreserving, which is independent of the shear rate.
Floc size distribution is skewed and hence does not tend to follow the normal distribution. The lognormal distribution and Weibull distribution are commonly used to model skewed distributions; however, the physical origin of the distribution is not well understood. Brown and Wohletz (1995) derived the Weibull distribution with respect to the fragmentation process, in which a power law was used to describe the breakup of a single particle into smaller particles. The Weibull distribution has been widely used as particle size distribution for coarse grains (Fang et al., 1993;Kondolf and Adhikari, 2000). Previous studies of fiber pulp suspension in a flat channel (Huber et al., 2006) and activated sludge flocs (Li and Ganczarczyk, 1991) showed that Weibull distribution is the best descriptor for the floc size distribution. On the other hand, Kiss et al. (1999) developed a model for particle growth that predicts the lognormal particle distribution. They assumed the rate of change of the particle mass is proportional to the surface area, and the particle residence time in the active zone of particle interactions is lognormally distributed. Floc growth is due to collisions with other flocs, and the collision frequency is proportional to the surface area of the floc. A lognormal distribution of velocity fluctuations (Mouri et al., 2009) or dissipation rate (Yeung et al., 2006) that drive inter-particle collisions could also lead to the lognormal floc size distribution. Byun and Son (2020) applied a stochastic approach to model the size distribution of suspended flocs, in which the breakup process is modeled by a lognormal distribution. They showed the lognormal distribution is the best descriptor for the floc size distribution. Hosoda et al. (2011) showed that a stochastic process of halving followed by addition can yield a stationary lognormal distribution. For cohesive sediment flocs, this suggests the breakup of a large floc into two small flocs of equal size followed by the aggregation with another floc could lead to a lognormal floc size distribution. Overall, it is difficult to distinguish the lognormal and Weibull distribution in floc size distribution curves and hence the physical origin of the size distribution, which requires priori knowledge on both the particle-particle and particle-fluid interactions during flocculation.
In the mathematical approach, the aggregation and breakup processes are parameterized. The accuracy of the predictive cohesive sediment transport model strongly depends on the aggregation and breakup models. The two-phase Eulerian-Lagrangian model can resolve both the particle-particle and particle-fluid interactions and can provide the particle-level information on the aggregation and breakup processes. In Eulerian-Lagrangian two-phase models, the carrier fluid is modeled as the continuous phase and the particles are modeled as the dispersed phase (Balachandar and Eaton, 2010). In total, two approaches, namely, the particle-resolving approach (PR) and the point-particle approach (PP), have been developed and implemented to study cohesive sediment dynamics. In both approaches, the discrete element method (DEM) models the particle-particle interactions. Particles are modeled as soft spheres, allowing a small overlap when two particles collide. When one particle collides with another particle or floc, they may stick together. In DEM, the motion of an individual particle is tracked, along with the aggregation and breakup of flocs. Collisions among particles are modeled by the contact mechanics theory, such as Hookean or Hertzian contact models (Johnson, 1985). In the particle-resolving approach (Vowinckel et al., 2019), flows around individual particles are fully resolved. Due to the large computational cost, the particle resolving approach is often limited to systems with a few thousands of particles, which may not generate satisfactory statistics for flocculation dynamics. On the other hand, in the point-particle approach (Marshall, 2009;Zhou et al., 2010), hydrodynamic forces, such as drag force, lift force, and inertial force, on the particle are modeled. The point-particle approach can be implemented to millions of particles easily. However, the accuracy of the point-particle approach strongly depends on the hydrodynamic force models (Akiki et al., 2017). To investigate the flocculation processes in homogeneous isotropic turbulence, the point-particle approach is implemented to get better statistics of particle dynamics.
In this study, we investigate the floc size distribution in homogeneous isotropic turbulence using a two-phase Eulerian-Lagrangian model, in which the particle-particle interactions are modeled by the discrete element method. Due to the limit of computational resources, we focused on flocculation processes in high-energy estuaries or near-field river plumes with high turbulent shear rate, in which turbulence dictates the aggregation, breakup, and restructuring processes of flocs. We investigated how the primary particle properties affect the aggregation and breakup processes and hence the floc size distribution by varying the stickiness, stiffness, and size of the primary particle while keeping the turbulent shear rate the same. We focus on the physical origin of the floc size distribution and assess the performance of the lognormal distribution and the Weibull distribution at the equilibrium stage. This study is organized as follows. Methods are described in Section 2, including the adhesive contact model and the one-way coupling of the fluid and particle phases. Model validation and model setup are also presented in Section 2. Model results are presented in section 3 followed by the discussion in Section 4 and concluding remarks in Section 5.

Direct Numerical Simulation of Homogeneous Isotropic Turbulence
Turbulence is characterized by a wide range of length scales. Interactions between turbulent eddies of different length scales with flocs play a critical role in flocculation dynamics. The primary particles are smaller than the Kolmogorov length scale (Kolmogorov, 1941a,b) in this study. Although the flow around the particle is not resolved, all turbulent scales including the Kolmogorov scale and larger, are fully resolved. Thus, the present approach is the particle-unresolved direct numerical simulations (DNS).The homogeneous isotropic turbulence is implemented in this study, which is an idealized version of the realistic turbulence and a reasonable approximation of the turbulent flow away from bottom boundary. To generate homogeneous isotropic turbulence, the linear forcing model (Lundgren, 2003;Rosales and Meneveau, 2005) was implemented. Instead of applying forces only to lowwavenumber modes, a force proportional to velocity is introduced in the momentum equation in the form of αu. Because the volumetric sediment concentration is dilute (≪ 1%) and the dominant effect is that of the turbulent carrier flow on the particle dynamics, the one-way coupling approach is adopted, and the governing equations of the fluid phase are as follows: in which u is the fluid velocity, p is the pressure, ρ is the density of the fluid, ] is the kinematic viscosity of the fluid, and α is the linear forcing coefficient. The direct numerical simulations were conducted with the open source code Nek5000 (Fischer et al., 2008;Zwick and Balachandar, 2020), which uses a high-order spectral element method.

Discrete Element Method for Cohesive Sediment
To resolve particle-particle interactions, the sediment phase is modeled by using the discrete element method (DEM), in which motions of individual particles are tracked.
x is the position vector, v is the particle velocity vector, F is the force vector, and m is the mass of the particle. The subscript "i" is the particle label. The force on particle "i" is the sum of the collision force (F c ) between particle i and all other particles j, the hydrodynamic force (F f ), and the gravitational force (F g ) as F i = j F c,ij + F f,i + F g,i . I is the moment of inertia, ω is the angular velocity of the particle, and T is the torque on the particle. In this study, we coupled the CFD code nek5000 with the molecular dynamic code LAMMPS (Plimpton, 1995). The granular package in LAMMPS provides a variety of options for the normal, tangential, rolling, and twisting forces resulting from the contact between two particles, and hence is used to model the complex interactions among cohesive sediment particles. For soft clay particles, the Johnson-Kendall-Roberts (JKR) model is adopted.
where a is the radius of the contact zone and is related to the overlap δ according to where E is the Young's Modulus, R is the radius of the particle, and γ is the surface energy density. The overlap between particle "i" and particle "j" is given as δ = R i + R j − |x i − x j |. The JKR model allows for a tensile force beyond contact (δ < 0), up to a maximum of 3πγR. When two particles are not in contact initially, they will not experience this force until they come into contact (δ > 0), then as they move apart, they experience a tensile force up to 3πγR till they lose contact. This force can be used to define the yield strength of the floc. In addition, a viscoelastic damping force model is used.
where η D is the viscoelastic damping coefficient, and V n is the relative velocity along the direction of the vector n, which is the unit vector along the line connecting the centers of the two particles. The total normal force is the sum of the adhesive JKR and viscoelastic damping terms.
The Mindlin no-slip model (Mindlin, 1949) is used to compute the tangential force (F t ), as follows: where μ t is the friction coefficient, k t is the elastic constant for tangential contact, and ξ is the tangential displacement accumulated during the entire duration of the contact. The vector t is the unit vector in the relative tangential velocity direction. F t,d is the damping term for the tangential force, which follows the same general form as the normal damping force (Eq. 9) but uses the relative velocity along the direction of the tangential vector t. The normal force value F n0 used to compute the critical force is given as follows: The floc restructuring, in which particles change their relative positions while remaining connected, could also play an important role in the flocculation dynamics. Compaction of flocs by turbulent shear may occur with preferential floc structures. To account for floc restructuring, a rolling friction model of a pseduo-force formulation (Luding, 2008) was implemented. The rolling friction model allows the adjustment of rolling displacement of the contacting pair. The rolling pseudo-force is computed analogously to the tangential force, as follows: where k roll is the elastic constant for rolling, γ roll is the damping constant for rolling, ξ roll is the rolling displacement, and v roll is the relative rolling velocity (Wang et al., 2015). A Coulomb friction criterion truncates the rolling pseudo-force if it exceeds a critical value of where k is the direction of the pseudo-force. The rolling pseudo-force does not contribute to the total force on either particle, but it acts only to induce an equal and opposite torque on each particle.
T roll,j −T roll,i .

Hydrodynamic Force
The total hydrodynamic force on particle "i" is given as follows: where F d and F p are the quasi-steady force and stress-divergence force, respectively. The added-mass force is neglected in this study assuming the small particle Stokes number. The drag force F d on particle "i" is given as follows: where ρ is the fluid density and A πD 2 p /4 is the projected area of the spherical particle with D p as the diameter of the spherical particle. For very dilute flow with sediment concentration ϕ ≪ 0.1%, the standard drag coefficient C D for an individual particle is used, which is given as follows: where Re p = |u − v|D p /] is the particle Reynolds number. The stress-divergence force experienced by the particle is calculated as follows: where the pressure gradient and stress divergence are interpolated to the particle center. In the current formulation, only the hydrostatic pressure is used to calculate the force F p for simplicity. The buoyancy force due to the hydrostatic pressure is F p,i = −ρ f gV p,i , where g is the gravitational acceleration vector and V p πD 3 p /6 is the volume of the particle.

Model Setup and Model Validation
The aforementioned governing equations are solved in nondimensional forms. With the characteristic velocity scale U and length scale L, the non-dimensionalized variables are defined as follows: In homogeneous isotropic turbulence, the Reynolds number based on the Taylor microscale (λ) and the root mean square of turbulent velocity fluctuations (u rms ) is commonly used, which is defined as Re λ = λu rms /]. The Taylor microscale is computed by λ 15]/ϵ √ u rms , where ϵ is the viscous dissipation rate and computed from the simulation results. We varied the properties of primary particles, including the particle diameter (D p p D p /L), Young's modulus (Ep = E/ρU 2 ), viscoelastic damping coefficient (η p D η D /ρUL 2 ), and the surface energy density (γp = γ/ρU 2 L). Properties of primary particle used in Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 this study are summarized in Table 1. Coefficients in the tangential force and rolling friction models were kept the same for all simulations. The Young's modulus of soft clay particle is on the order of 1 MPa, which can significantly restrict the critical time step for DEM simulations. In practice, the Young's modulus used in the model is often several orders of magnitude smaller than the actual value to accelerate the computation. Tsuji et al. (1993) showed that stiffness can be reduced by orders of magnitude without altering the collisional behaviors of particles. In this study, the Young's modulus and the surface energy density are scaled down properly to make sure the relative importance of the elastic force and adhesive force is kept the same and the same floc structures can be reproduced. For simplicity, the superscript "*" in the non-dimensionalized variables are omitted in the analysis. The computational domain is a periodic box of size 8 × 8 × 8, and 16 elements of uniform size were used in each direction. A polynomial order P N = 8 was used within each element, which yields a total resolution of around 2.1 million grid points. The third order Adams-Bashforth method was used for time integration. A fixed time step was chosen in all simulations, which ensures the maximum Courant-Friedrichs-Lewy (CFL) number to remain around 0.2. The "3/2" rule was used for dealiasing. The PNPN-2 algorithm was applied (Maday et al., 1990;Fischer, 1997), in which pressure is solved on a coarser grid with polynomial order 6.
The DNS Model has been validated with previous DNS study by Rosales and Meneveau (2005) using the time-averaged energy spectrum. Due to the small Reynolds number used in both studies, there is no clear "−5/3" slope. The red solid curve represents the averaged non-dimensional energy spectra over cases with different Reynolds numbers. Our model results agree with the previous DNS study reasonably well ( Figure 1). The forcing coefficient (α = 0.033) and the viscosity (] = 5 × 10 -3 ) are kept the same for all cases. This gives the Reynolds number of 200 based on unit characteristic velocity and length scales, and Taylor Reynolds number (Re λ ) of 32. For homogeneous isotropic turbulent flow, there are no intrinsic characteristic scales for the mean flow. The Taylor Reynolds number is commonly used in homogeneous isotropic turbulence because it uses the fundamental length scale and velocity scale of turbulence to define the Reynolds number. To relate the idealized simulation conducted in this study with field condition, the Taylor Reynolds number can be used to obtain the turbulent shear rate when characteristic scales in dimensional form are given. The average Kolmogorov length scale is η = 0.049 and the average turbulent kinetic energy is k t = 0.12.
To make results relevant to geophysical or engineering applications, simulation results can be interpreted in the dimensional forms with given characteristic length scale and velocity scale. Due to the limitation of computational resources, the present study focuses on energetic environment with high turbulent shear only. For a characteristic length scale of L = 10 -3 m, the characteristic velocity scale is U = 0.2 m/s, based on the Reynolds number, and the particle diameter is D p = 12.8 μm for cases P 2 . The particles can be interpreted as the smallest clay-based aggregates, namely, flocculi. Flocculi seldom break down to the lowest-level primary particles even at the high turbulent shear and hence are the building blocks of large flocs. The turbulent shear rate based on the characteristic scales is 350 s −1 for all cases. The shear rate is higher than the values in most laboratory experiments and field observations; however, the model captures how turbulence affects flocculation dynamics. The Young's modulus for soft clay is in the range of 0.5 to 5MPa, and the Poisson ratio of clay is 0.3. The Young's modulus used in the simulation is between 6.4 and 16 kPa, which is reduced by two order of magnitude to accelerate the computation. In the JKR theory, the pull-off force to separate two particles is 3πγ JKR D p /2, which can be used to define the yield strength of the floc. The softness of particles does not affect the yield-strength of flocs directly and hence the aggregation of particles. Detailed measurement on the surface energy density γ is still lacking. In this study, the surface energy density used in the JKR model is in the range of 2 × 10 -4 to 1 × 10 -3 J/m 2 . The stickiness of the particle can be characterized by the adhesive number, which is defined as the ratio of the yield strength of flocs represented by the surface energy density γ to the turbulent kinetic energy Ad = γ/ ρk t D. Because the turbulent intensity remains the same in all cases, the average floc size increases with the adhesive number as expected when the primary particle is kept the same ( Table 1). In addition, the averaged floc size is almost three times greater than the Kolmogorov length scale for the cases with the largest adhesive number. However, for cases with relatively small adhesive number (case P 2 S 2 and P 2 S 3 ), the floc size is limited by the Kolmogorov length scale, and the average floc size is comparable to the Kolmogorov length scale.

Flow Visualization
The flow velocity field from case P 1 S 1 is shown in Figure 2A. The horizontal x − y plane is located at z 0 = 0. Only particles whose FIGURE 1 | Model validation with DNS simulation results from Rosales and Meneveau (2005).
Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 centers are located between z 0 − D p /2 and z 0 + D p /2 are projected to the x − y plane. The zoom in view is shown in Figure 2B to show the detailed structures of the flocs. Clusters of primary particles can be observed as flocs form. The Stokes numbers St sD p u rms * /18] are 0.2 and 0.13 for P 1 and P 2 cases with the specific gravity of s = 2.65, therefore we did not observe strong local preferential accumulation of particles.

Time Evolution of the Floc Size Distribution
All simulations are initialized with mono-dispersed spherical particles, which are uniformly distributed in the simulation domain. Initially, the particle velocity is set to zero. Due to collisions driven by turbulence, flocs start to form gradually and the floc size distribution evolves into an equilibrium distribution. In total, two contrasting cases, P 2 S 1 (softer and more sticky primary particles) and P 2 S 2 (more stiff and less sticky primary particles), are selected to investigate the dynamics. Figure 3 shows the time-evolution of the floc population (N n f ) from case P 2 S 1 , where the primary particles of size D p = 0.0128 are the stickiest (the adhesive number Ad is the greatest), and n f represents the number of primary particles consisting the floc. At the beginning stage (t = 2), a power law relation can be observed. With the formation of larger flocs (t = 1000-3000), the power-law relation can still be observed for small flocs, but the slope starts to decrease. The slope is significantly different from the beginning stage. At the early stage, small flocs form mainly due to aggregation, and the power-law relation can well capture the size distribution for those small flocs (Hunt, 1982). At the intermediate stage (t = 4000-6000), we observe the accumulation of intermediatesized flocs with n f between 20 and 90, which forms a staircase in the floc size distribution (for instance at t = 5000). At the late stage (t ≥ 7000), the population of floc (N n f ) with n f between 20 and 90 starts to decrease and a peak appears around n f = 95. N n f only changes slightly for relatively large flocs at t ≥ 10,000, suggesting the equilibrium floc size distribution is reached. The floc size distribution shows an asymmetric shape with respect to the peak of N n f at n f = 95 on the log-log plot. Figure 4 shows the time-evolution of floc size distribution from a contrasting case P 2 S 2 with less sticky and more stiff particles. Similarly, at the early stage (t ≤ 2000), the power-law relation between N n f and n f can be observed. However, there is no formation of the staircase-shaped structure at the intermediate stage, and the peak around n f = 16 shows up at a much earlier time and is evident for t ≥ 3000. In addition, the population (N n f ) of large flocs of size n f ≥ 16 does not change with time much. However, the depletions of small flocs of size n f ≤ 5 due to aggregation can still be observed at the late stage. The floc size distribution shows a more symmetric shape with respect to the peak at n f = 16 on the log-log plot compared to the case P 2 S 1 .
To further investigate the aggregation process of flocs, the time-evolution of the floc population for given floc sizes are shown in Figures 5, 6. For case P 2 S 1 , small flocs of size n f between 2 and 5 show a similar pattern ( Figure 5A). The population N n f first increases and reaches a peak, then it drops and approaches to an asymptotic value at the late stage when the equilibrium is reached. For relatively large flocs ( Figure 5B), the time of the first appearance of the floc of size n f increases with floc size n f , since the flocs are built gradually when the floc grows. For intermediate flocs of size 10 <n f < 40, we observed a similar pattern to small flocs, N n f increases, peaks, and then decreases and approaches to the asymptotic value. However, for large flocs of size n f ≥ 50, N n f increases to the peak value and then approaches to the asymptotic value. Continuous aggregation and breakup keep Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 occurring at the equilibrium stage. The oscillations in N n f at the late stage are due to the intermittent nature of turbulence. The time evolution of floc population N n f from case P 2 S 2 is shown in Figure 6. Similar patterns for small flocs of size n f ≤ 5 can be observed. However, large flocs of size n f ≥ 10 show a different pattern that N n f first increases and then approaches the asymptotic constant at the equilibrium stage. This is consistent with the time evolution of the floc size distribution that the peak appears in the early case P 2 S 2 and the population of large flocs does not change much with time at the late stage ( Figure 4). The time of the first appearance of large floc of n f ≥ 10 also increases with floc size, again suggesting the flocs grow gradually. Again, we observe oscillations of N n f at the equilibrium stage due to the intermittent nature of turbulence. The oscillation is much stronger for larger flocs because large flocs are more fragile and more susceptible to breakup by turbulent shear.  Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 7

Flocculation Dynamics: Breakup and Aggregation
To better understand how physicochemical properties of the primary particle (i.e, stickiness) affect the equilibrium floc size distribution, we examine the flocculation dynamics using the population balance equation as follows: where n (v, x, t) is the number density of flocs with volume (or size) v at time t and location x, W s is the floc settling velocity, u i is the fluid velocity component in the i-th direction, and κ is the sum of the molecular and turbulent diffusivity. On the right hand side of the equation, Q is the aggregation kernel and Γ is the breakup kernel. β is the fragmentation distribution, which describes the created number of daughter flocs of volume v after the breakage of a mother floc of volume v′. The aggregation kernel (Q) is a function of the collision frequency and collision efficiency. The collision frequency is a function of the turbulent shear rate and increases with the turbulent shear rate. The collision efficiency is defined as the rate of successful collisions resulting in the aggregation of flocs to the total number of collisions, which is a function of the properties of sediment particles.  Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 The first term (term I) on the right hand side represents the formation of a floc of volume v by the aggregation of two smaller flocs of volume v − v′ and v′. The second term (term II) on the right hand side represents the formation of a floc of volume v from the breakup of a larger floc of volume v′. In both term I and term II, a new floc of volume v is generated, and hence they are the construction terms. The third term (term III) represents the aggregation of a floc of volume v with another floc of volume v′ to form a larger floc of volume v + v′. The last term (term IV) represents the breakup of a floc of volume v. In both term III and term IV, a floc of volume v is consumed, and hence they are the destruction terms.
In Eq. 21, both aggregation and breakup processes in the population balance equation require parameterization, including the aggregation kernel, the breakup kernel, and the fragmentation distribution (Jeldres et al., 2015). By modeling the sediment phase as the dispersed phase using the Lagrangian framework, we can track the time-evolution of individual flocs of different size n f (or volume v = n f V p , with V p as the volume of the primary particle) to understand the aggregate and breakup processes at the particle scale. By comparing simulation results at two consecutive time instances, aggregation and breakup of flocs can be obtained in the time-driven Lagrangian model. Considering all the flocs of size n f = 16, we investigate the state of each one of them at a time interval Δt = 2 before. Most of the flocs of size n f = 16 have remained the same over this small time interval. Some flocs would have been of a smaller size (i.e., n f < 16) at the previous time (t − Δt) and have grown to flocs of n f = 16 due to aggregation, while some of the flocs would have been larger at the previous time and have reduced in size to n f = 16 due to breakup. We refer to the previous time floc size at t − Δt as the "prior-size". Figure 7A shows the probability density function (PDF) of the prior-size of flocs whose current size is n f = 16. Most of the flocs of prior-size n f = 16 that have remained the same without aggregation or breakup and are not included in the analysis. In Figure 7A, the blue circle symbols represent the source of n f = 16 flocs. The circles to the left of the dash line correspond to the PDF of smaller flocs aggregating and becoming n f = 16 floc, while circles to the right of the dash line correspond to the PDF of larger flocs breaking up and generating a daughter floc of size n f =16. These are terms I and II on the right hand side of (21).
In a similar manner, the red plus symbols represent the sink of n f = 16 flocs, i.e., they measure the PDF of what a floc of size n f = 16 floc becomes after a small time interval of Δt = 2. The pluses to the left of the dash line correspond to the PDF of smaller flocs that form from the breakup of n f = 16 flocs, while pluses to the right of the dash line correspond to the PDF of larger flocs that are formed by the aggregation of a floc of size n f = 16 with another floc (or other flocs). These are terms III and IV on the right hand side of (Eq. 21). The collapse of the two curves (circle sources and plus sinks) suggests a dynamic equilibrium with the balance between the aggregation and breakup processes. The PDF is almost uniform for small flocs. A peak is evident at n f ≈ 80. A power-law distribution of the PDF can be observed for the large flocs with n f > 80. For small flocs, we observe a drastic drop from the peak to n f ≈ 30, and the distribution is quite uniform for n f ≤ 30.
We carried out the same analysis for flocs of size n f = 32 ( Figure 7B), the power-law distribution is evident for large flocs of n f > 80. A uniform distribution can be observed between n f = 32 and 80. For smaller flocs of n f < 32, the distribution shows a minimum around n f ≈ 9 and peaks around n f = 32 and n f = 1 (primary particles). For floc of size n f = 64 ( Figure 7C), the power-law distribution is still evident for large flocs of n f > 80. In addition, a significant change of the slope for large flocs of size greater than n f = 150 can be identified. For small flocs of n f < 64, the distribution shows a uniform distribution between 20 and 60 and two peaks near n f = 64 and n f = 1 (primary particles). For case P 2 S 2 with less sticky and more stiff primary particles, the terms are plotted for flocs of size n f = 8, 16, and 24 ( Figure 8). The PDFs for flocs of different size are quite similar. It shows a power-law relation for large flocs of n f ≥ 16. For smaller flocs of n f < 16, the distribution exhibits a minimum at n f = 5 and two peaks at n f = 16 and n f = 1 (primary particle). Compared to the case P 2 S 1 with the stickies primary particles, the presence of the uniform distribution for intermediate-sized flocs (n f between 30 and 50 in case P 2 S 1 ) is not evident in case P 2 S 2 . In both cases, at the equilibrium stage when the breakup and aggregation processes balance with each other, simulation results show that the construction by aggregation is primarily balanced by the destruction by breakup and the construction by breakup is primarily balanced by the destruction by aggregation. At the microscopic level, each aggregation or breakup pathway is reversible and hence in a dynamic equilibrium. With given aggregation and breakup kernels, equilibrium solutions of the floc size distribution exist (Vigil, 2009), and the mathematical form of the equilibrium floc size distribution could be derived.

Equilibrium Floc Size Distribution
An equilibrium floc size distribution develops when the aggregation process balances with the breakup process. The floc size distribution is modeled as a function of n f , instead of the floc size D f , because the number density (n) in Equation (21) is expressed as a function of the floc volume v. In this study, flocs are consisted of slightly overlapping identical spheres, and the volume of a floc consisted of n f primary particles can be approximated by n f as v ≈ n f πD 3 p /6. To calculate the floc volume using the actual floc size D f requires a priori knowledge of the floc internal structure, which is difficult to obtain. Because of the asymmetry around the peak, we tested the floc size distribution against two widely used asymmetric distributions, namely, the lognormal distribution and the Weibull distribution. The lognormal probability distribution function is expressed as follows: where ln n f follows the normal distribution, μ is the mean, and σ 2 is the variance. The Weibull distribution is a special form of Gamma distribution with two parameters, namely, the scaling parameter λ and the shape parameter k.
The Weibull distribution interpolates between the exponential distribution and Rayleigh distribution. The shape parameter k affects the shape of the distribution rather than simply shifting or stretching it. Figure 9 shows the curve fitting for case P 2 S 1 (panel a), P 2 S 2 (panel b), and P 2 S 3 (panel c). For case P 2 S 1 with the stickiest primary particles, model results fit better with the Weibull distribution, while results from case P 2 S 2 with less sticky but more stiff primary particles fit better with the lognormal distribution. For case P 2 S 3 , neither lognormal nor Weibull distribution can fit the data for the entire range of the floc size n f . To assess the performance of different distributions, the Anderson-Darling (AD) test (Anderson and Darling, 1952) was conducted, which is based on the empirical cumulative distribution obtained from the sample. The AD test is commonly used to test if a sample of data comes from a population with a specific distribution. We used the significant level of α = 0.01, which is commonly used in statistical hypothesis test (Fisher, 1992). The results are summarized in Table 2, where E s is the sum of squared residual errors. The accepted hypothesis for each case is shown with "*" in Table 2, which means the optimal descriptor for the floc size distribution. For the case with less sticky primary particles (case P 2 S 2 ), the lognormal distribution fits better and for the cases with very sticky primary particles (case P 1 S 2 and P 2 S 1 ), the Weibull distribution fits better. The AD test rejects both lognormal and Weibull distribution hypothesis for most cases, suggesting neither lognormal nor Weibull distribution can accurately predict the floc size distribution. For instance, the lognormal distribution fits better for small flocs and also captures the peak more accurately in case P 2 S 3 ( Figure 9C), while the Weibull distribution fits better for large flocs (n f > 35). The adjusted coefficient of determination (R 2 adj ) was then used to evaluate the goodness of the fit (Ezekiel, 1930). However, the R 2 adj for both distributions are quite close, and hence it is difficult to distinguish the two distributions. Based on curve fitting results, the shape parameter k is around 2.5 for all cases, suggesting a similarity in the floc size distribution.
Floc size distribution from P 2 S 1 fits better with the Weibull distribution. To further investigate the floc size distribution from case P 2 S 1 , we plotted model results under log-log scale Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 ( Figure 10A) and semi-log scale ( Figure 10B). We averaged N n f from t = 11,000 to 12,000, when the equilibrium is achieved for large flocs. For flocs with size n f between 30 and 75, a power-law relation can be identified when the aggregation and breakup processes balance. For large flocs with n f greater than 110, we observe an exponential decay of N n f with n f . For case P 2 S 2 with lognormal distribution, we did not observe the power-law relation, and hence the results are not shown.

DISCUSSION
Different mathematical formulas of floc size distribution arise in the aggregation or breakage processes (Huber et al., 2006). Lognormal distribution has been observed in particle growth or coagulation processes (Smoluchowski, 1918;Friedlander and Wang, 1966), in which aggregation process dominates the dynamics. On the other hand, Weibull distribution has been commonly observed in the fragmentation process of large   particles (Brown and Wohletz, 1995). In the flocculation process, both the aggregation and fragmentation processes play an important role. At the equilibrium state, a floc can be constructed either by aggregation of smaller flocs or breakup from larger flocs. To further understand under which circumstances a lognormal or Weibull distribution performs better, we analyzed the dominant floc construction mechanisms at the equilibrium stage. Figure 11 shows the relative importance of construction by aggregation and construction by the breakup for flocs of size n f from the two contrasting cases P 2 S 1 and P 2 S 2 . In case P 2 S 1 ( Figure 11A), the majority of small flocs are constructed by breakup of larger flocs, while large flocs (n f > 75) are mainly constructed by aggregation as expected. The aggregation and breakup processes play equally important roles for flocs in the range of 30<n f < 75. The primary particles first aggregate into microflocs, the microflocs are quite resilient to turbulent shear and serve as the building blocks for larger flocs. For case P 2 S 2 ( Figure 11B), we observe a monotonic increase in relative importance for aggregation and a decrease for breakup process with respect to n f . In general, a large portion of flocs (n f < 75) in case P 2 S 1 are generated mainly from breakup of larger flocs. However, in case P 2 S 2 , breakup only controls the formation of a small portion of flocs with n f < 15, and the aggregation process dominates the formation of flocs for a  Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 wide range of flocs. The lognormal distribution (case P 2 S 2 ) develops when the large flocs grow gradually from primary particles without the intermediate stage (the formation of "microflocs"). This behavior has not been observed in laboratory experiments before, and future laboratory studies are required to test whether this is an artifact of the numerical model or an actual physical process. Mathematical approach of floc modeling that has gained much interest by modelers is the fractal representation of flocs (Kranenburg, 1994;Merckelbach, 2000;Graham and Manning, 2007). Fractal theory is dependent on the successive aggregations of self-similar flocs, thereby producing a structure that is independent of the scale (or scale invariant). This is similar to the order-of-aggregation theory (Krone, 1963). Following the fractal theory, simple power laws can be used to describe floc properties such as floc density and settling velocity, as well as the aggregation and breakup processes (Winterwerp, 1998). Although some studies suggest that individually some natural muddy flocs (particularly those with high organic contents) may not be fully fractal in structure (Zhang et al., 2018;Spencer et al., 2021), the wider examination of in situ floc populations shows that a fractal representation of flocs still has many merits Winterwerp et al., 2006). The fractal dimension d 0 used to characterize the floc structure is defined as where n f is the total number of the primary particles consisting the floc, D f is the floc size and D p is the primary particle diameter. The major axis length (longest axis) is used as floc size D f , which is obtained by the principal component analysis (PCA). In general, the fractal dimension (d 0 ) is 1 for chain-like flocs and 2 for flat plane-like flocs. Flocs with fractal dimension close to 3 have compact structure and spherical shape. Flocs with the same number of primary particles (n f ) can exhibit different structures, and the averaged fractal dimension (d 0 ) for flocs with the same n f is shown in Figure 12. For case P 2 S 1 , the fractal dimension first increases and reaches a constant value around 2.4 for flocs of size n f between 50 and 120, suggesting compact and similar floc structures. Fractal dimension d 0 decreases for large flocs of n f > 120. The increase in d 0 for small flocs is due to the limited configurations of floc structure by finite n f . For instance, the most compact structure for a dimer (aggregate consisted of exact two spheres) is a rod with a fractal dimension of 1, and the most compact structure for a trimer (aggregate consisted of exact three spheres) is an equilateral triangle of fractal dimension of 1.7. The "microflocs" in the range of 30<n f <75 ( Figure 11A) have relatively small fractal dimensions. The primary particles in case P 2 S 1 are the stickiest and the turbulent shear stress is less efficient to break these 'microflocs' at the scale. The decrease of d 0 for large flocs suggests that they are more porous. Similar trends can be observed for cases P 1 S 2 , in which Weibull distribution better describes the flocs size distribution.
In contrast, in case P 2 S 2 with the least sticky primary particles, the fractal dimension increases to the peak value around 1.95 and then drops for large flocs. The decrease in the fractal dimension for large flocs has also been observed by Khelifa and Hill (2006) and Maggi (2007). Our numerical results suggest the scaledependence of floc structure as the fractal dimension is not constant for the entire range of the floc size. A variable fractal dimension should be considered to characterize the flocs.
In general, for the cases with the same primary particle diameter, primary particles with small adhesive numbers (or surface energy density) lead to flocs with smaller fractal Frontiers in Earth Science | www.frontiersin.org July 2022 | Volume 10 | Article 815652 dimensions, and primary particles with large adhesive numbers generate more compact flocs with fractal dimension as large as 2.4 in case P 2 S 1 . Floc compaction by the breakage-regrowth and restructuring mechanisms are more pronounced in the case with the stickiest primary particle. By comparing case P 2 S 2 and P 2 S 3 , the stiffness of the particle also affects the floc structures. Primary particles with larger Young's modulus (case P 2 S 2 ) lead to flocs with smaller fractal dimension.

CONCLUSION
A two-phase Eulerian-Lagrangian framework was implemented to investigate the equilibrium floc size distribution of cohesive sediment in homogeneous isotropic turbulence. The primary particles are modeled as identical sticky soft spheres, and particle-particle interactions are modeled by the discrete element method. The adhesive contact JKR model was implemented to model cohesive sediment particles, which is a tensile force model with hysteretic effect. In the adhesive contact model, the pull-off force to break two particles apart scales with both the particle size and the surface energy density (i.e., the physicochemical properties of the primary particle). A series of numerical simulations were conducted by varying the size and properties of the primary particles. At the equilibrium state, the construction by breakup is balanced with the destruction by aggregation, and the construction by aggregation is balanced with the destruction by breakup. The equilibrium floc size distribution depends on primary particle properties, including the stiffness and the surface energy density. For cases with more sticky primary particles, the floc size distribution can be better described by the Weibull distribution with a shape parameter around 2.5. In addition, at the intermediate stage, a staircase structure develops in the floc size distribution. The primary particles first form the 'microflocs', which serve as the building blocks for large flocs. For the case with less sticky primary particles, the lognormal distribution performs better. Flocs grow gradually from primary particles without the intermediate stage of 'microflocs'.
By analyzing the construction mechanisms of flocs of different size, when the Weibull distribution develops, construction by breakup and construction by aggregation are of equal importance for the intermediate-sized flocs. The fractal dimension of large flocs then decreases with floc size, suggesting large "macroflocs" are more porous and fragile. For less sticky particles, the lognormal distribution develops, and the aggregation dominates the floc construction for a wide range of flocs. The fractal dimension of flocs first increases with floc size, reaches the peak value, and then decreases with the floc size. However, given the similarity between the lognormal and Weibull distributions and hence the difficulties in distinguishing between them in confidence, it is recommended to choose the floc size distribution and make interpretations in practice with caution.
Due to the limited computational resources, the current simulation focuses on the high-energy environments with large turbulent shear rate (350 s −1 in this study). Simulations with more particles (several millions to billions of particles) are therefore required for more realistic cohesive sediment transport studies in low-to moderate-energy environments. In addition, current model framework oversimplifies the hydrodynamic interactions among particles without the influence from neighboring particles. For cohesive sediment, the particle Reynolds number based on Stokes settling velocity is small, and hence the sheltering and blockage effects from neighboring particles could play an important role. The sheltering effects from neighboring particles lead to reduced hydrodynamic drag, and hence could affect the breakup processes. A more sophisticated efficient model that can accurately predict hydrodynamic interactions among a large amount of particles is required and will be the future work.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://figshare.com/ articles/dataset/FlocSizeDistribution/16619734.

AUTHOR CONTRIBUTIONS
MY: methodology; formal analysis; and writing-original draft, review, and editing. XY: supervision; conceptualization; methodology; formal analysis; writing-original draft, review, and editing. SB: conceptualization; methodology, and writing-review and editing. AM: writing-review and editing.

ACKNOWLEDGMENTS
All simulations were carried out on the HiperGator super computer at the University of Florida. Computer resources, technical expertise, and assistance provided by the HiperGator are gratefully acknowledged. The CFD model nek5000 is developed primarily by the Argonne National Laboratory at https://nek5000.mcs.anl.gov/, and the DEM model LAMMPS is developed by Sandia National Laboratories at https://lammps. sandia.gov/.