The 1DV two-phase model used in this study is a simplified version of Xu and Dong (2017b). Eq. 1 represents the momentum equation for both the fluid and solid phases in the horizontal direction, Eqs 2, 3 represent the continuity equations, and Eqs 4, 5 represent the momentum equations for both phases in the vertical direction. Note that Eq. 1 only applies to field measurements where the rate of change of mean horizontal flow velocity (∂U/∂t) and horizontal pressure (∂P/∂x) need to be considered. For modelling the simplified case of cohesive sediment settling vertically (i.e., within a settling column), the horizontal flow terms and terms involving horizontal gradients are omitted.

Within Eqs 1–6, U is the horizontal velocity for both phases (i.e., fluid phase is denoted with subscript f and the solid phase with subscript S), ρ_mix = α_sρ_s + α_fρ_f is the bulk density of the fluid-sediment mixture, α_s and α_f are the volume fractions of solid and fluid phase, ρ_s and ρ_f are the solid and fluid phase densities, respectively, _T is time, w_s and w_f are the floc settling velocities and fluid velocities, respectively, P is the pressure of mixture (with p_f corresponding to the fluid pressure), g is the gravitational acceleration, τ_v is the viscous shear stress of the mixture, and f_i is the momentum transfer between two phases. In this study, f_i is used to describe the drag force from the other phase (i.e., the drag force exerted on the fluid phase from the solid phase, or vice versa). The modified classical mixing length method is adopted to calculate turbulent eddy viscosity (υ_T) and eddy diffusivity (Γ_T):

where σ_T is the turbulent Prandtl-Schmidt number (usually specified as 0.7 or 1.0), κ is the Karman constant, F_ν and F_d are the correction coefficients for eddy viscosity and eddy diffusivity, respectively, to describe the buoyancy effects caused by suspended sediments, and h is the height of vertical water column. Here, the eddy viscosity is modified by the formulation presented by Busch (1973), while the Munk-Anderson formula is applied for the calculation of F_d:

where, Ri is the gradient Richardson number, defined as:

Here, we assume the shear stress for the solid and fluid phases are equal (Chauchat et al., 2013; Xu and Dong, 2017b), and is presented as:

where um=αfuf+αsus is the volume-averaged velocity and μ_mix = μ_f(1 + β_aα_s) is the augmented viscosity, where β_a is the amplification factor. With an increase of the solid fraction, the mixture goes through the transition from Newtonian to non-Newtonian fluid. To account for the non-Newtonian effects, the amplification factor β_a is specified as (Graham, 1981):

where d* is defined as non-dimensional inter-particle distance. From geometrical considerations, it is expressed as a function of sediment volumetric concentration d*=[1-(αs/αsmax)1/3]/(αs/αsmax)1/3, where αsmax=0.625 is the maximum solid volume of simple cubic packed spheres (Chauchat et al., 2013). The calculated viscosity from Eqs 12, 13 are suitable for sediment transport with large variation of sediment concentration, as the model results are consistent with results from both the classic formulae μ_mix = μ_f(1 + 2.5α_s) and μmix=μf9/8[(αsmax/αs)1/3-1]-1 for the dilute case (Einstein, 1905) and for the dense case (Frankel and Acrivos, 1967), respectively.

In considering the aggregation and break up of flocs, Chauchat et al. (2013) suggested that the drag force should be given from a macroscopic point of view for the two-phase model. As the inverse of water flow resistance can be measured using the permeability parameter K, here the generalized Darcy law is adopted to describe the drag force (Toorman, 1996):

Permeability K is usually applied only when the sediment concentration reaches the gelling concentration (Winterwerp and Van Kesteren, 2004). Based on the stress balance equation, Toorman (1999) also extended the permeability K to the cases of dilute sediment concentration, the sedimentation and consolidation processes giving a unified expression as:

where, W is the settling velocity including the hindered settling effects and specified as:

