Using an idealized width-averaged process-based model, the role of a mud pool on the bed and time-varying river discharge on the trapping of fine sediment is systematically investigated. For this purpose, a dynamically and physically motivated description of erodibility is presented, which relates the amount of sediment on the bed to the suspended sediment concentration (SSC). We can distinguish between two states: in the availability-limited state, the SSC is limited by the amount of erodible sediment at the bed. Over time, under constant forcing conditions, the estuary evolves to morphodynamic equilibrium. In the erosion-limited state, there is an abundant amount of sediment at the bed so that sediment pickup occurs at the maximum possible rate. The SSC is then limited by the local hydrodynamic conditions. In this state, the estuary keeps importing sediment, forming an erodible bottom pool that grows in time. These two states can be used to explain the response of an estuary to changing river discharge. Under availability-limited conditions, periods of high river discharge push estuarine turbidity maxima (ETMs) downstream, while drier periods allow ETMs to move upstream. However, when the estuary is in an erosion-limited state during low river discharge, a bottom pool is formed. When the discharge then increases, it takes time to deplete this pool, so that an ETM located over a bottom pool moves with a significant time lag relative to changes in the river discharge. Good qualitative agreement is found between model results and observations in the Scheldt Estuary of surface SSC using a representative year of discharge conditions.
The fine sediment distribution in estuaries is strongly influenced by the effects of climate change, such as accelerated sea level rise and intensified river discharge (e.g., Scavia et al. 2002; Robins et al. 2016; Achete et al. 2017), and human interventions, such as land reclamation, channel deepening, and channelization (e.g., de Jonge 1983; Winterwerp and Wang 2013; de Jonge et al. 2014; van Maren et al. 2015). In turn, the sediment distribution impacts the economic and ecological value of the system. For instance, a large (local) suspended sediment concentration (SSC) may contribute to the infilling of navigation channels and/or the siltation of harbors (van Maren et al. 2009), which likely results in an increase of dredging activities. Additionally, regions of high turbidity severely impact (local) light penetration and oxygen levels, resulting in a decrease of primary production (Cloern 1987; Talke et al. 2009; Liu and de Swart 2015).
Variations or changes in the sediment distribution occur at different time scales. For instance, at the decadal time scale channel deepening has caused some European estuaries to undergo a dramatic transition from a state with SSC up to a few hundred milligrams per liter to a hyperturbid state with maximum concentrations of tens of grams per liter. Examples include the Ems Estuary in Germany (Talke et al. 2009; Schuttelaars et al. 2013; de Jonge et al. 2014) and the Loire Estuary in France (Jalón-Rojas et al. 2016). At the annual and seasonal time scale, the river discharge regime induces variations in along-channel sediment distribution, where regions of high SSC, the estuarine turbidity maxima (ETMs), are usually most pronounced during low discharge periods (e.g., Castaing and Allen 1981; Woodruff et al. 2001; Lesourd et al. 2003).
Generally, the governing processes that facilitate these variations in sediment distribution consist of interacting vertical and horizontal processes. Vertical processes determine the ability of the water motion to bring into and keep the available sediment in suspension. Horizontal processes control the trapping of fine sediment by either redistributing the sediment already present in the system or by importing it from the sea, the upstream part of the river, or land. The relative importance of the various processes depends on the (time-varying) external forcings and the sediment availability (e.g., Burchard and Baumert 1998). For example, after a wet year in 1998, San Francisco Bay experienced a significant 36% decrease in SSC, which was hypothesized by Schoellhamer (2011) to be related to a depleted erodible sediment pool. Using a quantitative conceptual model, he confirmed that SSC can suddenly decrease when the threshold from an erosion-limited state to an availability-limited state is crossed, meaning that an erodible sediment pool is depleted. The terms erosion-limited and availability-limited are also referred to as Type II and Type I erosion, respectively (Mehta and Partheniades 1982). Type I erosion (cf. depth-limited, supply-limited, or availability-limited erosion) occurs if SSC is limited by the amount of sediment available for erosion. This occurs when the critical bed shear stress increases with depth into the bed, at the depth where the actual bed shear stress equals the critical bed shear stress. Type II erosion (cf. unlimited erosion) is related to an abundant amount of sediment available for erosion. Thus, SSC is limited by the maximum erosion rate induced by the hydrodynamic conditions.
The main aim of this study is to enhance our understanding of the trapping of fine sediment in estuaries under influence of various processes that act on different time scales. Specifically, we systematically investigate the role of erosion- and availability-limited conditions and time-varying external forcings. For this purpose, several types of process-based models can be employed ranging from exploratory to complex models (Murray 2003). Exploratory, or idealized, models typically aim at reproducing the major phenomena of the system under investigation by only taking into account a limited number of processes that are thought to be important. The key advantages of these models are their excellent ability to quickly investigate the sensitivity of the model outcomes to parameter variations and the possibility to systematically study physical processes in isolation. However, because of idealized nature of these models, the comparison with natural systems has to be qualitative and interpretation of the results requires careful consideration of the effects of the underlying assumptions. Examples of studies on sediment dynamics in estuaries using exploratory models are Burchard and Baumert (1998), Friedrichs et al. (1998), Chernetsky et al. (2010), de Jonge et al. (2014), Burchard et al. (2013), and Schulz and Umlauf (2016). On the other end of the modeling spectrum, complex models try to reproduce a natural system as closely as possible by implementing all known processes and state-of-the-art parameterizations. These types of models are the method of choice when a high level of accuracy is required. However, they are (presumably) computationally expensive and the results can be more difficult to interpret. Examples of the use of complex models in relation to sediment dynamics in estuaries and coastal seas are Wang (2002), Burchard et al. (2008), and van Maren et al. (2015).
Facilitating the aim of this study, we choose to employ the idealized width-averaged [two-dimensional vertical (2DV)] modeling framework iFlow v2.5 (Dijkstra et al. 2017). It is based on the model by Chernetsky et al. (2010) that was used to study the along-channel distribution of fine sediment in the Ems Estuary, in particular the locations of ETMs. An important assumption in their model is the existence of a morphodynamic equilibrium, that is, a divergence-free, tidally averaged sediment transport. While the approach by Chernetsky et al. (2010) has been successful in determining the physical mechanisms that constitute ETM dynamics in tidal estuaries, it has two drawbacks. First, the model assumes that the sediment concentration scales linearly with the total amount of available sediment, even when sediment is abundant. As a result, the SSC may exceed the carrying capacity of the flow. Second, the morphodynamic equilibrium approach cannot be used for changing external conditions (tides, river discharge), which may vary on time scales from days to years. The approach thus gives the sediment distribution for time-invariant external forcing.
To overcome these two drawbacks, we 1) introduce the concept of dynamic erodibility to set an upper limit for the SSC equal to the carrying capacity of the flow and 2) allow the external conditions and the amount of sediment available in the estuary to change with time. Both concepts are derived from first principles.
This paper is organized as follows. In sections 2 and 3 the adopted model and solution method will be outlined. This includes a brief review of the iFlow model (Dijkstra et al. 2017) and introduces the dynamic erodibility concept and the mass balance of the total amount of sediment in the estuary (both suspended and at the bed). Section 4 is devoted to model results with idealized river discharge time series to explain the concepts of availability- and erosion-limited conditions and to provide an interpretation framework. Subsequently, the sensitivity of the results to the dimensionless erosion parameter, which needs to be calibrated, is discussed in section 5. The model and interpretation framework are applied to the case study of the Scheldt Estuary in section 6. Finally, the main results are summarized in section 7.
2. Model description
The trapping of suspended sediments and its evolution on a long time scale will be investigated by extending the width-averaged, semi-analytical approach available in the iFlow modeling framework (Dijkstra et al. 2017). These extensions allow for the temporal evolution of the total sediment stock, which is defined as the sum of the amount of sediment in suspension and in the bottom pool available for erosion, thereby replacing the morphodynamic equilibrium concept with that of dynamic erodibility. In the following sections, the model equations for the water motion and sediment dynamics will be presented, with a special emphasis on the important extensions highlighted above. To aid the model description, a comprehensive overview of the model variables and parameters is presented in Table 1.
a. Model geometry and water motion
A schematized geometry of a single tidal channel of finite length L is considered, that is, no tributaries are allowed (Fig. 1). The width B(x) and bed level H(x) can vary gradually in the along-channel direction on a length scale comparable to the basin length. The water level ζ(x, t) horizontal velocity u(x, z, t), and vertical velocity w(x, z, t) follow from solving the Reynolds-averaged and width-averaged shallow water equations. We neglect the effects of Coriolis and assume that density variations are small compared to the average density, allowing for the Boussinesq approximation. We assume a hydrostatic balance and model the baroclinic pressure by a prescribed (i.e., diagnostic) time-independent salinity field that is vertically uniform. The eddy viscosity is assumed to be constant in depth and time, but is allowed to vary in the horizontal direction.
The seaward boundary is situated at x = 0, where a prescribed tidal water level ζe forces the water motion inside the basin. The water motion is only forced by M2 and M4 tidal constituents, resulting in
where , , and φ denote the amplitudes of the vertical M2 and M4 tide at the entrance and their relative phase, respectively. The quantity σ ~ 1.4 × 10−4 rad s−1 is the angular frequency of the M2 tide.
A weir is located at the landward end (x = L), where a time-dependent river discharge Q(t) is prescribed, such that
Here, denotes the tidal average and the minus sign reflects the fact that the river flow is pointing in the seaward direction. Note that condition (2) allows for a tidally averaged, slowly varying river discharge over the weir.
At the water surface, the usual kinematic boundary condition is adopted in conjunction with the assumption of vanishing shear stress. Moreover, the effect of wind shear stress is neglected. The bed is assumed to be impermeable, while the kinematic bed shear stress depends linearly on the velocity at the bottom (the so-called partial slip condition), that is
where the quantity sf is the friction parameter and assumed to be constant. Thus, the possible effects of sediment-induced stratification reducing bed shear stress are not incorporated in this model. The subscript denotes the derivative with respect to height. The vertical eddy viscosity Aυ is assumed to be vertically uniform and is related to the friction parameter and the local bottom depth according to (Dijkstra et al. 2017)
where H0 and Aυ0 = sfH0/2 denote the bottom depth and vertical viscosity at the entrance, respectively.
b. Sediment concentration
The width-averaged SSC c(x, z, t) evolves according to the advection–diffusion equation
where the subscripts and denote the derivative with respect to along-channel direction and time, respectively. Furthermore, ws is the settling velocity, which is assumed to be constant in space and time. The vertical diffusivity Kυ = AυσP, with σP the Prandtl–Schmidt number taken as 1 for simplicity. The horizontal eddy diffusivity coefficient Kh is constant.
At the seaward boundary, a depth and tidally averaged sediment concentration csea is prescribed. On the landward side, the tidally averaged sediment transport equals a prescribed fluvial import of sediment that may vary on a time scale that is long compared to the tidal period.
At the water surface, a zero sediment flux condition is applied. At the bed, the sediment flux follows from the difference between the instantaneous erosion E and deposition D. The deposition flux is modeled as D = wscbed, while the erosion flux reads (Dijkstra et al. 2018)
where cbed is the suspended sediment concentration just above the bed and is the amount of erodible sediment at the bed per unit area. The quantity is the potential erosion, that is, the maximum erosion flux given an abundant amount of erodible sediment at the bed. The first condition thus states that the actual erosion rate E is equal to the potential erosion when there is any sediment at the bed available for erosion. The second condition means that the erosion rate can only compensate for the deposition rate, if no sediment is available for erosion at the bed. For notational convenience, Eq. (6) is rewritten as
where denotes the instantaneous relative erodibility, hereinafter denoted as instantaneous erodibility, such that
Hence, varies between zero (no sediment at the bed and zero deposition) and one (sediment at the bed). The introduction of the instantaneous erodibility is the first novelty compared to earlier work from Chernetsky et al. (2010) and Dijkstra et al. (2017).
The potential erosion can be described by any existing erosion formulation and is usually expressed as , with M being an empirical constant, τb being the dynamic bed shear stress, and τc being the critical stress for erosion [see Sanford and Maa (2001) for an extensive discussion]. Among other definitions, in many studies τc is a function of erosion depth, representing a continuously increasing critical stress with depth or accounting for various bed layers (Mehta 2014, and references therein). As a consequence, both Type I and Type II erosion can be modeled. The depth dependency of the critical bed shear stress has been applied to several complex 3D numerical models dealing with sediment erodibility, such as the Zuidelijke NoordZee (ZUNO) model (van Kessel et al. 2011), the Community Sediment-Transport Modeling System (CSTMS) model (Warner et al. 2008), and the 3D Hydrodynamic Model for Applications at Regional Scale (MARS3D) model (Mengual et al. 2017). In this study, we take a simpler alternative approach to modeling Type I and Type II erosion. We assume a soft muddy layer (represented by ) on top of a sandy bed or consolidated mud layer (hereinafter called “sand”). Correspondingly, τc is zero for the muddy layer and τc ≫ τ for the sandy bed. Consolidation effects are not taken into account. Type I erosion is represented by the condition that , and thus . The condition that , that is, , represents Type II erosion, where the erosion rate is maximal based on the prevailing water motion. Following Huijts et al. (2006) and Chernetsky et al. (2010), we take , resulting in
In Eq. (9), ρs, ρ0, and ds denote the density of sediment and water and the grain size, respectively; is reduced gravity; and is the dimensionless erosion parameter, which needs to be calibrated. The dynamic bed shear stress τb(x, t) is defined as
where the partial slip condition (3) has been used in the last step.
c. Bottom pool evolution
Contrary to Chernetsky et al. (2010) and Dijkstra et al. (2017), we consider the temporal evolution of the sediment distribution for time-varying external conditions, which is the second novel aspect of this study relative to these previous publications using idealized models. With the time-varying water motion and sediment transport known, we can describe the time evolution of the active bottom pool. To this end, we use , of which the time evolution is governed by the difference between the (local) deposition D and erosion E that is
3. Solution method
a. Perturbation approach
Similar to Dijkstra et al. (2017) and Chernetsky et al. (2010), a perturbation approach is used to obtain an approximate solution to the full system of equations. This also enables the identification of relevant processes that underlie the global sediment balance. This approach involves 1) scaling of the governing equations, 2) asymptotic expansion of the physical variables, and 3) harmonic decomposition. Here, these steps are only shortly discussed; for a detailed description, see Dijkstra et al. (2017).
With the scaling method (step 1), the physical variables in the governing equations are made dimensionless using their typical scales. Consequently, the order of magnitude of the various terms in the governing equations is estimated. In the scaled equations a small parameter ε is identified, which is defined as the ratio of the prescribed amplitude of the M2 water level and the bottom depth at the entrance:
By relating the magnitude of all terms to ε, the order of magnitude of each term can be determined. It is further assumed that the ratio and the residual velocities due to river discharge and baroclinicity are all of order ε compared to the M2 tidal velocity. By expanding the physical variables into a series in ε (step 2), the full set of equations for the water motion and sediment transport is found to reduce to a coupled set of linear equations at different orders of ε, which are much easier to solve than the original problem.
The last step in the perturbation approach is the harmonic decomposition. All physical variables (water level, velocity, sediment concentration) vary both on the long subtidal and short tidal time scale. The leading-order linear equation for the water motion gives the dominant linear M2 tidal flow, while the first-order problem results in subtidal and M4 contributions due to nonlinear interactions of the leading-order tidal signal, the river discharge, density gradients, and the externally prescribed M4 tide. Neglecting contributions of second order and higher, the horizontal velocity and water level can be written in terms of harmonic components as
where subscripts denote the tidal constituent and superscripts denote the order of ε. The sediment concentration is written as
where now the residual and M4 contributions are the dominant leading-order quantities, while is of order ε. The quantities and denote temporal variations of leading- and first-order sediment concentration at lunar harmonic components and , respectively. Since these constituents do not result in contributions to the dominant tidally averaged sediment transport, they will be ignored hereinafter.
b. Tidally averaged erodibility
The sediment concentration depends on the harmonically decomposed instantaneous erodibility through Eq. (7). In the following, only the tidally averaged effect of on the long time scale is taken into account. The resulting tidally averaged erodibility, denoted by f(x, t), is defined as a weighted mean (for details, see appendix A),
where denotes the tidal average. Parameter is the local depth-integrated concentration at maximum instantaneous erodibility (i.e., ), defined by
Hence, can be seen as the subtidal carrying capacity, that is, the maximum amount of sediment that (on average) can be in suspension locally for given hydraulic conditions. Explicitly, a tidally averaged erodibility f =1 means the presence of an erodible, muddy bed layer throughout the entire tidal cycle, while f < 1 implies that the bed is sandy during at least part of the tidal cycle. Note that we do not account for sand transport or otherwise include sand dynamics in this model.
Now, using the tidally averaged erodibility f, we can write the erosion flux, Eq. (6), as
Thus, denotes SSC components that scale linearly with the erodibility f, that is, the subtidal and M4 components and a part of the M2 component. On the other hand, denotes SSC components that scale linearly with the longitudinal gradient of the erodibility fx, that is, a part of the M2 component. Note that and are completely determined by the water motion and sediment-related parameters. Using these expressions, the subtidal carrying capacity , defined in Eq. (17), can also be expressed as
c. Relation between erodibility and sediment stock
We now have an expression for the concentration in terms of f, but have not yet related the tidally averaged f to the amount of sediment at the bed. To do this, we assume that only the upper layer of the bed can exchange sediment with the water column through erosion and deposition, and any change of mass in this layer does not significantly change the bed level z = −H(x). Recall that we define the sediment stock as the total tidally averaged amount of sediment in the active bed layer and in the water column per unit area, that is
Defining the maximum leading-order deviation from [Eq. (17)] as
the ratio is a measure for the relative intertidal variability of the depth-integrated sediment concentration. Using and , an explicit functional relation can be derived between the erodibility f and the dimensionless stock (for details, see appendix A), resulting in the graphical representation shown in Fig. 2.
From this figure, it follows that for sufficiently small , such that
the stock is smaller than the approximated carrying capacity, that is, the amount of sediment that can be kept in suspension during almost the entire tidal cycle. In this case, the bed shear stress is almost always strong enough to erode all sediment from the active layer, resulting in a predominantly sandy bed. Since we assume that f only varies on the long time scale, this results in the approximation . For higher values of the dimensionless stock , such that
there will be mud in the active layer during some parts of the tidal cycle, resulting in a sublinear relation between f and [see Eq. (A8) in appendix A]. We denote this situation as a partly muddy bed. For large enough values of the dimensionless stock , , there will be mud at the bed throughout the entire tidal cycle. This means that the amount of sediment that can be eroded per unit time has reached a maximum and, by definition, f is equal to unity. This situation is considered a muddy bed.
We would like to stress that the dynamic erodibility concept developed here and the availability concept used in Friedrichs et al. (1998) and Chernetsky et al. (2010) are equivalent if the erodibility depends linearly on . Only if abundant amounts of sediments are available do the formulations differ, allowing for erosion-limited behavior when employing the dynamic erodibility formulation.
For future purposes, the notions of availability-limited (also called supply-limited, depth-limited, or Type I) and erosion-limited (Type II) circumstances are related to here. Availability is limiting if . In that case, the amount of sediment in suspension is uniquely related to the amount of sediment in the stock, and this relation exists because the bed is completely sandy during at least some part of the tidal cycle. The situation is referred to as erosion-limited since here the amount of sediment in suspension is at its maximum given the hydraulic conditions (carrying capacity). In the erosion-limited regime, there is no longer a unique relation between erodibility and the stock: all stock values above correspond to f = 1. Note that can vary along a tidal channel, and thus availability-limited and erosion-limited regions can be present simultaneously in an estuary.
d. Subtidal dynamics
To obtain an explicit expression for the unknown erodibility f(x, t), we consider the time evolution of the total tidally averaged amount of sediment. By integrating the concentration [Eq. (5)] over depth and width and using Eq. (21), the leading-order tidally averaged contribution to the bottom pool evolution [Eq. (11)] can be rewritten in terms of the total stock and the divergence of the subtidal sediment transport . The fact that the stock evolves on a time scale that is large compared to the M2 tidal period (see appendix B) justifies the use of the subtidal sediment transport. We thus have
The boundary conditions to Eq. (27) are specified in terms of the erodibility f(x, t). At the seaward boundary, the depth- and tidally averaged SSC csea is prescribed, resulting in fsea such that
At the landward boundary, the total net sediment transport equals the fluvial sediment input such that
where the minus sign on the right-hand side indicates that implies import of sediment from the riverine side.
In addition, an initial condition for f needs to be prescribed, with f being a nonnegative function that is smaller than one and obeys the boundary conditions. Equation (27) is solved for f by applying a backward Euler scheme for time integration, using a second-order upwind scheme for spatial discretization. This approach has second-order accuracy in space and first-order accuracy in time. The backward Euler scheme for time integration has been used as it provides better stability compared to second-order integration like Crank–Nicolson.
4. Model results
In this section, we present model results of the tidally averaged sediment distribution for different discharge time series. First, we consider a situation where there is hardly any sediment in the system initially and prescribe both high and low constant river discharges. This enables a thorough analysis of a system’s evolution under constant river discharge and under which conditions an equilibrium distribution of suspended sediment can occur. Explicitly, it allows an assessment in terms of the availability- and erosion-limited conditions discussed in section 3c. Next, we consider a discharge time series that smoothly alternates between high and low river discharge and focus on the behavior of the system during the transitions from high to low river discharge and vice versa.
For the experiments, we use parameter values representative for the Scheldt Estuary, which is described in more detail in section 6a. The values are listed in Table 2. Additionally, the along-channel profiles for the bed level H(x), width B(x), and salinity s(x) are taken from Dijkstra et al. (2017). We consider Q = 60 m3 s−1 and Q = 25 m3 s−1 as representative high and low river discharge, respectively, and assume no influx of fluvial sediment, that is, . Finally, we initially set the dimensionless erosion parameter to = 1 × 10−4.
a. Constant high river discharge
We start with a near-empty system by prescribing the initial condition finit = fsea(1 − x/L)20 and calculate the tidally averaged SSC after one year for a constant high river discharge (Fig. 3a). It follows that an ETM has formed around 100 km from the seaward side with a maximum concentration of approximately 160 mg L−1. The temporal evolution of the SSC toward this final distribution (Fig. 3b) shows that qualitatively the spatial sediment distribution has already formed within 25 days and that equilibrium is reached after approximately 125 days. Dijkstra et al. (2017) found that the dominant physical mechanisms constituting this ETM are the river flow, which transports sediment downstream, and sediment advection and tidal return flow (due to correlations between M4 velocity and M4 concentration) that mainly transport sediment upstream. The accompanying total sediment transport to form the ETM decreases in time to zero until equilibrium is reached (Fig. 3c). The total amount of imported sediment is only distributed throughout the water column (Fig. 3d, blue line), resulting in a depleted bottom pool during the entire simulation (red line). The final along-channel distribution of the erodibility (Fig. 3c, gray dashed line) is qualitatively similar to the distribution of the SSC (Fig. 3b). The dashed line in Fig. 3c indicates regions where , that is, predominantly sandy bed conditions. The maximum value of the along-channel erodibility for each time step, fmax(t), which is usually located in the vicinity of maximum , is always smaller than one (Fig. 3d, black line), indicating that the dynamics are availability-limited.
When availability-limited conditions apply throughout the entire estuary (Fig. 3), an equilibrium sediment distribution [; Eq. (25)] occurs if the residual sediment transport is divergence free (). This implies by virtue of the up-estuary boundary condition [Eq. (31)]. In the case discussed here, with , this means that equilibrium is approached when the net transport across the estuary entrance approaches zero (see Fig. 3c, cyan line). The corresponding equilibrium for is the morphodynamic equilibrium state discussed in Friedrichs et al. (1998) and Chernetsky et al. (2010). It can be proved ( appendix C) that for any a global availability-limited morphodynamic equilibrium is linearly stable, that is, any small perturbation of this equilibrium will disappear over time.
b. Constant low river discharge
Conducting the same experiment for a constant low river discharge, the final tidally averaged SSC distribution shows two ETMs (Fig. 4a). One ETM is very pronounced and located close to the weir [river kilometer (rkm) 150] with a maximum concentration of approximately 750 mg L−1, and the other ETM is less pronounced, located around rkm 110 with a maximum concentration of approximately 160 mg L−1. The latter ETM appears quite similar to the only ETM in the high river discharge case (see Fig. 3b). However, it turns out that the dominant importing mechanism for this ETM is velocity–depth asymmetry, with sediment advection having a secondary role [see Dijkstra et al. (2017) for a detailed description of the processes]. For the ETM close to the weir, both velocity–depth asymmetry and tidal return flow are important for sediment import. For both ETMs the river flow is still the dominant exporting mechanism. The suspended sediment distribution and mass show a similar evolution as for the high river discharge case and reach equilibrium after approximately 125 days (Fig. 4b, green line, and Fig. 4d, blue line).
In contrast to the high river discharge case, the amount of mass in the bottom pool (Fig. 4d, red line) starts to grow linearly in time after approximately 50 days. The onset of this linear growth coincides with the moment the maximum erodibility reaches a value of 1 (Fig. 4d, black line), indicating that erosion-limited conditions apply. The accumulating bottom pool in the erosion-limited region is supplied by an import of sediment from the seaward boundary that eventually reaches a constant value (Fig. 4c, cyan line). Notice that the final along-channel distribution of f (Fig. 4c, gray line) shows that predominantly sandy bed conditions (dashed gray line) coexist with (partly) muddy bed conditions (solid gray line). Seaward of the erosion-limited region f < 1 and thus availability-limited conditions apply. Therefore, sediment import becomes constant when an availability-limited equilibrium has set in (). On the landward side of the erosion-limited region f < 1 and , resulting in zero sediment transport. Note that if , the bottom pool in the erosion-limited region can additionally be supplied from the landward boundary.
To sum up, similar to the availability-limited equilibrium, the erosion-limited equilibrium is characterized by a tidally averaged suspended sediment distribution that does not vary in time (cf. blue lines in Figs. 3d and 4d). The amount of mass in the bottom pool (cf. red lines in Figs. 3d and 4d), however, grows in time.
c. Transition between high and low river discharge
To study the transition between high and low river discharge, we impose a discharge time series that smoothly alternates between high (Q = 60 m3 s−1) and low (Q = 25 m3 s−1) values (see Fig. 5a). The maximum gradients in the river discharge are indicated by the vertical dotted lines.
The simulation starts with the equilibrium SSC distribution corresponding to the high river discharge (see Fig. 3a), which shows a single ETM around rkm 100 (Fig. 5b). After the discharge has dropped below approximately 30 m3 s−1 around 225 days, a pronounced second ETM emerges close to rkm 150. When the river discharge transitions back to a high value, initially sediment concentrations significantly increase from rkm 80 to 155. Subsequently, the SSC values and along-channel distribution evolves back to those corresponding to the high discharge equilibrium.
The maximum along-channel value of the erodibility (Fig. 5c, black line) indicates that, in accordance with Fig. 3d, availability-limited conditions prevail everywhere in the basin at high discharges and no bottom pool has formed (Fig. 5c, red line). At low river discharges, erosion-limited conditions (f = 1) apply and we observe the characteristic linear growth of the bottom pool (see also Fig. 4d). Note that during the phases of constant high and low river discharge, the amount of SSC (Fig. 5c, blue line) reaches a constant value, indicating that both availability- and erosion-limited equilibrium situations occur.
Remarkably, the total amount of suspended sediment (Fig. 5c, blue line) shows a sharp peak when transitioning from low to high discharge, whereas from high to low discharge this is not the case. This is further explained in Fig. 6. For the transition from high to low discharge (Figs. 6a–e), the local erodibility increases (Fig. 6b) as the river discharge decreases (Fig. 6a). At the same time, the system adapts to the lower discharge by importing sediment over the seaward boundary (Fig. 6c) and depositing it at a second trapping location around x = 150 km (Fig. 6e). As soon as the new trapping location has reached the maximum possible value of the SSC (~235 days), and thus an erosion-limited condition sets in, the bottom pool starts to grow (Fig. 6d). After approximately 260 days an erosion-limited equilibrium is reached.
The transition from low to high river discharge (Figs. 6f–j) shows a different trend. During the period of low river discharge, the bottom pool has grown considerably and peaks around 540 days (cf. Figs. 6d and 6i). After that, as the discharge increases (Fig. 6f), the system adapts by transporting sediment from the bottom pool back to the downstream trapping location (Figs. 6i and 6j). Up to 565 days, this transport to the downstream trapping location is supported by sediment import over the seaward boundary (Fig. 6h). The total amount of suspended sediment peaks around 580 days, which coincides with the moment that the system turns from a state of local erosion limitation to a global state of availability limitation (i.e., fmax drops below unity). At this moment, the bottom pool is depleted and the excess amount of suspended sediment from the upstream trapping location is exported via the downstream trapping location to the sea. In short, it takes time to deplete the bottom pool when transitioning from erosion to availability-limited conditions.
5. Sensitivity to the erosion parameter
The dimensionless erosion parameter , via the potential erosion function [Eq. (7)], determines the amount of sediment that can be kept in suspension under the prevailing hydrodynamic conditions. Consequently, influences the occurrence of regions of erosion limitation: see Fig. 7a. In this figure, for the case of constant low river discharge (section 4b and Fig. 4), is varied between 0.05 × 10−4 and 1.1 × 10−4 and shows the corresponding along-channel erodibility f after one year. Confirming the results in Fig. 4, for = 1 × 10−4 (upper black dashed line) a single erosion-limited region occurs around rkm 150 (indicated by the area between solid light blue lines). Decreasing ultimately leads to the emergence of a second erosion-limited region around rkm 100–110. On the other hand, increasing further eventually results in the disappearance of the erosion-limited region around rkm 150 (not shown here). Notice that for most values the largest part of the bed of the estuary is predominantly sandy throughout the entire tidal cycle (area outside of the light blue dashed lines).
The corresponding scaled distribution and maximum value of along-channel tidally averaged surface SSCs are shown in Figs. 7b and 7c, respectively. In general, the sediment distribution follows the distribution of the erodibility: high surface SSCs are found where the erodibility approaches unity. Additionally, decreasing values result in decreasing values of the maximum along-channel surface SSC. Note that around a value of 0.2 × 10−4 there is a switch in the location where maximum concentrations are found. This is explained by the fact that the SSC is directly linked to the magnitude of the local bed shear stress and the erodibility via the erosion flux . As the bed shear stress is higher at the downstream trapping location than at the upstream one (not shown here), the maximum SSC is higher as well for equal f. Since the upstream trapping location already is erosion-limited with f = 1, the downstream trapping location gains importance compared to the upstream trapping location for decreasing values and consequent increasing f values.
The emergence of a second erosion-limited region around rkm 100–110 for values below approximately 0.1 × 10−4 (lower black dashed line) also shows an additional step in the sediment transport (Fig. 8). This additional sediment transport seaward of the most-downstream trapping location supplies the growth of the local bottom pool at that location. Any further seaward-located erosion-limited region will give rise to similar additional contributions to the sediment transport. All these additions are positive and add up to give either a further reduction of the down-estuary-directed fluvial sediment transport or an import of sediment at the seaward boundary.
6. Application to the Scheldt Estuary
In the previous sections, we have used an idealized model for the Scheldt Estuary forced by simplified time series for the river discharge. The model results give insight into the sediment dynamics and the role of availability- and erosion-limited conditions. Here, we will use this interpretation framework for a realistic river discharge time series that is based on data from the Scheldt Estuary. Our aim of the comparison is to reproduce trends and orders of magnitude of observed SSC to obtain a better understanding of the underlying physical mechanisms governing the sediment distribution under time-varying forcing.
a. The Scheldt Estuary
The Scheldt Estuary is a macrotidal estuary located on the border between Belgium and the Netherlands (see Fig. 9) and runs from Gentbrugge to Vlissingen. The seaward (Dutch) part is referred to as the Western Scheldt while the Belgian part is called the Sea Scheldt, which is divided in the Upper and Lower Sea Scheldt. The total length L of the estuary is approximately 160 km. At Vlissingen, the estuary is approximately 6 km wide and has an average depth of 15 m. The mean tidal range at the entrance is about 3.8 m and increases in the landward direction to reach a maximum of approximately 5.5 m at Antwerp (rkm 75; Vandenbruwaene et al. 2013). Landward from Antwerp, the tidal range decreases to ~2.7 m, where the estuary has converged to a width and depth of about 50 m and 3 m, respectively. Here, a weir is situated so that the tidal wave cannot propagate further upstream. At Melle, which is only slightly downstream from the weir, a discharge gauge is located that measures daily averaged runoff values (https://www.waterinfo.be). Median values in the period 1996–2016 range between 10 and 70 m3 s−1 (Fig. 10), while 25–75 percentile (PCTL) values fluctuate between 5 and 150 m3 s−1. The saltwater influence typically reaches up to the port of Antwerp and the salinity can be regarded as vertically well mixed.
Near-surface SSC in the Sea Scheldt have been monitored monthly along longitudinal transects from rkm 59 to 150 since 1996 and gathered in the OMES dataset (Maris and Meire 2016). These measurements are taken independent of tidal phase and thus cover many tidal conditions encountered in the Scheldt Estuary. From the data it follows that SSC are moderate, with near-surface concentrations only occasionally and locally (usually between 100 and 140 km) exceeding 200 mg L−1 (Fig. 11, blue lines and areas). A dumping area of dredged material is located at rkm 70–75, giving rise to unnatural elevated surface SSC (Fig. 11, gray areas).
The parameter values adopted for the current Scheldt model are taken from Dijkstra et al. (2017) and are summarized in Table 2. Width and depth profiles are obtained by fitting smooth functions through the 2013 data (Coen et al. 2015). Furthermore, it is assumed that the horizontal salinity profile is stationary and depth independent and can be described by a tangent hyperbolic function. The model has been calibrated against M2 and M4 water levels by optimization of a cost function that varies with the hydraulic roughness sf, which in turn determines the vertical eddy viscosity Aυ0 (Dijkstra et al. 2017). Following Coen et al. (2015), the settling velocity of the sediment is 2 mm s−1. The depth- and tide-averaged SSC concentration at the entrance, csea, is set to 40 mg L−1. The sediment input from the river is set to zero, that is, . Finally, the dimensionless erosion parameter is chosen to match observations of surface SSC as well as possible, resulting in = 2 × 10−5.
b. Results and discussion
Because SSC data are measured at random moments in the tidal phase only once every month, we define a representative year based on both the median value of daily discharge data and monthly SSC data from 1996 to 2016. Additionally, we use the 25th–75th PCTL values to highlight the variability of the data. Using the median river discharge time series (Fig. 10a), we start the simulation in morphodynamic equilibrium and run two consecutive years. The computed tidally averaged distribution of surface SSC for the second year and the comparison with observed monthly median SSC from the OMES dataset are depicted in Figs. 10b and 11, respectively. The variability of the tidally averaged model results is denoted by the leading-order M4 concentration [see Eq. (15)]. Apart from the dumping area, the model results match the observations remarkably well. Both the location of the ETM between rkm 100 and 120 and the magnitude of the surface SSC are within reasonable ranges. However, the flushing of the ETM, for example, from January to March, is less well captured by the model. This is, in part, due to the fact that by taking the median values of the river discharge data, major flushing events or periods that push the ETM downstream are filtered out. Additional model runs of individual years with a relatively high river discharge regime (not shown here) revealed that transient behavior of the ETM location is better captured.
During the simulation, two bottom pools are continuously present where the dynamics are erosion-limited (Fig. 12a, red areas). One is located close to the weir, between rkm 140 and 160, and one more downstream between rkm 100 and 120. Furthermore, upstream from rkm 70 the bottom composition is muddy during at least part of the tidal cycle. Downstream of rkm 70 the bottom is predominantly sandy.
The total mass of sediment suspended in the water column is almost constant (Fig. 12b, blue line), while the total mass in the two bottom pools fluctuates considerably under the influence of the river discharge (Fig. 12b, red and green line). The bottom pool close to the weir (green line) follows the discharge trend and empties during periods of relatively high river discharge, but never depletes, and fills up during periods of relatively low discharge. In turn, the interior bottom pool (red line) fills up during periods of high river discharge when sediment is pushed downstream and empties during periods of low river discharge, but at a much slower rate. The mass balance is closed by either an additional import or export of sediment through the seaward boundary (black line), which follows the same trend as the river discharge, that is, import during low discharge and export during high discharge.
7. Summary and conclusions
In this paper we have systematically investigated the role of erosion- and availability-limited conditions and time-varying river discharge controlling the trapping of fine sediment in estuaries. To facilitate the investigation we have used the idealized, process-based iFlow modeling framework (Dijkstra et al. 2017), where both vertical and temporal variations of vertical viscosity and diffusivity have been ignored and the bed shear stress is taken to scale proportional with the near-bed velocity (Lorentz linearization). We have utilized a dynamically and physically motivated description of erodibility that extends the availability concept that was used in similar models by Friedrichs et al. (1998) and Chernetsky et al. (2010). Hereby, the role of a critical shear stress for erosion is not considered.
The erodibility f(x, t) varies on the long subtidal time scale. It is related to the sediment stock , which measures the total amount of sediment in suspension and material in the active bed layer that is available for erosion. When is low, almost all sediment is in suspension. As is increased, the bed may become muddy during parts of the tidal cycle. These two situations are referred to as availability-limited conditions, which give concentrations that are limited by the amount of erodible sediment at the bed. In this state, sediment transport is in morphodynamic equilibrium. The SSC in availability-limited conditions tends to, but not necessarily, reach the morphodynamic equilibrium distribution that was already studied in Chernetsky et al. (2010). In the current contribution, it is added that such equilibria are linearly stable. For higher values of , the bed is muddy throughout the entire tidal cycle so that sediment pickup occurs at the maximum possible rate and erosion-limited conditions apply. The SSC is then limited by the local hydrodynamic conditions, while the mass of the bottom pool increases in time. This situation is also stable in that the SSC does not vary on a subtidal time scale. The growth of the bottom pool requires a net import of sediment through at least one of the boundaries.
In general, during periods of high discharge an ETM is pushed downstream, while during periods of low discharge an ETM moves upstream. For sufficiently low discharges, both availability- and erosion-limited regions will coexist within an estuary. Consequently, in erosion-limited regions an erodible bottom pool will form. When the river discharge increases again, an ETM located at an erosion-limited region maintains its high concentrations until the bottom pool is depleted and thus shows a time lag relative to changes in river discharge. At the same time, the formation of a new turbidity maximum elsewhere in the estuary requires sediment being transported to that location, which might result in an additional sediment import over the system boundaries.
It is further found that the instantaneous horizontal sediment distribution shows good qualitative agreement with the equilibrium distribution. Only when the system transitions from a state of local erosion limitation to global availability limitation is the equilibrium distribution not representative for the instantaneous distribution.
The model is applied to the Scheldt Estuary using a representative year based on median values of a 20-yr discharge and SSC dataset. Despite the adoption of some simplifying assumptions, model results show good qualitative agreement with observations of near-surface SSC. During periods of low river discharge the system tends to import sediment, while during periods of high river discharge it exports sediment.
Ronald Brouwer is funded by the VNSC (http://www.vnsc.eu) Agenda for the Future scientific research program through Contract 3109 6925 and Yoeri Dijkstra is funded through Contract 3110 6170. The authors wish to thank The Flemish Waterway for permission to use the OMES dataset. The authors thank Wouter Vandenbruwaene and Tom Cox for their constructive discussions regarding the OMES data set of the Scheldt Estuary. The authors thank Carl Friedrichs and an anonymous reviewer for their constructive comments.
Relationship between Erodibility and the Sediment Stock
In this appendix, a closed functional expression for the erodibility f(x) will be derived. The starting point is Eq. (21), which is repeated for clarity:
Here, and c may vary on the short (intertidal) time scale. Furthermore, the stock is a conserved quantity on the short time scale, since it can only change by divergence of the horizontal transport , which evolves on a long time scale (see appendix B). Next, we rewrite the depth integral of the sediment concentration in Eq. (A1) by considering only sediment erosion by the dominant M2 tide:
where is the instantaneous erodibility, is the subtidal carrying capacity [Eq. (17)], (n = 1, 2, …) are the even harmonic components of the carrying capacity [ is defined by Eq. (22)], and are the corresponding phases.
where is the relative stock at the bed and
is a nondimensional measure of the maximum amount of sediment in suspension at a given time. Alternatively, can be interpreted as the dimensionless sediment carrying capacity including temporal settling lag. Note that this capacity has a nonzero minimum, which means that for sufficiently low stock values all sediment will be in suspension for almost the entire tidal cycle. The last term in Eq. (A3) yields, by definition, f by tidal averaging, that is
To explicitly calculate f, we approximate with only the subtidal and M4 harmonic component (see Fig. A1, black line). Now, consider the case , for example, the blue line in Fig. A1. The approximated carrying capacity is larger than the stock. Hence, almost all material is in suspension and the erodibility follows from
We define this situation as a predominantly sandy bed.
On the other hand, if , for example, the green line in Fig. A1, the total amount of sediment is always larger than the carrying capacity. Therefore, is positive at any time during the tidal cycle and the erosion is maximum (), that is
This situation is considered as a muddy bed.
In the intermediate regime , for example, the red line in Fig. A1, the bed will be muddy only during parts of the tidal cycle, in the intervals and (red hatched areas), and the mud fraction will increase as increases. The tidally averaged erodibility f is then obtained as
where and P is the tidal period. The relation for this intermediate regime is then found to be given by
While Eq. (A8) looks somewhat involved, it merely represents a smooth sublinear increasing transition between a predominantly sandy bed () and a fully muddy bed (f = 1).
To summarize, is given by
It should be noted that and in general depend on location x and, thus, the transition points and between the regimes will vary throughout an estuary.
Estimate of the Sediment Stock Evolution Time Scale
In this appendix an estimate will be given for the time scale on which the amount of sediment in the stock varies. To this end, the net sediment transport is estimated as
where , , , and denote typical values for width, depth, tidal velocity, and sediment concentration, respectively. The along-channel variation of sediment transport is assumed to occur on the length scale L. Using these estimates in Eq. (25) yields
where is an estimate for the amount of sediment in the stock. Comparing the tidal time scale P (i.e., the tidal period) with the typical time scale at which varies, we find that
The term between parentheses is the ratio of depth-integrated concentration to the amount of sediment in the stock, which is typically of order one or less. Using the depth-integrated mass balance, it can be shown that (Schuttelaars and de Swart 2000). It thus follows that
and therefore the sediment stock evolves on a time scale that is large compared to the tidal time scale. This justifies the use of tidally averaged sediment transport when considering the evolution of the sediment stock in Eq. (25).
Linear Stability of the Availability-Limited Equilibrium
Here, it will be shown that the availability-limited equilibrium solution is a stable static solution to Eq. (C1), that is, any small perturbation of this solution will be damped in time. In the availability-limited regime there is a unique (and invertible) relation between f and the amount of sediment in the stock . The time evolution of [Eq. (27)] can then be re-expressed in terms of f alone as
An equilibrium sediment distribution (ft = 0) occurs if vanishes throughout the estuary, which implies by virtue of the up-estuary condition (31). This is the morphodynamic equilibrium state discussed in Friedrichs et al. (1998) and Chernetsky et al. (2010), albeit with a possibly nonzero tidally averaged sediment transport. The corresponding equilibrium erodibility function, hereinafter denoted by fmeq, is then obtained by solving , which results in
Next, consider the perturbed solution to Eq. (C1),
We then find that any small perturbation f′ obeys the following time evolution equation:
Note that I(t) is always positive, unless f′ vanishes everywhere. Also, F(x) is always negative so that the right-hand side of Eq. (C6) is negative, unless again f′ is identically zero. Hence, it follows from Eq. (C6) that dI/dt < 0, which implies that f′ → 0, and, therefore, any small perturbation of fmeq will damp out and the morphodynamic equilibrium is stable. Note that this argument holds for any relation that is a monotonically increasing function.
The stability problem can be used to estimate the adaptation time tadapt. To investigate this further, we write the perturbation f′(x, t) as
We the find that Eq. (C3) can be written as
where fμ(x) obeys boundary conditions (C4) and (C5). Equation (C8) is a so-called Sturm–Liouville eigenvalue problem. This means that there are infinitely many eigenvalues μ [and corresponding eigenfunctions fμ(x)], which obey Eqs. (C4), (C5), and (C8). Furthermore, the eigenfunctions form a complete set, implying that any arbitrary perturbation f′(x, t) is a superposition of these eigenfunctions. Moreover, the eigenvalues are real and can be ordered as follows
Consequently, there exists a smallest eigenvalue that, since the equilibrium is stable μ1, must be positive. Any eigenfunction is thus found to decay in time according ~exp(−μt). From Eq. (C9) it follows that the first eigenfunction is the slowest decaying one. Therefore, this eigenmode will dominate the temporal behavior toward the morphodynamic equilibrium. Hence, the inverse of μ1 can be viewed as the time scale on which the equilibrium will settle in, the so-called adaptation time scale tadapt. We thus have
Publisher's Note: This article was revised on 8 February 2019 to include an addition in the Acknowledgments section that was omitted when originally published.
The notation used in the erosion flux is slightly different from the one used in Friedrichs et al. (1998) and Chernetsky et al. (2010). In the two papers, the (possibly) spatially varying availability, denoted as a(x), is contained in the reference concentration . In the formulation used here, this spatial dependency (now captured by f) is not absorbed in , but made explicit in the factor f. Hence, to compare the results between the various papers, one can observe that , but only when a = 1 in Eq. (12) in Chernetsky et al. (2010).