Effects of discharge, wind, and tide on sedimentation in a recently restored tidal freshwater wetland

Sediment deposition is one of the key mechanisms to counteract the impact of sea level rise in tidal freshwater wetlands (TFWs). However, information about sediment deposition rates in TFWs is limited, especially for those located in the transition zone between the fluvially dominated and tidally dominated sections of a river delta where sedimentation rates are affected by the combined impact of river discharge, wind, and tides. Using a combined hydrodynamic–morphological model, we examined how hydrometeorological boundary conditions control sedimentation rates and patterns in a TFW located in the Rhine–Meuse estuary in the Netherlands. The modelling results show that net sedimentation rate increases with the magnitude of the river discharge, whereas stronger wind increasingly prevents sedimentation. Sediment trapping efficiency decreases for both increasing river discharge and wind magnitude. The impact of wind storms on the trapping efficiency becomes smaller for higher water discharge. The spatial sedimentation patterns are affected by all controls. Our study illustrates the importance of evaluating both the separate and the joint impact of discharge, wind, and tides when estimating sedimentation rates in a TFW affected by these controls. Such insights are relevant to design measures to reactivate the sedimentation process in these areas.

Yet little research has focused on TFWs located in the transition zone between the fluvially dominated and tidally dominated sections of a river delta (terms in italic as defined by Leonardi, Kolker, & Fagherazzi, 2015), where sedimentation rates are controlled by the combined impact of discharge, wind, and tide. The majority of previous studies has focused on cases where only one or two of these controls are relevant, for example, to study the relation between (1a) stationary water discharges and sediment deposition in a synthetic TFW (Nardin & Edmonds, 2014); (b) tidal ranges, sediment concentrations and bed-level changes in a TFW in the Scheldt estuary, Belgium/Netherlands (Temmerman, Govers, Meire, & Wartel, 2003a); and (c) wind waves and resuspension on a tidal mudflat in Willipa Bay, USA (Mariotti & Fagherazzi, 2013). However, a more thorough understanding of the combined impact of these hydrometeorological controls is essential to develop successful sedimentological restoration strategies in TFWs in the transition zone of a delta.
The objective of this study was to quantify and understand how sedimentation rates and patterns of mud and sand in TFW are affected by the interplay of river discharge, wind waves, and tide. To this end, we carried out numerical experiments using a hydrodynamic and sediment transport model of a recently restored, sparsely vegetated TFW in the south-western part of the Netherlands. We conducted 14 simulations with varying discharge, wind magnitude, and tidal conditions, and compared average surface accretion, trapping efficiency (defined as the proportion of the incoming sediment that is deposited or trapped in the area), and sedimentation patterns. This study concentrates on the first stage of renewed sedimentation of the TFW after the opening of the levees. Therefore, the effect of vegetation on the vertical mass balance (through increased sedimentation of suspended material and possible accretion due to production of organic material) is not considered.
This study took place within the framework of a larger project on the effects of restoring sedimentation in a former polder area, in which field measurements were carried out (water levels, flow velocities, turbidity, sediment concentrations, settling velocities, and sediment thickness). These measurements were used for the model set-up and calibration. The measurement programme and results will be described in a separate paper in preparation.