where, ϕ_f is the volumetric concentration of cohesive sediment flocs and ϕ_fmax is the maximum value. The fractal dimension is denoted as nf. In the right-hand side of Eq. 16, the first two terms represent the effects of buoyancy, viscosity, and wake on the settling process of sediment particles. The ϕ_fmax is introduced to describe that the settling velocity of sediment particles approaches zero when ϕ_f is approaching ϕ_fmax. Following Chauchat et al. (2013), the value of 0.85 is adopted for ϕ_fmax, while w₀ is the settling velocities of cohesive sediment flocs in the dilute case. To be consistent with the flocculation models adopted in this work, the settling velocities of cohesive sediment flocs are calculated based on fractal theory presented by Winterwerp (1998) as follows:

where, α₁ is a coefficient depending on the sphericity of cohesive sediment flocs, Re is the particle Reynold number and defined as R_e = w_sD/ν, with D being the representative sizes of flocs. The boundary condition for sediment concentration, which also serves as the bed erodibility, is specified by van der Ham and Winterwerp (2001):

where, τ_b is the bed shear stress, τ_cr is the critical bed shear stress for sediment erosion, and M is the erosion coefficient.

The two flocculation models outlined in section “Two-Phase Flow Model” relate to the time evolution of a representative floc size, while the two-phase model calculates the SSC. Therefore, the number density N of flocs may be introduced as an intermediate variable to link these models. The volumetric floc concentration ϕ_f can be linked to the number concentration of flocs N via the equation:

where, f_s is a floc shape factor. During the flocculation process, ϕ_f varies with floc size D and fractal dimension nf, and can be calculated from the sediment volumetric concentration α_s as follows:

where, ρ_floc is the density of flocs. According to fractal theory, the floc density can be presented as (Kranenburg, 1994):

From Eqs 23–25, the variable floc size D can therefore be determined if the sediment volumetric concentration α_s and the number concentration N of flocs are known. Furthermore, the settling velocity of cohesive sediment flocs in a dilute suspension w₀ can be calculated from Eq. 17. Therefore, the floc settling velocities can be linked to their number concentration N. As discussed above, the floc number concentration N also needs to be resolved. Here, following Winterwerp (2002), we propose the balance equation for number density as:

where Γ_T is the turbulent diffusion coefficient and F_N is the flocculation term. The two flocculation models (i.e., Model A and Model B, section “Cohesive Sediment Flocculation models”) are examined, in turn, by combining each with the 1DV two-phase model (see section “Two-Phase Flow Model”). These flocculation models are first-order differential equations for floc size D, while F_N is in the form of a first-order differential equation for number density N. As such, Eq. 19 can be rewritten as (see Winterwerp, 1998 for more details):

Based on Eq. 27, flocculation Model B has the form:

The flow chart in Figure 1 shows the coupling procedures between the flocculation models and the 1DV two-phase model. For each time step, the governing equations of the two-phase model (Eqs 1–6) are firstly solved to obtain the sediment concentration. This concentration, and other relevant parameters, are then input into the flocculation models to solve the number density equation (Eq. 26). From Eqs 23–25, information on the floc size D, fractal dimension nf and floc density ρ_floc is obtained. Based on fractal theory, the settling velocities w₀ of the cohesive sediment flocs are then calculated by Eq. 17. Finally, these settling velocities are used to determine the drag force closure for the 1DV two-phase model (Eqs 14–16).

As discussed above in the introduction, most previous laboratory experiments on cohesive sediment flocculation have been conducted under controlled, idealized conditions within mixing tanks with pre-determined constant sediment concentrations, turbulent shear rates and/or water salinities. However, under non-equilibrium conditions, where sediment concentrations vary in both time and space, flocculation processes become more complicated due to the relative influence of residence and flocculation times on the floc sizes generated.

For this reason, the 1DV two-phase model is applied to simulate recent grid-stirred flocculation experiments conducted by Cuthbertson et al. (2018) for pure kaolin clay suspensions within a vertical, grid-stirred settling column (details of the experimental arrangement are given in Cuthbertson et al., 2010 and Cuthbertson et al., 2018). In this case, the calculations are focused in the vertical direction, therefore horizontal flow terms and other terms involving horizontal gradients in the two-phase model (Eqs 1–6) are omitted. During individual experimental runs, a highly concentrated kaolin suspension was fed at a constant inflow rate via a peristaltic pump from an external mixing tank into the upper buffer mixing tank at the top of the main grid-stirred settling column section. Two counter-rotating mixing paddles within the buffer mixing tank generated an established circulation that diluted the kaolin suspension within the preset volume water (50 L) and gradually transferred the dilute clay suspension into the main column section via a gate opening. Time series measurements of sediment concentrations were collected using optical backscatter (OBS) probes located at 0.5 and 1.2 m above the bottom of the main column section. These OBS probes were calibrated over a wide range of pure kaolin clay suspensions (with mass concentrations ranging from C = 0–1 g.l^–1), and relationships were established between turbidity (NTU) and suspended sediment concentration (Cuthbertson et al., 2018). The time evolution of floc sizes was collected 0.4 m above the base of the column, via a macro-CCD camera (see Cuthbertson et al., 2018).

