*2.2. Workflow Components*

Sections 2.2.1–2.2.6 describe individual workflow components, each developed to quantify a different component of the sediment dispersal pathway, from delivery of sediment to the norhtern Gulf of Mexico from river discharge, to turbidity current transport in deep water.

#### 2.2.1. River Discharge Modeling Results and Observations

Few rivers that discharge into the Gulf are adequately gauged, with only the Mississippi and Pearl Rivers having associated sediment flux determinations. Therefore, a global *WBMsed* [31,32] was used to estimate daily discharge and sediment flux from rivers into the northern Gulf of Mexico. *WBMsed* combined the Water Balance Model (*WBM*) with the *BQART* and *Psi* models. Specifically, *BQART* simulates long-term (30+ years) average suspended sediment loads for a basin outlet and is based on individual upstream basin properties for each distributed pixel, including geographical, geological and human factors [33]. The *Psi* variability model resolves the suspended sediment flux on a daily time step from the long-term sediment flux estimated by *BQART*, able to capture the intra-annual and inter-annual variability observed in natural river systems [34]. Skill assessment of *WBMsed* was based on a comparison to daily USGS observations of water and sediment discharge [35]; and daily discharge predictions compared favorably to both ground-based gauging stations and satellite-based observations [31,36]. Sixteen rivers that discharge to the northern Gulf of Mexico (Figure 3A) were simulated using observed conditions for 1995–2011.

Freshwater discharge to the northern Gulf of Mexico peaks seasonally in February through March (Figure 3B). The combined Mississippi–Atchafalaya Rivers supply 81% of the freshwater discharge into the northern Gulf; other significant riverine sources are identified in Figure 3A. The combined Mississippi and Atchafalaya flow averages 20,874 m<sup>3</sup>/s with a standard deviation of 11,211 m<sup>3</sup>/s (USACE observations). The *WBMsed* estimate of the combined Atchafalaya–Mississippi discharge for the same period was 18,300 m<sup>3</sup>/s, with a standard deviation of 12,400 m<sup>3</sup>/s. The total predicted (1995–2011) discharge for all northern Gulf rivers was 22,800 m<sup>3</sup>/s, with a standard deviation of 15,400 m<sup>3</sup>/s. The merged discharge of non-Mississippi rivers was 4540 m<sup>3</sup>/s or 19% of the total flow into the northern Gulf of Mexico, with a standard deviation of 4730 m<sup>3</sup>/s (Figure 3C). The 17 y one-day high of these non-Mississippi sources was 30,000 m<sup>3</sup>/s (8 March 1998), highlighting a potential issue of studies that solely consider the Mississippi River discharge. On 13 November 1997, these non-Mississippi sources accounted for 66% of the total freshwater input into the Gulf (of 15,500 m<sup>3</sup>/s), a 17 y one-day maximum contribution.

The USGS observations (1995–2011) of Mississippi River sediment load indicate a mean value of 170 Mt/y. Sediment discharge to the northern Gulf of Mexico is seasonal but with peak loads of short duration (Figure 3D). On average, the Mississippi River supplies 88% of the fluvial sediment load to the northern Gulf, although its contribution varies from more than 99% to less than 15% on any given day. Based on *WBMsed*, the 14 non-Mississippi Rivers identified in Figure 3A supplied 16.2 Mt/y of sediment to the northern Gulf of Mexico during the same period.

**Figure 3.** (**A**) *WBMsed* model: color shows discharge rates for March 1, 2005; numbers identify river outlets: (1) Apalachee, (2) Apalachicola, (3) Conecuh, (4) Choctawhatchee, (5) Escambia, (6) Mobile, (7) Pascagoula, (8) Mississippi and Atchafalaya, (9) Pearl, (11) Grand, (12) Sulphur, (13) Sabine, (14) Neches, (15) Trinity, and (16) Conroe. (**B**) Total daily water discharge (m<sup>3</sup>/s) entering the northern Gulf of Mexico (*WBMsed* simulation). The 2011 flood season (partly shown) was amongs<sup>t</sup> the most devastating floods in the continental US history. (**C**) The percentage of non-Mississippi freshwater discharge entering the northern Gulf of Mexico study region (*WBMsed* simulations). (**D**) Total suspended load (Kt/d) entering the northern Gulf of Mexico (*WBMsed* simulations).

