The hyperturbid state of the water column in estuaries and rivers: the importance of hindered settling

Over the last few decades, some estuaries have undergone a transition to a hyperturbid state, characterised by suspended sediment concentrations of several grammes per litre averaged over the water column. To improve our understanding of this transition and of naturally hyperturbid estuaries, we systematically identify the processes allowing for high suspended sediment concentrations using a water column (1DV) model. Under a range of realistic forcing conditions, the state of the water column can be characterised by one of two equilibrium states. The first is an erosion-limited state, in which there still is sediment available for erosion at the bed. We find that this state only occurs with relatively low concentrations. The second is a supply-limited state, in which all erodable sediment is in suspension. The concentration in this state depends entirely on the amount of sediment in the system and can potentially be very high. We identify the conditions under which the state of the water column can jump from a low to a high concentration and identify hysteresis in the transition between the two states. The mechanism responsible for this hysteresis is hindered settling. It thus follows that hyperturbidity is only possible in a supply-limited state. From this observation we derive a necessary condition for an estuarine system to make the transition from low turbidity to hyperturbidity in a 1DV context. This is an important step towards understanding why some estuaries are hyperturbid and assessing the risk that particular estuaries may become hyperturbid in the future.


Introduction
Suspended sediment concentrations in estuarine turbidity maxima range between fairly low suspended sediment concentrations of a few tens of mg/L averaged over the water column, to hyperturbid levels of several 10 g/L (Uncles et al. 2002). Apart from large variations in sediment concentrations between estuaries, there is also a large temporal variability of the sediment trapping location. The timescale of such variations ranges from the tidal time-scale up to several years (Jalón-Rojas et al. 2016). Concerning the variability on a decadal time-scale, it has recently been highlighted that concentrations have increased significantly in some estuaries. Examples of such estuaries are the Ems River, where the average concentration at the surface has increased from 200 mg/L in the 1950s to over 1 g/L now (e.g. Talke et al. 2009;Schuttelaars et al. 2013;De Jonge et al. 2014), and the Loire River, where the average concentration near the surface has increased from around 500 mg/L in the 1970s to several grammes per litre now (Jalón-Rojas et al. 2016). This dramatic increase in suspended sediment concentration has a severe negative impact on light penetration and oxygen conditions, resulting in a strong reduction in primary production (e.g. Cloern 1987;Talke et al. 2009).
There are strong indications that this long-term increase in sediment concentration is related to ongoing human interventions, including removal of intertidal area, deepening of navigation channels and continuous dredging to maintain the channel at navigation depth. Several studies have made this connection by using the combined knowledge of measurements and model results (Chernetsky et al. 2010;De Jonge et al. 2014;Van Maren et al. 2015). A possible feedback mechanism that could result in such a transition was proposed by . According to their hypothesis, channel deepening results in tidal amplification and more import of fine sediment. This sediment leads to a reduction of the drag and therefore more tidal amplification. This describes a feedback process, which may induce a strong increase in concentration, in response to a relatively small deepening. However, to date, no study has been able to validate this hypothesis by modelling the transition to a hyperturbid state as a consequence of one or multiple human interventions. Neither are the physical mechanisms featured in the hypothesis fully understood. A better understanding of these processes is important to make well-founded management decisions in estuaries that face substantial channel deepening.
To model the transition to hyperturbidity, at least two requirements have to be met. Firstly, it requires the availability of fine sediment. This typically requires an import of sediment into the system from the sea, the upstream part of the river or from land. Secondly, once a sufficient amount of sediment is available for erosion, the water motion should be able to bring and keep this sediment in suspension.
In this paper, we will focus on the second requirement: the ability of the flow to bring and keep sediment in suspension. The aim of the paper is to gain a better understanding of the processes that allow for high suspended sediment concentrations in the water column. To this end, the physical processes relating to erosion and settling of sediment in a water column (1DV) model are systematically investigated. Although we make no explicit comparison to field data, our parametrisations for erosion and settling are (semi-)empirical formulations based on laboratory and field measurements. Furthermore, the sensitivity of the results to the choice of parametrisations, parameter values, and external forcing conditions is analysed, so that the results should be globally applicable.
A lot of theoretical and experimental work has been done on individual processes acting on sediment in the water column, especially related to high concentrations or fluid mud (see e.g. Mehta 2014, for a review). The combined effect of many of these processes on sediment concentrations in water column models has been studied by for example Winterwerp (2001) and Le Hir et al. (2001). These studies have mainly focussed on stratification of the water column and the formation of lutoclines, identifying sediment-induced turbulence damping and hindered settling as the most important processes governing the amount of stratification and the formation of fluid mud from particle settling. By including these processes, they show that water column models are able to reproduce much of the behaviour of suspended sediment stratification observed in estuaries. Given this knowledge, we will not focus much on lutoclines and the structure of the suspended sediment concentration in the water column. Rather, we will focus on the overall magnitude of the concentration near the bed, while accounting for stratification of the water column in our model. The water column model is introduced in Section 2, focussing on the formulations for erosion and settling. We will then discuss our main result in the context of stationary flows in Section 3 and extend this result to tidal flows in Section 4. We will qualitatively discuss the consequences of different model formulations for erosion and settling in Section 5. Finally, the main findings are summarised in Section 6.

