Erodibility data, along with a suite of sediment and water column parameters, were collected from cores over 160 times in the York River—Pamunkey River system (Figure 1) between 2005 and 2019. A GOMEX box core (surface area 625 cm²) was collected at each location, and then multiple 10-cm diameter sub-cores were collected by hand from the box core. These sub-cores were used for analysis of erodibility and for determining disaggregated sediment grain size, resilient pellet content, organic content, and water content. For this project, parameter values from only the top centimeter were used. Disaggregated fractions of sand (>63 μm), silt (63 to 4 μm), and clay (<4 μm) by mass were determined by sonification followed by a combination of wet sieving and pipette analysis. Starting in 2010, resilient fecal pellet abundance was determined by wet sieving untreated sediment following the methods of Kraatz (2013). This involved first using a 63-μm sieve to isolate sand plus resilient sand-sized mud pellets and then subtracting out the previously calculated disaggregated sand component. More details regarding steps in the field data collection and laboratory analysis are provided in Wright (2021).

Sediment erodibility at each location was evaluated using a Gust microcosm (e.g., Dickhudt et al., 2009). During erodibility analysis, a sediment core is topped with a specially designed cap that is equipped with a rotating disk. The disk spins at varying speeds to create known stresses on the sediment bed while site water is pumped through the core. The water pumped out of the core is analyzed for turbidity and then filtered for suspended solids, and the calibrated turbidity time series is used to determine the mass of sediment eroded at each applied stress (Figure 3A). After a 30-min initial interval at 0.01 Pa, the stress is increased stepwise from 0.02 to 0.56 Pa, each lasting 20 min. These stresses are within the range of natural shear stresses observed in the York River in response to typical tidal and fluvial processes (Figure 3B).

Water content of each sediment sample by mass was determined by calculating the difference in mass from wet weight and dry weight after being thoroughly dried in an oven at 103–105°C. Due to most samples having high contents of mud, the percent water of the mud portion [100 × water mass/(water + dry mud mass)] was used instead of percent water of the total sample, similar to analysis performed by Dickhudt et al. (2011). The logic behind this approach is that at high mud content (when interlocking sand grains are not supporting the weight of the sediment), the most relevant role of water is to indicate the compaction of the mud matrix, not the compaction of the sand. The samples used for percent water were then muffled at 550°C to determine percent organic based on loss on ignition.

In the majority of cases, a 11.5 cm by 2.5 cm rectangular sub-core was also collected for X-radiography, and the resulting X-ray images were manually sorted into two categories. One category included samples that had a clearly mottled fabric (Figure 4A), and the other category included those that were laminated (Figure 4B) or that had a distinct low-density deposit layer at the surface (Figure 4C). Because 1) there were individually far fewer cases that were either laminated or that had a detectable surface layer relative to those that displayed mottling up to the surface, and 2) the hypothesized meaning of lamination or layering was similar (i.e., recent net deposition without time for extensive bioturbation), the laminated and surface layer categories were combined into one “layered” category. “Layered” vs. “Not layered” (i.e., mottled) was then included as an explanatory variable in model runs.

It is important to note that on each sampling cruise, a Gust erodibility sub-core was typically collected from each of two consecutive box cores, and the two sub-cores were then each analyzed by the Gust microcosm. The two box cores were usually collected on the order of a few to 10 s of meters apart, as dictated by the boat’s movement around its anchor point, typically close to slack tide. Two sediment property sub-cores were taken (one from each box core), and sediment from corresponding depths from both cores were combined for an average for the location. Because of this, many of the Gust experiment data have independent erodibility measurements but identical grain size, percent organic, percent water, and pellet content for a single sampling station. Limitations of sampling approaches in this study with regards to the interpretation of statistical analyses are further addressed in the discussion section.

In the depth-limited erodibility extreme, which has been found to apply very well to the muddy bed of the York River estuary (Rinehimer et al., 2008; Cartwright et al., 2009; Dickhudt et al., 2011; Fall, 2012), the sediment’s critical stress profile (τ_c ) alone is sufficient to fully quantify its erodibility. As surficial sediment is removed during active bed erosion in estuaries like the York, the magnitude of τ_c at the newly exposed, underlying bed surface typically approaches the value of the bed stress (τ_b ) imposed by the overlying flow within ∼15 min (Sanford and Maa, 2001). This means that at the end of each of the 20–30 min periods of constant stress utilized by the Gust microcosm on cores from the York as τ_b was increased incrementally from 0.01 to 0.56 Pa, one could assume that τ_c ≈ τ_b . The corresponding mass of sediment eroded at each of these levels of τ_c ≈ τ_b was known from integrating the flow rate through microcosm times the recorded time-series of suspended sediment concentration. Thus, the vertical profile of τ_c as a function of eroded mass was successfully quantified.

The next step was to condense the information provided by the τ_c profile at each site (i.e., seven eroded masses at seven corresponding values of τ_b ) into just two values, one representing the overall magnitude of erodibility, and the other representing the shape of the τ_c profile. Each of these values was used as a response variable in subsequent statistical models. The first response variable was chosen to be the mass of sediment that had been eroded at a single representative bed stress common to all of the Gust microcosm experiments (and also characteristic to typical tidal flows in the York). The second response variable was obtained by using the dominant component of a principal component analysis to describe the erosion profile shape.