#### 2.2.2. Ocean Hydrodynamic and Wave Model

*ROMS* is a mature numerical framework that represents ocean dynamics over a wide range of spatial (coastal to basin) and temporal (days to inter-annual) scales. A three-dimensional, free-surface, terrain-following ocean model, *ROMS* resolves the primitive momentum and continuity equations modeling large-scale ocean circulation using the hydrostatic vertical momentum balance and Boussinesq approximation [37,38]. The dynamical kernel includes accurate and e fficient algorithms for time-stepping, advection, pressure gradient [38,39], subgrid-scale parameterizations to represent small-scale turbulent processes [40,41] and various bottom boundary layer formulations to determine the stress exerted on the flow by the sediment bed.

For our implementation, two nested *ROMS* grids were run. The larger-scale coarser model represented hydrodynamics over the entire Gulf of Mexico (Figure 4, black box; hereafter, Grid-g) and provided boundary conditions to a finer-scale model that calculated higher-resolution hydrodynamics and sediment transport (Figure 4, red box; hereafter, Grid-f). The initial and lateral boundary conditions for Grid-g were derived from the northwestern Atlantic *ROMS* 50-year solution (courtesy E. Curchitser, Rutgers U.) and the Simple Oceanic Data Assimilation (*SODA*) global reanalysis 50-year dataset (stored as 5-day averages). The annual and monthly temperature and salinity climatology for Grid-g were objectively analyzed from the 1998 World Ocean Atlas. The tidal amplitude and currents (S2, M2, K1, O1 semidiurnal and diurnal components) forcing for Grid-g's open boundaries were derived from the Oregon State University Tidal Prediction Software (*OTPS*). Both Grid-g and Grid-f included river runo ff, and their forcing atmospheric fields were obtained from the European Centre for Medium-Range Weather Forecasts (*ECMWF*) *ERA*-Interim, 3-hour dataset available since 1 January 1978, to the present.

**Figure 4.** Gulf of Mexico bathymetry showing grid domains from the full coarse grid (Grid-g, black box), and the northern nested grid (Grid-f, red box). Black dashed lines mark approximate tracks of Hurricanes Gustav and Ike. Locations of De Soto and Mississippi Canyons also noted.

The bathymetric grids used by *ROMS* were melded between *ETOPO2* (2 arcminute resolution) and 15arcsecond resolution for the shelf and canyons. Apparent in the bathymetry is the narrowing of the continental shelf near the bird-foot delta, and the presence of both the Mississippi and De Soto Canyons. *ROMS* has terrain-following vertical coordinates, preferred for modeling suspended sediment transport. The bathymetry was smoothed to suppress computational errors in the discretization of horizontal operators (pressure gradient, advection, and di ffusion) using a method [42] that allows constraints in the smoothing minimization like preserving the bathymetry in specific grid cells (e.g., on the continental

y

 g

plateau), maximal amplitude modification, desired slope and steepness (**r**-factor), land/sea masking, and preservation of volume.

Spatial and temporal wave data are required to parameterize bottom stress due to wave-current interactions, which affects seafloor sediment transport. *ROMS* requires several wind-induced wavefields to compute bottom stresses from the various bottom boundary layer sub-models available [43]. These fields include significant wave height, wave direction, surface wave period, bottom wave period, bottom orbital velocity, and wave energy dissipation rate. The required fields were processed from the NOAA/NCEP WaveWatch III®dataset (*WW3*; [44]). They were available at three-hour intervals on a grid having a ten arc-minute resolution, and driven by *GFDL-GFS* winds. The *WW3* data were processed from 1 January 2006, to 31 December 2012. Figure 5A,B show wave height and the period during Hurricane Gustav, which impacted the study area during the late summer, 2008. Near-bed wave orbital velocity and near-bed wave periods were estimated from the surface wave characteristics (calculated followng [45]). Figure 5C,D show a sample of bottom wave period and bottom orbital velocity during Hurricane Gustav.

