Modelling flocculation Towards an integration in large-scale sediment transport models

Despite recent advancement in the field of sediment transport, the integration of cohesive sediment properties in large-scale transport models remain a challenging task. In order to model adequately the change in particle size that occurs in different environmental conditions, flocculation models based on the so-called Population Balance Equations (PBE) are often used. These models have to be efficient enough to be implemented in numerical transport models, and as full PBE ’ s are time-expensive to run and depend on a huge amount of a-priori unknown parameters, simplifications have to be made. These simplifications comes unavoidably at the cost of properly accounting for the complex particle-particle and particle-fluid interactions. In order to stay as close as possible to the physical processes, we propose a different approach based on a logistic growth model that mimics the Particle Size Distribution (PSD) measured over time for all size classes. The parameters of the model can easily be found from laboratory measurements. In contrast to most models, the particle classes we propose are not defined by particle size, but in terms of mineral sediment composition. One class is composed of (unflocculated) mineral sediment particles, another of flocculated sediment particles and a third one of organic particles. The mass balance between classes and the way to obtain their corresponding average settling velocity are given. Mass balance and settling velocities are the required input parameters for all sediment transport models. The simplicity of the derived expressions, and their link with measurable variables, makes them good candidates for future implementation in such models.


Introduction
Large-scale numerical models, such as DELFT3D (Lesser et al., 2000), (Lesser et al., 2004), ECOMSED (Blumberg and Mellor, 1987), (Blumberg et al., 1996) and TELEMAC (Hervouet, 1991), (Le Normant, 2000), (Sherwood et al., 2018) enable to study the transport of sediment and from this derive their impact on the environment. In the regions where the sediment is fine-grained and prone to interact with other particles or the surrounding fluid, these models were until recently limited as they did not account properly for these interactions (Droppo, 2006).
Because of the interactions, the particles present in the water column will change size, and hence the concentrations of the different size classes will vary over time. This has important consequences in terms of sediment transport, as by aggregation (flocculation) for example, the concentration of a large size class will increase while the concentration of a small size class will decrease accordingly. Realizing that every size class has a different mineral sediment content, this implies that the spreading of mineral sediment within the water column will be very different from the case where no aggregation occur. Another consequence of flocculation is that a change in particle settling velocity will occur. Aggregates are principally composed of mineral sediment particles, organic matter and voids. As the density of organic matter and the voids (filled with water) are lower than the density of mineral sediment, for a same size, the settling velocity of a mineral particle will be much larger than the settling velocity of an organic floc. Therefore mineral sediment particles will be transported less far than organic flocs of same size under the same conditions.
The most commonly used model to describe the aggregation of colloidal particles in natural systems are the Population Balance Equations (PBE) (Maggi, 2009), (Mietta et al., 2011), (Shen et al., 2018), (Lai et al., 2018). This set of equations describe the time evolution of the concentration of particles in a size class, where a size class is defined as the ensemble of particles of a given size and density. A PBE equation for a size class k (there are as many PBE's as there are size classes) is given by the time derivative of the mass, volume or number concentration of particles in the size class. This derivative is function of the sum of four terms, representing the gain by aggregation (particles of smaller sizes aggregating to form a particle of size k), the gain by break-up (particles of higher classes breaking to form particles of size k), the loss by aggregation (particles of size k aggregating and therefore leaving class k) and loss by break-up (particles of size k breaking and therefore leaving class k). These terms depend in particular on four parameters: the collision efficiency, the collision frequency, the break-up rate and the break-up distribution function. In principle, each of these parameters is class-dependent.
In controlled experiments, where the aggregation of monodisperse colloidal spheres in the presence of salt is studied, the terms can be estimated and the PBE's then reduce to Smoluchowski' s set of equation. It has then been verified that theory and experiments are in good agreement (Cahill et al., 1987), (Fernández-Barbero et al., 1996), (Russel et al., 1991). It has also been verified that clay flocculated in the presence of salt can relatively well be modelled using PBE's, despite the particles being polydisperse and having complex surface properties (Ives, 1978).
In natural systems, and especially in coastal areas which are the topic of many researches, clay (inorganic) particles interact with organic matter (Avnimelech et al., 1982), (Eisma and Irion, 1993), (Guenther and Bozelli, 2004). Schematically, organic matter is composed of living organisms (algae, bacteria…) and organic debris (dead organisms or their remains, polyelectrolytes such as polysaccharides…). In the present study we will use as an example the flocculation of clay by a polyelectrolyte. For convenience, we will study the flocculation of clay in the presence of a commercially available cationic polyelectrolyte. Using a cationic polyelectrolyte has the advantage that the interaction between the negatively charged clay and positively charged polyelectrolyte is driven by a strong Coulombic attraction between the two, which is weakly depending on the characteristics of the solvent. In nature, most polysaccharides are negatively charged (Karunaratne, 2012) and therefore solvent properties like salinity are important parameters that influence the flocculation ability.
Flocculation in the presence of polyelectrolyte is dominated by the characteristics of the polyelectrolyte (and the solvent). A major issue in electrolyte-induced flocculation is that, in contrast to salt-induced flocculation, the formed flocs have a different relation to shear as, where the size of salt-created flocs follows the Kolmogorov microscale upon increase and decrease of shear (Mietta et al., 2009a(Mietta et al., , 2009b, polymeric flocs do not relate to the microscale in a simple way (Ibanez Sanz, 2018). Upon increasing shear, the flexible polymer chains will tend to coil and the floc's density will increase while its volume decreases (Adachi et al., 2012). When the shear is decreased, the polymer will not fully disentangle. Flocs in a natural environment are therefore more complicated to model as they do not solely depend on aggregation/break-up mechanisms, even though the detachment of small flocs and particles from a large floc under shear is not to be excluded.
In recent years, simplified PBE modelsto be incorporated in largescale transport modelshave been proposed to study the behavior of suspended particulate matter (SPM) in-situ. In particular a tri-modal flocculation model has been proposed (Shen et al., 2018) where three types of flocs are considered, termed microflocs (of size 3 μm and 15 μm), macroflocs (of size in the range 20-200 μm) and megaflocs (of size about 360 μm). Even though the tri-modal flocculation model has been shown to be an improvement on lesser modal PBE's, it is not certain, considering the role of organic matter on the flocculation process, that it captures the physics accurately. Moreover, even though the tri-modal model has successfully been implemented in sediment transport model TELEMAC, it still requires some computational efforts.
The aim of the present article is to propose a new approach, based on a generic formulation for floc size growth using logistic functions (see Section 2). The main advantage in using these functions is that they can easily be calibrated for in-situ conditions and parameterized using laboratory experiments. The number of adjustable parameters in the sediment transport model is hereby virtually zero in the case of laboratory experiments, since all the parameters can be assessed by independent measurements. In Section 3 the different classes of particles are defined and the sediment transport model accounting for mineral sediment mass transfer between classes is presented. In Section 4 two examples are given. The first one shows how the model can be parameterized from laboratory studies and the second how the model can be linked to a study that presents generalized formulations for the settling velocities of flocs. These expressions have been validated for the study of three Northern European estuaries.

Studying flocculation using a logistic growth model
In this section, we start with presenting the logistic growth model that we propose to use in replacement for the Population Balance Equation (PBE) that is generally used to describe flocculation. We explain why using a logistic growth model is a good alternative to the PBE in Subsection 2.3.

Logistic growth
The model proposed in the present article is adapted from the logistic growth theory derived by Verhulst in 1838 (Verhulst, 1845). The logistic growth model was originally meant to describe the time evolution of a population (the number of inhabitants in Belgium). This model has later been applied to other systems, and used for example to estimate the autocatalysis rate in chemical reactions or the growth of bacteria (Cunningham and Cunningham, 2002), (Tsai, 2005), (Chen et al., 2017).
Two terms, representing the growth and decay of a given population are given in the equation. We define N as the number of particles of a given sizein other words N represents the population of a given class whereby the class is defined as the ensemble of particles having the same size. The change in number of particles N in a class (the change in population) is given by where τ and a are constants, τ being a characteristic time and a the coefficient which can be related to the rate constant k = a/τ. For 1 ≫ aN the equation reduces to which represents the growth term, and has as solution where we have defined N(t = 0) = N 0 . The characteristic time of growth is given by τ.. This growth model is sometimes called a Malthusian growth model after Thomas R. Malthus who introduced it in 1798 (Malthus, 1986). On the other hand, for 1 ≪ aN and using k = a/τ the equation reduces to which is to be compared with the derivation done by Smoluchowski to represent the decay in population of a class of primary particles at the early stage of aggregation (Russel et al., 1991). The solution of this equation is given by The characteristic time of decay is given by τ decay = 1/(kN 0 ) = τ/(aN 0 ). The full solution of the equation (accounting for both growth and decay) is given by the sigmoid function where N ∞ = N(t → ∞) = 1/a. One can verify that the full solution reduces to the particular ones in their range of validity. We have τ decay /τ = N ∞ /N 0 . Eq.(6) is usually the model adopted to study the growth of bacteria in different bed layers (Chen et al., 2017).

Logistic function for flocculation
Other logistic functions can be defined, based on the same principle that one term controls the growth (or birth), symbolized by the function b(t) and another term the decay, symbolized by the function d(t).The general expression for the change in population over time is then given by The (positive) functions for birth and decay proposed in the present article are given by: It is emphasized that the names "birth" and "decay" are solely chosen because they appear with a plus and minus sign in the balance equation Eq. (7). One can easily show that considering only the birth function (d (t) = 0)) in fact gives as solution a sigmoid function similar to Eq. (6). The full solution of the change in population equation is given by In the case that t b = t d = τ the solution can also be written The right-hand side of Eq. (11) is a sigmoid function similar to Eq. (6).

Comparison with the Population Balance Equation
The Population Balance Equation (PBE) is defined as the one used commonly in sediment dynamics to describe the time evolution of a population of flocs (Mietta et al., 2009a(Mietta et al., , 2009b, (Mietta et al., 2011), (Shen et al., 2018), (Lai et al., 2018). The equation representing the population evolution in time then can be reduced to the sum of two terms (given in square brackets), the first one being representing the changes in population due to aggregation and the second one the changes due to break-up: The break-up rate is given by s i and break-up distribution function by γ i, j . The term involving break-up is related to the property of a single particle to be ableor notto break and is therefore proportional to the number of particles N i in a class. The term corresponding to aggregation on the other hand depends on the product of the number of particles in two classes, as it depends on the interaction between two particles. The collision efficiency α i, j is usually taken to be a constant, i.e. α i, j = α and the collision frequency due to shear can be estimated by (Russel et al., 1991) where G is the shear rate and L i is the diameter of a particle of class i. From a simple dimension analysis, it follows that the characteristic aggregation (birth) and break-up (decay) rates are proportional to These rates are to be compared with the ones found from a dimension analysis of Eqs. (8), (9): A major difference is found between the aggregation rate as k b does depend on class j in the PBE model. The collision frequency β i, j in the PBE model is derived under the assumption that the collision frequency is depending on the size of the particle that the particle of class i is colliding with (large particles collide more frequently than smaller particles). This expression for β i, j and the associated PBE model has originally been set-up and validated for hard spheres aggregation in an electrolyte (Russel et al., 1991). In the case of flocculation induced by organic matter, the considered particle of class i can be colliding with an organic particle of very anisotropic shape. For instance, it is well-know that in-situ flocculation is driven by polymeric substances (Van Leussen, 1988), (Ingerson, 1960). Polymeric chains are very flexible and their radius of gyration is a function of shear and solvent properties (like salinity). It is therefore likely that the expression for β i, j should be adapted in this case. The collision efficiency α i, j is similarly depending on the presence of organic matter and the way it interacts with the clay. The collision efficiencies between two sediment particles, two flocs or a floc with organic matter are very different, due to the interactions involved. The convention to take α i, j constant (α i, j = α) is consequently questionable.
The main advantage of the PBE model is nonetheless that conservation of mineral sediment mass is, by construction, obeyed: all particles in a class have a given particle size and are composed of a given number of primary mineral sediment particles. When a particle aggregates or breaks, its primary particles are redistributed in other classes. One practical disadvantage is that depending on the modes of aggregation and break-up a significant amount of classes have to be created. This is the reason why, in most models, break-up functions are chosen such that a particle can only break in a limited amount of classes. In earlier models (Mietta et al., 2009a(Mietta et al., , 2009b, binary break-up was therefore chosen (a particle breaks in two parts of same size), but this solution is clearly unrealistic.
In-situ, a driver for break-up is shear rate, and as the shear rate is usually changing gradually, most flocs will adapt likewise. It has been shown by earlier studies that in-situ floc sizes follow relatively well the Kolmogorov microscale: when the shear increases, the mean floc size will decrease by either erosion of small flocs from the outer part of the floc or by reconformation of the floc due to the coiling of organic matter (Verney et al., 2009).
The binary break-up condition has been relaxed in recent models, but this then leads to a significant amount of adjustable parameters. In (Shen et al., 2018) for example, where only 3 classes of particles are considered (called microflocs, macroflocs and megaflocs) the mass transfer of one class to another by break-up required as parameters: a breakup fitting parameter, the mass fraction of microflocs produced when a macrofloc breaks up, the mass fraction of microflocs produced when a megafloc breaks up, the mass fraction of the remaining megafloc when a larger megafloc breaks up, the number of generated macroflocs when a larger macrofloc breaks up and the number of produced macroflocs when a megafloc breaks up (6 adjustable parameters). All these parameters have to be calibrated with in-situ data and cannot easily be parameterized.
The advantage of the logistic growth approach is that each class is (mathematically) independent of another. There is therefore no need to set-up a given amount of classes to account for flocs that change size during the flocculation (birth)/break-up (decay) processes. In short, the standard PBE model, which has originally been set-up and validated for hard spheres aggregation in an electrolyte (Russel et al., 1991), does not offer a better physical insight than the model proposed here to describe the interaction between particles in-situ.
We discuss in the next section (Section 3) three main classes that could be considered in most situations, even though, as will be discussed in Section 4.2.2 there can be situations where one class might be sufficient. The interaction between particles is accounted for through the rates k b and k d which can implicitly depend on the presence and properties of other classes of particles. The rates k b and k d can be space and time dependent. These rates can be quantified by experiments performed in the laboratory, as will be shown in Section 4. By performing systematic series of measurements, k b and k d can in the future be parameterized as function of relevant variables, such as shear, salinity, organic matter type and concentration. A point of concern is that mass conservation should be accounted for between classes, as this is not automatically accounted for as in the PBE model. In the next section (Section 3) we will show how this can be done.

Defining classes and their properties
In contrast to previous work, we here will define particle classes based on particle composition (and not on size). Three major classes are defined. Each class can contain sub-classes and these sub-classes have to be chosen depending on the system investigated.

Class 1: fine mineral sediment
Class 1 contains only mineral sediment particles. If required, this class can be split into 3 fractions (clay, silt and sand). Both clay and silt particles can be part of flocculated structures which transport over significant distances, but the size and density of sand particles usually prevent these particles to be advected over large distances or to be part of flocs.
Class 1 is associated with a (time-dependent) total mineral mass concentration m 1 and an average settling velocity ω s, 1 .

Class 2: flocculated mineral sediment particles
Class 2 particles represent flocs. In contrast to other models, when Class 2 particles "break" (restructure) or aggregate they remain in Class 2: only the D50 of that class (which is usually also representative for the D50 of the whole PSD at steady-state) will change. Class 2 is associated with a (time-dependent) total mineral mass concentration m 2 and an average settling velocity ω s, 2 .
There are mainly two types of flocculation to be discussed: (a) salt-induced flocculation, which is occurring in coastal regions where fresh river water comes into contact with saline sea water. The aggregation between clay particles is driven by salt ions, thanks to the compression of their electric double layer and the role of Van der Waals attraction. The density of the flocs depend on the interaction between particles, and it has then been proven experimentally that the flocs such obtained have a fractal structure (Schaefer et al., 1984), (Russel et al., 1991). The density of salt-induced flocs can then be modelled using the relation where ρ floc is the density of the flocculated particle, ρ w the density of the suspending fluid (water), ρ p and D p the density and diameter of the primary particles that constitute the floc. The fractal dimension n depends on the type of aggregation, that can be either diffusion-limited aggregation (DLA, α ≃ 1) or reaction-limited aggregation (RLA, α < 1). The fractal dimension for DLA flocs is lower than for RLA flocs, as primary particles have to find an energetically favorable configuration to attach to a floc in the RLA case (Russel et al., 1991). It is for this case that the Population Balance Equation, and its simplest form (The Smoluchowski relations) has been developed and validated (Cahill et al., 1987), (Fernández-Barbero et al., 1996). Salt-induced flocculation/aggregation (often called coagulation) requires much longer time to reach a steady-state at given environmental conditions than most organicinduced flocculation. Salt-induced flocs follow the Kolmogorov microscale, and the aggregation/break-up mechanisms are reversible upon a change in shear (Mietta et al., 2009b).
(b) organic-based flocculation, which requires the interaction between Class 1 and organic particles. Organic based flocculation is dependent in majority on seasonality and salinity: different microorganisms live in fresh or salt water and the characteristics of the organic polymeric substances they produce are depending on ionic strength. Polyelectrolytes can for instance better coil in saline water due to the screening of their surface charge and this will influence the way they interact with clay particles. As most organic matter particles are overall negatively charged, they require cationic bridging to bind to the negatively charged clay. This is why flocculation is mainly diffusion limited in a saline environment. Due to their properties (flexible polymeric chains, affinity with specific surfaces) organic matter-induced flocculation results in larger and stronger flocs than salt-induced ones. The large radius of gyration of polymeric organic matter also ensure that this type of flocculation is very efficient. Due to the nature of organic particles, it is not possible to derive adequate aggregation and break-up parameters for a PBE model as discussed in Section 2.
It has been shown that the mean floc size in the turbidity maximum zone of the Yangtse estuary is corresponding to its steady-state size found in laboratory experiments, where the flocculation between mineral Yangtse sediment and algae was studied (Deng et al., 2019). A large study on flocs from three main European estuaries, discussed in more details in Section 4.2.2, shows by another arguments that the flocculation kinetics are usually fast and the system always near steady-state (Soulsby et al., 2013). These results are not so surprising as flocculation between mineral sediment and organic matter (especially polyelectrolytes) in a marine environment is usually only diffusion limited, and that the residence time of organic matter of low density is very high. In most situations therefore (when no major input of unflocculated mineral sediment is observed) the aquatic system is close to steady-state. In that case, sediment dynamics in the water column are driven by transport only, even though some care should be taken on the resuspension term, see Eq. (25).

Class 3: organic matter
Class 3 particles are composed of organic particles. Organic particles in situ can be of very different nature, shape, size and sticking properties, ranging from debris to living microorganisms. In (Deng et al., 2019), it was shown that algae cells could flocculated with themselves over time (this implies a change in size of Class 3 particles). In algae-rich regions (like in the surface waters in summer), these particles might be dominant over mineral-sediment rich particles.
The largest Class 3 particles (> 200 μm) have, due to their low density, a settling velocity of the order of 1 or 2 mm/s, independently of environmental conditions. These are particles that can be observed by in-situ systems such as under water microscopy or LISST (Safar et al., 2020). The LISST is a submersible laser-diffraction based particle size analyzer commonly used in-situ (Sequoia Scientific, 2011). Particles larger than 200 μm are usually only organic based. Particles of this size containing a significant amount of mineral sediment can only appear in the water column under strong hydrodynamic conditions (storms). In the 100-200 μm size range it becomes problematic to discriminate from LISST data between Class 2 flocs and Class 3 particles, therefore additional measurements (analysis of bottle samples) are then required.
Smaller Class 3 particles exist, but their properties (low density, transparency and shape) make them difficult or impossible (for EPS in particular) to be detected by current light scattering or video techniques.

Mass transfer between Class 1 and Class 2
From a sediment modelling perspective, only the transfer of mineral sediment mass is of importance. Defining a Class 3 particles containing no mineral sediment could therefore seem superfluous, but as this type of particles appear in the measured in-situ data, we will see how this class influences the parametrization and calibration of the model. In particular Class 3 particles can contribute to the turbidity of the water column. For the remainder of this section, we will refer to Class 3 as "flocculant".
In Fig. 1 a sketch representation of flocculation is given. The graphs schematically represent LISST data (% volume of size detected). The transfer of mineral sediment mass between Class 1 and Class 2 particles within a closed volume V tot (no particles can enter or escape this volume) can mathematically be expressed as: where m 1 = M 1 /V tot is the total mass of mineral sediment per unit of volume belonging to Class 1 particles (unflocculated particles), m 2 = M 2 /V tot is the total mass of mineral sediment per unit of volume belonging to Class 2 particles (flocculated particles). We have M = M 1 + M 2 where M is the total mass of (dry) mineral sediment in V tot . We emphasize that despite having the same units ( of components other than mineral sediment (ex: water, organic matter) in the total volume V k . The rate constant k m depends on variables such as shear G, m 1 , mass of flocculating agent m OM per volume and other possible variables such as salinity, pH, temperature, seasonality. The total mineral sediment mass within a given volume is The value of m is a constant in a closed vessel. In the experiments performed in the laboratory m is simply the mass of dry mineral sediment particles divided by the volume of water added to the 1 L jar used in the test.
The rate constant k m can be assessed from laboratory PSD experiments, where unflocculated clay is set in presence of flocculant. Class 1 particles have the particularity to obey the relation ρ 1 = ρ m, 1 (ρ o, 1 = 0) where the density ρ 1 is be taken equal to the density of a constant mineral sediment ρ s . As the density of Class 1 particles is constant, the change in volume V 1 of Class 1 is proportional to the mass concentration m 1 that left the class: All measured PSD's (in the laboratory or in-situ) are given in relative volume and therefore the volumes V k are given as which is the volume of Class k particles divided by the total volume (V 1 + V 2 + V 3 ) of particles detected in time (The % indicates that V k, % is usually given in % volume). The total relative volume (V 1 + V 2 + V 3 )/ V tot of particles detected is available from in-situ LISST, or by Optical Backscatter point Sensor (OBS). Note however that these measurements can be imprecise as (a) organic matter is relatively transparent for light scattering techniques and (b) organic-based flocs can be elongated, which has two major consequences. One is that the conversion from raw data into PSD, which is done using the assumption that particles are spherical, gives inaccurate results (Safar et al., 2020). The other is that the estimated volume of particles detected, which is also based on the sphericity assumption, similarly can give inadequate estimations. We define (The % indicates that V LISST, % is usually given in % volume). The total volume (V 1 + V 2 + V 3 ) varies in time. It follows that Fig. 1. Schematic view of the flocculation principle between Class 1 (mineral sediment, black bullets) and Class 3 (flocculant, black spaghetti) particles: when put in presence, Class 1 particles will flocculate to create Class 2 particles (grey circles). The average size of Class 2 particles during the flocculation process is depending on the flocculation kinetics. At steady-state the size of Class 2 flocs is representative for the D50 particle size distribution. The value of the D50 is highly depending on parameters such as shear rate, organic matter composition, salinity. Note that Class 3 particles cannot always been visualized correctly using light scattering devices. This is why their PSD is given in dashed curves. In the illustration, at steady-state Class 1 has been depleted: this occurs mainly at or above optimal flocculant/clay ratio. Below this optimal ratio, a proportion of Class 1 particles will remain unflocculated.

Special case where m 1 ≪ m 2
In the case where most of Class 1 particles have largely aggregated and Class 2 particles are dominant (m 1 ≪ m 2 ), it can be assumed that k m ≃ 0 (as V 1 ≃ 0). This does not imply that V LISST, % is constant as Class 2 particles can restructure and hence V LISST, % can be changing in time. This case will be discussed in more details in the next section.

Implementation in a sediment transport model
In the previous subsection, a mass transfer equation has been set-up that is valid in a closed volume V tot . This condition has now to be relaxed, as describing the transport of sediment is based on the transfer of mineral sediment mass from one region of space to another. This is done through the use of advection-diffusion equations. For simplicity, we present here a system whereby the water and fluid movement only occur in the vertical direction z, and where no horizontal advection is considered. The transport equations read: where D z is the eddy diffusivity (not to be confused with the class diameters D k , k = 1,2,3). Summing these two equations it is evident that no sediment mass is lost through the sink/source term k m . Note that in contrast to most models, the settling velocity of Class 2 particles ω s, 2 is a function of time and space as it is depending (through D k and ρ 2 ) on time The first term at the right-hand-side of Eq. (24) represents mass deposition flux S k of particles and the second term the mass erosion flux E k . The following equations are set-up for within the bed: The rate function k m, bed , which expressed the amount of Class 1 particles in the bed that becomes part of Class 2, is not necessarily positive: if Class 1 particles interact with organic matter on/in the bed, they will be part of Class 2 upon resuspension (k m, bed > 0) but if by decay of organic matter Class 1 particles are created then k m, bed < 0.
The settling flux can be expressed as whereby the bottom shear stress is defined as τ and the critical shear stress τ c, k is the stress at which particles of class k start to be eroded.
Note that when k m, bed ∕ = 0 we have m k, bed ∕ = m k (z = 0) and that hence a balance equation should be kept for m k, bed . The empirical exponent n is close to 1.5 for a sandy seabed and close to 1 for a fine sediment bed. For highly organic particles it is expected that τ c, k ≃ 0 as these particles are easily picked-up (in which case, Eq. (27) should be adapted in order to prevent to divide by zero). The curve displaying the dependence of τ c, k on the particle Reynolds number Re p = u k D k ρ w /η where u k is a characteristic particle velocity and η the water viscosity is known as the Shields curve. The dependence of τ c, k on Re p has been observed to obey the Shields curve for a large number of granular (non-cohesive) beds, irrespective of the size, shape, roughness of the grains and flow conditions (from laminar to turbulent). The robustness of the Shields curve has been explained by molecular dynamics simulations of granular beds, where the grain-grain interactions have been varied (Clark et al., 2017). The parameter C k (s − 1 ) depends on the contact forces between particles.
In the extreme case of a very cohesive bed, one can imagine that in good approximation C k ≃ 0 and/or τ/τ c, k ≪ 1. We also refer to (Mehta et al., 2014) for a discussion about the Krone deposition equation in the context of floc aggregation and (Malarkey et al., 2015), (Tolhurst et al., 2002) for the role of biological cohesion in subaqueous beds.

Model parametrization and calibration
In section relations were given to determine the mass transfer rates between Class 1 (mineral sediment) and Class 2 (flocs). In the first subsection, we show how these rates can be determined using laboratory experiments in closed vessels. In a second subsection we show how the model can be parameterized and calibrated using in-situ data.

Materials and methods
To illustrate the model parameterization, we present the results obtained while studying the flocculation of a clay suspension with a cationic polyelectrolyte referenced Zetag 7587 (BASF). The flocculant is composed of a copolymer of acrylamide and quaternary cationic monomer usually used for the conditioning of municipal and industrial substrates. The clay was purchased from the company VE-KA (The Netherlands) and referenced K-10000 and dispersed in demi-water with a conductivity of less than 0.005 mS/cm. The composition of the clay is predominantly quartz, calcite, anorthite and muscovite (Ibanez Sanz, 2018).
The particle size distribution (PSD) as function of time of this suspension was measured by static light scattering using a Malvern Mas-terSizer 2000. The Malvern Mastersiser not only gives PSD's but also relative concentrations. The dry clay particles were suspended in a 1 L jar and each suspension was stirred using a rectangular paddle (25 mm high, 75 mm wide) at 75 RPM, which corresponds to a shear rate of about 50 s − 1 (Logan, 2012). The suspension was pumped from the jar, through the MasterSizer with the help of a peristaltic pump at 20 RPM (± 150 s − 1 ) and then back into the jar. The paddle and pump speeds were chosen as low as possible to minimize the shear, but sufficient enough to prevent any settling of particles in the jar or the tubes. A full PSD could be recorded every 30 s. The PSD of the clay sample is given as the PSD at t = 0 s. The concentration of clay was chosen to be 0.7 g/L, which was the highest concentration at which the laser obscuration was in an acceptable range. At a time defined as t = 1 s flocculant was added to the clay suspension.

Results
The particle size distributions and time evolution of characteristic diameters are given in Fig. 2 A more detailed analysis of the size evolution as function of time in given in Appendix A, where the logistic growth model presented in Section 2 is used. It is shown that the model can be used to fit a very wide family of size evolutions, and that characteristic times for birth and decay can be obtained. The logistic growth model can be used to study more into detail the evolution of the model parameters as function of key variables such as salinity, shear stress and organic matter type and concentration.
The concentration reported by the light scattering device corresponds to the variable defined as V LISST, % , see Eq. (21). Strictly speaking the variable should have been named here V Mastersizer, % as LISST refers to the equipment used in-situ. Both devices are operating using the same principles but with different resolutions. At the beginning of the experiment M = 0 0 .7 g clay is added to V tot = 1L of water. The concentration estimated by the device can be estimated using the following relation: whereby m = M/V tot . The density of mineral clay is ρ 1 = ρ s = 2600 g/L.
This gives a value of V LISST, % = 0.0269% that compares quite well with the initial concentration reported by the MasterSizer of 0.026%. The transfer of mineral sediment mass between Class 1 and Class 2 and the associate rate of change can be estimated using The results are displayed in Fig. 3. It is clear that after that the polymer has been added to the suspension, an error has been introduced as some polymer is detected (leading to an increase in m 1 at t = 30 s). It has therefore been decided to start the evaluation of k m at 60 s. A rapid estimation gives that 0.7 g of clay should have been transferred to Class 2 in a time 0.7/k m (60 s) = 280 s, which is in agreement with the observed PSD evolution in Fig. 2.  The density of Class 2 flocs can be estimated using Eq. (35) (see the corresponding subsection for an explicit derivation): Following the practice adopted to estimate the density of flocs insitu, that V 2, % ≃ 1 (all volume detected by LISST belongs to Class 2 particles) and that hence m = m 2 the following estimate for density of Class 2 particles can be made: The results are shown in Fig. 4. The density calculated according to Eq. (31) is decreasing as function of time. This is easy to understand, as this model assumes that the constant mineral sediment mass m is, from start, contained in flocs. As V LISST, % increases in time (flocs are growing), this leads to a decrease in ρ 2, in− situ . Eq. (30) accounts for the fact that not all mineral sediment is contained in the flocs and that m 2 is increasing in time faster than the product (V LISST, % V 2, % ). This leads to an increase in ρ 2 = ρ 2, lab as a function of time. At the end of flocculation not all mineral sediment is part of Class 2, see Fig. 4, right panel, and this explains the difference between ρ 2 and ρ 2, in− situ .
The settling velocity of the different classes can be estimated using Stokes' relation (Stokes, 1851) For Class 1 this gives, using D 1, mean = 5.5 μm, ω s, 1 = 0.025 mm/s.
The settling velocity of Class 2 is given in Fig. 5. Different expressions for D 2 and ρ 2 have been used to evaluate ω s, 2 . The settling velocity ω s, 2 has been evaluated using D 2 = D50 and ρ 2 is evaluated using Eq. (30), whereas ω s, 2 (in situ) has been evaluated using D 2 = D50 and ρ 2 is evaluated using Eq. (31). The difference in densities results in a difference of about 8 mm/s in settling velocity. The settling velocity ω s, 2 (D 2 ) is evaluated using D 2 = D 2, mean and ρ 2 is evaluated using Eq. (30). The difference in size between ω s, 2 and ω s, 2 (D 2 ) results in a difference of about 6 mm/s. These differences in settling velocity illustrate two major points of concern in the frame of sediment modelling: (a) instrumentation and protocols should be available to correctly estimate size and density of particles, as otherwise this results in large errors in settling velocities and (b) it is important to know the relative amount of mineral sediment in each class. Indeed, if Class 3, as defined at the beginning of this section to be the class of particles of size >200 μm would be composed in a large part of non-mineral components, D 2, mean would be the characteristic size of Class 2 and not D50, which would also induce an error in estimating the settling velocities. In the present case, as discussed above, D50 is probably a good estimate for D2 since Class 3 is rapidly becoming part of Class 2.

Parameters for Class 1
The size of mineral sediment particles is, due to the nature of mineral sediment particles, constant and independent of environmental conditions. Their settling velocity is depending on the water density and viscosity, which is a function of temperature and salinity. Based on modelling practices (Lesser et al., 2000), (Lesser et al., 2004), it is reasonable to assume for this class a settling velocity of ω s, 1 = 0.25 − 0.5 mm/s, which corresponds to the settling velocities of solid mineral particles of density 2600 kg/m3 of size about 10 μm. In in-situ conditions, Class 1 particles (especially the clay fraction) are mainly present in the water column during storm events or dredging activities, when the bed is eroded. Because small organic particles can also have a clay and  silt size it could be difficult to attribute the relative proportion of microorganisms and mineral particles in this size range from LISST data, even though it is expected that during resuspension the majority of clay particles are mineral particles. The mass concentration of these particles can be assessed combining turbidity measurements and estimating the relative ratio between organic and inorganic particles from laboratory experiments on in-situ water samples. After removal of the organic matter, the mass of Class 1 particles can be assessed. This method is more reliable than using LISST data in this small particle size range, as it is known that the "tails" observed in PSD's found by LISST at the smallest and largest sizes are caused by diffractions that cannot be correlated to a proper size (Sequoia Scientific, 2011).

Parameters for Class 2
Class 2 particles have both a dynamic size and a dynamic settling velocity. The size of Class 2 particles can easily be estimated from the D50 measured by LISST, as this Class (in most situations) is containing the largest amount of particles. As stated in the previous subsection, this does not hold in cases where a large inflow of Class 1 particles is observed (close to the bed during storms or dredging activities). Outside these situations, LISST data can be used to calibrate the model and checking that D 2 = D 50 .
The settling velocity of Class 2 particles is dynamic, as flocs can grow, coil or break in space and time. The settling velocity of class 2 particles can be estimated using Eq. (32). A very simple estimation of the floc density in-situ can be made as follows. If we assume that the density of organic matter is very close to water, and that each floc contains a large volume of water, the density ρ floc of a floc is given by where V tot, % = V tot /(V 1 + V 2 + V 3 ) is the inverse of the total volume concentration V LISST, % measured by LISST (total volume of particles detected divided by the considered volume). Most authors assume that m 1 ≪ m 2 and therefore m = m 2 and that V 2, % ≃ 1 (all volume detected by LISST belongs to Class 2 particles). The total mass of mineral sediment per unit volume can then be obtained by OBS (calibrated through filtration or drying experiments in the laboratory) and hence m = m OBS . These assumptions leads to This method to estimate density has been adopted by several authors (Fettweis, 2008), (Many et al., 2016). By comparing the estimation of ρ floc to the predicted ρ floc obtained from laboratory studies, it is possible to confirm that the laboratory experiments have indeed accounted for the relevant variables when studying the dependence of ρ floc on these variables.
The mass concentration of Class 2 particles can then also be calibrated using the simple relation m 2 = m = m OBS ..

Microflocs and macroflocs within Class 2
The formulation for settling velocity for example can greatly benefit from different existing studies. An overview of different expressions for the settling velocities of flocs can be found in (Soulsby et al., 2013), (Shen et al., 2018). In (Soulsby et al., 2013), new formulations for the settling velocity and mass settling flux (the product of settling velocity and sediment concentration) are presented. This study is an extension of the work presented in (Manning and Dyer, 2007) and also provides an overview of alternative expressions used by other authors. The formulations have been tested against independent measurements made in two Northern European estuaries.
The formulation adopted by Soulsby et al. (Soulsby et al., 2013) is very close to the pragmatic approach used in the present article. The authors have defined, based on experimental evidence, two classes of particles: "microflocs" (in the size range  μm) and "macroflocs" (of size >160 μm). The lower boundary of 20 μm corresponds to the detection limit of the equipment. In the frame of the article presented here, these micro and macroflocs are part of Class 2 particles. The authors define the mean diameter of steady-state (equilibrium) micro-and macroflocs as The parameters α μ , α M and k are adjustable and dimensionless. Note that we have used m 2 as the representative mass concentration for the sake of model but that the authors define it as SPM concentration (i.e. m in our notation). This difference is of no consequence when m 1 ≪ m 2 = m. The Kolmogorov lengthscale is defined by where υ (m 2 /s) the kinematic viscosity defined by υ = η/ρ w where η is the water viscosity and ε is the energy dissipation rate per unit of mass of water. The viscous shear rate within the small eddies, G (s − 1 ), is related to the energy dissipation rate by: From which it follows that G is can be evaluated from the root-mean-square of the gradient in the turbulent velocity fluctuations. These are also related to the turbulent shear stress τ (Pa) through the shear velocity u * by where κ is the Karman constant (taken to be 0.40), and where h is the water depth and z is the height above the bed. Near the bed, (z ≪ h) one gets Z ≃ z. The found dependence of floc size on shear rate has been verified by multiple studies (Jarvis et al., 2005), (Manning and Dyer, 2007), (Mietta et al., 2009a(Mietta et al., , 2009b, . However it was also found that this relation mainly hold for salt-induced flocs. Polyelectrolyte-induced flocs are usually larger than the Kolmogorov microscale (Ibanez Sanz, 2018) and the average aggregate size can then be fitted using a relation whereby D macro, eq = A/(1 + BG γ ) where A and B are adjustable parameters. The parameter γ is found to be in the range [1-2], which is in agreement with the results of Bubakova et al. (Bubakova et al., 2013) who analyzed the shear rate dependency at steady-state of aggregates found in a water reservoir. The dependence of D macro, eq on shear might therefore require further investigation. For the density, the authors propose the following equations, to be compared with Eq. (16): In these equations, the densities are in kg/m 3 and the diameters in μm. These equations have been obtained by analyzing the data of the Tamar, Gironde and Dollard estuaries (Soulsby et al., 2013). Interestingly, the authors have found that neither the densities nor the relative ratio of micro and macroflocs depend on shear. As they also discuss, this suggests that the flocculation kinetics are usually fast and the system always near steady-state. This holds of course mainly for systems where no rapid and major changes are occurring. During dredging activities, for example, it is much more likely that flocculation kinetics are to be considered, and hence Class 1 particles should be accounted for. We here (as the authors do) will assume that m 1 = 0 and m = m 2 . Defining m macro as the mass mineral sediment within macroflocs, the authors found that m macro /m verifies whereby m should be in mg/L. This equation holds for 1 < m < 1174 mg/L. By extrapolation, the authors proposed that m macro /m = 0.1 for m < 1 mg/L and m macro /m = 1 for m > 1174 mg/L. Another equation was proposed in (Manning and Dyer, 2007), but one can verify that both equations are within the same order of magnitude for the whole range of SPM investigated.
In terms of the model detailed in Section 2, it follows that Class 2 particles are now split in two, according to m = m micro + m macro . The mass transfer equation to consider is given by The assumption is made that the rate constant is only depending on the total mass concentration m(t) which is changing in time by import or export of sediment from the considered volume. One finds that with dr/dm = 0.221/(m × ln (10)). The mass concentration m macro is increasing when m is increasing and decreasing when m is decreasing.
The assumption used relies on the fact that the flocs are mainly locally produced, and that the kinetics of reaction is rapid (in fact here we assume that the (de)flocculation is instantaneous with the change in m).
One could easily relax this assumption by considering that the steadystate (equilibrium) mass concentrations are given by m macro, eq = r × m, in which case we get whereby an estimation for the rate dm macro /dm macro, eq should be given and could a-priori be a function of time (among other parameters such as m OM , shear, salinity…).
The floc settling velocities ω s, macro and ω s, micro are evaluated starting from Stokes equation, but also accounts for the fact that observations have shown that the settling velocities of flocs were increasing with shear for small shear, before reaching a maximum and then decreasing with large shear. This behavior as function of shear at low shear is to be related to the discussion given at the beginning of this subsection, where we question the validity of Eq. (37) for the dependence of size on shear.
As ω s, micro/macro is only decreasing with shear when Stokes is used, the Stokes settling formula was multiplied by a correcting exponential term that accounts for the initial increase with shear: The parameters are g = 9.81 m/s 2 , u *, macro = 0.067 m/s, u *, micro = 0.025 m/s, d macro = 10 − 4 m and d micro = 10 − 5 m. We refer to (Soulsby et al., 2013) for a further discussion about how these formulations compare with in-situ data and how they compare with expressions proposed by other authors. Note that, as proposed in Section 3, the settling velocities are time-dependent.
As Stokes settling velocity, evaluated using Eqs. (37), (43) for the diameters and densities, does not account for the increase in settling velocities at small shears, it has been argued (Soulsby et al., 2013), (Winterwerp et al., 2006) that this might be due to the limited residence time during which flocs can form. If this is the case, flocculation dynamics should be accounted forwhich can be done following the procedure described in Section 3.
If, on the other hand, we follow the hypothesis used in the present section that flocculation kinetics are very fast another explanation might hold, if we consider flocs as being produced by organic matter. At small shear the flocs usually have a very open structure and the organic matter is loosely bound. This creates large flocs, with dangling ends, which are responsible for the fluctuations in particle size for large particles sizes, see Fig. 2. When these ends collapse on the flocs they hereby reduce slightly the mean floc size. At higher shear, the flocs, without a significant change in mean size, can have a different structure, whereby less water is present in the floc structure and replaced by densely coiled organic substances. The mean density of flocs is then not much different from its mean density at lower shear as the density of organic matter is close to the density of water, but this change in structure can result in a significant increase in settling velocity as function of shear. At higher shears, the steady-state size of flocs becomes a function of shear, as has been discussed above. This results in a decrease in settling velocity with increasing shear.

Conclusion
The main aim of the article was to propose a new flocculation model to be incorporated in a large-scale sediment transport model. In the first part of the article, we present a logistic growth model to follow the changes in particle sizes distributions as measured by laser diffraction techniques in the laboratory or in-situ. The time scales of the different processes (growth of flocs, their decay) can easily be found using the proposed model. This model, in contrast to the standard population balance equation with aggregation and break-up parameters, gives the possibility to have uncorrelated birth and decay rates. This feature is shown to be of importance when flocculation is induced by the presence of organic matter. Upon a reduction in volume of a floc (triggered for instance by the coiling of organic matter due to a change in shear), the ω s,macro = 0.095g mass of mineral sediment within the floc remains constant. This phenomenon cannot properly be accounted for using the classical PBE model as the model attributes a given mineral sediment mass to a floc of a given size. The flocculation model we propose to be incorporate in a sediment transport model is based on a simple mass transfer between two classes of particles: one class (Class 1) consisting of mineral sediment particles and one class (Class 2) consisting of mineral sediment flocculated with organic matter. Class 1 is defined by a total mineral sediment mass, associated with a constant settling velocity. Class 2 is also defined by a total mineral sediment mass, associated with a dynamic settling velocity. We discuss the time evolution of size and density of Class 2 particles. In the model description we attribute a single size (D50) and settling velocity to Class 2. We then proceed to show with an example that the model can, without problem, account for a subdivision of Class 2 in microflocs and macroflocs as has been proposed by other authors. Further studies should demonstrate the necessity to subdivide (or not) the different classes in sub-classes for the purpose of sediment transport modelling. The methods to assess the masses and settling velocities in each classes is discussed. These are needed to estimate the rates of mass, density and size changes, and validate the model. In particular, we show that only using LISST and OBS data to estimate mineral sediment mass and settling velocities can lead to significant errors when the system is not at steady-state. In order to better assess the settling velocity of flocs, additional measurements should be performed, such as video microscopy of flocs in a settling column at rest. The proportion of mineral sediment and organic content of flocs should also be better studied by laboratory analysis so that in the future a correlation can be made between floc size and mineral sediment content. At present, it is usually assumed that the mass of mineral sediment in a floc follows a fractal growth, see Eq. (16), which in theory is only valid for salt-induced flocculation in the absence of organic matter.
A third class (Class 3) is defined as being the class containing organic matter. This class cannot be accounted for in a mass balance set of equations, as the mass of organic matter cannot properly be defined due to large variation in organic matter type and quantity (from small polyelectrolyte to living microorganisms that are able to produce poly-electrolytes…). The largest particles from Class 3 are however visible in in-situ PSD data, and this is why this class is defined.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table 1
Fit values for the characteristic times.