Ocean DynamicsTheoretical, Computational and Observational Oceanography
© Springer-Verlag Berlin Heidelberg 2014

Modeling flocculation in a hypertidal estuary

Rafael Ramírez-Mendoza , Alejandro J. Souza  and Laurent O. Amoudry 
National Oceanography Centre, Joseph Proudman Building, 6 Brownlow Street, Liverpool, L3 5DA, UK
Rafael Ramírez-Mendoza (Corresponding author)
Alejandro J. Souza
Laurent O. Amoudry
Received: 20 March 2013Accepted: 29 November 2013Published online: 15 January 2014
Responsible Editor: Bob Chant
When fine particles are involved, cohesive properties of sediment can result in flocculation and significantly complicate sediment process studies. We combine data from field observations and state-of-the-art modeling to investigate and predict flocculation processes within a hypertidal estuary. The study site is the Welsh Channel located at the entrance of the Dee Estuary in Liverpool Bay. Field data consist of measurements from a fixed site deployment during 12–22 February 2008. Grain size, suspended sediment volume concentration, and current velocity were obtained hourly from moored instruments at 1.5 m above bed. Near-bottom water samples taken every hour from a research vessel are used to convert volume concentrations to mass concentrations for the moored measurements. We use the hydrodynamic model Proudman Oceanographic Laboratory Coastal Ocean Modelling System (POLCOMS) coupled with the turbulence model General Ocean Turbulence Model (GOTM) and a sediment module to obtain three-dimensional distributions of suspended particulate matter (SPM). Flocculation is identified by changes in grain size. Small flocs were found during flood and ebb periods—and correlate with strong currents—due to breakup, while coarse flocs were present during slack waters because of aggregation. A fractal number of 2.4 is found for the study site. Turbulent stresses and particle settling velocities are estimated and are found to be related via an exponential function. The result is a simple semiempirical formulation for the fall velocity of the particles solely depending on turbulent stresses. The formula is implemented in the full three-dimensional model to represent changes in particle size due to flocculation processes. Predictions from the model are in agreement with observations for both settling velocity and SPM. The SPM fortnight variability was reproduced by the model and the concentration peaks are almost in phase with those from field data.
Flocculation Modeling Liverpool Bay Hypertidal estuaries Suspended sediments Sediment transport
This article is part of the Topical Collection on Physics of Estuaries and Coastal Seas 2012

1 Introduction

Most estuaries are zones of significant human intervention. For example, many estuaries, ports, and navigation channels have to be maintained, which often lead to dredging and disposal of large amounts of sediment. Increased understanding on such anthropogenic forcing relies on a synergetic use of state-of-the-art, process-based, numerical models, and observations (Schuttelaars et al. 2012). Good predictive skill based on thorough understanding on the dynamical processes involved is necessary for informed decision making.
The fate of suspended material is not only important in terms of the overall sediment budget in estuaries, but also for primary productivity and pollution. Suspended particulate matter (SPM) can be divided between noncohesive and cohesive sediments depending on, amongst others, its mineral composition (Winterwerp and van Kesteren 2004). While good approximations of noncohesive sediment dynamics can be achieved (e.g., van Rijn 2007a), predictions of cohesive sediment behavior present significant challenges because of numerous processes varying both in time and in space (Winterwerp and van Kesteren 2004). In particular, a critical process in cohesive sediment studies is flocculation. In this study, this term is used to describe the combination of aggregation and breakup of flocs. This process, or combination of processes, has been observed and found to be important in a number of coastal areas and estuaries. Fugate and Friedrichs (2003) reported SPM concentrations of about 300 mg l−1 and intratidal grain size changes in the range of 10–300 μm in Chesapeake Bay, USA. On the Belgian coast, flocculation leads to particle sizes varying in the range of 10–120 μm between flood and ebb cycles with concentrations of 100 mg l−1(Fettweis et al. 2006).
Modeling flocculation and the induced variations in settling velocity is critical towards reproducing the dynamics of cohesive suspended sediment. Winterwerp et al. (2006) used a one-dimensional model coupled to a formulation for flocculation and compared the results against observations over one tidal period. In a formulation of van Rijn (2007b), the flocculation effect is included as a factor in the calculation of the suspended sediment fall velocity. More recently, Kombiadou and Krestenitis (2012) have used a particle tracking model to give three-dimensional results and compared numerical predictions to field data under very dilute conditions (concentrations lower than 7 mg l−1). Nevertheless, modeling studies that consider the flocculation process in three-dimensional baroclinic estuarine models are still scarce or are not fully tested (e.g., Amoudry and Souza 2011a). One of the key variables related to the flocculation process is the floc settling velocity (W s ) (Fugate and Friedrichs 2003), and it is often found to vary significantly in estuaries at intratidal time scales (e.g., van Leussen 1994; Winterwerp et al. 2006; Manning 2008; Yuan et al. 2008). The present study seeks to improve long-term predictions of SPM behavior by accounting for flocculation via the implementation of a semiempirical formulation for the settling of flocs in a three-dimensional baroclinic coastal ocean model.
The study site is the Dee Estuary, located in the Liverpool Bay (Fig. 1) in the eastern Irish Sea. The Dee Estuary is an energetic estuary with a tidal range reaching about 10 m and tidal currents reaching 1.2 m s−1 at spring tides. Two deep channels are present near the mouth: the Hilbre Channel in the east and the Welsh Channel in the west. Figure 2 shows near-bottom current magnitude for the Welsh channel and river discharge from the Dee. The overall estuarine circulation pattern is vertically sheared in the Hilbre Channel and horizontally sheared in the Welsh Channel (Bolaños et al. 2013), both with little influence of stormy periods on the respective pattern (Brown et al. 2014). The Dee Estuary consists of mixed sediments and cohesive behavior has been observed (Thurston 2009). We focus on combining state-of-the-art modeling with field observations for a calm period in February 2008 and only for the Welsh Channel. Comparison between observations and modeled suspended sediment concentration for the Hilbre Channel are provided by Amoudry et al. (2014).
Fig. 1
Site study location. a United Kingdom. Red square marks the location of the Liverpool Bay. b Liverpool Bay and, in red square, the location of the Dee Estuary. c Bathymetry in meters in the entrance to the Dee Estuary and Welsh and Hilbre Channels
Fig. 2
a Current speed measured in the Welsh Channel at 0.35 m above the sea bed. b Fresh water river discharge to the Dee Estuary
We present in the next Section 2 the description of the measurements used in this investigation followed by some characteristics of the numerical models and the parameters used in the simulations. Results are presented in Section 3 starting with SPM observations suggesting flocculation of particles at the study location. We also derive a semiempirical formulation for the settling velocity which is then implemented in the model, and finally a comparison between observations and numerical results is presented. In the discussion in Section 4, we comment on the dynamics of fine sediments, their description from a fractal point of view, the importance of the floc effective density, the role of the settling velocity as proposed in this investigation and compared with a widely used model, and the relationship between hydrodynamics and SPM in the Dee Estuary.