**Figure 5.** Top row: estimates of wave properties for September 1, 2008 during Hurricane Gustav. (**A**) wave height (m), and (**B**) wave period (s) from NOAA-NCEP *Wavewatch III*® model. Black box indicates the location of *ROMS* hydrodynamic Grid-f, and contours show Grid-f bathymetry (m). Bottom row: calculated estimates of (**C**) bottom wave period (s) and (**D**) bottom orbital velocity (m/s) (calculated following [45]). Transects show locations of flux calculations shown in Figure 9A (MW: Mississippi West); Figure 9B–E (MC: Mississippi Canyon), and Figure 9F–I (DC: De Soto Canyon).

#### 2.2.3. Spatial Seabed Datasets

The *dbSEABED* facility [46–48] supplied information on the spatial distributions of seabed sediment type based on interpolations of more than 10<sup>5</sup> individual data records gleaned from numerous published and unpublished sources. The database provides 0.01-degree resolution mappings of mean grain size (Figure 6A), as well as sorting (*Phi*), gravel, sand and mud fractions (%), exposure of rock (%), and sediment carbonate percent. Local patchiness results from the presence of deep cold-water coral banks [49], low-stand shelf-edge delta remnants [50], and methanogenic carbonate rock and

rubble [51]. The shelf areas have important occurrences of gravel, shell, and hard grounds colonized by skeletal-benthos [52].

