**Kelsey A. Fall, Courtney K. Harris, Carl T. Friedrichs, J. Paul Rinehimer and Christopher R. Sherwood**

**Abstract:** The Community Sediment Transport Modeling System (CSTMS) cohesive bed sub-model that accounts for erosion, deposition, consolidation, and swelling was implemented in a three-dimensional domain to represent the York River estuary, Virginia. The objectives of this paper are to (1) describe the application of the three-dimensional hydrodynamic York Cohesive Bed Model, (2) compare calculations to observations, and (3) investigate sensitivities of the cohesive bed sub-model to user-defined parameters. Model results for summer 2007 showed good agreement with tidal-phase averaged estimates of sediment concentration, bed stress, and current velocity derived from Acoustic Doppler Velocimeter (ADV) field measurements. An important step in implementing the cohesive bed model was specification of both the initial and equilibrium critical shear stress profiles, in addition to choosing other parameters like the consolidation and swelling timescales. This model promises to be a useful tool for investigating the fundamental controls on bed erodibility and settling velocity in the York River, a classical muddy estuary, provided that appropriate data exists to inform the choice of model parameters.

Reprinted from *J. Mar. Sci. Eng.* Cite as: Fall, K.A.; Harris, C.K.; Friedrichs, C.T.; Rinehimer, J.P.; Sherwood, C.R. Model Behavior and Sensitivity in an Application of the Cohesive Bed Component of the Community Sediment Transport Modeling System for the York River Estuary, VA, USA. *J. Mar. Sci. Eng.* **2014**, *2*, 413-436.

#### **1. Introduction**

Fine sediment transport in coastal and estuarine environments has significant physical, biological, and chemical ramifications. Mobilized sediments reduce water clarity, transport toxic materials, pathogens and nutrients, and fill navigational channels [1,2]. Bed erodibility and settling velocity are key parameters influencing fine sediment dynamics in coastal and estuarine environments. Bed erodibility controls the amount of sediment suspended while settling velocity influences how far it is transported [3–7]. Erodibility and settling velocity vary widely over time and space in estuaries, and these variations are closely related to sediment flux convergences and divergences at Estuarine Turbidity Maxima (ETMs) [8–12].

ETMs are regions of locally high suspended-sediment concentrations that often occur immediately landward of the salt limit, where convergence in near-bottom flow traps suspended sediment in the bottom layer [13–15]. Suspended sediment trapped in the bottom layer near the head of salt cannot be entrained into the upper layer due to damping of turbulent mixing by stratification [16]. Tidal asymmetries in vertical velocity and suspended-sediment profiles can also produce ETMs [9,17,18]. Although these processes tend to be most important near the head of salt, they can create ETMs in other areas where the channel geometry alters the salinity field, often called Secondary Turbidity Maxima or STMS [19–22]. Hydrodynamic forces and sediment and bed properties influence both ETMs and STMs, and because these factors vary with time and in all three spatial dimensions, they are difficult to study using field measurements alone. A three-dimensional model is helpful for evaluating and understanding these complex processes that influence sediment dynamics in muddy estuaries.

Erosion of sediment from the bed provides an important control on estuarine turbidity. Here we define erodibility (*İ*) as the asymptotic relationship between a steady, externally imposed bed stress (*IJb*), and total eroded mass from the seabed when *IJb* exceeds the critical stress (*IJc*) of the sediment surface. The critical shear stress, *IJc*, represents the stress at which motion or suspension of sediment first occurs [23]. Erosion rates have been estimated using many different formulas [24–27], but the simple linear Ariathurai-Partheniades erosion formulation has been assumed in many cases:

$$E = M(\mathfrak{r}\_b - \mathfrak{r}\_c) \tag{1}$$

Here, the erosion rate, *E* (kg m<sup>í</sup><sup>2</sup> s<sup>í</sup><sup>1</sup> ), varies linearly with the excess shear stress (*IJ<sup>b</sup>* í *IJc*) (Pa) according to the erosion rate parameter, *M* (kg m<sup>í</sup><sup>2</sup> s<sup>í</sup><sup>1</sup> Pa<sup>í</sup><sup>1</sup> m<sup>í</sup><sup>2</sup> ) [11,12,28–33]. Though developed for cohesive sediment, Equation (1) can also be applied to non-cohesive beds. For purely non-cohesive beds, *IJc* mostly depends on the grain size and density of individual particles, and typically increases with diameter [23]. For cohesive sediment, *IJ<sup>c</sup>* represents a bulk characteristic of the seabed, and may depend on grain size, porosity, organic content, and depositional history, and often increases with depth in the seabed and with time since deposition [24].

