Effects of mud supply on large-scale estuary morphology and development over centuries to millennia

Alluvial river estuaries consist largely of sand but are typically flanked by mudflats and salt marshes. The analogy with meandering rivers that are kept narrower than braided rivers by cohesive floodplain formation raises the question of how large-scale estuarine morphology and the late Holocene development of estuaries are affected by cohesive sediment. In this study we combine sand and mud transport processes and study their interaction effects on morphologically modelled estuaries on centennial to millennial timescales. The numerical modelling package Delft3D was applied in 2-DH starting from an idealised convergent estuary. The mixed sediment was modelled with an active layer and storage module with fluxes predicted by the Partheniades–Krone relations for mud and Engelund–Hansen for sand. The model was subjected to a range of idealised boundary conditions of tidal range, river discharge, waves and mud input. The model results show that mud is predominantly stored in mudflats on the side of the estuary. Marine mud supply only influences the mouth of the estuary, whereas fluvial mud is distributed along the whole estuary. Coastal waves stir up mud and remove the tendency to form muddy coastlines and the formation of mudflats in the downstream part of the estuary. Widening continues in estuaries with only sand, while mud supply leads to a narrower constant width and reduced channel and bar dynamics. This self-confinement eventually leads to a dynamic equilibrium in which lateral channel migration and mudflat expansion are balanced on average. However, for higher mud concentrations, higher discharge and low tidal amplitude, the estuary narrows and fills to become a tidal delta.


Introduction
Sandy river estuaries with continuously migrating channels and bars have great and often conflicting economic and ecological value.These estuaries are typically dominantly built of sand, but mud and salt marshes also form significant parts of these systems.Mud plays a critical role in ecological restoration measures and harbour maintenance, but it is rarely taken into account in numerical morphological models.Due to human interference, mud concentrations have increased far above the desired values in many estuaries (Winterwerp, 2011;Van Maren et al., 2016).Mud problems arise from pollutants attached to clay particles, mud deposits covering benthic species, rapidly siltating harbours and channels and changing hydrodynamic and morphodynamic conditions by higher resistance against erosion.This raises questions about the effects of mud on large-scale estuary morphology in natural alluvial systems as a control for cases with human interference.
In rivers, the formation of cohesive floodplains with mud and vegetation causes river channels to be narrower and deeper than in systems with only sand given otherwise equal conditions (Tal and Paola, 2007;Kleinhans, 2010;Van Dijk et al., 2013;Schuurman et al., 2016).This results from a dynamic balance between floodplain erosion by migration of channels and new floodplain formation by mud sedimentation and/or vegetation development.The effective cohesiveness may change an unconfined braided system into a dynamic self-confined meandering system or even a straight, laterally immobile channel without bars (Makaske et al., 2002;Kleinhans and Van den Berg, 2011).Here we inves-L.Braat et al.: Effects of mud supply on large-scale estuarine morphology tigate whether mud has similar effects on large-scale planforms that develop over centuries to millennia in estuaries.We especially need more knowledge about where mud deposits occur and how they influence the evolution of the estuary over long timescales.We first quantify mudflat properties in two Dutch estuaries and then review approaches to mud modelling.

Spatial pattern of mudflats in estuaries
In this study we use data from two Dutch estuaries, the Western Scheldt estuary and the Ems-Dollard estuary.The Western Scheldt is a mesotidal to macrotidal estuary with a semidiurnal tide and is located in the southwest of the Netherlands (Fig. 1f).The estuary has a tidal prism of 2×10 9 m 3 and maximum channel velocities are on the order of 1 − 1.5 m s −1 (Wang et al., 2002).The freshwater discharge is on average 120 m 3 s −1 from the Scheldt River.The Ems-Dollard is a mesotidal estuary with a semi-diurnal tide and is located at the most northern part of the border between Germany and the Netherlands (Fig. 1f).The estuary has a tidal prism of 1 × 10 9 m 3 and maximum channel velocities are on the order of 1 m s −1 (Dyer et al., 2000).Freshwater input comes from the Ems River with an average discharge of 80 m 3 s −1 .We use these estuaries because they are relatively well documented, although bed composition data are rather scarce compared to bed elevation scans.The disadvantage of data from a well-studied estuary is that anthropogenic influences are usually considerable, so we only look at the general patterns and properties of the mud.Here we combine independent measures of mud content in surficial sediment: (1) a bed sampling dataset of the Western Scheldt (Fig. 1a; McLaren, 1993McLaren, , 1994)), (2) probability of clay in the GeoTOP map (v1.3) of interpolated borehole data in the top 50 cm of the bed (TNO, 2016) (Fig. 1b and e) where clay is defined as more than 35 % lutum (< 2 µm) and less than 65 % silt (< 63 µm) (Vernes and Van Doorn, 2005), (3) yearly Western Scheldt ecotope maps of Rijkswaterstaat (2012), in particular the mud-rich areas above the low water level (Fig. 1c) that are based on aerial photographs, and (4) the sediment atlas of the Waddenzee (Rijkswaterstaat, 2009) drawn from bed sampling in 1989 (Van Heuvel, 1991), which includes the Ems-Dollard (Fig. 1d).
Data from the two estuaries indicate that mud deposits on the sides of the estuary that are then shielded from the strongest tidal flow (Fig. 1a-e).Large fractions of mud are also found on bars, which is in general agreement with the estuarine facies description of Dalrymple and Choi (2007).The hypsometric curves indicate that most of the mud is deposited on the intertidal areas (Fig. 1h and i), yet significant mud fractions are also found in channels.Additionally, larger mud fractions occur in the single-channel upper estuaries and cover a large part of the width of the estuary (Fig. 1a, d and  g).To summarise, 10-20 % of the lower estuary cross section is typically covered by mud with higher fractions up to about half the cross section in the single-channel upper estuary.

Past and novel modelling approaches for sand-mud mixtures
In past long-term morphological modelling of estuaries, sand and mud were always considered separately, partly because the interactions between sand and mud are complicated.Models used either sand (e.g.Van der Wegen et al., 2008) or sand and mud without interactive transport (e.g.Sanford, 2008).However, sand and mud interact, which affects the erodibility (see Van Ledden et al., 2004a, for review).Such interactions include dominant mud with some sand that behaves as mud, but for lower mud fractions there is mixed behaviour (Van Ledden et al., 2004a).In particular, mixed sediments increase erosion resistance and decrease erosion rates when the critical shear stress is exceeded compared to pure sand (e.g.Torfs, 1995;Mitchener and Torfs, 1996).This behaviour is highly sensitive to small amounts of mud, and the highest critical shear stresses for erosion occur with 30-50 wt % sand (e.g.Mitchener and Torfs, 1996).
Over the past decade, mixed sediments have been implemented in several modelling software packages (Van Ledden et al., 2004a;Waeles et al., 2007;Van Kessel et al., 2011;Le Hir et al., 2011;Dam et al., 2016).Long-term morphologic calculations are rare due to computer limitations and lack of spatially and temporally dense data of mud in the bed.For deltas, on the other hand, long-term morphologic development by numerical modelling (Edmonds and Slingerland, 2009;Caldwell and Edmonds, 2014;Burpee et al., 2015) showed large effects of mud on plan shapes, patterns and dynamics with fairly simplistic sediment transport processes.In particular, cohesion reduces the ability to re-erode, resulting in more stable bars and levees and longer and deeper channels.Physical experiments produced similar results for deltas (Hoyal and Sheets, 2009) and for river meandering (Van Dijk et al., 2013).However, the sensitivity of the numerical models to parameters such as erodibility and settling velocity indicate that the value of long-term modelling exercises with the current state of the art is to develop generalisations and trends rather than precise hindcasts and predictions of specific cases.
Past long-term morphological modelling studies of estuaries that did not include mud showed channel bar patterns that are similar to those in nature (Hibma et al., 2003;Van der Wegen and Roelvink, 2008;Van der Wegen et al., 2008;Dam et al., 2013).Cases in which boundaries eroded unhindered (Van der Wegen et al., 2008) developed towards a state of decreasing morphodynamic activity as size and depth continued to increase and morphodynamic equilibrium was not reached.Most models, however, including the few models with mud, assumed prescribed planform shapes with nonerodible boundaries (Lanzoni and Seminara, 2002;Hibma et al., 2003;Van der Wegen and Roelvink, 2008;Dam et al., 2013; Dam and Bliek, 2013) allowing equilibrium in some cases.However, to obtain a dynamic equilibrium of planform shape and dimensions in which bank erosion on average equals sedimentation, the formation of cohesive mudflats needs to be incorporated in models with erodible banks.Regardless of the fact that most natural estuaries are in disequilibrium as they continuously adapt to changing boundary conditions and anthropogenic influences, it is of interest to know whether these systems could develop a morphodynamic equilibrium and on which variables this depends most.
The objective of this research is to determine the effects of mud supply on equilibrium estuary shape and dynamics.This fills a gap in the literature by combining millenniumscale morphological modelling of estuaries and the effects of sand-mud interaction.We examine estuary formation from idealised initial conditions and a range of boundary condi-tions and run models for 2000 years in order to study tendencies towards dynamic equilibrium.We hypothesise that mud will settle into mudflats flanking the estuary that resist erosion and thus self-confine and narrow the estuary and reduce channel bar mobility and the braiding index.As a result we expect that self-formed estuaries develop a dynamic balance between bank erosion on the one hand and bar and mudflat sedimentation with resistant cohesive mud on the other hand.

