Data-Driven, Multi-Model Workflow Suggests Strong Influence from Hurricanes on the Generation of Turbidity Currents in the Gulf of Mexico

Abstract

Turbidity currents deliver sediment rapidly from the continental shelf to the slope and beyond; and can be triggered by processes such as shelf resuspension during oceanic storms; mass failure of slope deposits due to sediment- and wave-pressure loadings; and localized events that grow into sustained currents via self-amplifying ignition. Because these operate over multiple spatial and temporal scales, ranging from the eddy-scale to continental-scale; coupled numerical models that represent the full transport pathway have proved elusive though individual models have been developed to describe each of these processes. Toward a more holistic tool, a numerical workflow was developed to address pathways for sediment routing from terrestrial and coastal sources, across the continental shelf and ultimately down continental slope canyons of the northern Gulf of Mexico, where offshore infrastructure is susceptible to damage by turbidity currents. Workflow components included: (1) a calibrated simulator for fluvial discharge (Water Balance Model - Sediment; WBMsed); (2) domain grids for seabed sediment textures (dbSEABED); bathymetry, and channelization; (3) a simulator for ocean dynamics and resuspension (the Regional Ocean Modeling System; ROMS); (4) A simulator (HurriSlip) of seafloor failure and flow ignition; and (5) A Reynolds-averaged Navier–Stokes (RANS) turbidity current model (TURBINS). Model simulations explored physical oceanic conditions that might generate turbidity currents, and allowed the workflow to be tested for a year that included two hurricanes. Results showed that extreme storms were especially effective at delivering sediment from coastal source areas to the deep sea, at timescales that ranged from individual wave events (~hours), to the settling lag of fine sediment (~days).

1. Introduction

The Gulf of Mexico continental margin generates >1.7 million barrels of oil per day, through >3500 oil platforms. The northern Gulf of Mexico houses >45,000 km of underwater pipes that may be exposed to structural damage from extreme oceanic events. During the passage of a hurricane, storm waves can exceed 10 m in height, resuspending seafloor sediment and potentially liquefying the seafloor. Both of these mechanisms may induce sediment turbidity currents, and in fact, ~5% of the underwater petroleum pipes appear to be broken or damaged by sudden powerful turbidity currents (BOEM pers. comm. 2015). For example, in 2004, a large sediment failure in the wake of Hurricane Ivan toppled an oil platform offshore of the Gulf of Mexico and moved it ~0.17 km downslope, initiating oil and gas leaks at a water depth of 140 m [1]. Leakage from such offshore oil and gas infrastructure puts at risk about 40% of the USA’s coastal and estuarine wetlands, which are vital to recreation, agriculture, and a $1B/y seafood industry [2].

Turbidity currents are important transport mechanisms in submarine canyons [3,4], such as the Mississippi and the De Soto Canyons, which incise the continental slope offshore of the Mississippi Delta. Several processes have been shown to have the potential to generate turbidity currents, including physical oceanographic mechanisms. Internal wave breaking on the upper slope may mobilize seafloor sediment [5]. Wave-current interactions on continental shelves during large oceanic storms can initiate wave-supported gravity flows [6]. Continental slope deposits may experience sediment failure triggered by sediment loading and over-steepening, and aided by excess pore pressure brought on by ground accelerations [7,8]. Localized events may grow into sustained currents via a self-amplifying ‘ignition’ process with accelerating erosion and entrainment of sediment from the seafloor [9,10].

