Influence of asymmetric tidal mixing on sediment dynamics in a partially mixed estuary

1. Introduction

Estuarine turbulence is characterized by high rates of kinetic energy dissipation rate and strong vertical density stratification or buoyancy frequency with both extending to much higher values than in other natural environments (Geyer et al., 2008). The asymmetric variations of hydrodynamics and sediment transport during flood and ebb periods are ubiquitous in estuaries. Mechanistic understanding of the asymmetries of hydrodynamics and sediment dynamics is important because the characteristics of sediment distributions and transport have a profound influence on harbor silting (Jiang et al., 2013), pollutant transport (Tessier et al., 2011), as well as on geomorphological evolution (Luan et al., 2017).

Asymmetry of hydrodynamics during a tidal period includes several aspects, such as peak current asymmetry, slack water asymmetry, and asymmetric tidal mixing (ATM). The latter is controlled by tidal straining, which is the semidiurnal variation in thermohaline stratification induced by the interaction of the vertically sheared tidal currents with the along-estuary salinity gradient (Simpson et al., 1990). Field measurements in many partially mixed estuaries demonstrated that tidal straining typically promotes vertical mixing during the flood tide and stratification during ebb, such as in the Liverpool Bay (Simpson et al., 1990), Columbia River Estuary (Jay and Musiak, 1994), Hudson River Estuary (Scully and Geyer, 2012; Geyer and MacCready, 2014) and York River Estuary (Scully and Friedrichs, 2007). Similar results for semidiurnal variations of mixing and stratification were obtained from a numerical model of a partially mixed estuary (Cheng et al., 2011). However, reverse asymmetry of tidal mixing with enhanced stratification during flood and intensified mixing during the ebb phase have also been examined (Fugate et al., 2007). Similar observations have been made in estuaries with a side embayment freshwater source (Lacy and Monismith, 2001; Fram et al., 2007) or with a significant variation of lateral bathymetry (Cheng and Wilson, 2008; Li et al., 2015).

The residual flow induced by ATM could be solved analytically (Jiang et al., 2013) or decomposed from numerical results (Chen and De Swart, 2018; Cheng et al., 2019). The strength and vertical structure of ATM-induced flow are associated to distinct estuarine mixing types. According to the numerical study of Cheng et al. (2011, 2013), the ATM induced flow has a two layer structure similar to that of density driven flow in weakly stratified estuaries. In partially mixed and highly stratified estuaries, it tends to have a three-layer structure of residual flow, with landward flows near the surface and bottom, and seaward flow in the middle of the water column (Cheng et al., 2011). In addition, the strength of ATM induced flow is positively correlated to the degree of vertical mixing (Cheng et al., 2013). These arguments are consistent with observations (Stacey et al., 2010; Cheng et al., 2020) and results from analytical models (Jiang et al., 2013).

The variation of vertical mixing and stratification between flood and ebb periods not only affects vertical residual circulation, but also the suspended sediment load (Scully and Friedrichs, 2007; Dijkstra et al., 2017). Tidal straining of the horizontal density gradient causes a stronger pycnocline lower in the water column, which acts as a barrier for sediment resuspension and reduces eddy diffusivity by limiting turbulent length scales during the ebb period (Geyer and MacCready, 2014; Li et al., 2018). This process is accounted for by modifying the value of eddy diffusivity based on the gradient Richardson number (Friedrichs and Hamrick, 1996). The observation in the York River Estuary showed that, higher values of eddy viscosity would occur during the less stratified flood period due to tidal straining. As a result, more sediment become re-suspended during the flood period and the net sediment transport is directed landward, despite of a net seaward advective flux with approximately equal current magnitude and shear (Scully and Friedrichs, 2003). The reverse asymmetry of sediment diffusivity is also observed in some other estuaries, e.g., the Changjiang River Estuary, where the stratification was strengthened during the flood tide and weakened during the ebb tide, with the control of the advection of the salt wedge (Li et al., 2015). The asymmetry of diffusivity induced suspended sediment transport, which is a part of tidal pumping transport, is likely important for the fine sediment movement in partially mixed estuaries (Uncles et al., 2006; Burchard et al., 2018). For silt-sized sediment particles, the numerical results of Geyer et al. (1997) indicated a 20-fold increase in the sediment transportation rate with inclusion of the mixing/stratification effect on turbulent sediment diffusion. With a semi-analytical model (iFlow model), Dijkstra et al. (2017, 2019) showed that the time-varying eddy viscosity and sediment diffusivity may have a similar important contribution on both tidal flow (traditional ATM) and gravitational circulation (eddy viscosity-shear covariance, or ESCO).

In summary, the residual flow induced by asymmetric eddy viscosity has been widely observed and studied using numerical and analytical models. Observations also indicated that the time varying diffusivity may be important for the sediment distribution and transport. However, accurate simulation of suspended sediment concentration (SSC) are subject to limitations, such as temporal variations of bed roughness, settling velocity, eddy diffusivity, which constrain mechanistic understanding of the effects of ATM on sediment dynamics.

In this paper, a vertical one-dimensional idealized model is developed to examine the influence of ATM on SSC and sediment transport in estuaries. Although the rigorous model assumptions are challenged by the complex estuarine hydrodynamics and sediment dynamics, the idealized model could be solved analytically and facilitates mechanistic understanding by isolating the individual mechanisms. In this study, analytical solutions are obtained for the main component of tidal flow (M₂) and for the residual flow, as well as for the mean and tidally varying components of SSC. The main objectives of this study are: (1) to examine the vertical structure of ATM-induced SSC; (2) to analyze the time lag between ATM-induced mixing and corresponding changes in the vertical distribution of SSC; (3) to discuss the importance of ATM-induced changes in SSC on sediment transport. The paper is structured as follows: the model and the solution technique are introduced in Section 2; in Section 3, we apply the model to a specific case study (the York River Estuary) and compare model results with observations taken from literature; further analyses of model results are presented in Section 4, followed by a discussion in Section 5 and conclusions in Section 6.

2. Analytical model

2.1 Hydrodynamics

We describe the vertical (z, unit: m) distribution of the along-estuary (x, unit: m) velocity (u(z)) using the one-dimensional (vertical) shallow water equation in an Cartesian coordinate system (x, z):