Model equations
The water motion and sediment dynamics are described by the momentum and continuity equations and the sediment mass balance equation. We concentrate on the vertical profiles of the velocity and sediment concentration. To this end, we use a water column model similar to the one used in Winterwerp (2001). The focus on vertical processes, thereby ignoring horizontal gradients, is a reasonable approximation in many estuarine systems, since vertical exchange processes are typically dominant up to leading order (see, e.g. the scaling analysis 97 in the models of Chernetsky et al. 2010;Dijkstra et al. 2017).
In this section, we will focus on the physical formulations used in the model; for details on the numerical implementation, the reader is referred to Winterwerp and Van Kesteren (2004). We assume a hydrostatic flow in a water column with vertical coordinate z (positive-upwards) that varies between the bed at z = −H and a fixed surface at z = 0 (i.e. rigid lid approximation). The flow velocity u is assumed to be uni-directional and the effects of the Earth's rotation are neglected. The vertical density differences are assumed to be small compared to the actual density, allowing the use of the Boussinesq approximation. The resulting Reynolds-averaged momentum equation reads with boundary conditions where g is the acceleration of gravity, A ν is the eddy viscosity, τ is the bed shear stress and ρ 0 is a reference density for water. The subscripts z and t denote derivatives with respect to space and time respectively. The model is forced by a prescribed water level gradient ζ x , which is either constant or varies in time. The bottom boundary condition is further rewritten using the definition τ = ρu * |u * |, where u * is the bed shear velocity. The bed shear velocity follows from the flow velocity and bed roughness by assuming that the flow near the bed has logarithmic profile according to with roughness height z 0 and Von Kármán coefficient κ = 0.4. The eddy viscosity A ν is computed using the k − model and depends on the flow velocity profile and sedimentinduced buoyancy destruction (see Dijkstra et al. 2016 for details on the numerical implementation of the k − model used here). The sediment-induced vertical density gradient is related to the sediment concentration through a linear equation of state ρ = ρ w + c(1 − ρ w /ρ s ), where c is the sediment mass concentration and ρ w and ρ s are the densities of water and dry sediment respectively.
The sediment is assumed to consist of a single mud fraction, so that the sediment dynamics is described by with boundary conditions where w s is the settling velocity of sediment flocs and K ν is the vertical eddy diffusivity. The eddy diffusivity is related to the eddy viscosity through a constant Prandtl-Schmidt number as A ν = K ν σ ρ , with σ ρ = 2 (Van Maren et al. 2009). At the bed, sediment deposits at a rate D and erodes at rate E, which are described by The deposition rate is defined as the flux of sediment settling out from the water column on the bed. It therefore equals w s c bed , where c bed is the concentration suspended in the water column just above the bed. In order to specify the erosion rate, we first define the sediment availability. This is denoted by the sediment stock S (in kg per m 2 surface area), which is the sum of the mass of sediment per m 2 at the bed available for erosion, called S bed , and the mass of sediment in the water column, i.e. S = S bed + 0 −H c dz. Next we define the potential erosionÊ as the erosion rate provided that there is enough sediment at the bed available for erosion. The potential erosion is normally described by (semi-)empirical formulations, such as Eq. 9, which is referred to as Partheniades' formula (Kandiah 1974), based on experimental work by Partheniades (1962Partheniades ( , 1965. In this expression, M is an erosion parameter and τ c is the critical shear stress that needs to be exceeded for erosion. For simplicity, it is assumed M and τ c are uniform over the depth of the sediment on the bed. The consequences of relaxing these assumptions are explored qualitatively in Section 5.
When S bed equals zero, there is no sediment at the bed that can be eroded and the erosion rate cannot be equal to the potential erosion, unless the potential erosion equals zero. This simply follows from the principle of mass conservation, i.e. the amount of sediment at the bed cannot become negative. At maximum, the erosion rate can compensate for the deposition rate, so that deposited sediment is re-suspended immediately.
We will refer to conditions where S bed > 0 as erosion limited, as the erosive strength of the flow limits the maximum erosion rate in this state. In literature this is sometimes referred to as erosion rate limited. The condition S bed = 0 is referred to as supply limited, as it is the sediment supply that limits the erosion rate (e.g. Scully and Friedrichs 2007;Winterwerp et al. 2012). In literature, this is alternatively referred to as depth limited, expressing that the erosion has reached a depth below which sediments are too consolidated to be suspended given the present flow conditions. In this research, we will show that these two states lead to a clearly different behaviour of the water column. Which of the two states a water column is in depends dynamically on the flow and the parameters in the erosion model. In our model simulations, we prescribe the sediment stock S. Whether the model is in an erosion-or supply-limited state given this stock follows as a model result.
When the concentrations in the water column are high, interactions between sediment particles and the ambient water reduce the effective settling velocity. To account for these effects, we use the hindered settling formulation proposed by Richardson and Zaki (1954). Their formulation is based on the reasoning that highconcentration suspensions increase the drag exerted on particles by the water that flows through the narrow space between the settling particles. Their parametrisation reads where w s,0 is the settling velocity of a single sediment floc and φ the volumetric concentration of flocs defined as c/c gel . The gelling concentration c gel is the concentration at which the sediment flocs form a self-supporting network (i.e. fluid mud). The parameter m is determined from experiments and assumes values between 2.4 and 4.7 for coarse to fine particles. The value m = 5, often used in practice for fine sediment, is taken as the default value in this article. The effect of other choices for m will be discussed in Section 5.