As noted above, the first erodibility response variable was the mass of sediment that had been eroded at a specific applied bed stress. In this study, the cumulative eroded mass at 0.2 Pa was interpolated and used as the characteristic eroded mass value for each Gust erosion experiment (Figure 3A). Based on observations of τ_b from acoustic Doppler velocimeters (ADVs), previous investigators have identified τ_b = 0.2 Pa as a representative tidal stress for the York River estuary (Friedrichs et al., 2008; Fall, 2012; Kraatz, 2013). For example, Fall (2012) found that during the summer of 2007 (Figure 3B), τ_b during that time often reached ∼0.2 Pa during peak flood and ebb flows, with a few maximum stresses reaching ∼0.5–0.6 Pa. Similar magnitudes for τ_b can be found for each tide based on the analytical solution for the amplitude of cross-sectionally averaged tidal velocity (U) in a long, straight channel with weak friction (Friedrichs, 2010): U = ( R T / 2 ) ( g / 〈 h 〉 ) 0.5 (2) where R_T is the observed range between each low and high tide, g = 9.8 m/s², <h> is spatially- averaged water depth, and τ_b = C_D U ². For the York, <h> ≈ 6 m (Nichols et al., 1991; Friedrichs, 2009), and C_D ≈ 0.0012 based on the ADV measurements of Fall (2012) (Figure 3B).

In channelized tidal estuaries such as the York-Pamunkey system, the amplitude of tidal shear stress associated with a given tidal range is often relatively uniform for long distances along the main tidal channels. For example, observations by Cartwright et al. (2009) found a similar range of τ_b at both Gloucester Point and Claybank during several months of ADV deployments, and 3-D modeling by Fall et al. (2014) calculated relatively similar τ_b in channels along the York in both downstream and upstream locations. Thus, tidal range squared at a given time is expected to be approximately proportional to the amplitude of τ_b within channels along much of the system. Morphodynamics favors relatively uniform peak tidal velocity and τ_b along channelized tidal estuaries in general, and this pattern has been noted along additional systems such as the Delaware, Gironde, Hudson, James, Thames, and many others (Friedrichs, 1995; Lanzoni and Seminara, 2002; Olabarrita et al., 2018).

Erodibility was also quantified by the shape of the eroded mass profile. Profiles for each Gust experiment were compiled (Figure 5) and analyzed using principal component analysis (PCA) to determine the shape of the most common deviation for the mean profile. Eroded masses for τ_b = 0.01, 0.05, 0.1, 0.2, 0.3, 0.45, and 0.56 Pa were interpolated for each profile produced by the Gust experiments (Figure 5A). When profiles did not reach τ_b = 0.56 Pa, but the final experimental stress point exceeded 0.50 Pa, the last data point was linearly extrapolated. All Gust erosion profiles were divided by the mass eroded at τ_b = 0.56 Pa, resulting in a normalized erosion profile shape constrained between 0 and 1 (Figure 5B). The normalized values were then averaged across all the profiles, and the average normalized profile was subtracted from the individual normalized profiles to center the data around a global mean of zero (Figure 5C). Finally, the deviations from the average normalized profile were used as input to the PCA. Similar PCA methods for describing profile shapes have been used for the topography of tidal flats and are described in Bearman et al. (2010). Note that for this analysis, some Gust experiments did not have applied τ_b exceeding 0.45 Pa. These profiles were removed (see Wright et al., 2021 for a list of the removed profiles), which created a slightly smaller data subset for this response variable relative to the set for eroded mass at τ_b = 0.2 Pa. Figures 5D–F display the output of the PCA analysis, which is described in the results section.

Observed and predicted tidal ranges were downloaded from the National Oceanic and Atmospheric Administration (NOAA, 2021) for the tide gauge located at the Yorktown Coast Guard pier (see Figure 1) to provide a variable proportional to the amplitude of the tidal current. Tidal ranges are defined as the absolute difference between each high tide and its preceding low tide and vice versa; so, each day includes an average of 3.9 tidal ranges. The observed tidal ranges were then squared to transform that variable into one proportional to the amplitude of the tidal bed stress. Time series of tidal range squared prior to each box core collection were used to produce running averages of progressively increasing length, preceding the timing of each associated erodibility measurement. The preceding running average length that produced the highest correlation between past tidal range squared and eroded mass at 0.2 Pa was then used for calculating the tidal range squared explanatory variable. In addition to observed tidal ranges, the absolute difference between NOAA predicted and observed tidal ranges was recorded, which here is termed the water level anomaly. This variable was used to determine if meteorological events that drive relatively rapid changes in set-up or set-down of water level along the estuary influence sediment erodibility. The water level anomaly was likewise put through a series of preceding running averages of increasing length, and the case with the best correlation with eroded mass at 0.2 Pa was retained as the water level anomaly explanatory variable.