Methods
The methodology was to set up an idealised scenario loosely inspired by the Dyfi, i.e.Dovey, estuary in Wales and to vary the most relevant boundary conditions.These include mud concentration supplied at the upstream boundary, mud supplied at the coastal boundary, surface waves, river discharge and tidal amplitude.There is a host of other initial conditions, boundary conditions and other variables that can be tested, such as other tidal components and other initial valley shapes.For example, the application of certain tidal components can lead to a change in the import or export tendencies of tidal systems (Moore et al., 2009), as can river inflow (Guo et al., 2016).However, our aim is to isolate the effects of mud, which requires the simplest possible conditions without non-linear interactions between imposed tidal components.Furthermore, we tentatively assume that the model is sufficiently sophisticated to reproduce the general behaviour found in nature of the phenomena under investigation, which will be discussed later.We chose the Dovey estuary as inspiration because direct human influences are relatively low compared to the Western Scheldt and the Ems-Dollard.Even though the system is still very natural, there is enough information about bathymetry and hydrodynamic data to develop the model and complementary model studies (Brown and Davies, 2010).Furthermore, it is one of the sandy estuaries in the UK that is included in the dataset of Prandle et al. (2005) that we will use later in the discussion.

Numerical model description
We used the modelling package Delft3D version 4.01.00, which is a process-based modelling system that consists of several integrated modules (Lesser et al., 2004).This modelling system is state of the art, open source and widely used and tested.It includes the possibility to use both sand and mud in the calculations.The depth-averaged version of Delft3D with parameterisation of spiral flow was used to keep the computational time for long-term simulations below 1 month.Furthermore, we excluded the effect of salinity and temperature on the hydrodynamics, as it was assumed that the effect of density differences would be limited in 2-DH and in well-mixed shallow estuaries.Auxiliary tests in 3-D with five layers and salinity confirm the assumption of well-mixed conditions.Furthermore, the estuary Richardson number (as defined by Fischer, 1972) is 0.036 and the Rouse number is < 0.01, further supporting the assumption of a well-mixed estuary for salinity and suspended sediment.The effects of the Coriolis force, organisms and wind are ignored for generalisation and simplicity.Hydrodynamics were calculated by solving the depth-averaged shallow water equations (Eqs.1-3): where η is water level with respect to datum (m), h is water depth (m), u is depth-averaged velocity in the x direction (m s −1 ), v is depth-averaged velocity in the y direction (m s −1 ), g is gravitational acceleration (m s −2 ), C is the Chezy friction parameter (m 0.5 s −1 ) and v w is the eddy viscosity (m 2 s −1 ).
The SWAN module was used to implement the effect of short waves.We used two-way coupling between the flow and wave module with an interval of 3 h.At four stages during every tidal cycle, SWAN calculated the wave conditions from the current situation of the morphological model.The waves enhanced turbulence and bed shear stress by wavedriven currents in the morphological model.The sediment transport was only affected by the enhanced bed shear stress by the wave-current interaction and not by enhanced turbulence.
A recently developed module for mixed sediments incorporates the effect of bed composition on erosional behaviour and hence morphology (Van Kessel et al., 2011, 2012).This module is a partial implementation of Van Ledden (2001) and Jacobs et al. (2011) and tracks spatial and temporal bed composition for multiple grain sizes of sand and mud with erosional characteristics depending on bed composition.In this paper we only used one sand fraction and one mud fraction (Table 1) and applied a uniform roughness.
Cohesive sediment, i.e. mud, is defined as the mixture of the clay (< 2 µm) and silt (2-63 µm) fractions with cohesive behaviour caused mainly by physico-chemical forces between the clay particles.This cohesive behaviour causes complex processes that influence the erosion and deposition of sediments.In the model we distinguish two erosion modes.Above a critical mud content (p m,cr ) of the bed, cohesive particles cover sand particles so they are not in direct contact, which limits erosion for both sand and mud (Torfs, 1995(Torfs, , 1996)).Below this critical mud content, friction and gravity oppose sediment transport for sand.The critical mud content was chosen to be at a mass fraction of 0.4, which depends on site-specific silt-clay ratios because only the clay fraction is cohesive (McAnally and Mehta, 2001;Van Ledden et al., 2004a).This value is higher than found in flume experiments (0.1-0.2, Torfs 1995;0.05-0.15, Torfs 1996;0.02-0.15, Mitchener and Torfs 1996) but was based on the silt-clay ratios of Dutch tidal systems (0.25-0.5;Van Ledden et al., 2004b).
When the bed is defined as non-cohesive (p m < p m,cr ), a traditional sand transport equation was used.Here we chose the Engelund and Hansen transport equation (1967;Eq. 4): where q s is sediment transport (m 3 m −1 s −1 ), U is the magnitude of the flow velocity (m s −1 ), is the relative density (ρ s − ρ w )/ρ w and D 50 is the median grain size (m).This equation does not distinguish between suspended and bedload transport but considers total transport.The Partheniades-Krone formulation was used to calculate the erosion rate of mud (Partheniades, 1965, Eq. 5): where E m is the erosion flux of mud (kg m −2 s −1 ), M is the erosion parameter (kgm −2 s −1 ), S is the erosion or depositional step function, τ cr,e is critical shear stress for erosion (N m −2 ) and τ cw is the maximum bed shear stress due to currents and waves (N m −2 ).When the bed is cohesive (p m > p m,cr ), the mud and sand fluxes are proportional to the mud and sand fraction.The erosion rate of mud is calculated by the Partheniades-Krone formulation (Partheniades, 1965;Eq. 5) similar to the noncohesive regime.The erosion rate for sand, on the other hand, was based on the entrainment of mud because sand particles are included in the cohesive matrix (Eq.6).In this way sand can only be eroded when mud is eroded.Bedload transport was assumed to be zero in the cohesive regime.
The advection-diffusion equation further describes the suspended sediment following from the Partheniades-Krone formulation.Sand and mud behave independently in suspension and segregation will occur with low concentrations (Torfs, 1996).For simplicity we assumed a constant settling velocity of 0.25 mm s −1 for mud, ignoring the fact that settling velocity depends on flocculation influenced by concentration, residence time, salinity, pH, turbulence and biochemical effects (e.g.Mietta et al., 2009).The settling velocity is typical for fluvial mud (0.1-0.4 mm s −1 , Temmerman et al., 2003), which we supply in the default run, and is relatively low for marine mud.Deposition of mud is determined by the concentration, settling velocity and the step function similar to Eq. (5).However, many studies show that deposition is continuous, and a threshold for deposition is therefore absent (short review in Sanford, 2008).This is approached in the model by setting a very high critical shear stress for sedimentation (Table 1).
Divergence of sediment fluxes for bedload and the erosion-deposition difference for suspended sediment cause bed level changes.To track the mud and sand fractions in the bed, a bed module was used with a mixed Eulerian-Lagrangian approach (Van Kessel et al., 2011, 2012) similar to Le Hir et al. (2011) and Sanford (2008).An active Lagrangian layer of 10 cm was used in which sediment exchange occurs with the water column.This active layer had a constant thickness and moved through the vertical framework with bed aggradation and erosion.Below the active layer we used several vertically fixed Eulerian layers to store bed composition in the vertical (Table 2).The advantage of Eulerian bed layers is that artificial mixing by vertically moving layers is prevented.The advantage of a Lagrangian active layer is that the thickness is constant, which is desired because strong bed armouring is prevented and the thickness affects the timescales of the system (Van Kessel et al., 2012).
To speed up morphodynamic calculations, the bed level change in each time step calculated from the divergence of sediment fluxes and the erosion-deposition difference for suspended sediment was multiplied with a morphological factor of 400 (Table 2).Extensive studies showed reasonable results up to a morphological factor of 1000, though it is recommended not to go above 400 (Roelvink, 2006;Van der Wegen and Roelvink, 2008, Fig. A13).Using a morphological factor is an efficient way of speeding up long-term morphodynamic calculations that is widely used (Roelvink, , 2008;Le Hir et al., 2015;Dam et al., 2016).
When the water level changes during a tidal cycle, flooding and drying of the intertidal area occurs.To prevent complicated and time-consuming hydrodynamic calculations with very small water depths, a threshold is set for drying and flooding (Table 2).When the water depth is below this threshold the velocity is set to zero.Since the velocity in dry cells is zero, there is no sediment transport in dry cells, even when considerable erosion occurs in a wet cell next to it.Therefore, dry beach and bank erosion was implemented to drive lateral bed lowering.A user-defined factor (Table 2) determines the fraction of the erosion flux that is assigned to the adjacent dry cells.
The transverse bed slope effect is a very important parameter for bar dimensions and behaviour in morphological models that is often used as a calibration parameter (Schuurman et al., 2013).In estuary models the transverse bed slope effect is often set to be much stronger than the advised default settings to prevent unrealistically steep banks and narrow bars and channels from forming (Van der Wegen and Roelvink, 2012).The reason for this is unclear, but unravelling this is beyond the scope of the present paper so we use settings similar to earlier studies (Van der Wegen and Roelvink, 2012).We used the transverse bed slope predictor of Koch and Flokstra (1980) as extended by Talmon et al. (1995): where θ is the shields parameter and α = 0.2, which is much lower than the default of 1.5 for rivers, and β = 0.5.We chose the Engelund-Hansen transport formulation because other relations, in particular Van Rijn (1993), Van Rijn et al. (2004), andVan Rijn (2007), resulted in higher bars and much deeper and straight channels with sudden (up to 90 • ) sharp bends, which would require transverse bed slope parameters that differ by 2 orders of magnitude from the theoretical value in estuarine settings (Van der Wegen and Roelvink, 2012).Furthermore, changing bed slope parameters does not fix the channel pattern issues.For long-term morphological modelling, Engelund-Hansen produces more realistic morphologies.The disadvantage of our method is that the present code for sand-mud interaction with Engelund-Hansen does not yet incorporate a gradual transition in critical shear stress for erosion between the cohesive and non-cohesive regime.Additionally, mud would ideally erode proportionally with sand in the non-cohesive regime as sand erodes with mud in the cohesive regime, but this is not yet implemented for Engelund-Hansen and is therefore also not described in this method section.These issues are beyond the scope of the present paper and require further research and model code development.