Sediment stabilities and sediment dispersal patterns are strongly determined by the seabed geomorphology, especially the slope and curvature. To assist the modeling of the generation and then the fate of the turbidity currents, derivative layers were computed from the *SRTM30*+ bathymetry, including slope gradients, and location and dimensions of the channelizations. Channel locations are well-discriminated using integration methods such as contributing area; here, the *PyDEM* package [53] was used for such procedure (see background on Figure 6B). On those features, channel-floor dimensions like widths and gradients were computed using operations on the original gridded bathymetry [54], and their values were mainly supplied to the *RANS*/*TURBINS* turbidity current models.

**Figure 6.** (**A**) The background shows the region's bottom-sediment grain sizes, blue-yellow-red for the range 10 to -8 phi (clay-sand-cobble) from *dbSEABED*. The superimposed points show locations of modeled turbidity current ignitions for Hurricanes Gustav and Ike. Purple points near Mississippi Canyon and Delta mark locations that are particularly prone to ignitions based on density-stability and energy-balance (Knapp–Bagnold, 'KB') measures. The blue points show locations that may also be prone to ignitions. This analysis makes no determination on whether flows will persist over distances. (**B**) The background shows channelized structures (integrated contributing area) for the region, with dark traces marking the most pronounced channels. (Note: The shelf area indications of channeling are dominated by noise in this mapping.) The purple to blue points mark locations of modeled surficial mass failures for Hurricanes Gustav and Ike. The Factors of Safety (*FoS*) indicate the potential for wave-induced mass failure and hence possibly, turbidity currents: purple—very high, blue—high, pale blue—significant. Bathymetric contours shown with depth (m) labeled.

#### 2.2.4. Suspended Sediment Transport Model (*CSTMS*)

The Community Sediment Transport Modeling System (*CSTMS*) has been coupled to the *ROMS* hydrodynamic kernel to represent suspended and bed sediment using user-defined sediment classes; to date, most published *CSTMS* simulations use the non-cohesive routine (see [43]). Each sediment class has attributes of grain diameter, density, settling velocity, and an erosion rate parameter. These are specified in an input file and held constant for the model run. The erodibility of non-cohesive sediment depends on the critical shear stress for erosion (<sup>τ</sup>*cr*), specified for each sediment type in an input file. Suspended transport is estimated by assuming that each sediment class acts independently of the others, and travels along with the ambient current velocities, with the addition of the sediment class' settling velocity. The contribution of suspended sediment to water column density is included in the equation of state, and allows for gravitationally driven bottom-boundary layer flows [55,56]. Net exchange of sediment between the seabed and suspended load are estimated by assuming simultaneous erosion, and deposition via settling [43].

We implemented suspended sediment transport simulations for the northern Gulf of Mexico using *CSTMS* on the three-dimensional Grid-f (see Figure 4) from 1 October 2007, through 30 September 2008. This included periods of energetic waves, elevated fluvial discharge, and also Hurricanes Gustav and Ike (Figure 7). Transport and deposition was calculated for seven sediment classes, representing fluvial and seabed sources (Table 1). Fast- and slow-settling sediment was simulated for riverine sediment from the Mississippi, Atchafalaya, and Mobile Rivers. Discharge and sediment concentrations were derived from USGS gauges. Model calculations included estimates of suspended sediment concentration and flux at each location in the three-dimensional *ROMS* Grid-f, and sediment deposition and erosion at each of the horizontal grid points. While the model timestep was 20 s, the output data was saved at hourly intervals.

**Figure 7.** Observed Mississippi River discharge (USGS), and wave height (NOAA's NDBC Buoy #42889) during the modeled period. Hurricanes were in the Gulf of Mexico between 30 August–1 September (Gustav) and 10–13 September (Ike), 2008.

**Table 1.** Parameters for the suspended sediment transport model. Three sediment classes represented the initial seabed, two sediment classes were discharged by the Mississippi River, and two sediment classes were discharged by the Atchafalaya and Mobile rivers. Critical shear stress and settling velocity for these were based on previous studies [19,20].


#### 2.2.5. Turbidity Current Ignition Models

A package of one-dimensional, time-dependent, process-numerical modeling modules was used to investigate conditions for wave-induced sediment resuspension and mass wasting, which could potentially lead to turbidity current ignitions. Turbidity currents are known to be generated during events of intense sediment resuspension and mass failure, especially over sloping seafloor [57,58].

The inputs to the modeling package *HurriSlip* included a three-hourly spatially-gridded wave climate based on WaveWatch III ®data, surficial seabed material properties from *dbSEABED*, and slope calculations derived from the *SRTM30*+ bathymetry. Whereas the *CSTMS* suspension model uses the significant orbital velocity to calculate bed stresses, the implementation of *HurriSlip* relied on a more energetic member of the wave spectra ( *H1*/*10*) to represent resuspension by extreme waves. The predicted sediment failure and ignition events were passed to the *RANS*/*TURBINS* model-suite, which could simulate the subsequent turbidity current flows down the continental slope. For the predicted cases, the starting flow height, suspended sediment concentration and grain size, and flow velocity were provided. The focus of the work with *HurriSlip* was on the scale of 0-20 m above the seabed with a horizontal resolution of about 1 km.

Sediment resuspension sources: The primary sub-module, *SuspendiSlip*, tested for a likely distribution of turbidity flows arising from wave-induced resuspension of surficial bottom sediment. Most sediment suspension in the continental shelf is thought to be from wave activity during storms [19]. Under significant wave action, bottom-water layers hold significant suspended sediment and turbulent kinetic energy. The module computed several criteria about the ignition of flows. That is, the transformation from bottom waters having significant sediment loading and density to self-sustaining, downslope density-flows undergoing an auto suspension process [59], which allows them to travel for long distances at high speeds. The sediment-laden bottom-water layers were tested from a reference height corresponding to the wave boundary layer thickness up to a height of significant suspension in a Rouse profile. Sediment pickup was modeled using the excess-over-critical bed shear stress for the sediment, using di fferent formulations for muds [60] and sands [61]. Those published formulations focus on granular erosion at low velocities (mostly <0.5 m/s). However, fine sediments under extreme bed shear during storms are known to erode by bulk-failure [62,63]. The fine sediment erosion rates were capped at the values reported in the publications for the highest bed shear stresses to allow for this.

The bulk, densiometric, Richardson Number (*Ri*, non-dimensional) divides layers between subcritical (>1.0) and supercritical (<1.0) on the value of *Ri* = (*gRCh*)/*U*<sup>2</sup> , which depends on gravitational acceleration (*g*, m/s2), sediment grain immersed specific gravity ( *R*, non-dimensional), suspended sediment concentration ( *C*, ppm v/v), flow thickness ( *H*, m), and flow velocity ( *U*, m/s). Flows in supercritical disequilibrium are observed to form sustained turbidity currents [64]. The turbulence-supporting flow velocity of layers is reported, based on bottom orbital velocity and ambient currents. For the wave characteristics, Airy linear wave theory was employed. Water properties were not relevant to this calculation; the work of the gravity flow is based on density contrasts due to the suspended sediment.

The primary criterion for ignition was the Knapp–Bagnold criterion (Equation (14)b from [59]), which approximately relates the necessary energy balance (*US*)/*ws* > 1, formulated with the seabed gradient (*S*, non-dimensional) and the grain settling velocity ( *ws*, m/s). Note that other criteria involving sediment and water entrainment (Equation (16) from [59]) apply to later flow-stages and are less relevant to initial ignition. The modeled events which satisfied the criteria were logged with their associated parameters, and collated onto a mapping (Figure 6A).

Mass failure sources: Large-scale mass failure events are also known to yield or transform into turbidity currents that can travel much further and faster than the original failure structure or debris flow [65,66]. The *WaveSlip* submodule tested for a wave-induced mass failure of seabed sediment during storms based on the circular failure approach [67]. It proved an array of plausible failure arcs, depths, and footprints. Cyclic force moments for each wave period were combined with gravitational

moments, and those driving forces were balanced against resisting ones (e.g., gradient shear strength) to test for mass failure. The complex interplay between wave-induced pressures, the footprints of loadings, sub-bottom depths of possible failure arcs, and the gravity-driven and wave-driven moments was integrated into a seafloor Factor-of-Safety (*FoS*) where values less than unity imply instability—the situation of interest.

The possibilities of mass failure were explored for three particular seabed conditions: (i) static undrained conditions of intact shear strength; (ii) remolded shear strengths considering cyclic wave-induced shear strains in the bottom; and (iii) in the presence of liquefaction, especially in shallow waters under long-period surface-waves. For (i), the static shear strengths were calculated using look-up values for the Mississippi Delta area [67]. Remolded values (ii) were computed from those based on wave-induced strains (after [68]). They were scaled linearly against a full remolding to 30% of the intact strengths occurring at 15% cumulative strain. To assess liquefaction potential (iii), a dedicated submodule *LiquiSlip* compared results using previous analytical solutions i.e., [69,70]. Significant wave heights and periods (*Hs*, m; *Tp*, s) were assumed to hold for more than 100 wave cycles, and were extracted from *WaveWatch III*® data for each modeled location and time. (Cases of breaking waves were excluded from the analysis; see [70]). Required values for seabed porosity, cohesion, permeability, and relative density (after [71,72]) were calculated based on surface sediment type from *dbSEABED*. The sediment thickness, for which only sparse sub-bottom data exists, was assumed to be effectively infinite. Note that this assumption will not apply in areas <30 m water depth where a "basal, erosional unconformity" at approximately 10 m sub-bottom marks the presence of a firm foundation under Holocene sediments (see [73]). Our study excluded such shallow areas. Time-series of the essential parameters were plotted (not shown) for selected sites in the area to monitor the *WaveSlip* and *LiquiSlip* modeling components.

After the modeling, which took place through the approximately 9 million cell spatial-temporal domain of the project, events at the lowest slope-stability *FoS* were collated and plotted, culminating in the mapping of failure predicted events (Figure 6B). All modeled mass failure events were indicated as potential sites of associated turbidity current ignitions, and their details were passed to the *RANS*/*TURBINS* component.

