## Abstract

In this idealized numerical modeling study, the composition of residual sediment fluxes in energetic (e.g., weakly or periodically stratified) tidal estuaries is investigated by means of one-dimensional water column models, with some focus on the sediment availability. Scaling of the underlying dynamic equations shows dependence of the results on the Simpson number (relative strength of horizontal density gradient) and the Rouse number (relative settling velocity) as well as impacts of the Unsteadiness number (relative tidal frequency). Here, the parameter space given by the Simpson and Rouse numbers is mainly investigated. A simple analytical model based on the assumption of stationarity shows that for small Simpson and Rouse numbers sediment flux is down estuary and vice versa for large Simpson and Rouse numbers. A fully dynamic water column model coupled to a second-moment turbulence closure model allows to decompose the sediment flux profiles into contributions from the transport flux (product of subtidal velocity and sediment concentration profiles) and the fluctuation flux profiles (tidal covariance between current velocity and sediment concentration). Three different types of bottom sediment pools are distinguished to vary the sediment availability, by defining a time scale for complete sediment erosion. For short erosion times scales, the transport sediment flux may dominate, but for larger erosion time scales the fluctuation sediment flux largely dominates the tidal sediment flux. When quarter-diurnal components are added to the tidal forcing, up-estuary sediment fluxes are strongly increased for stronger and shorter flood tides and vice versa. The theoretical results are compared to field observations in a tidally energetic inlet.

## 1. Introduction

In tidal estuarine systems with net buoyancy input from land (rivers, groundwater) or in shallow coastal waters (net precipitation, net warming), subtidal landward sediment fluxes are common. On longer time scales, the estuarine sediment budgets are closed by a seaward sediment flux due to dispersion caused by the positive landward sediment gradient, or by net sedimentation due to increased estuarine sediment concentrations. In many estuaries, sediment reduction caused by bacterial decomposition of particulate organic matter or dredging may play a substantial role.

Estuarine turbidity maxima (ETMs) are one characteristic consequence of estuarine sediment fluxes. They often occur at the landward end of the salinity intrusion such as in the Columbia (Jay and Musiak 1994), Hudson (Geyer et al. 2001), Elbe (Kappenberg et al. 1995), and Weser (Wellershaus 1981) Rivers, and are supposed to be driven by processes related to the internal density gradient. In shallow tidally energetic estuaries, ETMs may also occur in the freshwater range upstream of the salinity intrusion such as in the Ems (Talke et al. 2009; Schuttelaars et al. 2013) and Gironde (Allen et al. 1980) Rivers, where they are believed to be driven by asymmetries in tidal energy. Landward sediment fluxes may also be responsible for the retention and accumulation of pollutants, which tend to adhere to suspended particulate matter (Kuncheva et al. 2000). Also, in regions of freshwater influence (ROFIs) (Simpson 1997), landward sediment fluxes are ubiquitous features such as in Liverpool Bay (Souza and Lane 2013).

Landward net sediment fluxes are believed to prevent losses of intertidal flats due to sea level rise in shallow tidally energetic coastal environments such as the Wadden Sea of the southeastern North Sea (Postma 1961, 1981). The net import of the organic fraction of the particulate matter into the Wadden Sea plays an important role for the nutrient supply with possibly far-reaching consequences for the filter feeder communities (van Beusekom and de Jonge 2002).

There is a variety of processes that contribute in driving subtidal landward sediment fluxes. In well-mixed shallow water, tidal velocities and bed stresses are typically larger than the ebb values leading to spatial (Postma 1954; van Straaten and Kuenen 1958) and temporal (Groen 1967) settling lag as well as to scour lag (Bartholdy 2000), processes leading to net onshore sediment transports.

In estuaries, the superposition of the tidally averaged external pressure gradient (because of the mean surface slope, pointing seawards with no vertical variation) and the tidally averaged internal pressure gradient (because of horizontal density gradients, pointing landwards with increasing magnitude toward the bed) results in near-bed landward and near-surface seaward advective transports (Hansen and Rattray 1965) driving advective landward sediment fluxes in regions of sufficiently strong density gradient forcing. The near-bed landward currents may be enhanced by tidal asymmetries in eddy viscosity (Jay and Musiak 1994) caused by tidal straining (Simpson et al. 1990) leading to destratification, enhancement of turbulence, and thus increased mixing during flood and vice versa during ebb. Further enhancement of near-bed landward currents may come from lateral circulation systematically varying between ebb and flood (Lerczak and Geyer 2004) or from seaward wind stress causing wind straining (Scully et al. 2005). In an idealized numerical study, Burchard et al. (2011) could show that—for tidally energetic estuaries with modest freshwater inflow destratifying during each tidal cycle—subtidal circulation resulting from tidal straining was strongest, but for partially mixed estuaries gravitational and lateral circulation processes were more important.

The process of tidal straining does also lead to mixing of suspended matter higher up into the water column during flood than ebb, with the consequence of further-enhanced landward sediment fluxes as a result of tidal pumping (Scully and Friedrichs 2007b). A related process supporting ETM formation by convergent sediment fluxes has been proposed by Geyer (1993): in the well-mixed regime upstream of the landward end of the salinity intrusion, sediment is mixed high up into the water column and thus transported seawards, but once it enters the stratified regime of the salt wedge, the particulate matter sinks down to the near-bottom region where it is subject to landward advective fluxes. In an idealized model study, Burchard and Baumert (1998) show the importance of tidal straining for ETM formation and conclude that tidal pumping and gravitational circulation are less important.

Scully and Friedrichs (2007b) discuss the additional hypothesis that the increased breakup of suspended particulate matter flocs—because of more energetic turbulence during flood (and thus smaller Kolmogorov scales) and increased flocculation during ebb—could lead to higher settling velocities of suspended matter during ebb, which would be a further process supporting landward sediment fluxes in estuarine waters subject to horizontal density gradients.