Many three-dimensional numerical models assume a constant *IJc*, even for cohesive sediment [7,34–36], which can produce satisfactory results when applied for short time scales, but neglects feedbacks between erodibility, erosion, and deposition that develop in response to events such as flood deposition, spring freshets, storm erosion, or biogenic seasonal variations [37–40]. A model of the York River estuary in Virginia that defined a constant value of *IJc* [35,41] was able to represent the STM but underestimated suspended-sediment concentrations, despite using 0.05 mm/s for settling velocity, significantly smaller than the values inferred from recent observations [8,10,42]. More recently, a York River model was implemented to estimate suspended-sediment concentrations associated with Hurricane Isabel [7]. This used Sanford's (2008) consolidation model [43] that allows *IJc* to vary in response to consolidation. The parameters in their bed consolidation model were based, when possible, on data from the York River, or otherwise from literature values, however the article did not discuss the sensitivity of calculations to these parameters.

In this study we applied a three-dimensional hydrodynamic and sediment-transport model to represent processes in the York River estuary, VA, USA, using the Community Sediment Transport Modeling System (CSTMS) with a cohesive bed sub-model [33]. The critical shear stress profile of the cohesive bed was estimated following the consolidation and swelling model presented in [43]. The objectives of this paper are to (1) describe the application of the three-dimensional hydrodynamic York cohesive bed model, (2) compare calculations to observations, and (3) investigate sensitivities of the cohesive bed sub-model to user-defined parameters.

#### **2. Study Site: York River Estuary, VA, USA**

The York River (Figure 1), a partially mixed microtidal estuary, spans 50 km from West Point to Gloucester Point [44]. Tidal currents approach 1 m/s at the water surface and dominate sediment resuspension [45]. The estuary contains a main channel (~5–20 m deep) and a secondary channel (~5 m deep), both dominated by mud, and bordered by well-developed sandy shoals (~2 m deep) [44,46]. The primary ETM occurs upstream near West Point. A secondary turbidity maximum (STM) occurs about 40 km landward of the river mouth in a mid-estuary transition region near the Intermediate site between the well-mixed zone upriver and more stratified zone seaward [20]. Physical processes dominate in the upper and middle York, where the benthic diversity, productivity and biomass are suppressed by the more intense tidal currents, greater range of salinities, and the ephemeral presence of the STM [44]. In the lower York, near Gloucester Point, the estuary turns toward the northwest, widens from about 3 km to 6 km, and the main channel depth increases to 20 m. As a result, the seabed is only disturbed during storms, and suspended-sediment concentrations are lower (10 s of mg/L), creating a more favorable environment for benthic organisms [45].

**Figure 1.** Map of York River estuary, southeastern VA, USA. The location of the Multidisciplinary Benthic Exchange Dynamics (MUDBED) Intermediate site is in green and MUDBED Biologically Dominated site is in blue.

The York River estuary has been studied as a part of the Multidisciplinary Benthic Exchange Dynamics (MUDBED) project [4], particularly near Clay Bank at the MUDBED Intermediate site, and near Gloucester Point at the MUDBED Biologically Dominated site (Figure 1). The seabed at both sites is dominated by mud and resilient biologically repackaged fecal pellets, with very little sand [47,48]. Individual mud particles and flocculated muds form the major fraction of suspended material in the York. Repackaged fecal pellets are suspended during periods of high bed stresses [42].

While data from the Biologically Dominated site revealed little variability in either settling velocity or erodibility, at the Intermediate Site long-term observations from the MUDBED benthic Acoustic Doppler Velocimeter (ADV) tripod [4,10] and Gust microcosm experiments [11] showed temporal variability in sediment settling velocity and bed erodibility. During the late winter and spring (February–May), erodibility at the Intermediate site generally exceeded that during the summer and fall, with the most erodible beds occurring in April and May [11]. Settling velocities based on ADV data, and bed erodibilities estimated by the Gust microcosm and ADV data were inversely correlated at the Intermediate site, with lower settling velocities observed during periods of increased erodibility, while settling velocities increased during periods of reduced erodibility [10].