| STUDY AREA
The study area is located in the eastern section of De Biesbosch National Park, a 9,000 ha TWF in the lower part of the Rhine-Meuse delta in the Netherlands. The study area comprises three former polders (Spiering, Kleine Noordwaard, and Maltha) and has a surface area of around 700 ha. It was depoldered in 2008 by the park authority (State Forestry Service or Staatsbosbeheer in Dutch) as part of an ongoing programme to reduce flood water levels by enlarging inundation areas and to restore former wetland areas. Not only is this area itself potentially threatened by future SLR and, as such, a relevant case, the depoldering also created an excellent research environment to study sedimentation processes due to the size of the area and the limited number of in-and outlets, which facilitated the establishment of water and sediment balances.
The embankment around the polder was opened at two locations: on the northern side along the river Nieuwe Merwede (a major Rhine branch) close to location g1 in Figure 1, and on the southern side along the Gat van de Noorderklip (location g4), a smaller branch that connects with the Hollands Diep estuary. The dominant flow direction through the study area is from the North to the South. The area consists of inundated flats (former grassland and arable fields), a man-made channel system connecting the northern and southern in-and outlets, and a vegetated island in the centre of Kleine Noordwaard, which was constructed using the material dug from the channels. The substrate in the study area consists of a clay layer on top of a thick layer of fluvialtidal splay sands (Kleinhans, Weerts, & Cohen, 2010). Artificial channels were dug through this clay layer into the sandy layer underneath.
The hydraulic regime in the study area is semidiurnal microtidal with an average tidal range of 0.2 to 0.4 m. Because the TFW is located in the backwater of the North Sea, water levels in the TFW are influenced by storm surges as a result of heavy westerly wind storms at sea and the operation settings of the Haringvliet barrier (between Hollands Diep and the North Sea). The water levels are also affected by the discharge of the Rivers Rhine and Meuse. Most of the time, the tidal flats are inundated with depths ranging from 0 to 50 cm. Complete exposure of the flats only occurs in summer (when river discharge is low) at low tide or during strong easterly winds. The wave climate within the area is characterized by local short waves generated by winds mainly coming from the west-south-west. The significant wave height during windstorm events was observed to grow up to 0.2 m, as a result of the relatively long fetch lengths across the inundated flats and the distinct lack of vegetation, especially in winter. However, the development of the waves is hampered by the low water depths that occur during low tides or low river discharge.
The suspended sediment concentration (SSC) in the Nieuwe Merwede typically varies from 10 to 40 mg/L during average flow conditions with estimated peak values of up to 140 mg/L during periods of high discharge in the River Rhine (Asselman, 2000;Asselman, Middelkoop, & van Dijk, 2003).
We used Delft3D (Lesser, Roelvink, van Kester, & Stelling, 2004) to model hydrodynamics, sediment transport, and bed-level changes in the study area. The following sections describe the set-up of the model domain and the modules for hydrodynamics, sediment transport, and morphology.

| Model domain and bathymetry
The computational grid covers the polders Spiering, Kleine Noordwaard, and Maltha ( Figure 1). The land boundary largely follows the highest point of the original embankment around the polders. The upstream and downstream boundaries were chosen to coincide with the locations of fixed monitoring stations. The resolution and grid orientation were defined to account for dominant flow directions and important features in bathymetry (e.g., channels, island, in-and outlets). This resulted in a curvilinear grid of 144 × 145 cells, with cell sizes varying from 5 to 30 m ( Figure 2).
The initial bathymetry of the model area was constructed using the 2003 version of the official Dutch DEM "AHN1" with a horizontal resolution of 5 × 5 m 2 (Van der Zon, 2013), supplemented by a local LIDAR DEM with a horizontal resolution of 1 × 1 m 2 from 2010 for the central island and other artificially elevated areas, and a 2011 multibeam echosounder dataset covering the channel system. All bathymetric data sets were provided by the National Water Authority (Rijkswaterstaat).

| Hydrodynamics
Delft3D-FLOW calculates water levels and water flow velocities for every computational time step on spherical or orthogonal curvilinear coordinates by solving the unsteady shallow water equations in two or three dimensions. For this study, we used a curvilinear grid and a computational time step of 30 s. Given the explorative character of this study, the small gradient over observed vertical sediment concentration profiles and the focus on large scale horizontal sediment gradients, we decided to use depth-averaged simulations (2DH) to speed up the simulations.
We used the third-generation short wave model SWAN (Booij, Ris, & Holthuijsen, 1999) to simulate the effect of wind-driven short waves on hydrodynamics, morphology, and transport of sand and mud through increase in bed shear stress and wave-induced momentum.
This model is available within Delft3D as Delft3D-WAVE and calculates a wave field on the basis of hydrodynamic conditions and wind data. Delft3D-WAVE was coupled dynamically to Delft3D-FLOW, with the wave field being updated every hour. The effect of bottom friction in the energy balance equation in SWAN was included using the JONSWAP method described by Hasselmann, Barnett, Bouws, and Carlson (1973), and the wave-induced bed shear stress in FLOW was included using the method described by Fredsøe (1984).
Discharge time series were imposed at the two upstream The hydrodynamic roughness was defined using Manning's friction coefficient, which was set as a uniform value for the entire area after an a priori sensitivity analysis that showed that simulated levels and flows  (Deltares, 2014), except for the parameters listed in Table 1; the values for these parameters were defined on the basis of field observations and expert judgement.

