1 Introduction

Estuarine sediment dynamics are a consequence of the advective and dispersive transport of suspended sediment originating from freshwater and marine sources, and the bed-water exchange processes of erosion and deposition. These sediment fate and transport processes are a function of the hydrodynamic forcings involved, and the local bathymetry. The primary hydrodynamic forcings typically include barotropic circulation induced by the tide, estuarine circulation, and river flow, i.e., freshwater inflow from the head-of-tide (Dyer 1997). Additional, primarily episodic, forcings that can influence sediment dynamics include barotropic circulation induced by coastal setup and setdown, wind-induced circulation and mixing, local wind-waves, and ocean waves (Dyer 1997). The latter three are mainly relevant for relatively wide estuaries, estuaries subject to the propagation of ocean waves, and/or estuaries with large tidal shallows. Since the estuary reported upon in this paper is relatively narrow (i.e., fetch-limited for wind-waves) and is not directly influenced by ocean waves, the impacts of winds and waves are not examined further in this paper. Given their relevance to the analyses in this paper, the impact of the primary forcings (barotropic, estuarine, riverine) on sediment dynamics is briefly summarized first followed by a discussion of the present research.

1.1 Barotropic effects

Non-linear tidal distortion in shallow water (Dronkers and Schönfeld 1959; Friedrichs and Aubrey 1988) and the spatially-variable and oscillatory nature of estuarine tidal currents give rise to several residual sediment transport mechanisms, grouped into what is referred to as tidal asymmetries and lag effects, respectively. Dronkers (2005) and Gatto et al. (2017) provide a detailed review of the various transport mechanisms summarized here. These transport mechanisms, as described below with respect to their impact on residual sediment transport, assume that suspended sediment concentrations (SSC) scale over tidal time-scales.

The impact of tidal distortion on sediment transport is apparent in a Eulerian frame of reference and is driven by asymmetry in peak currents and asymmetry in slack-water duration. Residual (fine) sediment transport follows the asymmetry in peak currents, with flood-dominant systems (peak flood currents greater than peak ebb currents) exhibiting net up-estuary transport and vice versa for ebb-dominant systems (Dronkers 1986; van de Kreeke and Robaczewska 1993). Slack-water asymmetry originates from differences in the deceleration and acceleration of currents in the transition from flood to ebb and vice versa. Residual sediment transport follows the asymmetry in slack water duration, with longer slack water duration at the end of flood than ebb associated with net up-estuary transport, and vice versa for longer slack water duration at the end of ebb than flood. In addition to asymmetries generated by tidal distortion, the interaction between certain principal astronomic tidal constituents can also lead to the development of tidal asymmetries and residual sediment transport (Hoitink et al. 2003). Lag effects are conceptualized in a Lagrangian frame of reference and refer to residual sediment transport induced by the oscillatory and spatially variable nature of tidal currents in combination with certain sediment transport parameters, namely the critical shear stress for erosion, τcr, and a corresponding threshold for deposition (Postma 1954, 1961; Van Straaten and Kuenen 1957, 1958). Lag effects can further be categorized into settling lag and scour lag; both generally promote up-estuary residual transport. In addition to the periodic barotropic circulation induced by the tide, episodic sub-tidal variations in sea-level (induced by coastal setup and setdown) also result in volume fluxes into and out of estuaries (Salas-Monreal and Valle-Levinson 2008). However, the importance of such events on sediment transport dynamics is not straight-forward, depending on the interaction of these events with the other forcings. The use of the term barotropic in the remainder of this paper is limited to only the tide.

The barotropic effects on sediment dynamics are also collectively referred to as tidal pumping in the literature (Geyer et al. 2001; Scully and Friedrichs 2007; Sommerfield and Wong 2011). Tidal pumping has been shown by these authors to be a significant mechanism responsible for the up-estuary transport of sediments. Tidal pumping has also been shown to be dependent on river flow and the spring-neap cycle. Tidal pumping can influence the formation of estuarine turbidity maxima (ETM), a zone of elevated SSC and enhanced sediment trapping, in the tidal freshwater regions of an estuary due to the convergence of up-estuary barotropic transport and down-estuary fluvial transport (Burchard et al. 2018).

1.2 Estuarine circulation

Estuarine circulation includes several processes such as baroclinic (gravitational) circulation, shear induced by river flow, lateral and longitudinal advection, eddy viscosity-shear covariance, etc. (Dijkstra et al. 2017) resulting in residual near-bottom currents directed up-estuary and near-surface currents directed down-estuary. In combination with a vertical gradient in SSC (typically increasing from surface to bottom of water column), estuarine circulation induces an up-estuary flux of sediment that extends up to the limit of salinity intrusion and resulting in the formation of ETMs co-located with the salt front (Festa and Hansen 1978; Burchard and Baumert 1998; Sanford et al. 2001; Burchard et al. 2018).

1.3 Fluvial effects

In addition to delivering sediment to the estuary, freshwater inflow varying seasonally or episodically in response to rainfall runoff events can impact sediment dynamics within the estuary by pushing the salt-front seaward, enhancing bed shear stress, and potentially causing erosion and export of sediments from the estuary (Ralston et al. 2013). In addition, the direction of residual sediment transport can depend on river flow, with net up-estuary flux of sediments during low-flow conditions and net down-estuary flux during high-flow conditions (Geyer et al. 2001). The additional freshwater flow during runoff events can also impact tidal asymmetries by altering the relative magnitudes of net ebb and flood currents (Winterwerp et al. 2017; Guo et al. 2014).

1.4 Present research

The gross sediment dynamics described above have been assessed both numerically (Gatto et al. 2017; Guo et al. 2014) and using empirical data from estuaries such as the Hudson River (Geyer et al. 2001), Lower Passaic River (Chant et al. 2011), Huangmaohai Estuary (Gong et al. 2014), Delaware River (McSweeney et al. 2016), and Wairoa River (Pritchard and Green 2017). These studies have focused on the role of barotropic and estuarine circulation in promoting up-estuary residual sediment transport and on the role of river flow in promoting down-estuary residual sediment transport. As such, sediment dynamics in these estuaries have been classified into an importing and an exporting regime. In some cases, the exporting regime is described as primarily a flushing event where sediments trapped within the ETM zone are exported during the above-average river flows associated with the spring freshet (Pritchard and Green 2017; Geyer et al. 2001). Less studied is the impact of extreme or relatively infrequent river flow events on estuarine sediment dynamics and morphological evolution. For instance, in the Hudson River estuary, Geyer et al. (2001) hypothesize the occurrence of major erosion events occurring at roughly decadal time-scales. Similarly, a model application by Ralston et al. (2013) calculated significant erosion within the lower Hudson River estuary during an extreme event. Nonetheless, empirical evidence and understanding of sediment dynamics during extreme events is incomplete and is a subject of ongoing research (Ouillon 2018).

Here, we present an analysis of sediment dynamics in a short, narrow, microtidal estuary over short-term and long-term time-scales. The key objectives of our study are to assess estuarine suspended sediment dynamics over the full range of hydrologic conditions including extreme river flow events and to relate suspended sediment dynamics to the response of the bed in the particular estuary that is the focus of the analyses presented herein. The findings from this study are subsequently examined in the context of the long-term morphological evolution of estuaries. Furthermore, reliance on any single empirical line of evidence in developing a conceptual picture of sediment dynamics in such systems may lead to a biased result due to factors such as limited availability of measurements, episodic variations in the behavior of the system, etc. Therefore, a secondary objective of our study was to use a multiple lines of evidence approach including five separate and independent metrics in order to account for the limitations and uncertainties inherent in any single line of evidence. The analysis involves (1) measurements of residual sediment transport from a moored deployment over a range of river flows, (2) along-channel water column measurements over a range of river flows, (3) measurements of morphological change over the full range of river flows, (4) measurements of sediment erodibility, and (5) the results of a numerical hydrodynamic model. The mooring data are first used to assess sediment dynamics and relevant forcings during a limited range of river flows. The conclusions inferred from the analysis of the mooring data are compared against the along-channel water column measurements and measurements of morphological change for an assessment of sediment dynamics during a larger range of river flows including extreme events. The sediment dynamics are interpreted with the aid of a numerical hydrodynamic model and related to measurements of sediment erodibility for an assessment of bed dynamics driving sediment transport. Finally, the results are synthesized into a conceptual picture of sediment dynamics in the estuary. The following sections provide an overview of the study area, the data used, the analytical procedures involved, followed by a discussion of the results.

2 Site overview