Model schematisation
The modelled domain is 30 by 15 km of which 10 by 15 km is sea area (Fig. 2).The grid is rectilinear with a resolution that varies between 50 by 80 and 125 by 230 m.Cell size increases from the initial estuary shape to the sides and offshore to increase resolution in regions of interest and to decrease computation time.The initial bathymetry is in the shape of an idealised funnel-shaped estuary.This exponential shape was also found in previous modelling research (Lanzoni and Seminara, 2002;Canestrelli et al., 2008;Lanzoni and D'Alpaos, 2015) and obtained from field data (Savenije, 2015).The estuary is 3 km wide at the mouth and decreases exponentially to a channel of 300 m width over 20 km.The bed level linearly increases in elevation from −2 at the mouth to 2 m at the upstream boundary and 2 to 3 m on dry land (Van der Wegen and Roelvink, 2008).The sea has a maximum depth of 15 m.Van der Wegen and Roelvink (2008) argued that initial bathymetry does not greatly affect the dynamic equilibrium shape because dry cell or bank erosion is allowed in the model and the model will therefore develop a self-formed (alluvial) estuary shape.However, initial shape affects the time needed to form the equilibrium planform shape and the size of the ebb delta in the absence of waves and littoral transport, which is the default situation in our idealised estuary.We therefore started with a funnel shape to save calculation time and decrease the size of the ebb tidal delta.The shape is given as where W mouth = 3000 m is the width of the estuary at the mouth , L b = 3362.6m is the e-folding distance over which the width of an exponential channel is reduced by a factor of e and x is distance from the mouth (m).The shapes of modelled estuaries are characterised by the funnel-shape parameter (Davies and Woodroffe, 2010) calculated as efolding length normalised by mouth width at that point in time (Eq.9).Lower values of the characteristic funnel length indicate stronger funnelling in the sense of more rapid narrowing from the mouth in the landward direction.In this way estuary shape is normalised by estuary size.Three open boundaries are used: two cross-shore water level boundaries and one upstream discharge boundary.At the water level boundaries an M2 tide is prescribed with a default tidal amplitude of 1.5 m and a phase difference of 3 • (6 min between the two cross-shore boundaries over 15 m), resulting in alongshore tidal wave propagation as for the Dovey estuary.The western boundary is closed because three open sea boundaries created instabilities in the corners of the model.The chosen tide is exactly cross-shore and therefore the closed boundary does not affect the tide.With these simple conditions tidal asymmetry in the estuary is entirely caused by basin geometry and river flow and not by prescribed overtides.For generalisation purposes and simplicity we exclude the known effects of imposed multiple tidal constituents on residual transport and morphology (Guo et al., 2016).Discharge is prescribed as a constant value through time of 100 m 3 s −1 over the inflow cross section.However, the partitioning of discharge between the upstream grid cells of the river at the boundary varies sinusoidally through time from one side to the other to simulate weak upstream "meandering" perturbations with a period of 500 years to trigger bars if the self-formed channel aspect ratio would allow bars (Schuurman et al., 2016).
In some model scenarios waves were imposed at all sea boundaries including the closed, offshore boundary parallel to the coast.Waves have a 6 s peak period and a significant wave height of 0.7 m.This is the highest possible significant wave height for which no coastal erosion occurs.The effects of the added waves keep mud in suspension in the sea area because of the enhanced bed shear stress (Eq.5), which prevents the formation of a muddy coastline and meanwhile causes no significant sand transport outside the estuary.The addition of waves prevents instabilities at the boundaries that were due to the deposition of marine mud (Fig. A1 in the Appendix).Due to the choice for Engelund-Hansen as sediment transport formulation, sand stirring is excluded.Only indirect sand transport effects occur because the enhanced bed shear stress influences the currents.
The initial bed composition in the entire domain is 100 % sand.In some scenarios mud was supplied to the estuary at the discharge boundary and/or the sea level boundaries.Mud was supplied as a constant concentration, which means that the mass of mud transported into the model depends on the hydrodynamic conditions.For sand supply we used equilibrium conditions at the boundaries, meaning that the capacity of the flow to transport sand (Eq.4) at the boundary determined the sand supply rate, which prevents erosion and deposition at the boundaries.

Scenarios
The model was run for 23 scenarios of boundary conditions on the same initial conditions.One run constituted 5 hydrodynamics years and 2000 morphological years and took about 20 days in real time on one desktop core.Multiple scenarios were computed at the same time, so parallel computing was not necessary.The runs with waves took much longer and were therefore terminated at 1250 years.
We varied fluvial mud input concentration to assess the primary effect on the shape and size of estuaries.The effects of the source of mud were tested by comparing scenarios with fluvial input, marine input, both marine and fluvial input at the same time or no mud input.We further examined the effects of waves, river discharge and tidal range.To assess the sensitivity of the model we varied uncertain numerical and process-related parameters: critical mud content for cohesive behaviour, active layer thickness and settling velocity (Table 3).The model scenarios were analysed by studying the bathymetric changes, mud deposits and geometry of the final bathymetry.These results are compared to each other.In the discussion we compare model results to the data of natural estuaries presented above.
About 100 pilot models led us to select the model settings and boundary conditions presented in this paper.For example, we evaluated different initial bathymetries to test their effect on time to equilibrium.We found that the model could both erode and fill the initial basins for otherwise equal conditions, meaning that the initial shape is only of limited influence on final equilibrium.Moreover, we found that an exponential shape close to the equilibrium size saves considerable computation time and reduces the size of the ebb delta, which then, in turn, has a smaller effect on the incoming tide.Pilot runs with alternative sediment transport formulations confirmed the findings of Van der Wegen et al. (2008) and led to the choice for the Engelund-Hansen transport equation and a transverse bed slope parameter of 0.2 (Table 2).Furthermore, pilot runs showed that initial random bed perturbation was unnecessary to trigger bar development.Finally, to test the assumption that the estuaries are well mixed, the default run was restarted in 3-D after 1200 years with salinity and five sigma layers.These results indicated well-mixed condi-L.Braat et al.: Effects of mud supply on large-scale estuarine morphology tions and some influence on sediment transport but limited influence on large-scale morphology.

Results
Here we first describe the general development towards equilibrium of the default scenario with fluvial-fed mudflats.Secondly, we study the hydrodynamics and sediment transport in more detail for the equilibrium condition of the default scenario.Then we describe and compare trends in all scenarios, focussing on mud supply, mud source, the effects of waves, river discharge and tidal amplitude.Finally the mud parameter sensitivity runs are presented.Figures with detailed results for scenarios with hydrodynamic variables are shown in the Appendix.

General development
The final morphology of the default scenario (run 01) after 2000 years with a fluvial mud supply concentration of 20 mg L −1 is a self-confining bar-built estuary flanked by mudflats with migrating channels and bars (Fig. 3d and i).The width of the final morphology decreases exponentially in the upstream direction, similar to the initial condition but with self-formed, freely erodible banks.
During the first stage of the development the mud enters the system by river discharge, which is rapidly distributed over the whole estuary within the first few years.The upstream part narrows immediately, while narrowing at the mouth starts after about 150 years and continues for roughly 700 years (Fig. 3b and c).After 200 years the sand within the estuary is redistributed and the ebb delta starts to form.The ebb delta continues to prograde as sand and mud are supplied constantly, whilst coastal sediment transport is absent.Since we are not interested in the evolution of the ebb tidal delta and littoral processes are not well modelled in this set-up, the area downstream of the coastline is excluded in further analyses.
Within the first 3 years the upstream part of the estuary starts meandering and the downstream part starts braiding.Meanders grow and migrate downstream, while bifurcations develop and chute cut-offs occur.Within 200 years, an initial bar pattern has developed throughout the estuary, and the channel pattern is characterised by mutually evasive ebb-and flood-dominated channels (Fig. 3b and c).The bars continue to migrate downstream throughout the simulation as an effect of the fluvial discharge and sediment input.After about 1000 years the bar channel pattern appears to have reached a dynamic equilibrium with channels approximately 4 m deep, bars elevated to the mean water level and a linear sloping bed level (Fig. 8a-d).
Morphodynamic equilibrium in which average bank erosion equals sedimentation is indicated by a net bed level change fluctuation around zero (Fig. 3j).This means that there is no net accumulation or erosion in the estuary, so there is no net import or export of sediment.Furthermore, we observe that the absolute bed level changes approach a constant value (Fig. 3e).Net bed level change is determined by summation of the elevation change of each cell multiplied by the area of the cell, while the absolute bed level change uses the absolute value of the change in elevation.The initial changes in which the estuary adapts to the boundary conditions (like width and depth adaptation) happen within centuries, while the dynamic behaviour of bars and channels continues throughout the simulation, as also shown by the constant non-zero value approached in Fig. 3e.If the mean of the absolute bed level changes approached zero, then the bathymetry would become fixed.It could be argued that the lowering in Fig. 3e indicates that a true equilibrium was not reached and will not be reached because the river continues to import sand and the ebb delta continues to grow in the near absence of littoral processes.However, for our purposes and timescale of interest, Fig. 3j indicates that the equilibrium planform geometry of the channels, bars and mudflats in the estuary was reached after about 500-1000 years.

Hydrodynamics and sediment transport
Tidal water levels and velocity vary along the estuary: at the seaward boundary the tide is a symmetrical M2, while further into the estuary the tide becomes asymmetrical (Fig. 4).At the mouth the water level rapidly increases from low to high water and slowly decreases from high to low water.There is no phase lag anywhere along the estuary between water level and velocity since slack water occurs exactly at high and low water.The tidal range decreases further into the estuary mainly by a decrease in the low water amplitude (Fig. 4a).Likewise, flood velocities are reduced, while ebb velocities remain about constant (Fig. 4b, e and f).Additionally, the duration of the ebb flow is longer than the flood (Fig. 4g), which is similar to the current Dovey estuary.In the Dovey, the flood phase takes 5 h and the ebb phase 7 h near Aberdyfi (Brown and Davies, 2009), which is similar to the green line in Fig. 4b.The ebb flow increases in relative velocity and duration further upstream because the contribution of the river increases in this direction.The tidal excursion length is a little over 6 km.The location of the 1 ppt isohaline is 5 km upstream of the estuary mouth during high tide, which was inferred from the 3-D restart of the default scenario with salinity.In the Dovey estuary the 1 ppt isohaline is around 7.5 km and also well mixed (Baas et al., 2008).We consider this to be in good agreement because the model discharge is larger and the spit is ignored.
These hydrodynamics might suggest that the estuary is still an exporting system and not in equilibrium.However, spatial variation is very important.We observed flood-dominant velocity amplitudes over shallow areas, like bars and mudflats, and ebb-dominant velocity amplitudes in the channels (Fig. 4h) so that flood discharge and ebb discharge balance except for the net river inflow.In more detail, flood-dominant No Larger settling velocity for mud velocity amplitudes occur especially above mudflats (Fig. 5) which typically occur at higher elevations (Fig. 7).We observe that the lower estuary evolves from a system with very strong ebb-dominant peak velocities to a system which is only slightly ebb dominant (Fig. 5, black squares).Both high and low areas show this trend, but high areas are typically less ebb dominant.When more mudflats build up in the system these areas change from ebb-to flood-dominant peak velocities over time.SPM (suspended particulate matter) levels reach the highest local concentrations of 45 mg L −1 between 2 and 4 h after high tide with a typical mean concentration similar to the input concentration of 20 mg L −1 .The mean SPM levels of the Dovey are comparable and estimated to be 32 mg L −1 (Painting et al., 2007).The typical non-dimensional shear stress (Shields number) of the model is 0.27.Over the whole model run, 4000 m 3 of sediment is imported into the estuary of which 7800 m 3 of mud is imported and 3800 m 3 of sand is exported.
In the final stage of the model, net sediment transport is in the ebb direction for bedload transport and suspended transport, i.e. of sand and mud (Fig. 4d).Notably, the amount that is transported through the mouth (solid line) is equal to the sediment input from the river (dotted line).This shows that there is no net deposition or erosion in the estuary in agreement with Fig. 3d, meaning that the estuary is in equilibrium.
The river discharge is about 7.5 % of the maximum tidal flood discharge and contributes to the ebb flow.Therefore the volume of water flowing through the mouth during ebb is always slightly higher than in the flood direction (Fig. 4c).Because the tidal prism is calculated as the product of veloc-ity, duration and cross-sectional area and because of the ebb dominance through duration and velocity amplitude asymmetry, the cross-sectional flow area is larger for the flood flow; otherwise there would be a large imbalance in tidal prism going in and out of the estuary.

Effects of mud supply
We will now compare other scenarios (runs 03, 10 and 09) with the default run (01).In most scenarios mud accumulates on the flanks of the estuary where the velocities are low and in the upper estuary where it covers a relatively large fraction of the width (Fig. 6f, g and h).Locally, mud accretes on bars that are rather stable (e.g.Fig. 6g and h; on the ebb delta).The initiation of mudflats proceeds through the positive feedback identified in the model description: once mud starts settling somewhere, the mud fraction in the bed rapidly increases beyond the critical mud fraction for muddominated behaviour.As a consequence, the critical shear stress for sand erosion equals the entrainment threshold of mud (Eq.6).The mud-dominated mixed sediment thus becomes more difficult to erode and more rapid aggradation of mud is likely to occur.
Migration rates of channels decrease considerably due to the addition of cohesive material (Fig. 9a-h).Bar splitting and merging related to chute cut-offs and avulsion are also reduced with increasing mud concentrations.In Fig. 9a-d channels move through a cross section at the mouth through time.The experiments with a larger supply show slower and less movement of the channels.For example, a large bar forms in the mouth after about 1100 years in the scenario with only sand (Fig. 9a) and in the scenario with a mud sup-.ply of 50 mg L −1 (Fig. 9d).In the run with mud, the bar is covered with mud and becomes fixed, while the large bar in the scenario with only sand migrates about 1 km.Absolute bed level changes also indicate that mud input decreases dynamics (Fig. 9y-II) because there is less bed level change per time step.
The mudflats have a strong effect on the final shape of the estuary (Fig. 6).Firstly, an increase in fluvial mud input concentration leads to stronger self-confinement of the estuary.By depositing mud on the sides of the estuary, the banks become more stable and limit (further) erosion due to increased critical shear stress.Self-confinement of estuaries is clearly observed when the models with mud supply are compared to the control run without mud (Fig. 6a).The runs with mud are narrower and have a smaller surface area due to filling of the initial bathymetry, while the sand run has expanded in size.Consequently, the braiding index lowers with increasing mud concentration (Fig. 8e-h).In contrast, estuarine surface area continues to increase over time for the control run with only sand (Fig. 9q).After the initial rapid change, the increase in area and width is linear, driven by dynamic channels and bars and is unhindered by bank stability.This suggests that there is no equilibrium shape under these conditions as is also reflected in the absolute and net bed level change (Fig. 9y and III).The absolute bed level change does not approach a constant value and the net bed level remains negative, demonstrating the sand-only estuary to be a continuously exporting system.
For estuaries with fluvial mud, higher concentrations lead to narrower (Figs.8i-l and 9i-l) and smaller (Figs.9q-t and 7) estuaries.Moreover, in some places the width of the estuaries with mud supply is narrower than the initial width, supporting the finding that the initial bathymetry is of limited influence because the system is able to fill and to expand (see methods).Furthermore, tidal bars become higher with increasing mud concentrations, which results in an increased average bed level (Fig. 8a-d).Mud is deposited almost nowhere in the channels and therefore does not limit bed erosion by cohesion (Fig. 7).As a result we infer that the shallower channels in increasingly muddy estuaries mainly result from the decrease in estuary width and concurrent reduction of intertidal area, tidal range and tidal currents (Fig. 8).
With larger mud concentrations, a larger area of the estuary is covered with mudflats (Fig. 8m-p).The mud cover maps (Fig. 6e-h) indicate that although the distribution of the mud is quite similar for different concentrations, the overall mud cover over the estuary length increases with mud input concentration (Fig. 8m-p).In general, more mud leads to wider mudflats on the side and seems more likely to deposit mud on mid-channel bars.The maximum fraction of inter-tidal area shifts from the middle estuary to the lower estuary for increasing mud concentration.At the same time mud increasingly deposits on lower elevations as seen in the strong increase in cumulative area just above the mean water level (Fig. 7).
The estuaries with mud are shorter than the estuary with only sand.The length of the estuary is defined as the distance between the mouth and the limit of tidal influence (where tidal range reaches zero in Fig. 8q-t).Estuaries are shorter for scenarios with mud compared to the scenario with only sand, but mud concentration seems of limited impact.Mud supply concentration has some effect on the funnelling of the estuary, although the temporal variation in the funnel shape is less for higher concentrations (Fig. 9m-p).In general, funnel-shape strength first decreases and then increases again.This has to do with the order of widening and narrowing of different parts of the estuary.The width of the mouth always decreases at the start of the run, but after about 400 years that can change into widening or continue (Fig. 9i-l).This narrowing at the start increases the funnel parameter.Upstream, the estuary width initially increases, which also decreases funnelling.The e-folding length of the scenario without mud supply is not the shortest compared to the other scenarios, but the run without mud results in a larger estuary.For a longer convergence length, the funnelling parameter can be the same if the estuarine mouth is bigger.

Difference between fluvial and marine mud supply
Estuaries develop very differently when mud is imported from the sea rather than from the river under the assumption that the mud characteristics are the same (runs 01, 02 and 04; Figs.A1-A3).This scenario could, for example, occur when upstream mud supply is obstructed by the construction of a dam, but it is of higher importance in our understanding of sediment provenance.For marine mud, the mudflats form only in the lower estuary up to 9 km upstream from the mouth because mud can only occur in regions where there is significant flood flow to transport the mud upstream.For fluvial mud, on the other hand, mudflats form along the entire length of the estuary.Mud supply from both boundaries simply has a combined effect with mud distributed along the whole estuary and the highest mud cover near the mouth.Estuaries are narrower with fluvial mud supply compared to the marine mud supply and the sand-only control run.In the case of marine mud supply, the estuary decreases in width near the mouth, but upstream width and bed level are similar to the estuary without cohesive sediment.In the first 500 years the width at 1 km from the mouth decreases and is partly taken in by mudflats but returns to the initial width after 2000 years.On the other hand, estuary width increases at 4 and 8 km from the mouth.For the scenario with both mud from the sea and the river, the estuary mouth is narrower than for only marine or fluvial mud.
The total estuary area continues to increase for the scenario with marine mud supply because the upper estuary widens similar to the run without mud.Likewise, the estuary does not confine itself by cohesion and therefore does not reach equilibrium.The estuaries that include fluvial mud supply eventually reach a constant area over time and do reach equilibrium.
The estuary with fluvial mud supply shows stronger funnelling due to more narrowing between 5 and 10 km than between 0 and 5 km from the mouth.The estuary with marine mud supply shows an opposite trend with stronger narrowing near the mouth, leading to a lower convergence.The length of the estuary is shorter for scenarios that include fluvial mud supply, and the tidal range and flood velocity along the estuary decrease faster.
Not all mud settles in the estuary, but a lot is transported out of the estuary or never enters the estuary from the seaward boundary.Part of this mud is deposited at the coastline and part is transported out of the model domain.Because mud is supplied as a concentration depending on the discharge, a much larger volume of mud is supplied to the system when the sea supplies mud.This large volume of mud causes significant deposition at the coast and affects morphology at the mouth.We consider this an artefact due to the lack of littoral processes.

Effects of hydrological boundary conditions: river, tide and waves
Changes in the boundary conditions in the form of tidal amplitude and discharge did not seem to alter the location of the mudflats but only the size (runs 01, 22, 07, 08, 21, 20, 05 and 06; Figs.A4-A9).More and larger mudflats formed with higher discharges.An optimum in mudflat size occurred for increasing tidal amplitude.With lower amplitudes there is less intertidal area and therefore less space for mudflats, and with higher amplitudes the higher velocities prevent deposition.There is a balance between the tidal flow and fluvial flow into the estuary.When the river discharge becomes larger, tidal damping occurs under the influence of increased river discharge by friction (Horrevoets et al., 2004).Therefore, the limit of tidal influence is further downstream, decreasing the tidal prism and therefore tidal velocity.This means that the excess width can be filled until the appropriate width-depth ratio of the river to this point.When the river has less influence, the tidal intrusion is larger with higher velocities.This balance influences the morphology: relatively stronger tidal influences lead to larger estuaries when the more river-dominated estuaries decrease in size, fill up and eventually evolve into deltas.More specifically, no river discharge leads to large tidal meandering channels in the lower estuary with a filled upper estuary.On the other hand, larger discharges lead to a transi- tion from filled estuaries to a delta.In addition to a change in the river-tide balance, increased discharge means more sediment input at the equilibrium boundary condition used for sand.Mud is also supplied as a concentration, and mud volumes therefore increase with higher discharges.As a result the system rapidly expands the ebb tidal delta, fills the estuary and transforms into a delta for the highest discharges.This means that the balance between fluvial discharge and sediment supply and the tide and tidal sediment export is changed.
The estuary shape scales with discharge, but size does not.Lower discharge leads to stronger funnelling of the estuary.On the other hand, size hardly changes with discharge despite the fact that larger discharges result in more vigorous channel migration and faster dynamics.We only observe a sudden transition in size from estuary to delta between a discharge of 100 and 150 m 3 s −1 .Adversely, tidal amplitude has a strong effect on the size of the estuary.In fact, a tidal amplitude of less than 2 m leads to closure of the estuary and formation of a muddy delta.The larger flow velocities with higher tidal range keep mud in suspension so that less mud settles in the estuary, in turn leading to less self-confinement.Systems with lower tidal amplitudes are therefore more likely to develop deltas rather than equilibrium estuaries.Further it is observed that larger tidal ranges lead to larger tidal meanders and bigger channels.An additional effect is that higher velocities due to increased tidal amplitude cause enhanced shifting of the channels, which prevents the settling of mud on the bars sufficiently to change the erosional behaviour and prevent the positive feedback of mud from having an influence.This effect is also caused by waves.
Waves (runs 29, 27, 28 and 25; Figs.A10-A12) prevent mud deposition at the coastline, prevent instabilities in the sea area and cause widening of the mouth.This especially leads to a limited influence of marine mud supply, though it is supplied 5 km further upstream with waves.For example, the run with marine mud supply and waves is very similar to the run without mud supply with waves.Because of the widening of the mouth by waves, tidal range, water levels and flow velocities increase, especially flood velocities.Additionally, widening at the mouth leads to a very strong funnel shape.Due to the waves there is generally little mud cover in the www.earth-surf-dynam.net/5/617/2017/Earth Surf.Dynam., 5, 617-652, 2017 lower part of the estuary.In nature, waves form spits and may even largely close off estuaries, but this does not occur in the model because the effects of short waves on the sediment dynamics is limited to the stirring of sediment.The results therefore strictly apply only to estuaries with limited wave influence and to inner estuaries more generally where wave action is limited.

Effects of sediment transport parameters
The sensitivity to active layer thickness (run 13 and 14), assessed by doubling and halving the active layer, did not lead to different large-scale trends in mudflat formation and estuary shape and dimensions.A different active layer thickness leads to a different pattern, but the large-scale characteristics of the pattern are the same.Likewise, the critical mud fraction (run 11 and 12) that determines cohesive and noncohesive behaviour had no significant effect on large-scale morphology.Initially there are slightly more dynamics in the run with the higher critical mud fraction, but this effect can be disregarded after some time.On the other hand, the order of magnitude of the settling velocity (run 15 and 16) had a considerable effect: a 10 times slower settling velocity re-sulted in an estuary with more similar geometry to the run with only sand, while 10 times higher settling velocities developed a delta due to larger sedimentation rates.This means that similar trends can probably be found for lower settling velocities with higher mud concentrations or by the addition of biotic effects on apparent cohesion.Furthermore, an increased tidal range with higher settling velocities might show similar results to the current settling velocity and tidal range.We predict that changing mud characteristics, such as settling velocity, erosion parameter and critical shear stress for erosion, would not affect the general trends and conclusions but would lead to slightly different equilibria.
We did not test the combined effect of changing the proportions of clay and silt, whereby the settling velocity and critical shear stress for erosion would probably be inversely correlated and have opposite effects, reducing the effects of these parameters.Additionally, we ignore consolidation, which especially affects the layer thickness and erosion characteristics of mud layers.With this in mind, we expect that the migration of deep channels eroding deep, old mud layers is overestimated.Additionally, we assume that the time in which thick mudflats develop is also overestimated and the critical shear stress of very recently deposited mud in reality is also overestimated due to the fluff characteristics of mud when it is still submerged.To summarise, we expect that most uncertainties are related to timescale, but we do not expect large differences in the general pattern and trends.

Discussion
The most important findings from the results are summarised in Fig. 10.Mud supply leads to self-confinement (Fig. 10d, blue) of the estuary by the development of mudflats on the sides (Fig. 10d, brown).We observed that larger mud supply concentrations leads to narrowing and filling of the estuary towards a dynamic equilibrium, while the estuary without mud supply continued to widen and grow in size (Fig. 10d and g, blue).Furthermore, we observe that mud raises the bed level, decreases the length, increases mudflat size, decreases dynamics and increases funnelling (Fig. 10a).Marine mud supply causes the development of a muddy coast and in this model only influences the mouth of the estuary.However, these effects might be overemphasised due to uncertainties in wave transport and chosen settling velocities.Narrowing of the mouth strongly decreases the funnelling of the estuary but is of little influence when waves are included.In scenarios with larger fluvial mud supply, larger flow discharge plus fluvial mud supply and lower tidal amplitude sediment filled the initial estuary shape and a delta developed (Fig. 10).By this we mean that the deltaic channels had only negligible tidal flow and were much smaller than the initial estuary.These results suggest a rather sharp transition from a narrow equilibrium estuary with significant tidal action to an extending river-dominated delta.10), 20 (default, 01) and 50 mg L −1 (09).(a-d) Minimum, mean and maximum bed elevation, high and low water level and minimum and maximum initial bed level, (e-h) braiding index and (i-l) estuary width defined as the initial width, maximum reach over the whole scenario run, the width of wet cells in the model, width defined by a threshold value that is used to mask the cells that are around the dry-wet cell threshold.(m-p) Intertidal area and mud cover as a percentage of the total area, (q-t) tidal range and (u-x) peak ebb and flood velocities.