2 Methodology

We employ a modeling suite which couples hydrodynamic and turbulence models and which has been applied to a variety of shelf sea domains (e.g., Holt and Umlauf 2008; Polton et al. 2011; Bolaños et al. 2013). A sediment module is included in the hydrodynamic model (Amoudry and Souza 2011b; Amoudry et al. 2014). Field measurements were obtained in the Dee Estuary from the deployment of instruments near the bed and samples taken from a research vessel (Bolaños and Souza 2010). The samples were used to calibrate the instruments and lead to a formulation for the flocculation process which in turn was implemented in the sediment module.

2.1 Observations

A mooring was deployed in the Welsh Channel (Fig. 1) from 12 February to 13 March 2008 and was equipped with a Laser In Situ Scattering and Transmissometry (LISST) and an acoustic Doppler velocimeter (ADV) amongst other instruments (for a more detailed description of the instrument deployment, see Bolanõs and Souza 2010). From the LISST, volume concentrations for each 32 grain sizes from 2.7 to 460.1 μm were obtained. The operational principles of the instrument are given by Agrawal and Pottsmit (2000). The LISST recorded one sample every 40 s during a 20-min period every hour which were then averaged to obtain hourly data. Summation of volume concentrations from all grain sizes produces a time series of the total volume concentration. Current velocity and pressure were measured with the ADV at 16 Hz during 20 min every hour. The fast sampling rate allows the calculation of accurate turbulent stresses (τ) using
$$ \tau=\rho_{w} \sqrt{{\overline{u'w'}}^{2}+{\overline{v'w'}}^{2}} $$
where ρ w is fluid density, and u , v , and w are velocity fluctuations. Both instruments recorded from 12 February to 13 March, 2008. However, wave analysis from the ADV showed the presence of stormy waves from 22 February onwards (Bolaños and Souza 2010), and the analysis was stopped at this date to avoid the effect of waves. This study only intends to investigate flocculation processes under current dominated conditions. Although a final hourly data sampling rate may seem to be low, it still allows to investigate processes occurring at a quarter-diurnal frequency.
Near-bottom water samples using a rossete mounted in a CTD were taken every hour in the Welsh Channel on 12–13 February 2008, and SPM mass concentrations were obtained from these samples using gravimetric technique (e.g., Jones et al. 1998). The relationship between mass concentration from water samples and volume concentration is widely used because of the capability of the LISST to capture the total volumetric concentration of suspended particles (e.g., Fugate and Friedrichs 2002; Badewien et al. 2009). In fact, Fugate and Friedrichs (2002) tried to improve the calibration, taking into account the variability of the grain size with no significant results. Therefore, concentration behavior between volumetric concentration C v and mass concentration C m seems to be qualitatively similar. The regression analysis between volume and mass concentration is shown in Fig. 3. Two data points are clearly outside the trend (triangles in Fig. 3) and may be the result of several factors: measurement errors by the instrument, water sample handling, occasional events at time of measurements, and sample taking. These data were not taken into account in the analysis. The minimum values were 354.1 μl l−1 corresponding to 18.5 g m−3, while maximum values were 1,887.4 μl l−1 with 138.1 g m−3. The resulting linear relationship is
$$ C_{\mathrm{m}}=0.07C_{\mathrm{v}}-0.65 $$
with a coefficient of determination r 2 = 0. 88 and 14 degrees of freedom. This relationship is used hereafter to convert the LISST volume concentrations into mass concentrations.
Fig. 3
LISST Volume concentrations (C v) and mass concentrations (C m) from water samples (blue points). Data outside the trend, marked also with triangles, are not included in the regression analysis. The line is the result of the linear relationship between the two variables

2.2 Computation of settling velocities

If most of suspended particles are flocs, the effective density ρ e is related to mass and volume concentrations according to the following relation (Mikkelsen and Pejrup 2001):
$$ \rho_{\mathrm{e}}=\frac{C_{\mathrm{m}}}{C_{\mathrm{v}}} $$
Therefore, following from the regression analysis in the previous section (Eq. 2), ρ e can be obtained from C m , and C v for each LISST measurement.In addition, a time series of median grain sizes (D 50) can be calculated from the LISST data for 32 grain size classes. Using the time series for ρ e and D 50, settling velocities (W s ) are calculated for the period 12–22 February, using the Stokes law as a first approximation.
$$ W_{\mathrm{s}}=\frac{g\rho_{\mathrm{e}}{D_{50}}^{2}}{18\mu} $$
where g and μ respectively, are the gravity acceleration and the dynamic viscosity. The validity of this assumption will be discussed when results are presented. Finally, a regression analysis of turbulent stresses against settling velocities is used to find a relation between these two quantities at the study site. The resulting data distribution leads to the formulation implemented in the sediment module (Fig. 5).

2.3 Numerical models

The hydrodynamic model Proudman Oceanographic Laboratory Coastal Ocean Modelling System (POLCOMS) solves the incompressible, hydrostatic, Boussinesq equations of motion. It uses an Arakawa B-grid and sigma coordinates for the horizontal and vertical discretizations, respectively, and employs a time-splitting technique to obtain barotropic and baroclinic components. In the Arakawa B-grid, elevation and scalar quantities are evaluated at each grid cell corner, while velocity and flux components are evaluated at the center of each grid cell (Arakawa 1972). The equations for depth-averaged currents and free surface elevation are solved following a forward-time centered-space technique, while advective terms use a piecewise parabolic method (more details by Holt and James 2001). Delhez et al. (2004) found similar results of POLCOMS and several other models in a comparison study. The good performance of the model applied to the Dee Estuary has been reported by Bolaños et al. (2013). The turbulence model GOTM—General Ocean Turbulence Model (Burchard et al. 1999; www.​gotm.​net)—is used with a kε closure scheme. Holt and Umlauf (2008) carried out a comparison of the use of coupled POLCOMS-GOTM with observations. The sediment transport module included in POLCOMS (e.g., Amoudry and Souza 2011b; Amoudry et al. 2014) consists of an advection-diffusion equation for suspended sediment. The critical erosion bed shear stress and erosion constant are user-defined, and for this study, the sediment is eroded from the bed. The effect of different turbulent closure schemes on sediment transport using POLCOMS, GOTM, and the sediment module for idealized scenarios was carried out by Amoudry and Souza (2011b). However, sediment modeling in POLCOMS, as well as in other coastal ocean modeling systems, has typically been restricted to constant settling velocities. In the present study, we implement a variable settling velocity (varying in time and space): W s is empirically related to the turbulent stress.