The data presented in this paper come from the Lower Passaic River (LPR), a tidal estuary that is part of New York Harbor (Fig. 1). The LPR stretches approximately 28 km long from its mouth in Newark Bay at approximately River Mile (RM) 0.5, to the head-of-tide at Dundee Dam (approximately RM 17.5). Newark Bay is connected to New York Harbor and Raritan Bay (and the Atlantic Ocean) via the tidal inlets Kill van Kull and Arthur Kill, respectively. The width of the LPR ranges from approximately 600 m at its mouth, declining to about 200 m at RM 1.4, 150 m at RM 4.2, 120 m at RM 6.7, 90 m at RM 10.2, and 75 m at RM 13.5, i.e., locations about 2.25 km, 6.75 km, 10.75 km, 16.5 km, and 21.75 km, respectively, from the mouth of the river and relevant to the analysis presented here. Typical water depths along the thalweg in the LPR currently range between 5 and 7 m (with respect to mean sea level) and up to 17 m in the navigation channels in Newark Bay and Kill van Kull (Chant et al. 2011; Sommerfield and Chant 2010). The LPR is characterized by the remnants of a former navigation channel that was last dredged several decades ago and is no longer actively maintained at design depths (which ranged from about 10 m near the mouth of the river to about 4 m in the upper reaches; Chant et al. 2011). The LPR has relatively few sub-tidal shallows or tidal wetlands, features that have a bearing on the hydrodynamics of the system. The sediment substrate in the LPR is composed of predominantly fine sediments (particle diameters less than 63 μm; Moffatt & Nichol and Deltares 2019). Therefore, the sediment dynamics presented here relate primarily to fine sediments. The LPR is the subject of ongoing environmental cleanup and restoration studies; the data presented here were collected as part of this process.

Fig. 1
figure 1

Location map of the Lower Passaic River along with the locations of the in situ moorings

The hydrodynamic forcings within the LPR include the tides, estuarine circulation, and river flow. Semi-diurnal tides (period of 12.42 h, corresponding to the dominant semi-diurnal M2 constituent) entering Newark Bay through the Kill van Kull and Arthur Kill propagate to the LPR and the head-of-tide at Dundee Dam, forming an almost standing wave, with maximum currents typically occurring around mid-tide (Mathew and Winterwerp 2017). The tidal range varies from 0.9 to 2.1 m from neap to spring; the corresponding flow rates due to barotropic circulation (estimated from current measurements at RM 1.4) range approximately 150 to 300 m3/s (averaged over the half tidal cycle). In comparison, the annual average river flow over Dundee Dam is only about 34 m3/s (a few minor tributaries contribute approximately an additional 15% freshwater). Based on an extensive dataset of measurements, Chant et al. (2011) found salinity intrusion within the LPR during periods of low to average river flow (resulting in a partially mixed water column), with the saline water flushed out of the LPR during periods of high river flows. The extent of salinity intrusion, as indicated by the location of the salt front, is a function of the tidal phase, river flow, spring-neap cycle, as well as offshore mean water level fluctuations. Chant et al. (2011) also show an ETM co-located with the salt front, with its location and average SSC a function of the river flow.

3 Materials and methods

Several lines of evidence are presented in this paper. These include (1) fixed mooring measurements of SSC, currents, salinity, and water depth; (2) along-channel shipboard measurements of SSC, salinity, and water depth; (3) multi-beam bathymetry data; (4) measurements of sediment erodibility; and (5) a numerical hydrodynamic model. These data sources are described next followed by a discussion of the analytical methods used in this paper.

3.1 Water column: mooring data

The mooring data presented in this paper was collected during a period of several months (October 10, 2009 to December 16, 2009, and March 22, 2010 to July 24, 2010) at several locations within the LPR (shown in Fig. 1). The mooring locations span a range of salinity regimes during this period, ranging from freshwater tidal at RM 13.5, mostly freshwater tidal at RM 10.2, alternating freshwater and brackish at RM 6.7, and mostly brackish at RM 1.4 (with the exception of 2 days in the 2010 deployment when the salt front was located further seaward). The deployment also spanned a range of river flows, from below-average flows of about 5 m3/s to above-average flows up to 280 m3/s.

The deployment included moored (1) Acoustic Doppler Current Profilers (ADCP), (2) Conductivity-Temperature-Depth (CTD) sensors, and (3) Optical Backscatter (OBS) sensors performing in situ measurements every 12 min. The ADCPs were deployed in the bottom-mounted, upward-facing configuration and measured the depth-profile of flow velocity and echo intensity. The CTD and OBS sensors were deployed floating 0.9 m below the water surface and fixed 0.9 m above the bed for measurements of surface and bottom salinity, temperature, water depth, and turbidity. Water samples were regularly collected at the mooring locations and measured for SSC and related to turbidity measured by OBS, and to acoustic back-scatter (ABS; calculated from echo intensity following the methods of Deines 1999, and Wall et al. 2006). The resulting turbidity-SSC and ABS-SSC relationships were applied to the continuous time-series measurements of turbidity and ABS to estimate time-series of SSC at the mooring locations. The analysis presented in this paper relies on the ABS-estimated SSC time-series since it provides data on time-variable vertical profiles of SSC and primarily use data from the 2009 deployment. For reasons not well understood, ABS-SSC relationships for the 2010 deployment required separate regressions (yielding somewhat poorer correlation) for the relatively high river flow periods in the first half of the deployment and for the relatively low river flow periods in the second half of the deployment. Therefore, the data from the 2010 deployment are used in a limited manner as described subsequently.

Since the ADCP sensors were mounted on a tripod placed on the sediment bed, a fraction of the water column near the bed was not measured. Similarly, a fraction near the surface of the water column was not measured due to interference and binning artifacts. Both velocity and ABS-estimated SSC in these unmeasured zones were estimated by different methods as described next. Velocity in the unmeasured near-surface zone was estimated by assuming that fluid shear stress decreases linearly from measured values to zero at the surface. Fluid shear stress is calculated as

$$ \tau ={\mu}_t\frac{du}{dz} $$
(1)

where τ = fluid shear stress, μt = turbulent eddy viscosity, u = turbulent flow velocity, and z = vertical coordinate (z = 0 at the sediment-water interface). Velocity in the unmeasured near-bottom zone was estimated assuming a logarithmic profile:

$$ {\overline{u}}_{z,t}=\frac{u_{\ast }}{\kappa}\ln \left(\frac{z}{z_0}\right) $$
(2)

where \( \overline{u} \) = the turbulence-mean velocity, u = the bottom friction velocity, κ = 0.4 = the von Karman constant, and z0 = bottom roughness length = 0.4 mm, taken from a hydrodynamic modeling study of the LPR (HydroQual 2008). SSC in the unmeasured near-bottom and near-surface zones was extrapolated assuming that the vertical SSC profile follows the Rouse distribution (Van Rijn 1984):

$$ \frac{c}{c_a}={\left[\frac{a\left(1-\frac{z}{h_t}\right)}{z\left(1-\frac{a}{h_t}\right)}\right]}^{\beta } $$
(3)

where c = SSC measured at level z, ca = SSC at reference height a, ht = instantaneous water depth, and β = the Rouse number. β was estimated by least-squares fitting of the measured instantaneous SSC profiles.

Due to tidal variations, the measured profiles include a variable number of constant thickness ADCP bins with data over time. In order to assist with subsequent data analysis, the velocity and SSC profiles were converted to a sigma (σ) coordinate system which results in profiles with a constant number of layers but of variable thickness over time. The σ coordinate system is defined as

$$ \sigma =\frac{z-\eta }{H+\eta } $$
(4)

where η  = the instantaneous water level with respect to the reference height H. The instantaneous profiles were interpolated to a 20-layer σ grid.

3.2 Water column: shipboard data

The shipboard data span a wide range of river flows and include measurements of salinity and SSC (using CTD and OBS casts) over depth in the water column at several locations along the LPR and extending into Newark Bay. Some of this data is presented in Chant et al. (2011); the data presented here were collected by the same authors. Data during a below-average river flow condition of 8 m3/s on June 23, 2005 and an extreme event on March 16, 2010 with river flow of about 450 m3/s (corresponding to a return period of 25 years) are presented subsequently.

3.3 Bathymetry data

The bathymetry data consists of a series of multi-beam surveys performed in September 2007, November 2008, June 2010, October 2011, and September 2012. The freshwater inflow from Dundee Dam during this 5-year period is shown in Fig. 2 in relation to the annual average river flow rate. The river flow ranged from a low of about 1 m3/s in October 2007, to highs of 450 m3/s in March 2010 and March 2011 (return period of 25 years), and 700 m3/s in August 2011 (return period of 90 years). The surveys extended from the mouth of the river to RM 14.5 (about 23.3 km from the mouth), with data from individual surveys mapped to a 1.5 m by 1.5 m resolution grid. The various surveys were referenced to the same horizontal datum (North American Datum 1983, New Jersey State Plane) and vertical datum (National Geodetic Vertical Datum of 1929). River flow in the intervening periods between surveys varied, with certain periods (2007 to 2008 and 2011 to 2012) encompassing predominantly low-flow periods (i.e., no events greater than about 200 m3/s—return period of 2 years), and other periods (2008 to 2010 and 2010 to 2011) encompassing events greater than 200 m3/s (the relevance of the 200 m3/s threshold is discussed subsequently). Morphological change during these periods was calculated by performing bathymetric differencing of consecutive surveys.

Fig. 2
figure 2

Time-series of measured river flow at Dundee Dam during the periods encompassing the water column data and the various bathymetry surveys in the LPR in relation to the annual average river flow (red dashed line), and two morphologically relevant river flow thresholds discussed in this paper (blue dashed lines)

3.4 Sediment erodibility data