Dimensionless erosion parameter
It is illustrative to consider a stationary flow and assume there is an abundant sediment supply (i.e. S = ∞, S bed = ∞). Under these assumptions, u t = 0, c t = 0 and E =Ê.
Integrating the sediment concentration (4) between z = −H and 0, it follows that deposition balances the potential sediment erosion, i.e. D =Ê. Substituting (7) for D and Eq. 9 forÊ, together with Expression (10) for hindered settling and using φ = c/c gel , the condition D =Ê can be written as which relates the near-bed concentration φ bed (on the lefthand side) to a quantity that is defined as the dimensionless erosion parameterẼ. In the results presented in the next sections, we will use the dimensionless erosion parameter to characterise the possible equilibrium states of the water column and show that this parameter is also useful in a context of non-stationary flows.

Abundant sediment supply
We first assume a stationary flow with an abundant supply of sediments, so that Eq. 11 definingẼ holds. This expression is a non-linear algebraic equation for φ bed that only depends on the bed shear stress τ and a number of model parameters. The bed shear stress follows by vertically integrating the momentum balance (1) and using the corresponding boundary conditions. It then follows that τ = ρ 0 gH ζ x . As the model is forced by ζ x , τ is known a priori. Hence, φ bed can be obtained by resolving nearbed processes, without solving for the entire water column. Thus, φ bed does not depend on turbulence in the water column. Figure 1 shows the resulting equilibrium near-bed concentration φ bed as a function ofẼ. We distinguish between two regions in this graph. In region I (0 < E <Ẽ crit ), two solutions for the near-bed concentration exist. However, linear stability analysis shows that only the solution depicted by the solid line is a stable solution, meaning that a system close to the equilibrium state evolves towards its equilibrium over time. The near-bed 16. These concentrations can be high (up to 16 g/l assuming a gelling concentration of 100 g/l), but are significantly smaller than the gelling concentration. The dashed line depicts an unstable solution, meaning that a system close to, but not exactly in this state, will move away from it over time.
The behaviour of the solution if it is not in equilibrium is best seen from the depth-integrated concentration (4) 0 Left of the equilibrium curve D > E and the depthintegrated concentration decreases over time. Conversely, on the right of the equilibrium curve, E > D and the concentration increases over time. Correspondingly, if the near-bed concentration is above the unstable equilibrium (dashed line), the system will continue to erode sediment and, within this model, the concentration will increase without bound.
In region II (Ẽ >Ẽ crit ) we are on the right of the equilibrium curve and the concentration continues to increase, regardless of the initial concentration. This is possible because of the effects of hindered settling: as the concentration increases due to erosion, hindered settling leads to a decrease in the settling velocity and therefore in the rate of sediment deposition. As a consequence, the net erosion rate of the system increases leading to a further increase of the concentration and reduction of the deposition.
The change in model behaviour as one moves from region I into II is known mathematically as a bifurcation. This bifurcation exists, because the left-hand side of Eq. 11, which represents a dimensionless deposition flux w s c w s,0 c gel , has a maximum. If φ is small, the deposition flux increases with φ, as the settling velocity is virtually constant and the sediment concentration increases. However, if φ is large, hindered settling has such a strong effect in decreasing the settling velocity, that the total deposition flux decreases with increasing φ.