where the x-axis points from land to sea and the z-axis points vertically upward, with z = 0 at the undisturbed water surface (Fig. 1). u is the along-estuary velocity (unit: m/s), t is time (unit: s), $ g $≈9.81 m/s² is gravitational acceleration, $ {\rho }_{0} $≈1 020 kg/m³ is the reference water density and $ \partial \eta /\partial x $ represents the water level gradient (unit: 1). The along-estuary density gradient ($ \partial \rho /\partial x $) is prescribed and assumed to be constant (Chernetsky et al., 2010). The vertical eddy viscosity coefficient A (unit: m²/s) is assumed to be constant over depth and varies during the tidal cycle as a cosine function (Cheng et al., 2010),

Figure  1.  Site map of the York River Estuary showing the location of field observations. A sketch map of the model geometry attached in the top-right corner, where $ z=\eta $ is the elevation of the water surface and $ z=-H $ is the bed. The horizontal dashed line indicates the average sea level (z = 0), while the solid line shows surface elevation during a tidal period of length P.

DownLoad: Full-Size Img PowerPoint

where $ {A}_{0} $ represents a constant (tidal mean) part and $ \left|{A}_{2}\right| $ is the amplitude of the tidal varying part, which is proportional to and an order smaller than the mean part, ($ \left|{A}_{2}\right|=\alpha {A}_{0} $). ω (unit: s⁻¹) is the M₂ tidal frequency and $ \theta $ (unit: 1) is the phase angle of $ {A}_{2} $. It should be noted that eddy viscosity is affected by the magnitude of the tidal velocity and by the vertical density stratification. Under the influence of the former, the eddy viscosity would vary at the M₄ frequency, since it is high during both flood maxima and ebb maxima, hence the frequency doubled. However, many observations and numerical studies found a typical semi-diurnal variation (M₂ frequency) of eddy viscosity over a tidal cycle, which results from the influence of density stratification (Stacey et al., 2008; Cheng et al., 2020). Therefore, we only considered variations at the M₂ frequency in this model. The mean part of vertical eddy viscosity coefficient ($ {A}_{0} $) is parameterized using the formulation of Munk and Anderson (1948) with $ {A}_{0}= {c}_{V} U\left(H/2\right){(1+10Ri)}^{-1/2} $. In this expression, $ {c}_{V} $ is an empirical coefficient, $ Ri=g\left(\Delta \rho /\rho \right)\left(H/2\right)/{U}^{2} $ is the depth-averaged bulk Richardson number, $ H $ is the water depth, $ \Delta \rho $ is the residual bed to surface density difference and U is a typical value of the current velocity, which is assumed to be equal to the amplitude of the tidal velocity. An external semidiurnal tidal discharge $ {q}_{2} $ (unit: m²/s) and a constant river discharge $ {q}_{0} $ are imposed,

The imposed semidiurnal tidal discharge $ {q}_{2}={q}_{2a}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega t\right) $ has an amplitude ($ {q}_{2a} $) and with the M₂ tidal frequency ($ \omega $). The kinematic boundary conditions at the sea surface is applied,

At the riverbed $ z=-H $, partial slip and impermeability condition are imposed (Chernetsky et al., 2010),

where the parameter s is the constant partial slip parameter.

2.2 Sediment dynamics

SSC (c, unit: kg/m³) is described by the one-dimensional (vertical) advection-diffusivity equation,

The settling velocity ($ {w}_{s} $) is assumed to be constant, as its variation over a tidal cycle is small for estuaries with relative low SSC (Manning and Bass, 2006). The sediment vertical eddy diffusivity (K, unit: m²/s) varies during the tidal cycle with the same scaling coefficient as eddy viscosity ($ {K}_{2}=\alpha \left|{K}_{0}\right| $) and the same phase value $ \theta $. Here we assume a same proportional coefficient ($ \alpha $) and phase lag value ($ \theta $) for eddy diffusivity (K) and eddy viscosity (A), since both are influenced by the tidal varying structure of mixing and stratification. Strong (weak) stratification during ebb (flood) indicates both low (high) eddy viscosity and eddy diffusivity. Hence,

with $ {K}_{0}={c}_{V}{U}_{2}\left(H/2\right){(1+3.33Ri)}^{-1/2} $ denoting the mean component of sediment vertical eddy diffusivity. At the water surface, the diffusive and settling sediment flux balance,

where both the diffusive flux $\left( K\dfrac{\partial c}{\partial z}\right) $ and settling flux ($ {w}_{s}c $) are positive since we use a positive value of settling velocity, although it should be negative in our coordinate system. A normal flux condition is imposed at the bottom,

where $ {c}_{a} $ is the dimensionless reference sediment concentration and it is parameterized as (Dyer, 1997)

In this expression, ${\rho }_{s}\approx 2\;650\;\mathrm{kg}/ {\mathrm{m}}^{3}$ is the sediment density, $p\approx 0.4$ is the porosity of the sediment, a indicates the amount of sediment at the bed available for resuspension, $ {\tau }_{c} $ is the threshold bed shear stress, and $ {\tau }_{b} $ is the bed shear stress, which defined as

where $ |\cdot | $ denotes the absolute value.

2.3 Analytical solutions

The governing equations of flow and SSC are solved analytically with perturbation techniques, in which the physical components are divided into several order of magnitudes (Huijts et al., 2006). During perturbation analysis, the governing equations are non-dimensionalized firstly. Then we introduced a small parameter ($ \varepsilon \ll 1 $) to indicate the order of each terms in the governing equations. In this process, we assumed that the values of some terms are significant smaller (or an order smaller) than other terms (e.g., $ \varepsilon =\dfrac{{u}_{0}}{{u}_{2}}\ll 1 $ means the residual flow $ {u}_{0} $ is an order lower than the tidal flow $ {u}_{2} $). All model assumptions are summarized in Table 1. They have been used in related studies, and are validated using observations for typical tidal estuaries (Huijts et al., 2006; Jiang et al., 2013).

Table  1.  Assumptions of the one-dimensional analytical model