Comparison to real estuaries
Model conditions fall within the parameter space of natural estuaries (Fig. 12; Table 4; Prandle et al., 2005;Leuven et al., 2016).The model has typically larger discharges than the small UK estuaries, but discharge and tidal amplitude falls well within the range of estuaries worldwide.
Several aspects of the bar patterns are further indications that the numerical models reproduce important emergent phenomena of real estuaries.For example, we observe ebband flood-dominated channels that are unique for tidal systems (Van Veen, 1950;Ahnert, 1960).Typical bar dimensions obtained from the models are in good agreement with natural estuaries from a large dataset (Leuven et al., 2016); for example, tidal bar length is approximately 7 times the partitioned bar width (maximum bar width divided by barb channels).Furthermore, bar length approximates the local width of the estuary.Bars without mud are generally longer and wider for this model study and relative to the local estu-ary width.The bars in the models are also slightly bigger with marine mud supply rather than for fluvial mud supply.The braiding index is strongly related to estuary width as found for natural estuaries (Leuven et al., 2016) and in agreement with the relation between tendencies to form floodplains in rivers and the resulting relation between channel aspect ratio and bar pattern (Kleinhans, 2010;Kleinhans and Van den Berg, 2011;Schuurman et al., 2016).
The completed model runs show that mudflat characteristics and behaviour are broadly comparable to natural estuaries.Spatial trends in the field data, shown earlier (Fig. 1), generally agree well with the model results.We observe similar depositional areas of mud on the sides of the estuaries in the form of mudflats (Figs.1a-e and 6e-h).In the centre of the lower estuary there is little mud compared to the mudflats on the sides.However, some mud is observed on some of the bars in the Western Scheldt (e.g.Fig. 1c) as in some model scenarios (Fig. 6h).Comparison of the observed and modelled hypsometries (Figs. 7 and 1g and h) shows that mud www.earth-surf-dynam.net/5/617/2017/Earth Surf.Dynam., 5, 617-652, 2017 is deposited at comparable elevations, mostly at intertidal areas and more specifically around the mean water level.We observe a strong increase in mudflats with the strongest increase is cumulative area.
The fluvial mud scenarios have relatively large fractions of width covered by mudflats in the upper estuary as in the single-channel upper estuaries in the Netherlands.Indeed, most mud is deposited in the middle and upper estuary where the estuary consists of only one channel.This is also clearly observed in the McLaren (1993McLaren ( , 1994) ) dataset of the Western Scheldt (Fig. 1a).The tidal river contained more mudflats than the lower estuary (Fig. 1f).Note that Fig. 8 underestimates the modelled mudflat surface shown in Fig. 6 because many cells are inactive in the computation because they increased in elevation.
Typically in the model, marine mud does not settle much or far in the estuary.This is not what is observed in the Western Scheldt.Verlaan (2000) studied the marine versus fluvial distribution of mud through the estuary.He found a sharp increase in mud fraction in the bed between Lillo and Saeftinghe from 10 to 75 %, which is far upstream in the narrow single-channel system.This might be a consequence of the assumption that settling velocities for fluvial and marine mud are the same, while the settling velocities of marine mud are typically significantly higher due to flocculation.Marine macrofloc settling rates might be as high as a few mm s −1   (Leuven et al., 2016).(Mietta et al., 2009;Leussen, 2011).It is also a likely possibility that the Western Scheldt is not comparable to our modelled system considering marine mud deposits because the salinity intrusion of the Dovey and Western Scheldt is incomparable.Mud deposition data from the Dovey estuary are unavailable although mudflats and muddy marshes are easily observable on aerial imagery (Leuven et al., 2017).
In the model we observe sharp transitions between areas without mud in the bed (< 10 %) and areas with very high mud fractions (70−100 %).This is also observed in the Western Scheldt according to Van Ledden (2002).More gradual transitions of mud are expected for w s ×c/M 1, where w s is fall velocity, c is concentration and M is the erosion parameter (Van Ledden, 2002).All the model scenarios have ratios below 1, which is in agreement with conditions in the Western Scheldt and probably in agreement with conditions in the Dovey given the clearly observable sand-mud transitions on imagery.
L. Braat et al.: Effects mud supply on large-scale estuarine morphology In the Western Scheldt the fluvial mud supply varies between 100 × 10 6 and 300 × 10 6 kg yr −1 at the Rupple mouth (Taverniers, 2000).In the model the mud input is 63 × 10 6 kg yr −1 .The mean discharge of the Scheldt, about 120 m 3 s −1 , is about 20 % higher than the default model scenario, while the sediment input is at least 60 % higher.This higher mud load might explain why the Western Scheldt has more mud deposits.In the field case, mudflats occur more on bars than on the sides compared to the models.We partly attribute this to the embankment and limited space to form mudflats and partly to spatially and temporally varying mud characteristics in the Western Scheldt.
The default scenario shows that the velocity amplitudes are flood dominant in shallow areas and ebb dominant in the channels (Fig. 4h).This is in general agreement with most earlier findings about tidal asymmetry (e.g.Speer and Aubrey, 1985;Friedrichs and Aubrey, 1988;Wang et al., 2002;Moore et al., 2009) including model studies on the Dovey (Robins and Davies, 2010;Brown andDavies, 2009, 2010).Our findings generalise these earlier trends because the estuary is self-formed.Several bathymetries tested in previous research are strongly simplified or arbitrarily chosen and might not represent a realistic state of an estuary, meaning that flood or ebb dominance could be the result of the imposed combination of initial conditions and boundary conditions.In contrast with our results, these case studies found higher flood peak velocities upstream (Brown and Davies, 2010;Robins and Davies, 2010).This is attributed to more intertidal area upstream that promotes flood dominance Moore et al. (2009); Brown and Davies (2010); Robins and Davies (2010).The default model showed stronger ebbdominant peak velocities in the landward part (Fig. 4g), which is caused by a higher river discharge in our model that causes ebb asymmetry.
Over time, the model evolved from a net exporting system to a dynamic equilibrium with balanced import and export (Fig. 3e and j).As more intertidal area and mudflats formed in the estuary, these areas gradually transformed from ebb-to flood-dominant peak velocities (Fig. 5).The mudflats particularly show much stronger flood-dominant peak velocities and a faster change over time than the intertidal area in general.This is because mudflats are significantly higher and have an elevation near the high water level, while typical sandy shoals only have a maximum height between the low and mean water level.This matches well with the sediment budget of the model that shows a net import of sediment resulting from mud import and sand export.This trend is also observed, most likely for the same reason, in the Western Scheldt on the basis of separate sand and mud balances (Cleveringa and Dam, 2013).Mud trapping is very efficient as the import is significant even though the duration asymmetry and peak velocity asymmetry are ebb dominant in most of the estuary.This again shows that the spatial variation in ebb and flood asymmetry is very important for understanding whether the estuary will grow or fill.Moreover, representa-tion of tidal asymmetry by width-averaged velocity ratios are insufficient and misleading in the presence of significant mud deposits.Due to mud deposition, the elevation of intertidal flats increases, which is therefore essential to change an estuary from exporting to importing or towards an equilibrium system.
Even though the tidal asymmetries in the model are comparable to many estuaries, waves are largely simplified.Different processes caused by waves promote flood dominance (e.g.Bertin et al., 2009;Nahon et al., 2012;Wargula et al., 2014).We expect that the inclusion of more wave processes on sediment transport would lead to faster development towards equilibrium by stimulating flood-directed transport.If the waves are very strong, we expect filling of the estuary by generation of a spit, and the estuary might never have been ebb-dominant in the first place.However, in the absence of waves, the continuous enlargement of estuaries with only sand might be as expected.