Observations indicated that the temporal variability in bed erodibility and settling velocity at the Intermediate site can be attributed to the presence or absence of the STM [4,11,49]. At times of high river flow, stratification of the water column increased, and slowly settling, easily erodible flocculated muds became trapped as the STM formed. The STM migrates along the middle reach of the estuary over pools of easily erodible, muddy sediment [11,20], similar to "mud reaches" observed in other estuaries including the Hudson [37] and the Weser [50]. During periods of elevated river flow, high suspended-sediment concentrations associated with the STM, in combination with increased salinity stratification, dampened near-bed turbulence and prevented the suspension of faster-settling biologically repackaged material [49]. When river discharge dropped, stratification decreased, allowing fine sediment to disperse through the water column, reducing sediment-induced stratification. This permitted higher bed stresses that could suspend the less erodible, faster settling biologically repackaged pellets in addition to slowly settling muddy flocs, so that lower bulk erodibility and higher bulk settling velocities were observed at the site.

The York River estuary is a cohesive and highly dynamic sedimentary environment subject to the presence of a seasonal STM [20]. Erosion and deposition occur alternately, resulting in temporal and spatial variability in seabed *IJc* profiles [11,47]. Therefore, a cohesive bed model that accounted for time varying *IJc* profiles would best represent the York and could be based on observed seabed profiles.

#### **3. Methods**

A three-dimensional representation of the York River estuary was developed using the Regional Ocean Modeling System (ROMS) v3.1 [51,52]. ROMS solves the hydrostatic Reynolds-averaged Navier-Stokes equations on a curvilinear orthogonal grid with vertical stretched terrain-following coordinates. ROMS provides a choice of boundary conditions, advection schemes, turbulence parameterizations, and sub-models for sediment transport [53]. Our implementation modified a previously described model [54] by adding wind forcing and updating parameters. We used ROMS options for the Mellor and Yamada level 2.5 turbulence closure model [55], a third-order upstream advection method for momentum, and the MPDATA advection method for tracers such as sediment [56]. Bottom drag was parameterized with a logarithmic current velocity profile [57] with a bed roughness of 0.005 cm based on observed bed stresses (discussed below). The ROMS sediment model of [53] was included but modified to account for consolidation and swelling processes following [43] and [33] as described below. The model run simulated two months after first completing month spin-up, and the model performance was evaluated using values obtained from ADV observations. Additionally, we investigated sensitivities of the cohesive bed model to two user-defined parameters.

#### *3.1. Cohesive Bed Model*

The CSTMS cohesive sediment bed model is based on [43] and represents consolidation and swelling processes, and an increase of *IJc* with depth in the sediment, which is typically observed in muddy seabeds. A brief summary of the model implementation in ROMS follows, see [33,58] for more detail. The model represents the seabed at each horizontal grid point as a stack of several bed layers, and stores a value of the critical shear stress for erosion, *IJc* (*m*), at the top of each bed layer, with *m* (kg m<sup>í</sup><sup>2</sup> ) being the sediment mass overlying that layer. Under erosional conditions, surficial sediment is entrained into suspension, and the *IJ<sup>c</sup>* at the bed surface increases as sediment having a higher *IJ<sup>c</sup>* becomes exposed. During times of deposition, the model assumes a low critical stress, *IJc* = 0.01 Pa, for newly deposited material, creating an easily erodible layer at the bed surface.

Relaxation equations account for the temporal effects of consolidation and swelling on *IJc* [43]. The bed sediment is assumed to have an equilibrium critical stress profile *IJceq* (*m*). At the end of every model time-step, the instantaneous critical stress profile, *IJc* (*m*), is nudged toward *IJceq* (*m*) to simulate consolidation and swelling according to:

$$\begin{aligned} \frac{\widehat{\mathcal{C}}}{\widehat{\mathcal{C}}t} \left[ \boldsymbol{\tau}\_{\boldsymbol{c}} \left( \boldsymbol{m} \right) \right] &= \frac{1}{T\_{\boldsymbol{c}}} \Big[ \boldsymbol{\tau}\_{\boldsymbol{c}aq} \left( \boldsymbol{m} \right) - \boldsymbol{\tau}\_{\boldsymbol{c}} \left( \boldsymbol{m} \right) \Big]; \text{ if } \boldsymbol{\tau}\_{\boldsymbol{c}aq} \left( \boldsymbol{m} \right) \geq \boldsymbol{\tau}\_{\boldsymbol{c}} \left( \boldsymbol{m} \right);\\ \frac{\widehat{\mathcal{C}}}{\widehat{\mathcal{C}}t} \Big[ \boldsymbol{\tau}\_{\boldsymbol{c}} \left( \boldsymbol{m} \right) \Big] &= -\frac{1}{T\_{\boldsymbol{s}}} \Big[ \boldsymbol{\tau}\_{\boldsymbol{c}aq} \left( \boldsymbol{m} \right) - \boldsymbol{\tau}\_{\boldsymbol{c}} \left( \boldsymbol{m} \right) \Big]; \text{ if } \boldsymbol{\tau}\_{\boldsymbol{c}aq} \left( \boldsymbol{m} \right) < \boldsymbol{\tau}\_{\boldsymbol{c}} \left( \boldsymbol{m} \right); \end{aligned} \tag{2}$$