Limited sediment availability
By limiting the sediment stock, the suspended sediment concentrations remain finite for all time. The equilibria are shown in Fig. 2. The equilibrium branches in Section 3.1 are still present and are shown in blue. Added to this are two new (orange) branches which result from setting a finite stock S of 15 (lower branch) and 150 kg/m 2 (upper branch). On these branches, all available sediment is suspended in the water column, preventing a further increase in concentration. This equilibrium is therefore referred to as supply limited, corresponding to the definition of supply limited given in Section 2 for the erosion rate. Naturally, a larger value of the stock leads to larger values of φ bed . Regardless of the stock, the supply-limited solution exists for all values ofẼ in region II, so that the concentration is bounded for all model settings. It also extends into region I, where it ceases to exist when it crosses the erosion-limited solution. At this point, the flow velocity and effect of hindered settling are no longer sufficient to keep all the available sediment in suspension and the concentration reduces towards the stable blue branch. On the blue branch, some of the available sediment remains on the bed and the suspended sediment concentration is restricted by the ability of the flow to erode and keep sediment in suspension. This is referred to as an erosion-limited equilibrium, also corresponding to the definition of erosion limited in Section 2.
Contrary to the erosion-limited case, the supply-limited solution cannot be computed from the near-bed balance  Fig. 1, with the addition of supplylimited branches for two different values of the stock, indicated by the orange lines. The lower orange line corresponds to S = 15 kg/m 2 , the upper orange line corresponds to S = 150 kg/m 2 . On these lines, all available sediment is suspended in the water column. The arrows illustrate how a system can jump from one branch to another for increasing E (from point a to b) and decreasingẼ (from b to c and further) E = D. This is because E reduces to D in this state (see the definition of erosion (8)) and the condition E = D reduces to the trivial w s c bed = w s c bed . In this case, the restricting condition reads S = 0 −H c dz, i.e. all sediment is suspended in the water column. This means that the near-bed concentration depends on the distribution of the sediment over the water column and thus on both the turbulence profile and hindered settling. Hence, the near-bed concentration follows from the water column model. The parameters for this model are assigned default values, which are presented in Table 1. The model uses 400 numerical grid cells and a time step of 5 s. We have conducted 40 model experiments with different values for the water level gradient. We start with the largest water level gradient and run the model for a sufficiently long time, until a stationary state is attained. Next, the water level gradient is decreased and the model is run again, using the result of the previous experiment as initial condition. This ensures that the solution remains on the supply-limited branch provided this branch still exists.
The supply-limited solutions in Fig. 2 are the orange down-sloping curves. The near-bed concentration decreases with increasingẼ because we vary ζ x in the model experiments. This not only leads to a variation inẼ via the bed shear stress, but also to a variation in the rate of mixing of the water column. As a result, a largerẼ coincides with a more well-mixed water column. As the same amount of sediment is in suspension for all solutions on the branch, the near-bed concentration decreases asẼ increases. This is also illustrated by the vertical sediment profiles plotted in Fig. 3, where the profiles from light to dark colours are obtained by increasingẼ from 0.015 to 0.15.
The distinction between erosion-and supply-limited equilibrium states has profound implications for hyperturbidity. Near-bed concentrations associated with hyperturbidity are of the order of the gelling concentration. However, the maximum near-bed concentration in a stable erosionlimited state is only approximately φ bed = 0.16 according to this model, one order of magnitude smaller than the gelling concentration. In order for a system to be hyperturbid, it therefore needs to be supply limited. A water column can only be supply limited with φ bed > 0.16 if it satisfies the Clear-water settling velocity 2 · 10 −3 m/s criterionẼ >Ẽ crit , or has satisfied this criterion at some point in its history and has stayed on a supply-limited equilibrium branch with high concentrations (see below). Figure 2 also has profound implications for the ways in which a transition between branches occurs. Consider a system on the erosion-limited stable equilibrium branch (e.g. at point a in the figure). Next, consider some change to, for example, the bed shear stress such thatẼ increases beyondẼ crit . The system will then evolve over time to a supply-limited equilibrium, e.g. at point b in the figure. If there is a sufficient amount of sediment available, this leads to a catastrophic increase in concentration. If the change toẼ is completely reversed, the potential erosion is still larger than the deposition. Therefore, all sediment remains suspended in the water column and the system remains on the supply-limited branch, with even higher concentrations near the bed (point c in the figure). A transition back to the erosion-limited equilibrium branch will only happen ifẼ is further reduced to below the point where the supply-limited branch starts. The erosion and supply-limited branches thus describe hysteresis between the transition from low to high concentrations and back.