Assumption	Explanation
$ \dfrac{\eta }{H}=O\left(\varepsilon\right) $	water level fluctuation is an order of magnitude smaller than water depth
$ \left(\dfrac{U_0}{U_2},\dfrac{U_d}{U_2},\dfrac{U_q}{U_2}\right)=O\left(\varepsilon\right) $	residual currents are an order of magnitude smaller than tidal currents
$ \left(\dfrac{A}{\omega H^2},\dfrac{K}{\omega H^2}\right)=O\left(1\right) $	frictional force can affect the whole water column
$ \left(\dfrac{A_2}{A_0},\dfrac{K_2}{K_0}\right)=O\left(\varepsilon\right) $	the fluctuating component of eddy viscosity (diffusivity) is an order of magnitude smaller than the mean component
$ \dfrac{w_s}{\omega H}=O\left(1\right) $	length scale of sediment settling during a tidal cycle is comparable to water depth
Note: Here, $ \left({U}_{0},{U}_{2}\right) $ represent typical values of the mean velocity and M2 tidal velocity, respectively. The parameters ${U}_{d}=g{H}^{3} \Big(\dfrac{\partial \rho }{\partial x}\Big)/ \left(48\rho A\right)$ and $ {U}_{q}=2q/\left(3H\right) $ are the typical velocity scales for residual currents driven by an along-estuary density gradient and fresh water discharge. The parameter $ \varepsilon \ll 1 $ indicates the magnitude of the difference to higher-order physical components.

 | Show Table

DownLoad: CSV

With scaling analysis and perturbation analysis, the governing equations and boundary conditions can be decomposed into separate equations for tidal flow ($ {u}_{2} $), residual flow ($ {u}_{0}={u}_{0\rho }+ {u}_{0a}+{u}_{0q} $), tidally-averaged SSC ($ {c}_{0} $), M₂ SSC variations ($ {c}_{2}= {c}_{2t}+{c}_{2a} $) and M₄ SSC variations ($ {c}_{4} $), which can be solved analytically. Here, $ {c}_{2t} $ is driven by the rectification and bottom shear stress from residual flow, while $ {c}_{2a} $ depicts the part of sediment vertical diffusion induced by the time-varying diffusivity $\Big({\rm{i.e.}},\; {K}_{2}\dfrac{{\partial }^{2}{c}_{0}}{\partial {z}^{2}}\Big)$. At last, the vertical profile of sediment flux per unit width ($ F={F}_{0}+{F}_{2}={F}_{0\rho }+{F}_{0a}+{F}_{0q}+{F}_{2t}+{F}_{2a} $, unit: kg/(m²·s)) is obtained as the mean product $ \left(uc\right) $ over a tidal cycle (P). The sediment transport per unit width ($ T={T}_{0}+{T}_{2}={T}_{0\rho }+{T}_{0a}+ {T}_{0q}+{T}_{2t}+{T}_{2a} $, unit: kg/(m·s)) is obtained by integrating the vertical profile of sediment flux over depth. All symbols introduced above are described in Table 2. As for the subscripts used during this process, numbers indicate the fluctuation frequency, while letters (if exist) indicate the physical driving force. Therefore, $ \left({u}_{0},{u}_{2}\right) $ denote residual flow and M₂ tidal flow, while $ ({u}_{0a}, {u}_{0q}, {u}_{0\rho }) $ denote residual flow driven by asymmetric tidal mixing, by fresh water discharge, and by the along-estuary density gradient, respectively. A detailed description of the decomposition is described in Appendix A.

Table  2.  Definition of symbols used in the model description

Symbol	Definition	Symbol	Definition
$ \omega $	M2 tidal frequency	$ {w}_{s} $	settling velocity
$ P $	M2 tidal period	$ c $	SSC
$ Ri $	Richardson number	$ {c}_{0} $	tidally mean SSC (dominant order)
$ s $	partially slip parameter	$ {c}_{4} $	M4 component SSC (dominant order)
$ \partial {\rho }/\partial {x} $	along-estuary density gradient	$ {c}_{2} $	M2 component SSC (first order)
$ {{q}}_{0} $	fresh water discharge	$ {c}_{2t} $	bottom shear stress induced $ {c}_{2} $
$ {{q}}_{2} $	M2 tidal discharge	$ {c}_{2a} $	ATM induced $ {c}_{2} $
$ {u} $	along-estuary velocity	$ {c}_{a} $	reference concentration
$ {{u}}_{2} $	M2 tidal velocity (dominant order)	($ A $, $ K $)	eddy viscosity and eddy diffusivity
$ {{u}}_{0} $	tidally mean velocity (first order)	($ {A}_{0} $, $ {K}_{0} $)	tidally mean component of A and K
$ {{u}}_{0{a}} $	ATM induced u0	($ {A}_{2} $, $ {K}_{2} $)	M2 component of A and K
$ {{u}}_{0{\rho }} $	$ \partial \rho /\partial x $ induced $ {u}_{0} $	($ \left|{A}_{2}\right| $, $ \left|{K}_{2}\right| $)	amplitude of A2 and K2
$ {{u}}_{0{q}} $	fresh water discharge induced u0	$ \theta $	phase of A2 or K2
Note: ATM: asymmetric tidal mixing; SSC: suspended sediment concentration.

 | Show Table

DownLoad: CSV

3. Field data

To compare the analytical model to observations, we used field data from the York River Estuary, Virginia, reported in a study by Scully and Friedrichs (2003). The York River Estuary is a sub-estuary of the west Chesapeake Bay. Most of the estuary has a main channel with maximum depth reaching 10 m and is flanked by two shallow shoals. The typical tidal range is about 0.7 m, with a typical tidal current reaching 1.0 m/s during spring. The average river discharge is about 31 m³/s. The estuary is partially mixed, but tidal asymmetries in salinity stratification were found to play an important role in vertical momentum transfer and sediment suspension, and to have a great influence on sediment transport (Lin and Kuo, 2001). The tidal velocity and SSC were recorded with an acoustic Doppler current profiler and an optical backscatter sensor for about 25 h in April 1999. Detailed information of the observation can be found in Scully and Friedrichs (2003). The parameter values applied in the analytical model are either taken from Scully and Friedrichs (2003, 2007), or deduced from the study results of these two studies (Table 3). For example, the ratio $ \left|{K}_{2}\right|/{K}_{0} $ (i.e., $ \alpha =0.3 $) and the phase lag (i.e., $ \theta =1.4\pi $) are estimated from the variation of eddy viscosity depicted in the reference study (Scully and Friedrichs, 2003, 2007). In addition, the selected model input parameters (marked with asterisks in Table 3) were adjusted to obtain best agreement of the model results with observation.

Table  3.  Model input parameter values representative for a point in the York River Estuary, USA during April 1999 (referenced from Scully and Friedrichs (2003, 2007))

