Cohesive and mixed sediment in the Regional Ocean Modeling 1 System (ROMS v3.6) implemented in the Coupled Ocean 2 Atmosphere Wave Sediment-Transport Modeling System (COAWST r1179)

. We describe and demonstrate algorithms for treating cohesive and mixed sediment that have been added to the 15 Regional Ocean Modeling System (ROMS version 3.6), as implemented in the Coupled Ocean Atmosphere Wave Sediment- 16 Transport Modeling System (COAWST Subversion repository revision 1179). These include: floc dynamics (aggregation 17 and disaggregation in the water column); changes in floc characteristics in the seabed; erosion and deposition of cohesive 18 and mixed (combination of cohesive and non-cohesive) sediment; and biodiffusive mixing of bed sediment. These routines 19 supplement existing non-cohesive sediment modules, thereby increasing our ability to model fine-grained and mixed- 20 sediment environments. Additionally, we describe changes to the sediment bed-layering scheme that improve the fidelity of 21 the modeled stratigraphic record. Finally, we provide examples of these modules implemented in idealized test cases and a 22 realistic application. components spectral irradiance growth model and a model for treating the effects of aquatic vegetation on waves and currents 2017). The present paper describes new components that model processes associated with cohesive sediment (mud) and mixtures of sand and mud. These include aggregation and disaggregation of flocs in the water column, sediment exchange with a cohesive bed where erosion is limited by a bulk critical shear stress parameter that increases with burial depth, and tracking stratigraphic changes in response to deposition, erosion, and biodiffusive mixing. Our goal is to demonstrate that the algorithms reproduce some of the important behaviors that distinguish cohesive sedimentary environments from sandy ones, and to demonstrate their utility for modeling muddy environments. The model processes are presented and discussed in Section 2. Additional details of the model implementation and their use in ROMS are presented in the Supplement. Examples of model behavior are presented in Section 3, and a realistic application in the York River Estuary is presented in Section 4. Discussion and Conclusions are in Sections 5 and 6. throughout the model calculations. Seafloor properties that describe the condition of the sediment surface are stored with spatial dimensions that correspond to the horizontal model domain. Seafloor properties include representative values (geometric means) of sediment properties in the top layer, including grain size, critical shear stress for erosion, settling velocity, and density; and properties of the sediment surface, such as ripple height, ripple wavelength, and bottom roughness. Seabed properties (i.e. stratigraphy) are tracked at each horizontal location and in each layer in the bed. The number of layers used to represent seabed properties is specified as input and remains constant throughout the model run. The mass of each sediment class, bulk porosity, and average sediment age is stored for each bed layer. The layer thickness, which is calculated from porosity and the mass and sediment density for each class is stored for convenience, as is the depth to the bottom of each layer. Additional information for bulk critical shear stress is stored if the cohesive sediment formulation is being used.