#### 2.2.6. *RANS*/*TURBINS:* a RANS Sediment Gravity Flow Model

*TURBINS* [74,75] solves the incompressible Navier–Stokes equations in the Boussinesq limit with a convection-diffusion equation for the sediment concentration of small, polydisperse particles whose density significantly exceeds fluid density [76,77]. As a three-dimensional, time-dependent model, *TURBINS* provides spatially and temporally resolved information about the turbulent velocity and sediment concentration fields, conversion of potential into kinetic energy, and the dissipation of this kinetic energy neglecting the effects of rotation. The dispersed phase is assumed to be sufficiently dilute so that the momentum equation governs the two-way coupling between the fluid and particles; the effect of particle loading in the continuity equation is neglected, as are particle interactions such as hindered settling. Particles are assumed to have an aerodynamic response time much smaller than typical fluid flow time scales [78]. Hence, the particle velocity is given by the sum of the fluid velocity and the constant settling velocity. Polydisperse distributions are implemented by considering different particle size classes, each assigned a settling velocity, and contributing to the overall fluid density distribution. Though there is a potential for non-Newtonian dynamics in the dense suspension region near the seafloor, *TURBINS* includes Newtonian fluid dynamics enabling the erosion and resuspension boundary conditions used within the gravity flow module.

An empirical formula to represent the resuspension flux of sediment into the current [79] has been used to estimate erosion in low Reynolds number simulation of turbidity currents [10]. A variation of this was implemented in the non-hydrostatic *RANS*/*TURBINS* code. The sediment flux due to erosion was introduced into the current as a diffusive flux from the bottom wall.

$$-\frac{1}{\text{Sc }\text{Re}} \frac{\partial \mathcal{C}}{\partial \eta} = u\_s E\_s \tag{1}$$

where *c* is the non-dimensional concentration of the sediment, η is the coordinate along the direction normal to the boundary, *us* is the settling velocity, *Es* is the resuspension flux, *Sc* is the Schmidt number, and *Re* is the Reynolds number. Based on [79], the resuspension flux, *Es*, was evaluated using

$$E\_s = \frac{1}{\mathcal{C}\_0} \frac{aZ^5}{1 + \frac{a}{0.3}Z^5};\tag{2}$$

with *a* = 1.3 × <sup>10</sup>−7, *C0* is the initial volume fraction of the sediment, and *Z* is the erosion parameter. A maximum of *0.3*/*C0* caps the resuspension flux. The erosion parameter, *Z,* is calculated as