Extension to tidal flows
The instantaneous bed shear stress in tidal flows varies from its maximum at peak ebb or flood to virtually zero at slack tide. As a result, the near-bed sediment concentration varies over the tidal cycle. To demonstrate this, we force the water column model by a tidally varying water level gradient consisting only of an M 2 tidal constituent: ζ x = |ζ x | cos(ωt), where |ζ x | denotes the amplitude of the forcing and ω is the angular frequency of the M 2 tide. Figure 4 shows model results for |ζ x | = 10 −5 (left) and 2.1 · 10 −5 (right), an infinite sediment supply and all other variables having their default values (Table 1). The figure shows the temporal evolution of the near-bed volumetric concentration as a function of the dimensionless erosion parameter during half a tidal cycle from slack tide to peak tide and back to slack tide. The solid and dashed blue lines in the figure indicate the stationary equilibrium and are identical to the solution in Fig. 1 obtained for stationary forcing conditions. In Fig. 4a, the tidal forcing is small, leading to a small concentration and small amplitude of the dimensionless erosion parameter. The tidal signal is always close to the solid line that indicates the stationary equilibrium. Therefore, it is concluded that the change of the concentration during the tidal period is small enough to adjust to the stationary equilibrium at all time. Figure 4b shows the temporal variation of the near-bed concentration for a larger tidal forcing. The concentration now deviates significantly from the stationary equilibrium. During accelerating tide, the concentration increases, but remains below the stationary equilibrium concentration. Conversely, during most of the decelerating tide, the concentration decreases, but remains above the equilibrium concentration. At peak tide, the maximum value ofẼ is attained andẼ exceedsẼ crit , allowing for an unbounded increase of the concentration. After peak tide,Ẽ still exceedsẼ crit for some time resulting in a further increase of the near-bed concentration, even though the bed shear stress already decreases. The increase of the concentration only stops when the pair (φ bed ,Ẽ) crosses the unstable stationary equilibrium (dashed line), well after peak tide. Within the remainder of the time until slack tide, the concentration reduces to almost zero. Therefore, there is no net change of the concentration over the tidal cycle. Figure 4c, d visualises the same experiments in a different way, by plotting the effective settling velocity versus t/T , where T is one tidal period. For the weak tidal forcing, the variation in settling velocity is weak and in phase with the tidal velocity near the bed. For the stronger tidal forcing, the variation in settling velocity is stronger and we again see that the lowest settling velocity (highest concentration) is attained after peak tidal velocity.
Even though the concentration varies over the tidal cycle, we are mainly interested in the average behaviour over a tidal cycle. Therefore, we will further consider the tidally averaged concentration. It is not possible to find a relation between the tidally averaged near-bed concentration φ bed and the tidally averaged dimensionless erosion parameter Ẽ only on the basis of near-bed processes, as was the case in stationary flows (see Appendix for details). Therefore, the water column model has to be solved for different tidal forcing conditions. In these experiments, we force the water motion by a prescribed tidally varying water level gradient, an infinite stock and use the default parameter values (Table 1). Each Fig. 4 Top: instantaneous values of φ bed andẼ for half a period of an M 2 tidal wave for two forcing strengths. This is compared to the stationary equilibrium branch (blue line). At slack tide,Ẽ = 0, whileẼ is at its maximum for peak flood and ebb velocities. Bottom: the same simulations but now visualising the effective settling velocity versus dimensionless time over one M 2 tidal cycle experiment for one value of |ζ x | is allowed a sufficiently long spin-up time to attain periodic conditions, resulting in a single equilibrium value of Ẽ and φ bed . The branch of stable equilibrium states is constructed by repeating this procedure for various water level gradient amplitudes. The resulting branch is plotted as the green solid line in Fig. 5a. For small values of Ẽ this branch corresponds closely to the stationary equilibrium branch (solid blue line). This is because the concentrations are small and are almost equal to the stationary equilibrium concentration throughout every instance of the tidal cycle (cf. Fig. 4a). For increasing Ẽ , the tidal and stationary branches start to differ. For these settings it is found that Ẽ crit is smaller in the tidal case than in the stationary case. To see why this is possible, consider some value of Ẽ . The instantaneous value ofẼ is smaller than Ẽ during approximately half of the tidal cycle and larger during the remainder of the tidal cycle. If Ẽ is sufficiently large, the instantaneousẼ also exceeds the instantaneousẼ crit for some smaller time interval of the tidal cycle. If, during this time interval, the concentration increases so much that the concentration cannot decrease back during the remainder of the tidal cycle, the average concentration keeps increasing over time and no equilibrium condition can be found. This means we are in domain II (cf. Fig. 1 for stationary flows) and we must have passed the critical point. It depends on the model settings how longẼ has to exceedẼ crit before this happens. We will elaborate on this below.
The unstable equilibrium branch is, by definition, never obtained using time-integration models. Nevertheless, this branch can be inferred from the model by introducing a limited sediment availability. As demonstrated in Fig. 2 for stationary flows, a limited sediment availability results in a new branch of stable supply-limited equilibrium solutions. This branch exists at all values ofẼ on the right side of the equilibrium curve. By running the model for various values of the sediment availability and water level gradient amplitude, we find various supply-limited branches. The left-most endpoints of these supply-limited branches mark the location where they cross the unstable equilibrium. This way, the unstable equilibrium solution can be reconstructed. The resulting line of unstable equilibria is plotted as the dashed green line in Fig. 5a. This is constructed from the result of over 4000 model simulations with the M 2 water level gradient amplitude |ζ | x varying between 0 and 5·10 −5 and the sediment stock varying between 5 and 400 kg/m 2 . The result for a sediment stock of 15 kg/m 2 is plotted as the orange supply-limited equilibrium line. Figure 5b shows the erosion-limited equilibrium branches and a supply-limited branches (S = 100 kg/m 2 ) for a different set of parameter values. We have decreased the erosion parameter M to 2 · 10 −4 g m/(L s) and increased the range of |ζ x | between 0 and 5 · 10 −4 . Compared to the default value of M, a larger |ζ x | and thus a higher degree of vertical mixing is required to suspend the same amount of sediment. This leads to a more uniform distribution of the sediment over the water column. Therefore, a much higher stock is required to find a supply-limited branch at similar values of φ bed compared to the results shown in Fig. 5a. The value of Ẽ crit in Fig. 5b is almost the same as the instantaneous critical value, and significantly larger than the critical value in Fig. 5a. As this concerns the tidal average in a symmetric tide, this implies that the instantaneousẼ exceeds the stationary critical value for almost half the tidal cycle. However, in this time, the concentration does not increase so much that it cannot return to low values around slack tide. In other words, the timescale for the concentration to change is longer. This can be explained from a rough scaling analysis of the depthintegrated version of Eq. 4. After applying the boundary conditions, this equation reads LetC denote a typical depth-averaged concentration and T denote a time-scale at which the average concentration changes, then HC T ∼ −w s c bed +Ê, i.e.
The depth-averaged concentration scales with the near-bed concentration and with some increasing function of the eddy diffusivity. In other words, an increase in the eddy diffusivity at the same near-bed concentration, leads to an increase in the depth-averaged concentration. Therefore, at the same near-bed concentration and potential erosion, a higher degree of mixing of the water column leads to a higher value ofC and thus a larger time-scale for the concentration to change. For an even larger tidal forcing and smaller erosion parameter, Ẽ crit can exceedẼ crit . The steepness of the tidal equilibrium curve in Fig. 5b is fairly similar to that of the stationary curve, indicating that sediment-induced turbulence damping is not important. This is because the tidal forcing is so strong that only very high concentrations cause stratification. However, at such concentrations, the fall velocity decreases due to hindered settling, therefore leading to a more uniform distribution of sediment over the water column and reducing stratification. This phenomenon has for example been observed in the Yellow River (Van Maren et al. 2009).
Though the horizontal axis of Fig. 5 shows the dimensionless erosion parameter, the model results where obtained by varying the water level amplitude. A change in the water level gradient amplitude not only affectsẼ through the bed shear stress, but also affects the rate of mixing of the water column. A variation of another parameter, such as the erosion parameter only affectsẼ and not the rate of mixing. Figure 6 again shows the stable and unstable equilibrium solution branches, but now obtained by varying the erosion parameter M at a fixed value of |ζ x | of 10 −4 . The orange supply-limited branch corresponds to S = 100 kg/m 2 . Similar to Fig. 5b, Ẽ crit is close to the stationary critical value, because the tidal forcing is relatively large and the time-scale for the concentration to change is relatively long. The main difference with the figures found before is the shape of the supply-limited branches. These are more flat or even slightly upward sloping. This is because a change inẼ now does not coincide with enhanced or decreased mixing of the water column; the vertical distribution of sediment remains the same at different values ofẼ.