Three datasets from the laboratory settling column experiments, denoted Cases 1–3, are used for validation of the 1DV two-phase flocculation model. In the model simulations of the settling column cases, the temporal variation of clay concentration at the upper model boundary is determined by specifying (i) the clay input conditions (see Table 1) and (ii) the upper buffer tank volume and specified mass transfer rate of clay from the buffer tank to the main column (i.e., through calibration with time series clay concentrations measured within the column by the OBS probes). The initial floc size of the clay suspension is set as the primary clay particle size d = 2.0 μm, which is a regarded as a conservative value as it assumes no flocculation occurs in the buffer tank. The turbulent shear rate G adopted in the simulations, and representing the turbulence intensity, is set as a constant value for each case (see Table 1). These represent the average shear rate values within the central flow region between the oscillating grid pairs (Cuthbertson et al., 2010), which vary depending on the grid oscillation stroke and frequency (for the fixed grid arrangement). The resulting zero-mean shear turbulence fields are demonstrated to be quasi-homogeneous and near-isotropic within the central flow region between the oscillating grid pairs (i.e., away from the grids themselves). To determine the mean fractal dimension within the two flocculation, according to Cuthbertson et al. (2010), for floc sizes of pure clay smaller than 100, the majority of the fractal dimensions lie in the range of 1.7 ≤ nf ≤ 2.0 (see Figure 11 in Cuthbertson et al., 2010). As in a normal distribution the probability P(μ−3σ < x < μ + 3σ) is larger than 99%, the standard deviation is estimated as (2.0-1.7)/6 = 0.05. The mean fractal dimension is adopted as 1.85. The constant fractal dimension in Model A is thus specified as 1.85, while the specific flocculation parameters adopted in Models A and B for the three experimental cases considered are summarized in Table 2.

To ensure a rational comparison between the two flocculation models (i.e., Model A and B) for the reproduced time evolution of clay flocs in the settling column experiments, it is necessary to establish the baseline model parameters through calibration. Here, the flocculation model coefficients are obtained by adopting the same final equilibrium floc size generated under the same fixed sediment concentration and turbulent shear rate for each data set. Under these steady-state conditions, the ratio between the aggregation and break-up parameters, kA′ and kB′, can be determined. Secondly, the value for kA′ is selected (i.e., so is the value of kB′, because the ratio of these two kA′ and kB′ has been determined) to fit best to the initial flocculation rate.

The time series measurements and model results of sediment concentration at 0.5 and 1.2 m above the bottom of the main grid-stirred settling column section are shown in Figure 2 for Case 1 (Table 1). Using these measured time series of sediment concentration to calibrate the upper clay input boundary condition (where t = 0 refers to start of the sediment input into the column), the two flocculation models are capable of reproducing the vertical profiles of sediment concentration. In the experimental data, the measured concentrations at 0.5 and 1.2 m converge around 12,000 s (Figure 2), with the results from both Models A and B converging around 13,000 s. It is also noted that a smaller vertical gradient of sediment concentration was obtained by Model B than that of Model A before convergence. Similar trends were also obtained in the model simulations of Cases 2 and 3 (Table 1).