River discharge and salinity were also compiled for use as explanatory variables. Changing discharge and salinity are expected to serve as proxies for the movement of temporary, highly erodibility mud deposits which migrate, up- and down-river in concert with the ETM and STM. Daily river discharge data were downloaded from the United States Geological Survey (USGS, 2021) and summed for the two main tributaries of the York River: the Pamunkey and Mattaponi. An approach similar the treatment to tidal range and water level anomaly was used to create an optimal running average of river discharge prior to sampling of eroded mass. Salinity data for erodibility samples collected at Claybank and Gloucester Point were downloaded from the Virginia Estuarine and Coastal Observing System (VECOS, 2021) stations (see Figure 1) that are also located near Claybank and Gloucester Point. Because short-term salinity variation was high due to tidal advection, the previous 25 h before sampling (the approximate length of two M₂ tidal constituent cycles) was averaged for each set of VECOS salinity readings and used for data analysis. For sample sites not located in close proximity to the VECOS buoys (only 26 out of ≥160 erodibility measurements), depth-averaged salinity collected at the coring site at the time of sampling was used as a best available estimate of the recent characteristic salinity.

Different data subsets with increasingly more common attributes (and smaller sample size) were created and were explored with multiple statistical methods (Table 1). The entire data set included 167 erodibility profiles collected at various locations (see Figure 1) along the York River system (which here includes the lower Pamunkey). The means, standard deviations, and ranges for the entire data set are provided in Table 2. Two of the 167 data points were found to be extreme outliers in all models over the course of model exploration and were removed from further model-fitting. [A complete listing of all data points and their individual properties are provided in Wright et al. (2021)]. The first data subset analyzed statistically included all the remaining observations (n = 165). Most samples were located below the ETM except those in the Pamunkey River. The next data subset removed the Pamunkey River samples and only contained the fully “estuarine” sites (i.e., sites located seaward of the transition to fresh water) (n = 158), the logic being that distinct physical processes (such as density-driven estuarine circulation) occur in brackish conditions.

The two sampling locations that were the most consistently visited over the 15-year sampling period were Claybank and Gloucester Point. Thus, these two locations were combined within one subset (n = 136). Previous studies (Dickhudt et al., 2009, 2011; Fall, 2012) focused on these two areas because of their contrasting degrees of physical bed disturbance and deposition frequency as revealed in X-rays. Due to a large percentage of samples being within the Claybank region alone, these samples were also analyzed separately (n = 110). The logic behind the focus on Claybank and Gloucester Point and then on Claybank alone is the larger number of quasi-replicates available to strengthen detection of meaningful statistical trends relative to other inherent environmental “noise.” Most, but not all, of the Claybank and Gloucester Point samples include X-ray classifications. Because near surface layering could be an important driver of erodibility, both the Claybank and Gloucester Point subset and the Claybank subset were divided further, with each subset containing only samples that included an X-ray classification (n = 104 and n = 82, respectively).

Many samples from 2010 to 2019 also included data for fecal pellet abundance. Kraatz (2013) hypothesized that pellet abundance could alter erodibility, and therefore various fecal pellet subsets were initially included in the analysis. However, none of the statistical models considered produced a linear dependence of erodibility on fecal abundance in a manner that was consistent with the simultaneous linear effects of other variables. Thus, fecal pellets were dropped from further modelfitting. Possible non-linear associations between fecal pellets and erodibility are revisited in Section Variables That Could Also Have High Impact but Were Not Included.

Data set familiarization techniques included histograms and scatter plot matrices. Histograms of eroded mass at 0.2 Pa appeared log normal, and eroded mass at 0.2 Pa was log-transformed for subsequent analysis. River discharge was also log-transformed because it is common practice in the hydrological literature (e.g., Iddrisu et al., 2017). All variables were then linearly standardized by subtracting their means and dividing by their standard deviations. To allow standardization, “mottled” and “layered” X-rays were first assigned values of 0 and 1, respectively, before being standardized. Scatter plot matrices were then used to examine the size and sign of correlations between the explanatory and response variables and to assess co-linearity between the explanatory variables themselves (Figure 6). Any two explanatory parameters with high correlation values were monitored if included in the same model, and variance inflation factor (vif) values were evaluated for each model as a test for severe collinearity. If strong collinearity existed (common threshold of vif values >5) (James et al., 2013), each explanatory variable was tested separately in similar models and only one of the pair was retained.

Two different model sets (Table 3) were used to evaluate the two erodibility response variables identified for this study: 1) eroded mass at τ_b = 0.2 Pa and 2) the shape of the erosion vs. τ_b profile as determined by its principal shape component. The goal of the first grouping of explanatory variables (the “Eroded Mass Model Set”) was to determine the relative importance of various sediment properties and hydrodynamic conditions in determining erodibility at a representative, commonly occurring bed stress within the York River. The goal of the second variable grouping (the “Erosion Shape Model Set”) was to determine which sediment and hydrodynamic factors may be important in determining the overall shape of the τ_c versus eroded mass profile. Explanatory variables that were considered for this second model set were identical to those in the previous model set, but also include eroded mass at 0.2 Pa.

Multiple linear regressions were run for each data subset within each model set (see Table 3 for model framework). The general equation for multiple linear regression is: y = β 0 + β 1 X 1 + β 2 X 2 + ⋯ + β q X q + ε (3) where y is the predicted response variable, β ₀ is the model intercept, X ₁ though X_q represent individual explanatory variables, and β ₁ through β_q represent the corresponding best-fit regression coefficients. Lastly, ϵ is the model residual, i.e., the component of y not reproduced by the other terms on the right-hand-side. Due to the data being standardized prior to model formulation, the absolute value of the β coefficients can be ranked to show which explanatory variables resulted in the largest change in the response variable. For instance, for a β coefficient of 0.3, for every increase of one standard deviation of the explanatory variable, the response variable will increase by 0.3 standard deviations.