Transition from estuary to delta
The parameter space of Prandle et al. (2005) suggests that tides and river flow are sufficient conditions to explain the bathymetry of an estuary, with longer tidal reaches with larger river inflow (Fig. 12).This trend is not reproduced in the idealised model scenarios that typically have a tidal reach of 5-15 km in length but plot far above the line of 20 km in Fig. 12.Likewise, the trend is not clear in the dataset (Fig. 3 in Prandle et al., 2005).Rather, we observe the opposite trend: shorter estuaries or even deltas form with larger river discharges, and longer estuaries form in higher tidal ranges.Possibly, longer estuaries form for larger total flow from the combination of tide and river.We found much stronger effects of mud supply, suggesting that the tide-discharge parameter space needs to be extended with sediment supply.
As the model runs cover transgressive and regressive trends as effects of tides, river, waves and sediment supply on morphology we attempted to position the results in the traditional ternary classification diagrams for deltas of Galloway (1975).An expanded version of this classification system includes all coastal environments in which larger river influence leads to delta development and low or absent river influence leads to lagoons, strandplains and tidal flats (Dalrymple et al., 1992;Boyd et al., 1992).Qualitatively our results also show that for higher river discharge the estuarine system transitions to a deltaic system (Fig. 10c-i) by filling of the estuary.Note that the width did not decrease because a small tidal basin north of the river mouth affected the automated calculation of the width of the system (Fig. A4a).We also observed a transition to deltas when the tidal range was decreased (Fig. 10b-h) so that the relative power of the river increases in qualitative agreement with the classification diagram.
However, the most important findings of this research are more difficult to relate to these diagrams.We found that an increase in mud supply concentration leads to confining and filling of the initial estuary shape (Fig. 10a-g), leading to a decrease in total area and width at the mouth, while the mud-covered area and mudflat width at the mouth increased and is more delta-like.Orton and Reading (1993) found that smaller grain size leads to narrower channels in deltas and a tendency to avulse rather than have migrating channels.We observe similar behaviour in the model scenarios but here this is related not merely to grain size but to the supply rate.Alternatively, Dalrymple et al. (1992) and Boyd et al. (1992) developed a classification system with a fourth dimension based on the evolution of coastal systems by defining it as a progradating or transgressive system on the basis of sea level rise and sediment supply.This system disregards the possibility of an equilibrium without progradation and without transgression through combinations of sediment supply but otherwise similar hydrodynamic conditions.The models with different fluvial mud supply concentrations lead to distinct different morphologies but would plot on the same coordinates in these diagrams.Additionally, sea level rise is an ambiguous and qualitative variable in their conceptual figure because it affects the hydrodynamic conditions of the primary ternary diagram.To conclude, the model results for estuaries qualitatively fit in the ternary plots of Dalrymple et al. (1992) and Boyd et al. (1992) for deltas when sea level rise is ignored and sediment supply is considered the only variable on the fourth axis.

