The Influence of Temperature on the Bulk Settling of Cohesive Sediment in Still Water with the Lattice Boltzmann Method

Flocculation is very common and significant for cohesive sediment in coastal areas, and the influence of temperature on it cannot be neglected. The Lattice Boltzmann Method (LBM), combined with the extended Derjaguin-Landau-Verwey-Overbeek (XDLVO) theory, which considers the micro-interaction forces between particles, was applied to simulate the settling and flocculation processes of cohesive sediment under various temperature conditions. The floc size, floc volume, suspended sediment concentration (SSC), and settling velocities were analyzed. The analyses revealed that with increasing temperature, both the mean floc diameter and floc volume grew, while the maximum floc diameter initially increased and then slightly decreased with its peak at 10 ◦C. During settling, the SSC change rate was exponentially related to the SSC, with an optimal fitting index of 0.3. The LBM sediment settling velocity was also compared with some formulas and physical model tests; the comparison results consistently demonstrated that the LBM was reasonable for modeling the bulk settling of cohesive sediment. Further discussions illustrated that the cohesive sediment is more difficult to flocculate at low temperatures due to the low aggregation frequency, while at high temperatures, some large flocs broke easily due to the effect of the short-distance force and macro force.


Introduction
Cohesive sediment is an important type of sediment in the near shore zone.As a result of its own physical properties and external environments, cohesive sediment can form flocs easily.Temperature is one of many factors that influence flocculation [1][2][3].Particularly in middle-and high-latitude areas, the water temperature varies drastically: it may be close to 0 • C in winter and exceed 30 • C in summer [4][5][6][7].
Thus far, much research has been conducted to study the influence of temperature on cohesive sediment.Lau [8] studied temperature's effects on the settling velocity and deposition of cohesive sediment in an annular channel of distilled water that was housed in a temperature-controlled chamber.He claimed that, as temperature increased, the repulsive forces between particles increased bulk settling and flocculation process of cohesive sediment under different temperature conditions.The simulation results were analyzed to reveal the influence of temperature and its mechanism.
The paper is organized as follows: the numerical approach is described in Section 2; the computational conditions are obtained in Section 3; in Section 4, the effects of the temperature on the flocculation and settling of illite are indicated; a discussion and the shortcomings of this research are included in Section 5; and we present our conclusions in Section 6.

Lattice Boltzmann Method
The LBM is a relatively new numerical technique for modeling physical system responses.The LBM originated from Lattice Gas Automata.The basic concept of the LBM is to represent fluid as a particle distribution function located at each lattice node.Fluid particles move to neighboring nodes at discrete time steps, colliding with other fluid particles.In the LBM approximation, the fluid is described by a density distribution function f i (x,t), which describes the particle status at a lattice location x at time t with the discrete velocity e i .The Boltzmann equation is used to solve the collision-induced evolution of the fluid particle, and the equation can be written as: where the subscript i represents the directions in which the particle may move.The D3Q19 topology [30] is used in this study, which is a three-dimensional cubic lattice with 19 velocity directions, as illustrated in Figure 1.Ω i (f i ) is the collision operator; in Nguyen and Ladd's model [30], Ω i (f i ) can be written as follows: where f j eq is the local equilibrium function, and f j neq is the non-equilibrium function, i.e., f neq j = f j − f eq j .The hydrodynamic parameters, such as the mass density ρ and momentum ρu, are the function of the distribution function f and discrete velocity e i , which can be written as follows: The density and momentum should satisfy the mass conservation and momentum conservation.A suitable equilibrium distribution form of the D3Q19 topology can be written as follows [30]: where c s = c 2 l /3 is the speed of sound; c l = ∆x/∆t is the lattice speed, in which ∆x is the lattice space size and ∆t is the time step; and |e i | equals 0 (i = 0), 1(i = 1, 2, . . ., 6) and √ 2(i = 7, 8, . . ., 18).The coefficients of the three speeds a e i are 0 (i = 0), 1/18 (i = 1, 2, . . ., 6) and 1/36 (i = 7, 8, . . ., 18).

