Sediment transport in estuaries and the formation of estuarine turbidity maxima (ETM) highly depend on the ability of suspended particulate matter (SPM) to flocculate into larger aggregates. While most literature focuses on the small-scale impact of biological flocculants on the formation of larger aggregates, the influence of the flocculation process on large-scale estuarine SPM profiles is still largely unknown. In this paper, we study the impact of flocculation of SPM on the formation of ETM. For this, a semianalytical width-integrated model called iFlow is utilized and extended by a flocculation model. Starting from a complex one-class flocculation model, we show that flocculation may be described as a linear relation between settling velocity and suspended sediment concentration to capture its leading-order effect on the ETM formation. The model is applied to a winter case in the Scheldt estuary (Belgium, Netherlands) and calibrated to a unique, long-term, two-dimensional set of turbidity (cf. SPM) observations. First, model results with and without the effect of flocculation are compared, showing that the spatial and temporal variations of the settling velocity due to flocculation are essential to reproduce the observed magnitude of the suspended sediment concentrations and its dependence on river discharge. Second, flocculation results in tidally averaged land-inward sediment transport. Third, we conduct a sensitivity analysis of the freshwater discharge and floc breakup parameter, which shows that flocculation can cause additional estuarine turbidity maxima and can prevent flushing of the ETM for high freshwater inflow.
Estuaries often contain regions in which the concentration of suspended particulate matter (SPM) is larger than landward and seaward of these regions. These regions are called estuarine turbidity maxima (ETM) and are caused by fine sediments being trapped as a result of the complex interactions of the water motion and sediment dynamics. Examples of such transport mechanisms are, among others, related to tidal asymmetries in water motion, water density gradients, and transport mechanisms related to temporal variability in settling velocity of SPM [see, e.g., Burchard et al. (2018) for a review]. One mechanism resulting in temporal and spatial variability of settling velocity is flocculation. Flocculation results in the aggregation and breakup of cohesive SPM, thus changing the floc size and hence the settling velocity.
Two main classes of models describing flocculation can be distinguished, the Lagrangian flocculation (LF) models and the population balance equation (PBE) models (Lai et al. 2018). PBE flocculation models typically consist of multiple discrete size classes and compute the evolution of the number of flocs within each class over time (Smoluchowski 1918; Verney et al. 2011). Although PBE models can calculate the evolution of particle size distribution over time, they come with high computational costs. In contrast, the LF models are dynamic models that resolve the dynamics of particle size and settling velocity by computing a balance between floc aggregation and floc breakup (e.g., Winterwerp 2002). Floc aggregation and breakup are implemented using empirical formulations, which depend on both abiotic and biotic factors (Dyer 1989; van Leussen 1994; Lai et al. 2018).
A first important abiotic factor that impacts flocculation is the suspended sediment concentration. Both in situ measurements (Pejrup and Mikkelsen 2010) and laboratory experiments (Tran et al. 2018) show that floc size typically increases with increased suspended sediment concentration. Another important abiotic driver is turbulence, which both promotes aggregation through enhanced mixing and breakup by increasing the shear stresses on the flocs. Turbulence promoting breakup is observed in various estuaries, such as the Yangtze estuary (Guo et al. 2017) and the Scheldt estuary (Manning et al. 2007; Schwarz et al. 2017). Also, laboratory experiments show that turbulence can decrease the averaged floc size (Mietta et al. 2009). A third abiotic environmental condition that impacts flocculation is salinity. Edzwald et al. (1974) conducted laboratory experiments in which they showed that salinity promotes aggregation of clay particles and very little salt (~5 ppt) is already sufficient to reach a maximum impact of salinity on flocculation. In contrast, Eisma et al. (1980) and Verney et al. (2009) found that salt flocculation is not a crucial factor in the Rhine estuary and Seine estuary, respectively. Biotic characteristics that impact flocculation are the organic content that directly alters the differential density and structure of flocs (Kranenburg 1994; van Leussen 1994). Furthermore, organic content impacts, for example, the floc strength and collision efficiency (Winterwerp and van Kesteren 2004), averaged floc size (Mietta et al. 2009), and floc breakup (Alldredge et al. 1990). Finally, in situ observations show a correlation between chlorophyll a and flocculation efficiency (Verney et al. 2009) and sticky biotic substances (i.e., transparent exopolymer particles) and floc strength (Fettweis et al. 2014).
Although temporal and spatial variations in settling velocity, which is related to floc size, are recognized as potentially important mechanisms for sediment trapping, most models for estuarine sediment transport use either a constant settling velocity or empirical flocculation relationships. Only a few studies included full floc dynamics in a sediment transport model (e.g., Ditschke and Markofsky 2008; Xu et al. 2010; Sherwood et al. 2018; Shen et al. 2018), but the importance of flocculation on sediment transport mechanisms cannot be clearly identified from these model results.
This study aims to model and gain an understanding of the impact of varying settling velocity on large-scale sediment transport and the corresponding development of ETM. As a first step, we consider estuaries in which the suspended sediment concentration is lower than 1 g L−1, allowing us to focus on flocculation process and neglect hindered settling effects, which also result in a temporally and spatially varying settling velocity.
To be able to gain an understanding of the underlying mechanisms that result in sediment trapping and quantify the relative importance of flocculation on ETM formation, we extend the width-averaged, hydrodynamics and sediment transport model known as “iFlow” (Dijkstra et al. 2017) to include flocculation processes. The iFlow model is specifically geared toward gaining an understanding of the water motion, sediment transport, and trapping in tidally dominated estuaries and allows for extensive sensitivity analysis. The model is extended with the one-class LF model of Winterwerp (2002), in which the impact on flocculation of suspended sediment concentration, turbulent shear, and floc size and structure are directly parameterized in the floc aggregation and breakup terms. Other effects, such as salinity and biotic factors, are indirectly included by fitting to observations.
We apply the extended iFlow model to the Scheldt estuary. The Scheldt estuary is located in Belgium and the Netherlands and debouches into the North Sea. Being a well-documented estuary (Meire et al. 2005), there have been intensive monitoring campaigns (e.g., Maris and Meire 2016) and modeling experiments. Moreover, the iFlow model without flocculation has already successfully been applied to the Scheldt estuary (Brouwer et al. 2018; Dijkstra et al. 2019a). The flocculation process is expected to be important in the Scheldt estuary (Peters 1972; Gourgue et al. 2013).
In this paper, we extend the iFlow model to include the effect of flocculation in section 2. In section 3, we calibrate the model to a unique, long-term turbidity (cf. SPM) dataset of the Scheldt estuary. In section 4, we apply the coupled flocculation–sediment transport model to the Scheldt estuary and study the effect of flocculation on the large-scale sediment distribution and the underlying sediment transport mechanisms.
In this section, we first present the sediment transport model (i.e., iFlow) and flocculation model of Winterwerp (2002). Next, we use scaling and perturbation methods to simplify the coupling of the two models. Then, we describe the theoretical impact of flocculation on sediment transport. Last, we introduce the (numerical) implementation of our model.
The iFlow model (Dijkstra et al. 2017) was developed to obtain the width-averaged water motion, sediment transport, and trapping in a tidally dominated estuarine system by solving the width-averaged shallow water equations, suspended sediment concentration equation, and a dynamic equation for the erodibility of the bed (Brouwer et al. 2018). The water motion is assumed to be forced by an M2 and M4 tidal signal at the mouth, while a river discharge is prescribed at the upstream tributaries. It is assumed that the M2 tidal signal dominates over the M4 and riverine signal. We decompose the water motion into a longitudinal and vertical velocity component, denoted by u and w. The suspended sediment concentration follows from an advection–diffusion equation and is forced by a constant suspended sediment concentration at the mouth and prescribed inflow of sediment at the upstream tributaries. The bathymetry and width are approximated by smooth polynomials, thereby neglecting the impact of small-scale gradients, as illustrated in Fig. 1. A diagnostic longitudinal salinity gradient is prescribed, consisting of a depth- and time-independent sigmoid function.
By assuming that the dynamics consist of strong tidal signals beside a subtidal part, the iFlow model solves the problem in frequency space for subtidal, M2, and M4 components. The water surface elevation is assumed small compared to the subtidal water depth, resulting in a small dimensionless parameter, given by the ratio of water surface elevation and water depth. Using a scaling procedure and relating the various dimensionless numbers to this small parameter allows for identifying the relative importance of individual processes (e.g., river, tide, density forcing, advection). Next, a perturbation approach is employed to solve the equations at different order. For further details on the iFlow model, we refer the reader to Dijkstra et al. (2017).
Using this information and truncating frequencies larger than M4 (i.e., M6, M8, etc. are neglected and, hence, not discussed), the leading-order water motion u0 and suspended sediment concentration c0 consist of an M2 tidal signal and subtidal and M4 tidal components, respectively (Chernetsky et al. 2010; Dijkstra et al. 2017), denoted by u02, c00, and c04 (in which the first superscript denotes the order and the second one denotes the frequency) and defined as
with x and z being the longitudinal and vertical coordinates and t being the time. At first order, the water motion u1 and suspended sediment concentration c1 consist of a subtidal and M4 tidal signal and an M2 tidal signal, respectively (Chernetsky et al. 2010; Dijkstra et al. 2017), denoted by u10, u14, and c12 and defined as
The net sediment transport , i.e., the sediment transport averaged over a tidal period (denoted by angle brackets), is given by the sum of the advective and diffusive sediment transport integrated over the cross section:
where B is the local channel width, R is the reference level of the water surface elevation, H(x) is the local water depth, and Kh is the horizontal eddy diffusivity coefficient. Moreover, we used that only subtidal–subtidal, M2–M2, and M4–M4 tidal interaction terms result in net sediment transport, which is a direct result of the orthogonality of the complex exponential functions.
In iFlow, parameters, such as the total freshwater discharge Q, horizontal eddy diffusivity coefficient Kh, and settling velocity , are usually prescribed as constant parameters [although this is not needed; see Dijkstra et al. (2017)], allowing for a quick solution procedure. To include flocculation in the model, the settling velocity ws has to be coupled to the suspended sediment concentration c, which makes the model nonlinear, requiring an iterative solution procedure (see section 2e).
b. Flocculation model
The Winterwerp flocculation model employed in this paper reads (Winterwerp 2002) as
where N is the number of flocs per unit volume, ws is the settling velocity, Kυ is the vertical eddy diffusivity coefficient, G is the shear rate, Df and Dp are the floc and primary particle size, and and kB are the empirical aggregation and breakup parameters. To acquire Eq. (6), we assumed a fractal dimension nf = 2 and we set model calibration parameters q = 0.5 and p = 1, following Winterwerp (2002).
The aggregation parameter and breakup parameters kB define the efficiency of the flocculation process and represent the intrinsic flocculation kinetics of the system of interest. Laboratory flocculation experiments show that floc kinetics depend, among others, on salinity (Peters 1972; Edzwald et al. 1974). To include a salinity dependence in our model, we allow the flocculation kinetics to depend on the depth- and tidally averaged salinity profile, resulting in and kB being functions of the longitudinal coordinate x. In section 2c, we show that net sediment transport only depends on the ratio of and kB and not on the individual parameters. From the experiments carried out by Edzwald et al. (1974), we therefore postulate the following dependence of the ratio of and kB on salinity S:
with S1 and S2 being empirical parameters that are calibrated to a salinity S(x) dataset. For more details, the reader is referred to appendix A.
Following Pejrup and Mikkelsen (2010), the shear rate G reads as
where κ′ is the von Kármán constant, ν is the kinematic viscosity, is the relative water depth, and ⟨u*⟩ is the subtidal friction velocity:
where z1 is the distance in the logarithmic layer above the bed, is the bed roughness coefficient, u0 is the leading-order, longitudinal flow velocity [see Eq. (1)], and is the subtidal, first-order, longitudinal, river-induced flow velocity [see Eq. (3)]. For simplicity, we only consider the impact of subtidal shear velocity u* at leading order. The first-order, subtidal, riverine contribution is also included (i.e., ) because this constituent is dominant at the upstream border. Apart from neglecting tidal variations of G, we also approximate G by its value at z1 = H/2 and use this as a proxy for the whole water column. Working with a depth-averaged value of G instead of a midpoint value has no major impact on our results and conclusion because vertical variations in G are relatively small (typically <10% in our case study) in the logarithmic layer. Although turbulence is an important driver of flocculation, we do not focus on this variable in our case study because we assume a tide-dominated system with low stratification. This results in a shear rate which has a similar order of magnitude within the ETM.
The Winterwerp model depends on the number of flocs per unit volume N, floc size Df, and settling velocity ws. These three variables are not independent. We rewrite Eq. (6) such that it solely depends on the settling velocity ws. First, we express the number of flocs per unit volume N in terms of the floc size Df (Kranenburg 1994; Winterwerp and van Kesteren 2004):
with fs being the floc shape factor and ρs being the floc density.
Next, to obtain an equation for the settling velocity, we use the generalized Stokes formulation to translate floc size Df to settling velocity ws (e.g., Winterwerp and van Kesteren 2004):
where ρw is the water density, g is the gravitational acceleration, and μ is the viscosity of water. We assumed spherical flocs (fs = π/6), a fractal dimension nf = 2, and that the floc Reynolds number Ref= wsDf/ν ≪ 1.
with γ defined in Table 1. In the following section, we approximate this equation and write the settling velocity ws as an explicit function of the suspended sediment concentration c.
c. Ordering of the flocculation model
1) Leading order
At leading order, the flocculation model reduces to the balance between floc aggregation and floc breakup:
This can be interpreted as if the flocculation process is instantaneous and local: there is no inertia and transport. As a result, the settling velocity scales linearly with the suspended sediment concentration:
where β and κ are defined in Table 1 and depend on several parameters including the shear rate and ratio between and kB. Using the definitions of β and κ (Table 1), we see that the βκ term is equal to the settling velocity ws,min corresponding to the settling velocity of primary particles. Indeed, it is equal to the Stokes formulation for massive (nf = 3) primary particles with primary particle size Dp. This ensures that for c → 0, . Moreover, our leading-order result in Eq. (14) is equivalent to the equilibrium floc size De presented in Winterwerp (2002):
The positive correlation of ws with suspended sediment concentration in the leading-order result complies with observations. Indeed, Pejrup and Mikkelsen (2010), assuming a power relationship between settling velocity and suspended sediment concentration, found exponents ranging between 0.47 and 2.9 with an average of 1.29, using 18 different tidal systems based on real measurements. Recently, Tran et al. (2018) studied the influence of SPM concentration on floc size in laboratory mixing tanks. They found a linear relationship between floc size and SPM concentration, which again results in the linear relationship between settling velocity and SPM concentration as given in Eq. (14).
2) First order
At first order, our scaling procedure shows that local inertia (cf. ∂tN term), settling, and vertical diffusion become significant [see Eq. (C5)]. Advection and horizontal diffusion are still negligible at this order. The solution then reads as
where τ is defined in Table 1. Apart from a contribution to that again scales linearly in c (βc1 term), we obtain additional terms when compared to the leading-order result. Because of inertial effects, settling, and vertical diffusion, we obtain an additional contribution that is constant in time and one that varies on the M4 time scale, no additional M2 contribution is generated [see Eqs. (D30) and (D31)]. As we show below, we are only interested in the M2 contribution because this is the only first-order contribution that is important for net sediment transport. The only M2 signal in the first-order settling velocity is related to the first-order balance between aggregation and breakup and reads as
d. Impact of flocculation on net sediment transport
In the previous section, we discussed the impact of suspended sediment concentration c on settling velocity ws, which is a direct consequence of flocculation. In the following, we present the influence of settling velocity ws on suspended sediment concentration c. As flocculation alters the leading and first-order suspended sediment concentration c0 and c1, the corresponding net sediment transport from Eq. (5) is also impacted. We rewrite Eq. (5) as
and examine the relation of , , and with the settling velocity ws.
A careful analysis shows that at leading order, both the settling velocity and suspended sediment concentration consist of a subtidal and M4 tidal signal. This can be inferred from the following argument: if at leading order, the (truncated) suspended sediment concentration c0 consists of only a subtidal and M4 tidal signal as found for a time-independent ws [see Eq. (2)], the settling velocity also gets an M4 tidal component [see Eq. (14)]. This M4 signal impacts the leading-order suspended sediment concentration c0 through the nonlinear interaction term in the leading-order differential equation for c0 [see Eq. (B13)]. Interestingly, the nonlinear interaction of c0 and results in an adjustment of the subtidal and M4 tidal signal in c0, that is, c00 and c04, and generates an infinite number of additional tidal frequencies in c0 because of its nonlinearity:
As stated before, frequencies larger than M4 are truncated.
Additionally, the subtidal and M4 tidal signal in impacts the M2 signals in the first-order differential equation for c1 through the forcing term [see Eq. (B16)], which again results in an adjustment of net sediment transport through :
Consequently, at leading order both the subtidal and M4 tidal signal in have an impact on net sediment transport through its impact on and and .
At first order, only the M2 tidal signal in , i.e., , results in net sediment transport through the forcing term in the first-order differential equation for c1 [see Eq. (B16)]:
e. Solution method
iFlow contains a numerical sediment module that solves for suspended sediment concentration c numerically in x and z using a grid of 200 cells in the x direction and 50 cells in the z direction. To do so, the module requires a leading-order settling velocity and first-order settling velocity as input. A new module is added to iFlow containing the explicit analytical expression of and [Eqs. (14) and (17)]. The coupling of ws to c results in a ws which is a function of time and space. Because the coupling is nonlinear, we solve the coupling iteratively using the Picard method. We start with an initial condition for and as input to the sediment module to acquire a suspended sediment concentration at leading order c0 and first order c1. Next we use c0 and c1 to recalculate and using Eqs. (14) and (17). We use the recalculated and as new input to the sediment module. We repeat this procedures until the solution converges using the following stop criterion:
in which is the leading-order settling velocity of the previous iteration. The iterative procedure is schematized in Fig. 2.
3. The Scheldt: Observations and calibration
We apply our coupled flocculation–sediment transport model to a winter case (i.e., January–March) in the Scheldt estuary. To do so, we calibrate the model to a unique, long-term turbidity (cf. SPM) dataset for 2015 until 2017 measured within the environmental monitoring program called Onderzoek Milieu Effecten Sigmaplan (OMES). We divide the turbidity dataset in winter and summer because the turbidity data shows a strong seasonality in the Scheldt estuary (Maris and Meire 2016). In this section, we present the Scheldt estuary, the long-term turbidity dataset, and the calibration method.
a. The Scheldt
The Scheldt estuary (Fig. 3) is a funnel-shaped estuary of approximately 160 km long, which flows through Belgium and debouches into the North Sea near Vlissingen (Netherlands). The Scheldt estuary has a relatively low averaged total freshwater discharge (172 m3 s−1 in our winter case) and can be considered a tide-dominated estuary (Meire et al. 2005; Waterinfo.be 2019). The main tributaries are the upper Scheldt (upstream boundary; ~34% of total river discharge), the Rupel (at 95 km from the mouth; ~54% of total river discharge), and the Dender (at 123 km from the mouth; ~12% of total river discharge) (Waterinfo.be 2019). Parameter values used in our winter case study in the Scheldt estuary are summarized in Table 2. These parameter values follow from Brouwer et al. (2018) if not mentioned explicitly in the text. The total river discharge Q, erosion parameter M, and floc breakup coefficient kB follow from calibration to the OMES SPM dataset presented in the following section.
b. OMES data
We calibrate our coupled flocculation–sediment transport model to a unique dataset representing a typical winter situation of SPM distribution. In the frame of the ongoing OMES monitoring, vertical turbidity profiles were obtained using conductivity–temperature–depth (CTD) and turbidity casts from aboard a ship at 16 fixed locations spread over the Belgian part of the Scheldt estuary during monthly or biweekly campaigns (Fig. 3). The campaigns were conducted independently of the tidal phase and spring neap tide. To obtain a typical winter distribution, we temporally averaged the January–March data from 2015 to 2017.
Turbidity was measured at various depths using an optical backscatter point sensor (OBS) of RBR type XR420 CTD+. During each campaign, the sensor was calibrated using a Formazine solution standard (Maris and Meire 2016). The water body was vertically profiled using a minimal sampling frequency of 10−1 s−1.
Simultaneously, two water samples were collected at each location using a Niskin bottle at approximately half the water depth and at the water surface. SPM concentrations were gravimetrically determined after filtration in the laboratory. During each campaign, 16 × 2 SPM water samples were collected, resulting in 32 SPM estimates. To translate turbidity to SPM, we used these 32 samples to apply a linear data fit. We assumed that the relation between turbidity and SPM is location (and time) independent, that is, equal for every location within one campaign. Calibration of turbidity to SPM resulted in a depth profile of SPM at 16 fixed locations in the Flemish part of the Scheldt estuary.
The winter case covers 11 measuring campaigns. We excluded data at the water surface and river bed because here we typically have distortions due to, for example, air bubbles and high turbidity, respectively. To do so, we manually excluded depths at which we did not measure 8 of 11 measuring campaigns. We averaged the SPM concentration observations of the 11 measuring campaigns at each depth. We assumed this averaged value approximates the residual SPM concentration following, for example, Dijkstra et al. (2019a) and Cox et al. (2019). This assumption is valid when the number of estimates is sufficiently high, so averaging of the periodic temporal variability of the SPM concentration is negligible compared to the magnitude of the residual SPM concentration. Here, we assumed that the campaigns are randomly distributed within the tidal phase and spring neap tide, which was shown to be valid for OMES SPM sampling between 1995 and 2015 (Vandenbruwaene et al. 2016). Furthermore, we assumed that the SPM distribution within the time frame of our winter case is fixed on an estuarine scale (i.e., the scale of the iFlow model in the longitudinal direction). So local effects due to, for example, temporal variations in river discharge can be neglected. Given a system-averaged standard deviation relative to the time-averaged SPM concentration of 0.43, we conclude that these assumptions are acceptable and that the number of estimates is sufficiently high.
On average, each turbidity profile consists of 20 individual observations, with an average distance between two consecutive observations of 0.36 m. We vertically interpolated the data to enable a comparison of the subtidal SPM model output at the vertical model grid cells in our calibration. This results in a total number of 606 (x, z) locations at which we compared data and model output. This number results from the fact that we measured at 16 stations and have 50 vertical model grid cells (total of 800 grid cells) but excluded data near the water surface and river bed. The average number of data points n for each location (x, z) is 9.4.
We calibrated the linearized bed roughness coefficient such that the modeled M2 water levels match observations following Dijkstra et al. (2017). In our coupled flocculation–sediment transport model, we have three additional calibration parameters:
floc breakup parameter kB,
erosion parameter M, and
river discharge Q.
The iFlow model approximates the freshwater inflow by a subtidal discharge. In reality, discharge shows a significant temporal variation and observed average discharge is not always representative for average sediment transport. For example, the standard deviation of the averaged river discharge corresponding to our winter case is equal to 92 m3 s−1 (Waterinfo.be 2019). To improve the correspondence between model output and observed sediment distribution, we included Q as a calibration parameter.
We run our model for various values of the calibration parameters kB ∈ [2400, …, 8000] s1/2 m−2, M ∈ [0.1, …, 4] × 10−3 s m−1, and Q ∈ [172, …, 340] m3 s−1 and compare the model result with the measurements for each setting to find the best parameter values. The choice of the range in kB and M is based on scaling (see Table C1 in appendix C) and observations (Zhu et al. 2017), respectively. The range of Q represents from the 60th until the 90th percentile of the observations corresponding to our winter case (Waterinfo.be 2019).
To objectively compare the model output to the measurements, we construct a cost function, which is based on a statistical two-tailed t test. For each SPM data point location (x, z), we compute whether the following test statistic:
has a t distribution with n − 1 degrees of freedom. If found to be true, we accept that the model output and data are statistically equal at location (x, z). Here, ⟨cmodel(x, z)⟩ and ⟨cmeas(x, z)⟩ are the model subtidal suspended sediment concentration output and time-averaged measured SPM concentration at location (x, z), respectively, std is the standard deviation of the measured SPM concentration at location (x, z), and n is the total number of data points at location (x, z). We choose a significance level α = 0.05. Last, we define a cost function value as
where Nequal are the number of locations at which the suspended sediment concentration model output and data are statistically equal and Nloc is the total number of locations for which we have data. Consequently, the cost function output is a value between 1 (mismatch of model output and data) and 0 (perfect match of model output and data).
Evaluated by our cost function, we obtain an optimal parameter set kB = 5600 s1/2 m−2, M = 0.003 36 s m−1, and Q = 233 m3 s−1. The cost function output for various kB, M, and Q is presented in Fig. 4. We acquire a higher discharge of 233 m3 s−1 in winter when compared with the time-averaged value of 172 m3 s−1. This larger discharge is to be expected because of the asymmetry in flushing and buildup of an ETM: the impact of high river discharges on suspended sediment concentration is larger than the impact of low river discharges. Indeed, Brouwer et al. (2018) showed that the time scale of building up an ETM is two orders of magnitude larger (~102 days) than the flushing of an ETM (~100 days).
To analyze the relative importance of flocculation on large-scale sediment transport in the Scheldt estuary, we compare the modeled suspended sediment distribution and the different terms contributing to net sediment transport with and without flocculation. Next, we apply a sensitivity analysis of the calibration parameters Q and kB when including flocculation and compare the results to those without flocculation.
a. Impact of flocculation on large-scale suspended sediment distribution
Application of the model to the winter case in the Scheldt estuary and model calibration results in the model output presented in Fig. 5.
Figure 5a shows the result of the subtidal suspended sediment concentration with flocculation. Figure 5b shows the long-term, time-averaged suspended sediment concentration dataset for winter in 2015–17. These observations agree with SPM samples simultaneously taken at the water surface, which are presented in Cox et al. (2019). When we include flocculation (Fig. 5a), the model output complies both quantitatively and qualitatively with observations (Fig. 5b). The location and intrusion of the ETM are well captured. The model shows typical suspended sediment concentrations of 100–300 mg L−1, which agrees with observations.
Figure 5c shows the model output of the leading-order settling velocity corresponding to the case with flocculation. Qualitatively this figure corresponds to the suspended sediment concentration in Fig. 5a, which is to be expected because we showed that the leading-order settling velocity scales linearly to the suspended sediment concentration [see Eq. (14)]. The distribution of is only (slightly) altered by the imposed longitudinal variation in shear rate G and aggregation parameter [see Eqs. (7) and (8)]. Hence, the dependence of on c0 is the most important factor affecting .
Figure 5d shows the subtidal suspended sediment concentration without including flocculation, using a constant ws and the same river discharge as used in Fig. 5a. In the simplest model case in which we excluded flocculation, a constant ws of 3 mm s−1 yields the best comparison with observations (Fig. 5b). This value corresponds to the value of the settling velocity found near the ETM (Fig. 5c). When we exclude flocculation (Fig. 5d), the model output qualitatively agrees with observations; we acquire an ETM at 60–80 km from the mouth and obtain more sediment near the bed than at the water surface. However, quantitatively the suspended sediment concentrations are up to a factor three too low. Furthermore, the intrusion of the ETM is too small. Note that our results are different from Dijkstra et al. (2019a) and Brouwer et al. (2018), since they only considered much lower discharges than what is considered here. Because we are interested in the individual impact of flocculation on the suspended sediment distribution, we did not recalibrate parameters M and Q for the case excluding flocculation. When these parameters are recalibrated, M = 0.004 s m−1 and Q = 172 m3 s−1, resulting in a cost function value of 0.87. The modeled suspended sediment concentrations are still too low, although they are slightly higher than the previous case due to a decrease of Q. A further increase of M does not result in an increase of the suspended sediment concentrations because we are in a suspended sediment supply-limited condition (see Brouwer et al. 2018); that is, all sediment is eroded from the bed during a tidal cycle. This implies that the amount of sediment trapped, and not the erosion from the bed, is the limiting factor for the observed suspended sediment concentrations. Sediment trapping is related to net sediment transport processes, which are clearly less efficient in importing sediment when flocculation is not taken into account (see next section).
b. Impact of flocculation on net sediment transport
In this section, we focus on the underlying sediment-transport mechanisms that result in sediment transport and the buildup of an ETM presented in Fig. 5. For this, we introduce the advective net transport capacity T, which equals the net redistribution of a (longitudinally) uniform layer of sediment on the bed (Dijkstra et al. 2019b). The net transport capacity only depends on hydrodynamic and sediment features and not on the location of the ETM. Therefore, it is a suitable measure to analyze the impact of various mechanisms on the total net sediment transport (Dijkstra et al. 2019a,b).
Figure 6 shows the six most important transport mechanisms resulting in net transport capacity with (Fig. 6a) and without (Fig. 6b) flocculation together with the total net transport capacity Ttotal. This total net transport capacity Ttotal is the sum of all separate net transport capacities. When Ttotal = 0 and Ttotal changes sign from positive to negative (in land-inward direction), we have a convergence point. Here, we expect sediment to accumulate and thus the formation of an ETM. When comparing Figs. 6a and 6b, we observe that the convergence point at approximately 80 km in the case with flocculation (Fig. 6a) shifted seaward to approximately 70 km when flocculation is not considered (Fig. 6b). In both cases, a significant decrease in total net transport capacity is present at 95 km, which is due to the Rupel tributary, resulting in downstream net sediment transport. Figure 6a also shows that the additional mechanism related to flocculation, in Eq. (21), is always positive and consequently results in land-inward net sediment transport. Furthermore, it is the dominant mechanism at the mouth.
By comparing Figs. 6a and 6b, it is clear that the contributions of the other transport mechanisms changed when flocculation is included. This is due to the spatial and temporal (M4) variations in the leading-order settling velocity , which were absent in the case without flocculation. As mentioned earlier, this alteration of the subtidal and M4 tidal signal in impacts all contributions to the net sediment transport [see Eqs. (19) and (20)].
c. Sensitivity analysis of floc breakup parameter kB and total freshwater discharge Q
A sensitivity analysis of the floc breakup parameter kB and freshwater discharge Q allows us to compute the impact of the magnitude of kB and Q on the formation of ETM. To stress the importance of spatially varying ws, we show the sensitivity of ETM formation to kB with flocculation only varying in space, not in time. We also compare the case with and without flocculation to calculate the impact of space and time dependence in ws on the emergence of ETM. Finally, we vary the constant settling velocity in the case without flocculation to assess the sensitivity of the ETM characteristics for this parameter.
Figure 7 shows the sensitivity of varying kB (Q = 233 m3 s−1 and M = 0.003 36 s m−1) to depth-averaged subtidal suspended sediment concentration . In Fig. 7a the time dependence in ws is included in the model, whereas in Fig. 7b we excluded the temporal dependence in ws (i.e., and ) and only considered the spatial dependence in ws (i.e., ). By excluding the temporal dependence in ws, we remove the transport mechanism depicted in Fig. 6a by the solid black line and the impact of the M4 tidal signal in on net sediment transport [see Eqs. (19)–(21)].
Figure 7 shows that the magnitude of kB not only determines the location and intensity of ETM, but also the number of ETM. In Fig. 7a, for kB values larger than a critical value of 8000 s1/2 m−2, no ETM is found. For values between 4000 and 8000 s1/2 m−2, one ETM is found at a fixed location. For values between 2400 and 4000 s1/2 m−2, two ETM are found, with the most downstream ETM moving downstream with decreasing kB. This behavior is strongly linked to the temporal variability of settling velocities introduced by the flocculation process. Indeed, Fig. 7b shows the impact on ETM formation when only spatial variability in ws is incorporated. In that case, ETM formation only occurs at significantly smaller kB values (≈3000 s1/2 m−2), while multiple ETM are not observed in the domain. Moreover, modeled suspended sediment concentrations are much lower when ignoring temporal variations in ws.
Figure 8 shows the impact of varying river discharge Q (M = 0.003 36 s m−1) on the depth-averaged subtidal suspended sediment concentration with flocculation (Fig. 8a, kB = 5600 s1/2 m−2) and without flocculation (Fig. 8b; = 3 mm s−1). We keep Q at typical winter values. Similar to the sensitivity in kB, the magnitude of Q not only determines the location and intensity of ETM, but also the number of ETM. When flocculation is incorporated (Fig. 8a), we have one ETM for values between approximately 200 and 340 m3 s−1, with suspended sediment concentrations that compare well with observations (Fig. 5b). For values between approximately 170 and 200 m3 s−1, we find an additional ETM at approximately 100–120 km. This second ETM has been observed in the Scheldt estuary when discharges are low (Cox et al. 2019). An ETM appearing in the tidal freshwater region during lower discharges was observed in various estuaries (e.g., Uncles et al. 2006; Allen et al. 1980). Both ETM move downstream with increasing Q. In the model with fixed settling velocity, only a single ETM is found, which shifts again downstream with increasing Q (Fig. 8b). Moreover, the suspended sediment concentrations obtained are low. Consequently, flocculation is important to obtain quantitatively realistic results for relatively large and observed freshwater discharges Q.
Figure 9 shows the impact of varying constant settling velocity ws on in case effects of flocculation are not considered (Q = 233 m3 s−1, M = 0.003 36 s m−1). The magnitude of the settling velocity determines the existence of an ETM. We need a minimal constant settling velocity of approximately 1.8 mm s−1 to obtain one ETM. The optimal choice of the constant settling velocity is approximately 3 mm s−1, resulting in the best comparison with observations. A settling velocity = 3 mm s−1 corresponds to a typical value near the ETM (Fig. 5c). Moreover, the depth-averaged subtidal suspended sediment concentrations are relatively low compared to the data, for all . At = 3 mm s−1, although the locations of the ETM are equal, the suspended sediment concentrations are still approximately one-third of those found with time-independent flocculation (Fig. 7b), stressing the importance of spatially varying ws. We cannot resolve the modeled low concentrations (cf. Fig. 8b) by altering . Consequently, flocculation is required to acquire quantitatively realistic concentrations for large freshwater discharges in winter.
Our results show that flocculation promotes land-inward net sediment transport in the Scheldt. Relative to the iFlow model without flocculation, the model with flocculation is less sensitive to freshwater discharge. As a result, we still obtain an ETM for high freshwater discharges, which complies with observations (Cox et al. 2019). In addition, the iFlow model without flocculation results in depth-averaged suspended sediment concentrations that are typically too low for large freshwater discharge, stressing the importance of a temporal and spatial dependency of the settling velocity ws caused by flocculation.
a. Indications of spatial and temporal variations of the settling velocity in the Scheldt estuary
Observations in the Scheldt estuary confirm both this temporal dependence in ws, which generates net sediment transport and the modeled magnitude of ws. Although direct settling velocity measurements are scarce, Manning et al. (2007) measured an averaged settling velocity of macroflocs of 3.9 mm s−1 in the ETM of the Scheldt at the entrance of Deurganckdok, approximately 62 km from Vlissingen. When we compare these measured settling velocities to the model output, we find that they more or less comply (Fig. 5c). Besides direct settling velocity measurements in the Scheldt estuary, Fettweis and Baeye (2015) observed a strong M2 and M4 tidal signal in floc size and settling velocity in the Southern North Sea, which complies with our model results [Eqs. (14) and (17)]. Furthermore, Schwarz et al. (2017) measured in situ settling velocities over one tidal cycle in the main channel of the Sieperda March (near the Dutch–Belgian border in Fig. 3) using both the Stokes formulation and the Reynolds-flux method. They also found a distinct M4 tidal signal in settling velocity and floc size, which again agrees with our model results.
b. Flocculation resulting in land-inward net sediment transport
We showed that flocculation alters the sediment transport mechanisms and generates an additional and important sediment transport mechanism (Figs. 6a,b). Our method enabled a systematic analysis showing that flocculation, and more specifically, the M2 tidal signal in the settling velocity, results in land-inward net sediment transport (solid black line in Fig. 6a). Our results comply with Winterwerp (2011) who showed that including flocculation results in land-inward net sediment transport in the Ems River using a 1D vertical model. Xu et al. (2010) also concluded that flocculation results in land-inward sediment transport and promotes accumulation of suspended sediment by performing an idealized 2D model study in the Upper Chesapeake Bay. In general, the direction of net sediment transport depends on the tidal phase difference between flow velocity, suspended sediment concentration, and settling velocity, which might differ for various estuaries and is affected by both the spatial and temporal variations of the settling velocity caused by flocculation.
To estimate a condition for which flocculation results in land-inward net sediment transport, we only consider the M2 tidal signal in the settling velocity , which we showed can have a significant impact on net sediment transport. When the longitudinal water velocity u (in land-inward direction) is in phase with the suspended sediment concentration generated by , we have the condition for maximum land-inward sediment transport related to the flocculation contribution (see Fig. 10).
To meet this condition, it is found that has to be shifted by a −π/2 phase relative to u, meaning that the M2 tidal signal in the settling velocity peaks at slack tide from flood to ebb (see Fig. 10). This condition holds assuming that the water column is fairly well mixed and and c04 are relatively small in comparison with c12 and c00. This condition is modified in other cases. Our reasoning is generalized in appendix E for c04 not relatively small when compared with c00 and general (i.e., not restricted to maximal) land-inward sediment transport.
c. The potential impact of flocculation on sediment transport in the Scheldt
Our model results show that the floc breakup parameter kB determines the position, intensity, and existence of ETM (Fig. 7) which complies with Xu et al. (2010) who found that particle stickiness (cf. kB) can have a significant influence on sediment trapping, especially when the stickiness is small (cf. kB large). As mentioned earlier, various factors such as salinity (Edzwald et al. 1974) and biotic sticky substances (Fettweis et al. 2014) alter kB. In our case study, the direct impact of salinity on flocculation and resulting sediment transport is relatively small. This complies with Einstein and Krone (1962), who concluded that variations of salinity in most of an estuary have only a small effect on the bond strength. Moreover, the experiments of Edzwald et al. (1974) show a relatively low impact in comparison with the minimal aggregation coefficient (cf. kB) and saturation at relatively low salinity. Salinity might also have an indirect impact on flocculation through, for example, the production of biotic sticky substances (e.g., Alldredge et al. 1993; Bar-Zeev et al. 2015), which is only implicitly included through calibration of kB.
Over the last two decades, the water quality has significantly improved in the Scheldt estuary because of the implementation of wastewater treatment in Brussels in 2006 (Brion et al. 2015; Cox et al. 2019). Cox et al. (2019) suggested that changes in water quality might result in a large-scale impact on suspended sediment distribution through its influence on flocculation. Our flocculation model can be used to quantify the impact of water quality improvement and verify the hypothesis of Cox et al. (2019), and supplies a tool to systematically assess the mechanisms that explain the changes to the sediment transport. The main unknown for this is the relation between water quality and kB, and more research is needed to establish practical relations for this.
d. Parameter variations in the Lagrangian flocculation model
To determine the individual impact of flocculation on the large-scale suspended sediment distribution, we compared a sediment transport model using a constant settling velocity and using the complex Lagrangian flocculation model of Winterwerp (2002). As a first step, we kept most parameters fixed in the Lagrangian flocculation model following Winterwerp (2002).
Over the last decade, various authors extended this flocculation model. For example, Maggi (2009) extended the model by separating the floc volume in a mineral and biomass fraction. More recently, Xu and Dong (2017) were able to model floc size dynamics more accurately when assuming that the fractal dimension nf follows a normal distribution. As a final example, Kuprenas et al. (2018) implemented a dependence of floc size Df and shear rate G in the calibration parameter q to correct for an overestimation of floc size at large shear rates (order of 50 s−1).
In further research, the same scaling and perturbation procedure presented in this paper can be applied to extended versions of the Lagrangian model of Winterwerp. For example, by assuming constant parameter values (e.g., nf = 2 and q = 0.5), we found a linear relationship between the settling velocity and suspended sediment concentration. Different parameter values might result in a nonlinear relationship. An example of such a relation is given by van Leussen (1994), who proposed an at equilibrium relation between settling velocity, suspended sediment concentration, and shear rate: , in which a, b, K, and m are empirically determined parameters. Such a formulation is similar to ours in that it depends on c and G but with a somewhat different relation, which may lead to a modification of the results. But, further discussion of this is out of the scope of the present paper.
In this paper, we coupled a flocculation model and a sediment transport model to acquire insight into the impact of flocculation on large-scale sediment transport in a tide-dominated estuary. The combination of this flocculation model and the iFlow model allowed us to identify the relative importance of individual processes for water flow and sediment transport (e.g., river, tide, density forcing, advection). We employed a perturbation approach to gain insight into the highly complex coupled equations that describe flocculation. The result reveals a simple linear relationship between ws and c, both at leading and first order.
We applied our framework to a winter case in the Scheldt estuary. We showed that the spatial and temporal variations of ws due to flocculation are essential to reproduce observed suspended sediment concentrations. We were able to identify the impact of flocculation on individual transport mechanisms. We showed that flocculation alters most of the dominant transport mechanisms (e.g., river and tidal return flow) by introducing both a spatial and temporal dependence in ws.
To further investigate the impact of flocculation on large-scale sediment transport, we carried out a sensitivity analysis in which we showed that the magnitude of the floc breakup parameters kB and total freshwater discharge Q determine the existence, intensity, and number of ETM.
Although we found examples in literature of several other estuaries where flocculation is thought to promote sediment import, we were able to show that flocculation can also theoretically promote sediment export, depending on the phase difference between the tidal flow and suspended sediment concentration. Therefore, the effect of flocculation on the large-scale sediment transport in other estuaries needs to be assessed carefully based on local conditions.
We thank Tom Maris and De Vlaamse Waterweg NV for providing the OMES turbidity dataset, Tom Van Engeland for approval to use the ScheldeData package to generate the illustration of the Scheldt estuary (Fig. 3), and Emma Boegborn for her critical review of this paper. The iFlow, version 2.9, code with the flocculation extension and an input file example are available through Github (https://github.com/YoeriDijkstra/iFlow). Author Horemans is an SB Ph.D. fellow at FWO (1S36518N), with support of the University Foundation of Belgium.
Following Warner et al. (2005), we fit the observed salinity data to the following postulated salinity distribution of our winter case in the Scheldt estuary:
with ssea being the salinity boundary condition at the mouth and and being further undefined calibration parameters. Figure A1a shows the data fit. Edzwald et al. (1974) calculated the salinity dependence for stability value α′, which is defined as
where ϕ is the volumetric concentration, G is the shear rate, and N is the number concentration of flocs. Winterwerp and van Kesteren (2004) showed that
with fs being the floc shape factor and Df being the floc size. Consequently, in identification with Eq. (6) and assuming the floc shape is spherical (fs = π/6), we have
The latter relation allows us to use the salinity dependence of α′ for . We postulate the following salinity S dependence:
and we fit the latter function to the averaged data presented in Edzwald et al. (1974), which is shown in Fig. A1b. We acquire αmin = 0.0904, S2 = 4.085‰, and S3 = 0.034 01. We showed that net sediment transport only depends on the ratio of and kB and not on the individual parameters. We can thus decide to make solely a function of salinity S and keep kB fixed without restricting the generality of our results:
Sediment Equations with Ordering
where c is the suspended sediment concentration, u and w are the water velocity in the x and z directions, ws is the settling velocity, Kυ is the constant vertical eddy diffusivity coefficient, and Kh is the constant horizontal diffusivity coefficient. At the water surface, we require that no sediment particles enter or leave the domain:
with R being the reference level and ξ being the water surface elevation. At the bed, we require
where D and E are the deposition and erosion of sediment defined as
where M is the erosion parameter, τb is the shear stress at the bed, and f(a) is the erodibility.
To apply ordering and perturbation theory, we first have to compute the order of each term in the differential equation and boundary conditions. For this, we use typical scales for each variables following Chernetsky et al. (2010):
where the tilde denotes dimensionless variables, () is the order of magnitude, σ is the M2 tidal angular frequency, U is the typical scale of the horizontal velocity of the M2 tide, AM2 is the M2 tidal amplitude at seaward side, H0,Ems is the water depth at the mouth, Ws is the typical settling velocity scale, υ is the typical vertical eddy diffusivity coefficient scale, and LEms is the length of the Ems estuary. The typical scales following Chernetsky et al. (2010) are summarized in Table B1.
The boundary conditions in dimensionless form read as
where R is the water surface reference level, is the dimensionless water surface elevation, is the erosion scale, and C is the suspended sediment concentration scale.
To apply perturbation theory, we write the solutions for u, w, ξ, and c as a power series of a small parameter ε with (1):
where u0, w0, ξ0, and c0 are assumed to be of leading order; u1, w1, ξ1, and c1 are of order ε, and so on. Using the scaling, we acquire at leading order
and at first order
Scaling Analysis and Perturbation Theory
Similarly, we scale the Winterwerp flocculation model shown in Eq. (12):
To apply perturbation theory, we write the solution of ws as a power series of a small parameter ε with (1) following Chernetsky et al. (2010):
where is assumed to be of leading order, is of order ε, and so on. Using the scaling from Eq. (C1), we acquire at leading order in dimensional form the balance between floc aggregation and floc breakup:
with κ defined in Table 1. We added the second-order term −kBG1/2 from the breakup term to ensure that, in the limit for c0 → 0, the floc size is equal to the primary particle size: Df → Dp, which corresponds to a settling velocity of massive (nf = 3) primary particles:
We can add this (small) term to our leading-order result and leave out the βκ term in our higher-order calculations without restricting the generality of our results. The first-order equation yields
Calculation of Tidal Components of wS
In iFlow, we work with tidal components instead of time series. Therefore, we project the solutions of and on the subtidal, M2, and M4 tidal components. The amplitudes of the terms that are linear in c0 and c1 are trivial, for example,
with and (sub)tidal amplitudes.
The projection of the second term in Eq. (C6) is nontrivial. In general, an amplitude An of a periodic function func(t) can be obtained by
with T being the period of function func(t) and in which n denotes the (sub)tidal constituent (i.e., subtidal, M2, and M4). Substitution of yields
Application to the second term of the right-hand side in Eq. (C6), without the −2τ factor, and substitution of yields
We use the residue theorem
in which ζk is a pole of integrand func and m is the order of pole ζk. The integrand in Eq. (D5) has five poles:
Assuming , only the three poles and ζ = 0 lay inside the curvature of the integral and thus result in nonzero contribution (this assumption is not necessary to prove that the M2 constituent is equal to zero). Consequently, we have to calculate the residue for three poles to compute the integral in Eq. (D5). In the following, we show that the M2 tidal component is equal to zero by computing the residue for pole ζ = 0 and using the symmetry of the residue for the other two poles:
a. Pole ζ = 0
Calculating the residual for pole ζ = 0 yields
Consequently, the residue corresponding to the M2 constituent for pole ζ = 0 is equal to zero.
We compute the residue for pole :
and being the roots of the numerator:
The order m of the pole in the first term in Eq. (D15) is equal to 1. In the second term, the order is equal to 3. Consequently, the residues yield
in which we used that for pole the factor is replaced by to determine the symmetry of Z2. For n is odd (cf. M2, M6, M10, etc.), the solution in Eq. (D19) is asymmetric and consequently the sum of residuals is equal to zero. For n is even (cf. M0, M4, etc.) the solution is symmetric. Consequently, adding the solution for the pole ζ = 0, the final solution yields
which combined with the other trivial terms results in
Phase Requirements for Land-Inward Sediment Transport due to Flocculation
In this section, we demonstrate the phase conditions for which the M2 contribution of flocculation (solid black line in Fig. 6a) results in land-inward net sediment transport:
in which follows from Eq. (B16):
We assume that the water column is fairly well mixed, such that the vertical phase differences between u02 and c12 and c04 are negligible. We define the phases in polar coordinates as
So if Re(eiζ) > 0, that is, if −π/2 < ζ < π/2.
In the following, we link ζ to ϕ and ψ by integrating Eq. (E2):
in which we used the fact that only the right-hand side of Eq. (E2) generates an M2 tidal signal and
which results in
in which we used the fact that a −1 factor results in a π phase. Therefore, as a strong requirement for land-inward sediment transport due to the M2 contribution of flocculation, we require that
Likely, , so the first requirement may be leading.