The measured and modelled temporal variations in the root-mean-square (rms) floc sizes generated in the settling column at z = 0.4 m, where the floc size measurements were obtained, are shown in Figures 3A–C for Cases 1–3, respectively. The main feature of these measured data is that near quasi-equilibrium floc sizes are already attained within the column by the time the flocs are first detected in the floc viewing chamber within the lower part of the settling column (Cuthbertson et al., 2010). The corresponding 1DV two-phase flocculation model results indicate that Model B (i.e., multiple fractal dimension) provides far closer agreement with the measured time evolution of rms floc sizes, both in terms of the initial rapid flocculation and equilibrium floc size attained, while Model A significantly underpredicts the initial flocculation rate before reaching the same equilibrium floc size at a later elapsed time. Indeed, Model A is shown to be incapable of reproducing the measured temporal evolution with of floc sizes with the settling column no matter what combination of k′A and kB′ are used. Specifically, the root-mean-square errors (RMSE) of the calculated time series of floc sizes are 18.5 (Case 1), 26.6 (Case 2), and 26.9 (Case 3) for the results of Model A. While for the results of Model B, the RMSEs are 4.3 (Case 1), 3.5 (Case 2), and 5.7 (Case 3), respectively. By incorporating multiple fractal dimensions and thus variable yield strengths, Model B is able to capture better the temporal characteristics of the rapid initial floc size adjustment at earlier stages of the runs (i.e., T < 5000 s), after which the calculated floc sizes increase only slowly and approach the quasi-equilibrium floc size.

To further demonstrate the temporal evolution of clay flocs in the settling column simulations, vertical profiles of floc sizes calculated by Models A and B are compared in Figures 4A–C with the floc size measurements (at z = 0.4 m in the settling column) at three different elapsed times for Case 3 (Table 1). In addition, the corresponding calculated vertical distributions of SSC are compared with the measured OBS data obtained at the two elevations (z = 0.5 and 1.2 m) at the same elapsed times in Figures 4D–F (Note: as the vertical profiles of sediment concentration predicted by Models A and B are very similar, only Model B profiles are shown). In all cases, these reproduced concentration distributions show excellent agreement with the equivalent measured concentrations at the two elevations in the column.

During the experiments, the sediment feed into, and mixing within, the buffer mixing tank was assumed to keep cohesive sediments in a largely unflocculated state (Cuthbertson et al., 2018), resulting in very small floc sizes at the very top of the column over the experimental run duration. In this modelling study of the settling column experiments, the effects of the volumetric floc flux on the evolution of local floc size are fully accounted for by solving sediment concentration and floc number concentration simultaneously. With this imposed upper boundary condition, the predicted maximum floc sizes are shown to occur in the upper-middle part of the settling column [i.e., z/h = ∼0.7 at T = 2000 s and z/h = ∼0.8 at T = 5000 s (for Models A and B), Figures 4A,B]. This floc size maximum is therefore shown to occur at a different elevation within the settling column compared to the elevation of the maximum SSC (i.e., z/h = 1.0) at these elapsed times.

At the later stages of the experimental run (i.e., T = 15000 s, Figure 4C), a quasi-equilibrium value of the floc size is reproduced by Models A and B in the vertical column profile, which becomes consistent with the overall shape of the sediment concentration distribution in the column (Figure 4F). However, as indicated previously, during earlier stages of the experimental run (i.e., T = 5000 s, Figure 4B), only the floc size results from Model B agree well with the measured floc size data and is, thus, to only model capable of reproducing the rapid initial floc size adjustment within the settling column.

In terms of the validation of flocculation models against field data from estuarine sites, it is difficult to find complete and synchronous datasets that include all hydrodynamic conditions (i.e., flow velocities, bed shear stresses), suspended sediment concentrations (SSC), and the physical characteristics of flocs (i.e., floc sizes, settling velocities) generated throughout the tidal cycle. Van Der Ham et al. (2001) reported high frequency SSC measurements and flow velocities in the tidal channel of the Ems/Dollard estuary over a 24 h period. Within this measurement area, the horizontal gradients of SSC are known to be negligible, with horizontal and vertical salinity gradients also small when the river discharge is low (Van Der Ham et al., 2001), making this an appropriate field site for the application of the 1DV two-phase model developed in the current study. This data set alone, however, cannot provide full validation of the flocculation model as no corresponding floc information was available over the same time period covered by Van Der Ham et al. (2001).