Large-scale equilibrium of estuary shape and dimensions
Estuaries with fluvial mud supply evolve into large-scale morphodynamic equilibrium (where absolute bathymetry change is constant, Fig. 3c, net bathymetry change is zero, Fig. 3d, and net export equals import, Fig. 4d) with dynamic  2015) with perpetually expanding tidal basins in cohesionless sand.After a rapid adjustment of basin size and bar and channel pattern the experiments developed to near equilibrium but never attained equality of sediment import and export.Our scenario without discharge is similar to these experiments and shows the same evolution, including the rapid adjustment and continuous erosion in a low dynamic state (Fig. A6d-VI).In braided rivers, such unhindered bank erosion leads to a "threshold channel" (Parker, 1978) with an equilibrium width related to the upstream flow discharge and the threshold for sediment motion.This theory was earlier suggested to be valid for tidal basins (Kleinhans et al., 2015).However, unlike rivers, estuaries are not limited by discharge because tidal prism can continue to increase as the estuary enlarges, leading to a potentially positive feedback only limited by friction.In other words, estuaries may expand to much larger systems because the tidal prism adapts to the estuary size and flow velocities and entrainment rates will not decrease with basin size unless opposed by cohesion.This proved to be the case in the models with mud.From this we conclude that development to an equilibrium shape for estuaries requires some form of apparent cohesion from mud, from species with sediment-binding effects and from non-erodible valley walls.This explains why previous studies found large-scale equilibrium in estuaries: these imposed a fixed estuary shape and size in 1-D simulations (e.g.Lanzoni and Seminara, 2002;Schuttelaars and De Swart, 2000;Todeschini et al., 2008) or imposed non-erodible boundaries in 2-DH (e.g.Hibma et al., 2003;Van der Wegen and Roelvink, 2008).
The novel model applications and results open up possibilities to incorporate the effects of species on flow and sediment transport (Van Oorschot et al., 2015), in which species and species density depend on substrate and salinity, and to unravel the effects of initial conditions inherited from early Holocene systems from the effects of boundary conditions (Townend, 2005).