All possible combinations of the explanatory variables for each data subset were used in assessing Eq. 3, and the relative likelihood of any given model being the best among the available choices for that data subset was assessed using the Corrected Akaike Information Criterion (AICc). A lower AICc value indicates that a model is more likely to be best, and a model’s ΔAICc value is defined as its own AICc minus that of the model with the lowest AICc value in its model subset. In this study, all models with a ΔAICc value ≤2.0 units are considered comparable models with useful information regarding trends (Burnham and Anderson, 2004). Because ΔAICc values are relative within a given model subset, they cannot be used to compare the explanation of variance of response variables across model subsets. To remedy this, adjusted R ² values were included to describe how much of the data variation was explained by the combination of parameters in each model across sub-model categories.

Application of running means indicated that the average of the past 33 tidal range squared observations (8.5 days—i.e., centered around conditions 4.3 days prior) had the highest correlation with eroded mass at τ_b = 0.2 Pa. The best correlation between the water level anomaly and eroded mass was found to be over the past 27 cycles (7.0 days). For river discharge, the running average with the highest correlation with eroded mass at 0.2 Pa was over the previous 170 days (i.e., centered around conditions 2.8 months prior).

PCA for Gust profiles was completed to reduce the shape profile containing seven eroded mass values and seven stress values, to one shape score that could be used as a response variable in multiple linear regression. Principal component one (PC1) (i.e., the shape score) explained approximately 83% of the variation of shape within all data subsets (Figure 5D). Positive and negative values of PC1 made the normalized eroded mass profile less and more concave, respectively, than the average profile shape (Figures 5E,F). Negative shape scores corresponded to a geometrically increasing shape profile, such that each increase in stress resulted in a proportionately greater increase in eroded mass. In contrast, positive shape scores tended to indicate a more linear profile or, when most positive, a proportionally smaller increase in eroded mass with each increase in stress. PC1 was the only component whose scores were used as a response variable in the models reported here because PC1 explained ∼83% of the variability in erosion profile shape in all data subsets. PC2 described only ∼11–13% of variation in the erosion profile shapes, and PC3 accounted for just ∼3–4% (Figure 5D).

Initial assessment of relationships among variables was facilitated by correlation plot results (Figure 6) applied to each of the six data subsets listed in Table 1 and Table 3. (The X-ray layering parameter is not part of the All Sites data set in Figure 6 because X-rays were not collected at every site.) A complete set of correlation plot results for all six data subsets are included in Wright et al. (2021).

Despite its dominant role in the erodibility literature (e.g., Mehta, 2014), percent water of mud was consistently found to be poorly or nonsensically correlated with erodibility, and percent water was ultimately dropped from all multiple regression models for all six data subsets. For example, for the All Sites data subset, percent water of mud had the poorest correlation with eroded mass at 0.2 Pa of any variable, with r = −0.002 in the Eroded Mass Model Set (Figure 6). Across the other five data sets, the average magnitude of the correlation between percent water of mud and eroded mass was likewise quite poor, with a mean |r| of 0.111 (Wright et al., 2021). Furthermore, when included in multiple regressions with the other variables, the few models for eroded mass with ΔAICc < 2.0 that retained water content produced a negative β value for percent water (Wright, 2021). This is physically nonsensical given that there is a rich and long-established literature (see Mehta, 2014 and references within) demonstrating that erodibility for fine sediment should be negatively correlated with sediment bulk density (and thus positively correlated with percent water). Possible explanations for the unexpected result for the role of percent water are discussed in Section Limitations Associated With Sampling Approach and Resolution of Bed Properties.

Collinearity among explanatory variables was also assessed using correlation plot results. Among the All Sites explanatory variables that were further considered, salinity had the highest correlations with other variables: namely, with distance upriver (at r = −0.701) and with river discharge (at r = −0.532) within the Erosion Shape Model Set (Figure 6). This is not unexpected, since salinity along the York River estuary is predicted well by a nonlinear model based only on position along the estuary and past river discharge (Parrish et al., 2019). Thus, for the All Sites linear models, salinity was dropped from further consideration. However, for the other data subsets, salinity turned out to be a better predictor of the response variables than river flow, so salinity was utilized in its place. Across all data subsets, which include 452 explanatory variable correlations, only three other cases (not involving river discharge or salinity) had correlation values over 0.6. However, in the individual candidate and global models for those data subsets, vif factors were all below 5.0, so the associated variables were retained.

There were 33 candidate models identified (i.e., models with ΔAICc ≤2.0) within the six data subsets for the Eroded Mass Model Set. Overall, tidal range squared and salinity/river discharge were the most consistent explanatory variables throughout the entire model set, with these variables being retained in 100% of the models in which they were considered. Figure 7A shows average |β| values for each individual data subset, with each subset corresponding to a given bar color. [All individual β and adjusted R ² values for all 33 models are provided in Wright et al. (2021)]. The variables with the highest average |β| values across all of the Eroded Mass Model subsets included past tidal range squared (0.424), salinity (−0.294) / river discharge (0.252), and percent organic (0.226). Variables that were less important, but still had consistent trends in average β were past water level anomaly (0.171), percent sand (−0.143), distance upriver (0.136), percent clay of mud (−0.135), and presence of X-ray layering (0.108).