Most recent field studies on floc characterization in estuaries have tended to focus on relating measured floc sizes directly to their corresponding settling velocities. Within this context, our multiple fractal dimension model (Model B) has been validated against such field-based floc measurements (i.e., Khelifa and Hill, 2006 data) within Xu and Dong (2017a). Furthermore, previous field studies conducted by Dyer et al. (2000) in the Ems/Dollard estuary also provide measured floc sizes and settling velocities during the flood phase of the tidal cycle, that can be compared qualitatively with the current model predictions, albeit under different tidal conditions. These findings and comparisons are discussed in detail within section “Model Application to Field Measurements” of the paper.

For the current field application, the empirical coefficients for both flocculation models need to be calibrated again. Based on previous studies, we assume that the sediment density ρ_s = 2650 kg/m³ and the size of primary particles is set as d = 4 μm (Winterwerp, 1998). For direct comparison between the predictive capabilities of models A and B, the same baseline conditions need to be set. Thus, the two flocculation models are calibrated to achieve the same equilibrium floc size of 300 μm (Figure 5A), under the shear rate condition of G = 2 s^–1 and for a sediment concentration c = 0.3 kg/m³. These represent appropriate values for the typical field measurement conditions found in the tidal channel of Ems/Dollard estuary (Van Der Ham et al., 2001). It is noted that, in the modelling of the settling column experiments conducted with pure kaolin clay suspensions, the constant fractal dimension (in Model A) and mean fractal dimension (in Model B) remained unchanged with an increase in the floc size. By contrast, Khelifa and Hill (2006) collected more than 26 laboratory and field site measurements of flocs to assess the size-dependency of fractal dimensions; their results suggesting that the fractal dimension decreases with increasing floc size. Thus, Eq. 29 is adopted for the calculation of both the constant fractal dimension (in Model A) and the mean fractal dimension (in Model B) for their application to field measurements (see Figure 5B). It should be noted that for a given floc size D, the fractal dimension is unique in Model A, while, in Model B multiple fractal dimensions are adopted.

where, _α and _β are coefficients and specified using following boundary conditions:

where, nf_c is a characteristic fractal dimension when floc size D equals a characteristic floc size D_c. In Eq. 29, the fractal dimension takes the maximum value of 3 when floc size approaches the primary particle size d and takes a lower value nf_c when floc size approaches the characteristic floc size D_c. The value of nf_c = 2.0 is adopted when the characteristic floc size D_c → 300 μm, which is the typical value selected to calibrate the flocculation models (e.g., Winterwerp, 1998). As to the variance of the fractal dimension, a logarithmic function is found to be physically more realistic (Vahedi and Gorczyca, 2012):

where, α₂ is an empirical coefficient. When the floc size approaches primary particle size, the variance in fractal dimension is assumed to be zero (i.e., nf = 3.0), while for floc sizes approaching the characteristic size _Dc, it is set at 0.6 (Winterwerp, 1998) (i.e., 1.7 ≤ nf ≤ 2.3). Therefore, according to Eq. 31, _α2 can be determined as 0.0174. In contrast to the settling column experiments, the shear rate G = ε/ν within the tidal channel is no longer constant, instead varying with the tidal cycles. As such, Eqs 32, 33 are adopted to describe the turbulent kinetic energy k and dissipation ε, as follows:

where L is the Prandtl’s mixing length, C_D and C_u are set at 0.1925 and 0.09 (Rodi, 1980), respectively. The coefficients adopted in the two flocculation models are summarized in Table 3. In terms of the two-phase model, following van der Ham and Winterwerp (2001), the erosion rate for the cohesive sediment bed M = 1.5×10⁻⁸ m/s is selected. The critical shear stress for the cohesive sediment erosion τ_cr = 0.1 Pa, which is the averaged critical shear stress suggested by Kornman and De Deckere (1998), based on erosion studies conducted at an adjacent tidal flat in the Ems/Dollard estuary. The critical shear stress for the deposition is also specified as τ_b = 0.1 Pa, while a maximum depth-averaged sediment concentration C_max = 0.5 kg.m^–3 is applied in both models to account for the limited sediment availability from the bed (van der Ham and Winterwerp, 2001).