where *Tc* and *Ts* are timescales for consolidation and swelling processes, respectively. Reasonable ranges for *Ts* are suggested to be on the order of 100 days, while the *Tc* is expected to be much shorter at around one day [42]. Following [33], profiles for *IJceq* were based on power law fits to erodibility experiments performed by [48] on cores collected in September, 2007, and April, 2007 (Figure 2).

The model calculates erosion when bed shear stress (*IJb*) exceeds *IJc* at the sediment surface (Equation (1)). In its original form, the model allowed *M* to vary with depth as a function of solids volume fraction *ijs* [43]. As implemented here, however, both solids fraction and the erosion rate parameter were held constant with *ijs* = 0.1 and *M* = 1 × 10<sup>í</sup><sup>3</sup> kg/m2 s Pa. Previous analysis concluded that erosion in the York River was primarily depth-limited and hence depended more on the critical stress profile than on the erosion rate parameter [33]. The York River model calculated bed stress assuming a logarithmic current profile in the bottom boundary layer. Stresses were computed by ROMS subroutine *uv\_logdrag* from the velocity in the bottom cell and a user-defined hydraulic roughness parameter (*Z0*) [53,54]. Based on observations *Z0* was defined as a constant 0.005 cmab.

**Figure 2.** Profiles of *IJceq* obtained by power law fit to observations (as in [33]). Power law fit to observed (solid lines) critical stress profiles for April and September, 2007. Triangles and asterisks show observed data from [48].

*3.2. Implementation of Three-Dimensional York River Hydrodynamic Cohesive Bed Model* 

This describes the three-dimensional York River hydrodynamic cohesive bed model, originally developed by [55], and the more recent modifications used in this study. The "Standard Model" refers to the version that best matched observations from the Intermediate site on the York River. Later sections evaluate the sensitivity of calculations to changes in model parameters by comparing results from the Standard Model to those from other implementations.

The horizontal model grid of [41] adopted for this study had an average grid resolution of 170 m in the along-channel direction and 110 m in the cross-channel direction (Figure 3). The major tributaries of the York, the Mattaponi, and Pamunkey Rivers, were represented in the model with only one cell in the cross-channel direction due to their narrow width. The model extended about 60 km up-river from West Point, VA, USA to Hanover, VA, USA on the Mattaponi River, and Beulahville, VA, USA on the Pamunkey River. The vertical grid used stretched terrain-following coordinates, with 20 stretched grid layers, and increased resolution at the seafloor and water surface.

**Figure 3.** The York River estuary cohesive bed ROMS model grid used in in this study had 92 across channel cells and 334 along channel grid cells. Each square in the figure represents twenty-five model grid cells. This paper focuses on model results from the area around the Intermediate site (Clay Bank, VA, USA), which is marked with a black circle.

Velocity and salinity estimates from the final time step of a spin-up model were used to initialize the hydrodynamic model presented here. For this spin-up, ROMS was run for 60 days including a spring-neap tidal cycle with 0.2 m neap amplitude and 0.4 m spring amplitude. The spin-up model used a steady river input equal to the 60-year median freshwater flows of 42 m3 /s and 25 m3 /s from the Pamunkey and Mattaponi Rivers, respectively [59].

Calculations included two cohesive sediment classes having settling velocities of 2.4 mm/s and 0.8 mm/s. These values were based on observed settling velocities of the two dominate particle types in suspension, biologically repackaged pellets and flocculated muds [42,49]. The sediment bed was initialized having a total bed thickness of 1 m, with a uniform distribution of sediment types, containing 20% of the coarser sediment class and 80% of the finer sediment class. Bed *IJ<sup>c</sup>* throughout the grid was initialized as the power law fit used by [33] to represent the September profile, *IJcinit* (*m*) = 1.0 *m*0.62, while the *IJceq* profile was defined by the power law fit to the more erodbile April (more erodible) profile, *IJceq* (*m*) = 0.4 *m*0.55 (Figure 2). Model sensitivity tests, presented in Section 5.2, found that this configuration produced a better match between calculations of suspended sediment concentrations and observed values, compared to other configurations studied. A consolidation timescale of *Tc* = 1 day and a swelling timescale of *Ts* = 25 days were used. The model configuration, specifically the model simulation time, must be considered when defining *IJcinit*, *IJceq*, *Tc* and *Ts*, as discussed further in Section 6.2. In the runs presented here, the model simulation time of only 120 days included 90 days of spin-up time. With this short simulation time the model instantaneous *IJc* profiles were sensitive to *IJcinit* and *Ts*, because of the limited time for the bed to adjust towards *IJceq*.