| Sediment transport
Delft3D simulates suspended load transport of cohesive sediment fractions, suspended and bedload transport of noncohesive sediment fractions and the morphological changes that result from these processes. We defined one cohesive mud fraction and one noncohesive sand fraction. Transport of the noncohesive sediment fraction was modelled with the Van Rijn equation (Van Rijn, 1984). Uptake and settling of suspended sediment of the cohesive fraction were modelled with the Krone and Ariathurai-Partheniades formulations (Partheniades, 1965). The implementation of both transport formulae in the Delft3D framework is described by Lesser et al. (2004).
To obtain input SSC at the upstream boundary of the model, we used a sediment rating curve that was constructed following a procedure described by Asselman (2000), using discharges and SSCs from station Vuren, which is the closest river gauging station along the River Waal, located 31 km upstream of the study area. The SSC values estimated with the rating curve range from 20 mg/L for average river discharge (1500 m 3 /s) to 140 mg/L for extreme river discharge (6800 m 3 /s). The estimated SSC values for corresponding wetland inlet discharges (between 20 and 100 m 3 /s) agree well with SSCs at the inlet of the model area, which were measured using a calibrated turbidity sensor between July 2014 and April 2015.
The sediment densities of both fractions were left at default values (specific density of both mud and sand: 2650 kg/m 3 ; dry bed density mud: 500 kg/m 3 ; dry bed density sand: 1600 kg/m 3 ) as defined in Deltares (2013). On the basis of field observations, the effective settling velocity W s of the mud fraction was set at 0.04 mm/s, and the D 50 of the sand fraction at 200 μm.
The initial channel bed composition was modelled as one uniformly mixed layer with a spatially varying composition on the basis of field observations and geological maps of the area: 100% mud on the flats with a layer thickness of 2 cm, and 100% sand on the island and in the artificial channel system with a layer thickness of 3 m. The stiff polder clay layer underneath the mud layer was assumed to be non-erodible.

| Model calibration & validation
Model calibration and validation were carried out using a stepwise approach. The hydrodynamic model (including the SWAN wave module) was calibrated and validated first, subsequently, the bedload transport of the coarse fraction, and finally the suspended load transport and deposition. Because of limitations in data availability, different calibration periods were chosen for hydrodynamics and morphology.

| Hydrodynamic model
The hydrodynamic model was calibrated against observed water levels at three gauging stations in the case study area (points g2, g3, and g4 in Figure 1), using the root mean square error (RMSE) as optimization criterion. The Manning's roughness coefficient was used as calibration parameter with an a priori range of 0.01 to 0.05 s/m 1/3 . We selected the period between August 1, 2014 and December 1, 2014, as calibration period. August and September 2014 were relatively dry, apart from a few small discharge peaks in August. In late October 2014, there was a heavy windstorm event in combination with a small discharge peak. The calibration period ended with relatively dry and calm conditions. Discharge and water level series were derived from ADCP and diver measurements taken at the locations of the gauging stations. For more information on these time series, we refer to Van der Deijl (2015). Wind conditions (hourly values of average wind speed and direction during the last 10 min of every  The rest of the validation period was relatively dry with below-average discharges and no significant discharge peaks.
The SWAN model was not calibrated separately due to lack of quantitative data on wave characteristics. Instead, the performance was checked by comparing the significant wave heights for the calibration period with qualitative visual observations during field visits.

| Sediment transport model
Calibration of the sediment transport models was done at the level of 20 subareas (10 sections, see Figure 3, each further subdivided into a channel and a flat subsection), comprising the polder Kleine Noordwaard and polder Maltha. Table 2 lists the calibration parameters and their a priori value range. The specified ranges were based on a combination of available literature, expert judgement, and a priori sensitivity analysis. Calibration was carried out manually with the RMSE between measured and simulated accretion volumes (i.e., area-weighted accretion rates) in the 20 subareas over the entire calibration period as evaluation criterion.
The calibrated sediment transport model was evaluated using the Brier skill score (BSS; Sutherland, Peet, and Soulsby (2004), which is commonly used to evaluate the performance of a morphological model. We redefined the BSS in terms of bed volume changes rather than bed-level changes:  (2013)).