Model performance tended to increase as data subsets became more spatially focused and otherwise specific. Partial residual plots are shown in Figure 8 for representative 5-component models from the All Sites, Claybank & Gloucester Point, and Claybank with X-rays subsets. The slope of the line in each panel is the β value for the corresponding explanatory variable within the model. The β value for the most important explanatory variable (tidal range squared) (Figure 8) increased as the data sets became more localized, and the overall performance of the models with ΔAICc scores of 0 also tended to improve (Appendix Figure A1). As presented in detail in Wright et al. (2021), the multiple All Sites subset models had an average adjusted R ² of 0.392, which increased to an average of 0.464 in the Estuarine Sites subset. Furthermore, both Claybank & Gloucester Point subsets (with and without X-rays) had an average of 0.585, and both Claybank subsets resulted in an average adjusted R ² of 0.609.

The model set for profile shape score consisted of 39 candidate models throughout the six different data subsets. Throughout all data subsets, eroded mass at 0.2 Pa was by far the most prominent explanatory variable for all candidate models, being included in 100% of candidate cases and exhibiting an absolute value for β twice as large as any other variable (Figure 7B, Figure 9). In every model, the relationship between eroded mass and profile shape was positive, such that greater erodibility favored a straighter (or less concave-up) profile (see Figure 5F). Following eroded mass at |β| = 0.519 (averaged across all data sets), the next highest |β| values were for organic proportion at 0.249, water level anomaly at 0.182, tidal range squared at 0.171, and percent sand at 0.166. Six other variables had moderately low average |β| values, ranging between 0.136 and 0.058.

Given the strong positive relationship between profile shape score and eroded mass at 0.2 Pa, one might initially expect the remaining explanatory variables to show relationships with shape score such that the sign of their β values matched the signs relating them to eroded mass in the previous section. This was indeed the case for percent organic, water level anomaly, tidal range squared, river discharge, and salinity. However, the signs for the β values relating percent sand, percent clay of mud, X-ray layering, and distance upriver to profile shape score were all opposite to signs of the β values relating them to eroded mass.

As was the case for eroded mass at 0.2 Pa, the models for profile shape score also tended to increase in performance as data subsets became more spatially focused, although not as dramatically (Figure 9 and Figure A2). Partial residual plots are shown in Figure 9 for representative 6-component models from the All Sites, Claybank & Gloucester Point, and Claybank with X-rays subsets. The All Sites subset had an averaged adjusted R ² of 0.537, which slightly decreased in the Estuarine Sites subset to 0.517 (see Wright et al., 2021). However, adjusted R ² averages increased again in the Claybank & Gloucester Point subsets and the Claybank subsets to 0.581 and 0.703, respectively.

Using the top-performing (ΔAICc = 0) model from each model subset (Blue circles and blue dashed lines in Appendix Figure A1 and Figure A2), a validation exercise was performed to test the general stability the model fits. Each data subset was split into two halves, with one half used to calibrate new best-fit β values. The new β values were then used to predict the response variable for the validation half (red squares and red dashed lines in Appendix Figure A1 and Figure A2). To help ensure that each half encompassed a similar range of variable values, the calibration half consisted of the first and fourth quarters of each subset in time, and the validation set consisted of the second and third quarters. This was done because the lower- and upper-most estuary had mainly been sampled during the first and second halves of the 15-year time record, respectively.

The adjusted R ² values for the validation data fits were generally similar to the adjusted R ² values for the best model sets (Figure A1 and Figure A2), which demonstrates that the model fits were reasonably stable. In several cases, adjusted R ² was larger for the validation set than the corresponding full data set, indicating that some of the validation data subsets were less noisy than their corresponding full data sets. The regression lines for the validation sets (red dashed lines), however, were almost always somewhat farther from the 1:1 line relative to the full data sets (blue dashed lines), as would be expected.

Two of the explanatory variables with the strongest positive effects on eroded mass at 0.2 Pa across all data subsets highlight the influence of recent physical bed disturbance on bed erodibility. Past tidal range squared, the most influential variable of all, represents the role of recent tidal bed stress in repeatedly suspending and depositing mud, such that the newly disturbed sediment is unconsolidated and easier to subsequently entrain into the water column again. Past water level anomaly, or the absolute difference between observed and predicted tidal range, is interpreted here to represent the role of recent winds causing set-up or set-down of water level, which is also associated with winddriven currents that further physically disturb the seabed, enhancing erodibility.

Several explanatory variables also emphasize that recent deposition increases erodibility as measured by eroded mass at 0.2 Pa. When considering observations from a location that is seaward of the ETM, such as data from Claybank, a decrease in salinity indicates the continued progression of the pool of easily suspended sediment down the estuary and net deposition of mud. Thus, erodibility increases at Claybank as salinity decreases. An increase in percent organic matter in the upper 1 cm of the bed also suggests recent deposition of “fresher,” easier to resuspend muddy flocs. So, erodibility increases with percent organics. Likewise, the presence of layering at the surface as seen by X-radiographs indicates recent deposition, and layering was found to be positively correlated with erodibility. The past water level anomaly may also be associated with recent deposition events in the channels where most of the erodibility measurements were collected. Episodic wind waves associated with the water level anomaly might winnow mud from the shoals, which then might produce temporary deposition and higher erodibility in the neighboring channels.