The model was then run to represent May–August, 2007. Results from May to June were treated as spin-up for the sediment field, and our analysis centered on results from July-August, which coincided with ADV deployment. Freshwater discharges were specified using data from USGS gages 1,674,500 near Beulahville, VA, USA for the Mattaponi River and 1,673,000 near Hanover, VA, USA for the Pamunkey River. Suspended-sediment concentrations at these boundaries were set at 5 mg/L. The value chosen approximates typical suspended-sediment concentrations at the boundaries. Material from upriver sources generally remains trapped at the primary ETM at West Point and therefore has little influence on conditions downstream. Wind velocities were assumed to be spatially uniform, and hourly velocities measured at the York Coast Guard Meteorological station near Gloucester Point (see Figure 1, weather station) were used as model input.

The open boundaries at the mouth of the river where the York River meets Chesapeake Bay were specified as follows to account for tides, sediment fluxes, and the salinity gradient. For sea surface elevation, we used data from about 10 km upstream at the US Coast Guard pier (see Figure 3 for location). To improve the agreement between the modeled and observed time series of water elevation at the Coast Guard pier, the data was lagged by 1 h and scaled by a factor of 1.4 before being applied at the open boundary [54]. For the suspended-sediment open boundary condition, a zero-gradient condition was applied.

 Salinity in the interior of the model was calculated as a state variable, but near the open boundary the salinity gradient was specified as follows using an empirical relationship that assumed the salinity along the river fit the hyperbolic tangent function [57],

$$\overline{S}(\mathbf{x}) = \frac{S\_0}{2} \left[ 1 + \tanh\left(2 - \frac{X}{\beta}\right) \right] \tag{3}$$

where *S0* was the maximum salinity at the river mouth, *X* was the along-channel distance, and *ȕ* was a length scale for the salt intrusion, based on river flow. These were used to estimate the crosschannel averaged salinity, The along-channel gradient at the open boundary could then be specified from the derivative of Equation (3) with respect to *X*, *S x*( ).

$$\frac{\partial}{\partial \mathbf{x}} \mathbf{s}(\mathbf{l}, \boldsymbol{\chi}, z, t) = -\frac{S\_0}{2\beta} \mathbf{s} \text{sech}^2 \left( \tanh^{-1} \left( \frac{\mathfrak{L} \mathbf{s}\_1}{S\_0} - \mathbf{l} \right) \right) \tag{4}$$

where *s1* and were the salinity and the salinity gradient at the first interior grid point. The salinity at the boundary *s0* was then defined as:  1, , , *<sup>s</sup> yzt x* w w

$$
\Delta s\_0 = s\_1 - \frac{\widehat{\text{Cs}}}{\widehat{\text{c}} \mathfrak{x}} \Delta \mathfrak{x} \tag{5}
$$

where *ǻx* was the along channel grid cell width. Data from the Chesapeake Bay Program [60] were used to fit *ȕ* in Equation (3). Best fit values of *ȕ* were then regressed against a four-day running mean, four-day lagged river flow (*Q4*) resulting in These data fit the hyperbolic 4 E290 50. *Q*

tangent function (Equation (3)) with *r*<sup>2</sup> = 0.65. The exercise described above (Equations (3)–(5)) was used to set a flow-dependent open boundary condition for salinity, as has been done for a model of the Hudson River estuary [57].

#### *3.3. Acoustic Doppler Velocimeter (ADV) Observations*

Sontek ADV Ocean Probes have been maintained nearly continually at the MUDBED Intermediate site since 2006 (see Figure 1). The ADVs were attached to bottom-mounted tripods and measured three-dimensional water velocity (*u*,*v*,*w*) and acoustic backscatter at roughly 35 cm above the bed, see [10] for details of the instrumentation. This project used observations collected by a 5MHz ADV deployed from June to August 2007 to estimate burst-averaged bottom stress (*IJb*, via Reynolds averaging of turbulent velocity), suspended-sediment concentration (*C*, via acoustic backscatter calibrated by pump samples), and bulk sediment settling velocity (*WsBULK*, via an assumed Rouse balance between upward Reynolds flux and gravitational settling) [8,10,61].