Figure 1.
The possible velocity directions in the D3Q19 topology.The speed of discrete velocity |ei| equals 0 when particle keeps its original position after a time step, equals 1 when it moves to the faces of the cubic and 2 to the edges of the cubic.lij are the matrix elements of the linearized collision operator, which must satisfy the following eigenvalue equations: where i i e e is the traceless part of i i e e .The first two equations result from the conservation of mass and momentum, and the last two equations describe the isotropic relaxation of the stress tensor.The eigenvalues λ and λv are related to the shear viscosity η and bulk viscosity ηv, respectively, which are within the range of −2 < λ < 0, in which ( ) (1 )( : ) 3 is the non-equilibrium second moment, in which Considering an externally imposed force density F, the time evolution of the LBM includes an additional contribution ( , The LBM is well-suited for the specific problem of modeling solid particle suspensions because of its ability to solve particles movement with arbitrary shapes and complex geometries [30][31][32].

The Extended Derjaguin-Landau-Verwey-Overbeek Theory
The extended Derjaguin-Landau-Verwey-Overbeek (XDLVO) theory on the interactions of particles can be used to explain the flocculation of cohesive sediment particles [33,34].According to this theory, there are three short-distance forces acting between particles in aqueous environments (Figure 2), i.e., the Lifshitz-van der Waals attractive force F LW i-j, the electrostatic double-layer repulsive force F EL i-j, and the Lewis acid-base force F AB i-j.Each of them is the negative derivative of the corresponding potential to the distance.Among them, temperature only affects the electrostatic double-layer repulsive force by changing the Debye length κ −1 .For more details, please refer to [20,26,[31][32][33][34][35].l ij are the matrix elements of the linearized collision operator, which must satisfy the following eigenvalue equations: where e i e i is the traceless part of e i e i .The first two equations result from the conservation of mass and momentum, and the last two equations describe the isotropic relaxation of the stress tensor.
The eigenvalues λ and λ v are related to the shear viscosity η and bulk viscosity η v , respectively, which are within the range of −2 < λ < 0, in which η = −ρc 2 s (1/λ + 1/2) and η ν = −ρc 2 s [2/(3λ ν ) + 1/3].A three-parameter collision operator is used in the present study.The post-collision distribution f i * can be written in the same form as f i eq : where neq, * = (1 + λ) neq + 1 3 (1 + λ v )( neq : I)I and neq = − eq is the non-equilibrium second moment, in which eq = i e i e i f eq i = ρc 2 s I + ρuu.Considering an externally imposed force density F, the time evolution of the LBM includes an additional contribution F i (x, t): The LBM is well-suited for the specific problem of modeling solid particle suspensions because of its ability to solve particles movement with arbitrary shapes and complex geometries [30][31][32].

The Extended Derjaguin-Landau-Verwey-Overbeek Theory
The extended Derjaguin-Landau-Verwey-Overbeek (XDLVO) theory on the interactions of particles can be used to explain the flocculation of cohesive sediment particles [33,34].According to this theory, there are three short-distance forces acting between particles in aqueous environments (Figure 2), i.e., the Lifshitz-van der Waals attractive force F LW i-j , the electrostatic double-layer repulsive force F EL i-j , and the Lewis acid-base force F AB i-j .Each of them is the negative derivative of the corresponding potential to the distance.Among them, temperature only affects the electrostatic double-layer repulsive force by changing the Debye length κ −1 .For more details, please refer to [20,26,[31][32][33][34][35].The electrostatic double-layer repulsive force between spherical particles with radii of R1 and R2 can be written as follows [20]: where hij = rij − (R1 + R2) is the net distance between spheres, and rij is the distance between sphere centers; ε0 is the dielectric permittivity in a vacuum with a value of 8.854 × 10 −12 C 2 /(J•m); εr is the relative dielectric constant, and for water, it is 78.5; ψ0 is the surface potential of particles, which is very sensitive to salinity and pH but not sensitive to temperature when the temperature is below approximately 150 °C [36]; κ −1 is the Debye length, related to temperature, cation valence, and ion concentration, which can be written as follows: where 0 e is the element charge 1.6 × 10 −19 C; A N is Avogadro's number 6.022 × 10 23 /mol; k is the Boltzmann constant 1.38 × 10 −23 J/K; T is the absolute temperature with a unit of K; c is the cation concentration with a unit of mol/L; and z represents the cation valence and is dimensionless.

Criterion Distance of Flocculation
Yang et al. [37] mentioned that the aggregation sign of cohesive sediment should be the contact of the sliding surface.They suggested that the anions are distributed on the surface of particles, and adsorb the cations in the media.The zone with cations can be separated into two layers, the inner layer with a high density of cations, which is called the absorbed layer, and the outer layer with less cations, which is called the diffuse layer, as illustrated in Figure 3. Yang et al. [37] consider that when the surface distance is less than twice the thickness of the slipping layer, the particles are wrapped by a common slipping layer, and they form a new or enlarge an old floc.The electric potential decays exponentially at large distance with decay length given by the Debye length κ −1 , with ψ0 on the particle surface and ζ at the distance δ.Both ψ0 and ζ can be measured in experiments, thus the slipping layer thickness δ can be calculated as follows: The electrostatic double-layer repulsive force between spherical particles with radii of R 1 and R 2 can be written as follows [20]: where is the net distance between spheres, and r ij is the distance between sphere centers; ε 0 is the dielectric permittivity in a vacuum with a value of 8.854 × 10 −12 C 2 /(J•m); ε r is the relative dielectric constant, and for water, it is 78.5; ψ 0 is the surface potential of particles, which is very sensitive to salinity and pH but not sensitive to temperature when the temperature is below approximately 150 • C [36]; κ −1 is the Debye length, related to temperature, cation valence, and ion concentration, which can be written as follows: where e 0 is the element charge 1.6 × 10 −19 C; N A is Avogadro's number 6.022 × 10 23 /mol; k is the Boltzmann constant 1.38 × 10 −23 J/K; T is the absolute temperature with a unit of K; c is the cation concentration with a unit of mol/L; and z represents the cation valence and is dimensionless.

Criterion Distance of Flocculation
Yang et al. [37] mentioned that the aggregation sign of cohesive sediment should be the contact of the sliding surface.They suggested that the anions are distributed on the surface of particles, and adsorb the cations in the media.The zone with cations can be separated into two layers, the inner layer with a high density of cations, which is called the absorbed layer, and the outer layer with less cations, which is called the diffuse layer, as illustrated in Figure 3. Yang et al. [37] consider that when the surface distance is less than twice the thickness of the slipping layer, the particles are wrapped by a common slipping layer, and they form a new or enlarge an old floc.The electric potential decays exponentially at large distance with decay length given by the Debye length κ −1 , with ψ 0 on the particle surface and ζ at the distance δ.Both ψ 0 and ζ can be measured in experiments, thus the slipping layer thickness δ can be calculated as follows: Water 2019, 11, 945 where δ is the slipping layer thickness when the potential ψ equals the zeta-potential ζ.From Equations ( 9) and ( 10), a rising temperature will thicken δ.Thus, twice the value of δ was taken as the criterion distance for flocculation.
Water 2019, 11, 945 6 of 16 where δ is the slipping layer thickness when the potential ψ equals the zeta-potential ζ.From Equations ( 9) and ( 10), a rising temperature will thicken δ.Thus, twice the value of δ was taken as the criterion distance for flocculation.

Computational Conditions
Four cases were set up to study the influence of temperature on the bulk settlement of cohesive sediment.The temperatures were 5 °C, 10 °C, 20 °C, and 30 °C in the four cases.The domain calculation area was 0.288 mm × 6 mm × 0.240 mm, and the space precision was 2 μm, thus a total of 144 lu × 3000 lu × 120 lu (lattice unit).The surroundings were defined as periodic boundaries, and the upper and lower boundaries were set as solid no-slip boundaries.The initial sediment concentration was set to 1.3 kg/m 3 (volume concentration of 0.049%).Eight hundred sediment particles with diameters from 5 μm to 10 μm were scatted randomly in the whole calculation area.The calculation time step was set as 10 −6 s, and the total time was 20 s.In all cases, the cation was selected as Na + , with a valence of +1 and a concentration of 0.085 mol/L or salinity of 5 ppt.Illite was considered to be the clay mineral with a surface potential of −27.22 mV [38] as illite is widely distributed and its chemical and physical properties that needed in this study are easy to obtain as they are studied more extensively.The sediment density was ρs = 2650 kg/m 3 .The water density was ρ = 1000 kg/m 3 , without considering the influence of salinity and temperature as the error of the density at the values adopted will be less than 4% according to the sea water state equation proposed by UNESCO [39].The other parameters are listed in Table 1.

Floc Size and Floc Volume
Floc size and volume are two of the floc properties that are easy to obtain from the numerical results.Figure 4a shows the time series of the maximum and mean floc sizes.The final results of them  [37].This zone is comprised of an absorbed layer and a diffuse layer.The slipping layer thickness can be calculated by the potentials ψ 0 on the particle surface and ζ at the distance δ as the potential decays exponentially.When the particles are close enough and wrapped by a common slipping layer, they are considered as forming a new or enlarging an old floc.

Computational Conditions
Four cases were set up to study the influence of temperature on the bulk settlement of cohesive sediment.The temperatures were 5 • C, 10 • C, 20 • C, and 30 • C in the four cases.The domain calculation area was 0.288 mm × 6 mm × 0.240 mm, and the space precision was 2 µm, thus a total of 144 lu × 3000 lu × 120 lu (lattice unit).The surroundings were defined as periodic boundaries, and the upper and lower boundaries were set as solid no-slip boundaries.The initial sediment concentration was set to 1.3 kg/m 3 (volume concentration of 0.049%).Eight hundred sediment particles with diameters from 5 µm to 10 µm were scatted randomly in the whole calculation area.The calculation time step was set as 10 −6 s, and the total time was 20 s.In all cases, the cation was selected as Na + , with a valence of +1 and a concentration of 0.085 mol/L or salinity of 5 ppt.Illite was considered to be the clay mineral with a surface potential of −27.22 mV [38] as illite is widely distributed and its chemical and physical properties that needed in this study are easy to obtain as they are studied more extensively.The sediment density was ρ s = 2650 kg/m 3 .The water density was ρ = 1000 kg/m 3 , without considering the influence of salinity and temperature as the error of the density at the values adopted will be less than 4% according to the sea water state equation proposed by UNESCO [39].The other parameters are listed in Table 1.

Floc Size and Floc Volume
Floc size and volume are two of the floc properties that are easy to obtain from the numerical results.Figure 4a   At the end of the simulation, in the cases of 5 °C, 10 °C, 20 °C, and 30 °C, the maximum floc sizes were 20.0 μm, 24.0 μm, 22.3 μm, and 23.3 μm, respectively, with the peak value at 10 °C.In the cases of 20 °C and 30 °C, the maximum floc sizes fluctuated slightly in the latter phase (Figure 4a).
The mean floc size increased slightly but more stably than the maximum size did, from 15.2 μm at 5 °C to 15.5 μm at 10 °C, 15.8 μm at 20 °C, and 15.9 μm at 30 °C.The floc volume increased more obviously than the floc size.They were 3.8% at 5 °C, 4.0% at 10 °C, 4.3% at 20 °C, and 5.5% at 30 °C (Figure 4b).

Settling and Flocculation Process
The microscopic process of settling and flocculation might explain the preceding phenomenon.This process can be easily given out via the numerical simulations.Three particles numbered #681, #689, and #721 were selected to illustrate the settling process.In all four cases, they had formed or almost formed a floc.
The mean floc size increased slightly but more stably than the maximum size did, from 15.2 µm at 5 • C to 15.5 µm at 10 • C, 15.8 µm at 20 • C, and 15.9 µm at 30 • C. The floc volume increased more obviously than the floc size.They were 3.8% at 5 • C, 4.0% at 10 • C, 4.3% at 20 • C, and 5.5% at 30 • C (Figure 4b).

Settling and Flocculation Process
The microscopic process of settling and flocculation might explain the preceding phenomenon.This process can be easily given out via the numerical simulations.Three particles numbered #681, #689, and #721 were selected to illustrate the settling process.In all four cases, they had formed or almost formed a floc.
The corresponding diameters for particles #681, #689, and #721 were 6.00 µm, 7.16 µm, and 9.36 µm, respectively.The distance between the primary particles and their velocities at different temperatures are plotted in Figure 5.
This process can be easily given out via the numerical simulations.Three particles numbered #681, #689, and #721 were selected to illustrate the settling process.In all four cases, they had formed or almost formed a floc.
The corresponding diameters for particles #681, #689, and #721 were 6.00 μm, 7.16 μm, and 9.36 μm, respectively.The distance between the primary particles and their velocities at different temperatures are plotted in Figure 5.At the beginning, the large particle falling fast was on top of the small one with a low speed, resulting in a diminishing net distance between them before they met, as is common in all cases.At a low temperature (5 °C, Figure 5a), the three primary particles cannot get close enough to form a new floc, and fell individually.The velocities of the particles increased when they collided (12.5 s to 17.5 s), but finally they readjusted to their initial speed with the increasing distance between each other.At a medium temperature (10 °C and 20 °C, Figure 5a,b), the three particles formed a floc and settled together with a common speed higher than that of each primary particle; however, the time when the floc formed was earlier at 20 °C than that at 10 °C.In the case of high temperature (30 °C, Figure 5d), the small particle (#681) was first captured by the medium one (#689), forming an unstable floc, which was destroyed by the collision of the large particle (#721), and then flocculated with the large particle, leaving the medium particle settling individually with a low speed.

Suspended Sediment Concentration
Figure 6 illustrates the variation in SSC for each case every 5 s.It is shown that, in all four cases, SSCs above 4.0 mm continued to decline, while those at the bottom layers were increasing.They remained relatively stable in the 4.0-5.0mm layer, but a small difference can be seen between the cases: the SSC in this layer went up slightly in the 5 °C case, but had a slight reduction in the 30 °C case, and stayed nearly constant in the other two cases.At the beginning, the large particle falling fast was on top of the small one with a low speed, resulting in a diminishing net distance between them before they met, as is common in all cases.At a low temperature (5 • C, Figure 5a), the three primary particles cannot get close enough to form a new floc, and fell individually.The velocities of the particles increased when they collided (12.5 s to 17.5 s), but finally they readjusted to their initial speed with the increasing distance between each other.At a medium temperature (10 • C and 20 • C, Figure 5a,b), the three particles formed a floc and settled together with a common speed higher than that of each primary particle; however, the time when the floc formed was earlier at 20 • C than that at 10 • C. In the case of high temperature (30 • C, Figure 5d), the small particle (#681) was first captured by the medium one (#689), forming an unstable floc, which was destroyed by the collision of the large particle (#721), and then flocculated with the large particle, leaving the medium particle settling individually with a low speed.

Suspended Sediment Concentration
Figure 6 illustrates the variation in SSC for each case every 5 s.It is shown that, in all four cases, SSCs above 4.0 mm continued to decline, while those at the bottom layers were increasing.They remained relatively stable in the 4.0-5.0mm layer, but a small difference can be seen between the cases: the SSC in this layer went up slightly in the 5 • C case, but had a slight reduction in the 30 • C case, and stayed nearly constant in the other two cases.

Suspended Sediment Concentration
Figure 6 illustrates the variation in SSC for each case every 5 s.It is shown that, in all four cases, SSCs above 4.0 mm continued to decline, while those at the bottom layers were increasing.They remained relatively stable in the 4.0-5.0mm layer, but a small difference can be seen between the cases: the SSC in this layer went up slightly in the 5 °C case, but had a slight reduction in the 30 °C case, and stayed nearly constant in the other two cases.Chen and Shao [40] revealed that the SSC change rate in still water was generally fitted to a firstorder equation as follows: where kc is the attenuation coefficient, in which a large value indicates a faster incline.The time series curves of relative SSC (ratio of c(t) at time t to initial c0) of water depths shallower than 4.5 mm, in which layer the SSC is somewhat unchanged as shown in Figure 6, are draw in Figure 7a.As shown in Figure 7a, the fitting results of Equation (11) (the dash lines) were less superior in these cases, therefore, some trials of indexes of n for C in the right side had been done.As a result, 0.3 was the best fitted index (the solid lines in Figure 7a), thus yielding the following formula: The time when the SSC reduced to half of the initial SSC was taken as the half-life settlement period t1/2.From Equation ( 12), the half-life period t1/2 can be written as follows: The kc and half-life period t1/2 can be solved from the fitting curve in Figure 7a and they are shown in Figure 7b.The temperature increase had a positive effect on kc and a negative effect on t1/2, indicating a faster change in SSC at high temperatures than at low temperatures.Chen and Shao [40] revealed that the SSC change rate in still water was generally fitted to a first-order equation as follows: where k c is the attenuation coefficient, in which a large value indicates a faster incline.The time series curves of relative SSC (ratio of c(t) at time t to initial c 0 ) of water depths shallower than 4.5 mm, in which layer the SSC is somewhat unchanged as shown in Figure 6, are draw in Figure 7a.As shown in Figure 7a, the fitting results of Equation ( 11) (the dash lines) were less superior in these cases, therefore, some trials of indexes of n for C in the right side had been done.As a result, 0.3 was the best fitted index (the solid lines in Figure 7a), thus yielding the following formula: The time when the SSC reduced to half of the initial SSC was taken as the half-life settlement period t 1/2 .From Equation ( 12), the half-life period t 1/2 can be written as follows: The k c and half-life period t 1/2 can be solved from the fitting curve in Figure 7a and they are shown in Figure 7b.The temperature increase had a positive effect on k c and a negative effect on t 1/2, indicating a faster change in SSC at high temperatures than at low temperatures.

Sediment Settling Velocity
In numerical simulations, the bulk velocities can be statistically analyzed from microscope by averaging the velocities of individual particles and flocs, which are the directional output of the simulation, with a weight of each volume; however, they are commonly difficult to measure directly in physical experiments.As a result, the bulk velocity is often calculated from the SSC half-life period t1/2 and the settling distance [6].The bulk velocity can be expressed simply as follows: ( ) where ( ) u H is the time-averaged (from 0 to t1/2) bulk settling velocity at water depth H, in which depth the SSC changes slightly during settling, and that is why 4.5 mm (4.0-5.0 mm) is chosen in Figure 7a.The derivation of Equation ( 14) can refer to You's work [41].The water depth H and time t1/2 in Equation ( 14) can be obtained easily through macroscopic observation.
In Table 2, velocity 1 summarized the bulk velocity calculated from Equation ( 14) with the parameters in Figure 7b.The time-averaged statistical velocities during the period of 0-t1/2 in the 4.5 mm layer (4.0-5.0 mm) were listed in Table 2, as velocity 2. The results of the two methods were almost identical, with a small relative error of approximately 2%.Equation ( 14) can be taken as the connection between the microscopic statistical methods and physical test methods.

Sediment Settling Velocity
In numerical simulations, the bulk velocities can be statistically analyzed from microscope by averaging the velocities of individual particles and flocs, which are the directional output of the simulation, with a weight of each volume; however, they are commonly difficult to measure directly in physical experiments.As a result, the bulk velocity is often calculated from the SSC half-life period t 1/2 and the settling distance [6].The bulk velocity can be expressed simply as follows: where u(H) is the time-averaged (from 0 to t 1/2 ) bulk settling velocity at water depth H, in which depth the SSC changes slightly during settling, and that is why 4.5 mm (4.0-5.0 mm) is chosen in Figure 7a.The derivation of Equation ( 14) can refer to You's work [41].The water depth H and time t 1/2 in Equation ( 14) can be obtained easily through macroscopic observation.
In Table 2, velocity 1 summarized the bulk velocity calculated from Equation ( 14) with the parameters in Figure 7b.The time-averaged statistical velocities during the period of 0-t 1/2 in the 4.5 mm layer (4.0-5.0 mm) were listed in Table 2, as velocity 2. The results of the two methods were almost identical, with a small relative error of approximately 2%.Equation ( 14) can be taken as the connection between the microscopic statistical methods and physical test methods.
The floc velocities are also compared with the data of Xia et al. [12], Khelifa and Hill [16], Guo and He [13], Dyer and Manning [42], and Manning et al. [43] and the formulas of Manning et al. [43], Khelifa and Hill [16] and Winterwerp [44] in Figure 8a.The primary particle diameter in the simulations ranged from 5 µm to 10 µm, and the simulated floc settling velocities lie between Winterwerp's [44] lines with primary particle sizes of 1 µm and 20 µm, and most of them between Khelifa and Hill's [13] lines with the same primary sizes.They also are between Manning et al.'s [43] lines with an effective density of 160 kg/m 3 and 1600 kg/m 3 , but closer to the line of 1600 kg/m 3 , probably due to the small size of the flocs, so they have a larger effective density [16,42].The simulation results overlap with the in situ observation data of Xia et al. [12] for July 1999 and January 2000.It should be noted that the simulation velocities are each floc's speed, while that of other studies are bulk velocity, or statistical macro data of suspended samples.
Table 2. Bulk settling velocity of different methods (velocity 1 v 1 is derived from H and t 1/2 ; velocity 2 v 2 is the statistical result of all particles and flocs at water depths from 4.0 mm to 5.0 mm; ABS(1-2) |v 1 −v 2 | means the absolute difference between velocity 1 and velocity 2; the relative difference between velocity 1 and velocity 2).simulations ranged from 5 μm to 10 μm, and the simulated floc settling velocities lie between Winterwerp's [44] lines with primary particle sizes of 1 μm and 20 μm, and most of them between Khelifa and Hill's [13] lines with the same primary sizes.They also are between Manning et al.'s [43] lines with an effective density of 160 kg/m 3 and 1600 kg/m 3 , but closer to the line of 1600 kg/m 3 , probably due to the small size of the flocs, so they have a larger effective density [16,42].The simulation results overlap with the in situ observation data of Xia et al. [12] for July 1999 and January 2000.It should be noted that the simulation velocities are each floc's speed, while that of other studies are bulk velocity, or statistical macro data of suspended samples.Figure 8b illustrates the variation of average floc velocities in different temperature conditions, including the observational data from Xia et al. [12] and Guo and He [43], the experimental data from Wan et al. [6], and the simulation results of LBM.The settling speed differs greatly owing to factors other than temperature, such as SSC, labeled in the legend; however, all the studies reveal a trend of the settling velocity increasing with the water temperature.The velocity change rate with temperature of this numerical result is 0.0382 mm/s/°C-larger than the rate of Wan et al.'s [6] physical experiment results of low SSC (0.0248 mm/s/°C), smaller than that of Wan et al.'s [6] high SSC (0.0513 mm/s/°C) and Xia et al.'s [12] observation results (0.1031 mm/s/°C), and very close to Guo and He's [43] observation data (0.0364 mm/s/°C).

Discussion
The above results indicate that the LBM is reasonable for using in modeling bulk settling velocity and that some collisions led to aggregation, while others did not, or, on the contrary, to disaggregation.On referring to Kim and Stolzenbach's work [45], the ratio of the global collision number to the total particle number is defined as the collision frequency ηglo, the ratio of collisions that created new flocs or enlarged original flocs to the total collisions is the capture frequency ηcap, and the product of the two ratios is the aggregation frequency Eagg = ηglo × ηcap.The three parameters were analyzed and are illustrated in Figure 9.
As shown in Figure 9a, with increasing temperature, ηglo continued to grow up, while ηcap began to decline slightly after its peak at about 20 °C after an initial increase.Sterling et al. [46] stated that the collision frequency was due to three factors: Brownian motion, turbulence shear, and differential settling.In this simulation, the sediment particles were much larger than the molecule and thus the Brownian motion could be ignored; the turbulence effect was omissible in still water; therefore, the main factor was the differential settling.Figure 9b showed a positive correlation between ηglo and  [42], and Manning et al. [43]; theory results of Manning et al. [43], Winterwerp [44], and Khelifa and Hill [43]; and the simulation results of LBM under temperature conditions ranging from 5 • C to 30 • C; (b) floc velocities and temperature from the observational data of Xia et al. [12], Wan et al. [6], Guo and He [13], and this simulation's results.
Figure 8b illustrates the variation of average floc velocities in different temperature conditions, including the observational data from Xia et al. [12] and Guo and He [43], the experimental data from Wan et al. [6], and the simulation results of LBM.The settling speed differs greatly owing to factors other than temperature, such as SSC, labeled in the legend; however, all the studies reveal a trend of the settling velocity increasing with the water temperature.The velocity change rate with temperature of this numerical result is 0.0382 mm/s/ • C-larger than the rate of Wan et al.'s [6] physical experiment results of low SSC (0.0248 mm/s/ • C), smaller than that of Wan et al.'s [6] high SSC (0.0513 mm/s/ • C) and Xia et al.'s [12] observation results (0.1031 mm/s/ • C), and very close to Guo and He's [43] observation data (0.0364 mm/s/ • C).

Discussion
The above results indicate that the LBM is reasonable for using in modeling bulk settling velocity and that some collisions led to aggregation, while others did not, or, on the contrary, to disaggregation.On referring to Kim and Stolzenbach's work [45], the ratio of the global collision number to the total particle number is defined as the collision frequency η glo , the ratio of collisions that created new flocs or enlarged original flocs to the total collisions is the capture frequency η cap , and the product of the two ratios is the aggregation frequency E agg = η glo × η cap .The three parameters were analyzed and are illustrated in Figure 9.
Water 2019, 11, 945 12 of 16 settling velocity, which was consistent with the viewpoint of Sterling et al. [46].ηglo in the simulations ranges from 2.1% to 4.2%, very close to the results of Kim and Stolzenbach [45] (from 2.96% to 3.20%).
The potential reasons for the difference between the two results might lie in Kim and Stolzenbach's [45] presumption that their model neglected the repulsive colloidal interaction, which is very sensitive to temperature (in Section 2.2) but was not mentioned in their study, so they might not take the temperature into account in their work.
(a) (b) When the distance between two particles was less than 25 nm, the three short-distance forces were the most important.Among the forces, only the electrostatic double-layer repulsion force is related to the temperature.From Figure 10a, the repulsive force was higher at a high temperature or a short distance between two particles.Additionally, the rising temperature thickened the slipping layer (Table 1).Taken together, at a distance of twice the slipping layer thickness (the distance value of the left point of each curve in Figure 10), the particles must overcome a stronger repulsion force before flocculation at high temperature than at low temperature.Thus, the higher the temperature, the greater the capture frequency, as noted by Qiao et al.'s simulation of two primary particles [26].However, the situation was different when considering multiple particles.Also, taking the above three particles in the 30 °C case as an example, when the large particle (#721) met the floc formed by the other two particles, this particle seized the small one (#681) and a new floc was formed, with As shown in Figure 9a, with increasing temperature, η glo continued to grow up, while η cap began to decline slightly after its peak at about 20 • C after an initial increase.Sterling et al. [46] stated that the collision frequency was due to three factors: Brownian motion, turbulence shear, and differential settling.In this simulation, the sediment particles were much larger than the molecule and thus the Brownian motion could be ignored; the turbulence effect was omissible in still water; therefore, the main factor was the differential settling.Figure 9b showed a positive correlation between η glo and settling velocity, which was consistent with the viewpoint of Sterling et al. [46].η glo in the simulations ranges from 2.1% to 4.2%, very close to the results of Kim and Stolzenbach [45] (from 2.96% to 3.20%).The potential reasons for the difference between the two results might lie in Kim and Stolzenbach's [45] presumption that their model neglected the repulsive colloidal interaction, which is very sensitive to temperature (in Section 2.2) but was not mentioned in their study, so they might not take the temperature into account in their work.
When the distance between two particles was less than 25 nm, the three short-distance forces were the most important.Among the forces, only the electrostatic double-layer repulsion force is related to the temperature.From Figure 10a, the repulsive force was higher at a high temperature or a short distance between two particles.Additionally, the rising temperature thickened the slipping layer (Table 1).Taken together, at a distance of twice the slipping layer thickness (the distance value of the left point of each curve in Figure 10), the particles must overcome a stronger repulsion force before flocculation at high temperature than at low temperature.Thus, the higher the temperature, the greater the capture frequency, as noted by Qiao et al.'s simulation of two primary particles [26].
However, the situation was different when considering multiple particles.Also, taking the above three particles in the 30 • C case as an example, when the large particle (#721) met the floc formed by the other two particles, this particle seized the small one (#681) and a new floc was formed, with particle #689 settling individually.This collision produced a new floc, but disaggregated an old floc, giving no change in the number of flocs and particles; thus, η cap remained unchanged, while in the 20 • C case this collision made the η cap larger as it enlarged the floc.From the aspect of force, when a third particle was involved, the difference between the two double-layer forces acting on particle #681, F EL681-689 -F EL681-721 , was the same in magnitude as the particle gravity (Figure 10b), and the macroscopic force, including the gravity and hydrodynamic force, was highlighted.The higher the temperature, the more obvious the macro forces were.related to the temperature.From Figure 10a, the repulsive force was higher at a high temperature or a short distance between two particles.Additionally, the rising temperature thickened the slipping layer (Table 1).Taken together, at a distance of twice the slipping layer thickness (the distance value of the left point of each curve in Figure 10), the particles must overcome a stronger repulsion force before flocculation at high temperature than at low temperature.Thus, the higher the temperature, the greater the capture frequency, as noted by Qiao et al.'s simulation of two primary particles [26].However, the situation was different when considering multiple particles.Also, taking the above three particles in the 30 °C case as an example, when the large particle (#721) met the floc formed by the other two particles, this particle seized the small one (#681) and a new floc was formed, with This phenomenon reveals that, at high temperature, few large flocs broke into small ones because of collision, resulting in a fluctuation of time series and a decrease in the maximum floc size, as illustrated in Figure 4a.However, the mean floc size and total floc volume did not decrease (in Figure 4b) because the large flocs that were collapsed by collision had a very small proportion and were only decreased in size but rarely dispersed into primary particles.
A detailed force analysis of a large floc is more convincing but more complex when it is comprised of more primary particles because of its continuous changing spatial statuses and complicated interactive forces.Although the analysis of the floc of three primary particles is limited, it illustrates the flocculation formation and breakage at different temperatures.From the analysis, it is certain that, at high temperature, the macroscopic force becomes more obvious due to the small difference of short-distance forces between particles.The above analysis qualitatively explains why the fracture frequency increases at high temperatures, which is the main cause of the decrease in η cap .However, the reason is different for the low-temperature case, where the large repulsion force makes it more difficult for the particles to get close enough to form flocs, resulting in a small η cap .
The above results could be explained by the present LBM model, but there is still an obvious shortage in that the flocs in each case were not sufficient in quantity and size.This shortage is because the LBM requires a large number of grids to describe the solid-fluid boundary and micro properties of each sediment particle; thus, the computational cost of the LBM is extremely high.Each case in this study took three months on a supercomputer with 288 CPUs.A higher computational cost is expected if the sediment bulk properties are described in more detail.That is the reason for the less convincing results of the present study.Therefore, subsequent studies will involve further optimizing the calculation method and case design for better simulation results.In a natural environment, changes in other factors, such as water salinity, turbulence, and SSC, are unavoidable, so their influences on flocculation will be further studied by LBM in the future.

Conclusions
The bulk settling and flocculation process of cohesive sediment containing primary particles of sizes of 510.0 µm in still water at various temperatures was simulated via the LBM.The floc size and volume were analyzed, and the effect of temperature on bulk settling was studied based on the macroscopic SSC change and microscopic statistics of particle and floc settling velocity.The difference of the properties at different temperatures was explained by the formation of flocs, aggregation frequency, and forces between particles.The following conclusions were obtained: parameter collision operator is used in the present study.The post-collision distribution fi * can be written in the same form as fi eq :

Figure 1 .
Figure 1.The possible velocity directions in the D3Q19 topology.The speed of discrete velocity |e i | equals 0 when particle keeps its original position after a time step, equals 1 when it moves to the faces of the cubic and √ 2 to the edges of the cubic.

Figure 2 .
Figure 2. The XDLVO potentials of different net distances hij between particles.The Lifshitz-van der Waals potential is always attractive; the electrostatic double-layer potential is repulsive and the Lewis acid-base potential's sign depends on the properties of colloids.XDLVO potential is the summation of the aforementioned three potentials, resulting in potential wells and potential barriers.Particles whose interactive forces conquer the potential barrier can form a stable floc.

Figure 2 .
Figure 2. The XDLVO potentials of different net distances h ij between particles.The Lifshitz-van der Waals potential is always attractive; the electrostatic double-layer potential is repulsive and the Lewis acid-base potential's sign depends on the properties of colloids.XDLVO potential is the summation of the aforementioned three potentials, resulting in potential wells and potential barriers.Particles whose interactive forces conquer the potential barrier can form a stable floc.

Figure 3 .
Figure 3. Illustration of slipping layer thickness, equaling the thickness of the electrostatic double layer in Yang et al.'s study [37].This zone is comprised of an absorbed layer and a diffuse layer.The slipping layer thickness can be calculated by the potentials ψ0 on the particle surface and ζ at the distance δ as the potential decays exponentially.When the particles are close enough and wrapped by a common slipping layer, they are considered as forming a new or enlarging an old floc.

Figure 3 .
Figure 3. Illustration of slipping layer thickness, equaling the thickness of the electrostatic double layer in Yang et al.'s study[37].This zone is comprised of an absorbed layer and a diffuse layer.The slipping layer thickness can be calculated by the potentials ψ 0 on the particle surface and ζ at the distance δ as the potential decays exponentially.When the particles are close enough and wrapped by a common slipping layer, they are considered as forming a new or enlarging an old floc.
shows the time series of the maximum and mean floc sizes.The final results of them and the floc volume are shown in Figure 4b.Floc volume is a sum of the volumes of all the flocs, and indicates the content of floc in the suspended column.During the same settling period, both the floc sizes and floc volume increased with increasing temperature.Water 2019, 11, 945 7 of 16 and the floc volume are shown in Figure 4b.Floc volume is a sum of the volumes of all the flocs, and indicates the content of floc in the suspended column.During the same settling period, both the floc sizes and floc volume increased with increasing temperature.(a) (b)

Figure 4 .
Figure 4. Floc properties of each case.(a) Time histories of maximum floc size (solid symbols) and mean floc size (hollow symbols) for the 5 °C, 10 °C, 20 °C, and 30 °C cases; (b) The maximum floc diameter, mean floc diameter (the left axis, with unit of μm) and floc volume (the right axis, with unit of 10 −15 m 3 ) at the end of the simulation.In (b), the thin solid line and dashed line represent the line fit of maximum and mean floc diameters, respectively, and the thick solid line is the exponential fit of floc volume.

Figure 4 .
Figure 4. Floc properties of each case.(a) Time histories of maximum floc size (solid symbols) and mean floc size (hollow symbols) for the 5 • C, 10 • C, 20 • C, and 30 • C cases; (b) The maximum floc diameter, mean floc diameter (the left axis, with unit of µm) and floc volume (the right axis, with unit of 10 −15 m 3 ) at the end of the simulation.In (b), the thin solid line and dashed line represent the line fit of maximum and mean floc diameters, respectively, and the thick solid line is the exponential fit of floc volume.At the end of the simulation, in the cases of 5 • C, 10 • C, 20 • C, and 30 • C, the maximum floc sizes were 20.0 µm, 24.0 µm, 22.3 µm, and 23.3 µm, respectively, with the peak value at 10 • C. In the cases of 20 • C and 30• C, the maximum floc sizes fluctuated slightly in the latter phase (Figure4a).The mean floc size increased slightly but more stably than the maximum size did, from 15.2 µm at 5 • C to 15.5 µm at 10 • C, 15.8 µm at 20 • C, and 15.9 µm at 30 • C. The floc volume increased more obviously than the floc size.They were 3.8% at 5 • C, 4.0% at 10 • C, 4.3% at 20 • C, and 5.5% at 30 • C (Figure4b).

Figure 5 .
Figure 5.Time series of particle velocities and particles' net distance, i.e., h ij in Equation (8) and r ij − R i − R j in Figure 3. Solid symbols are net distances, and hollow symbols mean velocities.d1 represents the distance between particles #681 and #689; d2 represents that between particles #689 and #721; d3 represents that between particles #681 and #721; u1, u2, and u3 represent the settling velocities of particles #681, #689, and #721, respectively.The left axes are for the distance, and the right ones are for the velocities.(a) T = 5 • C; (b) T = 10 • C; (c) T = 20 • C; and (d) T = 30 • C.

Figure 6 .
Figure 6.SSC of each water depth in (a) 5 °C, (b) 10 °C, (c) 20 °C, and (d) 30 °C.The symbols between two labeled depths represent the SSC between those two depths.The SSC can be statistically obtained by the sediment weight in this zone.It was generally stable at depth of 4.5 mm, but there were some small differences between the cases.

Figure 6 .
Figure 6.SSC of each water depth in (a) 5 • C, (b) 10 • C, (c) 20 • C, and (d) 30• C. The symbols between two labeled depths represent the SSC between those two depths.The SSC can be statistically obtained by the sediment weight in this zone.It was generally stable at depth of 4.5 mm, but there were some small differences between the cases.

Figure 7 .
Figure 7.The SSC time series and its relative parameters.(a) The time series curves of relative SSC (ratio of c(t) at time t to initial c0) of a water system shallower than 4.5 mm; (b) The kc in Equation (12) and t1/2 in Equation (13).In (a), the symbols are the experiment results, the solid lines and dashed lines are the fitting curves for n = 0.3 and n = 1.0, respectively.In (b), the solid square and hollow triangles symbols are t1/2 and kc, both of which are deduced from the solid curves in (a).The solid line and dashed line in (b) are the linear fittings for t1/2 (t1/2 = −0.380T+ 21.366) and kc (kc = 0.120T + 2.661).

Figure 7 .
Figure 7.The SSC time series and its relative parameters.(a) The time series curves of relative SSC (ratio of c(t) at time t to initial c 0 ) of a water system shallower than 4.5 mm; (b) The k c in Equation (12) and t 1/2 in Equation (13).In (a), the symbols are the experiment results, the solid lines and dashed lines are the fitting curves for n = 0.3 and n = 1.0, respectively.In (b), the solid square and hollow triangles symbols are t 1/2 and k c , both of which are deduced from the solid curves in (a).The solid line and dashed line in (b) are the linear fittings for t 1/2 (t 1/2 = −0.380T+ 21.366) and k c (k c = 0.120T + 2.661).

Figure 8 .
Figure 8.The settling velocities of flocs of different diameter and temperature in many studies.(a) Floc diameter and settling velocities from the observational data of Xia et al. [12], Khelifa and Hill [43], Guo and He [13], Dyer and Manning[42], and Manning et al.[43]; theory results of Manning et al.[43], Winterwerp[44], and Khelifa and Hill[43]; and the simulation results of LBM under temperature conditions ranging from 5 °C to 30 °C; (b) floc velocities and temperature from the observational data of Xia et al.[12], Wan et al.[6],Guo and He [13], and this simulation's results.

Figure 8 .
Figure 8.The settling velocities of flocs of different diameter and temperature in many studies.(a) Floc diameter and settling velocities from the observational data of Xia et al. [12], Khelifa and Hill [43], Guo and He [13], Dyer and Manning [42], and Manning et al. [43]; theory results of Manning et al. [43], Winterwerp [44], and Khelifa and Hill [43]; and the simulation results of LBM under temperature conditions ranging from 5 • C to 30 • C; (b) floc velocities and temperature from the observational data of Xia et al. [12], Wan et al. [6], Guo and He [13], and this simulation's results.

Figure 9 .
Figure 9.The change of Collision frequency ηglo, capture frequency ηcap and aggregation frequency Eagg with (a) temperature and (b) settling velocity.

Figure 10 .
Figure 10.Comparisons of FEL and the gravity of the primary particle with diameters of 9.36 μm and efficient densities ρs-ρ of 1650 kg/m 3 .(a) The electrostatic double layer force between particles #681 and #689, F EL 681-689 of different particle net distances at temperatures of 5 °C, 10 °C, 20 °C, and 30 °C; (b) the difference in electrostatic double-layer forces acting on particle #681, F EL 681-689-F EL 681-721 with different particle net distances in the cases of 5 °C, 10 °C, 20 °C, and 30 °C.

Figure 9 .
Figure 9.The change of Collision frequency η glo , capture frequency η cap and aggregation frequency E agg with (a) temperature and (b) settling velocity.

Figure 10 .
Figure 10.Comparisons of FEL and the gravity of the primary particle with diameters of 9.36 μm and efficient densities ρs-ρ of 1650 kg/m 3 .(a) The electrostatic double layer force between particles #681 and #689, F EL 681-689 of different particle net distances at temperatures of 5 °C, 10 °C, 20 °C, and 30 °C; (b) the difference in electrostatic double-layer forces acting on particle #681, F EL 681-689-F EL 681-721 with different particle net distances in the cases of 5 °C, 10 °C, 20 °C, and 30 °C.

Figure 10 .
Figure 10.Comparisons of F EL and the gravity of the primary particle with diameters of 9.36 µm and efficient densities ρ s -ρ of 1650 kg/m 3 .(a) The electrostatic double layer force between particles #681 and #689, F EL 681-689 of different particle net distances at temperatures of 5 • C, 10 • C, 20 • C, and 30 • C; (b) the difference in electrostatic double-layer forces acting on particle #681, F EL 681-689 -F EL 681-721 with different particle net distances in the cases of 5 • C, 10 • C, 20 • C, and 30 • C.

Table 1 .
Parameters in each case.

Table 1 .
Parameters in each case.