Multiple explanatory variables highlight the role of time since deposition (i.e., consolidation time), the effect of erosion exposing older sediment, and/or the past winnowing of fines decreasing eroded mass at 0.2 Pa. At a given site below the ETM, an increase in salinity indicates the likely net movement of mud upstream and away from the site, exposing older sediment and decreasing erodibility. As “fresher” mud migrates away, one would expect lower percent organics and higher percent sand, trends which are each seen to be associated with reduced erodibility in all models in which they appear. Likewise, the negative association between mottled X-radiographs and eroded mass at 0.2 Pa indicates that older sediment (i.e., that which has had time to be fully bioturbated) is less erodible. A lower fraction of silts (as indicated by a higher clay fraction of mud and/or higher sand fraction overall) is associated with lower erodibility in the York. As previously recognized by many investigators, fine sediments tend to become harder to erode as their clay-sized grain content increases because of an increase in grain-to-grain cohesion as grains become smaller (e.g., Postma, 1967; Dade et al., 1992; Roberts et al., 1998; Mehta, 2014). In contrast, mixtures containing sand-sized sediments eventually become harder to erode as their sand content increases, because suspension of non-cohesive sediment is more difficult as grain size increases.

An important finding in this study with regards to consolidation of the uppermost seabed is that the strong increase in erodibility due to recent physical disturbance (as parameterized by tidal range squared and water level anomaly) only lasts several days. For both tidal range squared and water level anomaly, the average of the previous ∼7 to 8 days of observations—i.e., a time period centered around conditions ∼4 days prior—had the highest correlation with eroded mass at 0.2 Pa. This timescale of a few days is consistent with consolidation times scales reported for high water content muds in laboratory studies (Mehta, 2014 and references within), with those found to work well in modeling studies (Rinehimer et al., 2008; Sanford, 2008; Sherwood et al., 2018), and is also close to the time scale for the transition from spring to neap tides (as likewise noted by Mehta, 2014). Similar results have also been seen in the field by others. Based on a weekly time-series of five cores from the York River estuary, Kraatz (2013) previously found that eroded mass at 0.2 Pa was best correlated to the tidal range averaged over the previous 5 days, and Huang et al. (2020) found that it took approximately 7 days for sediment layers to consolidate after being deposited by tides in the Pearl River Estuary.

Throughout the Eroded Mass Model Set, the hydrodynamic explanatory variables had a statistically stronger relationship to eroded mass at 0.2 Pa than did the sediment property variables, according to average |β| values. However, hydrodynamic variables influence erodibility indirectly. Physically, erodibility is associated more directly with the properties of the bed, as seen in laboratory consolidation and erosion experiments focusing on specific sediment grain sizes and bulk densities (Roberts et al., 1998; Mehta, 2014). Nonetheless, for the Eroded Mass Model Set, the parameter with the highest average |β| value was past tidal range squared (0.424; Figure 7A). Past tidal range squared was present in all candidate models, and it had the highest |β| within all 33 cases. The mechanistically similar parameter of past water level anomaly was also present in every subset of models. These two hydrodynamic parameters, which are interpreted as proxies for physical bed disturbance and reduced consolidation, predicted erodibility more successfully than most directly measured sediment properties. The second most important variable of all, salinity at an average |β| of 0.294, is also a hydrodynamic variable, as is river discharge at third most important (average |β| = 0.226).

The overall importance of hydrodynamic variables may have ramifications for enhanced performance of bed erodibility routines within computationally-intensive, multiparameter numerical models, such as 3D hydrodynamic codes that include biogeochemistry. It is easier to accurately model hydrodynamic environmental variables than to reproduce centimeter or millimeter-scale sediment bed properties. By applying these empirical relationships between hydrodynamics and erodibility in the York or by developing similar relationships for other estuaries, modelers may be able to produce a more accurate representation of sediment bed critical stresses and associated eroded mass in given areas of an estuary by only having to include the effects of external variables such as recent velocities (or bed stresses), salinity, and/or river discharge. These parameters are likely more straightforward to model than having to precisely reproduce sediment properties throughout the estuary, especially when they are spatially heterogeneous. If it is possible to represent bed grain size and organic matter well in the estuarine numerical model, a multiple regression to predict erodibility could also use these variables to further improve erodibility results.