### *3.4. Standard Model Evaluation*

The Standard Model, run from July–August 2007, encompassed a time period for which ADV data and bed samples from the Intermediate site were available. This study especially focused on model behavior over individual tidal phases at the MUDBED Intermediate site. Model estimates from a three-by-three set of adjacent grid cells were averaged to represent the Intermediate site (Figure 3, black dot). These grid cells were chosen because they were within the area where the ADV tripod was deployed, these model cells had water depths similar to that of the observations (~5–6 m), and previous studies collected sediment cores from this location [48].

Data averaged over tidal phases collected by the ADV during this study period has been analyzed to examine the influence of the presence and departure of the STM on sediment dynamics, using the tidal cycles having the top 20% of the observed bed stresses [49]. Typically, bed stresses in the York estuary are low, (tidal max *IJb* < 1 Pa [62]). Our focus on the most energetic tidal cycles did bias our analysis towards spring-tide periods relative to neap-tide periods, but highlighted the study on times when significant sediment would be suspended from the bed. Tidal cycles were defined using the time series of ADV observed current speed (*U*), such that an individual ebb or individual flood formed a full cycle, with its beginning and end corresponding to times of minimum observed *U*. Observations over the tidal cycles were then interpolated to a common interval in the tidal phase to produce tidal-phase averaged values of *U*, *IJb*, *C* and *WsBULK*. No distinction was made between floods and ebbs, removing any flood-ebb asymmetry.

Model results from the Standard Model were phase averaged to obtain representative values of hydrodynamic parameters (*U*, *IJb*, *C*, and *WsBULK*). Model estimates of these values were interpolated to the height of the ADV measurements (~35 cmab) prior to the tidal-phase averaging, and model performance was assessed based on the degree to which model calculations of these values represented those derived from ADV observations. The observed estimates of *WsBULK* were obtained by assuming a Rouse balance, while *WsBULK* from model results were estimated as the average of the settling velocities of the two size classes (2.4 mm/s and 0.8 mm/s), weighted by their respective suspended-sediment

**44**

concentrations at the height of the ADV sensor. Though they used different methods, both the model-calculated *WsBULK* and ADV-based values provided a bulk settling velocity that accounts for all particle sizes present in suspension [63].

Modeled estimates of bed erodibility were calculated by evaluating the depth of the sediment bed that could be resuspended by bed stresses that exceed some threshold, following the definition presented by [11]. For example, the instantaneous profile of critical bed stress (*IJcinit*) estimated for a location in the model typically had a low value at the sediment-water interface and increased with depth in the bed (*m*). The depth at which the critical stress exceeded the threshold of, e.g., 0.2 Pa, would give a model estimate of the erodibility of the bed at *IJb* = 0.2 Pa. Because the (*IJcinit*) changed in response to erosion, deposition, consolidation, and swelling, the model estimates of erodibility varied in both time and space.

#### **4. Standard Model Results**

System behavior is illustrated during times of high (Figure 4) and low (Figure 5) sediment bed erodibility, showing daily-averaged bed stresses, near-bed suspended-sediment concentrations, erodibility at 0.2 Pa, and near-bed settling velocities throughout the York River. These time periods corresponded to a spring tide (Figure 4) and a neap tide (Figure 5). At both times, the upper and middle estuary had increased erodibility and suspended-sediment concentrations, compared to the lower estuary. Both times had similar bed stresses (Figures 4a and 5a), but when erodibility was higher (Figure 4), more material was suspended, having ~0.9 kg m<sup>í</sup><sup>2</sup> more material in suspension than during the less erodible situation in Figure 5. The increased concentrations were associated with lower near bed settling velocities, particularly in the middle estuary (Figure 4b,d). Ten days later, erodibilty was lower (Figure 5), and decreased suspended-sediment concentrations (~100 mg/L) and higher near-bed settling velocities characterized conditions in the middle estuary (Figure 5b,d). This result was consistent with field observations indicating that during periods of high erodibility, the suspended material was composed of more slowly settling flocculated muds [10,11,48]. The model was also able to capture observed trends in erodibility attributed to the spring/neap tidal cycle. Generally, the sediment bed is more erodible during spring tides and less erodible during neap tides [64].