Figure 6A presents the time series measurements of depth-averaged velocities and elevations over a 24 h period, while Figures 6B–D present corresponding measured velocities (red circles) and modelled velocity profiles (black lines, calculated by Model B) at three different elapsed times. These represent hydrodynamic conditions at (high) slack water (08:00, Figure 6B), 1 h after peak ebb flow (12:00, Figure 6C), and 1 h before peak flood flow (16:00, Figure 6D), respectively. The measured velocities (red circles) are obtained at elevations of 0.1, 0.4, and 1.0 m above the bed surface. As the relative height z/h is adopted for the vertical axis, and the overall water depth h varies over the measurement duration (i.e., see elevations in Figure 6A), the velocity measurements are located at different relative heights in the individual figures. The RMSEs of the calculated time series of velocities and shear stresses are 0.163 and 0.115, respectively, for the results of Model B. For Model A, RMSEs are 0.165 (velocities) and 0.114 (shear stresses). The predicted results from Model B compare very well with the measured data (Note: equivalent results from Model A are found to be very similar and, as such, are not shown here). The measured and predicted shear stresses (calculated by Model B) at 0.4 m above the bed are shown in Figure 7. Again, the equivalent results predicted by Model A are very similar (not plotted) and thus both models are capable of reproducing the velocity profiles and shear stresses during the different tidal phases.

The measured and modelled SSC timeseries at elevations of 0.3 and 1.4 m above the bed level, corresponding to the same period of measured depth-averaged velocities and elevations (Figure 8A), are shown in Figures 8B,C, respectively. The coefficients in both flocculation models are calibrated here to enable the model results to capture the peak SSC values at the 0.3 m elevation above the bed. As such, the SSC time series produced by both models (i.e., Models A and B) at 0.3 and 1.4 m above the bed are shown to broadly follow the temporal trend of the measured SSC data. The exception to this is around 04:00 and just after 16:00 in the SSC measurements at 0.3 m (Figure 8B), where there are abrupt increases in SSC values [note: similar, but less abrupt increases are also seen 03:00 and 16:00 in SSC measurements at 1.4 m (Figure 8C)]. These abrupt changes in SSC are explained as local increases in sediment availability (van der Ham and Winterwerp, 2001), while similar maximum SSC values occur during the flood and ebb tides, despite larger shear stresses being generated during the ebb tide (Figure 7B). When compared with the results of Model A, Model B showed better overall prediction and fit to these field measurements. The RMSE values for the SSC timeseries results from both models, when compared directly with the field measurements are calculated at the 0.3 and 1.4 m elevation above the bed as 0.296 and 0.177 (for Model A) and 0.223 and 0.130 (for Model B).

According to Van Der Ham et al. (2001), a main feature of the measured concentration data is a small vertical gradient in SSC values that suggests well-mixed conditions exist within the estuary (at least in terms of SSC). The results from Model B again appear to capture this feature best [e.g., at around 13:00 (i.e., ebb tidal phase), the difference in calculated SSC values at the 0.3 and 1.4 m elevations is 0.8 kg.m^–3 for Model A and 0.5 kg.m^–3 for Model B, see Figures 8B,C]. TheRMSEs for both models predictions of SSC are calculated when compared directly with the field measurements. At the 0.3 m elevation above the bed, the RMSE values for Model A and B predictions are 0.296 and 0.223, respectively, while at 1.4 m elevation, the corresponding RMSE values are 0.177 and 0.130, respectively. To further illustrate the vertical structure of physical properties predicted Model B, vertical profiles of SSC during both the slack (high) water period and subsequent peak (ebb) tidal velocity period are presented in Figures 9A,B. For slack water conditions, Model B results show lower SSC values remain in the upper part of the water column, with larger SSC gradients formed in the near-bed flow region (Figure 9A). By contrast, during the peak ebb tidal velocity period, the vertical distribution of SSC represents well mixed conditions (Figure 9B). Overall, the suspended sediment concentration profiles predicted by Model B match well the measured SSC data at 0.3 and 1.4 m elevations above the bed (i.e., black triangles, Figures 9A,B).