The Erosion Shape Model Set included an additional sediment property, namely eroded mass at 0.2 Pa, as an explanatory variable, which ended up having the largest influence on shape score by far. Eroded mass at 0.2 Pa had the highest average |β| in all the candidate models across all Erosion Shape Model subsets, with an average of 0.519 which was more than twice as large as |β| for the next most important erosion shape variable (Figure 7B). Higher eroded mass at 0.2 Pa tended to make the shape score more positive, which indicated a more linear eroded mass profile (blue PC1 > SD lines in Figures 5E,F). Conversely, beds with low erodibility tended to have a more strongly curved, concave-upward relationship between applied stress and total eroded mass (red PC1 < -SD lines in Figures 5E,F). Although they did not comment on this shape trend, Dickhudt et al. (2009) also observed more strongly curved profiles at low stress and an increase in profile linearity at higher erodibility. Dickhudt et al. (2009) fit a curve of the form τ_c ∼ m^b plus an offset to groupings of their York River erodibility profiles, where m is cumulative eroded mass. Values of b closer to 1 indicate a more strongly linear fit, whereas smaller b indicates a more concave fit. For their “low,” “transitional,” and “high” erodibility groups, Dickhudt et al. (2009) found b = 0.51, 0.65, and 0.75, respectively. An analogous fit to the lower erodibility (red), intermediate (black), and higher erodibility (blue) curves from the present study (Figure 5F) yields b = 0.29, 0.60, and 0.89, which is a consistent trend.

The low erodibility, concave-upward shape may represent a near-steady-state erodibility profile present in the absence of significant recent bed disturbance, such as immediately following neap tide. This shape is similar to the concave equilibrium eroded mass vs. applied stress profile assumed by modelers of muddy bed erodibility in the absence of disturbance (Rinehimer et al., 2008; Sanford, 2008; Sherwood et al., 2018). The transition towards a more linear erodibility profile in the presence of more frequent recent disturbance (and higher overall erodibility), such as immediately following spring tide, reflects a relatively larger increase in incremental eroded mass in the range of ∼0.2 to 0.4 Pa than at higher or lower applied stress (Figure 5F). This makes sense, since the frequency at which the τ_b = 0.2 to 0.4 Pa is reached in the York notably increases during spring tides (see Figure 3B). Stresses larger than 0.4 Pa are not reached very often, even on spring tides, and stresses less than 0.2 Pa are reached on nearly every tide, so spring-neap variability in bed disturbance does not change erodibility as dramatically at either of these τ_b end members.

Similar to the magnitude of erodibility, the erosion shape profile can potentially be predicted by numerical models in some cases based solely on hydrodynamic forcings. As discussed in Section Hydrodynamics Can Serve as a Proxy for Effects of Bed Properties on Erodibility, eroded mass for very muddy beds can be predicted relatively well based on recent hydrodynamic conditions. In turn, eroded mass is the most influential predictor of erosion shape, with |β| values twice that of any other explanatory variable. Thus, modelers can constrain both the magnitude and profile shape of the critical stress in very muddy sediments by focusing largely on the hydrodynamic history of a study area. In this respect, the findings of this study are highly consistent with the conceptual approach of Sanford (2008), who included a mud-only case for which changes in erodibility in time were predicted based only on a time history of tidally varying applied bed stress.

This study was largely an opportunistic analysis of an extensive sediment erodibility data set that was not collected for a single purpose. The samples were collected over 15 years at differing temporal and spatial scales for various projects. Due to clustered sampling in time and space associated with the past measurements, many standard statistical assumptions were not well adhered to. Thus, formal statistical significance levels such as p-values were not reported here for any of the multiple linear regressions. Over 75% of the samples used in this study were taken during the spring and summer months (March through August). Therefore, it is hard to definitively state that explanatory variables with seasonal changes like discharge, salinity, and water level anomaly follow the same predictive models all year round. For example, the average wind speed and freshwater discharge along the York are greater in winter than in summer by a factor of ∼1.3 and ∼3.2, respectively. Additionally, the data set was heavily skewed by the Claybank location, with over 65% of samples being in that subset. Given these sampling design limitations, the goal here was not to rigorously focus on the absolute level of the statistical significance of the results. Rather, the focus was on trends in the responses and the relative importance of one explanatory variable versus another.

In addition, the field methodology for subsampling sediment cores for water content analysis most likely did not sufficiently resolve the very surface of the sediment core. Based on a mean water content of 75.5% (Table 2), the upper 1 mm of the seabed contained, on average, 0.29 kg/m² of dry sediment. The average dry sediment mass eroded at 0.2 Pa was only 0.089 kg/m² (Table 2), and >90% of cases eroded less than 0.29 kg/m² (Wright et al., 2021). Thus, the Gust experiments, on average, eroded only ∼0.3 mm into the bed at 0.2 Pa, and >90% of Gust experiments eroded less than 1 mm at 0.2 Pa. Unfortunately, subsamples from box cores did not attempt to resolve the uppermost 1 mm, and instead the top 1 cm from sub-cores was homogenized. Quite possibly, the water content of the full 1 cm of the upper seabed was not sensibly correlated to the water content of the uppermost 0.3 to 1 mm. The same methodology was used in the Dickhudt et al. (2009, 2011) studies, and it may be a contributing factor to why they also found no meaningful relationships between the variation of water content of mud in the York River estuary and sediment erodibility. In contrast, sensible and meaningful correlations with eroded mass at 0.2 Pa were found for the clay content of mud, sand content, and organic content of the top 1 cm. A logical explanation is that vertical gradients in clay, sand, and organic content over the top 1 cm are not nearly as strong as the vertical gradients in water content.