The sediment erodibility data used to support the analyses presented here is based on a series of erodibility measurements performed on surficial sediment cores from several locations in the LPR (Chesapeake Biogeochemical Associates, CBA 2006; measurements also presented in Mathew and Winterwerp 2017). Briefly, the measurements consisted of shallow cores collected from the LPR and subject to erosion experiments using a Gust Microcosm device. The resulting data were analyzed to calculate a depth-profile of τcr.

3.5 Numerical hydrodynamic model

Though not explicitly necessary for the present analyses, we use a three-dimensional hydrodynamic model, as this was available and well-calibrated. This model was developed by HydroQual (2008) using the Estuarine, Coastal and Ocean Model (ECOM) framework and applied as part of the environmental restoration activities in the LPR (U.S. Environmental Protection Agency, US EPA 2016). The model was applied over a domain that includes the LPR, Hackensack River, Newark Bay, and extending to the ends of the Arthur Kill and Kill van Kull—roughly the spatial extent shown in Fig. 1 (excluding the Hudson River and New York Harbor). The model grid resolution in the LPR ranges from 5 cells across the river at the mouth, decreasing to 4 cells near RM 1.4, to 3 cells near RM 4.2, and 2 cells above RM 15.7 (about 25 km from the mouth of the river). The average grid resolution in the LPR is about 40 m wide and 180 m long with 10 vertical layers (in a sigma coordinate system). Boundary conditions for the model include the measured river flow entering from the head-of-tide in the LPR and tributaries, and the water level, salinity, and temperature at the Kill van Kull and Arthur Kill boundaries (specified using the results of a regional-scale hydrodynamic model also described in HydroQual 2008). The model also includes meteorological forcings (winds, air temperature, relative humidity, barometric pressure, shortwave solar radiation, and cloud cover). It was calibrated against measured water levels, currents, temperature, and salinity at several locations within the model domain as well as validated against measurements from the 2009 and 2010 moored deployments described previously. The calibrated model was used for an assessment of currents and bed shear stresses under various steady-state river flows.

Following standard assumptions for hydrodynamic interactions at the bottom boundary, the effective bottom roughness used in Eq. (2) was assumed to be composed of form-related and grain-related fractions (Van Rijn 1993). The grain-related roughness, calculated as a function of the surficial sediment texture, is considered to generate the skin friction relevant for erosion. Therefore, skin friction was calculated as

$$ {\tau}_{SF}=\rho {\left[\frac{{\overline{u}}_t\kappa }{\ln \left(\frac{h_t}{2{z}_{0G}}\right)}\right]}^2 $$
(5)

where τSF = skin friction, ρ = density of water, the overbar represents depth-averaging, and z0G = grain roughness height calculated as

$$ {z}_{0G}=\frac{k_S}{30}=\frac{3{D}_{90}}{30} $$
(6)

where kS = Nikuradse grain roughness (Van Rijn 1993), D90 = particle diameter representing the 90% cumulative percentile of the sediment grain size distribution. The D90 was calculated using surficial sediment grain size distribution measurements in the LPR.

3.6 Decomposition methods

The majority of the analyses presented here relates to sediment dynamics using the mooring data. The analysis involves the decomposition of suspended sediment flux (SSF) into components attributable to the primary hydrodynamic forcings. This was accomplished by first decomposing the measured flow rates into barotropic, estuarine, and residual components followed by calculation of SSF associated with these processes. The term residual flow in the context of the analyses presented in this paper refers to the depth- and tidally-integrated quantity which in this case is mainly the river (freshwater) flow and flow induced by episodic subtidal barotropic events. Due to the lack of cross-sectional coverage in the mooring data, unless otherwise noted, the flow rates and SSF discussed in the remainder of the text refer to their channel width-normalized equivalents.

3.6.1 Flow decomposition

The flow decomposition uses a combination of analytical formulations and harmonic analysis of the σ-transformed measured currents and flow rates. Given the co-variance of tidal water depths and currents, in order to perform a mass conservative decomposition, the decomposition is applied to flow rates instead of currents. The measured instantaneous depth-dependent flow rate qz, t is

$$ {q}_{z,t}={u}_{z,t}\Delta {z}_t $$
(7)

where uz, t is the measured instantaneous velocity for sigma layer z in the water column with directionality assigned positive during flood and negative during ebb, and Δzt is the instantaneous thickness of the corresponding sigma layer. This instantaneous flow rate represents the combination of several components—a high-frequency component associated with barotropic circulation (qz, T), a low-frequency component associated with estuarine circulation (qz, E), and a low-frequency residual component that in this case is associated with the river (freshwater) flow (qz, R). High-and low-frequency are relative to the tidal period T (12.42 h). Accordingly, the instantaneous depth-dependent flow rate is

$$ {q}_{z,t}={q}_{z,T}+{q}_{z,E}+{q}_{z,R} $$
(8)

Various approaches were tested to perform the flow decomposition described by Eq. 8. The approach used in the analysis presented herein was chosen primarily for its ability to estimate the estuarine circulation component at locations that are relatively dynamic with respect to the salt front (freshwater and brackish). Appendix 1 presents the various approaches, along with a comparative evaluation of the results from these approaches.

The flow components in Eq. 8 were calculated by first applying an analytical formulation to extract the estuarine circulation component, followed by harmonic analysis to separate the barotropic and residual components. Estuarine circulation is classically defined as the tidally averaged deviation of the velocity profile from the depth-averaged velocity (Dyer 1997):

$$ {u}_{z,E}=\left\langle {u}_{z,t}-{\overline{u}}_t\right\rangle $$
(9)

where the overbar represents depth-averaging, angled brackets represents tidal-averaging, and uz, E is the velocity component associated with estuarine circulation. However, this definition includes vertical shearing by the logarithmic profile in Eq. 2. Application of Eq. 9 leads to results such as estuarine circulation being calculated landward of the salt front (see Appendix 1) due to deviations from the depth-averaged velocity (e.g., the logarithmic velocity profile). Such artifacts are avoided by modifying Eq. 9 to incorporate the logarithmic velocity profile using analytical formulations for velocity profiles that include the effects of bottom roughness and the pressure gradient induced by the longitudinal density gradient. Accordingly, Eq. 9 is rewritten as

$$ {u}_{z, Tot}=\left\langle {u}_{z,t}-{\overline{u}}_t\right\rangle $$
(10)

where uz, Tot is the total measured residual vertical circulation which represents the effects of estuarine circulation uz, E and the logarithmic velocity profile uz, log. To obtain a best estimate of uz, E from the data, uz, Tot is corrected for uz, log. In the presence of baroclinic effects, the velocity profile differs from the logarithmic profile, thus affecting the effective u. The logarithmic contribution is therefore obtained from the first-order analytical velocity profile for shear flow under the influence of a longitudinal salinity gradient (Winterwerp et al. 2006). This analytical velocity is indicated with the symbol v, to distinguish from the measured value u:

$$ {v}_{z,t}=-\frac{v_{\ast }}{\kappa}\mathit{\ln}\left(\frac{z}{z_0}\right)+\frac{1}{2}\frac{\alpha g{h}_t}{\kappa {v}_{\ast }}\left(z-{z}_0\right)\frac{1}{\rho}\frac{\partial S}{\partial x}={v}_{z,t,\mathit{\log}}+{v}_{z,t, bcl} $$
(11)

where α  = 0.8, g is the gravitational constant, ∂S/∂x is the measured longitudinal salinity gradient, and vz, t represents the analytical depth-dependent velocity. The first term on the right-hand side of Eq. 11 is the barotropic logarithmic velocity profile vz, t, log, and the second term vz, t, bcl accounts for the contribution from the longitudinal density gradient. The analytical shear velocity, v, is calculated using Eq. 12 which represents the depth-integration of Eq. 11:

$$ {\overline{v}}_t=-\frac{v_{\ast }}{\kappa}\left[\mathit{\ln}\left(\frac{h_t}{z_0}\right)-1\right]+\frac{1}{4}\frac{\alpha g{h_t}^2}{\kappa {v}_{\ast }}\frac{1}{\rho}\frac{\partial S}{\partial x}={\overline{v}}_{t,\mathit{\log}}+{\overline{v}}_{t, bcl} $$
(12)

Substituting the measured \( {\overline{u}}_t \) for the analytical \( {\overline{v}}_t \) in Eq. 12 resolves the analytical shear velocity v which includes contributions from barotropic and baroclinic components. Subsequently, the barotropic logarithmic velocity profile vz, t, log is assessed using Eq. 11 and used to calculate the estuarine circulation component uz, E by modifying Eq. 10 to include a correction for the logarithmic velocity profile (estimated as \( {v}_{z,t,\mathit{\log}}-{\overline{v}}_{t,\mathit{\log}} \)):

$$ {u}_{z,E}=\left\langle {u}_{z,t}-{\overline{u}}_t-\left({v}_{z,t,\mathit{\log}}-{\overline{v}}_{t,\mathit{\log}}\right)\right\rangle =\left\langle {u}_{z,t}-{\overline{u}}_t-{v}_{z,t,\mathit{\log}}+{\overline{v}}_{t,\mathit{\log}}\right\rangle $$
(13)

Given the tidal variations in water depth, the flow rate associated with estuarine circulation is calculated by incorporating the instantaneous sigma layer thickness in Eq. 13:

$$ {q}_{z,E}=\left\langle \left({u}_{z,t}-{\overline{u}}_t-{v}_{z,t,\mathit{\log}}+{\overline{v}}_{t,\mathit{\log}}\right){\Delta z}_t\right\rangle $$
(14)

Given the inequality in the semi-diurnal tides, the tidal-period averaging was performed over two tidal cycles (2*12.42 h) using a centered moving-window scheme. Note that the estuarine component thus calculated does not explicitly meet the constraint:

$$ \int {q}_{z,E} dz=0 $$
(15)

However, review of the residual from depth-integration of the calculated qz, E as well as comparison of the calculated near-bottom qz, E to the measured near-bottom salinity (presented subsequently) serves as a check on the decomposition formulations.

The difference of qz, t and qz, E represents the advection term, qz, A, which is the sum of the barotropic and residual flow components, as in Eq. 16:

$$ {q}_{z,A}={q}_{z,t}-{q}_{z,E}={q}_{z,T}+{q}_{z,R} $$
(16)

The depth-dependent barotropic and residual flow components were separated using a 35-h low-pass filter applied to the constituent periods of the Fourier-transformed qz, A time-series. Constituents with periods less than 35 h were considered to represent the barotropic component qz, T and constituents with periods greater than 35 h were considered to represent the residual flow component qz, R. Therefore, the barotropic component represents primarily tidal transport, whereas the residual flow component can include river flow, low-frequency barotropic events such as storm surges, as well as the effect of lateral variations in flow due to presence of river bends, etc.

Such flow decomposition techniques have been applied by other authors (Uncles and Jordan 1979; Winterwerp 1982; Costa 1989; Dyer 1997; Jay et al. 1997; Siegle et al. 2009). Alternative approaches using signal processing techniques have also been developed by other authors; see Lerczak et al. (2006) and Chant et al. (2011). These approaches were also applied to the data presented here; Appendix 1 includes a comparison of results from these approaches to the flow decomposition formulations described by Eqs. 816.

3.6.2 Suspended sediment flux decomposition

The various flow components and SSC time-series data were used to calculate SSF. Although decomposition techniques have been applied by others to separate SSC into components equivalent to that resulting from flow decomposition (Geyer et al. 2001; Scully and Friedrichs 2007; Siegle et al. 2009; Sommerfield and Wong 2011; Chant et al. 2011; Becherer et al. 2016), it was not applied in the present analysis for two inter-related reasons (see Appendix 2 for additional description of both arguments). The first reason involves the fact that the decomposition results in fluctuating SSC components that are often negative. Although mathematically tractable, negative SSC components are physically meaningless. The second reason is that integrated over depth and the tidal cycle, the scale of interest for the analyses presented here, SSF calculated with and without SSC decomposition are identical. Therefore, the tide- and depth-integrated SSF, FX, associated with given flow component was calculated as

$$ {F}_X={\int}_0^T{\int}_0^{h_t}{q}_{z,X}{c}_{z,t} dzdt $$
(17)

where subscript X refers to the various flow components described previously (residual, estuarine, and barotropic circulation), and cz, t is the measured depth-dependent instantaneous SSC. The net SSF representing the integrated effect of the individual components was calculated as the sum of the SSF associated with the individual flow components.

The SSF associated with the various flow components is subject to some uncertainty originating from the fact that the SSF decomposition procedure considers the various flow components and SSC as independent variables. However, SSC has a boundary condition (at the river bed) that scales over tidal time-scales. In other words, erosion from the bed scales as a function of the net force (i.e., the bed skin friction generated by the sum of the individual flow components) which varies over the tidal period. As shown in the next section, because the current associated with river flow adds to and enhances the ebb tidal currents, SSC during a given ebb tide is greater during periods of high river flow than during periods of lower river flow. However, the SSF decomposition does not apportion the incremental SSC generated by the higher river flow entirely to the SSF associated with river flow; such attribution implies that SSF associated with barotropic circulation would be independent of river flow which is the theoretically expected result. Rather, the SSF decomposition formulation associates the barotropic ebb flow rate with the net SSC, potentially resulting in depth- and tidally average SSF associated with barotropic circulation directed down-estuary during periods of high river flow. This limitation is entirely related to the use of empirical SSC data and is expected to be most prominently apparent in the SSF associated with barotropic circulation. SSC (and SSF) can reliably be attributed to the individual flow components only by using a numerical model such as used by Gatto et al. (2017). Nonetheless, the empirical data and SSF decomposition presented here help inform sediment transport dynamics and the relevance of various transport processes; the exact magnitude of impact of given forcing on SSF may be somewhat different.

4 Results

The results of the decomposition formulations include time-series of flow rates and SSF which are examined for dependencies with the measured river flow and salinity over a limited range of hydrologic conditions. The findings of sediment dynamics inferred from the review of SSF are also compared against morphological trends measured in the bathymetry data. The bathymetry data are also used to assess sediment dynamics over the full range of hydrologic conditions. The bed-water exchange processes of erosion and deposition inferred from the SSF and bathymetry data are interpreted with the aid of the hydrodynamic model and reviewed in the context of measurements of the vertical profile of τCr in the bed. These results are described in the following sub-sections, starting with an overview of SSC dynamics in relation to the primary hydrodynamic forcings.

4.1 SSC dynamics

Figure 3 shows a detailed view of SSC at RM 6.7 over a 2-day window to illustrate various sediment transport processes and the different dynamics during low and high river flow conditions. In general, during the low river flow period (left column of Fig. 3), SSC is in phase with current speed (and salinity), implying local resuspension of bed sediments. Also, both currents and SSC attain higher magnitudes during flood than ebb, indicating advection and net SSF directed up-estuary over the tidal cycle. During low river flow, due to its morphology, the LPR is characterized by flood-dominance in currents, and consequently, SSC tends to be higher during flood than ebb. Similarly, acceleration/deceleration asymmetries are apparent in the velocity data, with acceleration at the start of flood seen to be faster than deceleration at the end of flood, and vice versa during ebb. The effect of acceleration/deceleration asymmetries is also apparent in SSC, with SSC increasing relatively rapidly at the start of flood than at the start of ebb. In other words, SSC persists at a relatively low value of about 20 mg/L for a longer duration around high-water slack than around low-water slack. Therefore, both tidal asymmetry mechanisms described previously are apparent in the data at RM 6.7 during low river flow conditions implying that depth- and tidally integrated net SSF are directed up-estuary during this period. As mentioned previously, the other category of transport mechanisms, namely lag effects, are apparent only in a Lagrangian frame of reference and are therefore not explicitly apparent in Fig. 3.

Fig. 3
figure 3

Time-series of a, e the Dundee Dam discharge and water level at RM 6.7, b, f depth-average currents, c, g near-bottom salinity, and d, h depth-average SSC at RM 6.7. Hatched regions indicate duration of flood currents, i.e., directed up-estuary. Left and right columns show data during low and high-flow periods, respectively

The trends during high river flow contrast with low river flow in several respects (right column of Fig. 3). Although SSC is in phase with current speed during high river flow, in contrast to low river flow, both SSC and currents attain higher magnitudes during ebb than flood, thus indicating advection from up-estuary potentially in combination with erosion. The additional freshwater during high river flow (in conjunction with a setdown event apparent during the second half of December 11) results in higher currents during ebb than during flood. At the same time, the salt front is pushed down-estuary because of the higher river flow. Furthermore, during certain flood tides such as during the latter half of December 11th, SSC does not exhibit a peak as seen during other flood periods. This is due to the relatively low peak flood currents during this period, indicating that bed shear stresses are too low to cause erosion. In addition to the afore-mentioned trends in currents and SSC, during high river flow conditions, ebb duration is longer than flood duration. These trends imply that depth- and tidally integrated net SSF are directed down-estuary during this period.

Figure 4 includes along-channel transects in the LPR and in Newark Bay showing the measured bathymetry, salinity, and SSC during a below-average river flow condition of 8 m3/s on June 23, 2005, and an extreme event on March 16, 2010 with river flow of about 450 m3/s (return period of 25 years). Both transects show evidence of a well-developed ETM in the vicinity of the salt front (nominally defined as the location of the 2 PSU isohaline). The salt front and ETM respond to river flow—during the low-flow event, both the salt front and the ETM are located at RM 7, with depth-average SSC of about 75 mg/L within the ETM whereas during the high-flow event, the salt front and ETM are pushed to the mouth of the LPR with much higher depth-average SSC of about 250 mg/L within the ETM. The presence of the ETM co-located with the salt front within the LPR is indicative of the relevance of estuarine circulation to sediment dynamics in the LPR. The response to river flow is apparent in the SSC landward of the ETM which ranges about 20–40 mg/L during the low-flow event but increases to about 120 mg/L during the high-flow event.

Fig. 4
figure 4

Along-channel transects showing the cross-sectionally averaged bathymetry within the LPR (measured in 2007) and within the navigation channel in Newark Bay (both shown with an offset of − 5 m for plotting purposes). Also included are salinity contours and SSC (shaded colors) measured during a low river flows of about 8 m3/s and b during high river flows of about 450 m3/s and plotted relative to the measured water depth at the time of measurement. Both transects measured during mid-ebb

The sediment dynamics and the residual sediment transport apparent in Figs. 3 and 4 are further explored using the flow and SSF decomposition methods.

