Critical dependence of morphodynamic models of fluvial and tidal systems on empirical downslope sediment transport

The morphological development of fluvial and tidal systems is forecast more and more frequently by models in scientific and engineering studies for decision making regarding climate change mitigation, flood control, navigation and engineering works. However, many existing morphodynamic models predict unrealistically high channel incision, which is often dampened by increased gravity-driven sediment transport on side-slopes by up to two orders of magnitude too high. Here we show that such arbitrary calibrations dramatically bias sediment dynamics, channel patterns, and rate of morphological change. For five different models bracketing a range of scales and environments, we found that it is impossible to calibrate a model on both sediment transport magnitude and morphology. Consequently, present calibration practice may cause an order magnitude error in either morphology or morphological change. We show how model design can be optimized for different applications. We discuss the major implications for model interpretation and a critical knowledge gap.

: Concept of the analytical model. The cross-section is three grid cells wide with a bed level difference between the middle grid cell and the surrounding cells as an initial perturbation. The numerical channel is also based on this concept. a) Definition of the flow velocity and transport vectors, channel width (W ), initial channel depth (h i ), bed level difference (dh), and transverse slope (dz/dy). b) The perturbation decays when transverse sediment transport is larger than the difference between incoming and outgoing sediment transport. The middle grid cell will accrete, while the surrounding cells will erode till the average bed level. c) The perturbation grows when transverse sediment transport is smaller than the difference between incoming and outgoing sediment transport. The middle grid cell will incise further, while the surrounding cells will accrete.
Supplementary Figure 3: The trend in equilibrium width-to-depth ratios with increasing sediment mobility, for three different depths of the initial perturbation (dh), resulting from the analytical model. Colors indicate the non-linearity of sediment transport (k). Solid lines indicate default slope effect (α I = 1.5), dashed lines indicate an increased slope effect (α I = 10). Width-to-depth ratios to the left of these lines will result in a decay of the initial perturbation, while ratios towards the right will result in a growth.
Supplementary Figure 4: Channelization factor resulting from the analytical model, plotted against the depth of the initial perturbation. Width-to-depth ratios lower than these lines will result in a decay of the initial perturbation, while higher ratios will result in a growth. of slope effect and sediment transport predictors. Maps on the horizontal axis have an equal slope effect, with slope effect increasing downwards. The α I is the input parameter of the method of Ikeda 2 , while the α K is the input parameter of the method of Koch and Flokstra 3 . The models in the first two columns were run with the VR sediment transport predictor, while the models in the last two columns were run with the EH predictor. The average sediment transport rates plotted in Supplementary Figure 10 were computed for all model runs over a cross-section at 20 km,    Difference in direction of sediment transport between the models with different slope parameterizations. This distribution shows the relative abundance of these differences for all grid cells in the model.

Supplementary Tables
Supplementary Table 1 shows the results of the literature inventory in typical geomorphology journals. For each study, the morphodynamic model and the magnitude of the slope factor is indicated, and whether this study mentions or discusses the effect on morphology of the slope parameter (yes = 1, no = 0). The environment that is modelled can be erosional (1), depending on a large-scale balance between erosion and deposition (2), or depositional (3). Models are included in this literature study when they either have an upstream river boundary (1), a river boundary as well as a downstream tidal boundary (2), or only a tidal boundary (3). Lastly, it is noted whether the model considers suspension (1) or treats all sediment transport as bedload (0).  VR makes a distinction between bed load and suspended load transport, by imposing a reference height, below which sediment transport is treated as bed load and everything above this height is treated as suspended load. Gravity only acts on the bed load, which is calculated as follows: . As a result, the sediment transport rate is related to flow velocity to the power of 3, which determines the non-linearity of the sediment transport predictor. However, since this predictor also includes a critical flow velocity, the relation between flow velocity and sediment transport will be more non-linear near the beginning of motion.