While the relative importance of these mechanisms in the northern Gulf of Mexico remains to be seen, evidence points to the potential for oceanic storms to mobilize sediment there, either during the passage of moderate storms [11] or more extreme events such as hurricanes [12]. Analysis of sediment deposits indicated that most (~75%) of the sediment budget of the Mississippi Canyon could be attributed to delivery during major hurricanes, likely through gravity-driven transport [13]. Several processes affect the seafloor during short-lived hurricane passages, including sediment mass failures, erosion, and suspension. For example, mudflows in the Mississippi Delta area, triggered by the 1969 Category 5 Hurricane Camille, destroyed the offshore platform SB-70B. The seafloor at a depth of about 90 m moved more than 1000 m downslope with soil flows up to 30 m in thickness [14]. Seafloor shear stresses from waves and currents of up to 1 N/m² were monitored at a depth of 90 m during the 2004 Category 5 Hurricane Ivan, reaching the critical shear stress for fine gravel [15]. The Ivan event lifted suspended sediment as high as 25 m in the water column and eroded the seafloor up to 0.30 m vertically over more than 500 km², thus removing hundreds of millions of tons of sediment with deposits at the shelf edge and upper slope [16], and additionally causing apparent damage to oil infrastructure [1]. Evidence of the effects of large storms at great depth in the Gulf of Mexico has been seen in conjunction with other hurricanes, such as Hurricane Georges in the Mississippi Canyon [12]; Hurricane Frederic in the De Soto Canyon [17]; and Hurricane Allen [18]. Rapid loading of sections of the seafloor locally enhances the prospects for gravitational slope failures, given the associated rapid increase in pore pressures and reduction in effective sediment strengths [7]. Process-based numerical modeling offers a way to study such ephemeral high-energy processes.

Studies of sediment dispersal on continental margins, including the northern Gulf of Mexico, have typically focused on an individual component of the transport path such as gravity-driven transfer via canyons, shelf resuspension, or flood plume dispersal. For example, numerical models for suspended sediment transport have been developed and applied to the northern Gulf of Mexico [19,20,21,22], but these types of suspended transport models have not been directly linked to turbidity current models. This paper describes a numerical capability to simulate the transport of sediment, from fluvial sources, to the continental shelf, the deeper continental slope, and ultimate depocenters. Accounting for these sediment transport pathways, and the hazards that they present, is a problem of multi-scale physics, ranging from continental-scale drainage basins that deliver sediment to the sea, to shelf-wide storm systems that mobilize and redistribute sediment, to small-scale turbulent motions that affect turbidity current generation and structure.

This paper describes a loosely coupled numerical workflow that has been developed to address land-sea pathways for sediment routing of terrestrial and coastal sources, across the continental shelf, and ultimately down the continental slope and canyons of the northern Gulf of Mexico. Few studies have attempted to integrate the various transport mechanisms into a single comprehensive framework, accounting for the multi-scale physics that are relevant to the full sediment transport pathway. The workflow was used to explore conditions that may trigger episodes of sediment transport onto the continental slope and to evaluate two hypotheses: (1) episodic sediment transport down a submarine canyon is fed by sediment input at the canyon head from wave and current resuspension, and (2) turbidity currents are triggered by failures near the shelf-slope break and are likely to pass into the canyons of the continental slope. Simulation results were based on oceanographic and meteorological conditions that could impact the generation of turbidity currents. The workflow (Figure 1) includes modules that:

(1) Simulate the fluvial delivery of water and sediment into the Gulf of Mexico with the Water Balance Model-Sediment (WBMsed) and as augmented by USGS (US Geological Survey) and USACE (US Army Corps of Engineers) gauged river data;

(2) Develop domain grids and bathymetry for ocean circulation and sediment transport models;

(3) Compute spatial griddings of seabed sediment texture from dbSEABED, and of topographic channelization from the bathymetry, for use in sediment transport and seabed failure models;

(4) Employ a high resolution (10 km) spectral wave action model (WaveWatch III®) driven by GFDL-GFS (Geophysical Fluid Dynamics Laboratory–Global Forecast System) winds for use in the ocean and sediment transport models;

(5) Calculate hourly-timescale ocean circulation at a spatial resolution of a few kilometers via the Regional Ocean Modeling System (ROMS) forced with ECMWF (European Centre for Medium-Range Weather Forecasts) ERA (ECMWF Re-Analysis) winds;

(6) Represent seafloor resuspension and transport at the same resolution as ROMS’ hydrodynamics using the Community Sediment Transport Modeling System (CSTMS);

(7) Apply seabed mass-failure and a sediment suspension model (HurriSlip) to determine failure and ignition locations, and the conditions to be used as input to the turbidity current model;

(8) Develop and deploy a Reynolds-averaged Navier–Stokes (RANS) model (TURBINS) to route sediment flows down the Gulf of Mexico slopes and canyons, providing estimates of bottom shear stress needed for ascertaining possible damage to offshore infrastructure.

2. Materials and Methods