4.2 Hydrodynamic model and sediment erodibility

The hydrodynamic model was used to perform a series of simulations under constant salinity (at the marine boundaries in the Arthur Kill and Kill van Kull), over a spring-neap cycle, and for various river flow rates. Applied freshwater flow rates at the head-of-tide ranged from 0 to 500 m3/s (representing an event with a return period slightly greater than 25 years). These simulations help understand the response of the LPR to river flow and the impact on sediment dynamics. Figures 5 and 6 show various metrics from these simulations, calculated using the spring-neap mean of various cross-sectionally averaged quantities. Figure 5 shows results for computed currents and salinity, and Fig. 6 shows results for computed skin friction. The impact of localized variations such as a reduction in cross-sectional area due to a rock outcrop immediately up-estuary of RM 8, and an increase in cross-sectional area due to a widening of the river at RM 4.2 are apparent in the calculated peak tidal currents and skin friction. Asymmetry in peak tidal currents was quantified as the ratio of the peak flood-current to peak ebb-current, with values greater than one denoting flood-dominance and values less than one denoting ebb-dominance. At any given location, the system is flood-dominant at low river flows, transitioning to ebb-dominance with increasing river flow. The river flow associated with the transition from flood- to ebb-dominant increases with distance down-estuary and the entire river exhibits ebb-dominance at river flow slightly greater than 50 m3/s. Salinity responds in a similar manner, with the salt front pushed seaward as river flow increases and pushed out of the river at flow greater than 200 m3/s.

Fig. 5
figure 5

Results of a numerical model showing impact of river flow rate on a peak flood-phase currents, b peak ebb-phase currents, c asymmetry in peak currents, and d salinity

Fig. 6
figure 6

Results of a numerical model showing impact of river flow rate on a peak flood-phase skin friction, b peak ebb-phase skin friction, and c asymmetry in peak skin friction

Skin friction at the bed-water interface shown in Fig. 6 responds in a similar manner as currents. The upper two panels in Fig. 6 also show the τcr measured on sediment cores. As shown in Mathew and Winterwerp (2017), typical values of τcr in the LPR range from about 0.04 Pa at the surface of the cores, increasing to 0.4 Pa at a depth of about 2–4 mm below the surface. This thin layer (2–4 mm thick) of easily erodible sediments at the surface of the cores was shown to be indicative of a pool of sediments (referred to as the fluff layer) that is resuspended every tidal cycle (once during flood and again during ebb) and redeposited around slack water. Mathew and Winterwerp (2017) also present arguments supporting the importance of the fluff layer for the net transport of fine sediments against the direction of residual (river) flow. This is also seen in the skin friction results presented in Fig. 6—during low river flows (nominally defined as 0–10 m3/s in this context), peak skin friction ranges between 0.04 and 0.4 Pa over nearly the entire length of the LPR. In other words, erosion is expected to be restricted to the fluff layer during such conditions. Only when bed skin friction exceeds about 0.4 Pa is erosion expected to extend to deeper depths. The results for peak ebb skin friction show that as river flow increases beyond 25 m3/s, progressively larger areas of the LPR experience skin friction greater than 0.4 Pa. At river flows greater than about 100 m3/s, areas landward of RM 8 experience skin friction greater than 0.4 Pa, and at river flows beyond 200 m3/s, the entire river experiences skin friction greater than 0.4 Pa. These relative comparisons of τcr for various sediment strata and the system response to increasing river flow are used to inform the interpretation of sediment dynamics and morphological changes from the analysis of SSF and bathymetry data. Note that the values of τcr presented here are associated with predominantly cohesive sediments. In particular, the sediment substrate in areas above RM 14.5 is predominantly composed of sand and gravel, and therefore, the values of τcr presented here are not representative of the sediments in those areas.

4.3 Flow decomposition

This section presents a validation of the decomposition procedure, by comparing two of the three calculated flow components against measured metrics.

4.3.1 Residual flow

The residual flow rate is estimated to represent primarily river flow; this was assessed by comparison against the measured freshwater inflow from the head-of-tide. The calculated depth-integrated residual flow rate at given location was extrapolated over the cross-section using the effective channel width (calculated as the ratio of cross-sectional area to water depth at the mooring location; both at mean sea level). Figure 7 shows the results of this comparison for the 2009 moored deployment, including quantitative metrics. The calculated flow rates at all locations tend to reproduce the general temporal trends seen in the measured river flow, with episodic high-flow events during the last week of October, and during the first half of December. However, the magnitudes differ—at the two seaward locations (RM 1.4 and RM 4.2), calculated and measured flow rates are more comparable during high-flow periods than during low-flow periods. Occasionally, during low-flow periods, the flow decomposition procedures also result in calculated residual (river) flow rate qz, R directed up-estuary which is a spurious result. The alternative approaches applied in Appendix 1 also give such spurious results suggesting that this is not an artifact related to the flow decomposition formulations. Rather, these artifacts at RM 1.4 and RM 4.2 likely relate to uncertainties in measurement, sampling location, cross-sectional averaging and possibly secondary flows and preferential flow paths induced by sharp bends in the river. Such errors are inherent in the estimation of residual terms from gross fluxes in tidal and estuarine settings (Jay et al. 1997). In contrast, calculated flow rates at the three up-estuary locations are more comparable to measured flow rates over the full range of flows. The relatively consistent tendency for over-prediction at these three locations may be an indication of cross-channel variations in currents not captured by the mooring data. It should be noted that the decomposition of sediment fluxes is performed on a channel width-normalized basis and is therefore not affected by potential artifacts associated with cross-channel variations. The comparison of the calculated and measured river flow rates at the three up-estuary locations is taken as a validation of the flow decomposition formulations, thus supporting its use in SSF decomposition.

Fig. 7
figure 7

Time-series comparison of measured river flow rate at the head-of-tide and calculated residual (river) flow rate at the various mooring locations for the 2009 moored deployment. Performance metrics shown include the bias (calculated as the difference of mean calculated and mean measured flow rate, with the mean values calculated over the duration of the deployment), and the relative bias (calculated as the ratio of bias to the mean measured flow rate). Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4

4.3.2 Estuarine exchange flow

The calculated flow rates associated with estuarine circulation were also validated, albeit in a qualitative manner, using measured near-bottom salinity data. Since estuarine circulation is expected only seaward of the salt front, the presence or absence of salinity at a given location is indicative of the presence or absence of estuarine circulation. Therefore, the presence or absence of calculated estuarine circulation is compared against the presence or absence of salinity at given location as a qualitative check on the flow decomposition formulations. Figure 8 shows the results of this comparison for the 2009 moored deployment. At RM 6.7, which is the most dynamic location with respect to the salt front, the calculated near-bottom flow component associated with estuarine circulation shows consistent temporal trends as near-bottom salinity. Estuarine circulation is seen to occur only when the salt front is located landward of RM 6.7, a result consistent with theoretical expectations. The other locations are less dynamic with respect to estuarine circulation, with estuarine circulation persisting at RM 1.4 and RM 4.2 over the entire deployment. Similarly, the salt front and estuarine circulation is seen to extend up to RM 10.2 for a few days preceding October 17. Although estuarine circulation is calculated during brief periods at RM 10.2 (October 24 and November 14), and RM 13.5 (October 25 through November 2), during periods when the salt front is located seaward of these stations, the magnitude of the calculated flow rate is small. These false signals are likely related to deviations from the theoretical logarithmic vertical profile for currents. The significantly higher flow rates associated with estuarine circulation at the other locations, and its co-dependence with measured salinity is taken as a validation of the flow decomposition methods, thus supporting its use for SSF decomposition. It should be noted that neglecting the correction for the logarithmic velocity profile results in artifacts such as estuarine circulation being calculated even at RM 13.5 which is located landward of the salt front during the entire deployment. This is further elaborated in Appendix 1.

Fig. 8
figure 8

Time-series comparison of measured near-bottom salinity and calculated near-bottom flow rate associated with estuarine circulation at the various mooring locations for the 2009 moored deployment. Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4

4.4 Suspended sediment fluxes

The results of the SSF decomposition include time-series of the net SSF and SSF associated with various flow components. These are evaluated for their dependency with the primary hydrodynamic forcings.

4.4.1 SSF associated with estuarine circulation

Figure 9 shows the time-series of SSF associated with estuarine circulation relative to the measured near-bottom salinity for the 2009 moored deployment. As with the estuarine circulation flow rates, SSF associated with estuarine circulation tends to be correlated with the salt front location and is directed up-estuary at locations seaward of the salt front. The only exceptions are during the neap tides of October 24–30 and November 22–28 at RMs 4.2 and 1.4. This is due to a combination of low tidal energy (i.e., lesser amplitude of tidal currents and therefore lesser sediment resuspension and lower near-bottom SSC), strong stratification (also seen in the relatively small intra-tidal fluctuations in near-bottom salinity) which dampens vertical mixing, and elevated river flow (on October 25 and 28). The runoff events likely resulted in additional SSC loadings to the river (from the head-of-tide as well as stormwater outfalls in the estuary), resulting in measurements of higher near-surface SSC than near-bottom SSC (also seen in the turbidity-estimated SSC; concentrations on the order of 10–20 mg/L). Although not definitively related to a runoff event, the period of November 22–28 also experiences strong salinity stratification, and higher near-surface SSC than near-bottom SSC. The net result is a negligible up-estuary or even down-estuary SSF despite estuarine circulation during these periods, especially at RM 1.4.