Supplementary
EH is a total load predictor (q t ), and unlike VR, it does not include a critical velocity or critical shear stress: where α = a calibration coefficient in the order of 1. Here, the sediment transport rate is related to flow velocity to the power of 5.
The general sediment transport predictor in Delft3D is based on the predictor of Meyer-Peter Mueller 74 : where b and c are user defined parameters, which determine the non-linearity of the sediment transport predictor and the addition of a critical sediment mobility. The sediment mobility θ, a dimensionless form of the bed shear stress, reads: When the magnitude of the bed load or total load sediment transport is calculated parallel to the flow velocity, the direction and magnitude of the transport vector is adjusted for bed slopes.
For transverse slopes, the two commonly used parameterizations are the predictor of Koch and Flokstra (KF) 3 (ISlope = 3) and Ikeda (IK) 2 (ISlope = 2). The main difference between both options is in the calculation of the transport vector (Fig. 2 in main text). For KF the direction of sediment transport is corrected for transverse gradients by rotating the transport vector based on the user-defined factors α K and β K : For IK an additional transport vector is calculated perpendicular to the flow direction, based on the input parameter α I : where q = sediment transport load [m 2 s −1 ] in the streamwise (s) or transverse (n) direction, As a result, the IK method increases the direction and total magnitude of sediment transport when a transverse slope is present, while for KF only the direction is changed. Another difference is that the IK method uses a critical shear stress, which is absent in the KF method. The default value of α I in Delft3D is set to 1.5, while the parameter α K is not defined in the model, but should be 1.5 according to Koch and Flokstra 3 . The method of calculating the sediment transport vector in both slope options therefore has major implications for calibrating models with the transverse slope parameter. By increasing the α I in the IK method by a factor of ten for example, the amount of downslope sediment transport is also increased by a factor of ten, which increases the total sediment transport significantly (Fig. 2 in main text). With the KF method sediment transport is not increased, but here, decreasing the α K to values reported in literature 12,13,65 could easily result in more downslope sediment transport than streamwise sediment transport.
It is possible to compare the effect on resulting morphology of using different slope predictors by requiring either the magnitude or the direction of transverse sediment transport to be equal.
When assuming an equal magnitude, the method of KF needs to be corrected for a given slope and sediment mobility. Using Equation 5 and Equation 6 with a β K of 0.5 it follows that: The resulting relation between α I and α K is plotted in Supplementary Figure. 1 for four combinations of transverse slope and sediment mobility. When assuming equal direction of sediment transport, it follows that: which is shown as the linear solution in Supplementary Figure 1.
Supplementary Note 2 To help identify the cause of the overdeepening of channels in numerical models, we compare the balance between incision and transverse sediment transport in a straight river channel in Delft3D with an analytical model of a channel cross-section with the same characteristics. Since we only consider a cross-section, the streamwise sediment transport is in balance with the constant flow conditions and the model does not account for deposition along a river reach.
The analytical model consists of three grid cells in cross-section, with an initial bed level difference between the middle cell and the surrounding cells, representing a disturbance that either decays or grows by incising further (Supplementary Figure 2). The aim of this model is to find the equilibrium width-to-depth ratio at which incision is equal to transverse sediment transport, and how this ratio depends on flow conditions, sediment transport processes, and size of the disturbance.
The model first calculates upstream flow characteristics and corresponding sediment transport rate based on the input parameters, which are a constant Chezy coefficient for friction (C), channel slope (S), grain size (D 50 ), the non-linearity of the sediment transport predictor (k), and a height difference (dh). We assume a constant specific discharge such that the relation between channel width (W ) and discharge (Q) is linear: The upstream flow velocity (u i ) and water depth (h i ) are calculated by iteration, using the following equations for flow velocity: The upstream sediment transport rate (q i ) is based on the same general sediment transport predictor as in Delft3D: Then, flow characteristics and sediment transport fluxes are calculated for the cross-section under consideration, based on the height difference between the middle grid cell (h 2 ) and the outer two grid cells (h 1 , h 3 ) (Supplementary Figure 2a). It is assumed that the average water depth at the cross-section is equal to the initial water depth, which leads to: The sediment transport rate for each cell is then calculated with Equations 4, 10 and 12, but with the specific water depths. The sediment transport rate towards the middle cell as a result of the transverse slope (q n ) is based on the method of Ikeda 2 : where β = transverse slope parameter, which is based on α I from Equation 6. The transverse slope is defined as the height difference between two cells divided by the width of one grid cell, which is the same method as in Delft3D.
A balance between incision and downslope sediment transport is assumed when the difference between the upstream sediment transport and the sediment transport rate for the middle grid cell is equal to the total downslope sediment transport: When the transverse sediment flux is larger, there is sedimentation and the perturbation will likely decay (Supplementary Figure 2b), while when the transverse sediment flux is smaller, the grid cell is incised and the perturbation will grow (Supplementary Figure 2c). Using Equations 4, where W eq = width of the channel when incision is equal to the transverse sediment transport.
The equilibrium width-to-depth ratio is now a function of the size of the disturbance, sediment mobility and the non-linearity of sediment transport. All other parameters influence this equilibrium by changing the sediment mobility. In further analyses we first assume a constant channel slope of 0.5 mm m −1 , a Chezy coefficient of 40 √ m/s −1 , a ratio between channel width and discharge of 12.5, and a grain size of 0.5 mm.
With increasing sediment mobility, the equilibrium width-to-depth ratio decreases exponentially (Supplementary Figure 3), which means that at higher sediment mobility a channel is more likely to incise. A higher non-linearity of the transport predictor causes a higher sediment transport rate, and therefore results in more incision and a lower equilibrium width-to-depth ratio at any sediment mobility. Increasing the transverse slope parameter has the opposite effect, since more sediment is transported downslope which counteracts incision. Increasing the depth of the initial perturbation also decreases the equilibrium width (Supplementary Figure 3), since deeper channels attract more flow and therefore need more downslope sediment transport to counteract this. However, this influence is less than changing the bed slope effect or the non-linearity.
To be able to show the effects of height of the perturbation and the other parameters that influence sediment mobility, the width-to-depth ratio is multiplied by the square root of the sed-iment mobility divided by the slope parameter, which is the ratio that describes the slopes of the graphs in Supplementary Figure 3. We call the resulting parameter the channelization factor, since it describes the balance between the tendency to enhance perturbations determined by the widthto-depth ratio, and the bed slope effect that counteracts incision. This balance thereby controls the formation of channels. As a result, Supplementary Figure 4a shows how models with varying slope effect and sediment mobility collapse when plotting this factor against height of the perturbation.
Again, a higher non-linearity of sediment transport results in a growth of the perturbation at lower width-to-depth ratios. Higher Chezy values, and thus lower friction, also results in a growth of the perturbation at lower width-to-depth ratios when increasing the depth of the perturbation, but less dramatically. However, negative perturbations, i.e. when the middle grid cell is higher than the surrounding cells, need higher width-to-depth ratios for the perturbation to grow. Increasing the channel slope or decreasing the ratio between discharge and channel width shows the same trend.
Since the analytical model identifies the equilibrium channelization factor, perturbations in numerical models plotted below this line should theoretically decay, while models plotted above the line should have growing perturbations (Supplementary Figure 5). With the default value for the slope effect (α I = 1.5), the VR models corresponded reasonably well with the analytical model, since the transition from a dampened system towards a channel where the perturbation grows is around the theoretical equilibrium line (Supplementary Figure 5a). There was no effect of the depth of the initial perturbation in the numerical model. However, with increased slope effect, the numerical models significantly deviated from the analytical model. Here, the numerical model with wider channels required a disproportionately larger slope effect to dampen the initial perturbation (more than 30 times higher than the default factor as opposed to 4 times the default in the analytical model). On the other hand, the initial perturbation in models with the models with the EH predictor immediately decayed (Supplementary Figure 5b), until the channel has a width-todepth ratio around 36, which is more than 15 times higher than the theoretical model. This behavior was very similar to that of the models with the general predictor (Supplementary Figure 5c,d).
Even when sediment transport is related to flow velocity to the power of 10, perturbations did not start to grow at a lower width-to-depth-ratio, while this was expected based on the analytical model. The two slope parameterizations differed only slightly and removing the critical sediment mobility from the generic transport predictor had no effect on equilibrium morphology. These results demonstrate a stronger tendency to incise in the numerical model with VR than expected from theory, and a weaker tendency to incise in numerical models with EH.
Supplementary Note 3 The fifth set of models is a detailed case study of the topographically forced Western Scheldt estuary in the Netherlands, to test the sensitivity of a calibrated model with two different slope parameterizations in comparison with measured bathymetry. This topographic forcing is typical for many natural and engineered systems and is important because it limits free bar and pattern formation, rendering models less sensitive in large-scale pattern to chosen parameterizations. Here, we focused on differences in local sediment transport dynamics in two model runs with different slope predictors that showed the same large-scale morphology in view of the need to predict sediment transport rates for fairway maintenance dredging. After 10 years of morphological development, these models reproduced the cumulative slope distributions that were closest to the actual morphology of the Western Scheldt that was used as input (Supple-25 mentary Figure 8c). The models had a strong slope effect, namely an α I of 30 and an α K of 0.05, which again shows that a higher than physical slope effect is needed when calibrating the Western Scheldt model on existing morphology.
While large-scale morphology is similar between both models after ten years ( Supplementary   Figure 8d), the dynamics differ in local sediment transport. The model with the IK method has higher bed load transport rates on steeper slopes, while the model with the KF method has higher transport rates on lower slopes (Supplementary Figure 11). Furthermore, there is a significant difference in direction of the transport vectors in more than half of all grid cells in the model (Supplementary Figure 11d), which is independent of slope. These differences in direction and magnitude imply locally channels can be orientated differently and location and speed of bank erosion will differ. For fairway maintenance dredging this means that predicted time scales can significantly differ when models are calibrated with a different slope parametrization on the same measured morphology.