Conclusions
The size and shape of alluvial river estuaries depend strongly on the supply of mud because this determines the mudflat formation that protects erodible estuarine boundaries against erosion.This was concluded from a series of idealised morphological model runs for medium-sized estuaries with sand and varying concentrations of mud, a range of tidal amplitudes and river discharges and limited littoral processes.Estuaries with mud supply may develop a dynamic morphological equilibrium.On the other hand, estuaries with only sand in the bed and banks expand perpetually with a positive feedback between tidal prism and sediment export.This means that freely developing estuaries self-confine their size and reduce channel and bar dynamics with increasing fluvial mud supply.Within centuries they attain a large-scale equilibrium with balanced sediment import and export.Higher mud supply concentrations result in shorter, shallower, narrower and generally smaller estuaries with increasing mudflat area and stronger funnelling that may develop into tidal deltas depending on the littoral conditions.Spatial patterns of mudflat development in estuaries depend strongly on whether the mud originates from the sea or the river: marine mud only influences the lower estuary with these model conditions, while fluvial mud deposits along the entire system in qualitative agreement with field data.The effect of marine mud supply is even smaller when waves are included, even though mud is transported further upstream.Tidal range and river discharge have opposing effects on the balance between mud deposition and erosion.For higher fluvial mud concentrations, relatively high river discharges and low tidal amplitudes, estuaries transition into prograding deltas.These general trends are similar to the effects of floodplain formation and erosion on the width and bar pattern in rivers.Earth Surf.Dynam., 5, 617-652, 2017 www.earth-surf-dynam.net/5/617/2017/www.earth-surf-dynam.net/5/617/2017/Earth Surf.Dynam., 5, 617-652, 2017  www.earth-surf-dynam.net/5/617/2017/Earth Surf.Dynam., 5, 617-652, 2017 Earth Surf.Dynam., 5, 617-652, 2017 www.earth-surf-dynam.net/5/617/2017/

Figure 1 .
Figure 1.Mud in the bed of the Western Scheldt and the Ems-Dollard.(a) Percentage of mud in the top 10 cm of the bed (McLaren, 1993, 1994), (b) GeoTOP map (v1.3) of probability of clay in the top 50 cm of the bed (TNO, 2016) and (c) an indicative morphodynamics map of the Western Scheldt (Rijkswaterstaat, 2012).(d) Fraction of mud in the top 10 cm of the bed (Van Heuvel, 1991; Rijkswaterstaat, 2009) and (e) GeoTOP map (v1.3) of probability of clay occurrence in the top 50 cm of the bed (TNO, 2016).(f) Surface mud distribution along the Western Scheldt from the three datasets.For the ecotope data only the low dynamics muddy class was used.(g-h) Cumulative and normalised hypsometric curves of surface area related to bed elevation.Plot includes the (cumulative and normalised) distribution of mud relative to the total area with reference to figure panels for the mud datasets.Dotted lines indicate high and low water levels during spring and neap tide at the mouth.

Figure 2 .
Figure 2. Initial bathymetry with model boundaries and cross section (red line) for analysis.Initial depth increases linearly upstream and width decreases exponentially (Eq.8).Coordinates are defined at the coastline with the channel centreline and mean sea level (MSL) as origin.

Figure 3 .
Figure 3. Results of the default scenario (01) with a fluvial mud supply of 20 mg L −1 , 3 m tidal range and 150 m 3 s −1 river discharge.(a-d) Bathymetry and (f-i) mud fraction in the top layer of the bed after 50, 150, 500 and 2000 years.(e) Morphodynamic activity expressed by absolute bed level change over time between the original coastline and the upstream boundary, (j) net bed level change over time between the coastline and the upstream boundary; positive is net accretion and negative is erosion.The ages of the maps (a-d and f-i) are indicated with blue dashed lines in panels (e) and (j).A video of the model is available on YouTube: https://youtu.be/HAeka4e2_PY

Figure 4 .
Figure 4. Hydrodynamics of the last day after 2000 years.The left panels show temporal variation in 1 day and the right panels show spatial variation along the estuary.(a) Water level, (b) streamwise flow velocity in the deepest channel with negative velocity towards the sea, (c) instantaneous discharge through the cross section and (d) cumulative sediment transport through the cross section showing no net difference between the upstream boundary and the coastline.(e) Maximum peak velocity for ebb and flood, (f) ebb and flood duration, (g) peak ebb and flood velocity ratio and (h) spatial pattern of peak velocity ratio showing flood-dominated shallow areas.The solid lines in panels (e) and (g) are based on streamwise velocity and the dotted line is based on velocity magnitude showing the effects of bends at 10 km.

628L.Figure 5 .
Figure5.Flood-ebb peak velocity ratio over time for the first 5 km of the default estuary (run 01) integrated over the total area (black squares), the area above mudflats (brown circles) and areas above (cyan triangles) and below (blue triangles) the −1 m bed level.

Figure 6 .
Figure 6.Effects of mud supply concentration (runs 03, 10, 01 and 09).The left column shows the final bathymetry of model runs after 2000 years and the right column shows mud fractions in the top layer of the bed.Runs with (a, e) 0, (b, f) 5, (c, g) 20 (default) and (d, h) 50 mg L −1 fluvial mud supply concentration.

Figure 7 .
Figure 7. Hypsometric curves of the final bathymetry after 2000 years.Curves indicate the cumulative area below a certain elevation.Dotted lines indicate the mud-covered area below this elevation.Runs 03, 10, 01 and 09 with 0, 5, 20 and 50 mg L −1 fluvial mud supply concentration.