Quantity	Symbol	Value
Water depth	H	7.0 m
Amplitude of tidal velocity	U2	0.6 m/s
Depth mean velocity	U0	0.045 m/s
Angular M2 tidal frequency	$ \omega $	1.4 × 10−4 s−1
Gravitational acceleration	g	9.81 m/s2
Reference density	$ {\rho }_{0} $	1 020 kg/m3
*Bed to surface density difference	$ \Delta \rho $	2.0 kg/m3
*Along-channel density gradient	$ \partial \rho /\partial x $	8.6 × 10−4 kg/m4
*Threshold bed shear stress	$ {\tau }_{c} $	0.002 kg/(m·s2)
*Settling velocity	ws	0.6 mm/s
*Ratio of $ \left|{K}_{2}\right|/{K}_{0} $ or $ \left|{A}_{2}\right|/{A}_{0} $	$ \alpha $	0.3
*Phase lag of $ {K}_{2} $ or $ {A}_{2} $	$ \theta $	$ 1.4\pi $
Note: The stars mark the parameters deduced from the above references and were modified during model validation.

 | Show Table

DownLoad: CSV

4. Results

4.1 Model application to York River Estuary

4.1.1   Hydrodynamics

The simulations of the York River Estuary for the period of observations reported in Scully and Friedrichs (2003) covered two semi-diurnal tidal cycles (the observation lasted 25 h) (Fig. 2). The current velocity during flood period (negative values) and ebb period (positive values) were nearly equal since the semi-diurnal tidal current was the dominant flow (Fig. 2a). In reality, the ebb current was slightly stronger in the upper water column, while the flood current was slightly stronger in the lower part of the water column (Fig. 2b). For the mean flow, the up-estuary velocity (negative values) in the lower part of the water column was mainly contributed by the flow component driven by the longitudinal density gradient ($ {u}_{0\rho } $), while the downstream velocity (positive values) in the upper water column was mainly contributed by freshwater discharge ($ {u}_{0q} $). The component of the mean flow induced by ATM ($ {u}_{0a} $) had a nearly identical vertical structure as $ {u}_{0\rho } $, which is consistent with the analytical results of Jiang et al. (2013) for the North Passage of the Changjiang River Estuary with a same $ \theta $ value. According to the study of Cheng et al. (2010), the vertical structure of $ {u}_{0a} $ is a function of $ \theta $. However, the relative contribution of $ {u}_{0a} $ is not as significant as reported in former studies (Stacey et al., 2010; Burchard and Hetland, 2010; Jiang et al., 2013), since the strength of $ {u}_{0a} $ is highly correlated to the degree of mixing and stratification (Cheng et al., 2013). According to the study of Scully and Friedrichs (2003), the York River Estuary is partially mixed, the ratio of $ \left|{K}_{2}\right|/{K}_{0} $ ($\alpha \approx 0.3$) is relative small. The horizontal shift between observed mean flow and model results can be attributed to the ignored components of residual flow in our model, such as Stokes return flow, which drives a residual current to one direction over the whole depth.

Figure  2.  Model results of the vertical-temporal structure of total flow velocity (u = u₀+u₂) (a) and of the vertical structure of the residual components (b) based on the parameters of an observation point in the York River Estuary (Table 1). The blue line shows the flow component driven by the longitudinal density gradient ($ {u}_{0\rho } $), the green line is contribution of asymmetric tidal mixing ($ {u}_{0a} $), and the red line is the contribution from fresh water discharge ($ {u}_{0q} $). The solid black line shows the corresponding modeled mean velocity and the dashed line shows the observed mean velocity from Scully and Friedrichs (2003).

DownLoad: Full-Size Img PowerPoint

4.1.2   Sediment dynamics

The model results of vertical-temporal variation of SSC over two tidal periods showed that, SSC attained maximal values during both flood and ebb maxima in flow magnitude (Fig. 2). Hence, dominant variations in SSC occurred at a frequency of 2ω, which is mainly accounted for by the M₄ component SSC (c_4t in Fig. 3e). The M₂ component of SSC can be driven by the rectification and bottom shear stress from residual flow (c_2t in Fig. 3c) and advection of time-varying diffusivity with vertical shear of mean SSC (i.e., $ {K}_{2}\dfrac{{\partial }^{2}{c}_{0}}{\partial {z}^{2}} $, c_2a in Fig. 3d). The component c_2t was positive (negative) during flood (ebb) period since the near bottom residual flow was landward (Fig. 2b). The component c_2a was positive during flood slack period since the phase lag of $ {K}_{2} $ was $ 1.4\pi $ (i.e., highest diffusivity occurs during flood slack). Both observation and model results showed that, the near bottom SSC was significantly higher during flood period than that during ebb period (Figs 3a, b). Our model decomposition results reveal that the higher near bottom SSC during flood could mainly be attributed to c_2t. High values of c_2a were concentrated in middle layer, which can be also deduced from the vertical distribution of amplitude of each component (Fig. 3f).

Figure  3.  Model results (a) and observations (b, taken from Scully and Friedrichs (2003)) for the vertical-temporal structure of suspended sediment concentration (SSC) over two M₂ tidal periods (~ 25 h). The line plots in panel b compares temporal variation of near bottom SSC (z = 7 m) from observations and model results. The dashed horizontal lines in panel a mark the depth range with observational data shown in panel b. The labels “F” and “E” in panel a and b denote “flood” and “ebb”, respectively. Panels c, d and e show individual components of M₂ component SSC induced by bottom shear stress (c_2t), M₂ component SSC induced by asymmetry of diffusivity (c_2a), and M₄ component SSC induced by bottom shear stress (c_4t), respectively. Panel f shows the vertical distribution of SSC amplitude of each component.

DownLoad: Full-Size Img PowerPoint

4.1.3   Sediment flux