| Sensitivity analysis
To analyse the impact of different types of hydrometeorological controls on net sedimentation quantities and patterns, a sensitivity analysis was carried out (Table 3). The analysis was carried out for the following boundary conditions: • Discharge events of different duration and magnitudes. DISCH1 has a discharge peak with a return period of 1 year (T1) and DISCH2 a peak with a return period of 50 years (T50). Average tidal conditions (M2) apply in both cases, and there is no wind.
DISCH0 is the reference scenario: a stationary discharge event without any wind and with average tidal conditions. Chbab (2012), who describe extreme value statistics of river discharge, wind conditions, and sea levels in the Rhine delta. Figure 5 shows two examples of 1D boundary conditions that were constructed on the basis of these statistics. Next, using the 1D model output as boundary  Manning's roughness coefficient was set at 0.025 s/m 1/3 . Figure 6 shows the observed and simulated water levels using the calibrated model for the inlet of the Kleine Noordwaard (measurement location G2). Table 4 Figure 1). In (a), the discharge boundary Q (location Lobith) has a constant value of 2300 m 3 /s and is purposefully omitted. In (b), the wind velocity has a constant value of 0 m/s for this run and is purposefully omitted  Noordwaard, where the sand is deposited on the bar in-between the bifurcating channels ( Figure 8). Close to the outlet of the system (between g3 and g4 in Figure 1), the high shear stresses in the converging water flow cause the channel bed material to become mobilized and to leave the area both through bedload and suspended transport.
Thus, close to the inlet, mostly an internal redistribution of sand occurs, whereas close to the outlet, there is a net loss of sand from the area. In the channel system, further away from the in-and outlet shear stresses are too low for mobilizing or transporting sand, leading to stable channels. Overall, bed-level changes in the channel system are dominated by deposition and erosion of sand and exhibit a strong correlation with the magnitude of the discharge event ( Figure 9).
Sedimentation of mud on the flats mostly takes place close to the channels due the large gradient in flow velocity there (Figure 9). However, higher discharges deposit the material farther away from the channels due to the larger water depth and consequent smaller gradient in flow velocities. Local topographic irregularities in the bottom surface also affect the deposition patterns: The former drainage ditches prove to be very good sediment traps, and former roads or small dikes may prevent sediment loaded water from flowing back to the channels after the highest water levels have passed. This is especially notable for the heaviest discharge event, which inundates the entire system. Gradual sediment depletion causes a small gradient in sediment deposition from the inlet (more sedimentation) to the outlet (less sedimentation). Deposition of the mud inside the channels occurs only to a very small extent: mainly in the dead-end channels on the eastern side of the system.
Larger discharges cause both more erosion of mud in the channels and more sedimentation on the flats (Figure 9). There is no direct relation between these two effects: The bed level of the flats increases mostly because of sedimentation of silts coming from the upstream boundary, whereas the sand that erodes from the channels stays in suspension and leaves the area through the downstream boundary.
Sedimentation of suspended sand occurs only on a specific part of the flats close to the post-confluence channel section.
The total amount of mud retained in the study area increases with the magnitude of the discharge peak and corresponding increased influx of sediment whereas the mud trapping efficiency decreases ( Figure 10). The reduced trapping efficiency is caused by the increased shear stresses during the high-discharge events, which also causes most of the fines to stay in suspension during their transport through the channels in the area. Note. ME = mean error (cm); RMSE = root mean square error (cm); NSE = Nash-Sutcliffe efficiency index (−).

FIGURE 7
Cumulative sedimentation/erosion (in 1000 m 3 ) per subsection (numbered 1 to 10) over the calibration period. "C" stands for "channel," "F" for "flat". Note that the measured values for the flats are based on the estimation of a uniform accretion rate of 0.5 cm per year