Figure 8 .
Figure 8. Hydrodynamics and morphology along estuaries with different mud supply concentrations after 2000 years.From left to right (columns): model with only sand (03) and mud supply concentrations of 5 (10), 20 (default, 01) and 50 mg L −1 (09).(a-d) Minimum, mean and maximum bed elevation, high and low water level and minimum and maximum initial bed level, (e-h) braiding index and (i-l) estuary width defined as the initial width, maximum reach over the whole scenario run, the width of wet cells in the model, width defined by a threshold value that is used to mask the cells that are around the dry-wet cell threshold.(m-p) Intertidal area and mud cover as a percentage of the total area, (q-t) tidal range and (u-x) peak ebb and flood velocities.

Figure 9 .
Figure 9. Hydrodynamics and morphodynamics over time for estuaries with different mud supply concentrations.From left to right (columns): model with only sand (03) and mud supply concentrations of 5 (10), 20 (default, 01) and 50 mg L −1 (09).(a-d) Bathymetry of the cross section at the mouth plotted over time, (e-h) mud fraction in the top layer of the cross section at the mouth, (i-l) estuary width at 1, 4 and 8 km from the mouth, (m-p) funnel-shape parameter, (q-t) estuarine surface area, (u-x) intertidal area and mud in the bed relative to the total area, (y-II) absolute bed level change and (III-VI) net bed level change.

Figure 10 .
Figure 10.Most important large-scale morphological parameters after 2000 years as a function of the varied boundary conditions: fluvial mud supply concentration, tidal range and fluvial discharge.(a-c) Funnel-shape parameter,(d-f) mouth width (in blue colours) and mudflat width (brown colours) at the mouth and (g-i) total area (blue colours) and mud-covered area (brown colours).The data indicated in light blue and light brown are more conservative estimates as high mudflats (higher than 0.5 m below the high water level) are masked from the estuary shape from which area, width and funnel shape are calculated.Light grey areas indicate models in the transition from estuary to delta.Dark grey indicates models that evolved into a delta.

Figure 11 .
Figure 11.(a) Bar length versus partitioned bar width and (b) bar length against local estuary width.Model results plot in the same range as the data of the natural estuaries (Leuven et al., 2016).

Figure A2 .
Figure A2.Hydrodynamics and morphology along estuaries with different mud sources after 2000 years.From left to right (columns): model with only sand (03), marine mud supply (02), supply from both boundaries (04) and fluvial supply (default, 01).(a-d)Minimum, mean and maximum bed elevation, high and low water level and minimum and maximum initial bed level, (e-h) braiding index and (i-l) estuary width defined as the initial width, maximum reach over the whole scenario run, the width of wet cells in the model, width defined by a threshold value that is used to mask the cells that are around the dry-wet cell threshold.(m-p) Intertidal area and mud cover as a percentage of the total area, (q-t) tidal range and (u-x) peak ebb and flood velocities.

Figure A3 .
Figure A3.Hydrodynamics and morphodynamics over time for estuaries with different mud sources.From left to right (columns): model with only sand (03), marine mud supply (02), supply from both boundaries (04) and fluvial supply (default, 01).(a-d) Bathymetry of the cross section at the mouth plotted over time, (e-h) mud fraction in the top layer of the cross section at the mouth,(i-l) estuary width at 1, 4 and 8 km from the mouth, (m-p) funnel-shape parameter, (q-t) estuarine surface area, (u-x) intertidal area and mud in the bed relative to the total area, (y-II) absolute bed level change and (III-VI) net bed level change.

Figure A4 .
Figure A4.Effects of river discharge (runs 08, 01, 07 and 22).The left column shows the final bathymetry of model runs after 2000 years and the right column shows mud fractions in the top layer of the bed.Run with (a, e) 150, (b, f) 100 (default), (c, g) 50 and (d, h) 0 m 3 s −1 river discharge.

Figure A5 .
Figure A5.Hydrodynamics and morphology along estuaries with different discharge after 2000 yr.From left to right (columns): model with river discharge of 150 (08), 100 (default, 01), 50 (07) and 0 m 3 s −1 (22).(a-d) Minimum, mean and maximum bed elevation, high and low water level and minimum and maximum initial bed level, (e-h) braiding index and(i-l) estuary width defined as the initial width, maximum reach over the whole scenario run, the width of wet cells in the model, width defined by a threshold value that is used to mask the cells that are around the dry-wet cell threshold.(m-p) Intertidal area and mud cover as a percentage of the total area, (q-t) tidal range and (u-x) peak ebb and flood velocities.

Figure A6 .
Figure A6.Hydrodynamics and morphodynamics over time for estuaries with different discharge.From left to right (columns): model with river discharge of 150 (08), 100 (default, 01), 50 (07) and 0 m 3 s −1 (22).(a-d) Bathymetry of the cross section at the mouth plotted over time,(e-h) mud fraction in the top layer of the cross section at the mouth, (i-l) estuary width at 1, 4 and 8 km from the mouth, (m-p) funnel-shape parameter, (q-t) estuarine surface area, (u-x) intertidal area and mud in the bed relative to the total area, (y-II) absolute bed level change and (III-VI) net bed level change.

Figure A7 .
Figure A7.Effects of tidal range (runs 06, 01, 05 and 20).The left column shows the final bathymetry of model runs after 2000 years and the right column shows mud fractions in the top layer of the bed.Run with (a, e) 4 m, (b, f) 3 m (default), (c, g) 2 m and (d, h) 1 m tidal range.

Figure A8 .
Figure A8.Hydrodynamics and morphology along estuaries with different tidal ranges after 2000 years.From left to right (columns): model with 4 (06), 3 (default, 01), 2 (05) and 1 m (20) tidal range.(a-d) Minimum, mean and maximum bed elevation, high and low water level and minimum and maximum initial bed level, (e-h) braiding index and(i-l) estuary width defined as the initial width, maximum reach over the whole scenario run, the width of wet cells in the model, width defined by a threshold value that is used to mask the cells that are around the dry-wet cell threshold.(m-p) Intertidal area and mud cover as a percentage of the total area, (q-t) tidal range and (u-x) peak ebb and flood velocities.

Figure A9 .
Figure A9.Hydrodynamics and morphodynamics over time for estuaries with different tidal ranges.From left to right (columns): model with 4 (06), 3 (default, 01), 2 (05) and 1 m (20) tidal range.(a-d) Bathymetry of the cross section at the mouth plotted over time,(eh) mud fraction in the top layer of the cross section at the mouth, (i-l) estuary width at 1, 4 and 8 km from the mouth, (m-p) funnel-shape parameter, (q-t) estuarine surface area, (u-x) intertidal area and mud in the bed relative to the total area, (y-II) absolute bed level change and (III-VI) net bed level change.

Figure A10 .
Figure A10.Effects of mud source in the presence of waves (run 28, 27, 25 and 29).The left column shows the final bathymetry of model runs after 1250 years and the right column shows mud fractions in the top layer of the bed.Run with (a, e) only sand, (b, f) marine mud input (default), (c, g) marine and fluvial mud input and (d, h) fluvial mud input.

Figure A11 .
Figure A11.Hydrodynamics and morphology along estuaries for different mud sources in the presence of waves after 2000 years.From left to right (columns): model with only sand (28), marine mud supply (27), supply from both boundaries (25) and fluvial supply (29).(a-d)Minimum, mean and maximum bed elevation, high and low water level and minimum and maximum initial bed level, (e-h) braiding index and (i-l) estuary width defined as the initial width, maximum reach over the whole scenario run, the width of wet cells in the model, width defined by a threshold value that is used to mask the cells that are around the dry-wet cell threshold.(m-p) Intertidal area and mud cover as a percentage of the total area, (q-t) tidal range and (u-x) peak ebb and flood velocities.

Figure A12 .
Figure A12.Hydrodynamics and morphodynamics over time for estuaries for different mud sources in the presence of waves.From left to right (columns): model with only sand (28), marine mud supply (27), supply from both boundaries (25) and fluvial supply (29).(ad) Bathymetry of the cross section at the mouth plotted over time, (e-h) mud fraction in the top layer of the cross section at the mouth,(i-l) estuary width at 1, 4 and 8 km from the mouth, (m-p) funnel-shape parameter, (q-t) estuarine surface area, (u-x) intertidal area and mud in the bed relative to the total area, (y-II) absolute bed level change and (III-VI) net bed level change.

Table 1 .
Sediment characteristics applied in the default model.Variation in settling velocity will be discussed later as one of the sensitivity parameters.

Table 2 .
Parameters for processes and numerics.

Table 3 .
Overview of all model scenarios and runs to examine sensitivity to mud-related parameters.Run Marine mud Fluvial mud Tidal amplitude Discharge P m crit Act lyr thickness Settling velocity Waves Prandle et al. (2005)ed against river discharge for realworld estuaries and modelled scenarios.Field data are used fromPrandle et al. (2005)for estuaries in the UK and several other sources for different estuaries over the world.Lines indicate estimations of estuarine length byPrandle et al. (2005)of 5, 10 and 20 km from left to right.