Can the Scheldt River Estuary become hyperturbid?

We investigate the hypothesis by Winterwerp and Wang (Ocean Dyn 63:1279–1292, 2013) that channel deepening in the Scheldt River Estuary could lead to a large increase in suspended sediment concentrations, with subsequent severe consequences to primary production and navigation. To this end, we use an idealised model to investigate the long-term development of the sediment concentration under the uncertainty of future changes in model parameter values and channel deepening. The water motion is calibrated to recent conditions after which the sediment concentration is validated against long-term observations and is subsequently tested for a wide range of parameter settings and deepening scenarios. We also investigate the effect of anthropogenic dumping of dredged sediments in the estuary on the sediment concentration. Deepening the channel, but keeping all other model parameters the same, we find lower long-term average sediment concentrations in most of the estuary. Thereby, our results suggest that deepening in the Scheldt alone cannot lead to high sediment concentrations, and we suggest to reject the investigated hypothesis. Further study of uncertain model parameters reveals that an increase of the erosion parameter by an order of magnitude allows for the development of high concentrations of several tens of grams per liter near the bed in narrow turbidity zones. It is unknown whether such an increase of the erosion parameter can happen in the future, which stresses the importance of further research into the factors that can lead to a change of this parameter.


Introduction
Suspended sediment dynamics is an important subject in the management of the Scheldt River Estuary, an estuary located in Belgium and Netherlands, see Fig. 1. Its importance is related to two factors. Firstly, primary production in the Scheldt estuary is to a large extent light limited due to suspended sediments (Kromkamp and Peene 1995). Secondly, dredging of fine sediments poses a significant cost in the maintenance of the navigation channel to the port of Antwerp (IMDC et al. 2013a). These issues related to fine sediments have become increasingly important, as the long-term suspended sediment concentration has gradually increased by several tens of milligrams per liter over the last three decades in the Lower Sea Scheldt (km 55-95) and locally in the dry season in the Upper Sea Scheldt (> km 95) (Vandenbruwaene et al. 2016;Maris et al. 2017). At the same time, many anthropogenic changes to the estuary have been made, including extensive deepening of the channel for the mining of sand and improvement of navigability, the construction of locks and harbour basins and the development of intertidal area (see Van Braeckel et al. (2006) and Jeuken et al. (2007) for an overview). Moreover, sewage treatment has improved, thereby affecting the organic content of sediments and hence sediment properties (Maris and Meire 2017).
It has been suggested by  that deepening in the long term may lead to the development of hyperturbid conditions in the Scheldt. According to their hypothesis, deepening causes a deformation of the tide that leads to an increasing import of fine sediment. The imported sediment leads to a reduction of the hydraulic drag, which supposedly leads to a further deformation of the Fig. 1 Map of the tidal part of the Scheldt River from Vlissingen to Ghent. The Scheldt River is separated into the Dutch Western Scheldt (< km 55) and the Belgian Sea Scheldt (> km 55). At the head of the estuary, the estuary is fed by water from the Upper Scheldt -Leie catchment. Two more tributaries are marked on the map in italics: the Dender and Rupel Rivers. Major locations where fine sediments are dredged and dumped in the Sea Scheldt are marked in green and red, respectively. Dredging and dumping also happens in the Western Scheldt but in smaller amounts, which are not considered in this study tide and more sediment import, hence leading to a dramatic increase of the sediment concentration. This hypothesis is mainly based on examples of the Ems (Germany) and Loire (France) Rivers, which have become hyperturbid following decades of substantial channel deepening. Despite regular and long-term monitoring of sediment concentrations in the Scheldt River since the 1990s, it remains unclear whether the observed long-term trends are part of a development towards hyperturbidity (Vandenbruwaene et al. 2016;Maris et al. 2017). Statistical analysis of trends in the observations is not only complicated by a high degree of natural variability but also by memory effects, by which the sediment concentration depends on both the hydrodynamic conditions of several months in the past (Brouwer et al. 2018) and recent sediment dredging and dumping. Sediments dredged in the port of Antwerp are dumped back into the channel of the estuary a few kilometres upstream at rates that far exceed fluvial sediment supply. Hence, variations in the observed sediment concentrations are strongly influenced by the spatial and temporal variability in anthropogenic sediment dredging and dumping (e.g. Depreiter et al. 2015). While such memory effects and the effects of dredging and dumping are included in models, only a few model projections have been made of the response of sediment concentrations to channel deepening (Van Kessel et al. 2008) and dumping strategy (Van Kessel and Vanlede 2010;IMDC et al. 2013b). Some of the main challenges for such modelling studies are the long timescales at which the sediment concentration varies and large uncertainty in model parameters.
The goal of this study is to investigate if the Scheldt River can become hyperturbid as a response to channel deepening given the uncertainty in model parameters values. Moreover, we aim to gain insight into the processes and parameters that are most important to the sediment dynamics in the Scheldt. It should be stressed that we do not aim to explicitly describe variability on weekly, seasonal or yearly timescales in the past or future. Rather, we want to qualitatively capture the most important underlying physics, which allows us to extrapolate the modelled trends to uncertain future scenarios and the corresponding long-term average behaviour. To this end, we use the iFlow model, which is a width-averaged, idealised process-based model (Dijkstra et al. 2017). The model is used to directly compute the long-term equilibrium water motion and sediment concentration given prescribed geometry and forcing conditions, thus quickly showing the long-term response of the estuary to changing depth and model parameters. As the model is fast, it allows for extensive study of parameter sensitivity.
The set-up of the iFlow model used in this study is based on the model used by Dijkstra et al. (2019) to simulate the transition to hyperturbid conditions in the Ems after channel deepening. In addition, several processes thought to be essential for the sediment dynamics in the Scheldt River are added (Section 2). The model is calibrated against water levels for conditions of the year 2010 and modelled sediment concentrations are presented and compared to the long-term averaged observations (Section 3). Furthermore, the sensitivity of the results to the sediment settling velocity and erosion parameter is systematically analysed. Next, in Section 4, the model is applied to configurations with smaller and larger depths, investigating the response to past and possible future large-scale deepening. The physical processes explaining the results are analysed in Section 5. These processes are discussed in the context of the processes that act in the Ems estuary and in context of the model uncertainty in Section 6. This chapter ends with a summary of the main conclusions in Section 7.

