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

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 OPEN ACCESS J. Mar. Sci. Eng. 2014, 2 414 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.


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][4][5][6][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][9][10][11][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][14][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][20][21][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 (τ b ), and total eroded mass from the seabed when τ b exceeds the critical stress (τ c ) of the sediment surface. The critical shear stress, τ c , represents the stress at which motion or suspension of sediment first occurs [23]. Erosion rates have been estimated using many different formulas [24][25][26][27], but the simple linear Ariathurai-Partheniades erosion formulation has been assumed in many cases: Here, the erosion rate, E (kg m −2 s −1 ), varies linearly with the excess shear stress (τ b − τ c ) (Pa) according to the erosion rate parameter, M (kg m −2 s −1 Pa −1 m −2 ) [11,12,[28][29][30][31][32][33]. Though developed for cohesive sediment, Equation (1) can also be applied to non-cohesive beds. For purely non-cohesive beds, τ c mostly depends on the grain size and density of individual particles, and typically increases with diameter [23]. For cohesive sediment, τ c 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 τ c , even for cohesive sediment [7,[34][35][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][38][39][40]. A model of the York River estuary in Virginia that defined a constant value of τ c [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 τ c 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.

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]. 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 τ c profiles [11,47]. Therefore, a cohesive bed model that accounted for time varying τ c profiles would best represent the York and could be based on observed seabed profiles.

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.

Cohesive Bed Model
The CSTMS cohesive sediment bed model is based on [43] and represents consolidation and swelling processes, and an increase of τ c 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, τ c (m), at the top of each bed layer, with m (kg m −2 ) being the sediment mass overlying that layer. Under erosional conditions, surficial sediment is entrained into suspension, and the τ c at the bed surface increases as sediment having a higher τ c becomes exposed. During times of deposition, the model assumes a low critical stress, τ c = 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 τ c [43]. The bed sediment is assumed to have an equilibrium critical stress profile τ ceq (m). At the end of every model time-step, the instantaneous critical stress profile, τ c (m), is nudged toward τ ceq (m) to simulate consolidation and swelling according to: (2) where T c and T s are timescales for consolidation and swelling processes, respectively. Reasonable ranges for T s are suggested to be on the order of 100 days, while the T c is expected to be much shorter at around one day [42]. Following [33], profiles for τ ceq were based on power law fits to erodibility experiments performed by [48] on cores collected in September, 2007, and April, 2007 ( Figure 2).

Figure 2.
Profiles of τ ceq 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].
The model calculates erosion when bed shear stress (τ b ) exceeds τ c 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 φ s [43]. As implemented here, however, both solids fraction and the erosion rate parameter were held constant with φ s = 0.1 and M = 1 × 10 −3 kg/m 2 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 (Z 0 ) [53,54]. Based on observations Z 0 was defined as a constant 0.005 cmab.

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. 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 τ c throughout the grid was initialized as the power law fit used by [33] to represent the September profile, τ cinit (m) = 1.0 m 0.62 , while the τ ceq profile was defined by the power law fit to the more erodbile April (more erodible) profile, τ ceq (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 T c = 1 day and a swelling timescale of T s = 25 days were used. The model configuration, specifically the model simulation time, must be considered when defining τ cinit , τ ceq , T c and T s , 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 τ c profiles were sensitive to τ cinit and T s , because of the limited time for the bed to adjust towards τ ceq .
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], where S 0 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 cross-channel where s 1 and were the salinity and the salinity gradient at the first interior grid point. The salinity at the boundary s 0 was then defined as: 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 (Q 4 ) resulting in These data fit the hyperbolic tangent function (Equation (3)) with r 2 = 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].

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 (τ b , via Reynolds averaging of turbulent velocity), suspended-sediment concentration (C, via acoustic backscatter calibrated by pump samples), and bulk sediment settling velocity (W sBULK , via an assumed Rouse balance between upward Reynolds flux and gravitational settling) [8,10,61].

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 τ b < 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, τ b , C and W sBULK . 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, τ b , C, and W sBULK ). 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 W sBULK were obtained by assuming a Rouse balance, while W sBULK 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 concentrations at the height of the ADV sensor. Though they used different methods, both the model-calculated W sBULK 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 (τ cinit ) 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 τ b = 0.2 Pa. Because the (τ cinit ) changed in response to erosion, deposition, consolidation, and swelling, the model estimates of erodibility varied in both time and space.

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 −2 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. The most reliable approximation of W sBULK 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 W sBULK increased from about 1 mm/s to 1.4 mm/s (Figure 7). More specifically, the tidal-phase averaged analysis showed that observed W sBULK increased from about 1 mm/s to 1.4 mm/s during the accelerating tide (Figure 7). The model represented the observed changes in W sBULK seen in the phase averaged analysis, with modeled W sBULK increasing from about 0.9 mm/s to 1.7 mm/s, though the model estimated slightly higher (~0.3 mm/s) W sBULK at peak bed stresses.