Role of hindered settling
In a tidal flow, the instantaneous sediment concentration changes constantly in the direction of the stationary equilibrium or keeps increasing if a stationary equilibrium does not exist. The essential mechanism explaining the bifurcation in the stationary case therefore also explains the bifurcation in the tidal case. This mechanism is hindered settling. Additional to hindered settling, inertia is important in tidal flows: the instantaneousẼ needs to exceedẼ crit for a sufficiently long time for the concentration to move away from the erosion-limited equilibrium.
In order to further support the conclusion that hindered settling is essential for the bifurcation, we again consider the depth-integrated sediment balance (12) without hindered settling. We assume a dynamic tidal equilibrium exists and take the tidal average. The result reads i.e. the tidally averaged near-bed concentration depends linearly on the dimensionless erosion parameter and does not depend on the rate of turbulent mixing of the water column. For comparison, the solution with hindered settling leads to the non-linear relation (1 − φ bed ) m φ bed = Ẽ (see Eq. 20 in Appendix). The solution without hindered settling is plotted in Fig. 7 (blue line). Clearly, there is no bifurcation. High near-bed concentrations of the order of the gelling point can nevertheless be attained, but only at values of Ẽ which are one order of magnitude larger than in the case with hindered settling (compareẼ crit = 0.067). The figure also shows the depth-and time-averaged concentration without hindered settling, with and without sediment-induced turbulence damping (red solid line and red dotted line respectively). In these experiments, all parameters have their default values (Table 1), the stock is infinite and the variation of Ẽ is obtained by varying Fig. 7 Relation between φ bed and the tidally averaged dimensionless erosion parameter if hindered settling is removed from the model. The blue line shows the near-bed concentration. The red line shows the depth-averaged concentration, compared to the depth-averaged concentration if both hindered settling and sediment-induced damping are omitted (red dotted line). The figure is obtained by varying the tidal M 2 amplitude of ζ x in the model |ζ x |. The figure shows that sediment-induced turbulence damping has a significant effect by reducing the depthaveraged concentration compared to the situation without turbulence damping. It also shows that turbulence damping alone does not lead to any bifurcations.

Discussion
In this section, we investigate the sensitivity of the model results to different parametrisations and parameter choices. In Section 5.1, we discuss several formulations and parameter choices related to hindered settling. In Section 5.2, we look at several formulations for erosion. Finally, the role of flocculation is discussed in Section 5.3. An exhaustive discussion of each available formulation is beyond the scope of this paper, but the influence of different parametrisations will be at least qualitatively discussed and the methodology to perform an in-depth analysis will be indicated where deemed necessary.