The iFlow model
The model used in this study is the iFlow model (Dijkstra et al. 2017), a width-averaged idealised process-based model for water motion and sediment dynamics in estuaries and tidal rivers. The model solves approximations of the width-averaged continuity, momentum and sediment balance equations, and hence it resolves the flow and sediment concentration in the along-channel direction x, vertical direction z and in time t. The model resolves a dynamic equilibrium state, which means that the water motion and sediment concentration may vary on the tidal timescale, but not on the subtidal timescale. The main model features and simplifications are discussed below.
The geometry of the model is represented by smooth fits of the measured width and depth, thereby ignoring the effect of small-scale bathymetric features on the estuary-scale dynamics. Additionally, it is assumed that the water surface elevation is small compared to the subtidal depth. These assumptions allow for the use of scaling and perturbation methods, which lead to systems of mathematical equations that can be solved at low computational costs and allow for making a decomposition of the water motion and sediment transport into contributions by individual physical forcing mechanisms. For a further detailed description of the basic model equations and solution methods, we refer to Dijkstra et al. (2017).
The effects of turbulence in the model are parametrised by using an eddy viscosity, eddy diffusivity and a quasiquadratic bed friction formulation (see Appendix A or Dijkstra et al. (2019) for a mathematical description). The eddy viscosity and diffusivity are assumed to be depth uniform and constant in time but depend on the tideaveraged, depth-averaged flow velocity and depth. The quasi-quadratic bed friction formulation is obtained by linearising the bed friction but allowing the linearised bed friction coefficient to depend on the tide-averaged, depthaveraged velocity. The model furthermore accounts for the effects of sediment-induced damping of turbulence, which represents the reduction in turbulent mixing due to vertical density stratification by sediment. This is accounted for by reducing the eddy viscosity and eddy diffusivity using tide-averaged, depth-averaged damping functions based on the formulations of Munk and Anderson (1948). The bed friction is reduced by using a tide-averaged damping function of the Richardson number, following e.g. Adams and Weatherly (1981). The turbulence model is calibrated using two coefficients: a dimensionless roughness height z * 0 that affects the eddy viscosity, eddy diffusivity and bed friction, and a background shear u z,min that controls the strength of sediment-induced turbulence damping.
Sediment is modelled assuming a single class of non-cohesive sediment. The settling velocity is assumed constant but may be reduced due to effects of hindered settling, which represents a decreased sediment settling velocity due to particle-particle interactions when sediment concentrations are high. At the bed, the subtidal stock of sediment available for erosion is accounted for by the model. Erosion or resuspension of sediment from the bed occurs through an erosion law inspired by the Partheniades formulation (e.g. Kandiah 1974) but without a critical shear stress: Here, E is the erosion rate, M is an erosion parameter (in s/m), τ is the bed shear stress and f is the erodibility, which is a measure for the amount of available sediment on the bed and takes a value between 0 and 1. In regions with f < 1, the low availability of sediment limits the erosion (supply limited). In regions where f = 1 is predicted, excess sediment is available in the bed, and erosion is limited instead by bottom stress (erosion limited). This will locally result in continual net deposition of sediment at the bed in equilibrium (see Section 5.1 for additional explanation of erodibility). Salinity is modelled diagnostically, by setting a prescribed salinity gradient that is fitted to observations. The salinity is assumed to be depth uniform and constant over the tidal cycle, approximating conditions in very well-mixed estuaries.
The model is forced by constant M 2 and M 4 tidal amplitudes and a subtidal depth-averaged sediment concentration at the mouth. Furthermore, constant river discharges are prescribed at the head of the estuary and at the confluences with the Dender and Rupel tributaries (see Fig. 1), and a discharge-dependent fluvial supply of sediments is prescribed at the discharge locations. The inflow of water and sediment from tributaries is a new feature in iFlow, and its mathematical implementation is described in Appendix B.
Together, the equations form a non-linear set that is solved iteratively using the procedure summarised in Fig. 2. The solution procedure consists of a combination of algebraic relations and numerical and analytical solution methods, where analytical methods are used for speed and accuracy whenever possible. The numerical computation of the water level and erodibility are done using a second-order finite difference method on an equidistant grid with 250 grid cells. All vertical profiles are computed analytically and do not require a numerical grid. In the time dimension, the model is solved in terms of a constant subtidal component and M 2 and M 4 tidal constituents (i.e. spectral method). Hence, time stepping routines are not necessary and there is no time step or spin-up time.