The results of the tidal-phase average analysis for current speed, concentration, and bed stresses from the model were compared to the ADV estimates from [49]. The model reproduced realistic patterns of all of the parameters over individual tidal cycles. Both model and observed bed stresses and concentration increased with current speed (Figure 6). Peak bed stresses and current speeds in the model had an offset lag compared to the observations (Figure 6a,c). This suggested that the modeled tides may have been slightly out of phase with the observed ones. Overall, the model had more skill at predicting tidal-phase averaged current speeds than the other values considered (Figure 6a). Only at peak flow did model current speeds show a slight bias, underestimating measured speed by less than 5 cm/s (about 12%). The model overestimated suspended-sediment concentrations during the accelerating and decelerating phases of the tide, but matched the concentrations at peak flow (Figure 6b). Bed stresses were overestimated in the model over the

entire tidal phase by as much as 0.05 Pa (~18%) (Figure 6c), but this bias was accepted because the modeled currents and sediment concentrations matched reasonably well.

**Figure 4.** Daily-averaged model estimates of (**a**) bed stress, (**b**) near bed suspended sediment concentration, (**c**) mass eroded at 0.2 Pa, and (**d**) near bed settling velocity during a spring tide when high erodibility was found mid-estuary (15 July 2007).

(**c**) Eroded Mass at 0.2 Pa (erodibility) (**d**) Near Bed Setting Velocity

(**a**) Bed Stress (**b**) Near Bed SSC

The most reliable approximation of *WsBULK* from ADV data was obtained during the accelerating portion of the tidal phase, when the assumptions in the Rouse balance were most likely to be appropriate [49]. During the accelerating phase of the tide, the tidal-phase averaged analysis showed that observed *WsBULK* increased from about 1 mm/s to 1.4 mm/s (Figure 7). More specifically, the tidal-phase averaged analysis showed that observed *WsBULK* increased from about 1 mm/s to 1.4 mm/s during the accelerating tide (Figure 7). The model represented the observed changes in *WsBULK* seen in the phase averaged analysis, with modeled *WsBULK* increasing from about 0.9 mm/s to 1.7 mm/s, though the model estimated slightly higher (~0.3 mm/s) *WsBULK* at peak bed stresses.

**Figure 5.** Daily-averaged model estimates of (**a**) bed stress, (**b**) near bed suspended sediment concentration, (**c**) mass eroded at 0.2 Pa and (**d**) near bed settling velocity during a neap tide when lower erodibility was found mid-estuary (25 July 2007).

**Figure 6.** ADV and model estimated (see legend) (**a**) current speed, (**b**) concentrations, and (**c**) bed stresses for the top 20% of tidal cycles with strongest bed stresses. Error bars denote ±1 standard error.

**Figure 7.** ADV and model estimated (see legend) bulk settling velocity during periods of increasing tidal velocity for the top 20% of tidal cycles with the highest bed stresses. Error bars denote ±1 standard error.

#### **5. Cohesive Bed Sub-Model Sensitivity**

Additional model runs were analyzed to investigate the sensitivity of the cohesive bed model to the swelling time (*Ts*), critical shear stress equilibrium profile (*IJceq*), and initial critical shear stress profile (*IJcinit*). The sensitivity to *Ts* (section 5.1) was evaluated by varying it from 2 to 50 days, bracketing the value of 25 days used in the Standard Model. To test the sensitivity to *IJceq* and *IJcinit* (section 5.2), models were run that defined *IJceq* = *IJcinit* using power-law fits based on the September or the April data (Figure 2). Recall that in the Standard Model, the equilibrium and initial critical shear stress profiles differed as *IJceq* was based on the April power law fit, while *IJcinit* was based on the September profile.

### *5.1. Sensitivity to Ts*

Results from model runs that used different values of *Ts* (specifically 2, 25, and 50 days) were compared to investigate the sensitivity of the instantaneous *IJc* profiles and resultant suspendedsediment concentrations to swelling timescale. The model configurations were identical to the Standard Model except for the swelling timescales to isolate the influence of *Ts*. As for the Standard Model, different critical stress profiles were used to initialize the model, and for the equilibrium state toward which the model calculations were nudged (*IJceq*). Bed *IJc* throughout the grid was initialized as the power law fit used by [33] to represent the September sediment bed profile while the *IJceq* profile was defined by the power-law fit found to represent the April (more erodible) profile (Figure 2). *Ts* = 25 days was used in the Standard Model described above.

Figure 8 shows the values of daily calculated instantaneous *IJc* profiles, along with the user-defined *IJceq* and *IJcinit* for these three models. Recall, the model simulations defined *IJceq IJcrinit*. When a short *Ts* of two days was defined (Figure 8a) the instantaneous *IJc* profiles rapidly became more erodible as the bed adjusted from the less erodible IJcinit to the more erodible *IJceq*. When a longer *Ts* of 50 days was defined (Figure 8c) the instantaneous *IJc* profiles showed minimal adjustment from *IJcini*t to *IJceq*, so that the bed remained more consolidated. In the Standard Model, where *Ts* was defined as 25 days, the instantaneous *IJc* profiles adjusted somewhat from the less erodible *IJcinit* to *IJceq*.