To model equilibrium conditions concerning the suspended sediment dynamics in tidal estuaries, the sediment availability (or erodibility) must vary along the estuary (Friedrichs et al. 1998). The influence of sediment availability on the ETM dynamics was already observed in the work of Lang et al. (1989). Consequently, Burchard and Baumert (1998) used a limited amount of sediment in their model to study the equilibrium dynamics of an ETM. In the analytical approach of Huijts et al. (2006) and Chernetsky et al. (2010), the inclusion of a spatially varying erosion coefficient is essential to describe the trapping of fine sediment in the lateral and longitudinal directions, respectively. Limited sediment availability is reflected by the observation that erodibility decreases with depth within the sediment layer and can be expressed as a critical shear stress increasing with depth (Sheng and Villaret 1989; Cartwright et al. 2009; Dickhudt et al. 2009). However, the influence of sediment availability on residual sediment fluxes in tidal estuaries has not been investigated in detail.

Therefore, the major motivation of this study is to systematically investigate the influence of sediment availability on residual sediment fluxes in tidally energetic (e.g., weakly or periodically stratified) estuaries, with a focus on estuarine circulation dynamics and sediment flux composition. To obtain a good coverage of the complex parameter space, and to include a number of processes potentially relevant for sediment transport under the assumption of horizontal homogeneity, a one-dimensional numerical water column model is used (see section 2 for the model equations and section 3 for model setup and the residual flow profile decomposition method). This approach has the further advantage that discretization errors may be kept small because of high numerical resolution. Because tidal asymmetries are pronounced in estuaries and tidal inlets, the sensitivity of sediment flux dynamics to the M_{4} tidal phase is studied in detail (see section 5a for the reference simulations with symmetric forcing and section 5b for the simulations with asymmetric tidal forcing). The drawback of such a one-dimensional approach is certainly that a number of relevant processes owing to spatial inhomogeneities and temporal irregularities is excluded, a fact that has to be kept in mind when interpreting the results of this study in the light of real world estuaries and tidal inlets (see the discussions in the field experiment section 6 and the conclusions section 7). Furthermore, the feedback of estuarine dynamics to horizontal buoyancy gradients is neglected with the consequence that this model study does not represent processes in permanently stratified estuaries. High Simpson numbers would lead to the known problem of run-away stratification for one-dimensional models (Blaise and Deleersnijder 2008).

## 2. Model equations

There are three major parameters governing the dynamics of a tidal flow in a horizontally homogeneous water column under the influence of a horizontal buoyancy gradient, the Simpson number

the Unsteadiness number

and the Rouse number

with the horizontal buoyancy gradient [∂_{x}*b*] (note that square brackets denote diagnostically prescribed quantities), the water depth *H*, the friction velocity scale (representing the square root of the mean tidal stress), the tidal frequency *ω*, the constant settling velocity *w _{s}* (positive for sinking), and the von Kármán constant

*κ*= 0.4. The conversion between the friction velocity scale and the current velocity scale

*U*

_{2}(representing the semidiurnal current velocity amplitude) is obtained by assuming a logarithmic velocity profile at maximum current speed with a bulk drag coefficient

see Burchard et al. (2011). In (4), is the nondimensional bottom roughness. The factor 2 in (4) accounts for the fact that represents the square root of the tidal mean stress (note that the average of the cosine squared function is ½).

The Simpson number represents the balance between the mixing force of vertical stress divergence and the stratifying force of the horizontal density gradient and is known to be of the order of unity for the transition between mixed and stratified conditions (Stacey et al. 2001). It should, however, be noted that the critical Simpson number for transition to permanently stratified conditions depends on the local bathymetry. It is known that for channelized estuaries additional processes such as lateral circulation (Lerczak and Geyer 2004; Scully et al. 2009; Burchard et al. 2011), see also section 6, provide additional forcing to estuarine circulation such that Simpson numbers already significantly lower than unity may lead to stratified conditions. Because such processes are excluded in this one-dimensional modeling study, transition to stratified conditions is expected for Si = *O*(1). The Unsteadiness number represents the balance between stress divergence and local acceleration, or, as an alternative interpretation, the ratio of the water depth to the tidal boundary layer height. For tidal flow, this ratio is typically much smaller than unity. It should be noted that there are arguments (Friedrichs 2010) to formulate the Unsteadiness number based on the velocity scale instead of , which would result in an Unsteadiness number for tidal flow of the order of unity. The Rouse number describes the balance between settling and upward turbulent transport of sediment. A further relevant parameter is the nondimensional bottom roughness length , because it determines the ratio between the flow velocity and the bottom friction velocity [see (4)].

Typical dimensional scales for energetic tidal channels are collected in Table 1 for a number of estuaries and tidal inlets. Based on these numbers, a tidal period of *T* = 44 714 s; a bottom roughness of ; a typical range of settling velocities *w _{s}* is from 1 × 10

^{−3}to 2 × 10

^{−3}m s

^{−1}[as observed by Scully and Friedrichs (2007b) for the flood–ebb variability in the York River estuary]; and Si, Un, and Ro are calculated with (1)–(3), using (4) to calculate the friction velocity scale from the tidal velocity amplitude scale

*U*

_{2}. For these estuaries and tidal inlets, which can be categorized as weakly-to-periodically stratified, the nondimensional parameter ranges are 0.11 ≤ Si ≤ 0.61, 0.01 ≤ Un ≤ 0.07, and 0.05 ≤ Ro ≤ 0.36, where the channel site of the York River estuary alone (Scully and Friedrichs 2007b) covers most of the range within its spring–neap cycle. It should be noted that the values for these estuaries are examples and may strongly vary with the spring–neap cycle, river runoff, temperature, and many other external parameters.