Cohesive Bed Sub-Model Sensitivity
Additional model runs were analyzed to investigate the sensitivity of the cohesive bed model to the swelling time (T s ), critical shear stress equilibrium profile (τ ceq ), and initial critical shear stress profile (τ cinit ). The sensitivity to T s (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 τ ceq and τ cinit (section 5.2), models were run that defined τ ceq = τ cinit 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 τ ceq was based on the April power law fit, while τ cinit was based on the September profile.

Sensitivity to T s
Results from model runs that used different values of T s (specifically 2, 25, and 50 days) were compared to investigate the sensitivity of the instantaneous τ c profiles and resultant suspended-sediment concentrations to swelling timescale. The model configurations were identical to the Standard Model except for the swelling timescales to isolate the influence of T s . 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 (τ ceq ). Bed τ c throughout the grid was initialized as the power law fit used by [33] to represent the September sediment bed profile while the τ ceq profile was defined by the power-law fit found to represent the April (more erodible) profile ( Figure 2). T s = 25 days was used in the Standard Model described above. Figure 8 shows the values of daily calculated instantaneous τ c profiles, along with the user-defined τ ceq and τ cinit for these three models. Recall, the model simulations defined τ ceq ≠ τ crinit . When a short T s of two days was defined ( Figure 8a) the instantaneous τ c profiles rapidly became more erodible as the bed adjusted from the less erodible τ cinit to the more erodible τ ceq . When a longer T s of 50 days was defined ( Figure 8c) the instantaneous τ c profiles showed minimal adjustment from τ cinit to τ ceq , so that the bed remained more consolidated. In the Standard Model, where T s was defined as 25 days, the instantaneous τ c profiles adjusted somewhat from the less erodible τ cinit to τ ceq .   . 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).