Other parametrisations of hindered settling
The existence of a limit point when the erosion parameter attains a critical value, depends crucially on hindered settling. Three parametrisations of hindered settling prevail in literature. The first is the formulation proposed by Richardson and Zaki (1954) (see Eq. 10), in which the parameter m is empirical and still needs to be chosen. Laboratory measurements indicate m is likely between 2 and 5. The second formulation is by Dankers and Winterwerp (2007). It reads Here, φ p is the volumetric concentration of primary particles, defined as c/ρ s , with ρ s the dry density of the soil. Dankers and Winterwerp (2007) propose m = 2 after a comparison to measurements. The third formulation is by Toorman and Berlamont (1993), in which we neglect the contribution related to consolidation and rewrite to our notation as where α is an empirical parameter. Mehta (2014) summarises data from three studies, with α ranging between 5.5 and 8. Here we use α = 7, but the results are qualitatively the same for all values in this range. Figure 8 shows the stable and unstable erosion-limited branches in a stationary flow for various values of m in Eq. 10, m = 2 in Eq. 14 and α = 7 in Eq. 15. All curves show the same characteristics as found before: the existence of a bifurcation point, where one stable and one unstable solution branch meet. The various curves only differ in the value of the critical dimensionless erosion parameter and the corresponding near-bed concentration.
As was explained in Section 3, the necessary condition for the existence of a bifurcating solution is the existence of a maximum of the settling flux. This is obtained for all hindered settling parametrisations studied here. For the Richardson and Zaki (1954) formulation this bifurcation point can be easily derived to occur at φ = 1 m+1 . The value ofẼ crit can also be expressed analytically in terms of m and is given by m m+1 m 1 m+1 .

Other descriptions of erosion
Although much experimental work has been done on erosion, there is no consensus on the best mathematical description of erosion (e.g Sanford and Maa 2001;Mehta 2014). One popular formulation for erosion can be generalised aŝ wherez indicates the depth of erosion relative to the bottom of the water column (z = −H ). Many variations on this model exist in which the time-and/or depth-dependence of the erosion parameter M or critical shear stress are omitted. The time-and depth-dependence of τ c is often associated with consolidation, but may also be related to biological influence (e.g. Le Hir et al. 2007). The erosion parameter M is often kept constant over time and depth, but may be allowed to vary related to sediment composition or consolidation (e.g. Sanford and Maa 2001). Often n = 1, M(z, t) = M and τ c (z, t) = τ c are chosen for practical reasons. In this case, this formulation is typically associated with well-consolidated soils (referred to as type II erosion). Another formulation for erosion readŝ This form is sometimes favoured for unconsolidated soils (referred to as type I erosion), where the strength varies with depth.
The results presented in this paper do not hold for all forms of the general erosion formulations (16) and (17) with arbitrary t or z dependence. It depends on the exact formulation if the same qualitative results hold. We will not go into the effect of choosing any particular form of these general erosion formulations, as there is not one standard form other than Partheniades' formulation, which is used throughout this paper. The purpose of this section is to show that the same methods applied in this paper can be applied to study many forms of the erosion formulation. This mainly requires a suitable redefinition ofẼ. We will illustrate how one may determineẼ and the existence of a bifurcation point in this case by a simple example.
Consider a stationary flow forced by a constant water level slope. Let the critical shear strength increase with depth, e.g. due to consolidation, and let the erosion be described by Eq. 16 with n = 1. The condition for the equilibrium near-bed concentration of Eq. 11 is modified to Reordering this expression yields where all terms depending on φ are on the left-hand side, including τ c , which depends on φ throughz. Within this formulation, the right-hand side should be considered as the dimensionless erosion parameter, i.e.Ẽ = M 0 w s,0 c gel τ . Equation 18 has a bifurcation point if the left-hand side of this expression has a maximum. Whether this is the case depends on the relation between the critical shear stress and φ bed . If a maximum to the left-hand side of Eq. 18 exists, the results of this paper also hold qualitatively for this erosion formulation. If there is no such maximum, the stable erosion-limited solution branch exists for all values ofẼ. High concentrations can then be attained for erosion-limited conditions and no hysteresis is observed in the transition between low-and high-concentration states.

Role of flocculation
Fine cohesive sediments are known to flocculate. The size and fall velocity of sediment flocs may differ from that of the primary sediment particles by several orders of magnitude. It is known that the degree of flocculation depends on the concentration of primary particles and shear stress in the water column. Therefore, the flocculation also varies with φ bed andẼ and is expected to modify the equilibrium curves in φ bed -Ẽ space.
Experiments on flocculation show that the fall velocity of the sediment increases due to flocculation processes, but only up to values of φ around 0.1 (see e.g Figure 5.3 of Winterwerp and Van Kesteren 2004, and references therein). For higher volumetric concentrations, hindered settling takes over and dominates the settling velocity. Therefore, it is expected that flocculation only modifies the shape of the erosion-limited equilibrium branch in the φ bed -Ẽ diagram. Taking w s,0 as the settling velocity of flocs atẼ = E crit , the unstable branch andẼ crit are hardly affected by flocculation. Flocculation is therefore not expected to have any important implications for the results found in this paper.