The vertical distribution of sediment flux is shown in Fig. 4. In correspondence to residual flows, the residual sediment flux induced by $ {u}_{0a} $ ($ {F}_{0a} $) and $ {u}_{0\rho } $ ($ {F}_{0\rho } $) is landward, while that induced by $ {u}_{0q} $ ($ {F}_{0q} $) is seaward. The tidal flux related to c_2t ($ {F}_{2t} $) is landward. With the similar explanation as for c_2a, the tidal flux related to c_2a ($ {F}_{2a} $) is seaward. Although the maximum value of c_2t (~ 40 g/m³, Fig. 3c) is significantly larger than that of c_2a (~ 28 g/m³, Fig. 3d), the corresponding sediment fluxes ($ {F}_{2t} $ and $ {F}_{2a} $) are of comparable magnitude, since that the higher values of c_2a mainly occurred in the middle or upper water column, where current velocity was high. As for the component F, both the observation and model results indicate the values of F is positive (seaward) in the upper water column and negative (landward) at lower depth. The positive seaward flux in the upper water column was mainly contributed by $ {F}_{2a} $, while the negative landward flux in the lower water column were contributed from mutual effect of $ {F}_{2a} $ and $ {F}_{0\rho } $.

Figure  4.  Vertical structure of sediment flux and its components based on the parameters for a point in the York River Estuary. $ {F}_{0a} $, $ {F}_{0\rho } $ and $ {F}_{0q} $ represent residual sediment flux related to asymmetry tidal eddy viscosity, horizontal density gradient, and river discharge. $ {F}_{2t} $ and $ {F}_{2a} $ represent tidal sediment flux related to bottom shear stress varying in M₂ frequency and asymmetry tidal diffusivity. The dashed black line (YR99) shows the observed sediment flux cited from Scully and Friedrichs (2003).

DownLoad: Full-Size Img PowerPoint

4.2 Effect of asymmetric tidal diffusivity on vertical sediment distribution

The model results in Section 4.1 showed that, the high values of ATM-induced variations in SSC (c_2a) is concentrated on the vertical middle depth with high tidal flow velocity, while the high values of SSC component induced by the near bottom shear stress of tidal flow (c_2t) occurs at the near bottom area with low tidal flow velocity. Thus the sediment flux of $ {F}_{2a} $ ($ {F}_{2a}=\left\langle{{u}_{2}{c}_{2a}}\right\rangle $) is considerable to that of $ {F}_{2t} $ (${F}_{{{2t}}}=\left\langle{{u}_{2}{c}_{{{2t}}}}\right\rangle$), even though the absolute value of ${c}_{{{2a}}}$ is smaller than that of ${c}_{{{2t}}}$ (Fig. 3). Many estuaries are well mixed during the flood period and stratified during the ebb period (typical ATM). The occurrence of reverse tidal asymmetry was attributed to a number of reasons, including a freshwater source from a side embayment (Fram et al., 2007), or lateral straining effect on deep-shoal system (Cheng et al., 2011). To consider both types of asymmetries, we varied the phase angle of K₂ ($ \theta $) between 0 and $ 2\pi $. We considered four specific cases in which $ \theta =[0, 0.5\pi ,\pi , 1.5\pi ] $ represent highest vertical diffusion coefficient occurs at the period of ebb maxima, ebb slack, flood maxima and flood slack, respectively (Fig. 5). The structure of c_2a shifts gradually in time when the value $ \theta $ varies from 0 to $ 2\pi $. A small phase lag (about 1 h) between maximum values of K₂ and highest concentration c_2a was caused by the time required for resuspension and vertical transport of sediment. In the following subsections, the vertical structure of SSC under the influence of ATM will be examined firstly in dependence on the phase lag. Following that, the relationship between the phase angle of K₂ and sediment flux will be discussed.

Figure  5.  Vertical-temporal structure of the asymmetric tidal mixing (ATM)-induced SSC (c_2a) during a single M₂ tidal period for different phase angles between tidal flow and sediment diffusivity ($ \theta $). Panels a, b, c and d represent $\theta =0,\theta =0.5\pi , $$ \theta =\pi \;\mathrm{a}\mathrm{n}\mathrm{d}\;\theta =1.5\pi$, respectively. The remaining parameters are shown in Table 1.

DownLoad: Full-Size Img PowerPoint

In a tidal estuary, the amplitude of SSC variation decreases with height above the bed (i.e., amplitude attenuation). In addition, the response of SSC to the tidal current is not instantaneous, and the SSC variation lags behind the velocity (i.e., phase lag $ \phi $) due to adjustment of the balance between upward diffusion and downward settling of sediment particles. In the study of Yu et al. (2011), the vertical distribution of the phase lag $ \phi \left(z\right) $ between variations in SSC and tidal flow was analytically derived as

In their study, this analytical solution assumed a constant value of K, i.e., the component c_2a resulting from ATM was neglected. Besides, the bottom SSC was prescribed rather than driven by bed shear stress. However, our model results showed that, the time-varying sediment diffusivity have a significant influence on sediment vertical distribution (Fig. 5). Furthermore, the variation of bed shear stress highly depends on the residual flow. Here, we discuss the vertical phase lag based on two residual flow scenarios, i.e., (i) $ \partial \rho /\partial x=0,{U}_{0}=0.10\;{\rm{m/s}} $ and (ii) $ \partial \rho /\partial x=1\times {10}^{-3}\;{\rm{kg}}/{{\rm{m}}}^{4},{U}_{0}=0 $. The first scenario represents a strongly stratified estuary, while the second scenario represents a fully mixed estuary. In these two scenarios, the near bottom shear stress induced by residual flow, as well as the c_2t component (induced by net bottom shear stress), are nearly equivalent. However, the maximum value of near bottom shear stress occurs at flood period in scenario (i), while it occurs at ebb period in scenario (ii). It will be shown that the influence of c_2a on these two scenarios is quite different.

In agreement with the solution proposed by Yu et al. (2011), the tidal-driven component SSC predicted by our model showed an increasing phase lag for increasing distance above the bed (Fig. 6a). Except for the near-surface region, there was good agreement between both our model results and Yu et al. (2011) (gray line compares to black line in Fig. 6). Taking the effect of the ATM-induced SSC (c_2a) into account, the vertical profile of the phase lag of c₂ (c₂ = c_2a + c_2t) changes dramatically in dependence on the phase of sediment diffusivity ($ \theta $).

Figure  6.  Vertical structure of phase lag of suspended sediment concentration in response to diffusivity phase $ \theta $ in case of: $\partial \rho /\partial x=0,{U}_{0}=0.15\;\mathrm{m}/{\mathrm{s}}$ (scenario z (i)) (a) and $\partial \rho /\partial x=1\times {10}^{-3}\;\mathrm{k}\mathrm{g}/{\mathrm{m}}^{4},{U}_{0}=0$ (scenario (ii)) (b). The black line represents the analytical result of Yu et al. (2011). The gray line shows the model results of c_2t, i.e., without consideration of asymmetric tidal mixing. The remaining 4 colored lines represent the result of $ \phi \left(z\right) $ with consideration of c_2a for four $ \theta $ values ($ 0, 0.5\pi ,\pi , 1.5\pi $).