Using [∂_{x}*b*], *H*, , and *ω* to make the physical variables dimensionless, the nondimensional water column dynamic equations for the along-current velocity {with dimensional velocity *u*, the friction velocity scale , and the buoyancy , where dimensional buoyancy *b* = −*g*(*ρ* − *ρ*_{0})/*ρ*_{0}, gravitational acceleration is *g*, density is *ρ*, and reference density is *ρ*_{0}} read as

with nondimensional time (dimensional time *t*) and nondimensional vertical coordinate (dimensional vertical coordinate *z*) [see, e.g., Burchard (2009)]. In (5) and (6), and are the nondimensional eddy viscosity and eddy diffusivity, respectively, and is a time-dependent barotropic pressure gradient constructed in a way that

with the nondimensional depth-averaged tidal velocity amplitude and the nondimensional residual runoff velocity , see Burchard (2009). Here, and are nondimensional amplitude and relative phase of an M_{4} tidal component, which is added for sensitivity studies. Stratification is calculated as a result of the vertical buoyancy gradient:

Analogously, a dynamic equation for the dimensionless sediment concentration can be derived [see, e.g., Burchard and Baumert (1998) for a dimensional version]:

with the Rouse number defined in (3). The sediment concentration has been made dimensionless by the tidally and vertically averaged sediment concentration (which can only be calculated a posteriori):

with the tidal period *T* = 2*π*/*ω*.

At the bottom and the surface, zero-flux conditions are applied for the buoyancy . For velocity, a zero-flux condition is applied at the surface and a no-slip condition at the bed . For the sediment concentration, the classical sedimentation–erosion formulation (Krone 1962) is used at the bottom, while at the surface a zero-flux condition is applied (see appendix for details).

To limit sediment availability as it is observed in real estuaries (Friedrichs et al. 1998; Dickhudt et al. 2009), a bottom pool of sediment is tracked, which may be emptied with the consequence that erosion may cease also at high-bottom shear stresses. It is shown in the appendix how this bottom pool can be related to a depth-varying critical shear stress (expressed as a depth-varying critical bottom friction velocity).

We furthermore introduce a nondimensional sediment erosion time scale , which indicates how much of a tidal period it would take to empty the bottom sediment pool when applying a tidal-mean bed stress (see appendix for details). In this study three cases will be investigated, (bottom pool is empty almost all the time), (bottom pool is empty during some parts of the tidal period), and (bottom pool is never empty).

## 3. Methods

### a. Model setup

To consider nonstationary processes that involve tidally varying eddy viscosity and diffusivity, a one-dimensional numerical model is applied to solve the dynamic (5), (6), and (9). This model uses a two-equation turbulence closure model with an algebraic second-moment closure, as described in detail by Umlauf and Burchard (2005). To prevent unrealistic collapse of turbulence in regions of high stratification, the macro-length scale of turbulence was limited to the Ozmidov-length scale, and a minimum TKE threshold value of *k*_{min} = 10^{−5} J kg^{−1} was set. This procedure parameterizes unresolved background mixing because of small-scale internal waves, see Umlauf and Burchard (2005) for details. Application of this parameterization prevents the eddy viscosity from falling below values of *O*(10^{−6} m s^{−2}).

The application to tidal estuarine flow is described in Burchard (2009) and Burchard and Hetland (2010). For too-strong horizontal buoyancy gradient forcing, that is, for a too-large an Si number, the flow showed a run-away stratification because of the effect of the constant horizontal buoyancy gradient [∂_{x}*b*], which in real estuaries is adapting in a way that periodic states may also exist under high-Si conditions. Therefore, also high-Si situations will be excluded here, because their study requires full-3D simulations.

After having prescribed for each simulation values of *H* = 10 m, *ω* = 1.4 × 10^{−4} s^{−1} (the M_{2} tidal frequency), (according to the law of the wall equivalent to a 1-m drag coefficient of 2.8 × 10^{−3}), Un = 0.05, Si (variable), and Ro (variable); the friction velocity scale, the horizontal buoyancy gradient, and the settling velocity can be obtained from these variables by combining (1), (2), and (3):

The tidal velocity amplitude *U*_{2} = 1.25 m s^{−1} is then obtained by (4) using = 0.028 m s^{−1} from (11). All equations are solved in dimensional form and the results are presented as nondimensional values.

The numerical resolution is carried out using 200 equidistant vertical layers and the tidal cycle is resolved by 1000 equidistant time steps. The sediment concentration is initialized with an empty bottom sediment pool .

For convenience, all numerically simulated velocity data shown in this study are displayed as normalized by the tidal velocity amplitude *U*_{2} instead of the friction velocity scale :

### b. Residual flow profile decomposition

The residual dynamics of this reference scenario are analyzed using the decomposition method suggested by Burchard and Hetland (2010). This method allows decomposition of the residual flow profile into components driven by distinct processes:

where angle brackets denote tidal averages with

for any *z*-dependent quantity *X* with the time-dependent deviation from tidal averages *X*′. In (13), all variables are functions of *z* and the subscripts denote the processes leading to the residual flow profile components (in angle brackets, the major variable is given on which the process depends):

*s*: covariance of eddy viscosity and vertical shear (depending on and ), which is the tidal straining circulation term;*g*: gravitational circulation (depending on [∂_{x}*b*] and ); and*r*: runoff circulation (depending on and ).

In real estuaries, a number of other processes will drive or influence estuarine circulation such as lateral bathymetry variation, channel curvature, Earth's rotation, wind stress, nonlinear tides, etc., which are, however, not considered in the present study.

Assuming that the sediment concentrations are horizontally homogeneous, the longitudinal tidal-mean sediment flux profile can be decomposed as follows, using (13):

showing that the total sediment flux profile is composed of a fluctuation sediment flux (covariance between current velocity and sediment concentration) and a number of flux profiles due to tidal-mean advective transports of tidal-mean sediment concentration. In the stationary analytical scenario presented in the following section 4, only the latter two sediment flux profiles due to gravitational and runoff circulation have been considered.

## 4. Analytical solutions of steady-state conditions

The following stationary eddy viscosity and diffusivity profiles allow for analytical steady-state solutions for the velocity and the sediment profile:

with the nondimensional distance from the bed . It should be noted that the effective turbulent Prandtl number in this analytical test scenario is . The factor is introduced to obtain the well-known Rouse solution. Burchard and Hetland (2010) derived an analytical stationary solution of (5) using the eddy viscosity of (16):

with the integration constant

In (17), the first term is an estuarine exchange flow with zero vertical mean due to a balance of friction, and (barotropic and baroclinic) pressure gradient, and the second term is a logarithmic velocity profile, where the vertical average equals the runoff velocity . Note that exchange flow due to tidal straining is not included in this steady-state solution. A typical exchange flow resulting from (17) is shown in Fig. 1a.

With the parabolic eddy diffusivity profile (16), see Fig. 1b, and a zero-flux bottom boundary condition for the sediment concentration, that is, , the Rouse sediment profile is obtained as stationary solution of (9):

(see Fig. 1c), with the nondimensional bottom concentration . Note that because of the zero-flux bottom boundary condition, we obtain for with (A2):

The sediment profile is made here nondimensional according to (10), such that has a vertical mean of unity. Because no analytical expression can be found for the vertical integral of (19), needs to be obtained from *c* by numerical integration. Figure 1d shows the resulting analytical profile for the advective sediment flux for one parameter combination of Si, Ro, , and . With this parameter combination, the down-estuary sediment flux due to runoff is about balanced by the up-estuary sediment flux due to the exchange flow.

It should be noted that the analytical solutions for a constant eddy viscosity and eddy viscosity would result in the well-known third-order polynomials for the exchange flow (Hansen and Rattray 1965) and an exponential profile for the sediment concentration.

Results for vertically integrated sediment flux as function of Si and Ro are shown in Fig. 2. The sediment flux increases linearly with Si, which follows from (17), and has a maximum for Rouse numbers slightly below unity, with the optimum Rouse number slightly decreasing with increasing Simpson number. Down-estuary sediment flux occurs for low Si numbers (weak estuarine circulation) and low Ro numbers (small increase of sediment concentration toward the bottom).

## 5. Numerical simulations

### a. Reference simulations

The reference simulations are carried out with a pure M_{2} barotropic tidal forcing plus a weak river runoff, with high-resolution coverage of the Ro–Si parameter space, neglecting sediment cohesiveness and sediment density. Note that internal nonlinear dynamics due to time-dependent eddy viscosity will generate higher harmonics in the current velocity profiles the vertical mean of which will vanish. For the interaction with the bed, the sedimentation–erosion formulation (A1) is used with a zero critical shear stress , meaning that the sediment erosion flux is proportional to the shear stress. Three different erosion time scales as defined in (A6) will be used to reproduce unlimited bed sediment supply , limited bed sediment supply , or almost no bed sediment supply . Effects of variations in the Unsteadiness number, the bottom roughness, and the residual runoff are briefly discussed.

Figure 3 shows tidally resolved profiles for , , *Ã _{υ}*, and for Si = 0.8, Ro = 0.2, , and , which generates a typical strain-induced periodic stratification (SIPS) situation (note that for channelized estuaries the transition between SIPS and permanently stratified conditions may occur at a lower Si, see section 2). The water column is significantly stratified after ebb (for [∂

_{x}

*b*] = 6.6 × 10

^{−6}s

^{−2}and

*H*= 10 m, a maximum stratification from top to bottom of Δ

*ρ*= 1.67 kg m

^{−3}or is obtained) and destratified after flood. Eddy viscosity peaks at flood to a value that is about 75% higher during flood than during ebb. The eddy diffusivity that is driving the vertical turbulent salt and sediment transport is of similar shape because the turbulent Prandtl number is on the order of unity. The sediment dynamics does strongly depend on the availability of sediments from the bed (Figs. 3d–f). For unlimited sediment supply (Fig. 3d), bed sediment concentrations are highest during full flood and ebb, because then bottom shear stresses are at maximum values and consequently erosion is highest. For almost no sediment supply from the bed (Fig. 3f), near-bed sediment concentrations peak during early ebb and flood, because for fully developed tidal currents, sediments are mixed up into the water column, leading to relatively low near-bottom concentrations. The case for partially limited bed sediment supply (Fig. 3e) lies somehow in between these two extremes.

To give a better insight into the effect of the sediment bottom pool on the sediment dynamics in the water column, time series of unscaled vertically integrated sediment concentrations in the water column are shown for three different initial nondimensional sediment concentrations:

(equivalent to and , i.e., almost always empty bottom pool)

(equivalent to and , i.e., partially empty bottom pool)

(equivalent to and , i.e., never empty bottom pool).

Figure 4a shows that the erosional behavior after slack tide at and are identical for the three cases as long as the bottom pool is not emptied.

However, during decelerating tides the amount of sediment in the water column needs more time to settle for higher sediment concentrations as a result of temporal settling lag effects. To eliminate this time lag effect, and as consistency check, the three simulations are repeated using only for the sediment equation a 10-fold smaller Unsteadiness number which is equivalent to running the sediment dynamics on a 10-fold larger time scale with the consequence that the sediment dynamics is in quasi equilibrium. The result is shown in Fig. 4b: note also that the settling dynamics of the three scenarios is identical as long as the bottom sediment pool is not exhausted.

Figure 5 shows for well-mixed flow with Si = 0.0, that flood stresses are significantly weaker than ebb stresses, because of the effect of the runoff. In contrast to this, maximum bed stresses during flood are marginally larger than maximum bed stresses during ebb for Si = 0.6 and Si = 0.8. This is a result of a competition between down-estuary residual runoff and up-estuary residual near-bottom flow owing to estuarine circulation. Because eddy viscosity scales with and *H*, the much higher flood eddy viscosities in the water column cannot be explained by asymmetric bed stresses. Instead, this asymmetry is generated by tidal straining destabilizing the water column during flood and thus producing an upward buoyancy flux, see, for example, Burchard (2009). In addition, the flood is significantly longer than the ebb near the bed because of the fact that already before the reversal of the depth-mean flow after ebb when turbulent mixing is small as a result of vertical stratification, the horizontal density gradient moves near-bottom water upstream. This effect is not visible at the end of flood because vertical mixing is still high.

The decomposition of the sediment flux profiles for Si = 0.8 and Ro = 0.2 is shown in Fig. 6. As expected for this small runoff estuarine circulation generates significant up-estuary near-bottom residual circulation because up-estuary tidal straining circulation and gravitational circulation are larger than down-estuary river runoff (Fig. 6a).

The tidal mean sediment profiles (Figs. 6b–d) have similar shapes as the Rouse profile shown in Fig. 1, with the profile for the short erosion time scale (Fig. 6d) being less steep than the profile for unlimited sediment supply (Fig. 6b).

Figures 6e–g show for all three erosion time scales that the total sediment flux is dominated by the fluctuation sediment flux in the upper part of the water column, which is always up estuary. The reason for this is that during flood, as a consequence of the much larger flood eddy diffusivity, a significantly higher amount of sediments is transported up into the water column than during ebb. For the short erosion time scale (Fig. 6g), near-bottom fluctuation sediment flux is negative because near-bed sediment concentrations are smaller during flood than during ebb, which is a result of the ceased sediment flux from the bed and the vertically more homogeneous flood tide sediment profiles. For the other extreme of unlimited sediment supply from the bed (Fig. 6e), fluctuation sediment flux is highest during full flood, which in turn generated highest near-bed sediment concentrations during flood and a significantly positive residual sediment flux also near the bed. The sediment flux profiles for the partially limited bed erosion (Fig. 6f) is somewhat in between these two extremes.

These profiles of the fluctuation sediment flux are already captured taking only the correlation between the M_{2} tidal velocity and concentration into account, independent of the erosion time scale. To get the amount of sediment transport correct using tidal constituents, it is enough to take the first three tidal components into account (i.e., the M_{2}, M_{4}, and M_{6} tidal constituents): this results for all erosion time scales on an error of less than one percent, see Fig. 7. From this figure, it follows that only in case of a short erosion time scale (i.e., when the sediment pool is emptied), an approximation with an M_{2} tide only would result in a large error. For the other erosion time scales, only considering the M_{2} tidal constituent would result in a good reproduction of both the transport profiles and the depth-integrated amount of sediment transported (relative error of ≈3%).

The sediment flux profiles due to contributions from mean advective transports do only slightly vary between the case with different erosion time scales, because they are all based on the same residual flow profiles (Fig. 6a) and on only slightly different residual sediment profiles (Figs. 6b–d). The contributions from gravitational circulation and tidal straining are up estuary near the bed and down estuary near the surface, with a dominance of the up-estuary part, whereas the sediment flux due to runoff is down estuary over the whole water column.

The vertically integrated sediment fluxes for the flux profiles shown in the figure (Figs. 6e–g) are quantified in Table 2. For shorter erosion time scales, the total sediment flux decreases as well as the contribution from the fluctuation sediment flux . However, the advective contributions from gravitational circulation and tidal straining circulation are slightly increasing for shorter erosion time scales because the tidal mean sediment profiles are steeper for that case (see discussion above). For depleted bed sediment supply , the residual sediment flux due to straining is of almost the same magnitude as the fluctuation flux.

A more continuous picture is given by the Ro–Si parameter space study given in Fig. 8 showing tidally and vertically integrated sediment fluxes and their decomposition, resulting from 26 × 41 individual model simulations. The maximum Si number chosen here was 0.8, because for Si > Si_{c} with the critical Simpson number Si = 0.9 run-away stratification occurs for this one-dimensional modeling study.

For the total sediment flux (Figs. 8a,d,g), the principle picture is that for low Si and for low Ro numbers the residual sediment flux is down estuary, independent of the erosion time scale . With this, the result from the analytical solution presented in Fig. 2 is qualitatively reproduced, although the latter neglects the fluctuation sediment flux and the straining sediment flux. For intermediate Si numbers, scenarios with unlimited bed sediment supply show a maximum down-estuary sediment transport dominated by the fluctuation flux for Ro numbers around 0.1. This can be explained in a way that for Ro ≪ 0.1 the fluctuation flux is more inefficient because sediments will be homogeneously distributed in the water column and the exchange with the bed is suppressed because the sedimentation flux is reduced because it is proportional to Ro, see (A2). For Ro ≫ 0.1 when exchange with the bed is strongly enhanced, the bed stress asymmetry will favor up-estuary fluctuation sediment fluxes for Si > 0.4, when flood bottom stress is stronger than ebb bottom stress. For an intermediate erosion time scale (Figs. 8d–f), bed sediment supply is not limiting for Ro > 0.35 such that for this case the sediment fluxes are similar to the fully unlimited case with . The strong variation of the total sediment fluxes with Si, Ro, and indicates that in real estuaries a strong spring–neap variation is expected (varying Si and Ro), as well as a significant sediment sorting (including various sediment classes with different settling velocity and therefore different Ro) and spatial variability (variation in and other parameters).

The sensitivity of the results for changes in the Unsteadiness number is such that the threshold Simpson number for permanent stratification Si_{c} is lowered for decreased Un and vice versa. This can be interpreted in a way that a higher relative tidal frequency (equivalent to a higher Un) allows for more effective vertical mixing in estuaries, a finding which has already been discussed by Burchard (2009).

Increasing the roughness length for Un = 0.01 decreases the critical Simpson number and vice versa. This is mainly an effect of the increase of the friction velocity scale with the roughness length, see (4) which decreases the Simpson number (and the Unsteadiness number) for a given horizontal buoyancy gradient.

Runoff variations affect the sediment fluxes in a way that the threshold between down- and up-estuary sediment flux increases about linearly with the runoff velocity. For zero runoff, all components of the sediment flux are essentially up estuary.

### b. Effects of *M*_{4} barotropic forcing

_{4}

The following scenarios with an additional M_{4} barotropic tidal forcing according to (7) are carried out with an amplitude of [which seems to be a typical value for tidal inlets, see, e.g., Davies et al. (1997)] and four different phases :

transition from ebb to flood shorter than from flood to ebb;

stronger but shorter flood than ebb;

transition from ebb to flood longer than from flood to ebb; and

stronger but shorter ebb than flood.

The resulting time series of bottom friction velocity that are shown in Fig. 9 reflect the differences in flood and ebb length and velocity maximum, respectively. The additional asymmetry introduced by the M_{4} tidal forcing has major impacts on tidally averaged velocity profiles and sediment transports. The major factor in determining the tidally averaged velocity profiles seems to be the maximum bed stress of the flood relative to the ebb (Fig. 10). For unlimited sediment supply from the bed , the total sediment flux is dominated by the fluctuation transport because of erosion asymmetry between ebb and flood, which is subject to the asymmetric bottom shear stress (Fig. 9). The stronger and shorter flood (Fig. 11f) results in a strong up-estuary fluctuation flux, and the stronger and shorter ebb (Fig. 11l) results in a strong down-estuary fluctuation flux. In contrast to that, similar ebb and flood intensities (Figs. 11c,i) result in much weaker sediment fluxes, and are comparable to the pure M_{2} case shown in Fig. 5g.

Differences in advective sediment fluxes are best studied for the scenarios with limited bed sediment supply, where those are of comparable magnitude with the fluctuation fluxes. The shorter flood allows for less flood straining and therefore less convectively driven vertical mixing with the consequence of smaller tidally averaged eddy viscosity and eddy diffusivity. Thus exchange flows are more pronounced for this case (Fig. 10b), and tidal mean vertical sediment stratification is stronger than for equally long ebb and flood (not shown). The longer flood case generates the strongest tidal mean vertical mixing and thus the smallest exchange flow (Fig. 10d). This impact on the tidally averaged velocity and sediment profiles results in an advective up-estuary tidal mean sediment transport , which is twice as large in spatial amplitude for the short flood as compared to the case with the short ebb , compare Fig. 11d with Fig. 11j and Fig. 11e with Fig. 11k.

The total sediment flux for the full Ro–Si parameter space is shown in Fig. 12. In the case of unlimited sediment supply, the sediment transport for the more intense flood (Fig. 12f) is up estuary except for low Rouse numbers where sediments are so uniformly distributed in the vertical that the residual runoff dominates. In contrast, sediment flux is down estuary for all Rouse numbers when the ebb is stronger than the flood (Fig. 12l). The picture is slightly more complex for variations in relative lengths between transitions from ebb to flood and transitions from flood to ebb (Figs. 12c and 12i). At low Rouse numbers (Ro < 0.4) the short transition from flood to ebb generates significant down-estuary sediment fluxes whereas the short transition from ebb to flood generates some up-estuary sediment transport, as already seen in Fig. 11. This can be explained by the fact that the shorter transition between ebb and flood results (because of limited settling velocity) in a higher amount of sediment in the water column at slack after ebb, which is then transported upstream by the new bottom current, which has already reversed toward the up-estuary direction. This needs the turbulence near the bottom to decay faster than the sediment concentration decreases, which is given for relatively low Rouse numbers. This also explains why for Rouse numbers higher than 0.4 in the two scenarios with different lengths of transitions from ebb to flood behave in a similar way, but behave significantly different for low Rouse numbers. Despite these differences, both cases with a similar ebb and flood intensity qualitatively show a picture that is comparable to the pure M_{2} case shown in Fig. 6g: low Si or low Ro numbers represent down-estuary sediment fluxes, and sediment fluxes become more up-estuary oriented when increasing any of the two numbers.

Similar trends as discussed for the unlimited sediment supply are visible for the limited sediment supply scenarios, however to a much weaker extent, because of the decreased importance of the fluctuation flux. Also here, a stronger flood favors up-estuary fluxes and a strong ebb favors down-estuary fluxes.

## 6. Analysis of field data

The theory presented above including the significance of fluctuating sediment fluxes is approximately verified in a single spot of the parameter space considered here. For this purpose, observations in the tidally energetic Wadden Sea of the southeastern North Sea carried out by Becherer et al. (2011) are analyzed. The measurements were made in the Lister Deep on 15 May 2008, a curved tidal channel of up to 40-m depth connecting the Sylt–Rømø Bight in the northern part of the Wadden Sea with the North Sea. The dominant semidiurnal tide had a velocity amplitude of *U*_{2} = 0.7 m s^{−1}, and as a result of curvature effects of the tidal channel the residual velocity at the site of the measurements with a water depth of *H* = 14.8 m was directed upstream with *U _{r}* = 0.14 m s

^{−1}, leading to a dominance of the flood currents. Because of differential heating in spring and net precipitation and runoff from land, a horizontal buoyancy gradient of approximately ∂

_{x}

*b*= 1.4 × 10

^{−6}s

^{−1}was observed. Figure 13 shows observed sediment concentration profiles. Flood profiles (boldface lines) show a stronger vertical homogeneity than ebb profiles (thin lines) whereas the strongest vertical gradients in sediment concentration are visible for the slack tide profiles (dashed lines) when turbulent mixing reached a minimum. Based on a comparison to simulated sediment profiles (Fig. 14), the settling velocity of the sediment is estimated as

*w*= 2 × 10

_{s}^{−4}m s

^{−1}. The bed roughness length has been estimated as . With a friction velocity scale of = 0.028 s

^{−1}[note that this value deviates from the value presented by Becherer et al. (2011), because there the maximum friction velocity was used, whereas here we use the rms friction velocity], this leads to the following nondimensional parameters governing the dynamics: Si = 0.38, Un = 0.07, Ro = 0.02, , and (see also Table 1). For more details of the observations, see Becherer et al. (2011).

The observations are made nondimensional in the same way as presented in section 2 and were related to the relative tidal phase . For the sediment scaling, it was assumed that fine sediments were not able to settle out of the water column in this tidally energetic channel dominated by migrating sand waves (Becherer et al. 2011). To eliminate effects of lateral sediment advection that led to a fluctuation of the sediment contents in the water column by about a factor of two during the tidal cycle (Fig. 13), all sediment profiles are scaled to give a vertical mean of unity. The results of , , and are shown in Figs. 14a–c. Model simulations using the same nondimensional parameters are presented in Figs. 14d–f.

Because of the upstream-directed residual current of the considered scenario, the flood velocity is larger than the ebb velocity. The tidal straining leads furthermore to flood velocity profiles that vertically are more homogeneous than the ebb velocities. This effect is visualized by the black lines in Figs. 14a and 14d marking the lowest position in the water column where 75% of the maximum current speed in the water column is reached. This position is significantly lower during flood than ebb, for both the observations and the model simulations. Tidal straining is also visible in the observed and simulated buoyancy profiles (Figs. 14b,e). During and after full flood buoyancy tends to be well-mixed or even unstably stratified (less buoyant water overlaying more buoyant water), whereas during ebb a background stratification remains [for the observations, see also Becherer et al. (2011)]. This effect is more pronounced in the observations indicating that other processes such as lateral density straining (Lerczak and Geyer 2004; Scully et al. 2009; Burchard et al. 2011) is supporting restratification processes here. This is specifically evident during the slack tide after flood , where the observations show significant stratification, in contrast to the numerical simulations which ignore lateral straining effects.

Although the scaled tidal phase–averaged observed sediment data (Fig. 14f) still show significant sediment patchiness, it is evident that vertically the profiles are almost homogeneous during flood in the bulk of the water column and significantly more stratified during ebb. This difference between sediment profiles during flood and ebb is in agreement with the model simulations shown in Fig. 14f. During flood (*u*′ > 0), *c*′ < 0 is observed near the bottom and *c*′ > 0 is observed near the surface. During ebb (*u*′ < 0), the opposite is observed, with *c*′ > 0 near the bottom and *c*′ < 0 near the surface. Averaged over the tidal cycle, *u*′*c*′ < 0 near the bottom and *u*′*c*′ > 0 near the surface, a feature that is visible for the fluctuation sediment flux for example in Figs. 6f and 6g.

## 7. Conclusions

The results of this idealized numerical modeling study show that estuarine circulation is generally not the most important process driving sediment transport in tidally energetic (e.g., weakly stratified or periodically stratified) estuaries. Only for the simple case of a vanishing sediment bottom pool and symmetric tidal forcing (no M_{4} contribution to the external tidal forcing) transport due to estuarine circulation (here referred to as transport flux) was comparable to the sediment flux due to the covariance between current velocity and sediment concentration [here referred to as fluctuation flux, which is the same as the tidal pumping flux used by Scully and Friedrichs (2007b)], see Fig. 8. Whenever a bottom sediment pool is present or a significant M_{4} component in tidal forcing is present, the fluctuation sediment flux dominates the total residual sediment flux in the estuary. This result is contrary to the conventional view of the general dominance of transport sediment flux as supported by the analytical solution (ignoring tidal fluctuations) laid out in section 4.

This result could be confirmed for a large parameter space spanned by the Rouse (settling velocity relative to the rms friction velocity) and Simpson (horizontal density gradient relative to the stress divergence scale) numbers. Interestingly, the pattern of up- and down-estuary sediment flux as a function of Ro and Si generally follows the same pattern as predicted by a simple steady-state model with parabolic eddy viscosity and eddy diffusivity: for low Ro and Si, the sediment flux is down estuary (because of the weak tidal asymmetry and weak estuarine circulation and the down-estuary depth-mean flow owing to river runoff). Increasing Si increases the up-estuary components of the sediment flux, and the higher the Rouse number, the stronger the vertical sediment gradients and, as a consequence, the stronger the fluctuation flux (because of increased sediment concentration tidal fluctuations *c*′). Additionally, relatively high tidal mean sediment concentrations as a result of higher settling velocity result in differential advection of sediment owing to estuarine circulation (with a net up-estuary vertical integral).

Even for the case of a significant M_{4} component in the tidal forcing, this general picture (down-estuary sediment flux for low Ro and Si and up-estuary flux for high Ro and Si) is visible in most cases. For large bottom sediment pools, the sediment flux is generally up estuary for a stronger but shorter flood than ebb and vice versa. In contrast, in absence of a bottom sediment pool, total sediment fluxes are not highly sensitive to the M_{4} tidal component, because a large part of the sediment flux is due to the transport flux where the two factors residual velocity profile and tidal mean sediment profile are relatively insensitive to the M_{4} tidal component.

Although the general dependence of estuarine sediment fluxes on Si and Ro appears to be intuitive, the composition of sediment fluxes is complex, as seen, for example, in Figs. 6 and 11. While the vertical integrals of the sediment flux are generally similar for the total sediment flux (red line) and the fluctuation sediment flux (blue line), the shape of the profiles may vary significantly. For a small pool of erodible sediments, fluctuation flux and total flux may even have a different sign near the bed (Figs. 6f,g). This effect, however, is not present any more for tidal asymmetries as a result of M_{4} tidal forcing (Fig. 11).

This study is limited to weakly or periodically stratified estuaries. For mostly permanently stratified estuaries such as for the Columbia River (Jay and Musiak 1994), the Hudson River (Geyer et al. 2001), the Elbe River (Kappenberg et al. 1995), and the Weser River (Wellershaus 1981), it can be assumed that residual sediment transport due estuarine circulation plays a more important role, for two reasons: (i) the estuarine circulation is stronger and (ii) the sediment is largely trapped in the bottom boundary layer capped by the permanent pycnocline. An investigation of processes driving residual sediment transports in permanently stratified estuaries requires at least two- (longitudinal–vertical) or three-dimensional models, making such a study more demanding.

This present study ignores various processes that have the potential to even increase the importance of the fluctuation sediment flux. For cohesive sediments, flocculation processes (smaller flocs and thus smaller settling velocity during flood because of higher turbulence and smaller Kolmogorov micro scale) are suspected to quantitatively influence residual sediment fluxes in estuaries (Scully and Friedrichs 2007b). This process may be complicated because sediment stickiness strongly depends on biogeochemical factors such as exopolymer concentration (Schartau et al. 2007). In estuaries with high sediment concentrations, the sediment stratification itself has the potential to damp small-scale turbulence, with the potential to increase tidal sediment concentration asymmetries and thus the importance of the fluctuation flux (Winterwerp and van Kessel 2003). Also, the bed composition adds complexity to the generation of subtidal sediment fluxes in tidal estuaries. Sediments may be composed of mixtures of sand and mud with the consequence that the sand fraction results in retention of mud particles (van Kessel et al. 2011). Also, the presence of benthic biota may change the erodibility of sediments significantly (Borsje et al. 2008).

Given this complex dependency of estuarine sediment fluxes on physical and biogeochemical processes that are not well understood, it may be questioned if a predictive numerical model of estuarine sediment fluxes can be constructed. Emerging societal questions as for the susceptibility of estuarine ecosystems to global climate change may be hard to answer, given the sensitivity of many sedimentary processes on climate parameters such as temperature.

It should finally be noted that the present study is based on water column models, which may represent one certain location in an estuary. In real estuaries where the governing nondimensional parameters may spatially vary over a wide range, the total subtidal sediment flux is composed of a complex interplay of local processes.

## Acknowledgments

This study was partially carried out during a sabbatical visit of Hans Burchard at the Woods Hole Oceanographic Institution. Project funding was provided by the German Research Foundation (DFG) in the framework of the Project ECOWS (Role of Estuarine Circulation for Transport of Suspended Particulate Matter in the Wadden Sea, BU 1199/11) and by the German Federal Ministry of Research and Education in the framework of the Project PACE [The future of the Wadden Sea sediment fluxes: still keeping pace with sea level rise? (FKZ 03F0634A)].

### APPENDIX

#### Bottom Sediment Flux and Erosion Time Scale

For the sediment concentration, the classical sedimentation–erosion formulation (Krone 1962) is used at the bottom:

with the nondimensional sedimentation and erosion fluxes, and , respectively,

the nondimensional bottom friction velocity , and the critical friction velocity , which may be a function of the erosion depth , to reproduce effects of sediment erodibility caused by, for example, compaction or biological processes (Cartwright et al. 2009; Dickhudt et al. 2009). Note that while we ignore hydrodynamic changes in water depth caused by morphodynamic processes, the erosion depth might be essential to be tracked (Sheng and Villaret 1989). Furthermore, is the constant in time nondimensional erosion parameter. For the critical friction velocity, we use here

with the maximum depth of erodibility . In practical terms, this critical friction velocity criterion is applied such that a bottom pool of erodible sediment is tracked, and the erosion ceases when *B* approaches zero:

with the nondimensional bottom sediment concentration . Generalized forms of (A4) may be derived for variable bottom sediment concentrations and smooth (and potentially time dependent) critical friction velocity profiles as a function of the actual or average bottom friction velocity. Effectively, a dynamic equation for the sediment pool is calculated:

The variations of between the different numerical experiments could be obtained by varying the dimensional value of *α _{e}* or by varying the initial sediment concentration in the water column (which has a direct effect on the resulting value of ). Numerical test calculations showed that both ways give the same nondimensional results. For stationary solutions, will result, and for all periodic solutions, will be obtained. In the case of unlimited sediment supply and zero critical shear stress, holds, because the nondimensional tidal mean stress equals unity. As shown in (A6), the inverse of gives a nondimensional erosion time scale that can be normalized to give the fraction of a tidal period needed to erode the bottom sediment pool for the case that the average amount of water column sediment is in the bottom pool, that is, :

where *T _{e}* is the dimensional and is the nondimensional (normalized by the tidal period 2

*π*/

*ω*) time scale for erosion. With this, is the nondimensional time scale needed to erode a typical amount of sediments from the bottom pool into the water column, if initially all sediments would rest in the bottom pool. For , it would be expected that the bottom pool would be emptied during the tidal cycle.

## REFERENCES

*Coastal and Estuarine Environments: Sedimentology, Geomorphology and Geoarchaeology,*K. Pye and J. R. L. Allen, Eds., Geological Society of London, 13–30.

*Geophys. Res. Lett.,*

**38,**L17611, doi:10.1029/2011GL049005.

*Ocean Dyn.,*

**60,**1219–1241, doi:10.1007/s10236-010-0329-8.

_{2}and M

_{4}tide on the northwest European continental shelf

*Contemporary Issues in Estuarine Physics,*A. Valle-Levinson, Ed., Cambridge University Press, 27–61.

*Physics of Estuaries and Coastal Seas,*J. Dronkers and M. B. A. M. Scheffers, Eds., Balkema, 315–327.

*Ocean Coastal Manage.,*

*J. Geophys. Res.,*

**112,**C07028, doi:10.1029/2006JC003784.

**6,**27–31.

*Cont. Shelf Res.,*

**31,**S124–S134.