2.4 Numerical setup

For this study, sigma coordinates with ten equidistant vertical levels and 180-m horizontal grid size are used. Higher vertical resolution was tested, but the numerical results were not significantly modified. The overall modeling system is forced at the free surface by meteorological data obtained from a mesoscale model of the National Weather Service of the UK. Values for temperature, salinity, elevations, and barotropic tidal currents are specified hourly at lateral open boundaries. Values for the stability functions in the closure scheme of GOTM are 0.52 and 0.73, respectively. In the sediment module, the erosion rate was set to 1. 25 × 10−5 kg m−2 s−2 with a critical stress value of 0.18 Pa. Field observations show the presence of waves from 22 February 2008 in the Welsh Channel as mentioned in Section 2.1. Therefore, the models are set to run from January and February 2008 with January used to spin up hydrodynamics and turbulence. The first 10 days of February are used to spin up the sediment module.

3 Results

3.1 Suspended sediment observations

Figure 4 shows the time series of LISST observations. Figure 4a shows the changes in median grain size in relation with tidal phase. The change in tidal amplitude from neap to spring tides is shown by the tidal amplitude increasing from about 4 m to about 8 m. Minimum values of D 50 of about 50 μm are found at the end of the record during spring tides, while maximums of more than 200 μm were commonly found during neap tides. Similar to the tide, a fortnightly modulation seems to be present in D 50, but, conversely to the tide, the amplitude is large during neap tides and small during spring tides. Quarter-diurnal frequency cycles are also present in D 50 with maximums at slack water and minimums during maximum flood and ebb. During neap tides, higher values were achieved after ebb than after flood. Mass concentration for each grain size is plotted in Fig. 4b in logarithmic scale. The maximum observed concentration value was about 9 g m−3 for a sediment of 122- μm size and occurred on 21 February at 1600 hours. High concentrations of sediment of size about 100 μm were present during flood and ebb periods. High concentrations of coarse particles, about 350 μm, were recorded during slack waters, although magnitudes were lower than for flood and ebb. This behavior was observed even during neap tides when suspended sediment concentrations were overall smaller.
Fig. 4
LISST observations. a Depth (blue line) and median grain size D50 (green line). b Contours of mass calibrated concentrations in logarithmic scale (log10 C m for each grain size during the study period

3.2 Settling velocities

In Fig. 5, turbulent stresses (τ) obtained from ADV and settling velocities calculated using Stokes law are presented together with the adjusted curve. The maximum particle Reynolds number is on the order of 10−1, and our assumption of using the Stokes law is reasonable. Some low values of turbulent stresses correspond to a wide scatter for settling velocities (higher than 1.5 mm s−1) which seem to be outside the trend and were not taken into account in the study. The maximum value of turbulent stress is about 5 Pa and corresponds to a minimum settling velocity of 0.1 mm s−1. The minimum value of turbulent stress of about 0.025 Pa coincides with settling velocity about 0.8 mm s−1. However, there is a wide range of settling velocities for a fixed value of turbulent stress. The range of settling velocity increases as turbulent stresses diminishes from τ ∼ 5 to 0.3 Pa, while the scattering of W s covers most of the range from 0.2 to more than 1.5 mm s−1.
Fig. 5
Turbulent stresses from ADV and settling velocities obtained with Stokes using sediment samples and LISST data. The solid line is the adjusted curve to the exponential formulation shown at the upper right corner. Data with W s higher than 1.5 were not taken into account
We aim to represent the observed trend by an expression that can easily be implemented into our three-dimensional model. To that end, we seek a formula that will satisfy the following mathematical constraints: finite and nonzero settling velocities as stress becomes either infinitely large or small, continuity, and derivability. The simplest mathematical expression that will match these criteria follows the following form:
$$ W_{\mathrm{s}} = W_{\mathrm{s},0} + W_{\mathrm{s,m}} \: e^{- \tau / \tau_{cf}} $$
where W s, 0 represents the settling velocity of smallest flocs, W s, m is such that W s, 0 + W s, m is the maximum settling corresponding to minimum stress, and τ c f a critical stress value denoting how quickly the observed settling rates decrease with increasing turbulent stress. Further discussion on this expression is provided in Section 4.4. Fitting the above expression to the observations leads to the following values: W s, 0 = 0. 25 mm s−1; W s, m = 0. 89 mm s−1; and τ c f = 0. 19 Pa. The maximum settling velocity is then 1.14 mm s−1.

3.3 Model observations comparison

A comparison between model results and observations is presented in Fig. 6. The water depth (Fig. 6a) shows that the study period mostly covers neap tide conditions. Tidal amplitude during the study period approximately ranges from 5 to 9 m. Differences between model and observations are small over the entire period, and the semidiurnal tide is well predicted. The hydrodynamics of the Dee Estuary are also reasonably well described as shown by Bolaños et al. (2013) and Amoudry et al. (2014).
Fig. 6
Comparison between observations (blue points and lines) and model results (red lines). a Sea level, observations from LISST pressure sensor. b Turbulent stresses from ADV current recordings and modeled values. c Settling velocities computed with Stokes law using data from LISST and filtered samples. Modeled settling velocities obtained using the semiempirical formulation. d Total observed mass concentrations from LISST and predicted concentrations including flocculation effect
The turbulent stresses are qualitatively well predicted by the model (Fig. 6b) although quantitatively under-predicted; this is discussed by Bolaños et al. (2013). As expected, the neap to spring evolution, the quarter-diurnal variability with maximums during flood and ebb periods and the corresponding asymmetry (i.e., high stress during ebb and low stress during flood) are all reproduced by the model. Quantitatively, the model performs well in most of flood periods, although for ebbs, the predicted magnitude is lower than observations with bigger relative differences during neap tides. Maximum turbulent stress found in observations is 5.1 Pa, while the model maximum is 2.8 Pa, both occurring during spring tides at 1400 hours on day 21. The determination coefficient between calculated and predicted turbulent stresses is r 2 = 0. 75 with 95 % of confidence interval.
As expected, settling velocities obtained with the model are in good agreement with observations, which is a direct consequence of prescribing the settling rate in the model following Eq. 5. Nevertheless, there are still some discrepancies on the high and low values. Observed settling velocities (Fig. 6c) show large intratidal variations during neap tides (15 to 18 February) with values ranging from approximately 0.25 to 2.1 mm s−1. As the tidal range increases after 18 February, these intratidal variations in settling rates are reduced and W s remains between 0.09 and 0.3 mm s−1. In comparison, the modeled settling rates are restricted by Eq. 5, and the resulting values range from 0.25 to 1.14 mm s−1. Modeled fall velocities also have a lower range during spring tides, therefore reproducing the observed trend. The quarter-diurnal variations are qualitatively well captured by the model, even though maximum settling values are underestimated around neap, and asymmetry in settling rate between high water and low water is not reproduced. It is important to note that maximum settling velocities, both modeled and observed, coincide with slack water and minimum stress values. Such discrepancy on the maximum settling rates would therefore not be caused by the disagreement in maximum stress values (Fig. 6b). Maximum observed settling velocity of 2.12 mm s−1 was found at 1400 hours in day 17, while model maximum is 1.12 at 1330 hours in day 16, both maximums during neap tides. Minimum observed value of ∼ 0.09 mm s−1 was found at 1500 hours in day 21, whilst model minimum is 0.25 mm s−1 at 1400 hours in day 21, both during spring tides.
The model reproduces reasonably well the general variability of SPM concentration during the study period (Fig. 6d). Lowest concentrations are found during neaps and highest concentrations in spring tide. However, during most of the first 3 days of the record, the model underestimates observed concentrations. Maximum values are 123.2 and 148.5 g m−3 for the observed and model concentrations occurring on the last ebb cycle in day 21. The quarter-diurnal variability is captured by the model for the entire record, and tidal asymmetry is reproduced for most of the period. During the first 3 days, a tidal asymmetry is obtained in the model, but the maximums are found during flood, while they are found during ebb cycles in the observations. From the end of day 17, observed ebb peaks are higher than flood, and the model qualitatively reproduces this difference of maximum concentrations even though it overestimates its magnitude. There is a timing difference between the maximum flood concentrations for which the model values are about 1 h late.

4 Discussion of cohesive sediment dynamics in the Dee Estuary

4.1 Floc size variability

In this work, we use the changes in grain size to identify the flocculation process at the study site. In particular, the changes in grain size showed in Fig. 4 are useful to see the effect of currents on flocs. There are differences in floc maxima and minima between spring and neap tides. During neaps, the turbulent stresses are lower, and the breakup process is not as strong as during spring tides. As a consequence, the minimum median grain size is about 75 μm around neap tides, while it is 50 μm around spring tides. This spring-neap grain size variability has also been found by several authors (e.g., Fettweis et al. 2006; Bartholomä et al. 2009). Another feature is a difference in the range of grain sizes of ∼ 160 μm during neaps and ∼ 50 μm during springs. Therefore, the aggregation process is more important during neaps.
There is, in particular during neap tides, an asymmetry in maximum floc sizes between low and high water: bigger flocs are found at low water. We already mentioned that this is not due to asymmetry in turbulent stress values between ebb and flood as these two asymmetries occur at different phases of the tide. Instead, a possible explanation could be that flocs of different origins are advected during ebb and flood: during ebb, fine particles would be advected from the salt marshes or the upper estuary (e.g., Amoudry et al. 2014), while during flood, the observed flocs would be of marine provenance. The different origins of particles may result in different settling velocities for similar low stress values due to different sediment characteristics. For increased tidal range (i.e., not during neap tides), the magnitude of turbulent stress is higher (Fig. 6b), which leads to reduced asymmetry in W s because there is not enough time with low enough turbulent stresses for the sediments to aggregate.Under strong currents at flood and ebb, the bed shear stress increases, and sediment is taken into suspension. Small flocs are resuspended, and coarse flocs are subjected to breakup, both of which increase fine sediment concentration, as seen in Fig. 4. The presence of coarse sediment in low concentrations during flood and ebb may be the result of the same resuspension process or flocculation by interparticle collision due to high concentration of suspended sediments. At slack water, the turbulent stresses diminishes, and the aggregation process dominates and increases the concentration of big flocs. This behavior has been found in several places with similar hydrodynamic characteristics (e.g., Jago et al. 2006; Badewien et al. 2009; Bartholomä 2009).
Badewien et al. (2009) found changes in grain size from approx 20 to 400 μm, and Jago et al. (2006) found approx 50 to 150 μm of median grain size. Krivtsov et al. (2008) measured grain size in Liverpool Bay and found highest concentrations of particles bigger than 200 μm which is in agreement with the size of big flocs encountered in this study. These features have been found by Jago et al. (2007) and Krivtsov et al. (2008) in sites near the area of this study. Flocs are bigger on neap tides due to longer residence time in the water column and low turbulent stresses. On the other hand, during spring tides, turbulent stresses prevents the formation of big flocs.

4.2 Treatment of flocs as fractal structures

The treatment of flocs as fractal structures is now well established and relatively common both in modeling and observational studies (e.g., Kranenburg 1994; Winterwerp 1998; Jago et al. 2006; Verney et al. 2011; Kombiadou and Krestenitis 2012). The basis is that flocs have a complex structure, which is related to their density and yield strength, and which can be modeled using fractal theory applied to cohesive sediments. This establishes that primary particles and flocs are related by a fractal dimension or number in the range of 1–3 (Kranenburg 1994). Following from experiments, a fractal dimension close to one corresponds to very fragile flocs and a dimension of three to compact aggregates (e.g., Kranenburg 1994; Winterwerp 1998). Field measurements have shown that different values could be found, and the difficulty to obtain observational data to calculate the real fractal number has led some authors to assume a number of two (e.g., Jago et al. 2006; Kombiadou and Krestenitis 2012). However, precise determination of the floc fractal number is important because it also provides an indication of the strength of cohesive forces for the flocs of the study site and their tendency to aggregate or break up.
According to the fractal theory applied to flocs (Kranenburg 1994), mass concentration is related to volume concentration and particle size as follows:
$$ C_{\mathrm{m}}=\rho_{\mathrm{s}} C_{\mathrm{v}}\left(\frac{D_{\mathrm{p}}}{D_{50}} \right)^{3-n_{\mathrm{f}}} $$
where ρ s is sediment density, D p is primary particle diameter, and n f is the fractal number. This equation can be rearranged into a linear relationship between known quantities, i.e., mass and volume concentrations and median grain size:
$$ C_{\mathrm{m}}=\rho_{\mathrm{s}} {D_{\mathrm{p}}}^{3-n_{\mathrm{f}}} \left(\frac{C_{\mathrm{v}}}{D_{50}^{3-n_{\mathrm{f}}}}\right) $$
which can then be used to estimate the n f value for the study site. To that end, Eq. 7 is used with different values of n f , assuming here that the fractal number remains constant during the study period. The best linear relationship with a determination coefficient r 2 of 0.9 was found for a fractal number of 2.4 and is shown in Fig. 7. Units of mass concentration are showed in kilograms per cubic meter to be consistent with units of density. The corresponding value for the factor \(\rho _{\mathrm {s}}\cdot {D_{\mathrm {p}}}^{3-n_{\mathrm {f}}}\) is 0.0008. The floc fractal number of 2.4 obtained from observations in the present study is in the range expected (1–3).
Fig. 7
Regression to find the fractal number for the study site based on measured quantities

4.3 Floc effective density

In the present study, the values of effective density used in Eq. 4 are implicitly inferred by the conversion of LISST volume concentrations into mass concentrations. These results are in the range reported by Mikkelsen and Pejrup (2001) for different sites. The effective density has been related to the floc size via a range of formulations (e.g., McCave 1984; Gibbs 1985; Khelifa and Hill 2006; Maggi 2007). However, an important aspect of our observations is the very small dynamical range for both D 50 and ρ e : D 50 varies between 50 and 200 μm, and ρ e varies between 46 and 67 kg m−3. This makes it very difficult to compare our results with expressions which have been derived based on several data sets covering a very wide dynamical range for both D 50 and ρ e . Nevertheless, it has to be noted that the fractal number calculated in the previous Section (n f = 2. 4) implies that we have ρ e ∝ (D 50)−0.6, which is qualitatively similar to other expressions such as those of McCave (1984), i.e., ρ e proportional to the negative power of floc size.

4.4 Parameterization of flocculation via settling velocity formulation

The settling velocities obtained from Eq. 4 are in general agreement with values reported by several authors. For example, Winterwerp et al. (2006) and Manning et al. (2010) found W s between 0.1 and 2 mm s−1 for grain sizes between 50 and 200 μm, while van Rijn (2005) reported on settling velocities between 0.06 and 1 mm s−1 for the Mersey Estuary, which is a nearby estuary (approximately 10 to 20 km to east of the Dee Estuary) with similar environmental conditions.
A relationship between grain size and shear stress has conceptually been described by Dyer (1989). For low shear stress, aggregation occurs due to interparticle collisions until a threshold size is achieved. This results in floc size increasing with low shear stress until a maximum size is achieved, and breakup of flocs becomes more important, and floc size then decreases with increasing shear stress. This behavior has been found for turbulent stresses and floc settling velocity (instead of shear stress and grain size) using a settling chamber (e.g., Winterwerp et al. 2006; Manning 2008). Winterwerp et al. (2006), numerically, found maximum or critical settling velocity at intermediate values of turbulent stress. However, if the column water increases, higher values of the maximum settling velocity are found at lower values of turbulent stress. At turbulent stress values very close to zero, there is no maximum settling velocity as its value would be infinite (equilibrium settling velocity).
In the present investigation, we found from observations that the highest settling velocities (or largest flocs) correspond to the lowest values of turbulent stress. However, the big flocs are the result of aggregation as turbulent stress is diminishing rather increasing, i.e., the grain size changes from microflocs to macroflocs and not from primary particles. In addition, the water depth of about 14 m enables the floc growth to achieve an equilibrium settling velocity (Winterwerp et al. 2006). These results have also been found in experiments with mud (e.g., van Leussen 1994; Verney et al. 2011). Another feature found in this study at low values of turbulent stresses is the high variability of settling velocities which has also been reported in experiments with mud by Winterwerp (1998) and by Verney et al. (2011) for high variability of floc median diameter at low values of shear rate. This may be due to a high sensitivity of W s to the environment conditions (Winterwerp 1998) or to a longer aggregation time because of the low collision rates (Verney et al. 2011). The observations also showed settling velocities decreasing as turbulent stresses increase which coincides with the behavior described by Dyer (1989) after the grain size threshold at low values of shear rate. In fact, the relationship between shear rate and maximum settling velocities found by van Leussen (1994) is in qualitatively good agreement with the semiempirical formulation for W s proposed in this study.
Several relations have been proposed to obtain the settling velocities of flocs with different degrees of complexity and variables involved since the simple power relation with particle diameter to the Stokes law modification which include the fractal dimension proposed by Winterwerp (1998). However, as formulations become more complete and dependent on multiple variables, requirements on field data increase and the implementation in a numerical model becomes more complex. The flocculation process is affected by turbulent stresses, differential settling, Brownian motion, concentration, salinity, and biology factors (Fugate and Friedrichs 2003). For estuaries and coastal zones, Brownian motion and differential settling are negligible (Winterwerp 2002). Salinity is important if there are values lower than three PSU which would be at the river-estuary transition zone. In this study, we assume that the biology effects are negligible in winter, and we take flocculation dependent on turbulent stresses and concentration only. In the Dee Estuary, turbulence plays an important role in the changes of grain size (see Fig. 5 in Amoudry and Souza 2011a) and therefore in the settling velocity. Thurston (2009) found strong relationships of grain size with the Kolmogorov microscale, which can be related to the turbulent stress itself. The proposed fall velocity formula in this investigation is based on the behavior of observed quantities and depends only on turbulent stresses computed simultaneously with suspended sediment calculations. The formulation describes high settling velocities or big flocs formed because of the aggregation of smaller flocs during periods of low turbulent stresses or slack waters. As current increases, turbulent stresses also increase, and the breakup process begins to disrupt the flocs and their size diminishes. If turbulent stresses continue increase, flocs continue breaking apart until a minimum size is achieved. Also, the formulation has mathematical advantages, i.e., near-zero values of turbulent stress result in a settling velocity close to a defined value ( 0. 25 + 0. 89) and very high turbulent stresses result in a constant value of 0.25. These features also have physical meaning because flocs cannot either aggregate to an indefinite size or be breaking apart beyond a minimum size. We thus consider that the formula is suitable as a first approximation to model the flocculation process.
In the present study, settling velocities are calculated following Eq. 4, which implies (a) validity of the Stokes law (briefly discussed in Section 3.2) and (b) knowledge on both ρ e and D 50. In our case, D 50 is directly obtained from the LISST measurements and ρ e inferred by the calibration of LISST volume concentration into mass concentration. Other approaches could be employed for obtaining the effective density. As mentioned in the previous section, a number of expressions can be used to relate the effective density to the floc size and could therefore be used to calculate the settling velocity from our observational data. Comparisons with three models (McCave 1984; Winterwerp 1998; Khelifa and Hill 2006) have been carried out. Reasonable settling velocity values were only obtained following Winterwerp (1998), and only this case is discussed further. According to fractal theory, effective density is as follows:
$$ \rho_{\mathrm{e}} = (\rho_{\mathrm{s}}-\rho_{\mathrm{w}}) \left( \frac{D_{\mathrm{p}}}{D_{50}} \right)^{3-n_{\mathrm{f}}} $$
$$ \frac{\rho_{\mathrm e}}{{D_{50}}^{n_{\mathrm f}-3}}=(\rho_{\mathrm s}-\rho_{\mathrm w}){D_{\mathrm p}}^{3-n_{\mathrm{f}}} $$
where the left-hand side can be calculated using our data calibration set of filtered water samples for C m , corresponding C v , and D 50 from LISST and n f = 2. 4 found in Section 4.2. Averaging all samples, we found
$$ (\rho_{\mathrm{s}}-\rho_{\mathrm{w}}){D_{\mathrm{p}}}^{3-n_{\mathrm{f}}}=0.2662 $$
From Winterwerp (1998), the equation for W s can be written as follows:
$$ W_{\mathrm{s}} = (\rho_{\mathrm{s}}-\rho_{\mathrm{w}}){D_{\mathrm{p}}}^{3-n_{\mathrm{f}}} {D_{50}}^{1.4}\left[\frac{g}{18\mu(1+0.15Re^{0.687})}\right] $$
$$ W_{\mathrm{s}}=0.2662{D_{50}}^{1.4}\left[\frac{g}{18\mu(1+0.15Re^{0.687})}\right] $$
where α = β = 1 as flocs are assumed to have spherical shape. The comparison between settling velocities using Eqs. 4 and 12 is shown in Fig. 8. As expected, results are similar because Eq. 12 is a modification of the Stoke’s law. The largest differences are found for the highest values of W s which is due to the highest variability of the big flocs at low turbulent stresses.
Fig. 8
Comparison of models to obtain floc settling velocity based on the observations of this study. Red line denotes using Stokes law taking into account that Re < 1 for this investigation (i.e., using Eq. 4). Blue line denotes using that of Winterwerp (1998) in which effective density and median grain size avoid the assumptions of primary particle size and effective density (i.e., using Eq. 12)

4.5 Coupled sediments and hydrodynamics in the Dee Estuary

Turbulent stresses are of similar order of magnitude as some previously reported (e.g., Winterwerp et al. 2006; Baugh and Manning 2007). Based on quantities obtained from observations, the turbulent stress asymmetries have a clear relationship with corresponding minimums of settling velocities. Higher values of τ during ebb lead to a more efficient breakup, and smaller flocs are present with slightly lower settling velocities than during flood period. During neap tides, the model reproduces very well this feature with similar minimums of calculated and modeled W s values. During spring tides, maximum turbulent stresses are well predicted, but corresponding minimums of modeled settling velocities are limited by the formulation. This result implies that some processes may either not well reproduced or may not be taken into account in the model. Calculated settling velocity maxima coincided with low values of turbulent stresses, and the highest W s magnitudes happened during neap tides and seems to be related to the asymmetries between previous flood or ebb periods. This could be due to an enhancement of floc growth in successive periods of low turbulent stress and more availability of suspended sediment after ebb by either advection from the estuary or local resuspension. Modeled settling velocities have maximum values with minimum turbulent stresses and highest W s magnitudes during neap tides. However, there is no clear asymmetry as in observed data.
The stress asymmetries are due to barotropic and baroclinic processes. However, processes like strain-induced periodic stratification and internal tidal asymmetry (tidal straining) can play an important role in the estuarine circulation (Burchard and Hetland 2010). Giddings et al. (2011)found the variability of advection and straining during flood and ebb strongly related to stratification, mixing, and longitudinal dispersion. According to Simpson et al. (2005), tidal straining has been found to control the structure and intensity of turbulent stresses on both ebb and flood with higher values in the lower half of the water column. Also, stress asymmetries are also related with the sediment dynamics.The internal tidal asymmetry has been mentioned as responsible of the estuarine turbidity maximum in the Columbia River (Jay and Musiak 1994). Burchard et al. (2008) found an important contribution of tidal straining to net sediment transport in the Wadden Sea. Tidal current asymmetry with ebb dominance and no river discharge has been reported as the main factor for sediment transport (Bartholomä et al. 2009). In the Dee Estuary, there is a contribution from the river, higher than that of Simpson et al. (2005), and strong tidal currents. This seems to enhance stratification and enable asymmetries in turbulence production (Bolaños et al. 2013). Nevertheless, Bolaños et al. (2013) show that the POLCOMS-GOTM model predicts accurately barotropic processes, but not the baroclinic part, which leads to an underestimation of the current tidal asymmetries.
There are time lags between observed and modeled results. For turbulent stresses (Fig. 6b), the time difference is mainly during flood periods with modeled values delayed with respect to observations. Modeled settling velocities lag observations by about 1 h throughout the study period for both ebbs and floods (Fig. 6c). Modeled suspended sediment concentrations display a similar lag mostly during floods around neap (Fig. 6d). These results could be related to the time lags of the hydrodynamic model which can be seen in a comparison of modeled and observed density anomaly in Fig. 9. Observations show stratification at the end of ebb periods (stronger during the first), while modeled stratification is weak and the response is delayed. Then, stratification seems to affect the timing of predicted settling velocities and concentrations observed mainly with comparison at maximum values against equivalent variables computed from observations.
Fig. 9
Welsh Channel model and observations during two tidal periods. a Measured sea level. b Vertical distribution of observed density anomaly (kilograms per cubic meter). c Vertical distribution of modeled density anomaly (kilograms per cubic meter). Vertical axes in b and c are in \(\sigma \)-coordinates used in the models
Another difference in predicted and calculated W s and SPM is the magnitudes during spring tides. In this cycle, the floc concentration may be an important factor, while during neap cycles, the timing factor is most significant. Also, settling velocities in the second half of the record have a flat shape when turbulent stresses have minimum values. This is a result of the minimum limit established by the settling velocity formulation, and the comparison suggests a different minimum value. This factors can be taken into account to improve modeled predictions and, in general, the performance of the models (POLCOMS, GOTM, and the sediment module accounting for flocculation) is good for the time scale at which were applied with also the three-dimensional results.

5 Conclusions

A simple, semiempirical formulation for the settling velocity of flocs is proposed. The formula has the following features: (1) It is based on variables obtained directly from observations with instruments widely used and relatively easy to deploy in the field, (2) few intrinsic assumptions are made, mainly that suspended particles are flocs and that volume concentration measurements are for spheric equivalent diameter particles, (3) it depends solely on turbulent stress as the main factor for the flocculation process, (4) it has both physical and mathematical meaning, and (5) it is easy to implement in a numerical model because it depends only on a model-based variable τ. In the present study, the semiempirical formulation was implemented in a sediment module and state-of-the-art three-dimensional hydrodynamic and turbulence models. Observations showed aggregation of flocs during slack water and breakup during flood and ebb periods. Flocculation was also found as a result of advection of fine material during ebb periods followed by flocculation at slack water. Predicted settling velocities and concentrations are in reasonable agreement with calculated values from observations. However, the model fails to predict quantitatively the asymmetry in turbulent stresses and settling velocities. Modeled SPM concentrations follow the observed behavior with differences in timing and maximum values. The time difference is due to a lag on the salt wedge in the model results. These differences suggest that factors like sediment consolidation and turbulence advection may play an important role for the sediment transport in the Welsh Channel during this calm winter period.
We want to acknowledge the organizations of Mexico and the UK, which provided the funding and opportunity for this investigation to be carried out: Consejo Nacional de Ciencia y Tecnología (CONACyT, through the PhD scholarship ID 212026), Centro de Investigación Científica y de Educación Superior de Ensenada (CICESE), National Environmental Research Council (NERC, via the iCOASST project, grant NE/J005444/1), and the National Oceanography Centre (NOC). We also want to thank the anonymous reviewers for their important comments to improve this manuscript.
Agrawal Y, Pottsmit H (2000) Instruments for particle size and settling velocity observations in sediment transport. Mar Geol 168:89–114CrossRef
Amoudry LO, Souza AJ (2011a) Deterministic coastal morphological and sediment transport modeling: a review and discussion. Rev Geophys 49:1–21. rG2002. doi:10.​1029/​2010RG000341 CrossRef
Amoudry LO, Souza AJ (2011b) Impact of sediment induced stratification and turbulence closures on sediment transport and morphological modelling. Cont Shelf Res 31:912–928CrossRef
Amoudry LO, Ramírez-Mendoza R, Souza AJ, Brown JM (2014) Modelling-based assessment of suspended sediment dynamics in a hypertidal estuarine channel. Ocean Dyn (in press)
Arakawa A (1972) Design of the UCLA general circulation model. Tech rep. University of California, Los Angeles
Badewien TH, Zimmer E, Bartholomä A, Reuter R (2009) Towards continuous long-term measurements of suspended particulate matter (spm) in turbid coastal waters. Ocean Dyn 59:227–238CrossRef
Bartholomä A, Kubicki A, Badewien TH, Flemming BW (2009) Suspended sediment transport in the German Wadden Sea–seasonal variations and extreme events. Ocean Dyn 59:213–225CrossRef
Baugh JV, Manning AJ (2007) An assessment of a new settling velocity parameterisation for cohesive sediment transport modeling. Cont Shelf Res 27:1835–1855CrossRef
Bolaños R, Souza AJ (2010) Measuring hydrodynamics and sediment transport processes in the Dee estuary. Earth Syst Sci Data 2:157–165CrossRef
Bolaños R, Brown JM, Amoudry LO, Souza AJ (2013) Tidal, riverine and wind influences on the circulation of a macrotidal estuary. J Phys Oceanogr 43:29–50. doi:10.​1175/​JPO-D-11-0156.​1 CrossRef
Brown JM, Bolaños R, Souza AJ (2014) Controls on the mid-term estuarine residuals: circulation and elevation. Ocean Dyn (in press)
Burchard H, Hetland RD (2010) Quantifying the contributions of tidal straining and gravitational circulation to residual circulation in periodically stratified tidal estuaries. J Phys Oceanogr 40:1243–1262. doi:10.​1175/​2010JPO4270.​1 CrossRef
Burchard H, Bolding K, Villareal MR (1999) GOTM, a general ocean turbulence model. Theory, implementation and test cases. Tech rep 18745. European Commission, Brussels, pp 103
Burchard H, Föser G, Staneva JV, Badewien TH, Riethmüller R (2008) Impact of density gradients on net sediment transport into the Wadden Sea. J Phys Oceanogr 38:566–587. doi:10.​1175/​2007JPO3796.​1 CrossRef
Delhez E, Damm P, de Goede E, de Kok J, Dumas F, Gerritsen H, Jones J, Ozer J, Pohlmann T, Rasch P, Skogen M, Proctor R (2004) Variability of shelf-seas hydrodynamic models: lessons from the nomads2 projects. J Mar Syst 45:39–53CrossRef
Dyer KR (1989) Sediment processes in estuaries: future research requirements. J Geophys Res 94(C10):14,327–14,339CrossRef
Fettweis M, Francken F, Pison V, Van den Eynde D (2006) Suspended particulate matter dynamics and aggregate sizes in a high turbidity area. Mar Geol 235:63–74CrossRef
Fugate D, Friedrichs C (2002) Determining concentration and fall velocity of estuarine particle populations using ADV, OBS and LISST. Cont Shelf Res 22:1867–1886CrossRef
Fugate D, Friedrichs C (2003) Controls on suspended aggregate size in partially mixed estuaries. Estuar Coast Shelf Sci 58:389–404CrossRef
Gibbs RJ (1985) Estuarine flocs: their size, settling velocity and density. J Geophys Res 90(C2):3249–3251CrossRef
Giddings SN, Fong DA, Monismith SG (2011) Role of straining and advection in the intratidal evolution of stratification, vertical mixing, and longitudinal dispersion of a shallow, macrotidal, salt wedge estuary. J Geophys Res 116:1–20. doi:10.​1029/​2010JC006482 CrossRef
Holt J, Umlauf L (2008) Modelling the tidal mixing fronts and seasonal stratification of the Northwest European Continental shelf. Cont Shelf Res 28:887–903CrossRef
Holt JT, James ID (2001) An s coordinate density evolving model of the Northwest European Continental shelf 1, model description and density structure. J Geophys Res 160(C7):14,015–14,034CrossRef
Jago CF, Jones SE, Sykes P, Rippeth T (2006) Temporal variation of suspended particulate matter and turbulence in a high energy, tide-stirred, coastal sea: relative contributions of resuspension and disaggregation. Cont Shelf Res 26:2019–2028CrossRef
Jago CF, Kennaway GM, Novarino G, Jones SE (2007) Size and settling velocity of suspended flocs during a Phaeocystis bloom in the tidally stirred Irish Sea, NW European Shelf. Mar Ecol Prog Ser 345:51–62CrossRef
Jay DA, Musiak D (1994) Particle trapping in estuarine tidal flows. J Geophys Res 99:20,445–20,461CrossRef
Jones SE, Jago CF, Bale AJ, Chapman D, Howland RJM, Jackson J (1998) Aggregation and resuspension of suspended particulate matter at a seasonally stratified site in the southern North Sea: physical and biologic controls. Cont Shelf Res 18(11):1283–1309CrossRef
Khelifa A, Hill PS (2006) Models for effective density and settling velocity of flocs. J Hydraul Res 44(3):390–401CrossRef
Kombiadou K, Krestenitis NY (2012) Fine sediment transport model for river influenced microtidal shelf seas with application to the Thermaikos Gulf. Cont Shelf Res 36:41–62CrossRef
Kranenburg C (1994) The fractal structure of cohesive aggregates. Estuar Coast Shelf Sci 39:451–460
Krivtsov V, Howarth MJ, Jones SE, Souza AJ, Jago CF (2008) Monitoring and modelling of the Irish Sea and Liverpool Bay: an overview and SPM case study. Ecol Model 212:37–52CrossRef
van Leussen W (1994) Estuarine macroflocs and their role in fine-grained sediment transport. PhD thesis. University of Utrecht, The Netherlands
Maggi F (2007) Variable fractal dimension: a major control for floc structure and flocculation kinematics of suspended cohesive sediment. J Geophys Res 112(C07012):1–12
Manning A, Baugh J, Spearman J, Whitehouse R (2010) Flocculation settling characteristics of mud: sand mixtures. Ocean Dyn 60:237–253CrossRef
Manning AJ (2008) The development of algorithms to parameterise the mass settling flux of flocculated estuarine sediments. In: Kusuda T, Yamanishi H, Spearman J, Gailani J (eds)Sediments and ecohydraulics: INTERCOH 2005. Elsevier, Amsterdam, pp 193–210CrossRef
McCave IN (1984) Size spectra and aggregation of suspended particles in the deep ocean. Deep Sea Res 31(4):329–352CrossRef
Mikkelsen OA, Pejrup M (2001) The use of a LISST-100 laser particle sizer for in situ estimates of floc size, density and settling velocity. Geo-Mar Lett 20:187–195CrossRef
Polton J, Palmer M, Howarth M (2011) Physical and dynamical oceanography of Liverpool Bay. Ocean Dyn 61:1421–1439CrossRef
van Rijn L (2005) Principles of sedimentation and erosion engineering in rivers, estuaries and coastal seas. Aqua, Amsterdam
van Rijn LC (2007a) Unified view of sediment transport by currents and waves. I: initiation of motion, bed roughness, and bed load transport. J Hydraul Eng 133:649–667CrossRef
van Rijn LC (2007b) Unified view of sediment transport by currents and waves. II: suspended transport. J Hydraul Eng 133:668–689CrossRef
Schuttelaars M, de Jonge V, Chernetsky A (2012) Improving the predictive power when modelling physical effects of human interventions in estuarine systems. Ocean Coast Manag 1–33. doi:10.​1016/​j.​ocecoaman.​2012.​05.​009
Simpson JH, Williams E, Brasseur LH, Brubaker JM (2005) The impact of tidal straining on the cycle of turbulence in a partially stratified estuary. Cont Shelf Res 25:51–64. doi:10.​1016/​j.​csr.​2004.​08.​003 CrossRef
Thurston W (2009) Turbulence as a mediator of processes in a macrotidal estuary. PhD thesis. School of Earth and Environment, University of Leeds, Leeds
Verney R, Lafite R, Brun-Cottan J, Le Hir P (2011) Behaviour of a floc population during a tidal cycle: laboratory experiments and numerical modelling. Cont Shelf Res 31:S64–S83CrossRef
Winterwerp JC (1998) A simple model for turbulence induced flocculation of cohesive sediment. J Hydraul Res 36:309–326CrossRef
Winterwerp JC (2002) On the flocculation and settling velocity of estuarine mud. Cont Shelf Res 22:1339–1360CrossRef
Winterwerp JC, van Kesteren WGM (2004) Introduction to the physics of cohesive sediment in the marine environment. 1st edn. Elsevier, Amsterdam
Winterwerp JC, Manning AJ, Martens C, de Mulder T, Vanlede J (2006) A heuristic formula for turbulence-induced flocculation of cohesive sediment. Estuar Coast Shelf Sci 68:195–207CrossRef
Yuan Y, Wei H, Zhao L, Yiang W (2008) Observations of sediment resuspension and settling off the mouth of Jiaozhou Bay, Yellow Sea. Cont Shelf Res 28:2630–2643CrossRef