Sensitivity to τ ceq and τ cinit
Power law fits to Gust microcosm data for the York River were used to define two depth-dependent critical shear stress profiles, τ c (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 τ c (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 (τ cinit ) and to define the equilibrium toward which instantaneous profiles (τ c ) were nudged (τ ceq 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 τ cinit and τ ceq , while the second used the April profiles for both. Like the Standard Model, both assumed timescales for consolidation and swelling as T c = 1 day and T s = 25 days, but recall that the Standard Model used the April, more erodible, profile as τ ceq and the September, less erodible, profile for model initialization τ cinit . A second set of these sensitivity tests used identical parameterization as the first two, except T s = 50 days. The limits on computational time precluded doing additional tests such as a model using T s = 2 days. Figure 10 shows the calculated instantaneous τ c profiles for the first two sensitivity tests, compared to the user-defined profiles of τ ceq and τ cinit . When τ ceq and τ cinit were defined using the September fit, the instantaneous τ c profiles showed slight adjustment back to τ ceq (Figure 10a). This was not the case when τ ceq and τ cinit 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 τ c could not readjust to the equilibrium profile, τ ceq , but developed into an environment dominated by easily erodible, newly deposited sediment throughout the model run.  (colored lines) shown with the profile used to parameterize these models (black lines). The equilibrium profile (τ ceq ) and initial bed profile (τ crinit ) were defined based on (a) September and (b) April profiles. Note: τ ceq = τ crinit for these model runs. The black arrows show the direction that Equation (2) will nudge the instantaneous τ c toward the user-defined τ ceq .
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 T s , concentrations estimated by the model simulation that defined τ ceq and τ cinit 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.

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

Summary of York River 3-D Hydrodynamic Cohesive Bed Model Performance
Implementation of the cohesive bed model in conjunction with the MUDBED field study provided a unique opportunity to evaluate model behavior both in terms of the system-wide feedbacks between bed erodibility and sediment concentration, and more specifically in comparison to field data from the MUDBED Intermediate site. Analyses of model results from July-August 2007 showed that the model represented trends in the spatial and temporal variability of bed erodibility and settling velocity observed in the York River estuary (Figures 4 and 5). Figure 11. Phase averaged concentration observed using the ADV, and estimated by model runs that defined (τ ceq ) and (τ crinit ) using the power law fits to the September and April profiles with swelling times of 50 and 25 days (see legend).
Ongoing analysis of ADV data collected at the MUDBED Intermediate site during the summer of 2007 revealed that the variability in settling velocity and bed erodibility appear to fall into two regimes having distinct hydrodynamic and sediment characteristics [49]. Regime 1 represented periods dominated by easily suspended, slow settling flocculated muds, while Regime 2 represented periods strongly influenced by less easily suspended, faster settling, biologically repackaged fecal pellets mixed with flocs [49]. The analysis presented in this paper focused on a time classified as Regime 2.
Results from the three-dimensional York River Hydrodynamic Cohesive Bed model showed good agreement to field data from the observational period defined as Regime 2 [49], with only slight differences between model estimated and ADV observed bed stresses and bulk settling velocities (~18%). Though the model was biased towards higher bed stresses ( Figure 6c) and higher settling velocities at peak flow, both model estimated and observed bulk settling velocities increased along with bed stress (Figure 7). This implied that differential resuspension of sediment types having the range of settling velocities considered in the model (0.8 to 2.4 mm/s) could account for the changes to settling velocity seen during a tidal cycle, without invoking aggregation and disaggregation processes. Looking more closely, the model overestimated bulk settling velocity as flow accelerated, with the biggest difference occurring at peak flow, which can be attributed to two things. First, the model only included two settling classes (0.8 mm/s and 2.4 mm/s), while suspended material in the York River estuary consists of many more particle types having a wide range of settling velocities [42,65]. Second, the model overestimated bed stresses throughout the tidal phase ( Figure 6), which may have contributed to the overestimated concentrations (Figure 6b) as well as bulk settling velocities (Figure 7). Higher bed stresses were able to suspend a greater mass of material, including a greater amount of faster settling particles, which resulted in a higher estimate of bulk settling velocity [42]. Additionally the higher bed stresses would increase eddy diffusivity, mixing the more slowly settling material higher in the water column, and resulting in a greater fraction of faster settling material near the bed.
Though not shown in this paper, the model had less skill when applied to the period defined as Regime 1, overestimating the observed modeled bed stresses by 50%, while underestimating suspended-sediment concentrations by 17%. Field data imply that sediment-induced stratification more strongly influenced the system during Regime 1 than during Regime 2, and that stratification from both salinity and suspended-sediment gradients suppressed bed stresses during Regime 1 [49,61]. ROMS has the capacity to account for sediment induced stratification, but the process was neglected in this implementation. Future work will evaluate whether the turbulence closures in ROMS can represent sediment-induced stratification at the relevant vertical scales for the York River estuary, and whether including this process this improves the model skill during Regime 1.

Cohesive Bed Model Sensitivities
The cohesive bed model sensitivity tests done in this study explored the impact on estimated instantaneous critical stress profiles (τ c ) and suspended-sediment concentration to user-defined cohesive bed swelling time (T s ) and the initial and equilibrium τ c profiles (τ cinit and τ ceq ). We did not examine sensitivity to the consolidation timescale, because it seemed better constrained in the literature than the swelling timescale [33,43].
Use of a short T s allowed the instantaneous τ c profiles to rapidly adjust to the equilibrium profile for a case where the system was less erodible than the equilibrium (Figure 8a). As T s increased to 50 days, the adjustment of the instantaneous τ c profiles became negligible (Figure 8c). Use of a shorter T s resulted in a more erodible bed, which increased suspended-sediment concentrations in the water column ( Figure 9). The sensitivity to T s should be considered in the context of the model configuration, which for the Standard Model used a less erodible bed for initialization than for the equilibrium profile. The model simulation time, 120 days, included 90 days before the time period analyzed. When swelling timescale was long compared to the spin-up time, the model would be especially sensitive to the seabed initialization (i.e., τ cinit ). For example, using a T s of 50 days meant that the instantaneous τ c profiles would be especially sensitive to the initialization over the timescale considered by the model.
For the case where the model's seabed was initialized to be especially erodible, however, positive feedbacks between sediment erodibility and concentrations led to the development of an even more erodible bed, and model estimates were then insensitive to T s ( Figure 10). This sensitivity test found that the model estimated instantaneous τ c profiles and suspended-sediment concentrations were less sensitive to T s than the profiles that were used to define τ cinit and τ ceq (Figures 10 and 11). We expect that the model estimates for this case would be more sensitive to T c than T s , but this was not evaluated directly, in part because the literature provides guidance for the consolidation timescales in the literature than the swelling timescales.
An important step in model implementation for the cohesive bed model was specification of both the initial and equilibrium critical stress profiles, and the consolidation and swelling timescales. The cohesive bed model is still fairly new, and a consensus for how the initial and equilibrium profiles should be defined has not been established. It is clear, however, that for a case where measured critical shear stress profiles exist, the model parameters, including both the initial and equilibrium profiles, should be defined so that instantaneous shear stress profiles usually remain within the realm of measured profiles. Additionally, the choice of erosion rate parameter (M) probably would influence model behavior and have implications for the suitability of different values of these parameters, but the value of M = 1 × 10 −3 kg/(m 2 s P) was not modified for this implementation from the value originally used by [54]. Each sensitivity test required significant computer time (a few days on our cluster), and full three-dimensional model runs were needed to evaluate system behavior under various parameter sets. For that reason, an exhaustive set of parameter choices was not explored, and we were limited to the cases discussed in section 5.2.
The availability of field observations from the MUDBED Intermediate site were valuable for constraining reasonable seabed profiles of critical shear stress, and water column estimates of bulk settling velocity and sediment concentration. In cases when field observations are unavailable, model estimates may reproduce realistic processes and trends, but should be analyzed with caution. Although this study offered much insight on implementing the cohesive bed model in a realistic environment, ongoing work is needed to further understand the use of the cohesive bed model in a complex system.

Conclusions
This study described the use of the CSTSMS cohesive bed model that accounts for erosion, deposition, consolidation and swelling in a three-dimensional hydrodynamic model to represent the York estuary, VA, USA. Results from a time period when sediment-induced stratification seemed to be unimportant showed good agreement with observations of sediment concentration, bed stress, and current velocity. The sensitivities of the cohesive model were explored by evaluating the impact on calculations of changes to the swelling timescale and the initial and equilibrium critical stress profiles, and model estimates were found to be sensitive to these parameters. The Standard Model was developed by choosing a swelling timescale and initial and equilibrium critical stress profiles for the York River three-dimensional model that showed good agreement with field observations of bed stress, current velocity, and sediment concentration using a tidal-phase averaged analysis. The model was able to capture observed spatial and temporal trends in settling velocity, and unlike previous models that did not allow for a depth varying τ c this model was able to adequately represent feedbacks between erosion and deposition in a highly dynamic cohesive environment. The cohesive bed model is still fairly new, and more work is needed to determine how the initial and equilibrium profiles should be defined, especially in cases where observations are limited. However, results from this study demonstrate that 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. by the U.S. Government. This is contribution number 3363 of the Virginia Institute of Marine Science, College of William & Mary.

K.
A. Fall completed model runs and data analysis. She drafted the manuscript, which the co-authors then edited. C.K. Harris supervised the model runs and analysis. C.T. Friedrichs supervised the data analysis. J.P. Rinehimer initially implemented the model for the York River. C.R. Sherwood oversaw changes to the model source code to include the cohesive bed formulation.