*6.1. 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 (*IJceq*) and (*IJcrinit*) 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.

#### *6.2. Cohesive Bed Model Sensitivities*

The cohesive bed model sensitivity tests done in this study explored the impact on estimated instantaneous critical stress profiles (*IJc*) and suspended-sediment concentration to user-defined cohesive bed swelling time (*Ts*) and the initial and equilibrium *IJc* profiles (*IJcinit* and *IJceq*). 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 *Ts* allowed the instantaneous *IJc* profiles to rapidly adjust to the equilibrium profile for a case where the system was less erodible than the equilibrium (Figure 8a). As *Ts* increased to 50 days, the adjustment of the instantaneous *IJc* profiles became negligible (Figure 8c). Use of a shorter *Ts* resulted in a more erodible bed, which increased suspended-sediment concentrations in the water column (Figure 9). The sensitivity to *Ts* 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.*, *IJcinit*). For example, using a *Ts* of 50 days meant that the instantaneous IJc 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 *Ts* (Figure 10). This sensitivity test found that the model estimated instantaneous *IJc* profiles and suspended-sediment concentrations were less sensitive to *Ts* than the profiles that were used to define *IJcinit* and *IJceq* (Figures 10 and 11). We expect that the model estimates for this case would be more sensitive to *Tc* than *Ts*, 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<sup>í</sup><sup>3</sup> kg/(m2 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.

#### **7. 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 *IJ<sup>c</sup>* 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.

#### **Acknowledgments**

Funding by the National Science Foundation (OCE-1061781 and OCE-0536572) supported Fall, Harris, Friedrichs, and Rinehimer. Lindsey Kraatz (formerly VIMS) contributed to some of the model runs. The authors appreciate discussions with Linda Schaffner, Marjorie Friedrichs (both VIMS) and Larry Sanford (U. Maryland Center for Environmental Science), who aided model development. Comments from Patrick Dickhudt (U.S. Geological Survey) improved an earlier version of this paper. Use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government. This is contribution number 3363 of the Virginia Institute of Marine Science, College of William & Mary.

#### **Author Contributions**

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.

#### **Conflicts of Interest**

The authors declare no conflicts of interest.