Motivation 28
Fine cohesive sediment (mud) is present in almost every coastal environment, and influences water clarity, benthic habitats, 29 shoaling of harbors and channels, storage and transport of nutrients and contaminants, and morphologic evolution of 30 wetlands, deltas, estuaries, and muddy continental shelves (Winterwerp and van Kesteren, 2004; Edmonds and Slingerland, 31

Previous Modeling Efforts 54
Amoudry and Souza (2011) surveyed regional-scale sediment-transport and morphology models, and found that one of the 55 shortcomings was the treatment of cohesive-and mixed-sediment models. The water-column behavior of cohesive sediment 56 (e.g., flocculation and disaggregation, and settling) and the consolidation of settling particles to form a cohesive bed has been 57 modeled mostly with one-dimensional vertical (1DV) models or with empirical formulae that allow particle settling velocity 58 to vary as a function of salinity (Ralston et al., 2012) or suspended-sediment concentration (e.g., Mehta, 1986;Lick et al., 59 1993;Van Leussen, 1994;Lumborg and Windelin, 2003;Lumborg, 2005;and Lumborg and Pejrup, 2005). The primary 60 dynamical effect of flocculation is to increase settling velocities, thereby increasing the mass settling flux. Soulsby et al. 61 (2013) reviewed methods for estimating floc settling velocities and proposed a new formulation that depends primarily on 62 turbulence shear and instantaneous suspended-sediment concentration. Spearman et al. (2011) noted that adjustments to 63 settling velocity (e.g., Manning and Dyer, 2007) were able to successfully reproduce floc settling in one-dimensional estuary 64 modeling applications. However, these approaches do not allow analysis of other characteristics of the suspended particle 65 field, such as acoustic and optical properties or geochemical properties (water content and surface area). Full floc dynamics 66 have been incorporated in only a few coastal hydrodynamics and sediment-transport models. Winterwerp (2002) 67 incorporated his floc model (Winterwerp, 1999)

Floc Model 143
We implemented the floc model FLOCMOD (Verney et al., 2011) in ROMS to model changes in settling velocity and 144 particle size caused by aggregation and disaggregation. Flocculation is represented as a local process that moves mass among 145 the floc classes within each model grid cell during a ROMS baroclinic time step. Subsequent advection and mixing of floc 146 particles is performed along with other tracers (heat, salt, sand, biogeochemical constituents). FLOCMOD is a population 147 model (Smoluchowski, 1917) based on a finite number of size classes with representative floc diameters Df (m). The model 148 requires a relationship between floc size and floc density f ρ (kg/m 3 ) that is related to the primary disaggregated particle 149 diameter Dp (m) and density s ρ (kg/m 3 ) through a fractal dimension nf (dimensionless; Kranenburg, 1994)  conserving suspended mass in that cell, similar to the way that reaction terms are included in biogeochemical models (for 157 example, Fennel et al., 2006). FLOCMOD simulates aggregation from two-particle collisions caused by either shear or 158 differential settling, and disaggregation caused by turbulence shear and/or collisions. The rate of change in the number 159 concentration N(k) (m -3 ) of particles in the k th floc class is controlled by a coupled set of k of differential equations 160 where G and L terms (m -3 s -1 ) represent gain and loss of mass by the three processes denoted by subscripts: a (aggregation), 162 bs (breakup caused by shear), and bc (breakup caused by collisions). Equations 2 are integrated explicitly using adjustable 163 time steps that may be as long as the baroclinic model time step, but are decreased automatically when necessary to ensure 164 stability and maintain positive particle number concentrations. Particle number concentrations N(k) are related to suspended 165 mass concentrations Cm(k) (kg/m 3 ) via the volume and density of individual flocs. The aggregation and disaggregation terms 166 (Verney et al., 2011) both depend on local rates of turbulence shear, which are calculated from the turbulence submodel in 167 The floc model introduces several parameters (see Supplement), some of which have been evaluated by Verney et al. (2011). 169 These parameters are specified by the user. The equilibrium floc size depends on the ratio of aggregation to breakup 170 parameters, and the rate of floc formation and destruction depends on their magnitudes (Winterwerp, 1999;2002). The 171 diameter, settling velocity, density, critical stress for erosion, and critical stress for deposition (described below) are required 172 inputs for each sediment class, both cohesive and non-cohesive (see Supplement). The present implementation requires a 173 fractal relationship between floc diameter and floc density (Kranenburg, 1994) deposited as flocs can be eroded as denser, more angular aggregates (Stone et al., 2008). However, we find little guidance for 198 constraining this process. We therefore have implemented deflocculation, a simple process that stipulates an equilibrium 199 cohesive size-class distribution and an associated relaxation time scale. The time-varying size-class distribution in the bed 200 tends toward the user-specified equilibrium distribution while conserving mass (see Supplement). This influences the 201 amount of material in classes that are available for resuspension when a cohesive bed is eroded. Example cases presented 202 below demonstrate the effect of this process and the associated time scale on floc distributions both in the bed and in the 203 water column. 204

Stratigraphy 205
Stratigraphy serves two functions in the model as conditions change and sediment is added or removed from the bed: (1) to 206 represent the mixture of sediment available at the sediment-water interface for use in bedload transport, sediment 207 resuspension, and roughness calculations; and (2) to record the depositional history of sediment. Bookkeeping methods for 208 tracking and recording stratigraphy must conserve sediment mass and must accurately record and preserve age, porosity, and 209 other bulk properties that apply to each layer. Ideally, a layer could be produced for each time step in which deposition 210 occurs, and a layer could be removed when cumulative erosion exceeds layer thickness. In practice, the design of many 211 models is subject to computational constraints that limit resolution to a finite and relatively small number of layers. In 212 ROMS, this number is declared at the beginning of the model run and cannot change. Thus, when deposition requires a new 213 layer, or when erosion removes a layer, other layers must be split or merged so that the total number of layers remains 214 unchanged. Where and when this is done determines the fidelity and utility of the modeled stratigraphic record. Some 215 models have used a constant layer thickness (Harris and Wiberg, 2001); others (for example, ECOMSED) define layers as 216 isochrons deposited within a fixed time interval (HydroQual, Inc., 2004). Our approach is most similar to that described by 217 Le Hir et al. (2011) in that we allow mixing of deposited material into the top layer, and require a minimum thickness of 218 newly formed layers, merging the bottom layers when a new layer is formed. Likewise, the bottom layer is split when 219 erosion or thickening of the active layer, discussed below, reduces the number of layers. The sequence of layer calculations 220 is described in detail in the Supplement. 221 A key component of the bed model is the active layer (Hirano, 1971), which is the thin (usually mm-scale), top-most layer of 222 the seabed that participates in exchanges of sediment with the overlying water. During each model time step, deposition and 223 erosion may contribute or remove mass from the active layer. Any stratigraphy in the active layer is lost by instantaneous 224 mixing (Merkel and Klopmann, 2012), but this is consistent with the original concept of Hirano (1971)

Bulk Critical Shear Stress for Erosion for Cohesive Sediment 229
An important difference between cohesive and non-cohesive sediment behavior is that the erodibility of cohesive sediment is 230 treated primarily as a bulk property of the bed, whereas the erodibility of non-cohesive sediment is treated as the property of 231 individual sediment classes. The erodibility of cohesive sediment often decreases with depth in the bed, resulting in depth-232 limited erosion (Type 1 behavior according to Sanford and Maa, 2001). When the cohesive bed module is used, the 233 erodibility of cohesive beds depends on the bulk critical shear stress for erosion cb τ (Pa), which is a property of the bed

237
There is no generally accepted physically based model for determining cb τ from bed properties such as particle size, 238 mineralogy, and porosity. We adopted Sanford's (2008) heuristic approach based on the concept that the bulk critical shear 239 stress profile tends toward an equilibrium profile that depends on depth in the seabed ( equilibrium value, while Ts is the time scale for swelling and is applied when the instantaneous profile is less erodible than 255 the equilibrium value. The consolidation time scale is usually chosen to be much shorter than the one associated with 256 swelling (Sanford, 2008). New sediment deposited to the surface layer is assigned a bulk critical shear stress that may either 257 be (1) held constant at a low value (Rinehimer et al. 2008), or (2) set at the instantaneous bed shear stress of the flow. 258

Mixed Sediment 259
Mixed-sediment processes occur when both cohesive and non-cohesive sediment are present, and are typically sensitive to 260 the proportion of mud. Beds with very low mud content (<3%; Mitchener and Torfs, 1996) behave as non-cohesive 261 sediment: erodibility is determined by particle critical shear stress, which is an intrinsic characteristic of each particle class. 262 Non-cohesive beds may be winnowed and armored by selective erosion of the finer fraction. In contrast, beds with more than 263 content behave according to bulk properties that, in reality, depend on porosity, mineralogy, organic content, age, burial 265 depth, etc., but that, in the model, are characterized by the bulk critical shear stress for erosion. Mixed beds in the model 266 have low to moderate mud content (3% to 30%, subject to user specification) and their critical shear stress in the model is a 267 weighted combination of cohesive and non-cohesive values determined by the cohesive-behavior parameter Pc, which ranges 268 from 0 (non-cohesive) to 1 (cohesive; see Supplement). This approach allows fine material (e.g., clay) to be easily 269 resuspended when Pc is low and only a small fraction of mud is present in an otherwise sandy bed, and it limits the flux to 270 the amount available in the active mixed layer. It also allows non-cohesive silt or fine sand embedded in an otherwise muddy 271 bed to be resuspended during bulk erosion events when Pc is high, and it provides a simple and smooth transition between 272 these behaviors. The thickness of the active mixed layer is calculated as the thicker of the cohesive and non-cohesive 273 estimates. Figure 2 illustrates mixed-bed behavior as the mud (in this case, clay-sized) fraction fc increases for a constant 274 bottom stress of 0.12 Pa. At low fc, Pc is zero (Figure 2a), and clay and silt are easily eroded (high relative flux rates out of 275 the bed; Figure 2c) because the particle critical shear stress for non-cohesive behavior of these fine particles is low ( Figure  276 2b). The relative flux rates in Figure 2b are normalized by the fractional amount of each class and the erosion-rate 277 coefficient; the actual erosional fluxes for clay content would be low at Pc = 0 because of the low clay content in the bed. As 278 fc increases and the bed becomes more cohesive, relative erosion flux rates decline. When fc exceeds a critical value (0.2 in 279 the example shown in Figure 2), the bed is completely cohesive and erosion fluxes are determined by bulk critical shear 280 stress for erosion of cohesive sediment cb τ .

Bed Mixing and Stratigraphy 282
Mixing of bed properties in sediment can be caused by benthic fauna (ingestion, defecation, or motion such as burrowing) or 283 circulation of porewater, and tends to smooth gradients in stratigraphy and move material vertically in sediment. The model 284 (e.g., Boudreau, 1997) assumes that mixing is a one-dimensional vertical diffusive process and neglects non-local and lateral 285 mixing processes: 286 where Cv is the volume concentration of a conservative property (e.g., fractional concentration of sediment classes or 289 porosity), Db is a (bio)diffusion coefficient (m 2 /s) that may vary with depth in the bed (see below), and z (m) is depth in the 290 bed (zero at the sediment-water interface, positive downward). We have discretized equation (5) using the varying bed 291 thicknesses and solve it at each baroclinic time step using an implicit method that is stable and accurate (See Supplement). 292 Biodiffusivity is generally expected to decrease with depth in the sediment (Swift et al., 1994;1996), but is often assumed to 293 be uniform near the sediment-water interface. The typical depth of uniform mixing, based on worldwide estimates using 294 radionuclide profiles from cores, is 9.8±4.5 cm (Boudreau, 1994). Rates of biodiffusion estimated from profiles of excess 295 234 Th on a muddy mid-shelf deposit off Palos Verdes (California, USA) varied from ~2 cm 2 /yr to ~80 cm 2 /yr (Wheatcroft 296 and Martin, 1996;Sherwood et al., 2002) and values from the literature range from 0.01 -100 cm 2 /yr (Boudreau, 1997; 297 Lecroart et al., 2010). The depth-dependent biodiffusion rate profile in the model must be specified for each horizontal grid 298 cell using a generalized shape described in the Supplement. 299 Representation of seabed properties, i.e. the stratigraphy, has been modified slightly from the framework presented in 300 During deposition, the new algorithm prevents the second layer from becoming thicker than a user-specifed value, which 305 results in thinner layers that can record changes in sediment composition inherited from the active layer as materials settle. The following cases demonstrate the cohesive-sediment processes included in ROMS, explore model sensitivity to 312 parameters, and provide candidates for inter-model comparisons. 313

Floc Model 314
Tests using a quasi one-dimensional vertical implementation of ROMS were conducted to verify that the floc model was 315 implemented correctly and to gain some insight into model behavior under typical coastal conditions.

Comparison to equilibrium floc size 336
Simulations were conducted to further evaluate the ROMS implementation of FLOCMOD by comparing modeled 337 equilibrium floc sizes to equilibrium floc sizes predicted by Winterwerp (2006). He argued that, in steady conditions, 338 equilibrium floc sizes are determined by the fractal dimension nf, ratio of aggregation rates and breakup rates, concentration 339 C (kg/m 3 ), and turbulence shear rate G (s -1 ). The equilibrium median floc size D50 (m) is given by 340 where kA and kB are aggregation and breakup coefficients, respectively (Winterwerp, 1998

Non-cohesive bed simulation 421
A non-cohesive bed simulation was used to demonstrate the generation and preservation of sand and silt stratigraphy during 422 a resuspension and settling event ( Figure 5). The model was forced with two stress events ~ 1.5 d apart and lasting 1.5 d and 423 1 d respectively. Four sediment classes, representing particles with nominal diameters of 4, 30, 62.5, and 140 µm, particle 424 critical shear stresses of 0.05, 0.05, 0.1, and 0.1 Pa, and settling velocities of 0.1, 0.6, 2, and 8 mm s -1 were used. Although 425 the diameters of the first two sediment classes corresponded to mud, all sediment classes in this experiment were treated as 426 non-cohesive material. The initial sediment bed contained 41 layers, each 1 mm thick, and each holding equal fractions 427 (25%) of the four sediment classes. New sediment layers were constrained to be no more than 1 mm thick. 428 The first, larger stress event (maximum b τ = 1 Pa; Figure 5b), eroded 1.2 cm of bed, and recruited additional fine sediment 429 from the active layer, which extended 2 cm below the initial sediment surface (Figure 5d). The finer fractions dominated the 430 suspended sediment in the water column, which contained only a small fraction of the coarsest sand (Figure 5a). When the 431 stress subsided, coarser sediment deposited first, while finer material remained suspended, producing graded bedding above 432 the 2-cm limit of initial disturbance (Figure 5d). The second stress pulse eroded the bed down to 1 cm but only resuspended 433 minimal amounts of the 140-µm sand. Deposition resumed after the second pulse subsided and, at the end of the simulation, 434 some mud remained in the water column (Figure 5a), leaving the bed with net erosion of 5 mm (Figure 5d). The finest 435 material (4 µm) remained mostly in suspension after five days. The final thickness of the bottom five layers was smaller than 436 their initial value (1 mm), because, to maintain a constant number of bed layers, the deepest layer was split each time a 437 surface layer was formed during deposition. The two stress pulses affected sediment texture down to 2 cm. Above this level, 438 almost all of the finest class was winnowed, and remained mostly in suspension while the other classes settled to the bed, so 439 that the upper bed layers developed a fining-upward storm layer. The bottom portion of the storm layer (1 -2 cm depth) was 440 a lag layer comprised of the two coarsest classes, both because these resisted erosion and because the sand that did erode 441 settled to the bed quickly when shear stress decreased. 442

Mixed bed simulation 443
This case examined the stratigraphic consequences of cohesive behavior resulting from a single bottom-stress event (Figure  444 6). The model configuration was similar to the previous example. The same sediment classes were used, but the two finest

Biodiffusion simulations 466
We validated the numerical performance of the biodiffusion algorithms using two analytical test cases with a realistic range 467 of parameters. The implicit numerical solution is unconditionally stable and conserves mass to within 10 -8 %, but the 468 accuracy depends on time step, gradients in biodiffusivity, and bed thickness. Typical RMS differences in the fractional 469 amount of sediment in a particular class between the numerical solutions and the analytical solutions ranged from 10 -2 to 10 -470 6 . We found that, for modeled beds 5 m thick, solutions improved as layer thickness decreased from 50 to 5 cm, but beyond 471 that, higher resolution did not substantially improve the solution. Even in the worst case, where the numerical solution was 472 off by 1%, it was much more precise than our estimates of biodiffusivity coefficients. 473 Four cases are presented to demonstrate bed mixing (Figure 7). The first two used the same configuration as in the non-474 cohesive (Figures 5, 7a) and mixed-bed simulations (Figures 6d, 7b). The second two were identical to the mixed-bed case 475 except that biodiffusive mixing was enabled. The biodiffusivity profile used was similar to that proposed for the mid-shelf 476 deposit offshore of Palos Verdes, CA (Sherwood et al., 2002) that had a constant diffusivity Dbs from the sediment-water 477 interface down to 2 mm, an exponential decrease between 2 mm and 8 mm, and a linear decrease to zero at 1 cm depth. 478 jump in the fractional distribution between mostly sandy layers and predominantly muddy layers produced in cases that 492 neglected mixing (Figure 7a,b) and the smooth gradient produced by the strong mixing case (Figure 7c). 493

Estuarine Turbidity Maxima 494
High concentrations of suspended sediment often occur near the salt front in estuaries, forming estuary turbidity maxima 495 (ETM). We present ETM test cases that simulated sediment transport in a two-dimensional (longitudinal and vertical) salt-496 wedge estuary with tidal and riverine forcing. The cases investigated the formation of cohesive deposits beneath the ETM 497 with and without floc dynamics. The first case, without floc dynamics but with a mixed bed, is presented here. The second 498 case, presented below, adds floc dynamics. The model was forced with a 12-hour tidal oscillation modulated with a 14-day 499 spring-neap cycle. The idealized estuary was 100-km long with a sloping bottom 4 m deep at the head of the estuary and 10 500 m deep at the mouth (Figure 8a). In all cases, the simulations were run for twenty tidal cycles. Two non-cohesive sediment 501 classes (180-and 250-µm diameter) were represented with equal initial bed fractions (50% of each). One cohesive fraction 502 (37 µm, ρf = 1200 kg/m 3 , ws = 0.13 mm/s) was included, with an initial uniform suspended-sediment concentration of 1 503 kg/m 3 . The bed was initialized without any cohesive sediment, so it initially behaved non-cohesively. Later in the simulation, 504 bed behavior became mixed as suspended mud settled and was incorporated into the initially sandy bed. The chosen 505 equilibrium bulk critical shear stress profile (Equation 3) had slope = 5 ln(kg/m2) and offset = 2 ln(kg/m 2 ), with a minimum 506 value of 0.05 Pa and a maximum of 2.2 Pa. The time scale for consolidation was set to Tc=8 hours (Sanford, 2008; During the simulations, salinity and suspended-sediment field evolved into dynamic equilibria that were repeated over 509 consecutive tides. An estuarine turbidity maximum (ETM) developed between 10 km and 60 km from the mouth of the 510 estuary (Figure 8a) in the salt wedge generated by gravitational circulation and tidal straining (Burchard and Baumert, 1998;511 MacCready and Geyer, 2001). Elevated suspended-sediment concentrations ranging from 0.7 to 2.05 kg/m 3 occupied most of 512 the bottom layer and extended to mid-depth. All of the suspended material was in the 37-µm class (Figure 8a). 513 The second case was identical, except that it included floc dynamics. in the bottom layer (bottom 5% of the water column). The second layer (5 -10% of the water column) had concentrations 520 about half of the bottom layer. The bed sediment response for the two cases also differed. In the no-floc case, the ETM 521 deposit was slightly thinner, located closer to the mouth, and varied less from slack to flood (Figure 8c). Floc dynamics 522 created large tidal variations in the size of bed material (Figure 8d), which ranged up to 600 µm as flocs deposited during 523 slack, and decreased to 37 µm as flocs were resuspended during flood. The behavior in the unflocculated case was less 524 intuitive. Over the course of the simulation, enough fine material accumulated beneath the ETM to cause the bed to behave 525 cohesively, but the top, active layer remained mostly non-cohesive. During flood tide, bottom stresses were sufficient to 526 resuspend the non-cohesive 70 µm material, leaving the cohesive 37 µm material on the bed. Thus, in both cases, the bed 527 became finer during period of higher stress, but for different reasons. The two cases highlight the model-dependent changes 528 in location (driven primarily by settling velocities) and size distributions (driven by floc dynamics) of the ETM. 529 We next expanded the numerical experiment, using six floc cases to elucidate the effects of floc dynamics in the idealized 530 estuary ( Table 1). The two-dimensional model domain was the same as the ETM case described above. Three types of floc 531 behavior in the seabed were investigated: (1) no changes in size distribution occurred in the bed; (2) the bed deflocculation 532 process was invoked, which nudged all cohesive sediment into the 20-µm class over a long time scale (50 hours); and (3) the 533 bed deflocculation process was invoked with a short time scale (5 hours). Additionally, three other combinations of 534 aggregation ( α ) and disaggregation ( β ) rates were used with the slow deflocculation rate to explore floc processes in the 535 water column (Table 1). The following six metrics were compared at the location of the maximum depth-mean suspended-536 sediment concentration (SSC): depth-mean SSC; maximum SSC; median size (D50); 12-h mean of the D50; depth-mean 537 settling velocity ws; and depth-mean ws averaged over a 12-h tidal period (Table 1). The median size and mean settling 538 velocities were weighted by the mass in each class. Also listed in Table 1 are the locus of the maximum deposition, the 539 thickness at that location, and the median size of deposited material at that location. 540 Mean SSC in the ETM did not vary significantly among the floc cases, but the maximum SSC (located lower in the water 541 column) increased when the ratio of aggregation rate / disaggregation rate / α β was higher, which led to larger, faster-542 settling flocs. Among the four cases (3 -6) with slow deflocculation rates in the bed, settling velocities, maximum SSC, and 543 floc size covaried. The locus of maximum deposition of ETM material was insensitive to the deflocculation algorithms 544 (cases 1 -3), and most sensitive to the overall floc rates. The range of ETM locations is listed in Table 6 to highlight the 545 cases where ETM location varied. The case with lowest floc rates (case 5) produced the farthest upriver deposit, with the 546 most variation in the location of the maximum. The case with the highest settling velocities (case 6) produced deposits 547 closest to the estuary mouth. Overall, the simulated ETM was more sensitive to changes in floc parameters than to prescribed 548 behavior of the floc population in the seabed (deflocculation), and the greatest effect of varying floc dynamics was the 549 vertical location of the ETM, which was controlled by floc size and settling velocity. River channels consists of fine-grained material. We found that it was important to modify the sediment bed layering 561 management scheme, as discussed in section 5 below, to resolve the high gradients in bed erodibility evident in the sediment 562 bed model (i.e. Fall et all 2014) and data (i.e. Dickhudt et al. 2009Dickhudt et al. , 2011. 563 In this implementation, sediment deposited to the bed provided an easily erodible layer with an assumed low critical stress, τc 564 = 0.05 Pa. The modeled sediment bed erodibility and suspended-sediment concentrations both were found to be sensitive to 565 parameterization of the equilibrium critical stress profile, and to the consolidation and swelling timescales used (Fall et  The initially uniform bed evolved during the 60-day model run, developing areas of high sediment erodibility along the 582 shoals of the estuary and channel flanks (Figure 10a). In general, sediment was removed from the main channel, which 583 developed reduced erodibility (Figure 10a). At the Gloucester Point site, the initial bed evolved to become less erodible, with 584 a critical shear stress at the seabed that exceeded the equilibrium values specified for the model (Figure 10a). Conversely, at 585 the Clay Bank field site, conditions were variable in space. Sediment deposited on the shoal area, which evolved to enhanced 586 erodibility (Figure 10a). Within the channel, however, the equilibrium critical stress for erosion was often exceeded, 587 resulting in a strongly eroded sediment bed having larger values of critical shear at the sediment surface ( Figure 10a). 588 Resuspension and transport also changed the spatial distribution of sediment classes, with the erosional areas retaining only 589 the coarser, faster-settling classes, while depositional areas retained finer-grained, slower-settling particles (Figure 10 b, c). 590 These patterns, with coarse lag layers and reduced erodibility in the channels relative to the shoals, are consistent with the 591 known grain size distributions and properties of the York River Estuary. 592 5 Discussion 593 version of ROMS, which provides a framework for realistic two-way nested models with forcing from meteorology and 596 waves. ROMS includes options for several turbulence sub-models (e.g., k ε − , k ω − , Mellor-Yamada) and wave-current 597 bottom-boundary layer sub-models that allow us to calculate fields of shear velocity G. Implementation of FLOCMOD in 598 this framework provides a platform for numerical experiments and real-world applications of a full-featured floc model. 599 The primary role of the floc model is to simulate the dynamical response of particle settling velocities to spatial and temporal 600 variations in shear and suspended-sediment concentrations. This can also be achieved with simpler and computationally 601 more efficient parameterization in many applications. What are the advantages of the complex and much slower model 602 implemented here? There are several. The floc model provides fields of particles with dynamically varying density and 603 number of primary particles, which allow calculation of the acoustic and optical responses of the particle fields. In turn, this 604 allows direct comparison with field measurements of light attenuation, optical backscatterance, and acoustic backscatterance, 605 the de facto proxies for suspended-sediment concentration. This also allows calculation of derived properties in the water 606 column, including light penetration and diver visibility. Finally, the modeled particle properties can be used in geochemical 607 calculations that require estimates of particle radius, porosity, and reactive surface area. Depending on the application, this 608 additional information may justify the computational expense of the floc model. 609 The cohesive bed model provides a heuristic but demonstrably useful tool for representing muddy and mixed beds. The 610 cohesive bed framework captures the most important aspects of muddy environment: limitations on erosion caused by 611 increased bed strength with depth in the sediment, and changes toward user-defined equilibrium conditions as deposited (or 612 eroded) beds age. The physical processes of self-compaction and associated changes in porosity and bed strength are not 613 modeled, but the framework of particle-class and bed-layer variables are designed to accommodate a compaction algorithm. 614 The equilibrium profile method implemented here adds little computational expense, but allows the model to represent 615 depth-limited erosion, a key property of many cohesive beds. 616 Modeling stratigraphy effectively is challenging. Although conserving sediment mass among a fixed number of layers is 617 straightforward, it has proven difficult to devise a robust and efficient method that records relevant stratigraphic events in a 618 modeled sediment bed over the wide range of conditions that occur in coastal domains. For both sediment transport and 619 sediment bed geochemistry (i.e. Moriarty et al. 2017), it can be important for the sediment bed model to achieve its highest 620 vertical resolution near the sediment -water interface, but the original ROMS sediment bed model did not meet that goal 621 when the sediment bed was subject to frequent or repeated cycles of erosion. The modifications we have made to the bed-622 layer management have improved the fidelity with which we can record stratigraphic events in the model layers, particularly 623 at the sediment -water interface. Inclusion of biodiffusive mixing is important for environments where biological activity is 624 rapid, compared with sedimentation or physical reworking. Additionally, for problems of sediment geochemistry, it is 625 important to account for mixing of both particulate matter and porewater. Expansion of the ROMS sediment bed model to 626 (2008) has been added, allowing the erodibility of the sediment bed to evolve in response to the erosional and depositional 638 history. Mixing between bed layers has been implemented as biodiffusion using a user-specified diffusion coefficient profile. 639 In addition, the sediment bed layering routine has been modified so that bed layers maintain a high resolution near the 640 sediment water interface, as demonstrated by both our idealized and realistic case studies presented here. The paper presents 641 results of model runs that test and demonstrate these new features and to show their application to real-world systems. The 642 authors encourage the coastal modeling community to use, evaluate, and improve upon the new routines. 643 Code and Data Availability 644 The algorithms described here have been implemented in ROMS (version 3.6) distributed with the Coupled Ocean 645 Atmosphere Waves Sediment-Transport Modeling System (COAWST, Subversion repository revision number 1179). 646 COAWST is an open-source community modeling system with a Subversion source-control system maintained by John C.            Figure 10.