Necessity of hindered settling to reproduce high concentrations
The theory developed in this work helps to explain why some estuaries and coastal systems can support very high suspended sediment concentrations, while the bed shear stress is not particularly high. Here, we demonstrate this using two examples. The first is of the Amazon shelf, with concentrations of more than 100 g/l (Kineke and Sternberg 1995), while the average bed shear stresses are estimated by model studies to be of the order of 1 Pa (Gabioux et al. 2005). The second example is the Ems River, where concentrations of over 30 g/l have been observed (Talke et al. 2009), while bed shear stresses are up to 2 Pa (Van Maren et al. 2015). We compute the order of magnitude of the tide-averaged dimensionless erosion parameter for these systems, using its definition (11). The erosion parameter M is a calibration parameter in most models. The critical bed shear stress τ c can be measured, but has a significant measurement uncertainty and natural variability. However, both parameters have a well-established range of typical values. We estimate M to be of the order of 10 −4 − 10 −3 g/l m/s and τ c in the range 10 −1 −10 1 Pa (see, e.g. Partheniades 2010). We assume a typical order of magnitude of the settling velocity of 1 mm/s and c gel of 100 g/l. For the Amazon shelf and the Ems River, this means that typical values of Ẽ are 10 −3 to 10 −1 .
If the effects of hindered settling were not included in the model, the mean value ofẼ equals the mean nearbed volumetric concentration (also see Eq. 13). This means that φ can at maximum attain concentrations of around 10 −1 , corresponding to 10 g/l with the current choice of the gelling concentration. This is much lower than the concentrations observed on the Amazon shelf and also lower than the concentrations observed in the Ems River. Without hindered settling the observed concentrations can therefore not be reproduced using realistic parameter settings. In a model that includes the effects of hindered settling, the only requirement for obtaining high concentrations is that Ẽ exceeds Ẽ crit , which is of the order of magnitude of E crit = 0.067. This value is within the range of typical values Apart from requiring a sufficiently high Ẽ , high concentrations also require a large supply of sediment. In estuaries like the Ems, where the sediment mainly enters at the mouth of the estuary, this requires a strong trapping of sediment. This trapping depends on the along-channel dynamics of the system and has not been investigated in this study. The role of hindered settling on the trapping behaviour of systems like the Ems thus still needs to be investigated.

Conclusions
Over the last several decades, some estuaries have undergone a transition from low or average turbidity levels to a so called hyperturbid state. In order to improve our understanding of this transition, we have systematically identified the essential processes that allow for high suspended sediment concentrations in an estuarine water column. We adopt a framework where we characterise the dynamic equilibrium state of the water column as either erosion or supply limited. In the erosion-limited state, the concentration is bounded by the ability of the flow to erode and keep sediment in suspension. In the supply-limited state, all available sediment is suspended in the water column, so that the maximum concentration is bounded by the amount of sediment available in the system.
The behaviour of the water column in these equilibrium states can be described in terms of a time-averaged dimensionless erosion parameter Ẽ , which parametrises the rate of erosion. If Ẽ is above some threshold Ẽ crit , only the supply-limited equilibrium exists and the actual concentration entirely depends on the amount of sediment in the system. This potentially allows for extremely high concentrations. For values of Ẽ below the threshold, both types of equilibria may exist. Provided the sediment supply is sufficiently large, this implies that the flow can be in one of two states: an erosion-limited state with relatively low concentrations and a supply-limited solution with relatively high concentrations. The actual state of the system depends on its history. If a system is in the low-concentration, erosion-limited state, it can only jump to the supply-limited state if Ẽ increases beyond the threshold value. This might result in a dramatic increase in concentration. The jump back from the supply limited to the erosion-limited state only happens for Ẽ significantly smaller than the threshold, so that there is hysteresis in the transition between the two equilibrium states.
The threshold value of Ẽ and the hysteresis in the transition between states exist because of the effects of hindered settling. These phenomena therefore disappear if hindered settling is not taken into account, strongly altering the behaviour of the model. Hindered settling allows for a positive feedback loop that might lead to a catastrophic increase in concentration if Ẽ exceeds the threshold. An investigation into the different formulations and parameter values parametrising hindered settling shows that this behaviour is a robust characteristic of all known formulations for hindered settling. Sedimentinduced damping of turbulence is often thought to be an essential mechanism promoting a sudden transition between states of low and high sediment concentration, but is found to play no essential role in the behaviour identified here. However, it is important for the vertical distribution of sediment over the water column and quantitative determination of the threshold value for Ẽ in tidal flows.
The high suspended sediment concentrations observed in some estuaries are only possible in a supply-limited state. Therefore, it can be concluded that hyperturbid systems should at least satisfy the criterion Ẽ > Ẽ crit or have satisfied this criterion at some point in the past. This observation is an important step towards the understanding of why some systems have become hyperturbid and the assessment of the risk that a particular non-hyperturbid system may become hyperturbid in the future.