Fig. 9
figure 9

Time-series comparison of measured near-bottom salinity and calculated SSF associated with estuarine circulation for the 2009 moored deployment. Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4. Fluxes calculated by integrating over depth and over time (two tidal cycles) using a moving-window scheme. Positive and negative values indicate fluxes directed up-estuary and down-estuary, respectively

4.4.2 SSF associated with barotropic circulation

Figure 10 shows the calculated SSF associated with barotropic circulation as a function of the measured freshwater flow rate at Dundee Dam for the 2009 moored deployment. The comparison indicates SSF dynamics that are dependent on river flow, with low river flows generally associated with up-estuary SSF and high river flows associated with down-estuary SSF. The up-estuary SSF at low river flows can be attributed to the various barotropic processes described previously—lag effects and tidal asymmetries. Two of these processes can be seen in the results for the zero river flow case shown in Figure 5—the decrease in tidal currents with distance up-estuary, and flood dominant tidal currents. Theoretically, the barotropic lag effects and tidal asymmetries, and by extension their role in promoting net up-estuary SSF, are independent of river flow. However, as explained previously, the lack of reliable decomposition of SSC results in the apparent dependency of SSF with river flow and specifically, the down-estuary SSF at high river flows. Nonetheless, as apparent from the up-estuary SSF at low river flows, Fig. 10 shows that in the LPR as a whole (landward of RM 1.4), lag effects and tidal asymmetries influence sediment dynamics and induce up-estuary flux of sediments. Although we cannot distinguish between lag effects and tidal asymmetry, the former is probably more important in transporting sediment up-estuary because up-estuary transport is measured only at low river flows (below erosion threshold of parent bed), and Fig. 3 suggests that only a limited amount of sediment is mobilized in a short period during the flood tide.

Fig. 10
figure 10

Calculated SSF associated with barotropic processes as a function of the measured river flow rate at the head-of-tide for the 2009 moored deployment. Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4. Fluxes calculated by integrating over depth and over time (two tidal cycles) using a fixed-window scheme. Positive and negative values indicate fluxes directed up-estuary and down-estuary, respectively

4.4.3 SSF associated with residual flow

Figure 11 shows the calculated SSF associated with the residual flow as a function of the measured freshwater flow rate at the head-of-tide at Dundee Dam for the 2009 moored deployment. As expected, net SSF associated with the residual flow follows the direction of the residual current, i.e., river flow, and is directed down-estuary at all locations. The only exceptions are net up-estuary SSF at RMs 4.2 and 1.4 during a handful of tidal cycles during low river flows. These are associated with the artifacts associated with the estimated river flow rate at these locations discussed previously in association with Fig. 7.

Fig. 11
figure 11

Calculated SSF associated with the residual flow as a function of the measured river flow rate at the head-of-tide for the 2009 moored deployment. Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4. Fluxes calculated by integrating over depth and over time (two tidal cycles) using a fixed-window scheme. Positive and negative values indicate fluxes directed up-estuary and down-estuary, respectively

The relationship between SSF and river flow shows certain patterns that inform erosion and deposition patterns in the LPR. In general, as river flow increases, net down-estuary SSF increases due to a combination of additional sediment load from the head-of-tide as well as erosion within the LPR. Comparing the SSF at RM 13.5 to locations down-estuary shows the impact of erosion within the LPR. There is a general trend of increasing SSF from RM 13.5 to RM 6.7, indicating erosion within this reach, and a general trend of decreasing SSF from RM 6.7 to RM 1.4, indicating deposition within this reach. These patterns of erosion and deposition show a dynamic system, with spatially variable patterns of erosion and deposition that are dependent on river flow and are further elaborated upon in the following section. The SSF trends with river flow were also assessed using SSF calculated using the measured river flow instead of the calculated residual flow as an assessment of the uncertainty in the calculated residual flow rates presented in Fig. 7. The resulting trends were generally similar to Fig. 11 suggesting that the previously noted errors in the calculated residual flow do not impact the overall conclusions presented here.

4.4.4 Net SSF

The net SSF represents the integrated result of all the transport processes described previously (barotropic, estuarine, and fluvial). As seen in Fig. 12 for the 2009 moored deployment, the magnitude and direction of net SSF are strongly dependent on the river flow rate. With the exception of RM 13.5, areas landward of which are characterized by predominantly sandy sediments, the general trend is of net up-estuary SSF at low river flow, i.e., importing conditions, and a transition to net down-estuary SSF at higher river flows, i.e., exporting conditions. During low river flow, the up-estuary transport is driven by estuarine circulation, and barotropic lag effects and tidal asymmetries. The latter two transport mechanisms also persist landward of the salt front. Increasing river flow reduces and/or eliminates the processes responsible for up-estuary SSF (for instance, increasing river flow reverses the flood dominance in currents), increases ebb currents, and consequently, net SSF transport is directed down-estuary at high river flows. The river flow rate associated with the transition from net up-estuary to down-estuary SSF increases towards the mouth of the LPR.

Fig. 12
figure 12

Calculated net SSF (integrating all transport processes) as a function of measured river flow rate at the head-of-tide for the 2009 moored deployment. Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4. At locations seaward of RM 13.5, during above-average river flows, SSF is color-coded relative to the magnitude of SSF at the location immediately landward, with red and blue indicating erosion and deposition, respectively, in the intervening reach. Fluxes calculated by integrating over depth and over time (two tidal cycles) using a fixed-window scheme. Positive and negative values indicate fluxes directed up-estuary and down-estuary, respectively

The SSF in Fig. 12 during above-average river flows is also color-coded to indicate whether SSF at given location is greater or less than the SSF at the location immediately landward, thus indicating erosion or deposition in the intervening reach. The data show a general trend of erosion at the more landward locations (e.g., between RM 13.5 and RM 10.2) gradually transitioning to deposition at the more seaward locations (e.g., between RM 4.2 and RM 1.4) during above-average river flows. For context, using dry density of 0.5 MT/m3 (based on measurements in the LPR; see Moffatt & Nichol and Deltares 2019), 1 MT/m of erosion in the reach from RM 13.5 to RM 10.2 is equivalent to a depth of erosion of approximately 0.4 mm which is less than the estimated thickness of the fluff layer. The erosional response at the up-estuary locations can be interpreted using Fig. 6—barring a few localized exceptions seaward of RM 13.5 and the area above RM 15, skin friction is less than the threshold τcr of 0.4 Pa in the majority of the LPR for river flows up to 100 m3/s. Therefore, the more consolidated bed underneath the fluff layer is likely not eroded during such river flows, thus limiting scour depths and net down-estuary SSF. In other words, erosion is limited to the fluff layer in the majority of the LPR during river flows up to 100 m3/s, with some of the eroded material deposited at locations farther seaward. Furthermore, as seen in Fig. 6, peak ebb skin friction exceeds 0.4 Pa at the more seaward locations (e.g., in the vicinity of RM 1.4) only when river flow exceeds about 200 m3/s. Therefore, erosion is expected to be limited to the fluff layer for river flows up to about 200 m3/s at the more seaward reaches.

Figure 13 shows the tidally integrated net SSF for the 2010 moored deployment as a function of river flow. The general trends noted in Fig. 12 are apparent in this dataset as well, with the magnitude and direction of net SSF dependent on river flow, low river flow periods associated with net up-estuary SSF, and a transition to net down-estuary SSF at higher river flows. Note the higher range of river flow rates for this dataset, up to 280 m3/s, as compared to the 2009 moored deployment. Although the 2010 moored deployment includes river flows greater than 200 m3/s (threshold based on river-wide impacts as seen from the peak ebb skin friction in Fig. 6), sediment dynamics in this dataset are subject to certain qualifiers that are also relevant to one of the findings of this paper. In contrast to the erosional response noted in the more landward reaches during above-average river flows in the 2009 moored deployment, comparison of the net SSF during flow rates greater than 200 m3/s at the various locations shows relatively similar values between RM 13.5 and RM 4.2, with only RM 1.4 showing somewhat higher SSF. The time-history of river flow is relevant in interpreting sediment dynamics during river flows greater than 200 m3/s in this dataset. The 2010 deployment started on March 22, 2010 following a high river flow event of 450 m3/s (return period of 25 years) on March 16, 2010. Subsequently, river flow decreased to 170 m3/s on March 22, decreased further to 120 m3/s on March 28, before increasing to 280 m3/s on April 1. As elaborated in the following section, the 450-m3/s event on March 16 caused significant erosion (nominally defined as depths in excess of the fluff layer thickness) within the LPR which influenced the sediment dynamics of the system during the subsequent 280-m3/s event. Measurements of sediment erodibility in the LPR (Borrowman et al. 2006) show decreasing erodibility (τcr increasing) with depth in the bed (measurements up to about 40 cm deep). This suggests that the less erodible bed strata exposed as a consequence of erosion during the 450-m3/s event would not be expected to erode under the relatively lower skin friction of the following 280-m3/s event. The lack of significant erosion noted in Fig. 13 during flows greater than 200 m3/s is consistent with this explanation derived from empirical measurements of sediment erodibility. Furthermore, this pattern of significant erosion during the 450-m3/s event on March 16 and relatively little erosion during the subsequent 280-m3/s event was also reproduced by numerical sediment transport models (U.S. Environmental Protection Agency (US EPA) 2016; Moffatt & Nichol and Deltares 2019). Therefore, the lack of significant erosion in the river during flows greater than 200 m3/s in the 2010 moored deployment dataset is related to time-history of river flow rather than being indicative of sediment dynamics had these flow conditions occurred sufficiently removed in time from the preceding 450-m3/s event.