DownLoad: Full-Size Img PowerPoint

In the first residual flow scenario ($\partial \rho /\partial x=0,{U}_{0}= 0.15\;\mathrm{m}/{\mathrm{s}}$), bottom shear stress is higher during the ebb period compared to the flood period, thus the c_2t component is positive (negative) during ebb (flood). Considering the influence of c_2a, if $ \theta =0 $, maximum diffusivities occurred during ebb, thus c_2a is positive during ebb (Fig. 5a). c_2a did not change the vertical profile of $ \phi \left(z\right) $ too much (red lines compared to gray line in Fig. 6a). However, the occurrence of the maximum value of c_2a was gradually delayed for increasing $ \theta $ (Fig. 5), which results in a gradually increasing phase lag of c₂. The positive values of c_2a occur during the flood period in case of $ \theta =\pi $ (Fig. 5c), because of the higher sediment diffusivity. As a result, the vertical distribution of the phase lag of c₂ differs far strongly from that of c_2t (cyan lines compared to gray line in Fig. 6a). The temporal and vertical distribution of the M₂ component of SSC (c₂) under four $ \theta $ values ($ 0, 0.5\pi ,\pi , 1.5\pi $) in case of ($ \partial \rho /\partial x=0,{U}_{0}=0.15\;\mathrm{m}/{\mathrm{s}} $) are shown in Fig. 7, where the vertical distribution of the phase lag can easily be observed.

Figure  7.  Vertical structure of M₂ suspended sediment concentration (SSC) component (c₂) over a tidal period under four different values of phase angle of tidally varying sediment diffusivity ($ \theta $) in case of scenario (i), $\partial \rho /\partial x=0,{U}_{0}=0.15\;{\rm{m}}/{\rm{s }}$. Panels a, b, c and d show results for $\theta =0,\;\theta =0.5\pi ,\; \theta =\pi \;\mathrm{a}\mathrm{n}\mathrm{d}\;\theta =1.5\pi$, respectively.

DownLoad: Full-Size Img PowerPoint

In the second residual flow scenario ($ \partial \rho /\partial x=1\times {10}^{-3}\;\mathrm{k}\mathrm{g}/{\mathrm{m}}^{4}, {U}_{0}=0 $), the near bottom shear stress and the c_2t component were positive (negative) during the flood (ebb) period, since the residual shear stress induced by the density gradient is directed landward, thereby enhancing the flood current. In this residual case, if $ \theta =\pi $ (maximum diffusivities occur at flood), the c_2a component would only increase the value of the M₂ SSC component, but not change its phase, because both c_a and c_2t are highest during the flood period. Therefore, the zero phase lag would occur at $ \theta =\pi $ (Fig. 6b, cyan line). The temporal and vertical distributions of the M₂ component of SSC (c₂) for four $ \theta $ values $ \left[0, 0.5\pi ,\pi , 1.5\pi \right] $ for scenario (ii) are shown in Fig. 8.

Figure  8.  Vertical structure of M₂ suspended sediment concentration (SSC) component (c₂) over a tidal period under four different values of phase angle of tidally varying sediment diffusivity ($ \theta $) in case of scenario (ii), $ \partial \rho /\partial x=1\times {10}^{-3}\; {\rm{kg}}/{{\rm{m}}}^{4},\;{U}_{0}=0 $. Panels a, b, c and d show results for $ \theta =0, \; \theta =0.5\pi ,\;\theta =\pi \;\mathrm{and}\;\theta =1.5\pi $, respectively.

DownLoad: Full-Size Img PowerPoint

4.3 Effect of tidally varying diffusivity on vertical sediment transport

In tidal estuaries, the residual flow plays an important role in the process of both sediment resuspension and transport. On the one hand, the residual sediment transport is induced by the superposition of residual flow and tidal mean SSC. On the other hand, asymmetry of sediment transport would occur between flood and ebb periods under the influence of residual flow, and thus tidal transport is induced. River discharge and the horizontal density gradient are two important factors for residual flow in estuaries. According to the classification based on salinity structure (Dyer, 1997), high fresh water discharge would make the estuary stratified (e.g., Silver Bay, Alaska), while a strong tide would result in fully mixed conditions (e.g., James River Estuary). In this study, the ratio of the typical residual velocity induced by density gradient and fresh water discharge (γ$ =\left[{u}_{0\rho }\right]/\left[{u}_{0q}\right] $) is used to quantify the degree of vertical mixing. A higher value of γ indicates a highly mixed estuary while a lower value of γ indicates a strongly stratified estuary. Here square brackets ($ \left[\cdot \right] $) indicate the velocity scale of residual flow, which is derived from the analytical solutions of Huijts et al. (2009) to estimate the relative importance of the various residual flow mechanisms from estuarine parameters. The velocity scale of density-gradient-driven residual flow ($ \left[{u}_{0\rho }\right] $) and discharge-driven residual flow ($ \left[{u}_{0q}\right] $) reads,

The variation of tidally mean bottom shear stress $ \left\langle{{\tau }_{b}}\right\rangle $ over $ \gamma $ and $ \theta $ is shown in Fig. 9a. Generally, the value $ \left\langle{{\tau }_{b}}\right\rangle $ is flood-ward (negative) in estuaries with larger along-estuary density gradients (large $ \gamma $), while $ \left\langle{{\tau }_{b}}\right\rangle $ is ebb-ward (positive) in the estuaries with a relative small values of $ \gamma $ (larger scale of river discharge induced flow). It should be noted that, the value of bottom shear stress should be always positive due to its definition, here the positive (or negative) sign indicates the direction of $ \left\langle{{\tau }_{b}}\right\rangle $, i.e., positive means seaward, negative means landward. Following Fig. 9a, the bottom shear stress ($ {\tau }_{b} $) is higher in case of density driven flow than that of river discharge driven flow with an equal velocity length scale, as the vertical length scale (z) of the former is much smaller (the solid white line is below the dashed white line in Fig. 9a). In addition, there is a slight variation of $ \left\langle{{\tau }_{b}}\right\rangle $ when $ \theta $ increases from 0 to $ 2\pi $ due to the ATM generated residual flow ($ {u}_{0a} $). The vertical structure of $ {u}_{0a} $ is similar to that of $ {u}_{0\rho } $ for $ \theta =\pi $, which enhances the flood current over the near bottom area, thus the shear stress $ \left\langle{{\tau }_{b}}\right\rangle $ attains its maximum value in flood direction under this condition. On the contrary, the near bottom ebb current will be enhanced by $ {u}_{0a} $ during $ \theta =0 $ (Cheng et al., 2010), and the shear stress $ \left\langle{{\tau }_{b}}\right\rangle $ attains its maximum value in ebb direction.