Section 2.1 describes the northern Gulf of Mexico, where the model workflow was applied. Section 2.2 provides descriptions and methods for each component of the workflow, noting how the various components can interact with one another. Section 2.3 outlines the implementation of the suite of models used to evaluate sediment routing in the northern Gulf of Mexico.

2.1. Environmental Setting

The Mississippi River drains 41% of the continental United States before entering the northern Gulf of Mexico (Figure 2). The discharge of the Mississippi River is regulated so that approximately 70% of it enters the Gulf through its main Mississippi River channel, while the remaining 30% enters through the Atchafalaya River channel [23]. Average modern-day sediment loads of the Mississippi and Atchafalaya Rivers are 115 and 57 Mt/yr, respectively [23]. Sand is deposited near the river mouths while most of the remaining suspended silts and muds are dispersed more widely [19,24,25]. Rapid delta progradation during the Holocene has narrowed and steepened the continental shelf (~20 km wide, ~0.4° gradient). The Mississippi Canyon, which cuts into the continental slope to the west of the bird-foot delta has been implicated as a conduit for shelf sediment during large storms [12,26].

A fair amount is known about suspended sediment dispersal on the Gulf of Mexico continental shelf. Frontal systems that occur frequently during winter months can create energetic waves and currents that cause significant sediment transport [27,28]. Wave contributions dominate the bed stresses on the continental shelf offshore of the Mississippi Delta, but fairweather waves are typically capable of mobilizing the seabed only in the surf and nearshore zones [19]. During extreme oceanic storms, however, deep-water wave heights exceed 10 m, with nearshore waves east of the bird-foot delta reaching 9 m in 15 m of water during Hurricane Ivan [29]. Storm waves, either from moderate storms or intense but infrequent hurricanes, have been shown to mobilize sediment mass failures on the Mississippi River Delta Front at water depths of ~75 m [11]. Sediment trap data and allied mooring and camera data from deep-water locations (~1000 m) have indicated that frequent, small magnitude resuspension events driven by inertial currents contribute to sediment transport there [30]. Less is known, however, about the mechanisms that drive shelf–slope sediment exchange or transport down the continental slope or canyons.

2.2. Workflow Components