There may be additional factors relevant to moderately turbid, muddy tidal estuaries like the York that could make the correlation between water content over the top 1 cm and erodibility noisier than that observed in well-controlled lab experiments. As also noted by Dickhudt et al. (2009), the water content of mud over the top 1 cm in the York varies much less than other potential controls on erodibility. The standard deviation of water content divided by its mean was only 0.06 (Table 2). This means that even modest errors in quantifying water content could overwhelm the signal in the variance. In contrast, standard deviations divided by the mean for the bed properties with the most observed influence on erodibility, namely organic content and sand content, were 0.55 and 1.1, respectively, meaning that their signals were more likely to exceed their noise. Another complication when averaging over the top 1 cm may be that water content at this scale is positively correlated to percent clay of mud (Figure 6), possibly because high clay content reduces mud permeability (Zhang et al., 2015). Clay content of mud, in turn, was found to reduce erodibility. So, given the relatively small variation in water content in York River mud when averaged over the top cm of the bed, variations in clay content may at least partly counteract the influence of the mud’s water content.

There are variables that were not included in this analysis that could also influence sediment erodibility and produce some of the unexplained variance remaining in the models. These include local mean water depth, position across the channel, and bedform types. Mean water depth and across-channel position were considered but were not found to notably account for variance in erodibility in this study’s data set within the context of general linear models. This may be because so much of the data in this study were repeatedly collected at a few specific sites with relatively little sample-to-sample variation in depth or lateral location. Bedform types were not used because of a lack of observations corresponding with each sample, although mud furrows have been shown to potentially play a role in physical disturbance of the bed in the York River estuary (Dellapenna et al., 1998, 2003). In general, additional explanatory variables were not available because this data set dates back to 2005. The variables used in the final model formulations are either those that were collected according to lab protocol and have stayed consistent from 2005 to 2019 or ones that were accessible from online sources for each of the sediment cores. For example, the smaller number of pellet abundance samples (58 samples) versus most other bed properties (165 samples) stems in part from it not becoming a standard variable until later in the sampling record.

Within estuaries characterized by a main channel bordered by distinct shoals, such as the York River, lateral changes in sediment properties and hydrodynamic properties may especially influence sediment bed erodibility. The northeast and the southwest shoals of the York River tend to differ in terms of grain-size distribution and the mechanisms of sediment reworking (Hinchey, 2002; Dellapenna et al., 2003; Kniskern and Kuehl, 2003). For example, the northeast shoal of the river is sandier than the southwest shoal. Ideally, the inclusion of grain-size distribution and distance upriver might be able to adequately account for these observed differences. However, Huzzey (1988) reported lateral density gradients at Claybank, with homogenous lateral densities only present during maximum tidal currents. Huzzey (1988) also states that the cross-estuarine gradients could drive significant lateral water circulation in the area, which may likely be strong enough to preferentially erode sediment and deposit in different areas of a channel cross-section, especially in such a soft-bottomed area such as Claybank. The influence of lateral circulation and resulting changes in lateral transport convergence and erodibility have been suggested for the neighboring James River Estuary (Huijts et al., 2006). In the current study, only longitudinal estuary circulation was considered in the form of river discharge and salinity affecting the presence of the estuarine turbidity maxima.

Only linear relationships were reported here, but non-linear relationships were also considered during early stages of the analysis. Generalized additive models (GAMs) were initially included, but they tended to show exaggerated, oscillatory, and unrealistically non-linear relationships with many of the explanatory variables due to overfitting. When this tendency was compensated for by increased smoothing in the GAM settings, the results then provided insights very similar to multiple linear regression, suggesting that the added complexity of GAMs was not warranted. Nonetheless, at least one potential explanatory variable, namely resilient fecal pellets, seemed to have a systematically non-linear relationship with erodibility (Figure 10). In the pellet subset, there was an initial increase in eroded mass with pellet abundance at a low proportion of pellets, but then eroded mass began to decrease again at a higher proportion of pellets. Perhaps the initial increase in pellet abundance is acting as a proxy for decreasing clay content of mud, which would increase erodibility. Eventually, continued increases in pellet abundance may be a proxy for behavior similar to coarser non-cohesive sediment, such as sand, which would decrease erodibility once more. It could be beneficial to better explore these non-linear relationships in further investigations of this or other similar datasets.

The current study would have likely explained a greater amount of variance within the dataset had it included more accurate measurements for water content in the upper few mm of the seabed. The upper few mm of very soft muddy sediments, like those found in the York River estuary, can be difficult to sample. Perhaps higher resolution methodologies for measuring water content, such as resistivity profiling (Wheatcroft and Borgeld, 2000), fiber optical backscatter profiling (Hooshmand et al., 2015), or freeze coring followed by thin slicing (Harrison et al., 2016) could be used in conjunction with future sediment sampling to better determine the true fine-scale influence of this parameter on sediment erodibility.

Finally, future statistical studies of controls on bed erodibility in muddy estuaries would also benefit from improved field sampling in terms of experimental design. Choices for future erodibility coring sites should utilize aspects of random sampling in time and space to reduce potential autocorrelation of samples and better ensure that the true range of variability in explanatory and response variables is well represented. Nonetheless, the large sample size and temporal duration of the present data set still makes it a uniquely valuable resource for better understanding controls on estuarine sediment bed erodibility, as long as its inherent limitations are recognized.