Figure  9.  Tidally mean near bottom shear stress $ \left\langle{{\tau }_{b}}\right\rangle $ (a) and sediment flux (b–i) as a function of $ \theta $ and $ \gamma =\left[{u}_{0\rho }\right]/\left[{u}_{0q}\right] $. The dashed white line indicates a ratio $ \gamma =1 $. The solid white line indicates the value $ \gamma $ for which the average shear stress $ \left\langle{{\tau }_{b}}\right\rangle=0 $. The solid red line in panel b–i indicates values of $ \gamma $ for which the sediment flux is zero. The remaining boundary conditions can be found in Table 1.

DownLoad: Full-Size Img PowerPoint

The sediment transport per unit width of the estuary (i.e., the sediment flux) and its components are shown in Figs 9b–i as a function of $ \gamma $ and $ {\theta } $. Generally, the total sediment transport (Fig. 9i) is directed upstream ($ {T}_{{\rm{ssc}}} < 0 $) in estuaries with a strong density gradient and a weak river discharge ($ \gamma > 1 $), in which the net bottom shear stress is landward ($ \left\langle{{\tau }_{b}}\right\rangle < 0 $). Contrary, ${T}_{{\rm{ssc}}} > 0$ in the estuaries for which $ \gamma < 1 $ and $ \left\langle{{\tau }_{b}}\right\rangle > 0 $. However, there is a transitional region for ${\theta }=0.5\pi \;{\rm{to}}\; 1.2\pi$, for which the total sediment transport is directed upstream (negative ${T}_{{\rm{ss}}c}$) despite a net seaward directed bottom shear stress (positive $ \left\langle{{\tau }_{b}}\right\rangle $). Furthermore, seaward sediment transport can be expected in estuaries with a net upstream directed bottom shear stress for $\theta = 1.2\pi -2.0\pi$ or $\theta =0- 0.5\pi$. The former situation was actually observed in the York River Estuary (Scully and Friedrichs, 2003).

The total sediment transport is composed of residual sediment flux ($ {T}_{0} $, Fig. 9g) and the M₂ tidal transport ($ {T}_{2} $, Fig. 9h). The zero residual flux line ($ {T}_{0}=0 $, red solid line) locates at upward (upward means the direction of higher $ \gamma $ value) of zero bottom shear stress line ($ \left\langle{{\tau }_{b}}\right\rangle=0 $, white solid line) in Fig. 9g. The main reason could be that, in the estuaries of $\gamma \approx 1$, although the net bottom shear stress generated by $ {u}_{0\rho } $ is stronger than that associated with $ {u}_{0q} $, the residual sediment flux induced by $ {u}_{0q} $ ($ {T}_{0q} $) is much higher than that induced by $ {u}_{0\rho } $ ($ {T}_{0\rho } $) (Figs 9c, d). Although the total vertical profile of $ {u}_{0q} $ is directed seaward, half of the vertical profile of $ {u}_{0\rho } $ is landward (Fig. 2b). The residual transport induced by $ {u}_{0a} $ is shown in Fig. 9b, which is relative weak since that $ {u}_{0a} $ is relative weak.

The M₂ tidal transport can be decomposed into two components, $ {T}_{2a} $ and $ {T}_{2t} $ due to that the M₂ component SSC. c₂ is affected by bottom shear stress and by ATM, i.e., $ {T}_{2}={T}_{2a}+{T}_{2t} $. The component $ {T}_{2t} $ is driven by the interaction of u₂ and c_2t, in which the component c_2t is induced by the asymmetry of bed shear stress, thus the zero line of $ {T}_{2t} $ coincide completely with the zero line of $ \left\langle{{\tau }_{b}}\right\rangle $ (Fig. 9f). The transport component $ {T}_{2a} $ is independent from bottom shear stress or from the scale of residual flow, and only a function of $ \theta $. The value of c_2a is positive (negative) during flood (ebb) period when $\theta =0.3\pi -1.3\pi$ (see Figs 5b and c for $ \theta =0.5\pi $ and $\pi $ for example), thus the transport component $ {T}_{2a} $ is landward. In contrary, the value of c_2a is positive (negative) during ebb (flood) period at the remaining values of $ \theta $ (see Figs 5a and d for $ \theta =0$ and $\theta =2\pi$ for example), thus the transport $ {T}_{2a} $ is seaward. In addition, although the typical value of c_2a is smaller than that of c_2t, the contribution of $ {T}_{2a} $ is significantly larger than that of $ {T}_{2t} $. This is because the high values of c_2a are concentrated on the middle of the water column, where the velocity is relative high, while high values of c_2t mainly occur at the near bottom region, where the velocity is relative small. Furthermore, comparing Figs 9e, f and g shows that the component $ {T}_{2a} $ is the key factor for allowing the solid red line and white line to intersect, i.e., the inconsistently direction of net bottom shear stress and total sediment transport.

5. Discussion

5.1 Effect of ATM on estuarine sediment dynamics