Fig. 13
figure 13

Calculated net SSF (integrating all transport processes) as a function of measured river flow rate at the head-of-tide for the 2010 moored deployment. Comparisons at a RM 13.5, b RM 10.2, c RM 6.7, d RM 4.2, and e RM 1.4. Fluxes calculated by integrating over depth and over time (two tidal cycles) using a fixed-window scheme. Positive and negative values indicate fluxes directed up-estuary and down-estuary, respectively

Although the SSF from the mooring data do not allow for evaluation of events exhibiting signification erosion, the periodic bathymetry surveys encompassed periods with river flows up to 700 m3/s as well as below-average river flows, allowing for an evaluation of the river dynamics under importing, exporting, and erosional conditions. The along-channel transect during the 450-m3/s event on March 16, 2010 presented in Figure 4 also provides evidence of significant erosion and is discussed in the next section.

4.5 Sediment dynamics from morphological data

The series of consecutive multi-beam bathymetry surveys during 2007 through 2012 were analyzed to calculate morphological changes during the intervening periods. The calculated bathymetric changes between successive individual surveys were averaged cross-sectionally and longitudinally over 1-mile (1.6 km) intervals for an assessment of large-scale morphological changes. Figure 14 shows the results of this analysis, plotted as a profile of morphological changes along the longitudinal axis of the LPR.

Fig. 14
figure 14

Longitudinal profile of measured bathymetric changes over various survey periods laterally and longitudinally averaged over 1.6 km (1 mile) intervals. Positive values indicate net deposition and negative values indicate net erosion. Morphological change over the a 2007–2008, b 2008–2010, c 2010–2011, and d 2011–2012 survey periods (see also Fig. 2 for the hydrograph during these periods). Positive and negative values indicate deposition and erosion, respectively

The morphological changes show a correspondence with river flow that is similar to the SSF trends in the mooring data. For instance, river flow in the periods encompassed by the 2007–2008 and 2011–2012 surveys was relatively low compared to the periods encompassed by the 2008–2010 and 2010–2011 surveys. The former periods had only three events and one event, respectively, above 100 m3/s, and no events above 200 m3/s (see Fig. 2). However, the latter periods had several events above 100 m3/s as well as 200 m3/s, with 450-m3/s events in March 2010 and March 2011 and a 700-m3/s event in August 2011. The impact of river flow is seen in the morphological changes measured during these periods, with the 2007–2008 and 2011–2012 periods exhibiting a generally depositional signal, albeit with some localized erosion, whereas the high river flow periods (2008–2010 and 2010–2011) show erosion over much of the LPR. It should be noted that the erosion measured between RM 1–2 during the 2007–2008 period is mostly due to ship-induced scour; this reach of the LPR includes several active maritime terminals. With the exception of a few events with river flow between 100 and 200 m3/s during 2007–2008 and 2011–2012, the LPR mostly experienced river flows corresponding to importing regimes. Therefore, the net infilling during these periods is considered to be associated with the estuarine circulation and barotropic processes as seen in the analysis of SSF from the mooring data.

In contrast, the two high river flow periods (2008–2010 and 2010–2011) experienced events with river flow in excess of 200 m3/s when the consolidated bed underneath the fluff layer is expected to be subject to scour (as seen from Fig. 6). This helps explain the erosion noted in the 2008–2010 and 2010–2011 periods. There are certain localized signals such as (1) relatively lower erosion landward of RM 7 during 2010–2011 than during 2008–2010 which may be caused by armoring of the bed in these areas during the 2008–2010 period, (2) relatively large erosion between RM 11–12 during 2010–2011 which was caused by a temporary constriction of the river due to bridge construction, and (3) erosion extending farther seaward during 2010–2011 than during 2008–2010 due to the higher river flow during the former period. The overall pattern noted during these two high river flow periods is one of scour, with scour depths largely in excess of the fluff layer thickness, a result consistent with the estimates of skin friction from the numerical model and its comparison to the τcr for the fluff layer and the more consolidated bed underneath. Furthermore, the pattern of significant erosion during the 450-m3/s and 700-m3/s events was also reproduced by numerical sediment transport models (U.S. Environmental Protection Agency (US EPA) 2016; Moffatt & Nichol and Deltares 2019).

The erosion noted in the 2008–2010 period is also seen directly in the water column data. Specifically, the 450-m3/s event on March 16, 2010 was the largest event measured during the 2008–2010 period and the likely cause of the erosion measured in the bathymetry data. Measurements during this event are presented in Fig. 4b—both the salt front and ETM were pushed to mouth of the LPR, with elevated SSC within the ETM as well as landward of the ETM (as compared to the low river flow transect in Fig. 4a). Depth-average SSC landward of the ETM averaged approximately 120 mg/L whereas SSC at the head-of-tide (estimated using a rating curve) is only about 50 mg/L (Moffatt & Nichol and Deltares 2019) indicating significant erosion landward of the ETM. Two similar transects (data not shown) collected during other phases of the tide during this same event show similar trends. This is a direct line of evidence indicating scouring conditions at river flows greater than 200 m3/s and is consistent with the measured morphological changes. Both the 450-m3/s and 700-m3/s events in 2010 and 2011 represent extreme river flow conditions where the salt front is pushed to the mouth of the river and the entire river experiences skin friction greater than 0.4 Pa, with erosion of the more consolidated bed beneath the fluff layer.

Therefore, the erosion measured in the morphological data over the 2008–2010 and 2010–2011 periods is considered to have been caused by the relatively extreme 450-m3/s and 700-m3/s river flow events during these periods. Although not direct evidence linking the measured erosion to the extreme river flow events, the various lines of evidence discussed above support the notion of significant estuary-wide scour during river flows in excess of 200 m3/s. The findings from the analysis of the morphological changes are combined with the results of SSF dynamics from the mooring data and used to develop a conceptual picture of estuarine sediment dynamics in the following section.

5 Discussion

The results of the SSF decomposition allow for an assessment of the various hydrodynamic forcings on sediment transport in the LPR. Since the residual, i.e., river flow, is directed down-estuary, net sediment transport associated with river flow is also directed down-estuary. In contrast, barotropic and baroclinic processes promote up-estuary transport. During low river flow, lag effects and tidal asymmetries induced by barotropic circulation, and estuarine circulation induced by salinity intrusion and mixing typically lead to net up-estuary transport of sediment. Increasing river flow modifies both the barotropic and baroclinic transport pathways. The additional freshwater inflow pushes the salt front seaward, thus limiting the spatial extent of estuarine circulation as a transport mechanism. Similarly, increasing freshwater inflow reverses the flood-dominance in barotropic currents and also alters lag effects resulting in net SSF directed down-estuary during high river flows.

The primary finding of this paper is a conceptual picture of sediment dynamics, with the estuary following one of three modes—importing, flushing, or scouring. During low river flow conditions, tidal and estuarine processes dominate sediment transport and the system imports sediments from locations farther seaward (as well as from the head-of-tide) and is a net sink for sediments. In general terms, low river flow conditions represent an importing regime (regime I). Conversely, increasing river flow is associated with net down-estuary SSF at all locations, and a net export of sediment from the system. In other words, high river flow conditions represent an exporting regime, when the system is a net source for sediments and exports sediments to locations farther seaward. The exporting regime can be further divided into two conditions, distinguished by morphological impact—flushing conditions, thus referred to since only the easily erodible sediment in the fluff layer are expected to be eroded (regime II), and scouring conditions, thus referred to since the more consolidated bed layers underneath the fluff layer are eroded (regime III). In the LPR, the threshold between flushing and scouring regimes is estimated to be approximately 200 m3/s for the system as a whole, although locally, such as in the areas landward of RM 8, the threshold is lower. The threshold of 200 m3/s also represents conditions when the salt front is located near the mouth of the river.

Figure 15 shows a conceptual picture of the sediment dynamics in the LPR following the three regimes. The freshwater flow rates associated with the thresholds between regimes are spatially variable. For instance, as seen in Fig. 12, the threshold between regimes I and II is estimated to be about 10 m3/s at RM 10.2, and between 20 and 30 m3/s at RM 1.4, and as seen in Fig. 6b, the threshold between regimes II and III is estimated to be about 100 m3/s at RM 8, and about 200 m3/s at RM 1.4. The spatial variation in flow thresholds means that these regimes may not be synchronous along the length of the river, e.g., an event that produces regime II behavior at an up-estuary location may induce a regime I behavior at a down-estuary location. Comparison of these flow thresholds to the long-term (over a hundred years) hydrograph in the LPR shows the temporal prevalence of these regimes. Freshwater flow rate of 30 m3/s is exceeded about 40% of the time, 100 m3/s about 5% of the time, and 200 m3/s less than 1% of the time. In other words, importing regimes persist the majority of the time and exporting regimes, in particular, exporting-scouring regimes, occur only during a relatively small fraction of the time. However, as seen in Figs. 12 and 13, net down-estuary SSF at high river flows are orders-of-magnitude larger than the up-estuary SSF at low river flows, suggesting that even though exporting regimes occur a small fraction of the time, net export during such conditions is a significant fraction of the long-term sediment mass budget. As seen from the infill rate during 2011–2012 following the erosional events of 2008–2010 and 2010–2011, the recovery (infilling) time is estimated to be on the order of about 1 year or more (depending on the magnitude of the event). This implies that the impact of these extreme events persists for an extended period.