**Figure 8.** Daily instantaneous model calculated IJc profiles (colored lines) shown with user defined equilibrium (*IJceq*) and initial bed profiles (*IJcrinit*) for model runs with swelling times defined as (**a**) 2 days, (**b**) 25 days, and (**c**) 50 days. Note: *IJceq IJcrinit*. Black arrows show the directions that Equation (2) will nudge the instantaneous *IJ<sup>c</sup>* profiles toward *IJceq*.

Figure 9 compares the resulting tidal-phase average concentrations from these three models and the values derived from ADV observations. The concentrations increased as *Ts* decreased. Comparison of the three model simulations suggested that using a *Ts* of 25 days resulted in the most realistic estimates for suspended-sediment concentration, justifying the choice of *Ts* = 25 days in the Standard Model.

**Figure 9.** Phase-averaged concentration from ADV observations (black), and estimated by model using swelling times defined as 2 days (blue), 25 days (pink), and 50 days (blue).

#### *5.2. Sensitivity to IJceq and IJcinit*

Power law fits to Gust microcosm data for the York River were used to define two depth-dependent critical shear stress profiles, *IJc* (*m*), one based on cores collected in September 2007 and the second on cores collected in April 2007 [47] (Figure 2). Compared to the September data, the *IJc* (*m*) profiles from April 2007 represented a time that had a more erodible seafloor. For the model runs described in this section, the same profiles were used for model initialization (*IJcinit*) and to define the equilibrium toward which instantaneous profiles (*IJc*) were nudged (*IJceq* in Equation (2)). Results from two model simulations were compared to evaluate the sensitivity of model calculations to parameterization of the sediment bed's critical shear stress profile. One used the power law fit to the September profile for both *IJcinit* and *IJceq*, while the second used the April profiles for both. Like the Standard Model, both assumed timescales for consolidation and swelling as *Tc* = 1 day and *Ts* = 25 days, but recall that the Standard Model used the April, more erodible, profile as *IJceq* and the September, less erodible, profile for model initialization *IJcinit*. A second set of these sensitivity tests used identical parameterization as the first two, except *Ts* = 50 days. The limits on computational time precluded doing additional tests such as a model using *Ts* = 2 days.

Figure 10 shows the calculated instantaneous *IJ<sup>c</sup>* profiles for the first two sensitivity tests, compared to the user-defined profiles of *IJceq* and *IJcini*t. When *IJceq* and *IJcini*t were defined using the September fit, the instantaneous *IJc* profiles showed slight adjustment back to *IJceq* (Figure 10a). This was not the case when *IJceq* and *IJcinit* were based on the more erodible April seabed (Figure 10b). Instead, the increased erodibility of the initial sediment bed provided a pool of easily resuspended material that when redeposited at slack tide, further increased the erodibility of the bed. With the consolidation timescale of one day, the instantaneous *IJc* could not readjust to the equilibrium profile, *IJceq*, but developed into an environment dominated by easily erodible, newly deposited sediment throughout the model run.

**Figure 10.** Daily instantaneous model calculated *IJc* profiles (colored lines) shown with the profile used to parameterize these models (black lines). The equilibrium profile (*IJceq*) and initial bed profile (*IJcrinit*) were defined based on (**a**) September and (**b**) April profiles. Note: *IJceq* = *IJcrinit* for these model runs. The black arrows show the direction that Equation (2) will nudge the instantaneous *IJc* toward the user-defined *IJceq*.

The increased erodibility of the second case led to a positive feedback between bed erodibility and sediment concentrations. This resulted in much higher estimated suspended-sediment concentrations than either the Standard Model or the case using the September profiles for both the equilibrium and initial values, as seen in the tidal phase averaged values (Figure 11). This result was insensitive to the swelling timescale used, which was not surprising because the modeled bed was always freshly deposited and thus easily erodible. Regardless of the defined *Ts*, concentrations estimated by the model simulation that defined *IJceq* and *IJcinit* based on the April fit exceeded the observed values by about 80–300 mg/L, while the model that used the September profiles provided much better estimates.

### **6. Discussion**

This section first summarizes the model behavior for the York River case, and then discusses issues relating to model implementation and model parameters.