Description of the Scheldt River Estuary and model forcing
The Scheldt River Estuary is modelled as a single channel from the mouth at Vlissingen to the tidal weir and locks at Ghent, 161-km upstream. Tidal propagation into the tributaries (see Fig. 1) is not explicitly taken into account. The width of the Scheldt in the model is represented by the average of the width at the surface at high and low water and is fitted by a smooth function. The width-averaged depth of the Scheldt is derived by dividing the cross-sectional area by the width at high and low water, subsequently subtracting the average water level elevation at high and low water and then taking the average. The resulting depth is fitted using a smooth polynomial function, see Fig. 3. The procedure for deriving the depth is slightly different to that used by Brouwer et al. (2018) for the Scheldt and is used because, in contrast to the procedure of Brouwer et al. (2018), it can be repeated for the historical high and low water data used in Section 4. The model is forced by an M 2 and M 4 tide at the mouth, representing year-averaged tidal conditions, obtained using a complex demodulation analysis (e.g. Jalón-Rojas et al. 2016) on the 10-min resolution tidal elevation observations at Vlissingen for 2009. This yields an M 2 amplitude of 1.81 m and an M 4 amplitude of 0.16 m with a relative phase difference between the M 2 and M 4 tide of −4 degrees. Fresh water discharges into the Scheldt at three locations: at the upstream boundary from the Upper Scheldt -Leie system, at Dendermonde (km 123) from the Dender tributary and at Rupelmonde (km 95) from the Rupel tributary. The Rupel has several tributaries itself and its fresh water discharge equals the sum of the discharges of its tributaries. We use the average discharge for a year, summer (Jul-Sep) and winter (Jan-Mar) averaged over the years 1971-2017 (data from www.waterinfo.be), see Table 1.
Several methods can be used to estimate the fluvial sediment supply as a function of the river discharge. Vanlierde et al. (2016) present a simple regression model relating the instantaneous sediment concentration from 7hourly measurements to the instantaneous river discharge for data of 2015. However, the number of data points and quality of the fit are low, so that this method is not reliable for estimating the long-term sediment supply. Therefore, we choose to correlate the year-averaged estimated fluvial sediment supply and year-averaged river discharge of each tributary. The fluvial sediment supply is estimated using sediment concentration data obtained from weekly samples near the river bank between 2001 and 2015 (data from Vanlierde et al. (2014Vanlierde et al. ( , 2016 and Vandenbruwaene et al. 2017), see Fig. 4. A linear fit is chosen over more conventional power relations to prevent over-fitting of the small amount of data. The obtained relations for the fluvial sediment load (in kg/s) read A comparison of 7-hourly and weekly sediment measurements by Plancke et al. (2017) indicates that fits of Eqs. 2 and 4 likely underestimate the sediment supply per tributary by a factor 3 to 6 (3.5 for the entire estuary for 2016). However, as the 7-hourly data are only available for a few years, we base our fits on the weekly data. A sensitivity study (not presented here) showed that increasing the fluvial supply by a factor 3 to 6 does not notably change the results presented in this study. At the seaward boundary, a depth-and tide-averaged sediment concentration c sea = 0.06 kg/m 3 is prescribed, based on observations. The sediment settling velocity w s,0 (a) Bed level (b) Width

Dredging and dumping
Dredging of fine sediments predominantly takes place at eight sills, lock entrances and harbour basins located  Vandenbruwaene et al. (2017) between km 60 and 71, see the green dots in Fig. 1 (IMDC et al. 2013a). All of the dredged fine sediments are dumped back into the navigation channel a few kilometres upstream at Punt van Melsele (km 73) and Plaat van Boomke / Oosterweel (km 78), see the red dots in Fig. 1. Exceptions are the years 1990-2000, when 300,000 tons of dry fine sediment was removed from the estuary. Figure 5 shows the amount of dredging and dumping in the Western Scheldt (< km 55) and Sea Scheldt (> km 55) in tons of dry sediment per year, compared to the estimated fluvial fine sediment supply from all tributaries (data from IMDC et al. (2013a), Vandenbruwaene et al. (2016) and Barneveld et al. 2018). This shows that the sediment source due to dumping exceeds the fluvial supply by an order of magnitude.
As harbour basins, sills and shallow areas are not explicitly resolved by the model, sediment deposition is not fully resolved. Hence, if dredging were taken into account by imposing a sink in the model, sediment would be extracted while an insufficient amount is deposited according to the model, leading to negative sediment concentrations. Therefore we do not take dredging into account in the model, but we do consider dumping. Dumping of sediment is represented using continuous point sources at km 73 and 78 with rates of 60.5 and 98.5 kg/s, respectively, corresponding to the average dumping rate in the Sea Scheldt between 2001 and 2015. Dumping of sediment in the Western Scheldt is neglected as the dumping volumes are relatively small, especially when considering the much larger volume of the Western Scheldt.

Model calibration
The M 2 water level is calibrated to observations by varying the roughness parameter until a best fit to data is found; the optimal dimensionless roughness value is z * 0 = 5·10 −4 . The resulting M 2 and M 4 water level amplitude and phase and cross-sectionally averaged velocity amplitude are plotted in Fig. 6. The model results are compared to water level observations from 2009 and velocity observations from one 13-h cross-section measurement in 2009 (data from Clear-water settling velocity 2 mm/s c gel Gelling concentration 100 kg/m 3 Turbulence Velocity gradient for background turbulence production 0.03 1/s Rijkswaterstaat, HIC and Flanders Hydraulics Research). The data are analysed using complex demodulation after which the average amplitude and phase over the year is used. The overall observed patterns for the M 2 tidal amplitude are reproduced, but the M 4 tidal amplitude is overestimated locally by more than a factor two. The tidal phases of both the M 2 and M 4 tide are reproduced well.
The modelled M 2 tidal velocity shows two maxima near km 20 and km 120 and minima at the mouth, km 70 and at the landward boundary, where the tide vanishes. The same pattern is observed in the measurements and the overall magnitude of the modelled velocity corresponds to the measurements. The modelled M 4 velocity increases up to km 140 before it vanishes at the landward boundary. While a maximum in the M 4 velocity in the upstream part of the estuary is also found in the measurements, measured velocities are much lower than the modelled velocities, similar to what was found for the tidal elevation.

Sediment concentration compared to data
In order to verify the performance of the model, we compare the modelled sediment concentration in the 2010 case using default parameter settings and the year-averaged discharge ( Table 2) to long-term averages of sediment concentration observations. Observations of surface sediment concentrations were collected between 1990 and 2015 within the MWTL (Western Scheldt) and OMES (Sea Scheldt) programmes. The data was gathered bi-weekly to monthly, independent of the tidal conditions, by taking water samples (see Maris and Meire (2017) for details on the OMES programme). Observations of the depth-averaged sediment concentration were collected in the period 2001-2015 as part of the OMES programme and are based on pump samples at different depths with approximately equal coverage of the entire water column (Vandenbruwaene et al. 2016). Finally, we have included data from four permanent optical measurement stations at a depth roughly halfway the average water depth, hence roughly representing depth-averaged concentrations. Variations on timescales smaller than one M 2 tidal cycle have been filtered from the observations (see the caption of Fig. 7 for more information per station).
Modelled subtidal sediment concentrations for a yearaveraged discharge (Q = 110 m 3 /s) are plotted in Fig. 7. First, we compare the observations and model results When considering the model with dumping of sediment ( Fig. 7b), the main difference is in the depth-averaged concentration. The model results now display an ETM at km 85, approximately corresponding to the observations on both location and order of magnitude of the sediment concentration. The location of the ETM does not correspond exactly to the dumping location, indicating that sediment concentrations are not just elevated because of the dumping of sediment but because sediment is trapped by the flow some distance upstream from the dumping location.
As the highest sediment concentrations typically occur during the relatively dry summer months, we also verify the model results using the average summer river discharge (Q = 56 m 3 /s) to measurements taken in the summer (Jul-Sep), see Fig. 8. The modelled ETM at the surface moves a few kilometers upstream to km 120 with a concentration of 0.12 kg/m 3 without dumping (Fig. 8a) and 0.18 kg/m 3 with dumping ( Fig. 8b), corresponding to the large-scale trends observed in the measurements. The modelled surface ETM around 150 km decreases in magnitude to values that match OMES observations. Similar to the case with year-averaged discharge, the depth-averaged concentration observations are qualitatively reproduced in terms of the ETM location and magnitude only if dumping is included in the model. The depth-averaged ETM is located around km 100 with a concentration around 0.27 kg/m 3 .

Sensitivity to the settling velocity and erosion parameter
The characteristics of the sediment in the model are predominantly determined by the clear-water settling velocity w s,0 and erosion parameter M. Both parameters are highly uncertain and subject to natural variation. We therefore test the model sensitivity for settling velocities between 0.5 and For average discharge conditions in the ETM seaward of km 125 (Fig. 9a), the highest depth-averaged sediment concentrations are approximately 2 g/l and are found for   (Table 2), the coloured circles correspond to along-channel concentration profiles plotted in Fig. 10 M > 10 −2 s/m and w s,0 between 2 and 3 mm/s. The ETM landward of km 125 (Fig. 9b) shows the highest depthaveraged concentrations up to 9 g/l for a combination of M > 10 −2 s/m and large w s,0 . High settling velocities correspond to large vertical gradients in the sediment concentration. As a result, the maximum concentration near the bed is 22 g/l, which is significantly larger than the depth-averaged concentration. For summer discharge conditions, the ETM seaward of km 125 (Fig. 9c) shows the same behaviour for varying M and w s,0 as for average discharge conditions. However, depth-averaged concentrations are now up to approximately 1 g/l and therefore lower than for average discharge conditions. The ETM landward of km 125 (Fig. 9d) shows a different behaviour for varying M and w s,0 within the tested range. The highest depth-averaged concentrations of approximately 11 g/l are found for the largest values of M and w s,0 around 1-1.5 mm/s. Near the bed, maximum concentrations of up to 40 g/l are found in the ETM landward of km 125.
While the ETM seaward of km 125 is strongly affected by dumping of sediment, high concentrations can also be attained in this ETM if dumping is excluded from the model (not shown). This confirms the earlier observation that this ETM corresponds to a sediment trapping location, not just a plume of dumped sediment. Without dumping, concentrations in the ETM seaward of km 125 would be similar as in Fig. 9a and c for w s,0 approximately > 3 mm/s; dumping strongly affects the concentration in this ETM for w s,0 < 3 mm/s. The ETM landward of km 125 is not strongly affected by dumping for any combination of M and w s,0 .
To further illustrate the along-channel distribution of sediment for some of the cases with higher erosion parameter, Fig. 10 shows along-channel near-bed sediment concentrations for three situations, all with M = 0.1 s/m and different Q and w s,0 , as marked by the coloured circles in Fig. 9 (here near-bed concentration is defined as that found at the bottom boundary in the model). For all three plotted cases, the sediment concentration locally exceeds 10 g/l, concentrated around two ETM. Between these ETM concentrations are much lower with values around 100-200 mg/l. Even though concentrations are moderate in a large part of the estuary, we call these conditions hyperturbid, as the sediment concentrations are high over an along-channel distance of several tens of kilometers and have a visible

Water level and sediment concentration
Over the last decades, the Scheldt River Estuary has become deeper due to sand mining and channel deepening. Measurements and smoothed approximations of the widthaveraged depth in 1960 and 2010 are plotted in Fig. 11. Inspired by the along-channel pattern of deepening in the past, we define depth profiles of the form where H 1960 , H 2010 are the fitted depth profiles of 1960 and 2010 and α is a bed-profile parameter. For α = 0, we obtain the depth of 1960; for α = 1, we obtain the depth of 2010; and for α > 1, we obtain a depth larger than in 2010 by extrapolating the pattern of deepening between 1960 and 2010. We vary α between 0 and 2, keeping all other parameters the same as in the default experiment (Table 2). Model results show that channel deepening between 1960 and 2010 leads to an amplification of the M 2 water level, thereby confirming observations from measurements and earlier model studies Cai and Savenije 2013). Furthermore, we find a combination of amplification and damping of the M 4 water level with deepening. Figure 12a shows the M 2 and M 4 water level amplitude for α = 0, 1 and 2, together with observations of the tidal amplitude in 1960 and 2009. In order to compare the model result and measurements for 1960, the modelled water level for α = 0 is for a situation without sediment dumping, while dumping is taken into account for α > 0. Although the modelled M 2 tide is only calibrated for 2010 conditions, the M 2 tide for α = 0 shows good correspondence with the 1960 observations. As α increases, the M 2 tidal amplitude increases for all tested values of α. The M 4 tidal amplitude is overestimated compared to the measurements in both 1960 and 2010 conditions. Between 1960 and 2010, the measurements show only very minor changes in M 4 tidal amplitude, with amplification upstream from km 130 and reduction of the amplitude downstream from km 130. For α increasing from 1 to 2, the model does capture a trend similar to what was observed between 1960 and 2010, with increasing M 4 amplitude upstream from km 140 and decreasing amplitude elsewhere.
The maximum sediment concentrations become lower as a result of channel deepening. This is illustrated in Fig. 12b, which shows the depth-averaged sediment concentration as a function of x and the bed-profile parameter α for year-averaged discharge. Sediment dumping is now taken into account and has its default value for all α. For α = 0, the figure shows two ETM around km 80 and 115. As α increases (moving up along the vertical axis), the sediment concentrations in the ETM become lower. From approximately α > 0.9, the ETM at km 115 starts to disappear and is replaced by an ETM at km 150. As α increases further, the concentrations in this new ETM at km 150 also decrease. Repeating these model experiments without sediment dumping (not shown) yields similar results; however, the ETM at km 80 is much weaker. Hence, regardless of sediment dumping, the effect of deepening alone leads to an upstream shift of the ETM from km 115 to km 150 and lower maximum sediment concentrations.

Sensitivity to the settling velocity and erosion parameter
The effect of channel deepening for other values of the settling velocity w s,0 and erosion parameter M is generally consistent i.e. showing decreasing concentrations with

Fig. 13
Ratio of the maximum depth-averaged concentration for α = 2 (deepened) divided by that for α = 1 (2010) in the areas seaward and landward of km 125 for average and summer discharge conditions, plotted for a range of values of w s,0 and M. Blue colours denote that the maximum concentration is lower for α = 2 than for α = 1. Red colours denote that the maximum concentration is higher for α = 2 than for α = 1. The contour line indicates no change between the scenarios. The circles indicate the default parameter settings (Table 2) Increasing concentrations with deepening are found for a combination of a high erosion parameter and high settling velocity in the ETM seaward of km 125. For these settings, the estuary is already highly turbid for 2010 depth conditions (see Section 3.3), so the increasing concentration with deepening does not signify a transition from low to high sediment concentrations. In the ETM landward of km 125, increasing concentrations are only found for the year-averaged discharge case for settling velocities between approximately 1 and 2 mm/s and M > 10 −3 s/m. The increase is, however, less than a factor two and also does not indicate a transition from low to high sediment concentrations.

Analysis
In order to gain a better understanding of and confidence in the presented results, we investigate the physical processes underlying the sediment dynamics in the model. Following the approach taken by Dijkstra et al. (2019), the model results before and after deepening are analysed on the basis of two aspects: 1. along-channel suspended sediment transport; and 2. vertical resuspension These aspects are quantitatively expressed in terms of the transport capacity, erodibility and dimensionless erosion parameter, which are introduced in Section 5.1. Next, in Section 5.2, we analyse the sediment dynamics in the 2010 case. This is used to explain the sensitivity to the erosion parameter in Section 5.3 and to deepening in Section 5.4.

Transport capacity, erodibility and dimensionless erosion parameter
In order to analyse the vertical resuspension of sediment, we look closer at the formulation for erosion E in iFlow, which reads as Here, M is a prescribed erosion parameter, τ b is the bed shear stress and f is the tidally averaged erodibility. This erodibility indicates the tidally averaged amount of sediment on the bed. The erodibility is a number between 0 and 1, where 0 means that no sediment is available for erosion during the entire tidal cycle and 1 indicates that easily erodible sediment is available at the bed during the entire tidal cycle. A number between 0 and 1 indicates that sediment is available at the bed during some part of the tide. A formal mathematical definition is provided by Brouwer et al. (2018).
Using this erosion formulation, one can define a dimensionless erosion parameterẼ (see Dijkstra et al. (2018)), which expresses the ability of the flow to resuspend sediment from the bed. For our model, it is mathematically expressed as where the clear-water settling velocity w s,0 and gelling concentration c gel are constants in our model. Hence, the along-channel variation ofẼ expresses along-channel variations in the bed shear stress τ b . In order to analyse changes in sediment dynamics, one could directly investigate the changes in the sediment transport. The disadvantage of this is that the sediment transport is typically large near the ETM and small in areas with little sediment. Hence, when the ETM moves to a location where little sediment was available previously, this appears as a large change in the sediment transport. The changes in the sediment transport then simply reflect the changes in ETM location, not the changes in hydrodynamic forcing that caused the ETM to move. To circumvent this, we use the concept of transport capacity. A formal mathematical definition of the transport capacity is provided in by Chernetsky et al. (2010) or Dijkstra et al. (2019). More intuitively, the transport capacity is defined as the sediment transport that would occur if a uniform layer of sediment was added on the bed everywhere in the estuary, given the modelled hydrodynamic conditions (flow velocity, turbulence field) and sediment parameters (effective settling velocity, erosion parameter). It therefore indicates the tidally averaged redistribution of a uniform layer of sediment on the bed. Trapping of sediment, here defined as a local maximum of the erodibility, is indicated by a convergence of the transport capacity.
In iFlow, the transport capacity can be subdivided into various physical contributions. The most dominant contributions in the Scheldt River are the following.
-The external M 4 tide contribution is due to tidal asymmetry caused by the M 2 tide and M 4 tide entering the estuary at the mouth. The contribution to the M 4 tide is generated outside the estuary on the shallow shelf and propagates through the estuary, causing asymmetry in the velocity and sediment resuspension during ebb and flood and therefore net sediment transport. -The tidal return flow contribution is the transport capacity due to Stokes drift and the corresponding return flow. The Stokes drift is associated with sediment import. At least partly compensating this import, the return flow velocity contains a subtidal contribution which typically causes export of sediment. Additionally, the return flow velocity has an M 4 contribution, which may cause import or export of sediment, depending on the phase lag with the M 2 tide. -The velocity-depth asymmetry contribution is the transport capacity due to the asymmetry of the tide that is created because the depth is different during ebb and flood. This yields different velocity profiles during ebb and flood and hence asymmetric sediment resuspension and transport. Whether this effect is importing or exporting sediment depends on the phase difference between the M 2 velocity and surface elevation. -The sediment advection contribution represents the transport due to spatial settling lag (see e.g. Van Straaten andKuenen 1957, De Swart andZimmerman 2009). This contribution tends to transport sediment towards along-channel minima in the tidal velocity amplitude. -The river contribution consists of two parts: the riverinduced flushing of tidally resuspended sediment and the transport due to the tidal asymmetry caused by the tide-river interaction. Both contributions cause an export of sediment The above contributions to the sediment transport capacity contain correlations of subtidal velocities and sediment concentrations, as well as correlations of tidally varying velocities and sediment concentrations, known as tidal pumping (see e.g. Burchard et al. (2018) for a review). Tidal pumping can be further related to M 2 -M 4 tidal asymmetry, which is contained in the above contributions by the external M 4 tide, tidal return flow and velocity-depth asymmetry. Other contributions to tidal pumping are related to the asymmetry caused by the subtidal flow, contained in the river contribution, tidal return flow and velocitydepth asymmetry, and to spatial gradients, contained in the sediment advection contribution.

Analysis of the 2010 case
Analysis of the sediment transport capacity provides more insight into the ETM near km 80 and 150. Figure 14a shows the transport capacity for the 2010 case with average discharge (Q = 110 m 3 /s) with sediment dumping and is almost the same for the case without dumping. The total transport capacity (black line) shows the ability of the flow to transport sediment upstream (positive) or downstream (negative). Sharp jumps occur in the transport capacity at km 95 and 123, due to the inflow of fresh water from tributaries. The ETM correspond to the two convergence zones near km 80 and 150, indicated by the numbers in (a) Total transport capacity and decomposition of the transport capacity into its five most important contributions. The vertical markings and numbers denote trapping locations.
(c) Dimensionless erosion parameter. The reason why convergence of transport occurs, follows from the balance of dominant physical mechanisms, which are also shown in Fig. 14a. There are two dominant exporting (i.e. negative) contributions. The river discharge dominates the sediment export from the estuary for x > 80 km. For x < 80 km, sediment export is dominated by the sediment transport due to M 2 -M 4 tidal asymmetry that is related to the externally forced M 4 tide. Three contributions are important for import (i.e. positive): sediment advection (or spatial settling lag), tidal return flow and velocitydepth asymmetry. The latter two are associated with M 2 -M 4 asymmetry of the tide. Hence, sediment transport due to M 2 -M 4 tidal asymmetry is important, but not all contributions to this asymmetry lead to sediment import. The resulting combined transport by the M 2 -M 4 asymmetry is a small import of sediment for x > 40 km.
While the sediment transport capacity shows that the observed ETM near km 80 and 150 are results of sediment trapping, it cannot explain the ETM observed at the surface near km 115 (see Fig. 7b). This ETM is not directly related to a trapping zone but results from a large resuspension of sediment. This results from the combination of a sufficiently large availability of sediment, expressed by the erodibility (Fig. 14b) and a relatively large erosion, expressed by the tidally averaged dimensionless erosion parameter (Fig. 14c). The sediment available at the bed is suspended high up in the water column and is therefore observed as an ETM at the surface.

Analysis of the sensitivity to the erosion parameter
The erodibility helps to explain why the model results are sensitive to the erosion parameter, as was found in Fig. 9. Near the ETM at km 80 and 150, the erodibility (Fig. 14b) equals one. This means that sediment is always available on the bed at these locations. The maximum sediment concentration in these ETM is limited by the ability of the flow to resuspend sediment (i.e. erosion limited conditions). Hence, the sediment concentration at these locations increases when the erosion parameter is increased. Results indicate that erosion limited conditions prevail for average discharge conditions at these ETM locations for M up to 10 −2 s/m if the settling velocity exceeds 2 mm/s (not shown). This is over 10 times the erosion parameter used in this study. If dumping were not included in the model, less sediment would be available at the bed but erosion limited conditions still prevail for M up to 10 −3 s/m (for a settling velocity of 2 mm/s) to 10 −2 s/m (for a settling velocity of 4 mm/s). Therefore, regardless of dumping, the development of hyperturbid conditions in the Scheldt, within our model, is mainly controlled by the exchange of sediment between the water column and the bed, parametrised by the erosion parameter.

Analysis of the response to deepening
To analyse why deepening leads to lower sediment concentrations, the above analysis is repeated for the case α = 2 (i.e. deepening), and results are compared to the case α = 1 (i.e. 2010). The total sediment transport capacity (Fig. 15a) shows only minor changes due to deepening in the seaward half of the estuary. In addition, the locations of the trapping zones near km 80 and 150 (indicated by the numbers) change only little, moving slightly further apart. Only between the two trapping locations does the transport capacity change, leading to an increasing convergence of sediment near the trapping locations. Corresponding to this, the transport capacity diverges between the trapping zones, and less sediment is found between km 80 and 150.
The relatively minor changes in the transport capacity are caused by a mixed response of the underlying mechanisms to deepening (not shown). The sediment transport related to the external M 4 tide becomes more exporting. Additionally, the sediment transport related to tidal return flow becomes less importing between km 0 and 80. The transport due to other mechanisms does not change much between km 0 and 80. Between the ETM at km 80 and 150, the increasing divergence is found in the tidal return flow, as well as the velocity-depth asymmetry and spatial settling lag.
As the trapping locations do not change much due to deepening, the explanation for the lower maximum concentration follows from the erodibility and dimensionless erosion parameter. Near the ETM at km 115, increasing divergence of sediment transport capacity means that less sediment is available, which is expressed by a lower erodibility, see Fig. 15b. As a result, the sediment concentration in this ETM decreases with deepening. In the ETM near km 80 and 150, the erodibility remains equal to one, since the sediment transport still strongly converges in these (a) Total transport capacity for the 2010 case ( 1) and the case for 2.

Fig. 15
Transport capacity, erodibility and dimensionless erosion parameter for the cases α = 1 (2010) and α = 2 for average discharge (Q = 110 m 3 /s) with dumping areas. Therefore, the maximum concentration in these ETM is restricted by the dimensionless erosion parameter. The dimensionless erosion parameter (Fig. 15c) decreases after deepening near km 80. This is related to a decrease in the bed shear stress, which is caused by a decrease in the M 2 tidal velocity in response to deepening. This explains why the sediment concentration becomes lower near km 80 after deepening. Near km 150, the trapping zone moves slightly upstream (see above), where the dimensionless erosion parameter is smaller. This explains the lower sediment concentrations in this ETM.

Comparison with the Ems River Estuary
The main motivation to study whether hyperturbid conditions can develop in the Scheldt as a consequence of deepening is the development of hyperturbid conditions as a consequence of deepening in the Ems. However, we found that the effects of deepening in the Ems and Scheldt are different. These differences are explained below.
Using the same iFlow model as in this study, the observed transition to hyperturbid conditions in the Ems was qualitatively reproduced by Dijkstra et al. (2019). This was done by calibrating the model to a situation representing 1965 (before hyperturbid conditions developed) and then changing the depth to conditions representing 2005 (after hyperturbid conditions developed). It was concluded that the sediment dynamics in the Ems is supply limited i.e. the erodibility is smaller than 1, and sediment concentrations are restricted by the ability of the estuary to import sediment. It was found that the most important physical mechanisms responsible for import in the Ems are the M 2 -M 4 tidal asymmetries related to the externally generated M 4 tide and tidal return flow. With deepening, both mechanisms become more importing. The additional import of sediment leads to sediment-induced stratification, which leads to damping of turbulence. This in turn leads to a further increase of sediment import due to the M 2 -M 4 tidal asymmetry.
The sediment dynamics in the Scheldt behaves differently to the Ems in two aspects. Firstly, it is concluded in this study that the most intense ETM locations in the Scheldt are erosion limited. Hence, the maximum sediment concentration is not restricted by the ability of the estuary to import sediment but by the ability to resuspend sediment from the bed. With deepening, it is found that the bed shear stress decreases in the Scheldt in the ETM near km 80, explaining the lower sediment concentrations. A decrease in the bed shear stress with deepening is also found locally in the Ems but does not restrict the sediment concentration there.
Secondly, deepening does not lead to an increasing import by all sediment transport mechanisms. On the contrary, the transport related to the externally generated M 4 tide becomes more exporting, while transport related to the internally generated M 4 tide become less importing in a large part of the estuary. As deepening does not lead to an increasing sediment concentration, a feedback between sedimentinduced damping of turbulence and sediment import, as in the Ems, cannot develop in the Scheldt.
The underlying reasons why deepening lead to increasing sediment import in the Ems but not in in the Scheldt follow from a complex interplay between the non-linear generation of the M 4 tide and the phase difference between the M 2 and M 4 tide. Nevertheless, one of these reasons can be understood intuitively. At the mouth of the estuary (taking the mouth of the Ems River at Knock), the phase difference between the M 2 and M 4 tide is −4 degrees in the Scheldt and −172 degrees in the Ems. This is a difference of almost 180 degrees, explaining why the effects of the externally generated M 4 tide on transport are almost completely opposite in the Scheldt and the Ems. Deepening leads to amplification of the externally generated M 4 tide inside the estuary in both the Scheldt and Ems and hence to more export in the Scheldt and more import in the Ems.

Effect of model simplifications
As this model study is highly idealised, there are many physical processes that are not included and some processes that are not represented accurately. Nevertheless, as we have studied the sensitivity of the model to parameter variations and have investigated the most essential physical mechanisms in the model, it is possible to discuss the robustness of the results with respect to these model simplifications.
One of the main discrepancies between the model results and observations is the amplitude of the M 4 water level and the velocity (e.g. Fig. 6), which are overestimated by more than a factor 2 in the upstream part of the estuary. The reason for this overestimation is unknown and could potentially be related to an oversimplification of the geometry, such a neglecting tidal variations in estuary width (see e.g. Friedrichs and Aubrey 1994), or oversimplification of bed friction and eddy viscosity as a time-independent quantities, such that interactions of tidal constituents to create tidal damping is not fully included (see e.g. Godin 1999).
The effect of the overestimation of the M 4 tide has been tested in our model by artificially scaling to match the observations. Results show only slight variations in sediment concentration and the sensitivity studies of M and w s , such that the primary conclusions are not affected. This ad-hoc sensitivity test suggests that the mismatch between measured and modelled M 4 tide does not appreciably alter results. This relative unimportance of the error in the M 4 tide is because of several reasons. Firstly, since the estuary is erosion limited, the spatial variation of the M 2 tide (and modelling that correctly) is most important. Secondly, the phase difference between the MM 2 and M 4 tide determines the direction of the net transport, and this phase difference is captured correctly by the model. Finally, our conclusions consider the relative change of the concentration as a consequence of the relative change in dynamics after deepening. Hence, the absolute values of the M 4 tidal amplitude are less important as long as relative changes are captured. Hence, our conclusions seem robust to the errors in the modelled M 4 tide, but further study is strongly recommended to verify what contributions to the M 4 tide are missing in this model.
A process that is often considered to be important in the Scheldt is flocculation (e.g. Chen et al. 2005). Flocculation affects the settling velocity of sediment. In this study, we have shown that the results are robust for large changes in the settling velocity in the entire estuary. Hence, large-scale changes in the settling velocity due to changing flocculation properties do not affect our conclusions. The remaining uncertainty is related to spatial or temporal variations of the settling velocity, which are not taken into account.
As erosion is the most restricting process to the sediment concentration, and as it is found that higher sediment concentrations may occur in the Scheldt for larger values of the erosion parameter, the erosion formulation requires most direct attention of further research. The erosion formulation used in this study is based on Partheniades' formulation, which is also used in many state-of-the-art complex models. In this study, we have simplified this formulation by ignoring the critical shear stress and omitting tidal variations of the sediment availability at the bed. While these simplifications likely have important consequences for the quantitative results, they do not change the qualitative conclusions. The main source of uncertainty is the erosion parameter. It remains unknown whether the value of this parameter changes on the long timescale or as a response to deepening. It was identified in this study that hyperturbid conditions can occur in the Scheldt if the value of the erosion parameter is increased by one to two orders of magnitude. Further research is needed to investigate if this is possible.
Our conclusion is further supported by studies using complex models of the Scheldt. Using a depth-averaged Delft3D model, Van Kessel et al. (2008) investigated the effect of the second deepening campaign (1997)(1998) and found that this deepening should lead to lower suspended sediment concentrations. Using a three-dimensional model, Vandenbruwaene and Stark (2018) show that the tide in the estuary became less flood dominant due to deepening since the 1930s, also suggesting less sediment import due to deepening.

Conclusion
We have investigated the hypothesis of  that the Scheldt may become hyperturbid as a response to deepening. To this end, we have used the iFlow model to investigate the dynamic equilibrium sediment concentration in the Scheldt for a case representing conditions of 2010 and a range of cases with higher and lower bed levels, keeping all other parameters the same. In order to draw robust conclusions, all cases have been tested for a range of values of uncertain parameters and the physics underlying the sediment dynamics has been investigated.
For the conditions representing 2010, the modelled sediment concentrations qualitatively reproduce observed longterm average ETM locations and sediment concentration magnitudes. From the analysis, it is found that the most intense ETM locations in the Scheldt are erosion limited i.e. the maximum sediment concentration is restricted by the ability of the flow to resuspend sediment from the bed, not by the availability of sediment. Hence, model results are sensitive to the quantities that control the amount of resuspension, which are mainly the bed shear stress and an erosion parameter.
Deepening of the estuary in the model generally leads to lower maximum sediment concentrations in the Scheldt. When investigating the sensitivity to varying parameter values, some parameter settings were identified where the maximum concentration increases with deepening but such increase is minor and does not lead to the development of hyperturbid conditions. The analysis shows that the flow velocity and hence the bed shear stress at the ETM locations generally decrease with deepening. This results in a reduction of resuspension, which in turn results in lower sediment concentrations. Furthermore, deepening does not lead to a clear trend of increasing sediment import in the Scheldt. Overall, deepening leads to less import in the most seaward part of the estuary and more convergence of sediment around the ETM.
Based on these results, we suggest to reject the hypothesis of  that channel deepening alone may lead to development of hyperturbid conditions in the Scheldt. By combining the model results, sensitivity analysis and understanding of underlying processes, we argued that this is a robust conclusion, even though the model used is highly idealised. To further verify this conclusion, it is recommended to investigate some processes that are missing or inaccurately represented by the model. Firstly, this concerns the M 4 tide, which is overestimated in the model. Secondly, this concerns the parametrisation for erosion. High sediment concentrations were found in the model of the Scheldt when the erosion parameter is increased by one or two orders of magnitude compared to its calibrated value and it remains unknown whether such an increase of the erosion parameter could occur.
At the bed, we impose a non-linear friction law for the water motion by imposing that the bed shear stress equals ρ 0 s f u bed , where u bed is the along-channel velocity near the bed and s f is the partial slip parameter that depends on the depth-averaged velocity U as Here, c v,3 (z * 0 ) is a prescribed function of z * 0 and c D is a reduced drag function that depends on the sediment stratification near the bed according to Adams and Weatherly (1981), Friedrichs et al. (2000), and Wang (2002) where A = 5.5 is an empirical parameter and Rf bed is the flux Richardson number near the bed. We restrict Rf bed to a maximum value of 2 to ensure that the friction does not become much smaller than observed in laboratory studies. For sediment erosion and resuspension, the bed-shear stress τ b in Eq. 1 is parametrised as where the parameter s s equals s f (12) but assuming c D = 1. Sediment erosion is therefore not affected by reduced drag due to near-bed stratification.

Appendix B: New model additions related to tributaries
The discharge of water into the estuary by tributaries is added to the model as depth-integrated point sources by adding source terms to the depth-averaged continuity equation, i.e.
where B is the width, H is the depth below mean sea level (MSL) z = 0, R is the reference level above MSL, ζ is the surface elevation, u is the horizontal flow velocity and subscripts x and t indicate derivatives with respect to alongchannel distance and time, respectively. The source terms S Q,n represent the discharge of tributary n at location x n , δ denotes the Dirac delta function i.e. indicating a source at one point. These sources lead to an additional contribution to the first-order residual water motion (see Dijkstra et al. 2017) that can be analysed separately. Sources of sediment enter into the sediment concentration is computed using an equation for mass conservation and an equilibrium condition, which requires that the crosssectionally integrated sediment concentration does not vary on the subtidal timescale. This condition is equivalent to requiring In this equation, · denotes tidal averaging, c is the sediment concentration, K h is the horizontal eddy diffusivity, E denotes erosion or resuspension from the bed and D denotes deposition. Fluvial sources of sediment, dreding and dumping are added to the model by adding a source or sink of sediment S, which consists of pulses that equal the rate of sediment added (positive) or extracted (negative) from the system at the confluences, dredging and dumping locations. The model computes itself how the sources and sinks are distributed over the water column and the bed in such a way that the model remains in dynamic equilibrium.