Fig. 15
figure 15

Conceptual picture of estuarine sediment transport regimes—importing, exporting-flushing, and exporting-scouring. The gradient in colors indicates the increasing SSC and increasing potential for morphological change from regime I to regime III. The lines separating the regimes are inclined to indicate the generally decreasing river flow thresholds between regimes with distance up-estuary. The lines separating the regimes are jagged to indicate that local variations are possible around this general pattern of decreasing river flow threshold between regimes with distance up-estuary; for instance, in response to local changes in cross-sectional area. The lines separating the regimes are dashed to indicate that temporal variations are possible for the river flow threshold separating the regimes; for instance, in response to local erosion and deposition, dredging, etc

With the exception of regime III, the general conceptual picture of importing and exporting regimes described here is consistent with studies of other estuaries around the world. Based on analysis of SSF in the Hudson River, USA, Geyer et al. (2001) concluded net import of sediments during average flow conditions, and net export during elevated river flow (coincident with spring tides). Similar observations have also been reported in other estuaries (Huangmaohai Estuary, China by Gong et al. 2014; Delaware River, USA by McSweeney et al. 2016; and Wairoa River, New Zealand by Pritchard and Green 2017) as well as in the LPR by Chant et al. (2011). River discharge (possibly in combination with spring tides) is the main factor in determining the transition from importing to exporting regimes. Wider estuaries with significant inter-tidal areas or shallow tidal flats may also show cross-channel variations in importing and exporting regimes. For example, in the Delaware River, USA, McSweeney et al. (2016) show net sediment export across the entire cross-section during high river flow periods, but net import along the deeper portions of the channel and net export along the shallow flanks during low river flow periods. A similar result was also reported by Scully and Friedrichs (2007) in the York River, USA. The main contribution of this paper towards the understanding of estuarine sediment dynamics is the assessment of sediment dynamics over the full range of hydrologic conditions and over longer timescales which results in the delineation of the exporting regime into exporting-flushing and exporting-scouring regimes.

The finding of regime III is relevant in the context of the long-term morphological evolution of estuaries. Meade (1969) noted that the typical morphological behavior of estuaries in the Holocene epoch is a slow infill with sediment originating from fluvial and marine sources, eventually reaching a (dynamic) equilibrium between sediment supply and export. Further, the exporting conditions necessary for maintaining the long-term equilibrium are thought to be extreme events when the flow is directed seaward at all depths within the estuary and the salt front is pushed out of the estuary (Meade 1969). This general narrative was also proposed by Geyer et al. (2001) who found the short-term sediment import from fluvial and marine sources in the Hudson River estuary to be significantly higher than the long-term net sedimentation rate measured from radionuclide data. Since the Hudson River is at morphodynamic equilibrium, Geyer et al. (2001) hypothesized the occurrence of episodic events at decadal time-scales as the mechanism causing erosion and export of accumulated sediments, thus maintaining the long-term morphodynamic equilibrium and accounting for the difference between the short-term and long-term sedimentation rates. The hypothesis of episodic erosion in the lower Hudson River estuary was also verified by Ralston et al. (2013) who calculated significant erosion in model simulations of a river flow event with a return period of 60 years. The finding of regime III events in the LPR is consistent with the hypothesis of Meade (1969) and Geyer et al. (2001). Therefore, regime III events may be considered representative of the mechanism responsible for maintaining long-term estuarine morphodynamic equilibrium in such estuaries.

The conceptual picture outlined in Fig. 15 is also responsive to morphological changes. For instance, deepening of the system (due to dredging or erosion) or a reduction in water depths due to sedimentation may change the threshold river flow rates separating the three regimes (in addition, tidal propagation may also be impacted by morphological change). In other words, the morphodynamic response of the system at any given time is a function of the existing river flow rate, as well as the time-history of river flow rate and the morphological state of the system. This is seen in the lack of significant erosion at above-average river flows in the 2010 moored deployment data which is likely related to the armoring effects of the preceding 450-m3/s event on March 16, 2010. Another example relates to the long-term (decadal-scale) morphological evolution of the LPR relative to the short-term event-driven dynamics conceptualized by Fig. 15. As shown in Chant et al. (2011), the long-term large-scale morphological behavior of the LPR is an infilling trend following the dredging of the navigation channel in the late 1940s which deepened the LPR below its morphological equilibrium. However, despite this long-term infilling trend, short-term trends are more variable—high flows induce erosion as seen in the morphological change during 2008–2010 and 2010–2011, whereas low-average river-flows promote infilling as seen in the morphological change during 2007–2008 and 2011–2012, and in the mooring data. This implies that it is not possible to assess whether a river is in morphodynamic equilibrium on the basis of short-term measurements. Furthermore, the infilling trend noted in the mooring data under low-average river flows may also be altered by the time-history of river flow. For instance, the 2009 mooring data were collected during a prolonged period of infill (mostly regime I) following a regime III event in April 2007 (see Fig. 2). Measurements during similar river flows during other periods may show a slightly different picture, primarily for the river flow thresholds separating the three regimes. Therefore, the analysis of sediment dynamics in such systems also needs to consider the long-term response of the system (morphological status relative to the long-term equilibrium morphology), river flow at given time, and the time history of river flow in the system (for example, occurrence of regime III events which may promote regime I behavior during subsequent periods).

In addition to its morphodynamic relevance, the conceptual picture outlined in Fig. 15 also indicates the relevance of the fluff layer for estuarine sediment dynamics (note that this point is mainly relevant for transport associated with barotropic circulation). Sediment dynamics (erosion and deposition) are mostly restricted to the fluff layer during regimes I and II which persist about 95% of the time landward of RM 8, and about 99% of the time landward of RM 1.4. In other words, fluff layer dynamics dominate sediment dynamics in the system the majority of the time. Therefore, reproducing the spatial and temporal trends in SSC and SSF in a numerical model under such conditions is dependent on the appropriate parameterization of the erodibility of the fluff layer.

It should be noted that there is some uncertainty associated with the river flow thresholds characterizing the various regimes. The uncertainty derives from potential cross-channel variations in flow rates (the mooring data represents measurements at a single location in the cross-section), and variability in the ABS-SSC relationships used to estimate SSC, and subsequently, SSF at the mooring locations. Such artifacts mainly affect the river flow threshold separating regimes I and II rather than the conceptual definition of these regimes and relevant transport processes. Another uncertainty described previously relates to the use of empirical SSC measurements, and the limitation of the SSF decomposition approach in accurately ascribing SSF to the various flow components. However, the various lines of evidence (SSF dynamics, along-channel transects, morphological changes measured in the bathymetry surveys, sediment erodibility measurements, and the numerical model) tend to a single conceptual picture of sediment dynamics. Therefore, the overall conclusions from the analysis of sediment dynamics and the conceptual picture of sediment transport presented here are valid.

6 Conclusions

Estuarine sediment dynamics reflects a complex interplay between barotropic, estuarine, and fluvial forcings. The relative importance of these forcings can vary spatially along the length of the system as well as laterally (in relatively wide systems). Sediment dynamics during low river flow conditions are primarily influenced by the formation and persistence of a thin layer of easily erodible sediments termed the fluff layer, deposited around slack water and eroded during the following flood or ebb phase of the tide. In the LPR, for the system as a whole, fluff layer dynamics persist and are relevant for sediment dynamics up to about 99% of the time. The net direction of sediment transport during low river flow is typically up-estuary, induced by barotropic lag effects and tidal asymmetries, and estuarine circulation. Increasing river flow causes an increase in currents and shear stresses during ebb, resulting in net erosion and flushing of the sediments contained within the fluff layer (up to a certain river flow threshold, estimated as 200 m3/s for the LPR), and net export of sediments from the estuary. Only when river flow exceeds 200 m3/s is erosion expected to scour the more consolidated bed underneath the fluff layer and cause estuary-wide erosional impacts. Depending on the magnitude of the event, the impact of these extreme events may persist and influence sediment dynamics for an extended period (on the order of over a year for the events examined in this paper). The analyses presented in this study allow for the classification of a given estuary at any given time into one of three regimes—importing, exporting-flushing, and exporting-scouring, in decreasing order of temporal prevalence. The morphodynamic response and categorization of the system into one of these three regimes at any given time also depend on the long-term morphological progression, river flow rate at given time, and time-history of river flow. The exporting-scouring regime represents an important aspect of the hydro-sedimentological behavior of estuaries and likely represents the conditions responsible for maintaining long-term morphodynamic equilibrium in estuaries.