In the simulation of the grid-stirred settling column experiments with pure kaolin clay suspensions, the development of sediment concentration profiles within the column was shown to be well-represented by the 1DV two-phase model with either of the two flocculation models (i.e., with fixed or variable fractal dimensions) incorporated (Figure 2). By contrast, significant variability in the temporal development of rms floc sizes between the two models suggested that the adoption of a multiple fractal dimension approach (i.e., Model B) better replicated the floc size development in the settling column. Though direct measurements of floc settling velocities of kaolin clay generated in the grid-stirred settling column are not presented by Cuthbertson et al. (2018), they are shown in a previous study by Cuthbertson et al. (2010). Here, the measured kaolin clay flocs sizes and their corresponding settling velocities lay between two predicted settling rate curves (with fractal dimensions nf = 1.7 and nf = 2.3, see Figure 11 in Cuthbertson et al., 2010). The corresponding Model B results for floc sizes and settling velocities measured in the settling column tests were found to be consistent with this conclusion.

A conservative value of the initial floc size (D₀ = 2 μm) was adopted for simulations with both flocculation models. However, sensitivity of the model predictions to the specification of D₀ needs to be tested. Figure 10 shows sensitivity analysis runs of the predicted temporal development of the rms floc size for both flocculation models, where D₀ is set at 2, 5, and 10 μm. It is apparent that the different D₀ values influence floc development in both models, particularly during the initial stages of floc size evolution. Specifically, by increasing the initial floc size D₀ in Model B, the initial rapid floc size development occurs earlier, with the final floc adjustment to quasi-equilibrium floc sizes shown to converge for all D₀ after T = ∼5000 s. The initial rapid growth in floc sizes occurs as smaller flocs, with higher density and larger yield strengths, are more difficult to break up [i.e., with the aggregation term in Model B (Eq. 22) thus dominant]. Specification of larger D₀ therefore takes a shorter time to reach floc sizes where the break-up term in Model B (Eq. 22) becomes more important (i.e., represented by the change in gradient of the temporal floc size development) and floc sizes then adjust more gradually to their quasi-equilibrium floc size. For Model A, the effect of D₀ on the initial floc development is less consistent. Indeed, when D₀ is set at 10.0 μm (i.e., blue dotted line, Figure 10), the predicted rms floc size actually decreases initially before increasing steadily with time. This floc size reduction is due to the sediment concentration being initially very low in the column, resulting in low aggregation rates, while the initial break-up term for the D₀ = 10 μm flocs is higher (i.e., break-up > aggregation on right-hand side of Eq. 19). This initial reduction in rms floc size also means that convergence with the temporal floc size evolution for D₀ = 5 μm occurs significantly earlier than with the D₀ = 2 μm condition. As with Model B, once the temporal development of rms floc sizes have converged for all D₀ values (at T = ∼8000 s, Figure 10), the subsequent adjustment to the quasi-equilibrium floc size again also coincide.

To better explain the results of the settling column experiments by Model B, both flocculation models can be presented in the simplified general form F = A_f–B_f, where A_f and B_f represent the aggregation and break up terms, respectively. As indicated previously, smaller flocs with sizes approaching that of the primary particles (or small particle aggregates) have a denser structure (i.e., higher fractal dimension) and larger yield strength, making them more difficult to break up. For this particular condition, the turbulent stress μG is less than the floc yield strength τ_y, and the maximum fractal dimension nf_max [i.e., from μG=B1(Dd)2nfmax(D)/3D-2], is smaller than the value at which the flocs will break up. In other words, this indicates that the break up term B_f → 0 in Model B and, hence, the aggregation term will be dominant when floc sizes are small. This is the primary reason for the predicted rapid increase in floc size by Model B during the earlier stages of the runs.

Variations in the calculated floc sizes during slack water (Figure 11A) and ebb (Figure 11D) tidal phases indicate, firstly, a greater spatio-temporal variability in floc sizes is obtained with Model B than with Model A. This finding reflects (i) increased aggregation rates in Model B during slack water periods (i.e., floc sizes D up to ∼55 μm at z/h = 0.4, Figure 11A), compared to Model A (i.e., D up to ∼30 μm at z/h = 0.4, Figure 11A), and (ii) increased break up rates in Model B during high (ebb) tidal velocity periods (i.e., D ≈ 5–14 μm over z/h range, Figure 11D), compared to Model A (i.e., D ≈ 7–22 μm over z/h range, Figure 11D). This clearly suggests that the inclusion of variable fractal dimensions for all floc sizes (Model B) provides a more responsive flocculation model that better represents spatio-temporal floc evolution due to changing hydrodynamic conditions and SSC values (Figures 11C,F) within the tidal channel.