| Windstorm events
The analysed windstorms have a large impact on the net sedimentation/erosion patterns compared to the reference case: All storm scenarios lead to less sedimentation on the flats ( Figure 11 The rest of the resuspended mud stays mobile and leaves the study area through the downstream outlet with the stationary water discharge. This leads to a net reduction in mud trapping efficiency for all windstorm events, regardless magnitude and direction, when compared to the reference case DISCH0 without wind (Figures 10 and 12). The only exception is WIND4. This event has the lowest wind speed of all T1 windstorm events as well as a wind direction that is opposite to the flow direction, which causes it to have a less pronounced impact on the mud trapping efficiency than the other T1 windstorm events ( Figure 11).

| Tidal range
Varying the tidal range from average to neap or spring tide has little effect on the sedimentation/erosion patterns compared to both reference cases (Figure 13), even though these patterns seem to be strongly affected by tidal water level fluctuations. The net retention rates of mud compared to both reference cases remain also relatively unaffected ( Figure 14). Still, the net outflow of sand from the study area during the discharge event depends slightly on the tidal amplitude, with a decrease during neap tide (TIDE3) and an increase during spring tide (TIDE4; Annex). The average channel bed level changes accordingly ( Figure 9).

| Combined discharge -Windstorm events
A T1 windstorm coinciding with a discharge event (COMB1) mobilizes the initial mud layer on the flats in a similar fashion as during average discharge conditions (e.g., WIND1), albeit to a slightly lesser extent.
This is the result of the increased water depth in the system due to the discharge event, which makes it more difficult for the wind waves to reach the bed. The wind also hinders the settling on the flats of "new" sediment entering the system with the discharge wave, which results in a negligible increase of average bed level of the flats ( Figure 9). Part of the resuspended sediment settles in the deeper parts of channel system which, in combination with the bed erosion occurring in other parts during the passing of the discharge event, leads to an almost neutral channel bed level change (Figure 9).
Another part of the mobilized mud is redistributed over the flats.
The sedimentation/erosion pattern of COMB1 strongly resembles WIND1 (Figure 15), although the total amount of sedimentation on Average daily inflow and outflow of mud, and the fraction of mud retained in the study area (trapping efficiency) for alternative WIND directions (WIND1 to WIND4) and magnitude (WIND5) DISCH0 is included for reference the flats is slightly higher in the combined case because of the passing of a discharge wave.
The rest of the resuspended mud remains mobile and leaves the study area through the downstream outlet, together with the sediment that originates from the upstream boundary and that was unable to settle on the flats due to the wind conditions. This leads to a large reduction in mud trapping efficiency when compared to reference case DISCH1 ( Figure 16). A much smaller windstorm (COMB2) still leads to a reduction of sediment trapping efficiency compared to the base case DISCH1, albeit very small.
Trapping efficiencies of combined discharge-windstorm events depend on both parameters as follows: Given an average discharge regime, switching from a T1 windstorm to a T50 windstorm leads to a large reduction in trapping efficiency (Figure 17). Given a T50 discharge event however, the trapping efficiency is reduced by only a small amount when switching from a T1 to T50 windstorm. We hypothesize that the larger water depths during the T50 discharge event reduce the impact of the waves on the bed shear stresses, thereby effectively reducing the resuspension of mud. Furthermore, increased sedimentation occurs in those areas that are less affected by wind, such as the deep dead-end channel sections and parts on the lee side of the island.

| DISCUSSION
Our site-specific results show the sedimentation patterns and rates inside a flow-through TFW, and how these change in response to various combinations of hydrometeorological controls. Until now, most research focused on a single control on sedimentation processes in a wetland (e.g., Mariotti & Fagherazzi, 2013;Temmerman et al., 2003b); however, our results demonstrate that there are cases in which there is important interaction among the three controls (discharge, wind, and tide). In this section, we discuss the role of both separate and combined controls, first on the sediment balance and trapping efficiency and then on sedimentation patterns.
For our study area, the role of the separate controls on the sediment balance terms and the sediment trapping efficiency is as follows: River input is the major source of sediment of this flow-through wetland, and the discharge through the inlet determines net sediment deposition rates (positive correlation, due to increased amounts of With respect to sedimentation patterns in the study area, our results generally agree with previous research in similar areas (e.g., Delgado et al., 2013;Hupp et al., 2008;Mitsch et al., 2014;Temmerman et al., 2003b): The highest sedimentation rates are found close to the inlet of the wetland and close to the internal channel network. Topographic irregularities in the submerged terrain-particularly former polder drainage ditches and old embankments-also influence local sedimentation patterns. Larger discharge events cause a larger portion of the sediment to settle farther away from the inlet and the channels. Wind influences sedimentation patterns by resuspension and internal redistribution, usually resulting in a net transport of sediment from the flats towards the channels.
The impact of wind on bed shear stresses is caused especially by locally generated wind waves and not by wind-driven currents. This is also observed in tidal mudflats, for example, along the Westerscheldt estuary (Callaghan et al., 2010). Interestingly, we find an almost linear Average daily inflow and outflow of mud, and the fraction of mud retained in the study area (trapping efficiency) for combined discharge-windstorm events COMB1 and COMB2. DISCH1 and WIND1 are included for reference FIGURE 17 Relation between discharge magnitude (defined by the maximum 12-hr moving average water inflow from the upstream boundary) and sediment trapping efficiency for three different windstorm scenarios (all from SW). The discharge values of 82, 175, and 352 m 3 /s correspond to average discharge, T1 and T50 discharge, respectively (negative) correlation between wind speed and mud retention, for all wind directions (Figure 18), even though fetch lengths for the various wind directions show large differences due to the distinctly elongated shape of the study area along the NW-SE axis. As Mariotti and Fagherazzi (2013) pointed out for the case of a tidal mud flat in Willapa Bay, USA, different fetches only start to impact bed shear stress above a certain critical water depth. We hypothesize that water depths in our study area (generally between 0.1 and 1 m) are below this critical value.
The determination of this critical depth is an interesting topic for further research.
To answer the question whether this TFW will survive the impact of SLR, the next step will be to analyse long-term effects for different climate scenarios in which the frequency and duration of these and other events are incorporated in longer time series and are adapted in accordance with future climate scenarios. This step also requires that the effect of vegetation and possibly subsidence be accounted for in the models. Application to other TFWs in the world requires further extension of the analysis by evaluating the effects on sedimentation rates and patterns of for example wetland size, shape, position within the delta (distance to turbidity maximum), and vegetation.

| CONCLUSIONS
To gain insight in both the average wetland surface accretion rates and the spatial sediment distribution in a TFW and the role of various different controls, we developed a combined hydrodynamic, morphological and wave model of a TFW in the Netherlands and applied it to analyse sediment rates and patterns for various windstorm and discharge events under different tidal conditions. The main conclusions for this area are as follows: • The net sediment deposition rate inside the TFW increases with water discharge magnitude and associated increases in SSC of the inflowing water, decreases with windstorm magnitude, and is relatively unaffected by changes in tidal conditions (neap and spring).
• The trapping efficiency decreases with water discharge magnitude as a result of increased bed shear stresses. Windstorms during any discharge event reduce the trapping efficiency compared to the same discharge event without any wind. The actual reduction increases with wind velocity, depends on wind direction (highest for SW winds), and decreases for higher water inflow from the river.
• Sedimentation rates are highest close to the inlet of the wetland, and the channel system within the wetland. Local sedimentation patterns are affected by irregularities in the topography, particularly former polder drainage channels and old embankments.
• Regardless of wind direction, windstorms lead to (a) a net transport of sediment from the flats towards the channels, (b) a net transport from the downwind sections of the flats to other sections, and (c) an increased outflow of sediment from the study area.
• TFWs have the potential to trap large amounts of sediment, yet the actual deposition rate shows large variations depending on the interplay between discharge conditions, windstorms and tidal conditions. This interplay should be taken into account when predicting long-term sedimentation rates.
Results are found to be in line with findings from previous studies.
However, the specific location of this wetland in-between tidally and fluvially dominated areas makes it particularly important to consider the various controls or boundary conditions in combination. Followup research will focus on (a) the identification of critical thresholds of combined boundary conditions for sedimentation and erosion in TFWs and (b) the application of the model to quantify the long-term response of the TFW to SLR.

ACKNOWLEDGMENTS
This study was financed by the Dutch Technology Foundation STW (project no. 12431). We thank everybody who assisted us with the data collection, especially Dr. H. de Boois. We thank the people of Staatsbosbeheer and Rijkwaterstaat-WNZ for the provided data, logistic support, and assistance. We also thank two anonymous referees for their careful reviews and constructive comments.