$$\begin{aligned} Z &= 0.586 \frac{\mu\_{\text{t}}}{\mu\_{\text{s}}} Re\_p^{1.23} \text{ if } Re\_p \le 2.36; \\ Z &= \frac{\mu\_{\text{t}}}{\mu\_{\text{s}}} Re\_p^{0.6} \text{ if } Re\_p > 2.36; \end{aligned} \tag{3}$$

where *u\** is the shear velocity at the bottom wall,

$$\frac{\mu\_t}{\mu\_\*} = \frac{1}{\kappa} \log \left( \frac{\eta \mu\_\*}{\nu} \right) + B;\tag{4}$$

and *ut* is the tangential velocity at the first grid point off the wall, η is the wall-normal distance of the first grid point from the bottom wall, ν is kinematic viscosity, and constants κ = 0.41 and *B* = 5. Using *dp* as the particle diameter, ρ*p* as sediment density, ρ*0* as water density, and *g* as gravitational acceleration; the particle Reynolds number, *Rep,* is defined as:

$$Re\_p = \frac{d\_p \sqrt{gd\_p(\rho\_p - \rho\_0)/\rho\_0}}{\nu}.\tag{5}$$

As a proof-of-concept, *TURBINS* was used to represent a turbidity current generated by a lock-release (results in Section 3.3). The lock-release type simulation extended over a 21 km long domain in the streamwise (along the pathway) direction. In the vertical direction, the water depth varied from 130 m to 300 m. Dictated by a minimum resolution criterion of at least ten grid nodes over the current height, along with the condition that the grid spacing is similar in all directions, the simulation employed a grid spacing of 3 m in all directions. Consequently, the computational grid applied 7000 nodes in the streamwise direction, 100 nodes in the vertical direction, and 10 nodes in the spanwise direction. A time step of 0.6 s was used, based on a modified *CFL* (Courant–Friedrichs–Lewy) condition (*CFL* <0.5) involving both convective and viscous terms.