An interesting feature of the vertical distributions of floc sizes in Figure 11D (i.e., during the ebb tidal flow) is the general uniform and even a slight reduction in predicted floc size by both models from the water surface to the bed surface (i.e., z/h = 1 → 0). As, the largest turbulent shear rates are generated near bed, which result in smaller flocs, while the strong diffusion effects (ebb tidal flow) result in a more general uniform distribution. In this regard, the current model results for floc distributions within the water column are entirely consistent with previous observations by Guo et al. (2017) in the Yangtze river estuary. According to Guo et al. (2017), during the flood/ebb phase acceleration in tidal currents, the vertical distributions in measured mean floc sizes were relatively uniform (i.e., decreasing only slightly from the upper water column to the bed surface) and generally smaller than flocs generated under slack water conditions.

By contrast, Figure 11A shows that vertical floc size distributions are more variable during slack water (so as the settling velocity as shown in Figure 11B), with the largest floc sizes shown to occur at z/h = 0.4, and significant size reductions observed both in the water column above (approaching the water surface, z/h → 1) and below (approaching the bed surface, z/h → 0) this maximum. It is anticipated that this variable FSD occurs as larger flocs tend to settle out more quickly under quiescent conditions (i.e., during slack water) leaving only smaller flocs in the upper water column. The model results also show that average floc sizes are larger during the slack water than during the peak flood/ebb phases (and, thus the same with settling velocities as shown in Figures 11B,E), a trend that is again entirely consistent with the field measurements by Guo et al. (2017).

The floc sizes and effective floc densities are two key parameters that determine the settling velocities of flocs. Previous field studies in the Ems/Dollard estuary by Dyer et al. (2000) provided direct measurements of floc sizes and settling velocities during the flood phase of the tidal cycle (i.e., 2.13 and 0.14 h before HW). These can be compared, at least in a qualitative sense, with the current model predictions, albeit under different tidal conditions. Dyer et al. (2000) found that most smaller flocs measured in the estuary (d < 80 μm) had effective floc densities between 160 and 1600 kg.m^–3, with corresponding settling velocities between 0.01 and 1.0 mm.s^–1. By comparison, the calculated mean floc sizes within the Ems/Dollard estuary from the present modelling study (using Model B) during both peak flood/ebb phases and slack water periods are shown to typically vary between D = 10–60 μm, with effective densities between 160 and 1600 kg.m^–3 and settling velocities between 0.01 and 1.0 mm.s^–1. These values are therefore in broad agreement with the field measurements by Dyer et al. (2000) within the same estuary and provide further validation of the flocculation Model B with variable fractal dimensions.

Fettweis et al. (2006) also conducted field measurements of SSC, flow velocity, and floc size in the Belgian coastal zone and concluded that the Kolmogorov turbulent length scale was typically 3–10 times larger than the cohesive sediment flocs generated. Considering the field measurements from the Elms/Dollard estuary tidal channel used in the current study, the Model B predicted time series of average floc sizes at 0.4, 0.7, and 1.0 m elevations above the bed (i.e., equivalent to the elevations of the velocity measurements in Figures 6B–D) are shown in Figure 12, along with the predicted Kolmogorov scales at these elevations. It is shown that the predicted averaged floc sizes are generally significantly smaller than the Kolmogorov length scale, and only during periods of high SSC levels (i.e., on the flood and ebb phases, prior to slack water, Figures 8B,C) do we see significant floc growth (D ≈ 80–220 μm, Figure 12) at the three measurement elevations, which diminishes rapidly again at slack water, primarily due to floc settlement and the corresponding rapid reduction in SSC values (Figures 8B,C). Importantly, the corresponding Kolmogorov length scales at these elapsed times with high SSC values (and largest floc sizes) vary between about 400 and 720 μm, with the Kolmogorov length scale to peak floc size length ratio therefore varying between 3 and 5, in full accord with the findings of Fettweis et al. (2006). The current findings are also consistent with the assumption in the Winterwerp (1998) semi-empirical flocculation model (Eq. 19) that the Kolmogorov length represents the upper limit on the attainable equilibrium floc size generated under steady state conditions (i.e., constant concentration c and shear rate G).