Section 2.2.1, Section 2.2.2, Section 2.2.3, Section 2.2.4, Section 2.2.5 and Section 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³/s with a standard deviation of 11,211 m³/s (USACE observations). The WBMsed estimate of the combined Atchafalaya–Mississippi discharge for the same period was 18,300 m³/s, with a standard deviation of 12,400 m³/s. The total predicted (1995–2011) discharge for all northern Gulf rivers was 22,800 m³/s, with a standard deviation of 15,400 m³/s. The merged discharge of non-Mississippi rivers was 4540 m³/s or 19% of the total flow into the northern Gulf of Mexico, with a standard deviation of 4730 m³/s (Figure 3C). The 17 y one-day high of these non-Mississippi sources was 30,000 m³/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³/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.

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 efficient 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 (S₂, M₂, K₁, O₁ 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 runoff, 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.

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 diffusion) using a method [42] that allows constraints in the smoothing minimization like preserving the bathymetry in specific grid cells (e.g., on the continental 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.

2.2.3. Spatial Seabed Datasets

The dbSEABED facility [46,47,48] supplied information on the spatial distributions of seabed sediment type based on interpolations of more than 10⁵ 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.

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 (τ_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. 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].

Sediment Class	Source	Sediment Type	D (mm)	Tcr (Pa)	ws (mm/s)
1	Seabed	Mud	0.063	0.11	1.0
2		Sand	0.125	0.13	10.0
3		Gravel	10.0	10.0	70.0
4	Mississippi River	Small Mud	0.015	0.11	0.1
5		Larger Mud	0.063	0.11	1.0
6	Atchafalaya/Mobile Rivers	Small Mud	0.015	0.03	0.1
7		Larger Mud	0.063	0.03	1.0

). 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.

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 (H_1/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 different 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=\left(g R C h\right)∕U^2$, which depends on gravitational acceleration (g, m/s²), 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 $\left(US\right)∕w_{s }>1,$ formulated with the seabed gradient (S, non-dimensional) and the grain settling velocity (w_s, 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 (H_s, m; T_p, 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}{Sc Re}\frac{\partial c}{\partial \eta }=u_s{E}_s$$

(1)

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

$$E_s=\frac{1}{C_0}\frac{a{Z}^5}{1+\frac{a}{0.3}{Z}^5};$$

(2)

with a = 1.3 × 10⁻⁷, C₀ is the initial volume fraction of the sediment, and Z is the erosion parameter. A maximum of 0.3/C₀ caps the resuspension flux. The erosion parameter, Z, is calculated as

$$\begin{array}{c}Z=0.586\frac{u_{*}}{u_s}R{e}_p^{1.23}\mathrm{if}R{e}_p\le 2.36;\\ Z=\frac{u_{*}}{u_s}R{e}_p^{0.6}\mathrm{if}R{e}_p>2.36;\end{array}$$

(3)

where u_* is the shear velocity at the bottom wall,

$$\frac{u_t}{u_{*}}=\frac{1}{\kappa }\log \left(\frac{\eta u_{*}}{\nu }\right)+B;$$

(4)

and u_t 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 d_p as the particle diameter, ρ_p as sediment density, ρ₀ as water density, and g as gravitational acceleration; the particle Reynolds number, Re_p, is defined as:

$$R{e}_p=\frac{d_p\sqrt{g{d}_p\left({\rho }_p-{\rho }_0\right)/{\rho }_0}}{\nu }.$$

(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.

2.3. Modeling Approach

To account for the multi-scale physics of sediment delivery from rivers to the Gulf of Mexico, and subsequent mobilization by oceanic flows, the workflow was designed to operate as follows. Each model component was designed to deliver needed model inputs to the “downstream” models in a one-way coupling framework (Figure 1). Phasing of model development required coordination among the subject matter experts who developed various components of the workflow. The river discharge model (WBMsed) can provide values needed as input to the hydrodynamic ocean model. ROMS can use these discharges as point sources of freshwater and sediment, distributing the output from WBMsed for individual rivers onto the three-dimensional grid of the hydrodynamic model. For example, WBMsed provided Mississippi River discharges that were distributed to 39 Mississippi River discharge grid cells that were spread around the bird-foot delta of the ROMS Grid-f. Then, using input winds and open boundary conditions from the lower resolution Gulf model (Grid-g), the local model (Grid-f, see Figure 4) was used in the CSTMS to estimate the dispersal and deposition of sediment delivered from the rivers. The bed stresses calculated by ROMS accounted for wave-current bed shear stress, and along with WaveWatch III® data, could be employed by the HurriSlip modules to identify times and locations of sediment mass-failure and density-flow ignition. These events detected by HurriSlip, could be used to trigger a turbidity current calculation via RANS/TURBINS; which would also be informed by topographic gradients, the sediment properties from dbSEABED, and near-bed current velocities and sediment depositions calculated by the ROMS/CSTMS.

3. Results

The sections below describe model calculations from components of the workflow to demonstrate their capabilities.

3.1. Suspended Sediment Transport

The ROMS/CSTMS ocean model calculated current velocities, bed stresses, suspended sediment fluxes, and erosion/deposition from October 1, 2007, through September 30, 2008. CSTMS results for 2007–2008 indicated that the overall signature of sedimentation calculated from suspended sediment was deposition near fluvial sources, with patchy erosion and deposition elsewhere (Figure 8A). Sediment delivered by the Mobile River was largely retained within Mobile Bay. The Atchafalaya River sediment was deposited near the delta, but resuspension events on the inner shelf (depths < 30 m) created westward sediment transport along the coast. Mississippi River plumes more widely dispersed sediment around its bird-foot delta with some of the river load deposited in deeper water (>200 m). The model indicated that the deep sea experienced strong intermittent currents capable of mobilizing sediment, termed benthic storms [80].

Analysis of suspended sediment delivery to the continental slope indicated that about 70% resulted from delivery during low-intensity storms such as frontal systems, and fallout from the Mississippi River plume. The remaining 30% of the year-long delivery of suspended sediment to the continental slope occurred rapidly, during the days surrounding the passage of Hurricanes Gustav and Ike. This supports our first hypothesis, that episodic sediment transport down a submarine canyon is fed by sediment input from wave and current resuspension. Shelf erosion during non-hurricane times accounted for a small fraction of the cumulative erosion seen for the year (Figure 8B). The patchiness of erosion seen in the deep sea (Figure 8B) corresponded to the sediment texture assumed by the model (see Figure 6A). Hurricanes Gustav and Ike created widespread erosion on the shelf, and this material contributed disproportionately to sediment delivery from the shelf to the slope, compared to other resuspension events during the preceding eight months when elevated Mississippi River discharge also occurred (Figure 8). Bed shear stresses during the hurricanes were sufficient to suspend fine-grained sediment across the shelf break. Hurricanes Gustav and Ike produced distinct patterns of erosion and deposition (Figure 8C,D), mainly due to their differences in strength, duration, and storm track (see Figure 4). In general, Ike created higher bed stresses, sediment concentrations, and erosion; but in some locations, Gustav had more impact.

Suspended sediment fluxes along three cross-slope transects (locations shown on Figure 5D) were analyzed to evaluate the phasing and magnitude of hurricane-driven sediment delivery to the continental slope and beyond. The size of sediment flux generally decreased with water depth across the continental slope (Figure 9). While peak fluxes on the continental shelf coincided with the passage of hurricanes, there was often a lag of several days before suspended sediment reached deeper waters. Along the western continental slope, suspended sediment fluxes were larger for Hurricane Ike than Gustav and decreased in the deeper waters (Figure 9A). The peak suspended sediment fluxes at depth (~1300 m) occurred several days after the peak fluxes calculated for the shelf–slope break (~128 m). For the De Soto Canyon, the model estimated net downslope flux towards the south to southeast during the hurricanes (Figure 9F–I). Suspended sediment flux actually increased offshore there between depths of ~650–1115 m during and after the passage of the hurricanes, suggesting that suspended sediment transported over the sides of the canyon settled to the bottom boundary layer and contributed to the suspended sediment fluxes calculated within the canyon (Figure 9F–I).

Sediment fluxes down the Mississippi Canyon lagged behind the passage of the storms, being delayed by 1–5 days relative to the occurrence of peak wave energy on the shelf (Figure 9B–D). The model results showed that these lags corresponded to the time needed for nepheloid layers generated by cross-shelf transport of storm resuspension to be carried to, and settle into, continental slope depths. For example, Hurricane Gustav made landfall on 1 September 2008. For the Mississippi Canyon transect, the model indicated that Gustav created peak sediment fluxes on the outer continental shelf (water depth 98 m) around 2–3 September; while offshore, sediment fluxes did not peak until 4 September (688 m depth) and 8 September (1008 m depth) (Figure 9B,C). The distance along the Mississippi Canyon transect from the 98 m deep site to the 1008 m deep site is about 66 km, so the ~4.5 day lag in delivery to the 1008 m site can be explained by an average horizontal transport velocity of about 0.16 m/s. Similarly, vertical settling delays a storm’s impact on deep-sea locations. The fine sediment classes used in the model would settle about 10 or 100 m per day, so that fall out from nepheloid layers would require days to weeks to reach the near-bed continental slope and deeper. This process is illustrated using modeled suspended sediment concentrations along the Mississippi Canyon when Gustav was centered over the Lousiana shelf (Figure 10A), and five days later, the nepheloid layer was delivered to, and settled into, continental slope depths (Figure 10B).

3.2. Density Flow Ignitions

Results from the HurriSlip model suggest that during extreme storms, bed stresses are large enough to create conditions suitable for the ignition of turbidity currents from near-bottom layers of suspended sediment, especially in areas near the shelf break (Figure 6A). The modeling also suggests that small-thickness sediment mass-failure events, which may evolve into turbidity currents, are widespread around the shelf-slope transition under hurricane conditions (Figure 6B). There is some association between predicted ignitions’ locations, and the geomorphic channelizations of the upper continental slope (Figure 6B).

Sediment resuspension sources: The results on the resuspension of sediments into bottom waters (SuspendiSlip) indicated suspended sediment concentrations (SSC) during times of wave activity averaged ~300 ppm v/v, up to ~5000 ppm v/v (5% v/v) at levels 1 m above the bottom. During the storm events, in shoreface areas including at the delta front, some wave-induced bottom orbital velocities >4 m/s were indicated. At depths of 20–40 m this was reduced to >2 m/s. As modeled, wave-induced resuspensions occurred down to water depths of 189 m (at surface wave periods >13 s) in areas not sheltered from the storm wave effects.

Numerous density-flow ignition events were indicated. They were overwhelmingly in the bottom 1–2 m of the water column, but occasionally occupied water masses as thick as 8 m or more. Bulk densimetric Richardson Number values for the bottom flows ranged widely, but during the storm events were <<1.0 near-bottom i.e., were supercritical states susceptible to the onset of density flow [64]. The Knapp–Bagnold criterion discriminated events more closely and with the gravity influence of slope, identified locations of plausible density flow ignition (Figure 6A). There is some indication that suspension events in the waxing and waning of a storm are more likely to ignite because of the balance between densities and velocities.

Mass failure sources: In agreement with the extensive evidence of mass sediment failures in the region [11,81], the modeling indicated a potential for seafloor failures due to the combined effects of intense storm wave activity, shallow depth, and significant slope. The present prediction with WaveSlip, however, also extends over sandy areas not only the mudslide province at the Mississippi Delta front. There seems to be an increased potential for the failure to transform into turbidity current in sandy sediments [82].

Wave-induced liquefaction was predicted in the modeling for conditions of <30 m water depth, somewhat sandy sediments, surface wave wavelengths of >150 m, and significant wave heights of 10 m. Developed (residual and momentary) normal pore pressure increases to exceed normal overburden pressure were modeled down to subbottom depths of 10 m and more at some locations. In those circumstances, effective shear strength was reduced to near zero. The possibility of cyclic strain reduction of shear strengths was also investigated. However, the cumulative strains induced by waves, even during extreme events, were insufficient to produce significantly lowered (remoulded) shear strengths, the strains being at most of order 10⁻² cumulative (10⁻⁴ to 10⁻⁶ per cycle).

The circular-slip analyses indicated mass-failure instabilities (FoS <<1.0) over broad areas of sloping seafloor in the top 0.5 m of the seabed (Figure 6B). More deeply-seated failures, down to 20 m sub-bottom, were predicted at a small number of locations at about 30 m water depth. Still, all had a FoS >>2 and, therefore, apparently limited potential for actual failure. (They are also at the limits of the analysis in terms of wave-breaking and infinite sediment column.) HurriSlip results appear to suggest that without liquefaction or remoulding, probably very few significant wave-induced mass failures would occur in the region. However, the smaller occurrences which are also predicted, remain as candidates to release turbidity flows. They include particularly, many locales with a high likelihood of failure (FoS <<1) during storms, in seabed areas down to 100 m water depth, with a significant slope, and often near to the shelf edge.

3.3. Turbidity Currents

As a proof-of-concept for simulating turbidity currents in the northern Gulf of Mexico, we employed TURBINS for a lock-release type simulation for a narrow slice along a specific pathway in the failure location. Initially, the normalized sediment concentration, a proxy for excess density due to suspended particles, was set to one in the lock-region and zeroed elsewhere (Figure 11A). When the lock was released, a gravity-flow with a height of 30 m started to travel down the slope, so at 4 h past the lock release, the current was about 10 km downslope (Figure 11B). The current became diluted as a result of entrainment of ambient ocean water, and its height was reduced below the original height of 30 m as it traveled down the slope. Within about 8 h, the current had traveled 15 km down-slope (Figure 11C).

The velocities resulting from the momentum balance are shown in Figure 12 for different stages of the turbidity current. With time, as the turbidity current traveled downslope, thinned and became diluted; its velocities decreased (Figure 12). At 4 h post-ignition, the turbidity current had speeds exceeding 1 m/s, but by 8 h post-release the velocities were much lower. As the current traveled along the bed, it generated a counter-flowing current above that moved in the opposite direction (Figure 12A,B). The calculated velocity at the front of the turbidity current decreased from over 1 m/s to about 0.75 m/s over a period of 10 h (Figure 12C).

Figure 13 displays the modeled bed shear stress at two instants in time. A substantial level of bed shear stress exists along the current length as a result of the turbidity current created by the suspended sediment, which drives the flow. As the current decelerated, the bed shear stress value decreased. These levels of bed stress exceeded the critical shear stress levels for the seabed assumed by the suspended sediment transport model (~0.1 Pa, Table 1), indicating that the gravity flows could be auto suspending, though this process was neglected in this version of the modeling workflow.

4. Discussion

Our results offer estimates of northern Gulf of Mexico sediment delivery and oceanic transport conditions, including locations for gravity flow; and routing of riverine and shelf sediment into submarine canyons. With efforts such as these, that treat multiple time- and space-scales, modeling tools can be developed to deepen our understanding of how sediment is carried from riverine sources to various oceanic sinks. The challenges of integrating various modeling approaches across different spatial and temporal scales are substantial and require further research and code development. Both physical aspects (the implementation of erosion, resuspension of complex sediments into large-scale simulations), as well as numerical challenges (two-way coupling, temporal and spatial interpolation at the boundaries between models) require an additional community effort. The treatment of physical phase-transitions, such as between wave-supported suspended sediment flows and actual turbidity currents, requires more fundamental research.

The models developed for the workflow operated over a broad range of spatial and temporal scales. For example, the RANS/TURBINS model represented relatively thin (tens of meters) turbidity currents at higher temporal (<1 s) and spatial (~3 m) resolution than afforded by ROMS’ hydrodynamic and suspended sediment transport model. As a first step toward multi-scale modeling at the spatial level, our workflow follows sediment routing from the watershed scale via WBMsed, to the continental shelf scale via ROMS, to specific sediment gravity flows via the HurriSlip modules and RANS/TURBINS. Likewise, the processes encompassed in our workflow operate over a range of temporal scales, from that of hours for the TURBINS model, to the timescale of storm fluctuations for riverine delivery, flow ignition, and suspended transport. Changes in sediment transport that operate at seasonal and interannual timescales are likewise built into our workflow by using forcing functions for weather that represent variations in winds, precipitation, and air temperatures that operate at these timescales. Barriers in applying our methods to longer timescales (i.e., longer than decadal) include both computational limits, and difficulties in assuring that subtle biases in the models and their parameterizations do not cause the calculations to drift from realistic conditions.

The model workflow presented here is sequential, with limited two-way coupling. A fairly straightforward step is to link the riverine discharge model (WBMsed) to the oceanic ROMS and CSTMS models. It would facilitate studies aimed at quantifying oceanic dispersal of fluvial sediments for poorly gauged river systems [83]. Future efforts should explore a more direct model coupling between the suspended sediment transport and gravity flow mechanisms. Within this workflow, ROMS estimates the bed shear stresses, which can be used for the flow ignition model (HurriSlip). Locations of a slope failure can trigger simulation of a gravity current (e.g., Figure 11), which moves sediment downslope. More direct coupling between these modules would account for sedimentation via suspended sediment transport within the slope failure module, and for net erosion and deposition via gravity currents within the regional scale (ROMS) resuspension model.

Regional modeling in the northern Gulf of Mexico is not trivial. Sediment transport modeling requires high-spatial-resolution models to resolve the complex and steep bathymetry. The intense coastal circulation, eddy shedding from the Gulf Loop Current [84], and sporadic strong forcing from storms and hurricanes can affect sediment transport pathways across the continental shelf and slope. Therefore, a telescoping grid approach, from coarse (kilometers) to fine (10s of meters) horizontal scales, is required to obtain viable long-term (1–10 years) and affordable computations. Within our implementation, this was realized by using a low-resolution model for the entire Gulf, telescoping to finer-resolution for the region surrounding the bird-foot delta (Figure 3A). A similar approach has been employed to represent decadal-scale sediment transport in the northern Gulf of Mexico [22].

Joint modeling and field experiments are needed to develop reliable sediment transport models for the Gulf of Mexico continental shelf and slope. Sediment depositional data with which to compare the model calculations are severely lacking, especially at the spatial scales considered here. Recent observational efforts, some motivated by the response to the Deepwater Horizon event, have shown that sediment can be mobilized in deep Gulf of Mexico locations [30,85]. Many of our workflow’s sediment transport routines were based on parameterizations for other continental shelf systems, or on laboratory measurements. To improve and gain confidence in the models developed for this workflow requires allied field and modeling studies of sediment processes for the northern Gulf of Mexico continental shelf and slope. Because field sampling during and immediately after storm events is inherently challenging, coupled models that are consistent with observed transport processes and sedimentation are needed to characterize conditions during the extreme events most likely to lead to large sediment fluxes in the deep Gulf of Mexico, and which can damage offshore infrastructure (e.g., [1]).

5. Conclusions

A model-data workflow was developed to numerically represent sediment fluxes from fluvial sources on the inner continental shelf to the continental slope. The workflow is perhaps one of the more complex ever attempted for the problem of routing sediment from coastal sources to deep-sea sinks. The range of components (Figure 1) included: (1) database frameworks for sediment texture and bathymetry of the continental shelf and slope environments; (2) hydrology framework to simulate the discharge of water and sediment for multiple (fifteen) rivers geographically distributed along the northern Gulf of Mexico; (3) an ocean modeling framework that combined output from a spectral wave-action model with ocean circulation simulations, as driven by winds, tides and solar radiation; and tuned to the seafloor environments where bottom boundary layer dynamics can be sufficiently represented including the resuspension, transport and deposition of sediment; (4) a seafloor geotechnical modeling framework able to capture the strengthening and weakening of seafloor deposits, under both ambient ocean conditions, and high intensity, short-lived hurricanes; (5) a gravity flow generator able to determine the location(s) and sediment volume(s) displaced; and (6) a high-resolution CFD model able to simulate the development of a turbidity current, including the bottom shear stresses likely to impact offshore infrastructure. The immersed boundary RANS approach, in conjunction with the multiple successive streamwise modules, appears to be well suited to perform the Gulf of Mexico turbidity current simulations over the realistic length and time scales.

The workflow was exercised to explore the conditions that trigger episodes of sediment flux on the continental slope where gas and oil infrastructure exist. Several one-way nested grids from coarse to fine were developed to simulate the hydrodynamic circulation, sediment transport, sediment failure, sediment liquefaction, and turbidity currents in the northern Gulf of Mexico. A full Gulf of Mexico ROMS domain was run to provide boundary conditions to a higher-resolution grid that better resolved bathymetric features, river runoff, and sediment transport. Ocean hydrodynamic simulations covered the period from 1 January 2000, to 31 December 2005 (spinup), and from 1 January 2006, to 31 December 2012. It allowed us to characterize sediment transport scenarios during diverse forcing events (river discharge, storms, and multiple hurricanes). We ran focused suspended sediment transport solutions from 1 October 2007, to 30 September 2008, a time period that saw very active tropical storms and major hurricanes crossing the study area.

The suspended sediment model indicated that episodic suspended transport down the Mississippi and De Soto Canyons was fed principally by sediment fluxes generated by wave resuspension on the shelf. During the two hurricanes modeled (Ike and Gustav), suspended sediment fluxes were predominantly seaward in the vicinity of the Mississippi and De Soto Canyons. Peak suspended sediment fluxes coexisted with the occurrence of the highest wave-induced bed stresses on the continental shelf, but showed increasingly long delays relative to this timing with distance down the canyon or continental slope. While hurricane conditions only lasted for two brief episodes during the one-year model run, they accounted for about 30% of the sediment delivered from the continental shelf to the slope. Delivery of sediment directly from settling from the freshwater river plume at the canyon head or over the continental slope provided a more gradual source of sediment delivery for the study period from 1 October 2007, to 30 September 2008. Plume delivery and transport during moderate-intensity frontal passages accounted for 70% of the total sediment delivered to the continental slope during the study period.

The workflow applied a newly developed ignitions model, which was used to explore some particular mechanisms for creating turbidity currents as an additional, and perhaps the major, transportation of sediments to the slope and into channelized features there. Modeling of the flows explored physical constraints on the flow velocities and forces.

On the continental slope, turbidity currents can be triggered by slope failure when storm-driven supply forces accumulation of sediment in deeper water and steeper slopes. These appeared intense enough to both erode sediment along the path of the turbidity current and to damage offshore infrastructure. Modeling efforts in the future should explore more two-way coupling along with workflows such as developed here, and take advantage of observational methods for developing model parameterizations and confirming model estimates.