The vertical structure of tidal flow and SSC was modeled using a one-dimensional analytical model, with a focus on the asymmetric variation of the sediment diffusion coefficient during a tidal cycle. We analyzed the vertical profile of SSC in terms of the phase lag between temporal variations in SSC and tidal flow. With a one-dimensional analytical model, Yu et al. (2011) solved the vertical profile of the phase lag without consideration of ATM, and suggested that the phase lag of SSC increases linearly from the bottom to the surface (Fig. 6, dashed line). By considering the effects of ATM on sediment diffusion, our model showed that the vertical distribution of the phase lag is strongly influenced by ATM. In estuaries with typical ATM structure (higher diffusivity during flood and lower diffusion during ebb), and with residual flow that is dominated by fresh water discharge (e.g., during the flood season or neap tides), the near bottom SSC is be higher during ebb than that during flood due to higher near bottom shear stress. However, in the upper part of the water column, ATM causes an increased SSC during flood. This result is consistent with observations from the North Passage of the Changjiang River Estuary (Jiang et al., 2013). It has been found that the sediment transport associated with the M₂ tidal flow is landward in some areas, even though the fresh water discharge is quite strong during flood season, which may partly account for severe siltation in this area (Jiang et al., 2013; Wang et al., 2017). Similar observations were made in the York River Estuary (Scully and Friedrich, 2003), where the residual flow and near bottom shear stress was seaward due to river discharge, but SSC was higher during the flood period, bottom shear stress is reduced and the tidal-averaged transport is landward. These observation results are consistent with our model prediction (see Fig. 9h, seaward shear stress corresponds to landward sediment transport).

In estuaries with reverse ATM structure, i.e., higher eddy diffusion occurring during ebb, the sediment diffusion phase lag varies between $ 0 \;{\rm{and}}\; 0.5\mathrm{\pi } $. This kind of reverse ATM was observed in the two shallow shoals of upper Chesapeake Bay (Fugate et al., 2007), where the vertical distribution of SSC was more uniform during ebb than that during flood, indicating a higher diffusion during ebb and lower diffusion during ebb. The residual flow in the two shoals was seaward under the influence of river discharge according to Huijts et al. (2009). As a result, the SSC in the whole water column was much higher during ebb in response to higher diffusion and bottom shear stress, whereas much lower during flood in response to low diffusion and bottom shear stress. This is consistent with our model shown in Fig. 7a, in particular consistent with the modeled vertical distribution of the phase lag (Fig. 6a, $ \theta \approx 0 $). Similar reverse ATM structure and SSC variations were also observed in numerical simulations of an idealized estuary (Cheng et al., 2011) and in measurements at the South Channel in the Changjiang River Estuary (Yang et al., 2017).

5.2 Study limitations

The model in this paper focused on the vertical distribution of hydrodynamics and sediment dynamics, and provided mechanistic insight into the processes controlling sediment transport in estuaries and identification of their main environmental drivers. However, several factors, which might be important for the asymmetric tidal mixing of sediment diffusion, are not included in this model. First, and most important, our model is one dimensional, which neglects longitudinal and lateral variations and implies a straight channel shape of the estuary. Since we focused on the hydrodynamics and sediment dynamics and their vertical variation at a single location, longitudinal phenomena such as, turbidity maxima and salinity intrusions, could not be considered in this model. Furthermore, in a deep channel and shallow shoal system, the shoals may undergo a reverse tidal asymmetry process due to enhanced stratification by lateral straining during flood, while promoting vertical mixing during ebb (Cheng et al., 2011). Alternatively, the case of reverse tidal asymmetry was considered by the parameters, i.e., vertical eddy viscosity (A) and diffusivity (K). In reality, both parameters (A and K) are a function of salt transport (Cheng et al., 2010) and sediment distribution (Talke et al., 2009), and vary over depth (Chen and de Swart, 2018). Consideration of these dependencies prevents an analytical solution of the model; thus, we followed previous studies (e.g., Chernetsky et al., 2010) and considered them as fixed input parameters. Furthermore, the vertical velocity was ignored ($ w=0 $) due to the conservation law (no lateral process and no along-estuary variation of along-estuary velocity, i.e., $ v=0,\partial u/\partial x=0 $), which implies the vertical advection term ($ w\partial u/\partial z=0 $) was ignored, and the vertical movement of sediment only depends on diffusion and settling velocity. The vertical advection term was an order smaller than the inertia term ($ \partial u/\partial t $) according to the scaling analysis (Huijts et al., 2009). The effect of vertical velocity in sediment movement may result higher SSC during flood (upward tidal velocity) and lower SSC during ebb (downward tidal velocity). However, higher settling velocity can be expected due to flocculation process occurred during flood with higher salinity, while lower settling velocity may occur during ebb (Winterwerp, 2002; Dijkstra et al., 2019). These processes were represented by the asymmetry of sediment diffusivity in our model. At last, the sediment transport due to M₄ tidal flow was ignored, even though the M₄ component of SSC is one of the dominant parts of SSC, and M₄ tidal flow could be generated by nonlinear advection terms (Yang et al., 2014), or be amplified in a convergent channel (Friedrichs and Aubrey, 1994). The main reason for that is the lack of M₄ tidal discharge to solve for M₄ tidal flow. Fortunately, the M₄ tidal transport is relatively small as demonstrated in former studies (e.g., Jiang et al., 2013; Yang et al., 2014).

6. Conclusions

In this paper, we investigated the influence of time-varying diffusivity on sediment transport in a tidal estuary with an idealized analytical model. The one-dimensional (vertical) model considered M₂ tidal flow and a residual flow driven by along-estuary density gradient, fresh water discharge, and asymmetric tidal mixing. The suspended sediment concentration and sediment transport were balanced between settling and turbulent diffusion, and was decomposed into a tidal-averaged part and components varying with M₂ frequency and M₄ frequency. For the first time, our analytical model included SSC variations driven by time-varying diffusivity during a tidal cycle. As a case study, we applied the model to York River Estuary, where higher eddy diffusivity during flood and lower diffusivity during ebb were observed in a previous study. The model results agreed well with the observation. With the model, we analyzed the influence of asymmetric tidal mixing on sediment distribution and transport. The vertical distribution of the phase lag between SSC and tidal flow is strongly influenced by the time-varying diffusivity. The phase lag increases in the case of a typical ATM ($ \theta \approx \pi $) in an estuary with seaward tidally mean bottom shear stress, where the residual flow is dominated by river discharge. Contrary, the phase lag is reduced by ATM, if the tidally mean bottom shear stress is landward. Also sediment transport is highly dependent on the direction of bottom shear stress and ATM phase. In estuaries with typical ATM, higher diffusivity during flood results in enhanced landward tidal transport. Net landward sediment transport can be expected under these conditions, even if the net residual flow is seaward for relative high freshwater discharge (e.g., during flood season or neap tide period). On the contrary, in estuaries with relatively strong tides (e.g., during drought season or spring tides), reverse ATM may result in a net seaward sediment transport, even if the gravitational circulation induces a higher near bottom SSC during flood than that during ebb.