OCS Study

BOEM 2016-038

U.S. Department of the Interior

Bureau of Ocean Energy Management

Headquarters

Shelf-Slope Sediment Exchange in the 

Northern Gulf of Mexico: Application of 

Numerical Models for Extreme Events

Final Report
OCS Study

BOEM 2016-038

 

Published by

U.S. Department of the Interior

Bureau of Ocean Energy Management

Headquarters

Sterling, VA

March 2016

Shelf-Slope Sediment Exchange in the 

Northern Gulf of Mexico: Application of 

Numerical Models for Extreme Events

Authors

Hernan G. Arango

David J. Robertson

Courtney K. Harris

Justin J. Birchler

Tara A. Kniskern

J.P.M. Syvitski

Chris J. Jenkins

Eric Hutton

Eckart H. Meiburg

Senthil Radhakrishnan

Prepared under BOEM

M11AC00013

By

Rutgers, The State University of New Jersey
OCS Study

BOEM 2016-038

 

Published by

U.S. Department of the Interior

Bureau of Ocean Energy Management

Headquarters

Sterling, VA

March 2016

DISCLAIMER

Study concept, oversight, and funding were provided by the U.S. Department of the Interior, 

Bureau of Ocean Energy Management, Environmental Studies Program, Washington, DC, 

under Contract Number M 11AC00013. This report has been technically reviewed by BOEM, 

and has been approved for publication. The views and conclusions contained in this document 

are those of the authors and should not be interpreted as representing the opinions or policies 

of the US Government nor does mention of trade names or commercial products constitute 

endorsement or recommendation for use. 

REPORT AVAILABILITY

To download a PDF file of this report, go to the U.S. Department of the Interior, Bureau of 

Ocean Energy Management, Environmental Studies Program Information System website and 

search on OCS Study BOEM 2016-038.

This report can be viewed at select Federal Depository Libraries. It can also be obtained from 

the National Technical Information Service; the contact information is below.

U.S. Department of Commerce

National Technical Information Service

5301 Shawnee Rd.

Springfield, Virginia 22312

Phone: (703) 605-6000, 1(800) 553-6847

Fax: (703) 605-6900

Website: http://www.ntis.gov/

CITATION

Arango, H.G., D.J. Robertson, C.K. Harris, J.J. Birchler, T.A. Kniskern, J.P.M. Syvitski, C.J. 

Jenkins, E. Hutton, E. H. Meiburg, and S. Radhakrishnan. 2016. Shelf-Slope Sediment 

Exchange in the Northern Gulf of Mexico: Application of Numerical Models for Extreme 

Events. U.S. Dept. of the Interior, Bureau of Ocean Energy Management, Headquarters, 

Sterling, VA. OCS Study BOEM 2016-038. 116 p.

ABOUT THE COVER

The cover is a map of the study area identifying the locations of the bird-foot delta and 

Southwest Pass of the Mississippi River, Atchafalaya River and Bay, Wax Lake Delta, and 

Tarbert Landing (site of commonly used river gauge). Image is from MODIS on NASA’s Aqua 

satellite.
OCS Study

BOEM 2016-038

 

Published by

U.S. Department of the Interior

Bureau of Ocean Energy Management

Headquarters

Sterling, VA

March 2016

ACKNOWLEDGMENTS

The enthusiastic support of Dr. Guillermo Auad, BOEM project manager, is gratefully 

acknowledged. We would also like to acknowledge the Bureau of Ocean Energy Management 

(BOEM) and thank them for the opportunity to carry out this research and for providing the 

financial support.

PROJECT PERSONNEL

• Hernan G. Arango (Lead PI): ROMS model and nesting (Rutgers, The State 

University of New Jersey)

• David J. Robertson: Model compilation debugging and figure and report editing 

(Rutgers, The State University of New Jersey)

• Coutney K. Harris, Justin J. Birchler, and Tara A. Kniskern: Sediment modeling 

(Virginia Institute of Marine Science, College of William & Mary)

• J. P. M. Syvitski, Chris J. Jenkins, and Eric Hutton: Seabed failure and flow ignition 

modeling (INSTAAR, University of Colorado)

• Eckart H. Meiburg and Senthil Radhakrishnan: Turbity current modeling (University 

of California, Santa Barbara)
i

Contents

List of Figures ................................................................................................................ii

List of Tables .................................................................................................................iii

1. Northern Gulf of Mexico Background ......................................................................... 8

1.1 Economic Importance........................................................................................... 8

1.2 Sediment Fluxes Across the Continental Shelf and Slope .................................... 8

1.3 Geomorphological Framework.............................................................................. 9

1.4 Physical Oceanographic Setting ......................................................................... 12

1.5 Climatology......................................................................................................... 14

1.6 Sediment transport.............................................................................................. 15

2. Objectives and Hypothesis ....................................................................................... 15

3. Approach.................................................................................................................. 16

4. Numerical Models..................................................................................................... 17

4.1 Regional Ocean Modeling System (ROMS) ....................................................... 17

4.2 Community Sediment Transport Modeling System (CSTMS)............................. 18

4.3 Seabed Failure and Flow Ignition Models........................................................... 19

4.4 LES/RANS-TURBINS ......................................................................................... 35

5. Datasets ................................................................................................................... 54

5.1 Atmospheric Forcing Data .................................................................................. 54

5.2 Surface Wind-Induced Wave Data...................................................................... 55

5.3 Sea Surface Temperature Data.......................................................................... 56

5.4 River Discharge Simulation Data ........................................................................ 58

5.5 Bed Sediment (dbSEABED) ............................................................................... 64

5.6 Depositional Rates Data ..................................................................................... 65

5.7 Sediment Flow Pathways Analysis ..................................................................... 66

6. ROMS Simulations ................................................................................................... 67

7. CSTMS Simulations ................................................................................................. 79

8. HurriSlip Simulations ................................................................................................ 97

9. TURBINS Simulations .............................................................................................. 98

10. Conclusions and Recommendations.................................................................... 100

11. Appendix............................................................................................................... 102

12. References........................................................................................................... 103
ii

List of Figures

Figure 1.1: Study Area Map .................................................................................. 10

Figure 1.2: Views of Northern Continental Shelf ................................................... 11

Figure 1.3: Sediment Textures .............................................................................. 12

Figure 1.4: Tarbert Landing River Discharge ........................................................ 14

Figure 3.1: Project Workflow Schematic................................................................ 16

Figure 4.1: waveSlip Outputs ................................................................................ 22

Figure 4.2: Consolidation Properties ..................................................................... 25

Figure 4.3: Lithostratigraphy and in situ Cohesion ................................................ 26

Figure 4.4: Excess Pore Pressure......................................................................... 27

Figure 4.5: liquiSlip Output.................................................................................... 29

Figure 4.6: Study Area Gravel and Rock Maps ..................................................... 32

Figure 4.7: suspendiSlip Run-monitoring .............................................................. 34

Figure 4.8: Turbidity Currents Risk Map................................................................ 35

Figure 4.9: Gravity Current Over Flat Terrain........................................................ 42

Figure 4.10: Streamwise Velocity Fluctuations...................................................... 43

Figure 4.11: Immersed Boundary Grid Layout ...................................................... 44

Figure 4.12: Flow Grid........................................................................................... 45

Figure 4.13: Mean Streamwise Velocity Field and the Streamlines ...................... 46

Figure 4.14: Mean Streamwise Velocity at Selected Planes ................................. 47

Figure 4.15: Mean Wall-normal Velocity at Selected Planes................................. 47

Figure 4.16: Turbulent Kinetic Energy at Selected Planes .................................... 48

Figure 4.17: Predicted Velocity Profiles................................................................. 49

Figure 4.18: Streamwise Sediment Concentration at 600 Seconds ...................... 50

Figure 4.19: Streamwise Sediment Concentration at 1000 Seconds .................... 50

Figure 4.20: Streamwise Sediment Concentration at 1600 Seconds .................... 51

Figure 4.21: Sediment Concentration from Erosion at 0 Seconds......................... 53

Figure 4.22: Sediment Concentration from Erosion at 400 Seconds..................... 53

Figure 4.23: Sediment Concentration from Erosion at 600 Seconds..................... 53

Figure 4.24: Suspended Mass History .................................................................. 54

Figure 5.1: Hurricane Gustav Wave Height and Period......................................... 56

Figure 5.2: Hurricane Gustav Bottom Wave Period and Bottom Orbital Velocity .. 56

Figure 5.3: Grid-g Potential Temperature Solutions.............................................. 58

Figure 5.4: Modeled vs. Measured Water Discharge ............................................ 59

Figure 5.5: Modeled vs. Measure Sediment Load................................................. 60

Figure 5.6: Model Derived Sediment Load vs. Water Discharge........................... 61

Figure 5.7: Modeled River Output Map ................................................................. 62

Figure 5.8: Sum of All Rivers Flowing into the Study Area.................................... 63

Figure 5.9: Non-Mississippi Water Discharge ....................................................... 63
iii

Figure 5.10: Total Suspended Sediment Flux ....................................................... 64

Figure 5.11: Sediment Types Map ........................................................................ 65

Figure 5.12: Deposition Rates............................................................................... 66

Figure 5.13: Sediment Drainage Patterns ............................................................. 67

Figure 6.1: ROMS Grids Map................................................................................ 69

Figure 6.2: Grid-g to Grid-f Refinement Mesh Overlay.......................................... 70

Figure 6.3: Grid-f to Grid-b Refinement................................................................. 71

Figure 6.4: Mississippi Canyon Bathymetry .......................................................... 72

Figure 6.5: De Soto Canyon Bathymetry............................................................... 72

Figure 6.6: Initial and Lateral Boundary Conditions (SSH).................................... 73

Figure 6.7: One-way Nesting Sea-Surface Height ................................................ 75

Figure 6.8: 2008 Hurricane Track Map.................................................................. 76

Figure 6.9: Potential Temperature Latitudinal Cross-section................................. 77

Figure 6.10: Potential Temperature Longitudinal Cross-section............................ 78

Figure 7.1: De Soto Canyon Conditions................................................................ 82

Figure 7.2: Sediment Bed Fractions...................................................................... 83

Figure 7.3: Depth-averaged Currents.................................................................... 84

Figure 7.4: Final Sediment Distribution ................................................................. 85

Figure 7.5: Net Erosion and Deposition............................................................ 86-87

Figure 7.6: Cross-shelf and Slope Transects ........................................................ 88

Figure 7.7: Bed Shear Stresses ............................................................................ 89

Figure 7.8: Sediment Flux During Hurricanes Gustav and Ike ......................... 90-91

Figure 7.9: Depth-integrated Flux.......................................................................... 92

Figure 7.10: Suspended Sediment transect .......................................................... 93

Figure 7.11: Sediment Flux During Hurricane Gustav...................................... 94-95

Figure 7.12: Conditions During Hurricanes Gustav and Ike .................................. 96

Figure 9.1: RANS Pathway ................................................................................... 98

Figure 9.2: Sediment Concentration at 0 hours..................................................... 99

Figure 9.3: Sediment Concentration at 4 hours..................................................... 99

Figure 9.4: Sediment Concentration at 8 hours..................................................... 99

List of Tables

Table 4.1: CPU Run-time: Uniform Mesh.............................................................. 39

Table 4.2: CPU Run-time: Non-uniform Mesh....................................................... 40

Table 4.3: CPU Run-time: No Mesh...................................................................... 41

Table 5.1: River Gauging Sites ............................................................................. 60

Table 6.1: ROMS Grids Information ...................................................................... 69

Table 6.2: ROMS Configuration Options............................................................... 74

Table 7.1: Sediment Parameters........................................................................... 79
8

1. Northern Gulf of Mexico Background

The study area in the northern Gulf of Mexico (GOM) is described in terms of its economic 

importance, sediment fluxes, geomorphology, physical oceanography, climatology, and sediment 

transport.

1.1 Economic Importance

The Gulf of Mexico Outer Continental Shelf (GOM OCS) region is a mature offshore oil 

and gas production area with many producing leases generating more than 1.7 million of barrels 

of oil per day, through more than 3,500 oil platforms. Production on many of these leases is 

expected to continue into the next decade thus carrying into the future the associated risk of an 

oil spill. A large number of oil-rigs and other manmade structures, e.g., pipelines and platforms, 

are located on the continental slope where sudden and powerful turbidity currents are likely to 

occur. Over the last few years there has been a significant increase in the number of oil 

platforms operating in deep waters. The northern GOM currently has more than 28,000 miles of 

underwater pipes, which are all exposed to risks of various types of structural damage, most of 

them associated with extreme oceanic and atmospheric events. Of all damaged underwater 

pipes, about 5% are broken or damaged by sudden and violent cascading of sediments. The 

transport, burial, and impact of contaminated sediments within the GOM environment are of 

interest to many.

In addition to the oil and gas reserves of the GOM, Gulf wetlands are vital to recreational 

and agricultural interests, and the seafood industry (Stone and McBride, 1998). The State of 

Louisiana for example, contains about 40% of the nation’s coastal and estuarine wetlands. 

Additionally, sediment resuspension and redistribution affect nutrient cycles and productivity of 

the GOM continental shelf (McKee et al., 2004). The exchange of sediment between fluvial 

sources, near shore areas, the continental shelf, and outer shelf and slope are especially poorly 

constrained, both in terms of the magnitude and timing of fluxes.

1.2 Sediment Fluxes Across the Continental Shelf and Slope

While many studies have focused on the near field of the Mississippi and Atchafalya River 

sedimentary systems (i.e., see Wright and Nittrouer, 1995; McKee et al., 2004; Allison et al., 

2000), and on the continental slope of the northern Gulf (i.e., see Liu and Bryant, 2000), few 

studies consider the sediment routing from fluvial sources to offshore sinks (e.g. Coleman, 

1988). One under investigated aspect of the GOM dynamics is the transport of sediment away 

from fluvial sources and nearshore regions. This transport can be driven by either wave 

resuspension in response to the passage of cold fronts, or the collision of large northward 

propagating mesoscale eddies that periodically detach from the Loop Current. Even less well 

understood is the routing of sediment in response to extreme events such as hurricanes and 

unusually large discharges by the Mississippi and Atchafalaya rivers (e.g., May 2011). The 

existing knowledge on sediment transport in the northern GOM, under any oceanic or 

atmospheric condition, generally decreases in the offshore direction. Increased understanding of 

the mechanisms by which sediment is transported from sources such as the Mississippi and 

Atchafalaya Rivers across the continental shelf is needed because this material ultimately 

provides the source for catastrophic, downslope transport on the continental slope via slope 

failures and turbidity currents. Until a process-based understanding of sediment transport from 
9

source to sink is developed for the northern Gulf of Mexico, the frequency and triggering 

mechanisms for these catastrophic flows will remain unknown. Very little is known about 

sediment transport in the offshore areas of the continental shelf, slope, and rise regions. In 

particular, very little information exists on turbidity currents, their triggering conditions, speeds, 

transport, location, and geographical spread.

1.3 Geomorphological Framework

The largest river in North America, the Mississippi, drains 41% of the continental United 

States before entering the northern Gulf of Mexico. The Mississippi Delta, its associated 

wetlands, and barrier islands developed over geological timescales in response to continuous 

accumulation of fluvial sediment and reworking by waves and oceanographic currents. Over the 

past ~8,000 years the Mississippi River has built six delta complexes, with the most recent ones 

being the Balize (modern Mississippi bird-foot delta), and Atchafalaya-Wax Lake Deltas 

(Coleman et al., 1998). The Atchafalaya River has served as a distributary of the Mississippi for 

about 500 years (Fisk, 1952). Hydraulic control structures erected in 1963 regulate the 

distribution of freshwater to the main stem of the Mississippi River and the Atchafalaya, so that 

now approximately 70% of the combined discharge of Mississippi and Red Rivers enters the 

Gulf of Mexico through the modern Mississippi main stem, while 30% discharges through the 

Atchafalaya (Meade and Moody, 2010). The modern-day sediment loads of the Mississippi and 

Atchafalaya have been estimated to be, on average, 115 and 57 Mt/yr of sediment, respectively 

(Meade and Moody, 2010). Sand delivered both by the Atchafalaya and Mississippi Rivers 

deposits near the respective river mouths (Roberts, 1998; Wright and Coleman, 1974). However, 

about 70 – 80% of the sediment is finer-grained (silts and muds) suspended material (Wright and 

Nittrouer, 1995; see also Xu et al., 2011) that is dispersed more widely than the sands (Adams et 

al., 1987; Neil and Allison, 2005).

In response to high fluvial sediment discharge and relatively modest waves and tides, the 

Mississippi River has developed a river-dominated bird-foot shaped delta on the eastern 

Louisiana shelf (Fig. 1.1). High sediment supply has led to rapid delta progradation, developing 

a narrow and steep shelf (<20 km wide south of the delta, a gradient of ~0.4°; Fig. 1.1). The 

bird-foot delta provides the classic example of a “river-dominated” delta in the classification 

scheme of Coleman and Wright (1975) because the fluvial sediment supply overwhelms the 

oceanographic forces (waves and tides) that can act to redistribute coastal sediment. More 

recently, Walsh and Nittrouer (2009) classified the Mississippi / Atchafalaya systems as a 

Proximal Accumulation Delta because much of the sediment supplied by these rivers 

accumulates near the river mouths. The Atchafalaya River has been building the bird-foot 

Atchafalaya and Wax Lake deltas inside the shallow and semi-enclosed Atchafalaya Bay (Fig. 

1.1). Because only a small portion of sediment from the Atchafalaya River reaches deep water 

(Allison et al., 2000), the modern sediment supply to the middle and outer shelf is minimal and 

the shelf south of Atchafalaya Bay is broad, smooth, and gently sloped (~200 km wide, a 

gradient of 0.02°).

The continental shelf offshore of Louisiana and eastern Texas is, generally, somewhat 

gently sloped, as described for the shelf south of Atchafalaya Bay. Seaward of around the 60m 

isobath, however, the gradient steepens onto the continental slope, though a terrace does fringe 

most of the continental slope in the Gulf of Mexico. The continental slope is intersected by the 

Mississippi Canyon to the west of the bird-foot delta (Fig. 1.1). This feature has been implicated 
10

as a conduit for seaward flux of sediment, including during large storms such as Hurricanes 

Katrina and Rita (Goñi et al., 2007; Walsh et al., 2006). Salt-withdrawal mini-basins blanket the 

continental slope of the Louisiana – Texas shelf (Fig. 1.2), providing significant topographic 

complexity that has been argued to be important in the development and transport of turbidity 

currents there (Toniolo et al., 2006).

Figure 1.1: Map of study area identifying locations of bird-foot delta and Southwest Pass of 

Mississippi River, Atchafalaya River and Bay, Wax Lake Delta, Tarbert Landing (site of 

commonly used river gauge). Image from MODIS on NASA’s Aqua satellite.

Sediment texture within the northern Gulf of Mexico varies greatly (Fig. 1.3; Williams et 

al., 2006). While much of the shelf is muddy, the shoals offshore of Atchafalaya Bay are sandy. 

Also, Wright and Coleman (1974) found sand bars directly offshore of the bird-foot delta. Fluid 

muds have been the focus of recent studies on the inner shelf west of and offshore of Atchafalaya 

Bay, and have been shown to interact with surface gravity wave fields there (Sheremet et al., 

2003, 2005; Kineke et al., 2006). On the shelf, sediment texture and the availability of easily 

erodible mud seem to vary temporally. For example, Friedrichs et al. (2000) reported that a 

mobile layer of fluid mud (concentrations of 1 – 3 g/L) was observed at water depths of 15 – 20 

m during field operations in 1992 and 1993. Sediment cores taken during this experiment 

reported a “soupy mud layer” at the sediment-water interface (Friedrichs et al., 2000). 

Erodibility experiments, performed on over 100 sediment cores from the area offshore and west 
11

of the Mississippi birdfoot delta have found that erodibility offshore of the Mississippi delta 

exhibits spatial and temporal variability, probably in response to episodic delivery of fluvial 

sediment and wave-induced reworking (Xu et al. 2014). Other areas, such as along the 20-m 

isobaths appeared to be over-consolidated and difficult to erode (Xu, et al. 2014). It seems 

feasible that the erodibility of seabed sediment and the presence or absence of fluid muds depend 

upon feedbacks between themselves, and respond to sediment supply by the fluvial sources and 

the occurrence of energetic waves and wind-driven currents.

Mississippi Delta

Figure 1.2: View of northern continental slope of the 

Gulf of Mexico showing many crater-like features

associated with the movement of buried salt, driven by

massive loading of sediment eroded from the Rocky Mountains and delivered to the 

Louisiana coast by the Mississippi River. The majority of the mini-basins, as detailed in the 

right panels, are connected by turbidity current channels. Images courtesy of Lincoln 

Pratson.

Mississippi Canyon
12

Figure 1.3: (A) Sediment texture of study region showing bathymetry (vertical axis), crossshelf, and along-shelf location (number longitude/latitude points at 0.02 degree resolution each), 

and mud percentage (green color scale). (B) Distribution of sediment samples used to construct 

map in (A), data from Williams et al. (2006).

1.4 Physical Oceanographic Setting

The Mississippi River dispersal system is an example of a low oceanographic energy 

regime where resuspension is infrequent following initial deposition (Wright and Nittrouer 

1995). Walsh and Nittrouer (2009) assigned an average wave height of only 0.9m, and tide 

range of only 0.3m to the Mississippi continental margin, among the lower values of the global 

survey of rivers that they conducted. The Mississippi River mouth acts as a classical salt wedge 

estuary, whereby the freshwater plume overlies relatively undiluted seawater with sharp 

stratification between the two layers (Wright and Coleman 1971, 1974). Conditions important 

for sediment transport, however, cannot be characterized by the average behavior of the system, 

because resuspension and transport are often dominated by storm conditions when waves are 

energetic and wind-driven currents may also be strong.

Circulation on the Texas–Louisiana shelf has two distinct modes, with conditions during 

June – August differing from those during the rest of the year. Summer-time winds are generally 

from the west and thus upwelling favorable. This, combined with warm temperatures and large 

freshwater inputs from the Mississippi and Atchafalaya Rivers, stratifies the water column 

(Cochrane and Kelly, 1986; Cho et al., 1998). During non-summer months, winds blow 

predominately from the northeast, enhancing down-welling conditions and westward circulation 

(Rhodes et al., 1985). Large near-inertial variations over the Texas–Louisiana shelf have been 

attributed to high-frequency oscillations of the wind stress accompanying the passage of 

atmospheric fronts (Chen et al., 1996). Surface currents on the Texas – Louisiana shelf exhibit a

near-resonant response to the summer-time diurnal fluctuations in wind stress, which respond to 

the daily cycle of heating and cooling (DiMarco et al., 2000). Wind and baroclinic forces drive 
13

currents on the shelf, as the amplitudes of tidal fluctuations, are generally less than 0.4 m 

(Wright et al., 1997; DiMarco and Reid, 1998).

Observations for current velocities on the outer continental shelf and slope are limited. 

Moored measurements in 90 – 100 m of water show water column velocities to be weak, with 

means of < 5 cm/s directed toward the southwest, with some speeds reaching 15 cm/s (Halper 

and McGrail, 1988). These observations show suspended sediment concentrations higher when 

currents are directed offshore, indicating the advection of more turbid, shelf waters, past the 

sensors on the outer continental shelf (Halper and McGrail, 1988). Directly offshore of the 

Mississippi River bird-foot delta, a 60-day record of currents 1 meter and 5 meters above bed, 

show mean currents to be very weak (< 0.5 cm/s), with peak currents reaching ~8 cm/s (Adams 

et al., 1987). The circulation of the outer continental shelf and slope is impacted by the Loop 

Current, a warm, anticyclonic current (speeds of ~ 1 m/s) that enters the Gulf of Mexico near the 

Yucatan Peninsula and flows northward before turning to the southeast and exiting the Gulf 

south of Florida (see Hurlburt and Thompson, 1980; Oey et al., 2003). Anticyclonic eddies shed 

from the Loop Current often propagate northward, impacting currents and fluxes on the outer 

continental shelf. Eddy shedding is episodic, with events occurring every 3 to 17 months (Oey et 

al., 2003; Sturges and Leben, 2000).

Wave energy in the northern Gulf of Mexico is normally low compared to major oceans 

(Curray, 1960), though it is episodically high enough to suspend sediment during storms (Wright 

et al., 1997, Xu et al. 2011). Typical offshore waves have wave periods averaging from 3 to 8 s 

with heights rarely exceeding 1 m. Such waves are capable of mobilizing the seabed only in the 

surf and nearshore zones (Curray, 1960). Wave heights can reach 3 m in deep water, however, 

as has been reported from NDBC Buoy 42001 at a water depth of 3.3 km (see Kineke et al., 

2006). Comparison of these deepwater waves (heights ~ 3 m) to those observed in muddy sites 

of the Louisiana inner shelf (heights ~ 1.5 m) showed that waves can be strongly attenuated, not 

only by bathymetric shoaling, but also by interactions with fluidized mud on the inner shelf and 

nearshore (Sheremet and Stone, 2003; Kineke et al., 2006). During extreme events, however, 

wave energy can be much higher in the Gulf of Mexico. For example, deepwater wave heights 

exceeded 10 m during Hurricane Katrina in August, 2005 (NDBC Buoys 42001, 42003 both in 

about 3.3 km deep water), while nearshore waves recorded east of the bird-foot delta exceeded 

16 m (NDBC Buoy 42040 in 164 m of water), and 5.5 m (NDBC Buoy 42007, in 15 m of water 

offshore of Biloxi, MS; these wave heights reported by NOAA, 2006). Wave conditions can also 

be large enough to resuspend sediment during non-hurricane conditions. For example, during 

the so-called “Storm of the Century” in March, 1993, wave heights reached 3.5 m on the midshelf, and model estimates indicated that conditions were energetic enough to suspend sediment 

to water depths of ~30m (Xu et al. 2011). 

Wave contributions dominate the bed stresses on the continental shelf offshore of the 

Mississippi Delta (Xu et al. 2011), and therefore any estimates of suspended sediment transport 

require wave – generated bed stress. Because of the spatial and temporal complexity of the Gulf 

of Mexico’s wave field, a sediment transport model will require wave characteristics from a 

wave model capable of representing wave transformations from deepwater into the continental 

shelf.
14

1.5 Climatology

Seasonal wind patterns modify the circulation in coastal parts of the Gulf of Mexico. 

Seasonal variations in storm and freshwater discharge impact the sediment routing on the shelf. 

The most severe storms are hurricanes and tropical storms that occur during summer and early 

autumn, with about half of the Gulf Coast hurricanes making landfall in September (Landsea, 

1993). Hurricanes impact the Gulf via intense winds that drive strong currents and increase 

wave heights. Extratropical storms manifest with the passage of a cold front along the Gulf 

Coast. These tend to occur during the winter and early spring and typically bring energetic 

waves that promote sediment flux (Kineke et al., 2006). Arctic surges bring strong northeasterly 

winds and associated short-period waves, while migrating cyclones at first brought southerly 

winds, turning to weaker northwesterly winds after the front passed (Pepper and Stone, 2002). 

Migrating cyclones are associated with longer period waves and are therefore more likely to 

suspend sediment at greater water depths than Arctic surges. Migrating cyclones cause net 

sediment flux to the east, while the arctic surges cause net shoreward flux (Pepper and Stone, 

2002).

Mississippi River discharge, and therefore freshwater delivery to the Gulf of Mexico varies 

seasonally. Because the Mississippi drains such a large catchment, discharge tends to vary over 

relatively long timescales (weeks to months), with peak levels between mid-winter to late spring 

(Bratkovich et al., 1994). Freshwater discharge from the Mississippi and Red Rivers are 

combined and then re-divided to flow into the Atchafalaya and Mississippi Rivers, which then 

carry the freshwater to the coast. The normal diversion to the Atchafalaya is ~30% of the total 

(for details see below), though the fraction may be modified during extreme floods. The average 

discharge for the Mississippi River is about 23,500 m3

/s (Nittrouer et al., 2008). Recent flood 

years include 1993, where discharge peaked at 46,000 m3

/s (see Xu et al., 2011) and in May, 

2011 when the combined Mississippi and Atchafalaya River discharge reached 65,000 m3

/s (Fig. 

1.4).

Figure 1.4: Mississippi and Atchafalaya River discharge measured at Tarbert Landing (blue), 

and Simmsport (red) for 2010 – 2011. Data from Hydraulics Branch, CEMVN-ED-HD; New 

Orleans, http://www.mvn.usace.army.mil/Missions/Engineering/StageandHydrologicData.aspx.
15

1.6 Sediment transport

Observational and modeling efforts describing sediment transport on the continental shelf 

of the Gulf of Mexico have usually focused on the inner – to – middle shelf area, and highlight 

complex interactions between sediment availability and the hydrodynamic forces that mobilize 

the seabed. Sediment transport on the continental shelf of the northern Gulf of Mexico is 

dominated by suspended load and peaks during times of energetic waves and currents that 

combine to create large bed stresses. Observations from the outer continental shelf and slope (60 

– 100 m water depth) indicated that transport is dominated by storm events (Adams et al., 1987; 

Halper and McGrail 1988). Observations and estimates of bed stresses from much shallower 

locations (2 – 8 m depth) conclude that shear stresses are highest during storms, and cause the 

highest sediment fluxes (Pepper and Stone 2002; Jaffe et al., 1997). Sediment fluxes in the 

northern Gulf of Mexico may be heavily impacted by sediment supply from the large fluvial 

sources of the Mississippi and Atchafalaya Rivers. Based on observations from 15 – 20 m water 

depth, west of the Mississippi bird-foot delta, Wright et al. (1997) conclude that suspended 

sediment concentrations depend on local bed stresses and sediment availability. They observed 

high suspended sediment concentrations in the bottom boundary layer during times of low bed 

stresses and argued that these indicated the importance of sediment settling from the freshwater 

plumes and advecting to the field site.

Observations suggest the presence of fluid mud layers that attenuate wave heights at the 5-

m deep site on the inner shelf (Sheremet and Stone, 2003; Sheremet et al., 2005). The formation 

of fluid muds is attributed to liquefaction of the sediment bed by energetic waves, appear to 

prevail through a combination of turbulent mixing, stratification at a lutocline, and hindered 

settling (see Safak et al., 2010). These investigations demonstrate the importance of feedbacks 

between fluid mud, turbulence, and wave fields on the inner shelf, but did not establish the 

influence of cross-shelf flux of sediment, or linkages between inner shelf processes and the mid 

to outer continental shelf and continental slope. In the larger scheme of things, these fluid mud –

wave interactions appear to increase the retention of sediment on the inner shelf and aid in 

building the clinoform offshore of the Atchafalaya Bay (Jaramillo et al., 2008).

2. Objectives and Hypothesis

The primary goal of this project is to provide a comprehensive understanding of the factors 

that trigger significant across-isobaths sediment fluxes in the mid- and outer-continental shelf 

and slope in the northern Gulf of Mexico.

The objectives to accomplish the goal are:

i. To develop capability to numerically model sediment fluxes from inner continental shelf 

or fluvial sources to the continental slope;

ii. To explore the conditions that trigger episodes of sediment flux on the continental slope 

where gas and oil infrastructure exist;

iii. To study the sediment transport scenarios due to idealized extreme atmospheric forcing 

events (storms and hurricanes of various intensities and different paths);
16

iv. To study the sediment transport response to high discharge and ambient oceanic forcing 

(Loop Current eddy shedding, mesoscale and sub-mesoscale eddies, meanders, waves, 

etc).

The following hypotheses will be tested:

i. Episodic transport down features like Mississippi Canyon is fed by sediment input at the 

canyon head from either fluvial input or wave and current resuspension.

ii. On the continental slope, turbidity currents can be triggered by slope failure when storm 

driven supply causes sediment to accumulate in deep water.

3. Approach

The approach for achieving the project goal and objectives is via loosely coupled numerical 

simulations using hydrodynamic, sediment transport, and turbidity models. Several state-of-theart models were conceptually linked to address the various land-sea pathways for sediment 

routing from terrestrial and coastal sources, across the continental shelf and ultimately down the 

continental slope and canyon systems of the northern GOM. Model simulations were designed 

to represent a range of conditions and climatologies that might impact the generation of 

turbidity-current triggering mechanisms (i.e. extreme storms, high discharge, eddy shedding 

from the Loop Current, cumulative sediment deposition).

The simulation model/data workflow is illustrated in Figure 3.1:

Figure 3.1: Schematic workflow for the project in terms of model development and data 

generation and usage. The different component of the system is color coded for clarity.
17

It includes:

i. Model WMBsedV2 to simulate the flux of water and sediment discharged from rivers into 

the Gulf of Mexico with model calibration using USGS and USACE gauged river data;

ii. Develop domain grids and concomitant bathymetry (ETOPO2) for the oceanographic 

models;

iii. Augment these grids with realistic seabed sediment textures (dbSEABED);

iv. Employ a high resolution (10 km) spectral wave action model (WaveWatch III) as driven 

by GFDL - GFS winds;

v. Employ the Regional Ocean Modeling System (ROMS) for ocean circulation forced with 

ECMWF ERA winds;

vi. Model seafloor resuspension and suspended sediment transport using the Community 

Sediment Transport Modeling System (CSTMS) inside ROMS; 

vii. Develop and deploy a seabed failure and flow models (HurriSlip) to determine failure 

and flow ignition locations and conditions to be used as input in the turbidity current 

model; and

viii. Develop and deploy a Navier-Stokes Reynolds Averaged 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.

4. Numerical Models

A downscaling modeling approach using various types of numerical models is 

implemented to resolve the appropriate temporal and spatial scales associated with the 

circulation dynamics, sediment transport, and turbidity currents. The modeling framework 

includes ROMS, CSTMS, various seabed failure and flow ignition models, and LES (Large Eddy 

Simulating)/RANS (Reynolds Averaged Navier-Stokes)-TURBINS.

4.1 Regional Ocean Modeling System (ROMS)

ROMS is a mature numerical framework used routinely by the scientific and operational 

communities to study ocean dynamics over a wide range of spatial (coastal to basin) and 

temporal (days to seasons, inter-annual) scales. It is unique in that it is the only community 

framework including the adjoint-based analysis and prediction tools that are available in 

Numerical Weather Prediction (NWP), such as 4-dimensional variational data assimilation (4DVar), ensemble prediction, observation sensitivity and impact, adaptive sampling, and circulation 

stability and sensitivity analysis (Moore et al., 2011a,b,c). ROMS is freely distributed 

(http://www.myroms.org) to the Earth-modeling community and has thousands of users 

worldwide.

ROMS is a three-dimensional, free-surface, terrain-following ocean model that solves the 

Reynolds-averaged Navier-Stokes equations using the hydrostatic vertical momentum balance 

and Boussinesq approximation (Haidvogel et al. 2000, 2008; Shchepetkin and McWilliams, 

2005, 2009). The governing dynamical equations are discretized on a vertical coordinate that 

depends on the local water depth. The horizontal coordinates are orthogonal and curvilinear 
18

allowing Cartesian, spherical, and polar spatial discretization on an Arakawa C-grid. Its 

dynamical kernel includes accurate and efficient algorithms for time-stepping, advection, 

pressure gradient (Shchepetkin and McWilliams 2003, 2005), several subgridscale 

parameterizations (Durski et al., 2004; Warner et al., 2005) to represent small-scale turbulent 

processes at the dissipation level, and various bottom boundary layer formulations to determine 

the stress exerted on the flow by the bottom.

The momentum equations are solved using a split-explicit time-stepping scheme that 

requires special treatment and coupling between barotropic (fast) and baroclinic (slow) modes. 

A finite number of barotropic time steps, within each baroclinic time step, are carried out to 

evolve the free-surface and vertically integrated momentum equations. A cosine-shape time 

filter, centered at the new time level, is used for the averaging of the barotropic fields 

(Shchepetkin and McWilliams, 2005). In addition, the separated time stepping is constrained to 

exactly maintain both volume conservation and constancy preservation properties, which are 

needed for the tracer equations (Shchepetkin and McWilliams, 2005). By default, the 2D and 3D 

equations are time-discretized using a third-order predictor (Leap-frog) and corrector (AdamsMolton) scheme, which is very robust and stable. The enhanced stability of the scheme allows 

larger time steps, by a factor of about four, which more than offsets the increased cost of the 

predictor-corrector algorithm.

ROMS is a modular code written in F90/F95 and employs C pre-processing to activate the 

various physical and numerical options. The parallel framework is coarse-grained with both 

shared-memory (OpenMP) and distributed-memory (MPI) paradigms coexisting in the same 

code. Because of its construction, the parallelization of the adjoint is only available for MPI. 

Several coding standards have been established to facilitate model readability, maintenance, and 

portability. All the state model variables are dynamically allocated and passed as arguments to 

the computational routines via dereferenced pointer structures. All private arrays are automatic; 

their size is determined when the procedure is entered. This code structure facilitates 

computations over nested grids. There are three types of nesting capabilities in ROMS: (i) 

refinement grids which provide increased resolution (3:1 or 5:1) in a specific region; (ii) mosaics

which connect several grids along their edges, and (iii) composite grids which allow overlap 

regions of aligned and non-aligned grids (Warner et al., 2010).

4.2 Community Sediment Transport Modeling System (CSTMS)

The Community Sediment Transport Modeling System (CSTMS) was developed inside the 

ROMS kernel. The suspended and bed sediment concentrations can be represented using an 

unlimited number of user-defined cohesive (mud) and non-cohesive (sand) sediment classes, but 

to date, most published uses of the CSTMS have used the non-cohesive routine (see Warner et 

al., 2008). Each sediment class has attributes of grain diameter, density, settling velocity, and an

erodibility constant, which are specified in an input file and held constant for the model run. A 

multi-level bed framework tracks the distribution of every size class in each layer and stores bulk 

properties including layer thickness, porosity, and mass; allowing the computation of bed 

morphology and stratigraphy.

Typical implementations of the CSTMS use bed layers of ~mm to cm-scale thickness. 

Also tracked are bed-surface properties like active-layer thickness, ripple geometry, and bed 

roughness. Bedload transport can be calculated for mobile sediment classes in the top layer. The 
19

erodibility of non-cohesive sediment types depends on the critical shear stress for erosion (τcr), 

specified for each sediment type by the user in an input file. Suspended transport is estimated by 

assuming that each sediment class acts independently of the other classes, 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 can be included in the equation of 

state, and to the degree allowed by the model’s vertical resolution, account for gravitationally 

driven flows (see Bever, 2010; Chen et al., 2011). Net exchange of sediment between the seabed 

and suspended load are estimated by assuming simultaneous erosion and deposition via settling 

(Warner et al., 2008).

The sediment-modeling group at VIMS has actively participated in recent developments 

within the CSTMS, including a cohesive sediment model (see Rinehimer et al., 2008; Fall et al.

2014) and linking the sediment model with the biogeochemical model of Fennel (see, e.g. Fennel 

et al., 2011). A cohesive sediment formulation is now available within the CSTMS for which 

the critical shear stress for erosion of cohesive sediment is assumed to vary with time, and with 

depth in the bed, generally become less erodible (higher τcr) with age and with depth (see 

Sanford, 2008; Rinehimer et al., 2008). The instantaneous critical shear stress for a bed layer, τcr

(x,y,zbed,t), is nudged toward an assumed equilibrium value, τcr, eq (zbed), at a rate governed by 

user specified timescales of consolidation (for the case of τcr < τcr, eq) and swelling (for cases 

where τcr > τcr, eq), as described in Sanford (2008) and Rinehimer et al., (2008). Our previous 

implementations of the cohesive sediment module, for the York River Estuary, VA (Fall et al.

2014), relied on estimates of τcr, eq (zbed) obtained using Gust-microchamber erodibility 

measurements (see Gust and Muller, 1997; Thomsen and Gust, 2000; Dickhudt et al., 2009). 

Gust-microchamber erodibility measurements made during the summers of 2010 and 2011 were 

reported in Xu et al. (2014) and can be used to estimate τcr, eq (zbed) and characterize sediment 

erodibility in the northern Gulf of Mexico. Also relevant to the northern Gulf of Mexico, co-PI 

C. Harris is linking the CSTMS to biogeochemical models. Though not a focus of this project, 

links between the sediment transport model and biogeochemical model may be of interest in the 

future.

4.3 Seabed Failure and Flow Ignition Models

Several depth- and time-dependent (z, t) analytical models have been developed to

investigate the seabed response to extreme forcing events in the northern Gulf of Mexico, such as 

tropical storms and hurricanes. The focus is to determine when and where the conditions during

extreme events are suitable for initiating downslope sediment gravity flows. Information about 

seabed failure and flow ignition conditions can then be passed to other models, like TURBINS a

LES/RANS turbidity current model that would then carry the sediment mass as a suspension into 

deeper parts of the continental slope.

There are several known processes affecting the seafloor during hurricane passages

including sediment mass failures, erosion and suspension. Understandably, these extreme 

processes have rarely been observed directly but they can be interpreted from their impacts. 

Submarine mudflows in the Mississippi Delta area, triggered by the 1969 Category 5 Hurricane 

Camille, destroyed the offshore platform SB-70B. The sea floor at a depth of about 90m moved 

more than 1000m downslope with soil flows up to 30m in thickness (Sterling and Strohbeck 

1975). High seafloor shear stresses from waves and currents up to 1 N/m2

, and a critical shear 
20

stress for fine gravel (see Table 7.1 in Julien, 1998) were monitored at a depth of 90m during the 

2004 Category 5 Hurricane Ivan (Wijesekera et al., 2010). The Ivan event raised suspended 

sediment up to 25m in the water column and eroded the seafloor up to 30cm vertically over more 

than 500 km2 thus removing hundreds of millions of tons of sediment (Teague et al., 2006). The 

suspended sediment was most likely deposited at the shelf edge and upper slope of the 

continental margin. Rapid loading of sections of the seafloor from mudflow lobes (see Coleman 

et al., 1980) or deposition of suspended sediments locally enhance the prospects for gravitational 

slope failures given the associated rapid increase in pore pressures and reduction in effective 

sediment strengths (Guidroz 2009, unpubl.).

Observations of the processes affecting the seafloor sediment are very limited and likely do 

not represent the most extreme conditions. Process-based numerical modeling offers a way to

study these processes under extreme conditions. Several computational models are used to 

investigate sediment slope failure, liquefaction, and gravity flows ignition. They are used to 

bridge between the ROMS-CSTMS sediment transport and the TURBINS turbidity current 

models, and are in part based on guidance from the advanced sediment slope failure model 

within SedFlux (Hutton and Syvitski, 2004):

• hurriSlip: Offers the overarching framework from within the following modules can be 

called to simulate seabed sediment flow setup and gravity flow ignition model.

• waveSlip: Simulates the wave-induced, pressure-gradient mass failure of sediments 

(see Nodine et al., 2007).

• subSlip: Used for computing sediment slope stability parameters for a particular 

application.

• liquiSlip: Simulates sediment liquefaction and pore-pressure evolution during extreme 

events like tropical storms and hurricanes (see Jeng et al., 2006). suspendiSlip: 

Simulates the potential for near-bottom suspended sediment to initiate gravity flows

and tubidity currents, as described in Parker et al. (1986). It is a module within 

hurriSlip.

Model Descriptions

Each of the above turbidity current initiation modules employs a one dimensional (z), timedependent (t) model whose coordinates are established at grid-points from the large-scale ROMS 

grid (Grid-f, see below). The wave climate is based on the Wave-Watch III (WW3) model, also 

used as input for ROMS-CSTMS. Regional hydrodynamic flows are disregarded. The surficial 

seabed material properties are obtained from dbSEABED grids, along with slope calulations that 

are based on the bathymetry. Surface sediment types and depths are not observed to change 

significantly during the modeling.

Long-period waves are known to cause seafloor gravity failures in the Gulf of Mexico 

(Nodine et al., 2006). Measurements have demonstrated cyclic motions up to 0.5m in the 

shallow seabed under hurricane waves. Cyclic force moments are created during each wave 

period, and combine with gravitational moments to cause failure. Thus the WaveSlip module is 

based on the work of Nodine et al. (2007) and the impact of Hurricane Ivan where a limit 
21

equilibrium method was employed. One change from this earlier work was to replace a stepwise 

classification of wave type according to water depth (deep, intermediate, shallow and hence 

calculation of wavelength from period), with a continuous function following Hunt (1979).

The WaveSlip module data needs to include: (a) Neighborhood seabed topographic 

(terrain) statistics calculated using 3x3 samplings (see Burrough and McDonnell, 1998). These 

statistics are for spatial scales equivalent to twice the bathymetric grid cell size. (b) The time 

series of wave heights and periods (Hs, meters; Tp, seconds) extracted from the derived WW3 

wave data at each modeled location. This yields 3 hourly sampled records, or about 720 wave 

cycles under storm conditions. WaveSlip output parameters include water depth, slope, aspect 

and (de-trended) roughness. 

To assess the stability, a Factor of Safety (FoS) criterion is used, with limit equilibrium 

geometry (Fig. 7 in Nodine et al., 2007). The FoS is defined as

FF = MMrrrrrrrrrrrr⁄�MMgg gggg + MMww � (1)

with MMrrrrrrrrrrrr, MMgg gggg, and MMww as the moments of the resisting, gravity, and wave driving forces

respectively. A FoS<1 is indicates likely failure.

The weight of soil, calculated as 

WWgg gggg = γγ′(ππ RR2θθ − h�RR2 − h2) (2)

creates the gravity induced moment in the slope MMgg gggg = WWgg ggggaa. The variables are: γγ′

submerged unit weight of the soil, RR (m) radius of moment arm, h (m) height of the circle center 

above the seabed, and θθ (radians) angle of oscillation. The moment arm, aa (the ‘horizontal 

distance from the center of gravity of the soil weight to the center of the circle’), is calculated as 

aa = 2

3

RR �ssssss �

θθ

2

�

3

�

θθ

2 − ssssss �

θθ

2

� cccc �

θθ

2 � ��� sin ββ, (3)

with slope ββ (radians).

For the wave-induced moment component,

MMww = 2ppmm mm �kk−2ssssss(kkkk) − XX

kkcccc (kkkk)� (4)

where wavelength LL (m) is expressed as wave number (kk = 2ππ⁄LL), and the chord length (m) 

= RR �

θθ

2

�. The maximum wave-induced pressure over hydrostatic (Pa), ppmm mm , is obtained 

through linear wave theory (Wiegel, 1964; numerical approximations from Hunt 1979): 

ppmm mm = γγww

2 (HH⁄cccc h(kkkk)) (5)

with HH the wave height (m), water depth (m), γγww density of seawater (kg/m3).
22

These driving moments, MMww and MMgg gggg, are balanced against resisting moment based on 

geometry and the sediment shear strengths: 

SS = �cc0 − cczzγγ′RR �

θθ

2

�� RR + 2cczzγγ′RR2ssssss �

θθ

2

�, (6)

using cc0 as the shear strength at seabed surface, cczz as gradient in shear strength by effective 

stress, so that shear strength at level zz is cc = cc0 + cczzγγ′zz in the case with constant γγ′, else with an 

integrated solution. The resisting moment MMrrrrrrrrrrrr = SS .

Model results are realistic, though continued testing is certainly in order. The interplay 

between wavelengths, the footprint of the loadings, sub-bottom depths of possible failure arcs 

and the gravity-driven and wave-driven moments is complex. This is reflected in the FoS and 

parameter ratio values over time for a selected site (Fig. 4.1). 

Figure 4.1: Outputs monitoring the time-series results of the module waveSlip. The module 

depends on the ROMS/CSTMS wave field and bathymetric gridded data. In this case, no cyclic 

loading failures are expected with the local sediment (slightly stiff mud), gradient (6%) or depth 

(66m). (a) Wave height (Hs) and period (Tp) statistics; (b) calculated Factors of Safety (FoS), in 

this case all above 1.0; (c) maximum ratios of wave-induced moments versus total failuredriving moments for the site (are not a determinant of failure). Abbreviations: wv, wave induced, 

gv, gravity moments.
23

High energy sediment failure involving coherent materials are able to shed downslope, fast 

non-coherent (liquid) turbidity currents. However, Mohrig and Marr (2003) note with results 

from laboratory investigations that there is not yet enough understanding to formulate a more 

sophisticated model.

subSlip is developed to assess mass failures using information on the subbottom sediment 

properties. While dbSEABED holds data from many cores in the region, it was preferred to 

model the subbottom properties from first principles. subSlip is better applicable in less known 

areas such as the Northern Gulf of Mexico where local variations in seafloor properties make 

accurate prediction difficult. Important properties employed in subSlip are porosity (%) and 

shear strength (kPa) and derived quantities like density (kg/m3) and friction angles (degrees). 

The basic approach for shear strength is to have a downwards linear increase involving a 

surficial value C0 and a linear gradient Cdz. This approximation has been used elsewhere 

(Nodine et al., 2007). For this project a prediction framework of sediment physical properties 

was developed (Fig. 4.2), down to 100m below the seafloor. The parameters of interest are those 

employed for slope stability calculations, such as infinite slope models. Below is a brief 

description of these parameters:

a. The porosity profile is central to strength considerations (in absence of cementation). 

Porosity in a compacting column is controlled by coarse fraction and overburden stress. 

Near surface, sediments compact rapidly as documented in Boudreau and Bennett (1999). 

Mixed fine-and-coarse sediments generate a coarse-supported fabric that tends to restrict 

consolidation as they tighten and when the coarse fraction exceeds a certain threshold

(Winterwerp et al., 2012). Thus, Ocean Drilling (e.g., Cavin et al., 2000) and many other 

studies show that for depths down to a few hundred meters, sands compact little beyond 

the surface layers. Muds in that depth range compact appreciably from initial high 

porosity values as quantified by Bahr et al. (2001). For this project the Boudreau and 

Bennett (1999) relations were employed for the near surface (<3m). The Bahr et al.

(2001) and Cavin et al. (2000) relations were used at deeper levels. Prediction of the 

porosity profile is improved using the concepts of relative density (see Bardet 1998) and 

granular frameworks (van Ledden et al., 2004). In the shallow sediment column, the 

methods of Bahr et al. (2001) are not appropriate because the ultimate (very deep) 

porosity is not well defined. The relative density of a sediment column is defined as:

DDrr = eemm mm − ee�eemm mm − eemm mm , (7)

(x100 for percent) where eemm mm and eemm mm are the void ratios at loosest and densest states, 

and ee is the current void ratio. The eemm mm and eemm mm are estimated for surface sediment of 

known texture using the envelope of: (i) relations from Youd (1973) with moderately 

rough grain shapes, (ii) maximum and minimum porosities recorded by Jackson & 

Richardson (2007), (iii) local surface sediment values from dbSEABED.

b. Effective stress is approximately the integral of the weight of solids above a cell, 

adjusted for buoyancy (i.e., γγ′

, submerged weight). 
24

c. Cohesion for clays is defined by the effective stress as described in Hutton and Syvitski 

(2004, p. 370). The final result is prorated to all sediment by mud content with zero 

cohesion for the coarse sediment fraction. 

d. Effective Internal Friction Angle. Friction angles for surficial materials in the US Soils 

Classification (USCS) were obtained from various sources including Senneset & Janbu 

(1985), and were then synthesized into linear functions by grain size and effective stress. 

For in situ values with high mud content, the relation of Duncan et al. (1989) was used 

that relates the internal friction angle (φ’) to effective stress (σ’). See Figure 4.3.

e. Consolidation Excess Pore Pressure. Two methods were used to evaluate this 

parameter: (i) time-based, active settlement (ii) steady-state under steady small loading. 

For the first, the Gibson (1958) formulation for clays was implemented following 

Dimakis et al. (2000, Eq. 4 and their discussion), with Scipy Integrate evaluation of the 

integration. For the second method, the analytical solution of Viesca and Rice (2012, Eq. 

7) was employed. Both methods have advantages and disadvantages. The Gibson 

approach assumes a constant Coefficient of Consolidation (Cv) and is lithologyindependent but time-dependent. The Viesca and Rice method allows for varying Cv, can 

parameterize slopes, and seems compatible with the relatively modern net deposition 

rates of the shelf and upper slope. See Figure 4.4.

f. The transition between cohesive and granular behavior is assumed to be at ~8% mud 

fraction (5-13% range) and varies by water content due to the proximity of sand and silt 

grains affecting their interactions (van Ledden et al., 2004). Permeability is estimated 

using the Kozeny-Karman relation for sediments after Bardet (1998), depending on grain 

size-distributions and angularity, and void space parameters.

The physical properties data can be deployed to models via profiles at each 1D modeling 

site (Fig. 4.2) or where local geometry is important for larger-scale failures, via cross-sections 

with appropriate placement.

The possibility of over-consolidation remains an awkward issue, partly because it can be 

localized and may have several causes. The principal cause of over-consolidation is past 

overburden, leading to consolidation under a higher effective stress than existing today. Biologic 

factors (Kuo et al., 2010), wave and current shear hardening are also over-consolidation factors.

For this project over-consolidation is assumed not present until evidence shows otherwise.

Another challenge is defining the sub-bottom stratigraphic structure, which essentially is a 

matter of integrating core log and seismic profile information (e.g. Kindinger 1988), or running a 

stratigraphic simulator such as SedFlux (Hutton & Syvitski 2004). It is crucial for modules such 

as waveSlip to have the correct sediment textures for areas and layers, leading to correct shear 

strengths, etc. For some parameters, such as excess pore pressure the age of a sediment deposit 

is important (e.g., Gibson 1958).
25

Figure 4.2: Example of the modeling of shallow-subbottom consolidation properties for the 

sediment. This is for moderately sandy sediment. (Note: The parameters are plotted with a 

variety of units.)
26

Figure 4.3: (a) Plot of a simulated internal lithostratigraphy in terms of mud content from 

assignment of top-level sediment type using bathymetry, deposition rates and paleodepths. (b) 

corresponding calculated in situ cohesion (Cu’); to support models of large-scale slope stability. 

The axes are: (x) section cell, each 0.01 dg; (y) present water depth (m). The shoreline lies to 

right (north). Only zones adjacent to the seabed were computed, as shown.

a.

b.
27

Figure 4.4: Plots of (a) calculated excess pore pressure (‘uEppGibsn’ in MPa) by agesensitive, vertical Gibson (1958) method; (b) calculated excess pore pressure (‘uEppVr’ in MPa) 

by the steady-state, lithology-sensitive, slope-affected Viesca and Rice (2012) method. The xaxis is as in Figure 4.3.

a.

b.
28

liquiSlip. Liquefaction is the loss of internal sediment strength through a buildup of pore 

pressure. It may occur during cyclic deformation of marine sediments under storm waves (Finn 

et al., 1983), depending on wave characteristics, water depth and sediment thickness, and 

sediment geotechnical properties. Direct evidence of liquefaction has been observed in sidescan 

images of the submerged Mississippi Delta front (e.g., Prior et al., 1979).

The project module liquiSLIP is an implementation of the engineering practice methods in 

Jeng et al. (2007) that compute the cyclic and residual components of liquefaction under a 

modeled wave environment. The methods take into account the sediment character through 

shear modulus, Poissons ratio, permeability, and sediment thickness. Accordingly, the model is 

tuned to the varying sediment types across the project area, from fine muds to sands in thick 

accumulations to bodies near to and over bedrock.

Unfortunately the sub-bottom seabed structure is usually unknown, even in the relatively 

well-explored northern Gulf of Mexico. For this reason, and also due to different approaches in 

physics, several different formulations of the liquefaction problem were considered for the 

modeling.

a. wave-sediment interaction following Spierenburg (1987, 1989) over a rigid base after 

code in Verruijt (1995);

b. employing the classic analysis of Yamamoto et al. (1978) and de Wit & Kranenburg 

(1997); 

c. for sandy sediments (Finn et al., 1983);

d. boundary problem solution (Mei and Foda, 1981);

e. comparing pore pressures to earth pressures directly Geremew (2013), Jeng et al. (2006).

Analytical estimates for depth of liquefaction, loss of effective strength, and ratio of cyclic 

and residual pore pressures (PPoo oo, PPrrrrrr) are required. Jeng et al. (2006) proposed a simplified 

method whereby:

PPoo oo = PPbbeeee (−kk ) , and

PPrrrrrr = 2AA

vvλλ3�1−�

λλ

2

+1� ee(−λλ )� ,

(8)

where PPbb is the peak bottom pressure, and AA is a coefficient for the source term of Sumer and 

Fredsoe (2002) in the one dimensional consolidation equation. The term λλ relates to the surfacewave wavelength via empirical constants. The coefficient of consolidation is evaluated by 

cc = 2GGGG(1 − νν)

γγww(1 − 2νν) � (9)

using GG the shear modulus (Pa), the permeability (m2), vv Poissons ratio, and γγww the unit weight 

of pore fluid (N/m3). 

Loss of effective stress when PPrrrrrr = σσ′ marks a state of liquefaction. The effective stress is
29

σσ′ = (1+2KKoo)

3 γγ′zzLL (10)

where oo is the lateral earth pressure coefficient in a state of rest, and zzLL is the maximum 

liquefied depth (m), which can now be determined.

This approach to estimation of liquefaction is implemented in liquiSlip and covers cases of 

quasi-infinite sediment thickness, non-breaking waves of modest steepness, isotropic saturated 

soil, and fine sediments. Other modules in liquiSlip address somewhat different scenarios or 

physics.

From the present module implementation, this method (dashed lines, Fig. 4.5) predicts 

varying degrees of seabed liquefaction in the top 7m under extreme hurricane-grade waves in 

20m water depth in conditions similar to those of the Mississippi Delta front.

Figure 4.5: Output to monitor the results from the liquiSlip model, showing the seabed 

responses to wave conditions in terms of (log10) shear stresses and pore pressures. The 

vertical scale is depth in seabed normalized to sea-surface wave wavelengths. The plot colors 

represent different methods of calculation (see text body) for the shallow seabed pore 

pressures, shear stress, and overburden. The line patterns represent the different output 

parameters: panel a, solid – shear stress; panel b, solid – pore pressure; panel b, dash-dot -

earth pressure; panel b, dotted – developed residual pore pressure. In panel b, when the solid 

a b c
30

pore pressure exceeds earth pressure, liquefaction is indicated (Jeng et al., 2006). Panel c is 

for example only, showing model-run diagnostics.

suspendiSlip. Bottom-concentrated suspended sediments which are liable to initiate 

gravity flows from the shelf down onto the slope originate from two major sources: (i) riverine 

sediment discharges, and (ii) resuspension of the seafloor by currents and waves. 

Sediments of the region vary significantly from gravel through sands to muds and clays. 

Therefore, the ROMS data on waves and currents were applied separately to the spatially varying 

substrates from dbSEABED to calculate pickup. This was done in the module suspendiSlip 

using models of various applicability:

QQ = QQ exp�ββ�ττ − ττ �, Parchure and Mehta (1985) for muds, and

CC0 = 0.005ρρ θθrr

3, Nielsen (1986, 1992) for sands.

(11)

where QQ is sediment erosion rate (kg/m2s), QQ is the erosion without flow velocity, and ββ (~10 

m/N0.5) is a factor involving material and inversely dependent on temperature. The bottom shear 

stresses are ττ and ττ (ττ ), applied and threshold. 

The muds erosion rate is linked to concentrations through an expression for balance of the 

upwards turbulent entrainment and downwards settling out of sediment. The formulation of 

concentration at height zz above bed depends on erosion rate (QQ ) and settling velocity (ws, m/s): 

CCzz = QQ

ww � zzzzaa

�

−

(12)

where zz is the reference height above the bed (m, ~0.05), pp is the Rouse exponent pp = wwss

κκuu∗

, 

described using von Karman constant κκ (~0.41), the flow friction velocity uu∗ and particle settling 

velocity (ww ). An expression for settling velocity from Ferguson and Church (2004) was 

employed: 

ww = RRRRDD50

CC1νν + �(0.75CC2RRRRDD50

3 ) � (13)

.with νν the kinematic viscosity (1.83x10−6, m2/s) of seawater, RR is the submerged specific density 

(~1650, kg/m3

) of grains.

For the sands erosion rate relation, CC0 is the time-averaged sediment concentration 

(kg/m3

) at a reference level, ρρ is the sand particle density (kg/m3), and θθrr is the Shields 

Criterion value for rough (e.g., rippled) bed. The Shields Criterion (SC) has several forms: 

θθ2.5 = ff2.5UU0

2

2(ss − 1) DD50

� for grain roughness SC,

θθrr = θθ2.5

(1 − ππππ/λλ) � 2 for rough bed SC,

(14)

having the sediment median grain size (m) DD50 , gravitational acceleration (m/s2), bottom flow 

velocity UU0 (m/s), roughness (ripple) height , and wavelength λλ (m). The grain roughness wave 

friction factor is assessed as 
31

ff2.5 = 0.237 �AA0

5DD50 � �

−0.52

, (15)

where AA0 is the wave bottom excursion length.

The suspended concentration depends on a reference height and concentration: 

CCzz = CC0 �

zzaa(h−zz)

zz(h−zzaa)

�

− (16)

where zz is the reference height above the bed (m, ~0.05), h is the water depth (m) and pp is the 

Rouse exponent pp = wwss

κκuu∗

, described using von Kármán constant κκ (~0.41) and the flow friction 

velocity uu∗ (m/s).

The particle suspension criteria was pro-rated for the proportion of gravel and rock in areas, 

since these restrict the amount of sediment available. Pro-rating in this way will be extremely 

important in the Viosca Knolls – de Soto Canyon sectors where rock areas are widespread (Fig 

4.6; see also Gardner et al., 2007).

To model the distribution of turbidity flows in time and space over the area and through weather 

events, we adopted the methods of Parker (1982), Parker et al. (1986) and Pantin and Franklin 

(2009). They set out criteria for ignition of flows, that is, for transformation from sediment laden 

gravity flows to self-sustaining flows undergoing an autosuspension process that makes them 

capable of travelling for long distances at high speeds. In particular, we implemented the 

methods of Pantin and Franklin (2009, p.870), which employ a flow Richardson Number for 

detecting ignition. In a water column with suspended sediment, the increasing dilution of 

sediment upwards means that there will be an optimal thickness of turbid water for ignition of a 

flow. Accordingly, the basal most-turbid water layer is checked for the ignition (autosuspension) 

criterion, and the following are provided as inputs to the RANS/TURBINS model: flow height, 

suspended sediment concentration, and peak downslope velocity.
32

Figure 4.6: Maps of the regional distributions of (a) gravel and (b) rock - facies areas from 

which little suspended sediment will be available (data from dbSEABED). The deeper colors 

show where the greatest amount of gravel and rock lie, gray is zero, white unknown. Gravel 

includes shells and is likely to armor the seabed against erosion and resuspension.
33

The suspendiSlip basis at each modeling site is a profile of the suspended sediment 

concentrations (CC) and maximum oscillatory bottom velocity from waves (UU0) (Fig. 4.7). A first 

criterion is the densiometric Richardson number RR , calculated as: 

RR = RR h

UU0

2 � CCzz

h

zz=0

(17)

for each layer, integrated from bottom. The minimum RR and corresponding h are taken to 

represent the stability of the sediment-laden bottom waters at the site. For RR > 1 the layer is 

subcritical and deposition is likely to occur (Pantin and Franklin 2009). For RR < 1 the flow is 

supercritical and motion with mixing is likely. Implementation of the Pantin and Franklin (2009) 

criterion is complicated, involving numerical estimation methods, and description is left to their 

paper. The treatment in Parker et al. (1986) of turbidity current physics is not quite relevant to 

the case here, which is of a perched sediment-laden flow with only the potential to be gravityentrained. The gravity flow stage (if sustained) is for treatment in RANS/TURBINS.

The ignition results from suspendiSlip are passed to RANS/TURBINS via a data file. They 

included information in each case on the layer-averaged sediment concentrations (C), maximum 

bottom velocity (U0), thickness (h), the sediment texture (D50), and turbulent intensity (K; m2

/s2). 

The latter was scaled to friction velocity (u*) using the assumption of Parker et al. (1986, eq. 

26): = uu∗

2⁄αα where αα = ~0.1.

These tests in hurriSlip for flow ignition events are connected with the outputs from 

ROMS/CSTMS on the prevailing wave conditions. Taking a broad and early view of the results 

(Fig. 4.8) it appears that situations prone to gravity failure and sediment flow ignition are 

modeled as occurring quite widely over time on the Louisiana-Mississippi-Alabama shelf during 

hurricane events, depending on the paths of the storms. The plot is for one period of a month 

(September 2008), including two extreme events. This accords with the results of several studies 

in the area under such conditions (e.g., Nodine et al., 2006, 2007), and also with evidence of 

cross-shelf gravity flows during storms and flood discharges in other regions (Traykovski, et al., 

2000). Overall, about 640 separate sites have at least one, at least short-term, sediment flowignition during the wave records. The total number of at least momentary ignition events on a 3-

hourly time-stepping through one month was 20,000, the vast majority of which take place at the 

storm/hurricane climax.

The deepest water flow ignition site, based on the modeling, is at 135m water depth on the NE 

flank of the Mississippi Canyon, witnessed silt grain size and suspended sediment concentrations

of 5x10-5 v/v.

One caveat to this modeling exercise is that there is no allowance for the regional 

hydrodynamic flows. A protocol fixing for the reference level for concentrations is also needed, 

grain size differentiation in the suspendates is ignored, and wave-induced cyclic velocities 

should be better accounted for. Nevertheless, we provide a mapping of the modeled distribution 

of suspended sediment flow sites that may be prone to down-slope ignition of turbidity flows. 

They do not necessarily correspond with the areas of greatest sediment resuspension.
34

Figure 4.7: Example of the suspendiSlip run-monitoring display for one of the seabed locations 

where sediment flow ignition is predicted by the model. The upper panel shows the suspended 

sediment profile through the time period, with times of greater suspension (red). The middle 

panel shows variations of (i) layer densiometric Richardson Numbers (‘Ri’) and (ii) Bagnold 

Criterea (‘KBtst’), representing the balance of potential and kinetic energy in the flows. Turbid 

water masses where Ri<1 and KBtst>1 are capable of ignition. (The storm event here is 

marginal for ignition.) For each criterion, two methods of calculation were used (‘pfzRiA’, 

‘pfpRiA’ and ‘KBtstA’, ‘pfpKBtstA’), as described in the text. The lower panel shows the 

progress of the wave climate at the cell in terms of wave height (Hs, blue), wave period (Tp, 

red), and the water column suspended sediment (SSC ppm v/v, green) from seabed erosion. 

Some parameters were log-transformed for the purposes of display.
35

Figure 4.8: Modeled occurrences of near-bottom wave-supported sediment-laden waters that 

could have progressed to turbidity currents, in the period of September 2008. The criterion is a 

small (<1.0) Richardson Number. Red to gray signifies, Richardson number increasing from ~0 

to 2.

We expect the pattern to be strongly dependent on the path of a storm– among other factors. 

Depth contours are marked in meters. Note the generation of sediment laden flows on the 

shelf; only those near the shelf edge are germane to this project.

4.4 LES/RANS-TURBINS

Discretization

TURBINS solves the incompressible Navier-Stokes equations in the Boussinesq 

approximation, along with a convection-diffusion equation for the sediment concentration 

(Necker et al. 2002, 2005). The turbidity current is modeled as an incompressible fluid laden 

with small, polydisperse particles whose density is significantly higher than the fluid density. 

The dispersed phase is assumed to be sufficiently dilute so that the two-way coupling between 

the fluid and the particles will be primarily governed by the momentum equation, and the effect 

of the particle loading can be neglected in the continuity equation. This furthermore allows us to 

neglect interactions among the particles, such as hindered settling. In addition, the particles are 

assumed to have an aerodynamic response time much smaller than typical fluid flow time scales 

(Maxey and Riley, 1983). Hence, the particle velocity is given by the sum of the fluid velocity 

and the constant settling velocity. Particle field is treated in a Eulerian manner by solving a 

transport equation. Note that in this approach, individual particles are not tracked in a 

Lagrangian manner. Polydisperse distributions are implemented in a straightforward fashion by 
36

considering different particle families, each with its own settling velocity and contributing to the 

overall density distribution.

The governing equations for mass, momentum and sediment concentration are:

∂∂uujj

∂∂ jj

= 0 (18)

∂∂uu

∂∂∂∂ +

∂∂�uu uujj�

∂∂ jj

= − ∂∂

∂∂ +

1

RR

∂∂2uu

∂∂ jj∂∂ jj

+ cceegg (19)

∂∂

∂∂∂∂ +

∂∂�cc(uujj + uu eejj

gg)�

∂∂ jj

= 1

RR

∂∂2cc

∂∂ jj∂∂ jj

(20)

Here uu , pp, and cc refer to the velocity, pressure and sediment concentration, respectively, while

uu is the settling velocity of the particle, and the governing equations are in non-dimensional 

form. The Reynolds number (RR ) and Schmidt number (SS ) that appear in the above equations 

are defined as:

RR = UU

νν

(21)

SS = νν

αα

(22)

where UU is a characteristic velocity scale, LL is a characteristic length scale, νν is the kinematic 

viscosity and αα is the diffusivity coefficient of the sediment. 

The original TURBINS code (Nasr-Azadani and Meiburg 2011) used third order Total 

Variation Diminishing Runge–Kutta method for temporal discretization, and ENO for spatial 

discretization of the convective terms. This discretization yields accurate results for low 

Reynolds number laminar flows. However, this particular choice of discretization was found to 

be excessively dissipative for high Reynolds number turbulent flows. To obtain more accurate 

results, we implemented new non-dissipative discretization in the code. Specifically, we 

implemented energy-conserving second order scheme proposed by Morinishi et al. (1998). For 

the temporal discretization, we implemented a three-substep third-order accurate Runge-Kutta 

scheme by Spalart et al. (1991). This method determines the velocity field uu +1 at time ( + Δ ) 

from the known state uu at time by using three intermediate velocity fields uukk

uukk − uukk−1

Δ = ββkkDD(uukk−1) + ββkkDD(uukk) + γγkkCC(uukk−1) + γγkkCC(uukk−2) − 2ββkk∇∇ppkk (23)

for kk = 1, 2, 3. Here, uu0 (uukk−1 for kk = 1) is the starting point of the integration step (uu0 = uu ) 

and uu3 denotes the solution for the velocity field after one entire time step (uu3 = uu +1). Earlier 
37

stages of the velocity field (such as uukk−2 for kk = 1) are not required in this scheme, since ζζ1 is 

set to zero. The other eight coefficients are defined as follows:

ββ1 = 4

15 ββ2 = 115 ββ3 = 16

γγ1 = 8

15 γγ2 = 512 γγ3 = 34

ζζ2 = −17

60 ζζ3 = −512

(24)

In the temporal discretization, the operators D and C represent the spatial discretization of 

the diffusive and convective terms of the governing equations. The diffusive terms have been 

treated implicitly using a Crank-Nicolson scheme, which can be seen by the existence of the 

second term of the right side of the temporal discretization, whereas a second-order AdamsBashforth scheme is used for the convective terms. 

Due to the coupling between velocity and pressure field, the momentum equations are 

solved in time using a projection method which calculates an intermediate velocity field uu∗ based 

on the old pressure field ppkk−1 obtained by

uu∗ − uukk−1

Δ = RR − 2ββkk∇∇ppkk−1

(25)

with RR being the residual term that includes the convective and the diffusive terms. To account 

for the influence of the pressure gradient on the new velocity field uukk, the pressure correction is 

applied in a second step using 

uukk − uu∗

Δ = −2ββkk∇∇(ppkk − ppkk−1) (26)

The new pressure field ppkk, which is still unknown, is obtained by using the assumption of a 

solenoidal velocity field at any instance in time. For example, assuming the velocity field uukk is 

divergence free yields

∇∇2(ppkk − ppkk−1) = 1

2ββkkΔΔ ∇∇uu∗ (27)

The above Poisson can now be solved to evaluate ppkk, since ppkk−1 and uu∗ are known. 

Parallelization and Linear solver development

In the original TURBINS code, we used PETSc, a tool developed by Argonne National 

Laboratory, for various tasks including data structure definitions in the code, inter-processor 

communication of data and solving linear equations. To decrease the CPU time needed to run 

the code, we developed an in-house linear solver to solve the linear equation that arises while 

discretising the Poisson equation for pressure. 

The PETSc’s linear solver uses a Krylov subspace based iterative method with an algebraic 

multi-grid method as a preconditioner. The convergence of iterative methods depends on the

condition number of the matrix, which increases with grid-cell anisotropy. For highly 
38

anisotropic grid-cells needed to resolve wall-bounded flows, iterative methods take a large 

number of iterations and a much longer CPU time to converge. Direct solvers are preferable as 

the CPU time needed to solve the equations is independent of grid-cell anisotropy. When we use 

a uniform mesh along streamwise and spanwise directions, Fast Fourier or Fast Cosine 

transforms (FFT/FCT) can be used to reduce the three-dimensional Poisson equation into a series 

of one-dimensional Helmholtz equations. Discretizing these one-dimensional equations results 

in tri-diagonal system of linear equations, which can be solved using the Thomas algorithm.

Parallelization is achieved in the code by decomposing the domain along all three 

directions. The subdomains thus obtained by the decomposition contain a small range of the 

original nodes along each of the three directions. Since the FFT/FCT need all of the data along a 

given direction, the layout of the subdomains during the intermediate stages of solving the 

Poisson equation change. First the domain is decomposed along only the wall-normal and the 

streamwise directions so that each processor has all the data along the spanwise direction of a 

given x- and y-location. FFT/FCT is taken along the spanwise direction and then the domain is

decomposed along the spanwise and wall-normal directions. At this stage FFT/FCT is taken 

along the streamwise direction. Next the subdomain layout is changed so that the domain is 

decomposed along the streamwise and spanwise directions. The tri-diagonal system that relates 

the pressure along the wall-normal direction is solved using the Thomas algorithm. Next the 

same steps are performed but in reverse order. The domain is decomposed along spanwise and 

wall-normal directions. Inverse FFT/FCT is taken along the streamwise direction at this stage. 

Then, the layout is changed again so that domain is decomposed along streamwise and wallnormal directions. Now the pressure is obtained after taking inverse FFT along the spanwise 

direction. Finally, the original layout in which the domain is decomposed along all three 

directions is restored. Note that at each of the intermediate stages, when the layout of the 

domain decomposition is changed, the pressure or its derivatives are transferred from one 

processor to all the other processors using MPI commands.

Next we compare the performance of this Poisson solver with the one from PETSc’s linear 

solver. We simulated a periodic channel flow for the purpose of benchmarking our solver. We 

performed all the simulations with a grid that has 128 nodes along the streamwise direction, 128 

nodes along the wall-normal direction and 128 nodes along the spanwise direction. We always 

used a uniform mesh along streamwise and spanwise directions in all our test runs. We used 

both uniform and non-uniform meshes along the wall-normal direction. In actual calculations, 

we always have to use a non-uniform mesh to efficiently resolve the gradients along the wallnormal direction. We tested our solver by running the same mesh but with 4, 8 or 16 CPUs. We 

use the following quantities, SS and EE , as diagnostics for the performance of the solver. SS

compares the scalability of the solver. For a perfectly scalable solver, SS should be 2 with 8 

CPUs and 4 with 16 CPUs. EE measures the factor of computational savings achieved with the 

newly developed in-house solver as compared to PETSc’s solver. 

SS = CC tt ww h NNCCCC = 4

CC tt ww h hee ss ss ss ss

(28)

EE = CC tt ww h PP PP ss ss ss

CC tt ww h h ss ss ss aa ss NNCCCC

(29)
39

Table 4.1 summarizes the performance of both solvers with uniform mesh. As compared to 

PETSc’s solver, our in-house solver runs approximately 80-120 faster for uniform meshes. The 

PETSc’s solver scales slower than the desired linear speed-up with the increase in CPUs. A 

surprising fact is that our in-house solver has super-linear scalability, i.e., EE is higher than 2 

with 8 CPUs and is higher than 4 with 16 CPUs. Further analysis of the CPU time taken at 

various stages of our in-house solver reveals that the observed super-linear scalability is due to 

faster access of elements from the memory. For instance, the data the along the slowest varying 

index (k-index in kji format in programming language C) needs to be accessed repeatedly. This 

entails jumping (NY*NX) elements in the memory where NY and NX are the number of gridnodes used in the sub-domain by a given processor. When we increase the number of CPUs, the 

subdomain becomes smaller thereby decreasing NX and NY values. This results in smaller 

jumps in memory with a decrease in overall CPU time and a much better speed-up than the 

desired linear speed up. The grid spacing given below is non-dimensionalized by the height of 

the channel.

ΔΔ = ΔΔΔΔ = ΔΔ = 0.0234375

Solver NCPU Number of points 

per subdomain

CPU time 

in seconds

SS EE

PETSc 4 524288 15.99 1 N/A

PETSc 8 262144 8.18 1.76 N/A

PETSc 16 131072 4.64 3.45 N/A

In-house 4 524288 0.19 1 82.13

In-house 8 262144 0.086 2.25 94.55

In-house 16 131072 0.039 5.02 119.55

Table 4.1: CPU time in seconds taken in various runs by PETSc and in-house solver for 

uniform mesh.

Table 4.2 summarizes the performance of both solvers for non-uniform mesh along the 

wall-normal direction. Maximum aspect ratio of the grid-cell for the current mesh is 75. Our inhouse solver takes approximately the same CPU time as in the previous uniform mesh case. The 

small difference in CPU time that we observe here is perhaps due to other factors such as the 

difference in amount of communications that is taking place in the cluster due to other users 

programs etc. As expected, the iterative PETSc solver has considerably slowed down. Now our 

in-house code runs approximately 540-710 times faster than PETSc’s solver. Note that at higher 

Reynolds number, grid-cell aspect ratio of up to 100 needs to be used for efficient resolution of 

wall-normal gradients. In these cases, our in-house code would be even more efficient as 

compared to PETSc’s solver, as it would be faster by much more than the 710 times we observed 

here with aspect ratio of 75.

ΔΔ = 0.05 ΔΔ = 0.025

MM AA RR cc = ΔΔ

ΔΔΔΔmm mm

= 75

ΔΔ mm mm = 6.6668 × 10−4 ΔΔ mm mm = 0.037732
40

Solver NCPU Number of points 

per subdomain

CPU time 

in seconds

SS EE

PETSc 4 524288 102.69 1 N/A

PETSc 8 262144 54.14 1.95 N/A

PETSc 16 131072 27.81 3.69 N/A

In-house 4 524288 0.19 1 540.02

In-house 8 262144 0.086 2.21 629.67

In-house 16 131072 0.039 4.88 713.12

Table 4.2: CPU time in seconds taken in various runs by PETSc and in-house solver for nonuniform mesh.

In our original discretization of the momentum and scalar transport equations, we used 

semi-implicit treatment for all the diffusive terms. Semi-implicit treatment is preferred as this 

removes the time step restriction imposed by the viscous stability condition, which dictates that 

time step be proportional to the square of grid spacing when using explicit treatment. In wallbounded flows, grid spacing used along horizontal directions is often one or two orders of 

magnitude larger than the smallest grid spacing near the wall along the wall-normal direction. 

This implies that, when using explicit treatment, the time step dictated by horizontal diffusive 

terms is two to four-orders of magnitude larger than the time step dictated by the wall-normal 

diffusive term. In our current discretization, we are using semi-implicit treatment for the wallnormal diffusive term and explicit treatment for horizontal diffusive terms. This discretization 

results in a tri-diagonal system of equations instead of the original hepta-diagonal system of 

equations. With this discretization, the wall-normal diffusion term does not impose any limit on 

time step and the time step dictated by horizontal diffusion terms is much larger than the time 

step imposed by the convective CFL condition which effectively means that explicit treatment of 

horizontal diffusion terms do not impose any time step limit in our simulations.

Parallelization is achieved in the code by decomposing the domain along all three 

directions. Thus, the subdomains obtained by the decomposition contain a small range of the 

original nodes along each of the three directions. As the discretization of momentum and 

transport equations results in tri-diagonal equations along the wall-normal direction, the data 

along all the nodes in the wall-normal direction for a given x- and z location is needed to solve 

the tri-diagonal equations. After we have calculated the right hand side of the discretized 

equation for all the nodes in a given subdomain, the domain layout needs to be changed so that 

the domain is decomposed along x- and z-directions. Data from the original layout is transferred 

to the new layout using MPI commands. In the current layout, each processor has data on all the 

nodes along the wall-normal direction for a given x-and z location. The Thomas algorithm is 

employed to solve the tri-diagonal system of equations now. Once the velocity component or the 

concentration is obtained, the data in the current layout is transferred to the original layout using 

MPI commands.

Next we compare the performance of our newly developed linear solver with the one from 

PETSc’s linear solver. We simulated periodic channel flow for the purpose of benchmarking our 

solver. We performed all the simulations with a grid that has 128 nodes along the streamwise 

direction, 128 nodes along the wall-normal direction and 128 nodes along the spanwise direction. 

We used a uniform mesh along streamwise and spanwise directions in all our test runs. We used 

non-uniform mesh along the wall-normal direction to efficiently resolve the gradients along the 
41

wall-normal direction. We tested our solver by running the same mesh but with 8, 16 or 32 

CPUs. We use the quantities SS and EE as diagnostics for the performance of the solver. SS

compares the scalability of the solver. For a perfectly scalable solver, SS should be 2 with 16 

CPUs and 4 with 32 CPUs. EE measures the speed-up factor achieved with the newly developed 

in-house solver as compared to PETSc’s solver. 

SS = CC tt ww h NNCCCC = 8

CC tt ww h hee ss ss ss ss

(30)

EE = CC tt ww h PP PP ss ss ss

CC tt ww h h ss ss ss aa ss NNCCCC

(31)

Table 4.3 summarizes the performance of both solvers. As can be seen from the table, our inhouse code runs 45-65 times faster than PETSc’s solver. Both PETSc and our in-house solver 

show better than the expected linear speed up. With the newly developed in-house Poisson 

solver and the current linear solver for momentum and scalar transport equation, when a nonuniform mesh is used, our code (entire code for one time step advancement) runs about two 

orders of magnitude faster than compared to the previous version of the code.

ΔΔ = 0.05 ΔΔ = 0.025

MM AA RR cc = ΔΔ

ΔΔΔΔmm mm

= 75

ΔΔ mm mm = 6.6668 × 10−4 ΔΔ mm mm = 0.037732

Solver NCPU Number of points 

per subdomain

CPU time 

in seconds

SS EE

PETSc 8 262144 5.16 1 N/A

PETSc 16 131072 2.53 2.036 N/A

PETSc 32 65536 1.19 4.34 N/A

In-house 8 262144 0.117 1 44.03

In-house 16 131072 0.0409 2.86 61.86

In-house 32 65536 0.0182 6.43 65.2

Table 4.3: CPU time in seconds taken in various runs by PETSc and in-house solver.

Large-eddy Simulation (LES): model and representative results

As a first step to upscale the DNS code, we implemented a LES model in the code. More 

specifically we implemented a dynamic model that calculates the subgrid-scale stress (SGS) 

terms that appear in the governing equation, based on the resolved flow field. The SGS terms 

are

ττ – 1

3 ττkkkkδδ = −2νν SS̅ (32)

where νν is the SGS eddy viscosity, which is calculated using
42

νν = (CC ΔΔ)2|SS̅| (33)

Here SS̅ is the strain-rate tensor, |SS̅| is its magnitude and ΔΔ is the local grid spacing. The 

unknown coefficient CC is calculated using the dynamic model of Germano et al. (1991).

To evaluate the performance of the dynamic LES model for high Reynolds geophysical 

flows, LES simulations were performed of a gravity current propagating over flat terrain. We 

considered a lock-exchange configuration, in which heavy fluid (shown in red in Figure 4.9) is 

present on the left side of a partition and a light fluid (shown in blue in Figure 4.9) is present on 

the right side of the partition. As the partition between the two fluids is removed, the heavy fluid 

propagates to the right on the bottom terrain while the light fluid propagates to the left along the 

top surface. At a later time Kelvin-Helmholtz billows can be seen at the interface between the 

heavy and the light fluid, and lobe and cleft instability can be seen near the front of the current. 

Figure 4.10 shows the contours of streamwise velocity on a plane parallel to the wall. In a fully 

turbulent boundary layer, because of the presence of streamwise vortical structures, an 

alternating high and low speed fluid region is observed. The signature of this pattern is reflected 

in the final sediment deposit profiles. These high and low speed streaks can be seen in Figure 

4.10, indicating that LES was successful in accurately predicting the flow structures present in 

high Reynolds number fully turbulent flows.

Figure 4.9: Gravity current propagating over a flat terrain. The domain is 30 units long where 

the length is non-dimensionalized by the height of the domain. At this instant in time, the 

current occupies approximately 20 units.
43

Figure 4.10: Contours of streamwise velocity fluctuations. Low- and high-speed streaks 

expected in turbulent boundary layers can be seen. Here X denotes the streamwise direction

and Z denotes the spanwise direction. The streamwise velocity is non-dimensionalized by the 

buoyancy velocity.

Note that this simulation was performed at Reynolds number of 100,000, where the 

Reynolds number is based on buoyancy velocity and the height of the channel. This corresponds 

to a gravity current of height 1m and a current velocity 10 cm/sec. This LES calculation required 

270 million grid nodes and ran on several hundred CPUs over a few weeks. While the LES is 

capable of resolving fine flow details due to turbulence, it is still limited to small domain sizes 

due to its high computational resource requirements. 

Immersed Boundary Method

We implemented a new immersed boundary method that is appropriate for simulation of high 

Reynolds flows using either LES or RANS. Figure 4.11 shows the identification of various 

types of nodes in the new immersed boundary method. In this method, the first node just outside 

the solid object is treated as an immersed node (shown as filled triangles in Figure 4.11), all the 

other nodes that are outside of the solid object are treated as fluid nodes (shown as filled circles 

in Figure 4.11). Flow field values at an immersed node are obtained by interpolating the flow 

field from the nearest fluid nodes. Three different scenarios are shown in Figure 4.11. The flow 

field at node 1 and node 2 are interpolated to obtain the flow field at point 3, which is the 

intersection point of the normal from the boundary to the line joining node 1 and node 2. The 

flow field at point 3 and no-slip or no-flux condition at the boundary point B are then used to 

calculate the flow field at the immersed node. In LES and RANS simulations, we need to 

calculate the eddy viscosity at the immersed node as well. Since the eddy viscosity is zero at the 
44

wall, the eddy viscosity at the immersed node is obtained using the same interpolation procedure 

as the velocity field.

Figure 4.11: Grid layout showing fluid-solid interface treatment in the new immersed boundary 

method implementation. Shaded region: solid nodes, filled circle: fluid nodes, filled triangle: 

immersed nodes.

To validate the new immersed boundary method, we considered flow past a wavy wall. 

Direct Numerical Simulation (DNS) of this flow has been performed by Maas and Schumann 

(1996) using a curvilinear code that uses a body fitted grid. Figure 4.12 shows the outline of the 

geometry as well as the grid used in the LES calculation. The configuration consists of a channel 

whose bottom wall is wavy while the top wall is horizontal. In the simulation, periodic boundary 

conditions are used in the streamwise and spanwise directions. In the wall-normal direction, the 

no-slip condition is used at the top wall and the immersed boundary method is used to impose 

the no-slip condition at the bottom wall. The dimensions along the streamwise and the wallnormal direction are shown in Figure 4.12. The domain was 4 units wide in the spanwise the 

direction. The LES calculation used 400 X 130 X 64 grid nodes in the streamwise, wall-normal, 

and spanwise directions respectively.
45

Figure 4.12: Grid used in the calculation of flow past a wavy wall. Every other grid node is 

shown along both directions.

The grid spacing was uniform along the streamwise and the spanwise directions. A nonuniform grid was used along the wall-normal direction with grid nodes clustered near the top and 

bottom wavy wall to better resolve the boundary layer. Flow over the bottom wall experiences 

first an adverse pressure gradient on the diverging section and then experiences a favorable 

pressure gradient on the converging section. Figure 4.13 shows the mean streamwise velocity 

field and the streamlines over the first half of the domain. Due to the symmetry of the flow field 

about the centerline along the streamwise direction, the flow field over the first half of the 

domain is discussed. Streamlines on Figure 4.13 show that the flow separates on the diverging 

section and then reattaches on the converging section. The LES calculation predicts the 

separation at 0.31 and the reattachment at 1.22, which agrees quite well with the DNS prediction 

of the separation at 0.28 and the reattachment at 1.20. 
46

Figure 4.13: Mean streamwise velocity field and the streamlines.

Figures 4.14 and 4.15 show the mean streamwise and wall-normal velocity profiles at 

selected streamwise locations. The first location (x=0.2) shown is upstream of the separation 

point; here, both streamwise and wall-normal velocity show excellent agreement with the DNS 

data. The next two locations (x=0.6 and x=1.0) are inside the separation region, while the last 

location (x=1.4) is in the recovery region downstream of the reattachment point. At all these 

locations, the mean streamwise velocity shows excellent agreement with the DNS data. The 

mean wall-normal velocity has a small deviation from the DNS data in the near wall region. 

This small difference could be due to the differences in the size of the separation region 

predicted by LES and DNS.
47

Figure 4.14: Mean streamwise velocity at selected planes.

Figure 4.15: Mean wall-normal velocity at selected planes.

Figure 4.16 shows the turbulent kinetic energy profiles at the same locations discussed 

above. Note that there is an increase in turbulent kinetic energy in the separation region due to 

high values of velocity gradient present in this region. Overall there is a good agreement 

between the LES and the DNS predictions except for small differences in the recovery region 

(x=1.4).
48

Figure 4.16: Turbulent kinetic energy at selected planes.

Reynolds Averaged Navier-Stokes (RANS) model

As LES is an eddy resolving technique, it needs a computational grid with relatively fine 

spacing. This requirement limits the application of LES in terms of the Reynolds number, and to 

domains whose length along the longitudinal direction is less than 50 times the height of the 

current. For a current height of 10m, the domain length would be limited to less than 500m.

Since BOEM’s interest is in much larger domains than feasible with LES code, we 

upscaled the LES code using k-ɛ RANS model. RANS simulations can employ much coarser 

grid spacings than LES, enabling us to simulate much larger domain. For instance, for a current 

height of 10m, the RANS code can simulate a longitudinal domain of O(10KM) length. 

We implemented the non-hydrostatic k-ɛ RANS model in the LES code. As a validation of 

the RANS code, we ran channel flow driven by a constant streamwise direction at various 

Reynolds numbers. Figure 4.17 shows the mean velocity predicted by the k-ɛ RANS model. At 

all the Reynolds numbers, the mean velocity shows good agreement with the law of the wall 

indicating that k-ɛ RANS model is implemented accurately. Note that only the bottom 20-30% 

of the boundary layer (current height) is expected to obey the law of the wall. Current 

implementation of k-ɛ model uses wall functions approach so that coarser spacing can be 

employed along the wall-normal direction.
49

Figure 4.17: Velocity profile predicted by the k-ɛ RANS model at various Reynolds numbers.

Streamwise consecutive modules 

The RANS model enables us to simulate flows at any arbitrarily high Reynolds number. 

However, huge computational resources in terms of both CPUs and memory are needed to 

simulate large domains. Since, such computational resources are not available, we implemented 

inflow-outflow modules that simulate large domains by splitting them into multiple subdomains 

and transferring the flow field values at the boundaries to the adjacent subdomain. These inflowoutflow modules store the velocity, concentration, turbulent kinetic energy and turbulent 

dissipation data in an upstream domain. These data are then read by the downstream box, and 

used as the inflow conditions for the downstream box. 

To verify these modules, we first ran a simulation in a domain that is 2000m long in the 

streamwise direction and 200m tall in the vertical direction. Initially the first domain had only 

fresh water without any particles. At the beginning of the simulation, a current with a height of 

100m entered the domain. Data at 1000m from the inlet in the simulation were stored. These 

data were used as inflow conditions for a simulation performed on a second domain, which was

1000m long and 200m tall. Ideally, if these inflow-outflow modules work as expected, then the 

flow development in the second domain must be very similar the flow development in the last 

1000m of the first domain. 

Figures 4.18, 4.19 and 4.20 compare the sediment concentration in the first and the second 

simulations at t=600, 1000 and 1600 seconds. These figures show that the flow development in 

the second domain is very close to the flow observed in the last 1000m of the first domain. This 

indicates that inflow-outflow modules are working accurately. 
50

Figure 4.18: Contours of concentration at time = 600 seconds.

Figure 4.19: Contours of concentration at time = 1000 seconds.
51

Figure 4.20: Contours of concentration at time = 1600 seconds.

Erosion model

Garcia and Parker (1993) proposed an empirical formula to represent the resuspension flux 

of sediment into the current. Blanchette et al. (2005) has previously used this resuspension flux 

formula to model erosion in a low Reynolds number simulation of turbidity currents. We 

implemented a variation of this model in the non-hydrostatic RANS code. The sediment flux 

due to erosion is introduced into the current as a diffusive flux from the bottom wall.

− 1

SS

∂∂

∂∂ = uu EE (34)

where cc is the non-dimensional concentration of the sediment, is the coordinate along the 

normal direction to the boundary, uu is the settling velocity, EE is the resuspension flux, SS is the 

Schmidt number, and RR is the Reynolds number. The resuspension flux EE is evaluated using 

the following expression,

EE = 1

CC0

aaZZ5

1 + aa

0.3 ZZ5

(35)

Here aa = 1.3 × 10−7, CC0 is the initial volume fraction of the sediment and ZZ is the erosion 

parameter. The resuspension flux is also set a maximum value of 0.3/CC0. The erosion 

parameter ZZ is calculated as
52

ZZ = 0.586

uu∗

uu

RR 1.23 RR ≤ 2.36

ZZ = uu∗

uu

RR 0.6 RR > 2.36

(36)

In the above equation, uu∗ is the shear velocity at the bottom wall, and RR is the particle 

Reynolds number. The shear velocity uu∗ is calculated using the following equation.

uu

uu∗

= 1

κκ � uu∗νν � + BB (37)

Here uu 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 the kinematic viscosity, κκ = 0.41 and BB = 5. 

The particle Reynolds number is defined as:

RR =

√� �ρρ − ρρ0�

ρρ0 � �

νν

(38)

where is the diameter of particles, ρρ is the density of the sediment, ρρ0 is the density of the 

water, and is the acceleration due to gravity. 

We implemented this erosion model in our non-hydrostatic RANS code. We ran 

simulations in lock-exchange configurations to evaluate the performance of the erosion model. 

We considered a domain that is 4000m long and 200m tall. Initially the lock region extended 

over 400m along the streamwise direction, and 100m in the vertical direction. 

We ran two types of simulations in this domain. In the first, we assumed that particles 

deposited on the bottom wall without any erosion and resuspension. In the second simulation, 

we included erosional effects through Garcia and Parker’s model. Figures 4.21, 4.22 and 4.23

show the temporal evolution of the current for the case where the erosional effect is included. 

Initially the maximum value of the non-dimensional concentration is 1. However, at t=600 

seconds, the maximum value of the concentration reaches 1.03 near the bottom wall due to 

erosion and resuspension of the sediment from the bottom wall. Figure 4.24 compares the 

fraction of the initial mass of the particles that is still suspended in the current for the two cases. 

Note that the current in both cases is depositional in the net, as the suspended particle mass 

decreases with time. The case that included erosional effects has more suspended particles at 

any instant in time indicating that erosional model has resuspended some particles near the 

bottom wall. 
53

Figure 4.21: Contours of concentration at time = 0 seconds.

Figure 4.22: Contours of concentration at time = 400 seconds.

Figure 4.23: Contours of concentration at time = 600 seconds.
54

Figure 4.24: Temporal history of fraction of initial mass still suspended in the current.

5. Datasets

Several datasets were acquired and/or processed to configure, initialize, or force the 

hydrodynamic, sediment transport, and turbidity current models including atmospheric forcing, 

surface wind-induced wave data, satellite-derived sea surface temperature, river discharge, bed 

sediment properties, and high-resolution bathymetry for sediment flow pathway analysis.

5.1 Atmospheric Forcing Data

Surface atmospheric fluxes are required to force the ocean circulation models (ROMS and 

ROMS-CSTMS). Generally, ROMS requires surface vertical boundary conditions for momentum 

(surface wind stress), temperature (surface net heat flux), and salinity (surface net freshwater 

flux, E-P). Usually, the surface fluxes for passive tracers (like suspended sediment classes) are 

set to zero. There is a choice to process these surface vertical boundary conditions as direct 

kinematic fluxes (reading from an input NetCDF file) or compute them using an air-sea 

interaction boundary layer model. All of the bottom vertical fluxes for momentum and active 

and passive tracers are also set to zero, except for the case of sediment erosion and deposition.

In ROMS, the air-sea boundary layer model is based on the bulk parameterization of Fairall 

et al. (2003), which was adapted from the Coupled Ocean-Atmosphere Response Experiment 

(COARE) algorithms for the computation of surface fluxes of momentum, sensible heat, and 

latent heat. It uses similarity theory to compute the Monin-Obukhov stability parameters for 
55

turbulent fluxes of wind, heat, and moisture. If the bulk parameterization is activated in ROMS, 

different atmospheric forcing fields are required including: near surface (2m or 10m) zonal and 

meridional wind components, air temperature, solar shortwave radiation flux, downwelling 

longwave radiation flux, air relative humidity, air pressure, and rain fall. Otherwise, ROMS just 

needs to read and process the surface wind stress, net surface heat flux, surface shortwave 

radiation flux, and net freshwater flux. There are options to read these gridded atmospheric 

fields from input forcing NetCDF files or process from a coupled atmosphere-ocean model using 

the Model Coupling Toolkit (MCT) or the Earth System Modeling Framework (ESMF) libraries. 

In this project, the gridded atmospheric data is read and processed from input NetCDF files. 

Then, the atmospheric forcing fields are time-interpolated data snapshots at every ocean model 

time-step.

At the beginning of the project, we used the atmospheric fields from the CORE2 (Common 

Ocean Research Experiment) 1° x 1° global air-sea dataset. The atmospheric fields from this 

dataset are from multiple sources but the primary data is from NCEP 6-hourly reanalysis. 

However, this dataset was too coarse for our applications and we used higher spatial and 

temporal resolution atmospheric forcing from the ECMWF (European Centre for Medium-Range 

Weather Forecasts) ERA-Interim dataset. The ERA-Interim is a global reanalysis dataset 

available at 1/4° x 1/4° resolution every three-hours for surface file. It is archived from 1979 to 

the present. We processed the data for the full Gulf of Mexico from 1-Jan-2000 to 31-Aug-2013. 

The reanalysis dataset was generated using 4-Dimensional Variational (4D-Var) data 

assimilation based on a 2006 release for ECMWF IFS system (Dee et al., 2011). 

5.2 Surface Wind-Induced Wave Data

Wave energy plays an important role in the mobilization of seabed sediment in the northern 

Gulf of Mexico (see Section 1.4). During extreme events, large waves from hurricanes can 

resuspend seafloor sediment and / or liquefy the seafloor providing a viable mechanism for 

inducing sediment gravity flows, like turbidity currents. Therefore, accurate spatial and temporal 

wave data are required in the Bottom Boundary Layer (BBL) models used to parameterize 

bottom stress due to wave-current interactions, which affect seafloor sediment transport.

Several wind-induced wave fields are required in ROMS to compute bottom drag and 

bottom stresses from the various BBL sub-models available (Styles and Glenn, 2000; Blass et 

al., 2007; Warner et al., 2008). These fields include significant wave height, wave direction, 

surface wave period, bottom wave period, bottom orbital velocity, and wave energy dissipation 

rate. These wave properties are either read from an input NetCDF forcing file or computed from 

a two-way coupled wave model (like SWAN, Simulating Waves Nearshore) with ROMS (Warner 

et al., 2008). In this project, the required wave fields are processed from the NOAA/NCEP

WaveWatch III (WW3; Tolman and Chalikov, 1994) dataset available every three hours on a grid 

having 10-arc minute resolution, as driven by GFDL-GFS winds. The WW3 dataset was 

processed from 01-Jan-2006 to 31-Dec-2012. The web-site 

http://polar.ncep.noaa.gov/waves/wavewatch/wavewatch.shtml provides more detailed

information, including how to access the wave data. Figure 5.1 shows a sample of wave height 

and period during Hurricane Gustav, which impacted the study area during the late summer, 

2008.
56

Figure 5.1: Estimates of (a) wave height (m), and (b) wave period (s) from NOAA / NCEP 

Wavewatch III model for September 1, 2008 as Hurricane Gustav crossed the study area. 

Black box indicates approximate location of the ROMS hydrodynamic Grid-f.

While WW3 provides estimates of surface wave characteristics, the bottom boundary layer 

and sediment transport formulations require near-bed wave orbital velocity and near-bed wave 

period. These were estimated from the surface wave characteristics following Wiberg and 

Sherwood (2008). Figure 5.2 shows a sample of bottom wave period and bottom orbital velocity 

during Hurricane Gustav, which impacted the study area during late summer, 2008.

Figure 5.2: Estimates of (a) bottom wave period and (b) bottom orbital velocity from 

NOAA/NCEP WaveWatch III model following Wiberg and Sherwood (2008) and interpolated to 

the ROMS hydrodynamic Grid-f on September 1, 2008 as Hurricane Gustav crosses the study 

area.

5.3 Sea Surface Temperature Data

The sea surface temperature (SST) data is used in ROMS to improve the thermal surface 

boundary condition in the temperature vertical diffusion term. The advection-diffusion for the 

time evolution of the ocean potential temperature (T) is:

∂∂

∂∂∂∂ = Advection + Diffusion + F ss(z) (39)

(a) (b)

(a) (b)
57

where F ss(z) is a source term representing the depth-dependent heating by the solar radiation 

flux. In ROMS this term is parameterized by an absorption law (Paulson and Simpson, 1977) 

where the depth-dependent shortwave radiation flux decays exponentially according to the Jerlov 

water type. The surface vertical boundary condition for the diffusion term is:

� (zz)

∂∂

∂∂ �

zz=0

= QQ

ρρooCC

(40)

where (zz) is the vertical temperature diffusion coefficient that depends on depth, QQ is the 

surface net heat flux, ρρoo is a constant seawater density, and CC is the specific heat constant for 

seawater (3985 J Kg-1 Celsius-1). Following Barnier et al. (1995), the net heat flux on surface 

temperature (Tsur) can be described as the sum of two components

QQ ( ss ) = QQnet − �

∂∂QQ

∂∂ �

sst

(SSSS − ss ) (41)

where QQnet is the net heat flux from the atmospheric model and �

∂∂QQnnnnnn

aa

∂∂∂∂ �

sst

is the net heat flux 

sensitivity to SST data. Therefore, the second term in the above equation is a net heat flux 

correction proportional to the SST data and model surface temperature.

In this project, the SST was satellite-derived, monthly-averaged data obtained from NOAA 

CoastWatch blended product for the period of January 2005 to December 2013. These data are

available from their OpenDAP catalog oceanwatch.pfeg.noaa.gov/thredds/catalog.html and are

titled SST, Blended, Global, EXPERIMENTAL (/thredds/dodsC/satellite/BA/ssta/mday). The 

net heat flux correction eliminated the warming trends in the upper 50m temperature solutions 

during the summer. Figure 5.3 shows the simulated montly-averged temperature for June 2007 

(upper panel) and January 2008 (bottom panel). The summer and winter temperatures are within 

the range observed in the SST blended maps.
58

Figure 5.3: Monthly averaged 10m potential temperature (Celsius) solutions for Grid-g for 

summer (top) and winter (bottom).

5.4 River Discharge Simulation Data

Few streams and rivers flowing into the Gulf of Mexico (GOM) are adequately gauged. 

Table 5.1 outlines the available data. Except for the Mississippi and Pearl Rivers, sediment flux 

measurements are not available for GOM Rivers. To fill in this data gap a global Water Balance 

Model (WBM) was combined with the BQART and PSI models and simulations run to provide 

estimates of daily discharge and sediment-flux for the 20th century and first decade of the 21st 
59

century (Cohen et al., 2013) based on an extension of the WBMplus global hydrology model 

(Wisser et al., 2010), which is part of the FrAMES biogeochemical modeling framework 

(Wollheim et al., 2008). A surface water balance of non-irrigated areas forms the core of 

WBMplus with flow routing employing a cell-tree topology (e.g. HydroSHEDS) and using a 

semi-implicit finite difference solution to the Muskingum-Cunge equation via a diffusive wave 

solution to the St. Venant equations that neglect acceleration terms in the momentum equation. 

Bankfull discharge is monitored for each global pixel using the Yamazaki approach (Cohen et 

al., 2014). Irrigation is a function of irrigated land within each computational grid cell, with 

irrigational water demand obtained from four sources: (1) small reservoirs, (2) shallow

groundwater, (3) nearby rivers, and (4) unsustained deep aquifers.

The annual discharge for the Mississippi & Atchafalaya Rivers (USACE data 1995-2011) 

is about 20,874 m3/s, with a standard deviation of 11,211 m3/s. The WBMsed model predictions 

for the combined Mississippi and Atchafalaya Rivers (17 year simulation, 1995-2011) is 18,300 

m3/s, with a standard deviation of 12,400 m3/s. Figure 5.4 provides a day-to-day comparison of 

the model values to the measured USACE values for this 17-y period. Much of the scatter 

relates to timing issues. The model tends to introduce a lag in the water levels.

Figure 5.4: Comparison of the WBMsed modeled water discharge (Y-axis) with the USACE 

measured water discharge of the Mississippi River (before the diversion of the Mississippi and 

Atchafalaya Rivers) for 1995-2011.
60

Data

Source Station ID Data Range

Daily

Discharge

Daily

Load

Water 

Temp.

USGS Atchafalaya River

Morgan City, LA 7381600 1988- 2013 Yes

USGS Biloxi River

Wortham, MS 2481000 1952-2013 Yes

USGS

USAGE

Mississippi River

Tarbert Landing, MS 7295100 1975-2013 Yes Yes

USGS Mississippi River

Baton Rouge, LA 7374000 2004-2013 Yes Yes

USGS Pascagoula River

Graham Ferry, MS 2479310 1968-1989 Yes Yes

USGS Pearl River

Bogalusa, LA 2489500 1968-1989 Yes Yes

USGS Tombigbee River

Coffeeville, AL 2469761 1961-2013 Yes

USGS Wax Lake Outlet

Calumet, LA 7381590 1986-2013 Yes

Table 5.1: Northern Gulf of Mexico stream and river gauging sites.

Sediment flux predictions are based on the BQART and PSI models (Morehead et al., 

2003; Syvitski and Milliman, 2007) and take into account relief, lithology, basin area, discharge, 

local temperature, reservoir trapping, and land use practices. Model implementation involves 

matching input resolution (time and space) with model resolution for global air temperature and 

precipitation (e.g. GPCC/NCEP), and appropriate boundary conditions (e.g. soil parameters, crop 

land, vegetation, small reservoir fraction, irrigation parameters, ice cover, population). Figure 

5.5 provides a model verification example of the sediment load model skill where WBMsed

model results are compared with USGS observations. Other examples can be viewed in Cohen 

et al., 2013. Figure 5.6 provides an example of how WBMsed predictions of daily sediment load 

vary with WBMsed predictions of daily discharge.

Figure 5.5: Comparison of WBMsed modeled sediment loads of the Mississippi River (before 

the diversion of the Mississippi and Atchafalaya Rivers). The WBMsed model was used to 

provide the freshwater and sediment load boundary conditions for the project’s marine models.

10

100

1,000

10,000

10 100 1,000 10,000

WBMsed Sediment Load KT/d

USGS Observed Sediment Load KT/d
61

Figure 5.6: WBMsed-derived sediment load values compared with WBMsed-derived water 

discharge predictions (1995-2011) for the Mississippi River. Daily values are connected to 

show the level of hysteresis in the model between discharge and load.

The original Mississippi flow is diverted to the Atchafalaya River by 30% ±1% (USACE data 

1995-2011), with a one-day minimum of 25% and a one-day maximum of 45%. The WBMsed

model does not separate the Atchafalaya from the Mississippi discharge, but instead delivers the 

total along the Mississippi main stem.

The WBMsed model was used to simulate 15 (Fig. 5.7) Gulf of Mexico Rivers where the 

annual discharge of a river exceeded 30 m3/s, a threshold that acknowledges the limits of the 

model’s resolution. The WBMsed model was run at a daily time scale using observed 

climatology, and simulations were provided for the 1995-2011 period. Daily discharge 

predictions compare favorably to both ground-based gauging stations and satellite-based 

observations (Brakenridge et al., 2012; Fig. 5.4). Daily load predictions also compare favorably 

against daily observational data (Cohen et al., 2013; Fig. 5.5).

10

100

1000

10000

100000

1000000

1000 10000 100000

WBMsed Sediment Load KT/d

WBMsed Water Discharge (m3/s)
62

Figure 5.7: Simulation river outlets for the WBMsed model: (1) Apalachee River, (2) 

Apalachicola River, (3) Conecuh River, (4) Choctawhatchee River, (5) Escambia River, (6) 

Mobile River, (7) Pascagoula River, (8) Mississippi and Atchafalaya Rivers, (9) Pearl River, (11) 

Grand River, (12) Sulphur River, (13) Sabine River, (14) Neches River, (15) Trinity River, and 

(16) Conroe River.

Freshwater discharge to the Gulf of Mexico varies seasonally as seen in the sum of the 15 

rivers with flood peaks in February through March (Fig. 5.8). The Mississippi discharge ranges 

from a 17-y one day low of 2,100 m3

/s (WBMsed model) or 5,667 m3/s (USACE observation) to 

a one-day high of nearly 88,000 m3

/s (WBMsed model - spillway releases to avoid infrastructure 

losses downstream would reduce this theoretical upper limit) or 65,280 m3

/s (USACE 

observation). The total discharge for all Northern Gulf rivers as modeled by WBMsed has a 17-y 

average of 22,800 m3

/s, a standard deviation of 15,400 m3/s, a one day low of 2,330 m3/s and a 

17-y one day high of nearly 96,800 m3/s. In general the non-Mississippi merged discharges 

averaged 4,540 m3/s or 19% of the total flow into the northern Gulf of Mexico, with a standard 

deviation of 4,730 m3/s (Fig. 5.9). The 17-y one-day high of these non-Mississippi sources was 

nearly 30,000 m3/s (March 8, 1998). On Nov 13, 1997 these non-Mississippi sources accounted 

for 66% of the total freshwater input into the Gulf (of 15,500 m3/s), a 17-y one-day maximum 

contribution.
63

Figure 5.8: Total daily flux of water discharge (m3/s) entering our Gulf of Mexico study site. 

The 2011 flood season (partly shown) was amongst the most devastating floods in the 

continental US history.

Figure 5.9: Contribution of water discharge into the northern Gulf of Mexico from nonMississippi River sources. Simulation from WBMsed.
64

The USGS observations (1995-2011) of the sediment load of the Mississippi River (Tarbert 

Landing values reconstituted to take into account the amount of water diverted to the 

Atchafalaya River) showed a mean of 170 Mt/y, whereas the WBMsed model for the same period 

showed a larger mean of 289 Mt/y. The model also showed a greater variability (Fig. 5.6). On 

average the Mississippi River supplies 88% of the fluvial sediment load to the Gulf, although the 

range varies from more than 99% to less than 15% on any given day. The non-Mississippi River 

loads (summation of 14 rivers) supplied 16.2 Mt/y to the Gulf. As with water discharge (Fig. 

5.8), sediment discharge to the Gulf is very seasonal but with peak flows of shorter duration (Fig. 

5.10).

Figure 5.10: Total suspended sediment flux (kg/s) entering our BOEM Gulf of Mexico study 

site. The 2011 flood season was amongst the largest in US history carried much of the 

sediment mass.

5.5 Bed Sediment (dbSEABED)

The physical properties of bed sediments are a crucial control on the movement of 

sediments in shelf-slope settings. The grain size textures determine amount of sediment pickup 

under currents and waves, and particle settling rates in wave-supported and turbulent 

suspendates. The presence and quality of clays, and sortedness of sands determine the 

occurrence of mass failures. Rock and gravel areas stabilize loose sediments and may buttress 

them against slope failure. Bedrock, if shallow, determines the extent of the deformations of 

sediments under surface ocean waves. The database dbSEABED (database – Seabed; Goff et 

al., 2008) is a source of detailed information on these factors, especially for the Gulf of Mexico 

(GoMex) where the number of samplings over the decades has reached hundreds of thousands. 

The GoMex coverage for this project in dbSEABED was queried and gridded maps were 

computed to define the spatial variations of seabed type as they might affect the results of the 

CSTMS, hurriSlip, and TURBINS models.

The new grids have 0.01 degree resolution, deal with the seabed surface (0-0.3m subbottom 

reach), and cover these parameters: moment grain size and sorting (phi), gravel, sand, mud 

0

2000

4000

6000

8000

10000

12000

14000

16000

18000

20000

10/3/94 10/2/96 10/2/98 10/1/00 10/1/02 9/30/04 9/30/06 9/29/08 9/29/10 9/28/12

Sum of All Northern Gulf River loads (Kt/d)
65

fractions (%), exposure of unerodible rock (%), and sediment carbonate as a guide to biogenics 

(%). The grids are currently in ASCII ESRI raster format.

Turning to specific features of the area which could interfere with sediment flows, 

obviously, submarine canyons and gullies, such as the strongly incised Mississippi Canyon are of 

primary interest. But, upstanding rock ridges along the brim of the shelf in the area of the De

Soto Canyon will prevent (or divert) bottom-concentrated sediment flows generated by 

waves/currents on the shelf. They will also tend to stabilize the outer shelf against mass wasting. 

An example of the ridges is the Alabama Alps feature (Weaver et al., 2001). Areas of hard 

substrate are also commonplace on the upper slope west of the Mississippi Canyon. Coarse 

sandy and gravelly (including shell and coral) areas of the shelf are unlikely to yield much 

sediment for resuspension into wave-supported sediment flows. In Figure 5.11 rock areas are 

purple, and gravel is red.

Figure 5.11: Three-dimensional view of sediment types in the project area based on data from 

dbSEABED, showing rock (purple), gravel (red), sand (yellow), and mud (green). Patchiness 

between rock/gravel and mud in areas of the slope is due to the presence of deep cold-water 

coral banks (Brooke and Schroeder 2007), low-stand shelf-edge delta remnants (Gardner et al., 

2007), methanogenic carbonate rock and rubble (Cordes et al. 2009). The shelf areas have 

widespread occurrences of gravel, and hard grounds colonized by skeletal-benthos (Thompson 

et al. 1999).

5.6 Depositional Rates Data

The seafloor database also holds information on deposition rates as observed through 

isotopic and geochemical tracers and event layers. These data are important as a comparison 

with results from the numerical modeling. However, they also have a role in understanding the 

stratigraphic development of the shelf-slope sediments, which may be prone to failure. Figure 
66

5.12 illustrates some of the depositional rate data for the project area. As is usually the case, 

rates measured on a shorter time-frame are higher (Sadler 1981). The effect is exponential. 

Determinations on longer timescales incorporate more hiatus ‘waiting-time’ events.

Figure 5.12: Compilation of deposition rate data for the Mississippi shelf-slope region based on 

radiological, stratigraphic and biostratigraphic methods. The data are presented in the context 

of smaller rates for longer intervals – the ‘Sadler Effect’ (Sadler, 1981). Note that the timescales 

increase to the right. An example of the data incorporated into this synthesis is that of Winter et 

al., (2007).

5.7 Sediment Flow Pathways Analysis

Bottom-concentrated current/wave suspended sediment flows and gravity-induced mass 

flows are guided by the local bottom topography. Hollows (‘channels’) will be the preferred 

pathways especially if those courses have been incised by earlier flows. In terms of risks to 

infrastructure posed by sediment-laden flows - the subject of this study - the pathways indicate 

the geographic areas where this risk is expected to be greatest.

With sufficiently good bathymetry the channels can often be identified. In many areas 

their patterns are clear and lead to strong inferences on flow directions and destinations. In 

other areas the slopes are plain and only weakly dissected by channels. In almost all areas, 

however, flow directions can be identified.

The SRTM 30 Plus bathymetry for the area (1/120th degree, approx. 500m resolution; 

Sandwell et al., 2014) was analyzed for downslope channeling and directions using the TauDem 

software (Tarboton 2003; Tarboton and Baker 2008), co-opted from on-land applications in river 

network analysis (Jenkins, subm.).

Provinces of different drainage patterns are identifiable throughout the project area (Fig. 

5.13). Unsurprisingly, given the commonness of cross-shelf sediment flow phenomena, shelf 

areas show a channeling network, though it is often diverted by sand ridges (e.g. south of 

Terrebone parish, LA) or rock ridges (south of Mobile, AL). Upper slope areas have either 
67

dendritic or parallel downslope drainage. The latter is most apparent in areas with relatively 

featureless slopes. The Mississippi Canyon and sub-features are clearly delineated in the 

analyzed network. The channeling patterns change as slope varies from steep to gradual. The 

computed channelized pattern was rasterized and passed on to the RANS - TURBINS modeling 

group so that those models could incorporate the varying slopes, directions, and sediments (viz. 

critical shear stresses and bed grain-sizes) along the relevant pathways.

Figure 5.13: Topographically defined sediment drainage patterns (‘channels’) for the project 

area, computed with TauDem. Grey/red streams have least incised relief and are the most 

common, green/blue streams have incision up to 15m, measured on km scales. 

6. ROMS Simulations

Several ROMS model grids were built at various resolutions in the northern Gulf of 

Mexico. The basic information for each grid is summarized in Table 6.1 and their locations are 

shown in Figure 6.1. In order to facilitate the lateral open boundary conditions with the regional 

grids (a and f), a full Gulf of Mexico grid (g) was developed. Grid-g has a mean resolution of 

10.733 km and 14.239 km in the longitude and latitude directions, respectively. This is the 

coarser grid and data donor to finer resolution grids (a and f). Grid-f covers the same 

geographical region as Grid-a, but it has finer resolution and larger number of grid points (see 

Table 6.1). Grid-a was discarded early in the project and replaced with Grid-f to facilitate twoway nesting. Grid-f has a 1:5 refinement ratio from Grid-g (see Figure 6.2)
68

Grid-f extends from 92°59′W to 86°5′W and 26°6′N to 30°41′N. It has a mean resolution 

of 2.070 km and 2.746 km in the longitude and latitude directions, respectively. This is the 

medium resolution grid and data donor to other even finer resolution refinement grids. 

Refinement grids b and c are centered in the Mississippi River delta and canyon. Grid-b has a

refinement factor of 1:5 from Grid-f. Grid-c is obtained from b with a refinement factor of 1:5. 

Refinement grids d and e are centered in the De Soto canyon. Similarly, Grid-d has a refinement 

factor of 1:5 from grid f. Grid-e is obtained from d with a refinement factor of 1:5. Notice that 

in order to achieve a grid resolution around 100m, we need to have two layers of grid refinement 

for both the Mississippi River and De Soto canyons. This exposes the computational nature of 

the project. The De Soto grids were built to examine their viability for this project as BOEM has 

some interest in the De Soto canyon area. The De Soto canyon’s (grid e) extent is wider than the 

Mississippi River canyon (grid c) and substantially more grid points are required, as shown in 

Table 6.1. In light of this, we concentrated on the Mississippi River canyon for the turbidity 

currents and sediment transport partly because there is extensive data available in the Mississippi 

River canyon. Because of the computational expense, all the ROMS simulations were carried 

out on grids g and f.

The bathymetry for all grids is a melded product between ETOPO2 (2 minute resolution) 

and 15-second resolution in the finer grid to resolve the canyons. Figure 6.3a shows the 

bathymetry, mesh, and land/sea masking for Grid-f. Notice that the bird-foot delta is resolved at 

this resolution (~2 to 3 km). The Mississippi River bird-foot delta is well resolved in Grid-b, as 

shown in Figure 6.3b. Figures 6.4 and 6.5 show the fine resolution (~100m) refinement grids c

and e and their associated bathymetry, respectively.

The initial and boundary conditions for the full Gulf of Mexico model (Grid-g) were 

processed from the North West Atlantic (NWA) ROMS solution courtesy of Enrique Curchitser 

(Rutgers University). This is a 57 years solution from 1957-2012, which uses the SODA dataset 

as lateral boundary conditions. The Grid-g was initialized on 01-Jan-2000 and spun-up for 7-

years. The NWA sea surface height solution for 01-Jan-2000 is show in Figure 6.6. The strategy 

was to run Grid-g with initial conditions and daily lateral boundary conditions from the NWA 

solution until 01-Jan-2006. The primary focus was to fine tune ROMS model options, 

parameters, and external forcing. Table 6.2 shows a summary of options and main parameters. 

Once a satisfactory configuration was achieved, various circulation experiments were initialized 

on 01-Jan-2006 and run until 31-Dec-2012.

The atmospheric forcing fields are derived from the European Centre For Medium-Range 

Weather Forecasts (ECMWF) ERA-Interim, 3-hour dataset, as discussed in Section 5.1. The 

wind-induced wave forcing used in the bottom boundary layer model (BBL) is derived from the 

NOAA WaveWatch III dataset (see Section 5.2).
69

Grid

Number Points

Min/Max ∆x

(m)

Min/Max ∆y

(m)

Mean ∆x

(m)

Mean ∆y

(m)

Refinement

Factor

a (256x128) 2595 / 2688 3108 / 3220 2642 3165 -

b (325x275) 525 / 533 630 / 639 530 634 1:5 from f

c (825x500) 105 / 107 126 / 127 106 127 1:5 from b

d (400x365) 519 / 530 622 / 635 525 628 1:5 from f

e (1250x700) 105 / 106 126 / 127 105 126 1:5 from d

f (324x150) 2032 / 2106 2696 / 2794 2070 2746 1:5 from g

g (158x98) 10160 / 11236 13484 / 14905 10733 14239 -

Table 6.1: Information about northern Gulf of Mexico ROMS grids

Figure 6.1: Gulf of Mexico grids showing: (a) northern coarse parent grid, (b) Mississippi River 

delta nesting grid, (c) Mississippi River canyon nesting grid, (d) De Soto Canyon and coastal 

nesting grid, (e) De Soto canyon nesting grid, (f) northern nested grid, and (g) full coarse grid.
70

Figure 6.2: Full Gulf of Mexico coarse Grid-g (red box) and associated 1:5 refinement Grid-f 

(yellow box) for the Mississippi River delta. Lower panel shows a zoom of the grid meshes.
71

Figure 6.3: Mississippi River delta grids showing bird-foot delta: (a) coarser parent Grid-f and 

(b) fine resolution nested Grid-b with a refinement factor of 1:5 from Grid-f.

a.

b.
72

Figure 6.4: Mississippi canyon bathymetry (m) for nested grid c with an average horizontal 

resolution of the order of 100m.

Figure 6.5: De Soto canyon bathymetry (m) for nested Grid-d with an average horizontal 

resolution of the order of 100m.
73

Figure 6.6: Sea surface height (m) from the ROMS solution for the North Western Atlantic 

(NWA) used for initial and lateral boundary conditions (Courtesy of Enrique Curchitser, Rutgers 

University).

The tidal amplitude and currents (S2, M2, K1, O1 semidiurnal and diurnal components) 

forcing at the open boundaries were derived from the Oregon State University Tidal Prediction 

Software (OTPS) for the Gulf of Mexico. The model was forced with river runoff from various 

sources (see Section 5.4).

The modeling strategy was to run sequentially coarse and fine grids (g and f, respectively). 

The three-hour averaged solution of Grid-g was used as boundary conditions for Grid-f. The 

coupling between the two grids was one-way. These open boundary conditions worked very 

well as illustrated in Figure 6.7. A well-developed Loop Current and soon-to-be detached eddy 

can be seen. Notice that this eddy was centered on the southern lateral boundary of the finer –

resolution Grid-f. The sea-surface height contours were continuous between both grids making 

the transition almost invisible and without small-scale contamination. The lateral boundary 

exchange between donor and receiver grids was frequent and consistent. Both grids g and f were 

configured with the same algorithm options and forced with the same atmospheric and wave 

datasets.
74

ROMS Configuration

1. Third-order upstream-biased horizontal advection and fourth-order vertical advection 

for momentum (Shchepetkin and McWilliams, 2005, 2009; Haidvogel et al., 2008).

2. MPDATA (Multidimentional Positive Definite Advection Transport Algorithm) for active 

and passive tracers (Margolin and Smolarkiewicz, 1998).

3. Nonlinear equation of state for seawater (Jackett and McDougall, 1995).

4. Density Jacobian pressure algorithm with monotonized cubic polynomials (Shchepetkin 

and McWilliams, 2003).

5. Harmonic horizontal viscosity (Grid-g: 30 m2

/s; Grid-f: 10 m2/s) along constant Ssurfaces using the deviatory stress tensor (Wajsowicz, 1993).

6. Harmonic horizontal diffusion of active and passive tracers along geo-potential surfaces 

(Grid-g: 100 m2

/s; Grid-f: 10 m2/s).

7. Turbulent vertical mixing parameterization using the Generic Length Scale Method 

(Umlauf and Burchard, 2003) configured as a k-kl (Warner et al., 2005).

8. The bottom stresses were computed from the Bottom Boundary Layer (BBL) model 

(Warner et al., 2008) forced with 3-hour wave data from WaveWatch III.

9. Third-order predictor-corrector time-stepping (Shchepetkin and McWilliams, 2005, 

2009) with baroclinic time-steps of 150 s (Grid-g) and 30 s (Grid-f). 

10. Initial conditions from larger scale ROMS solution (see Figures 6.1 and 6.6).

11. One-way nesting of lateral boundary conditions from larger scale ROMS solutions (see 

Figure 6.7).

12. Atmospheric forcing from ECMWF ERA-Interim 3-hour dataset.

13. Tidal forcing from the OTPS dataset (Egbert and Erofeeva, 2002) for the Gulf of 

Mexico: S2, M2, K1, and O1 semidiurnal and diurnal components.

14. Surface net heat flux correction with monthly-averaged SST from the NOAA 

CoastWatch blended dataset (see Section 5.3).

15. River forcing from gauged and modeled data (see Section 5.4). The major river runoff 

was modeled with 13 (Grid-g) and 47 (Grid-f) point sources.

16. Simulation periods: Grid-g spin-up from 01-Jan-2000 to 31-Dec-2005, Grid-g

production runs from 01-Jan-2006 to 31-Dec-2012, and Grid-f sediment transport runs 

from 1-Oct-2007 to 30-Sep-2008 for selected strong hurricane season.

17. Output NetCDF files: daily instantaneous fields and 3-hour time-averaged fields. 

Solutions are available from

http://tds.marine.rutgers.edu/thredds/roms/gom/catalog.html.

Table 6.2: ROMS configuration options and main parameters for solutions with Grids g and f.
75

Figure 6.7: Simulated sea-surface height (m) seven days after initialization for a one-way 

nesting between: (a) coarse Grid-g and (b) fine Grid-f. The coarse to fine refinement ratio is 

1:5.

The grid-f solutions were run for the period from 1-Oct-2007 to 30-Sep-2008. The period 

for finer Grid-f was selected because the 2008 Atlantic hurricane season was very active in the 

Gulf of Mexico with several tropical storms (Edouard, Fay) and hurricanes (Dolly, Gustav, Ike). 

The 2008 Atlantic hurricane season tracks are shown in Figure 6.8 (from NOAA’s National 

Hurricane Center). The VIMS group ran Grid-f for the sediment transport simulations.

a

b
76

Figure 6.8: NOAA’s National Hurricane Center 2008 Atlantic hurricane season track map. Two 

named hurricanes impacted the study area during this time. Hurricane Gustav (Aug 25- Sep 4, 

2008; track 7) made landfall near Cocodrie, Louisiana on 1 September with maximum winds 

near 90kt (category 2); while Hurricane Ike (Sep 1-14, 2008; track 9) made landfall on 

September 13, 2008 on Galveston Island, Texas with maximum winds near 95 kt (category 2). 

The two storms had different storm tracks, with Gustav going over the Mississippi Canyon study 

area and Ike south of the study area.

The ROMS configuration scripts and input NetCDF files for Grids-f and -g can be 

downloaded from the following web site (www.myroms.org/arango/BOEM/ROMS). Several 

ROMS output NetCDF files and animations are available on this site for various solutions. In 

addition, an OpenDAP catalog (http://tds.marine.rutgers.edu/thredds/roms/gom/catalog.html) 

was set-up to store the current solutions for grids f and g.

It is beyond the scope of this report to show maps of the physical circulation solutions from 

grids g and f. The sediment transport solutions for Grid-f are discussed in Section 7. Figure 6.9

shows a monthly-averaged longitudinal cross-section along 25.5°N of potential temperature 

(colors) and sigma-t referenced to the surface (overlaid contours) for June 2007 (Fig. 6.9a) and 

January 2008 (Fig. 6.9b). The meandering in the upper 500m is due to Loop Current eddies and 

was more pronounced during the summer of 2007. Similarly, Figure 6.10 shows a monthlyaveraged latitudinal cross-section along 87.5°W of potential temperature (colors) and velocity 
77

magnitude (overlaid contours) for June 2007 (Fig. 6.10a) and January 2008 (Fig. 6.10b). Notice 

that the surface and sub-surface currents next to De Soto Canyon (~28°N) were stronger in 

winter 2008.

Figure 6.9: Monthly-averaged potential temperature (Celsius; colors) and sigma-0 (kg/m3; 

overlay; CI=0.5) longitudinal cross-section along 25.5°N for a Grid-g solution: (a) summer 2007 

and (b) winter 2008.

July 2007

January 2008

a

b
78

Figure 6.10: Monthly averaged sea potential temperature (Celsius) and velocity magnitude 

(m/s; overlay; CI=0.1) latitudinal cross-section along 87.5°W for a Grid-g solution: (a) summer 

2007 and (b) winter 2008.

July 2007

January 2008

a

b
79

7. CSTMS Simulations

To evaluate suspended sediment transport on the northern Gulf of Mexico, CSTMS 

routines were implemented within the three-dimensional Grid-f ROMS model for the time period 

of October 1, 2007 through September 30, 2008 (Fig. 6.2, lower panel). During this time frame, 

oceanographic conditions included episodic periods of energetic waves due to various types of 

extratropical storms, a period of elevated fluvial discharge, and two hurricanes, Gustav and Ike, 

during the late summer of 2008 (Fig. 7.1). Sediment transport was simulated using seven 

sediment class sizes, representing fluvial and seabed sources (Table 7.1). Two grain-size classes 

were represented in each fluvial source with duplicate classes used to distinguish sediment 

delivered by the Mississippi, Atchafalaya, and Mobile Rivers. Discharge data and rating curves, 

obtained from USGS personnel, were used to develop a time-series of sediment and freshwater 

discharge as described in Section 5.4. Three sediment classes represented a realistic, spatiallyvariable seabed interpolated to the model grid from dbSEABED data (Fig. 7.2, see Section 5.5). 

All other model input, including winds, waves, boundary conditions, tides, etc., derived from the 

hydrodynamic model setup outlined in Sections 5 and 6. Analyses of results that used different 

open boundary configurations revealed sensitivities in model estimates of currents and, therefore,

sediment transport (Fig. 7.3). The updated configuration had higher eastern and western 

elevations (on the order of 10 cm), and southern boundary elevations were nudged higher with 

time by a few cm. Boundary temperatures varied less over time than in the previous model 

setup. Similarly, salinities at the boundaries had a narrower range (less than 1) than the previous 

model configuration (greater than 5).

Class Source Sediment

Type D (mm) tcr (Pa) ws (mm/s)

1 Seabed Mud 0.063 0.11 1.0

2 Seabed Sand 0.125 0.13 10.0

3 Seabed Gravel 10.0 10.0 70.0

4 Mississippi River Micro-floc 0.015 0.11 0.1

5 Mississippi River Macro-floc 0.063 0.11 1.0

6 Atchafalaya/Mobile 

Rivers Micro-floc 0.015 0.03 0.1

7 Atchafalaya/Mobile 

Rivers Macro-floc 0.063 0.03 1.0

Table 7.1: Sediment parameters for the three-dimensional 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.
80

Throughout the year-long model run (October 2007 – September 2008), suspended 

transport and deposition of terrestrial sediments were limited along-shelf, such that the model 

indicated little comingling of material from the Mississippi, Atchafalaya, and Mobile Rivers 

(Fig. 7.4). Sediments delivered by the Mobile River were generally retained within Mobile Bay. 

Atchafalaya River sediments initially deposited proximally to the delta, and were subjected to 

frequent resuspension by waves and currents, and transported to the west, along the shelf. Crossshelf transport of Atchafalaya sediments was limited to depths of ~30 m. These results are 

consistent with results from Xu et al., 2011. Mississippi River plume dynamics largely 

influenced initial deposition, widely dispersing sediment a radial pattern around the birdfoot 

delta (Fig. 7.4A). A small portion of the Mississippi River load, however, was deposited beyond 

the shelf break during extratropical and tropical storms (> 200 m water depth). These various 

pathways have important implications for nutrient cycling, benthic habitats, and stratigraphic 

formation on the shelf, as well as implications for sediment supply to the deep sea (e.g. Bianchi 

et al., 2006).

Patterns of deposition and erosion resulting from wave- and current-resuspension over the 

simulated period revealed a high degree of spatial variability, extending from the coast and 

overlapping the edge of the shelf break (Fig 7.5). A closer inspection revealed that most erosion 

and deposition occurred during hurricanes Gustav and Ike, rather than during the preceding eight 

months, coinciding with a period of elevated Mississippi River discharge (Figs. 7.1 and 7.5D). 

Accordingly, the hurricane periods also contributed disproportionately to sediment flux to the 

slope Fig. 7.6), due to resuspension and cross-shelf transport. Periods of less elevated flux did 

occur between February and May 2008, as a result of plume dispersal combined with elevated 

waves during extratropical storms (7.6B). Buoyant dispersal during non-storm conditions 

contributed to sediment delivery to a lesser extent when currents drove the plume over the shelf 

break.

During less energetic wave events, such as winter storms, wave heights reached only about 

1 – 3 m, (Fig. 7.1), and were therefore energetic enough to suspend fine-grained sediment only in 

the shallow portions of the continental shelf (Fig. 7.7). Cross-shelf transport of the resuspended 

sediments during these events was generally limited because the sediments either rapidly settled 

to the seabed, or along-shelf currents dominated transport (Fig. 7.8 C and D). The shelf-break 

around the Mississippi River mouth, extending from the Mississippi Canyon to the DeSoto 

Canyon, received some sediment from the river plume during these winter storms when winds 

rotated from the south to the north, as cold fronts generally passed from west to east over the 

Gulf of Mexico. During south winds, the Mississippi River plume was constrained to the shelf 

by onshore currents, but as winds relaxed and then blew from the north, the plume began to 

extend out over the shelf break, resulting in cross-slope sediment delivery. The high frequency 

of the storms during the winter also contributed to enhanced delivery to the slope (Fig. 7.6 B).

Especially large cross-slope fluxes occurred during the two major hurricanes in the study 

period, Gustav and Ike (Figs. 7.6, 7.8, and 7.9). Because wave heights exceeded 3m during these 

events, bed stresses were sufficient to suspend fine-grained sediment at all continental shelf 

depths (Fig. 7.7). Transport across the shelf break and slope during hurricanes Gustav and Ike 

was a result of resuspension when elevated bed-shear stresses, sufficient to suspend fine grained 

sediment, reached water depths that overlapped the shelf break. Resuspended sediment was

transported downslope to depths of 1300 m at some locations (Figs. 7.6, 7.8, 7.9, and 7.10). 
81

Strong sediment flux episodes persisted for only a few hours, and the magnitude of flux 

generally decreased with water depth.

Hurricanes Gustav and Ike produced distinct patterns of transport and flux, largely due to 

difference in hurricane strength, duration, and storm track. Hurricane Ike persisted for a longer 

period, spending about twice as much time over the Gulf of Mexico as Gustav. Ike followed a 

more westerly path, making landfall near Galveston Texas, whereas the faster-moving Gustav

made landfall near the Mississippi delta (Fig. 6.8). Wave heights in the study area were larger 

for Ike than Gustav (Fig. 7.1). Consequently, lower sediment flux rates were observed during 

Hurricane Gustav, whereas hurricane Ike resulted in cross-slope transport at several locations 

(Figs. 7.8, 7.9, 7.11).

To better characterize suspended sediment fluxes across the shelf break, model calculations 

of depth-integrated flux were analyzed for five cross-slope transects (Fig. 7.6A). These data 

highlight the dominant role of the hurricanes in delivery large fluxes of suspended sediment from 

the continental shelf to the slope. Cross-slope fluxes were small during and after Hurricane 

Gustav west of the Mississippi Canyon, but created some of the largest sediment fluxes over the 

De Soto Canyon area (Figs. 7.9 and 7.11). In contrast, Hurricane Ike produced suspended 

sediment fluxes at all studied transects, but the magnitude of the suspended fluxes generally 

increased toward the west. During both storms, the suspended sediment flux was directed 

toward the southwest for water depths from 92 – 700 m, i.e. roughly across the Mississippi 

Canyon with a slight offshore component (Fig. 7.9). For the shallowest location analyzed (92 m) 

at the head of the Mississippi Canyon, the suspended sediment flux peak coincided 

approximately with the presence of storm waves on the shelf (Figs. 7.1, 7.9, 7.12). The peaks in 

sediment flux calculations down the Mississippi Canyon, however, lagged behind the passage of 

the storms, being delayed by 1 – 5 days relative to the occurrence of peak wave energy on the 

shelf at ~700 m water depth. Sediments resuspended on the shelf during the hurricanes were 

transported across the shelf break after the hurricanes began to cause relatively high resuspension 

on the west side of the Mississippi delta. As winds and along-shelf currents relaxed and changed 

direction, and sea-surface heights lowered, sediments were then transported down-slope, 

sometimes days after the hurricane (Figs. 7.10, and 7.12).

For the transect located in the DeSoto Canyon, however, the model estimated net 

downslope flux towards the south to southeast during the hurricanes (Fig. 7.9). The apparent 

increase in flux downslope (between 643 and 754 m) during and after the passage of hurricane 

Gustav (Aug25-Sep3) suggests possible contributions from suspended sediment transported over 

the sides of the canyon. Short periods (hours long) of transport apparently upslope, reflect the 

complex currents of the area, which occasionally form a gyre. Sediment transport in the DeSoto 

Canyon during Hurricane Ike appeared to begin at the 1115 m and then later at shallower depths. 

This occurred because a large suspension of sediment over the DeSoto Canyon was shifted to the 

east after Hurricane Ike passed through the area, and winds began blowing from the northwest.

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

et al., 2009); Hurricane Frederic in the DeSoto Canyon (Shay and Elsberry, 1987); Hurricane 

Allen (Brooks, 1983) and Hurricanes Katrina and Rita (Cole, 2010). Our model results indicate 

that the magnitude and direction of along-canyon fluxes is sensitive to storm intensity and storm 

track, as well as Mississippi River discharge.
82

Figure 7.1: Discharge (panel A in blue, m3 s-1), wave (panel A in red, m), surface wind 

momentum stress (panel B in black, Pa), and current (panel C, m s-1

) conditions during the 

modeled period at the DeSoto Canyon, -87.94, 29.11. In panel C, east (positive) and west 

(negative) currents in blue; whereas the red line denotes north (positive) and south (negative) 

currents. Discharge, wave, and wind inputs were based on data from the USGS, NOAA/NCEP 

WaveWatch III, and ECMWF; see Section 6). Hurricanes (H), tropical storms (T.S.) , and 

tropical depressions (T.D.) that impacted parts of the Gulf of Mexico denoted by gray bars in 

each panel: H. Noel (10/28-11/2/2007); T.S. Olga (12/9 - 12/14/2007); remnants of H. Arthur 

(6/5 – 6/7/2008); H. Dolly and T,S Christobal (7/20 -7/27/2008); T,S, Edouard (8/3 – 8/9/2008); 

T.D. Fay (8/18-8/26/2008); H. Gustav (8/25-9/4/2008); and H, Ike (9/1- -/14/2008).

A

B

C
83

Figure 7.2: Sediment bed fractions of (A) mud, (B) sand, (C) gravel used in the model. These 

correspond to sediment types 1, 2, and 3 in Table 7.1, respectively. Coastline shown in black; 

land in gray. White contours are bathymetry for depths of 20, 50, 100, 500, and 1000 m.

(A) Mud

(B) Sand

(C) Gravel

-92 -91 -90 -89 -88 -87

30

29

28

27

-92 -91 -90 -89 -88 -87

30

29

28

27

-92 -91 -90 -89 -88 -87

30

29

28

27

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4

0.6

0.8

1
84

Figure 7.3: Depth-averaged currents in the Mississippi (A) and DeSoto (B) canyons for two 

model runs that used different open boundary specifications. Period shown from October, 2007 

to September, 2008. The Mississippi Canyon location, -89.47, 28.39, is at 727 m depth, and the 

DeSoto Canyon location, -87.94, 29.11, is at 1052 m depth.
85

Figure 7.4: A) Final distribution calculated for sediment (cm) delivered from the A) Mississippi 

River B) Atchafalaya and Mobile River, and C) combined fluvial deposition, September 20, 

2008. Bathymetric contours (in black) drawn for depths of 10, and every 100 m up to 1500 m.

0.1 1 10 cm

A

B

C
86

Figure 7.5: A) Net erosion and deposition from suspended sediment transport for all sediment 

class sizes (see Table 7.1) calculated near the end of the simulation, September 20, 2008. 

Panel B represents the time-integrated deposition and erosion (cm) for the period prior to 

hurricanes Gustav and Ike, from Oct 1, 2007 up until August 25, 2008. Bathymetric contours (in 

black) drawn for depths of 10, and every 100 m up to 1500 m.

A

B
87

Figure 7.5 (cont.): Panels C and D represent time-integrated deposition and erosion (cm) 

during Hurricanes Gustav and Ike respectively. Bathymetric contours (in black) drawn for 

depths of 10, and every 100 m up to 1500 m.

D

C
88

Figure 7.6: Panel A shows locations of cross-shelf and slope transects. Data from the DeSoto 

Canyon (green circles), Mississippi Canyon (red circles), and western slope (black circles) are 

presented in Figure 7.9. Panel B shows time series depth-integrated horizontal flux (kg m-1 s-1

) 

for selected sites located east of the Mississippi Canyon (ME, blue circles), and the hurricane 

and storm record for the Gulf of Mexico. Note changes in scale in panel B. The bottom plot in B 

shows tropical storms (11 identified) in gray and extratropical storms (55 identified) in pink, 

which typically last 1-4 days. NOAA National Climatic Data Center Storm Events Data Base, 

Storm Prediction Center weather maps, and buoy 42040 winds used to compile storm list.

A

B
89

Figure 7.7: A) Time-averaged wave and current bed shear stresses (Pa) calculated over the 

simulated period and B) during Hurricane Ike, from September 4-20th, 2008. Bathymetric 

contours (in white) drawn for depths of 10, 100m, and every 100 m interval up to 1500 m.

A

B
90

Figure 7.8: Three-dimensional model estimates of time-averaged horizontal sediment flux 

during hurricanes Gustav and Ike combined (Aug 25-Sept 24, 2008). A) 1 m above the seabed 

(kg s-1

). B) depth-integrated (kg m-1 s-1). Red arrows depict flux direction and magnitude. 

Bathymetric contours (in black) drawn for depths of 10, and every 100 m up to 1500 m.

B

A
91

Figure 7.8 (cont.): Three-dimensional model estimates of time-averaged horizontal sediment 

flux during hurricanes Gustav and Ike combined (Aug 25-Sept 24, 2008). Panels C and D show 

time-averaged flux at 1 m above the seabed and depth integrated, respectively, during a period 

of high discharge, Jan 8 –Aug 24, 2008. Red arrows depict flux direction and magnitude. 

Bathymetric contours (in black) drawn for depths of 10, and every 100 m up to 1500 m.

C

D
92

Figure 7.9: See figure 7.6 for locations. Time-series of depth-integrated flux (kg m-2 s-1) down 

the western slope (A), the Mississippi Canyon (B), and the De Soto Canyon (C). Vector angles 

correspond to flux direction in accordance with map conventions. Text in B and C denotes grid 

locations (y,x), and water depth (m). Note difference in scale in figure (B), top panel.

A

B C
93

Figure 7.10: Panels A and B show suspended sediment along the Mississippi Canyon transect 

(see Fig. 7.6 for location) during and after Hurricane Gustav, 1 September 2008 and 6 

September 2008, respectively. The cross-slope distance is in km, while depth is in m.

A

B

Cs g m-3

Cs g m-3

0.001 0.01 0.1 1

0.001 0.01 0.1 1

0 20 40 60 80 100 120 
94

Figure 7.11: Three-dimensional model estimates of time-averaged sediment flux during 

hurricane Gustav: (1 m above the seabed; kg s-1 (A), depth-integrated; kg m -1 s-1 (B)), Aug 25-

Sept 3. Red arrows depict flux direction and magnitude. Bathymetric contours (in black) drawn 

for depths of 10, and every 100 m up to 1500 m.

B

A
95

Figure 7.11 (cont.): Three-dimensional model estimates of time-averaged sediment flux during 

hurricane Ike: (1 m above the seabed; kg s-1 (C), depth-integrated; kg m -1 s-1 (D)), Sep 4-Sept 

20, 2008. Red arrows depict flux direction and magnitude. Bathymetric contours (in black) 

drawn for depths of 10, and every 100 m up to 1500 m.

C

D
96

Figure 7.12: Conditions in the Mississippi Canyon during hurricanes Gustav and Ike. Surface 

momentum stresses (Pa) due to winds are shown in panel A, with arrows pointing in the 

compass direction that winds are blowing. Panels B, C, and D display surface elevation (m) in 

black, depth-averaged current (m s-1

) direction in blue arrows, and depth-integrated flux 

(kg m-1 s-1

) at 92 m, 688 m, and 1008 m depths. Locations in figure 7.6 A. Note changes in 

scale between panels.
97

8. HurriSlip Simulations

The overall results from the hurriSlip modeling of sediment failures and flows in the Gulf 

of Mexico under the extreme conditions of a hurricane can be described quantitatively by event 

likelihoods during storm passage, in terms of event locations relative to storm path, and also 

statistically relative to the model input parameters such as water depth.

Wave-Induced Mass failures. The module waveSlip for wave-induced mass failure, 

combined with the module liquiSlip for wave-induced weakening of the seafloor, concerns 

mainly outer shelf and upper slope locations. This is a consequence of the shelf-slope zonation 

of sediment types and the decay of wave-influence at the seabed as water depths increase. The 

mudflows observed on the Mississippi delta front (e.g., Sterling and Strohbeck 1975) are one 

manifestation of the processes being addressed in waveSlip.

For any particular area, the statistics of failure prediction by depth are affected by the 

geography of the storm relative to the area's bathymetry. This is apparent in spatial mappings of 

the modeled failure events. The modeling suggests that mass failures were likely at over 100 

locations in the area over 10 days during the passage of Ivan in the zone of the shelf-slope 

transition. Many of the failure surfaces are indicated at depths ~30m.

Wave-induced liquefaction appears to be very important to these scenarios. In many 

locations at storm peak, the effective shear strength of the seabed is decreased by up to 1kPa 

(from ~5kPa) due to the steady pore pressure buildup. Any transient pore pressure increases 

would be additional, further weakening the seabed. Indeed, the modeling shows that without the 

liquefaction, probably very few mass failures would occur in the region.

Although it is known that mass failures do shed turbidity currents (Mohrig and Marr 2003), 

our modeling did not address that phenomenon. 

Suspended Sediment Flows. The module suspendiSlip treated the possibility that bottomlayers of the water column that are heavy with suspended sediment might initiate turbidity 

currents to travel downslope from the shelf-slope transition zone. 

The modeling determined that during storm passage, wave-supported near-bottom 

suspended sediment flows would be generated very widely throughout the region, even from 

inner- to outer-shelf locations. This result conforms to observational evidence (e.g., Traykovski 

2000) and results from other modeling, for the widespread occurrence of cross-shelf flows of 

suspended sediment. The modeling also suggests that around storm climax, conditions suitable 

for ignition of turbid bottom flows from near-bottom layers of suspended sediment are 

continuously present near the shelf break. That is, flow initiations were indicated by the 

modeling at every modeled time-step. 

Over 100,000 occasions of flow ignition were determined for a 10-day period of the storm, 

with a spatial distribution as shown in Figure 4.8. These flows were identified using a threshold 

of a Richardson Number criterion. Therefore marginal and extreme cases are included in these 

counts. Looking at the statistics of how well the criteria were met, we see that only a very small 

number were extremely likely to become turbidity currents passing to the upper slope.
98

Our modeling suggests that sediment failure events are commonplace around the shelfslope transition under hurricane conditions. It extends - through time and across geography - the 

miscellaneous and scattered direct evidence for seabed failures and flows during hurricanes. The 

direct observations of failures and flows, such as the mudflows and some sand deposits on the 

upper slope could be employed in future studies to validate the discussed model-suite.

9. TURBINS Simulations

As a step towards simulating turbidity currents in the northern Gulf of Mexico, we 

performed a lock-exchange simulation using non-hydrostatic RANS, along a canyon, in a failure 

location identified by CU-Boulder team. This failure location corresponds to cell index of 

(252,44) in the grid for the project area, as shown in Figure 9.1. In the lock-exchange 

simulation, we used a domain 21km long along the streamwise (pathway) direction. In the 

vertical direction, the water depth varied from 130 m to 300m. We need at least 10 grid nodes 

per current height. The streamwise and spanwise grid spacings have to be similar to the vertical 

one. In the simulation, as the current height was 30m, a grid spacing of 3m was used in all three

directions. A time step of 0.6 seconds, calculated based on a modified CFL condition (CFL<0.5) 

involving both convective and viscous terms, was used.

Figure 9.1: Pathway used in the non-hydrostatic RANS simulation.

Figure 9.2 shows the initial condition representing this canyon in the Gulf of Mexico. 

Initially the concentration, which is a proxy for excess density due to suspended particles, is set 

to one in the lock-region and it is set to zero everywhere else. When the lock is released, a 

gravity current with a height of 30m starts to travel down the slope with a current velocity of 

approximately 1 m/sec. We ran the simulation for 11 hours in real time. Figures 9.3 and 9.4

show the concentrations calculated for 4 hours and 8 hours respectively after initiation of the 

current. As can be observed from the figures, the current is being diluted due to entrainment of 

the ambient ocean water, which does not contain any sediment. The current height decreases 

below the original height of 30m as it travels down the slope. Overall, the non-hydrostatic 

RANS simulation accurately captures the features of a gravity current travelling over such long 

distances. Work is underway to use the ROMS estimated near–bed current velocities and 

suspended sediment concentrations as input into the RANS model. Specifically, we are working 

Cell index:

I,J=252,44
99

to include the oceanographic conditions along the axis of De Soto Canyon during Hurricane Ike 

in September 2008, to demonstrate our ability to include oceanographic forcing within our 

turbidity current model.

Figure 9.2: Contours of concentration (proxy for excess density) at time = 0 hours. Here the 

sediment concentration is non-dimensionalized so that it lies between 0 and 1. A value of 0 

indicates clear water without any sediment.

Figure 9.3: Contours of concentration at time = 4 hours. Here the sediment concentration is 

non-dimensionalized so that it lies between 0 and 1. A value of 0 indicates clear water without 

any sediment.

Figure 9.4: Contours of concentration at time = 8 hours. Here the sediment concentration is 

non-dimensionalized so that it lies between 0 and 1. A value of 0 indicates clear water without 

any sediment.
100

10. Conclusions and Recommendations

The major conclusions drawn from this project are:

1. A workflow was developed to numerically model sediment fluxes from inner continental 

shelf or fluvial sources to the continental slope. All code is available, along with 

initializations, grids, and outputs. The model workflow is perhaps one of the more 

complex ever attempted. The range of models used include databases for sediment 

texture for the continental shelf and slope of the Gulf of Mexico; river discharge 

estimates; oceanographic conditions including hydrodynamic model estimates and 

suspended sediment transport; geotechnical models for flow ignition and slope failure; 

detailed turbidity current models. At present, the model workflow is sequential, without 

two-way coupling. Our results should therefore be considered first order, providing 

estimates of both location and volume of sediment made available for gravity flow 

transport into the deeper continental slope. Furthermore, sequential modeling is time 

consuming by nature. Model A often must be completed before Model B is started.

2. The challenges with regard to integrating the different modeling approaches across the 

different spatial and temporal scales are substantial, and they require further work, both 

on the research and the code development side. Both physical aspects (the 

implementation of erosion, resuspension of complex sediments into very large-scale 

simulation approaches) as well as numerical ones (one-way, two-way coupling, temporal 

and spatial interpolation at the boundaries between models etc.) will require additional 

work.

3. The immersed boundary/RANS approach in conjunction with the multiple successive 

streamwise modules appears to be well suited to perform Gulf of Mexico turbidity 

current simulations over realistic length and time scales.

4. 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 where 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 grid (g) was 

coupled with high-resolution grids to better resolve bathymetric features, river runoff, 

and sediment transport. The simulations over Grid-g cover the period from 01-Jan-2000 

to 31-Dec-2005 (spinoff) and from 01-Jan-2006 to 31-Dec-2012 (production runs). All 

the necessary atmospheric forcing fields, river discharge, and lateral boundary conditions 

were processed from 2000 to 2012. This allowed us to study the sediment transport 

scenarios due to idealized extreme atmospheric forcing events (storms and hurricanes of 

various intensities and different paths). We concentrated the sediment transport solutions 

(Grid-f) from 1-Oct-2007 to 30-Sep-2008 because of very active tropical storms and 

major hurricanes crossing the study area during this period.
101

5. The suspended sediment transport model indicates that episodic suspended transport 

down features like Mississippi and De Soto Canyons is fed principally by sediment fluxes

generated by wave resuspension. In particular, 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 was not a significant source of sediment delivery for the study period from 1-Oct2007 to 30-Sep-2008.

6. During the two Hurricanes modeled (Ike and Gustav), suspended sediment fluxes were 

predominantly offshore 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. 

7. On the continental slope, turbidity currents can be triggered by slope failure when storm 

driven supply forces accumulation sediment to deep water.

8. The weakest parts of our gravity flow model include the erosion and resuspension

boundary conditions used within the gravity flow module, and the potential for nonNewtonian dynamics in the dense suspension region close to the bed.

The main recommendations are:

1. The model workflow is perhaps one of the more complex ever attempted. As such, the 

workflow should be better exercised in ways beyond the original BOEM funding level.

2. 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, the eddy shedding from the Loop Current, 

and sporadic strong forcing from storms and hurricanes affects lateral open boundary 

conditions in such regional models. Therefore, a telescoping grid approach, from coarse 

(kilometers) to fine (meters) horizontal scales is required to obtain viable long-term (1-10 

years) and affordable circulations.

3. Joint modeling and field experiments are needed to develop reliable sediment transport 

models for the Gulf of Mexico continental shelf and slope. Our suspended sediment 

transport routines are based on parameterizations developed for other continental shelf 

systems. Water column and depositional data with which to compare the model 

calculations is severely lacking. Both to gain confidence in the models developed for this 

workflow, and to improve these models will require allied field studies of sediment 

processes for the northern Gulf of Mexico continental shelf and slope. Maximum benefit 

would be obtained via joint field studies and numerical modeling exercises. 
102

11. Appendix

Presentations based on research from “Shelf-Slope Sediment Exchange in the Northern 

Gulf of Mexico: Application of Numerical Models for Extreme Events”

2015

C.K. Harris, J. Birchler, T.A. Kniskern, D.R. Corbett. 2015. Linking sediment transport and 

radiosiotopes in a numerical model with application offshore of the Mississippi River, 

InterCoh Meeting, Leuven, Belgium, September, 2015.

Kniskern, T.A., C.K. Harris, H. Arango and D.R. Forrest. 2015. Storm-driven delivery of 

sediment to the continental slope: Numerical modeling for the northern Gulf of Mexico. 

Poster presented at the Community Surface Dynamics Modeling System (CSDMS) AllHands Meeting; Boulder, CO; May 2015. 

2014

Kniskern, T.A., C.K. Harris, and J.J. Birchler. 2014. Advances in sediment-transport modeling 

offshore of a fluvial source. Poster presented at the ChesMS14 (Chesapeake Modeling 

Symposium, 2014); Annapolis, MD; May 28, 2014. 

Kniskern, T.A., C.K. Harris, and J.J. Birchler. 2014. Advances in sediment-transport modeling 

offshore of a fluvial source. Poster presented at the Community Surface Dynamics 

Modeling System (CSDMS) All-Hands Meeting; Boulder, CO; May 2014. 

Birchler, J.J., C.K. Harris, H.G. Arango, E.W. Hutton, C. Jenkins, and J. Syvitski. 2014.

Delivery of sediment to the continental slope: numerical modeling tools for the northern 

Gulf of Mexico. EOS, Transactions of the American Geophysical Union Ocean Sciences 

Meeting Supplement, Abstract 16630, Honolulu, HI, February, 2014. 

2013

Birchler, J. and C.K. Harris. 2013. A numerical model to directly estimate sediment transport 

and short-lived radioisotopes: Application to the Mississippi shelf. Coastal and 

Estuarine Research Federation, November 6, 2013, San Diego, CA.

Birchler, J. and C.K. Harris. 2013. Sediment transport model that accounts for short-lived 

radioisotopes. Fourteenth International Conference on Estuarine and Coastal Modeling, 

November 4, 2013, San Diego, CA.

2012

Birchler, J.J. and C.K. Harris. 2012. A one-dimensional sediment-water column model 

incorporating radioisotope 7

Be. EOS, Transactions of the American Geophysical Union 

Ocean Sciences Meeting Supplement, Abstract B0898. Salt Lake City, Utah.
103

12. References

Adams, C.E., D.J.P. Swift, and J.M. Coleman, 1987: Sediment Transport on the Mississippi 

Prodelta Shelf, Winter 1984, Abstracts of Papers of the American Chemical Society, 194, 57–

ENVR.

Allison, M.A., G.C. Kineke, E.S. Gordon, and M.A. Goñi, 2000: Development and reworking of 

a seasonal flood deposit on the inner continental shelf off the Atchafalaya River, Cont. Shelf 

Res., 20(16), 2267–2294, doi:10.1016/S0278-4343(00)00070-4.

Bahr, D.B., E.W.H. Hutton, J.P.M. Syvitski, and L.F. Pratson, 2001: Exponential approximations 

to compacted sediment porosity profiles, Comput. & Geosci., 27(6), 691–700, 

doi:10.1016/S0098-3004(00)00140-0.

Bardet, J.P., 1998: Introduction to Computational Granular Mechanics, in Behavior of Granular 

Materials, B. Cambou, ed., International Centre for Mechanical Sciences, Springer-Verlag, 

Vienna, Austria, 79 pages.

Barnier, B., L. Siefridt, P. Marchesiello, 1995: Thermal forcing for a global ocean circulation 

model using a three-year climatology ECMWF analyses, J. Mar. Syst., 6(4), 363–380, 

doi:10.1016/0924-7963(94)00034-9.

Bever, A.J., 2010: Integrating space- and time-scales of sediment-transport for Poverty Bay, 

New Zealand, Ph.D. Dissertation, College of William and Mary, Virginia Institute of Marine 

Science.

Bianchi, T.S., M.A. Allison, E.A. Canuel, D.R. Corbett, B.A. McKee, T.P. Sampere, S.G. 

Wakeham, E. Waterson, 2006: Rapid export of organic matter to the Mississippi Canyon,

EOS Trans. AGU, 87(50), 565–573, doi:10.1029/2006EO500002.

Blaas, M., C. Dong, P. Marchesiello, J.C. McWilliams, and K.D. Stolzenbach, 2007: Sedimenttransport modeling on Southern California shelves: A ROMS case study, Cont. Shelf Res., 

27(6), 832–853, doi:10.1016/j.csr.2006.12003.

Blanchette, F., M. Strauss, E. Meiburg, B. Kneller, and M.E. Glinsky, 2005: High-resolution 

numerical simulations of resuspending gravity currents: Conditions for self-sustainment, J. 

Geophys. Res., 110, C12022, doi:10.1029/2005JC002927.

Boudreau, B.P., and R. H., Bennett, 1999: New Rheological and porosity equations for steadystate compaction, Am. J. Sci., 299(7-9), 517–528, doi:10.2475/ajs.299.7-9.517.

Brakenridge, GR, S. Cohen, A.J. Kettner, T. De Groeve, S.V. Nghiem, J.P.M. Syvitski, and B.M. 

Fekete, 2012: Calibration of satellite measurements of river discharge using a global 

hydrology model, J. Hydrology, 475, 123–136, doi:10.1016/j.jhydrol.2012.09.035.

Bratkovich, A., S. P. Diinnel, and D. A. Goolsby, 1994: Variability and prediction of freshwater 

and nitrate fluxes for the Louisiana-Texas shelf: Mississippi and Atchafalaya River source 

functions, Estuaries, 17(4), 766–778, doi:10.2307/1352746.
104

Brooke, S., and W. W. Schroeder. 2007: State of deep coral ecosystems in the Gulf of Mexico 

region. Texas to the Florida Straits. In S. E. Lumsden, T. F. Hourigan, A. W. Bruckner, and 

G. Dorr (Editors), The State of Deep Coral Ecosystems of the United States, 271–306. U.S. 

Department of Commerce, NOAA Technical Memorandum CRCP-3.

Brooks, D.A., 1983: The wake of Hurricane Allen in the western Gulf of Mexico, J. Phys. 

Oceanogr., 13(1), 117–129, doi:10.1175/1520-0485(1983)013<0117:TWOHAI>2.0.CO;2.

Burrough, P.A., and R. A. McDonnell, 1998: Principles of Geographical Information Systems. 

Spatial Information Systems and Geostatistics, Second Edition, Oxford University Press, 

Oxford U.K., 352 pages.

Cavin, A., M. Underwood, A.T. Fisher, and A. Johnston-Karas, 2000: Relations between textural 

characteristics and physical properties of sediments in northwestern Cascadia Basin, in

Fisher, A.T., Davis, E.E., and Escutia, C., eds., Proceedings of the Ocean Drilling Program, 

Scientific Results, v. 168: College Station, Texas (Ocean Drilling Program), p. 67–84 and 

CD-ROM.

Chen, C., R.O. Reid, and W.D. Nowlin Jr., 1996: Near-inertial oscillations over the TexasLouisiana shelf, J. Geophys. Res., 101(C2), 3509-3524, doi:10.1029/95JC03395.

Chen, X., Y. Li, X. Niu, M. Li, D. Chen, X. Yu, 2011: A general two-phase turbulent flow model 

applied to the study of sediment transport in open channels, Int. J. Multiphase Flow, 37(9), 

1099-1108, doi:10.1016/j.ijmultiphaseflow.2011.05.013.

Cho, K., R.O. Reid, and W.D. Nowlin Jr, 1998: Objectively mapped stream function fields on 

the Texas-Louisiana shelf based on 32 months of moored current meter data, J. Geophys. 

Res., 103(C5), 10377–10390, doi:10.1029/98JC00099.

Cochrane, J.D. and F.J. Kelly, 1986: Low-frequency circulation on the Texas-Louisiana 

continental shelf, J. Geophys. Res., 91(C9), 10645–10659, 10.1029/JC091iC09p10645.

Cohen, S., A.J. Kettner, and J.P.M. Syvitski, 2014: Global suspended sediment and water 

discharge dynamics between 1960 and 2010: Continental trends and intra-basin sensitivity,

Glob. Planet. Change, 115, 44–58, doi:10.1016/j.gloplacha.2014.01.011.

Cohen, S., A.J. Kettner, J.P.M. Syvitski, and B.M. Fekete, 2013: WBMsed: a distributed globalscale riverine sediment flux model: Model description and validation, Comput. Geosci., 53,

80–93, doi:10.1016/j.cageo.2011.08.011.

Cole, K., L., and S.F. DiMarco, 2010: Low-frequency variability of currents in the Deepwater 

Eastern Gulf of Mexico, U.S. Dept. of the Interior, Mineral Management Service, Gulf of 

Mexico OCS Region, New Orleans, LA. OCS Study MMS 2010-015, 136 pages.

Coleman, J.M. and L.D. Wright, 1975: Modern river deltas; variability of processes and sand 

bodies. M.L. Broussard (Editor). In: Deltas – Models for Exploration. Houston Geol. Soc., 

Houston, TX, p. 99–149.
105

Coleman, J.M., D.B. Prior, and L. E. Garrison, 1980: Subaqueous sediment instabilities in the 

offshore Mississippi River delta, Open File Rep 80-01, Bur. of Land. Manage., New Orleans, 

LA. 

Coleman, J.M., 1988: Dynamic Changes and Processes in the Mississippi River Delta, Geol. Soc. 

Am. Bull., 100(7), 999–1015, 10.1130/0016-7606(1988)100<0999:DCAPIT>2.3.CO;2.

Coleman, J.M., H.H. Roberts, and G.W. Stone, 1998: Mississippi River Delta: An Overview, J. 

Coastal Res., 14(3), 698–716.

Cordes, E.E., D.C. Bergquist, C.R. Fisher, 2009: Macro-ecology of Gulf of Mexico cold seeps,

Ann. Rev. Mar. Sci., 1, 143–168, doi:10.1146/annurev.marine.010908.163912.

Curray, J.R., 1960: Sediments and history of the Holocene transgression, continental shelf, 

northwest Gulf of Mexico. Recent Sediments, Northwest Gulf of Mexico, Amer. Assoc. 

Petrol. Geol. Bull., Tulsa, OK, 221–226.

De Wit, P.J. and C. Kranenburg, 1997: The wave-induced liquefaction of cohesive sediment 

beds, Estuar. Coast. Shelf Sci., 45(2), 261–271, doi:10.1006/ecss.1996.0184.

Dee, D.P., S.M. Uppala, A.J. Simmons, P. Berrisford, P. Poli, S. Kobayashi, U. Andrae, M.A. 

Balmaseda, G. Balsamo, P. Bauer, P. Bechtold, A.C.M. Beljaars, L. van de Berg, J. Bidlot, 

N. Bormann, C. Delsol, R. Dragani, M. Fuentes, A.J. Geer, L. Haimberger, S.B. Healy, H. 

Hersbach, E.V. Hólm, L. Isaksen, P. Kållberg, M. Köhler, M. Matricardi, A.P. McNally, 

B.M. Monge-Sanz, J.-J. Morcrette, B.-K. Park, C. Peubey, P. de Rosnay, C.Tavolato, J.N. 

Thépaut, and F. Vtart, 2011: The ERA-Interim reanalysis: configuration and performance of 

the data assimilation system, Q. J. R. Meteor. Soc., 137 (656), 553–597, doi:10.1002/qj.828.

Dickhudt, P.J., C.T. Friedrichs, L.C. Schaffner, and L.P. Sanford, 2009: Spatial and temporal 

variation in cohesive sediment erodibility in the York River estuary: a biologicallyinfluenced equilibrium modified by seasonal deposition, Mar. Geol., 267(3), 128–140, 

doi:10.1016/j.margeo.2009.09.009.

Dimakis, P., A. Elverhøi, K. Høeg, A., Solheim, C. B. Harbitz, J. S. Laberg, T. O. Vorren, and J. 

Marr, 2000: Submarine slope stability on high-latitude glaciated Svalbard-Barents Sea 

margin, Mar. Geol., 162(2-4), 303–316, doi:10.1016/S0025-3227(99)00076-6.

DiMarco, S.F. and R.O. Reid, 1998: Characterization of the principal tidal current constituents 

on the Texas-Louisiana shelf, J. Geophys. Res., 103(C2), 3093–3109, 

doi:10.1029/97JC03289.

DiMarco, S.F., M.K. Howard, and R.O. Reid, 2000: Seasonal variation of wind-driven diurnal 

cycling on the Texas-Louisiana continental shelf, Geophys. Res. Lett., 27(7), 1017–1020, 

doi:10.1029/1999GL010491.

Duncan, J.M., R.C. Horz, and T.L. Yang, 1989: Shear Strength Correlations for Geotechnical 

Engineering, Virginia Polytechnic Institute and State University, Blacksburg,VA, August, 

1989, 100 pages.
106

Durski, S.M., S.M. Glenn, and D.B. Haidvogel, 2004: Vertical mixing schemes in the coastal 

ocean: Comparison of the level 2.5 Mellor-Yamada scheme with an enhanced version of the 

K profile parameterization, J. Geophys. Res., 109(C1), C01015, doi:10.1029/2002JC001702.

Egbert, G.D. and S.Y. Erofeeva, 2002: Efficient inverse modeling of barotropic ocean tides, J. 

Atmos. and Oceanic Technol., 19(2), 183–204, doi:10.1175/1520-

0426(2002)019<0183:EIMOBO>2.0.CO;2.

Fairall, C.W, E.F. Bradley, J.E. Hare, A.A. Grachev, and J.B. Edson, 2003: Bulk 

Parameterization of Air-Sea Fluxes: Updates and Verification for the COARE Algorithm, J.

Clim., 16(4), 571–591, doi:10.1175/1520-0442(2003)016<0571:BPOASF>2.0.CO;2.

Fall, K.A., C.K. Harris, C.T. Friedrichs, J.P. Rinehimer, and C.R. Sherwood, 2014: Model 

behavior and sensitivity in an application of the cohesive bed component of the Community 

Sediment Transport Modeling System for the York River Estuary, VA, USA, J. Mar. Sci.

Eng., 2(2): 413–436; doi:10.3390/jmse2020413.

Fennel, K, R. Hetland, Y. Feng, and S. DiMarco, 2011: A coupled physical-biological model of 

the northern Gulf of Mexico shelf: Model description, validation and analysis of 

phytoplankton variability, Biogeosciences Discuss., 8(7), 1881–1899, doi:10.5194/bg-8-

1881-2011.

Ferguson, R.I., and M. Church, 2004: A simple universal equation for grain settling velocity, J.

Sediment. Res., 74(6), 933–937, doi:10.1306/051204740933.

Finn, W.D.L., R. Siddharthan, G.R. Martin, 1983: Response of seafloor to ocean waves, J. 

Geotech. Eng., 109(4), 556–572, doi:10.1061/(ASCE)0733-9410(1983)109:4(556).

Fisk, H.N., 1952: Geological Investigation of the Atchafalaya Basin and the Problem of 

Mississippi River Diversion, Vol. 1, US Army Corps of Engineers, Mississippi River 

Comm., Vicksburg, MS.

Friedrichs, C. T., L. D. Wright, D. A. Hepworth, and S. C. Kim, 2000: Bottom-boundary-layer 

processes associated with fine sediment accumulation in coastal seas and bays, Cont. Shelf 

Res., 20(7), 807–841, doi:10.1016/S0278-4343(00)00003-0.

Garcia, M. H., and G. Parker, 1993: Experiments on the entrainment of sediment into suspension 

by a dense bottom current, J. Geophys. Res., 98(C3), 4793–4807, doi:10.1029/92JC02404.

Gardner, J.V., B.R. Calder, J.E. Hughes Clarke, L.A. Mayer, G. Elston, and Y. Rzhanov, 2007:

Drowned shelf-edge deltas, barrier islands and related features along the outer continental 

shelf north of the head of De Soto Canyon, NE Gulf of Mexico, Geomorph., 89(3), 370–390, 

doi:10.1016/j.geomorph.2007.01.005.

Geremew, A.M., 2013: Pore-Water pressure development caused by wave-induced cyclic 

loading in deep porous formation, Int. J. Geomech., 13(1), 65–68, 

doi:10.1061/(ASCE)GM.1943-5622.0000189.

Germano, M., U. Piomelli, P. Moin, and W.H. Cabot, 1991: A dynamic subgrid-scale eddy 

viscosity model, Phys. Fluids, 3(7), 1760–1765, doi:10.1063/1.857955.
107

Gibson, R.E., 1958: The progress of consolidation in a clay layer increasing in thickness with 

time, Geotechnique, 8(4), 171–182, doi:10.1680/geot.1958.8.4.171.

Goff, J.A., C.J. Jenkins, and S.J. Williams, 2008: Seabed mapping and characterization of 

sediment variability using the usSEABED data base, Cont. Shelf Res., 28(4-5), 614–633, 

doi:10.1016/j.csr.2007.11.011.

Goñi, M.A., Y. Alleau, R. Corbett, J.P. Walsh, D. Mallinson, and M.A. Allison, 2007: The 

effects of hurricanes Katrina and Rita on the seabed of the Louisiana Shelf, Sedimentary 

Record, 5(1), 4–9.

Guidroz, W.S., 2009: Subaqueous, Hurricane-Initiated Shelf Failure Morphodynamics along the 

Mississippi River Delta Front, North-Central Gulf of Mexico, Unpubl. Ph.D. Thesis, 

Louisiana State University, 355 pages.

Gust, G. and V. Muller, 1997: Interfacial hydrodynamics and entrainment functions of currently 

used erosion devices, in Cohesive Sediments, Burt, Parker, and Watts, eds., John Wiley & 

Sons, West Sussex, U.K., 149–174.

Haidvogel, D.B., H.G. Arango, K. Hedstrom, A. Beckmann, P. Malanotte-Rizzoli, and A.F. 

Shchepetkin, 2000: Model Evaluation Experiments in the North Atlantic Basin: Simulations 

in Nonlinear Terrain-Following Coordinates, Dyn. Atmos. Oceans, 32(3-4): 239–281, 

doi:10.1016/S0377-0265(00)00049-X.

Haidvogel, D.B., H. Arango, W.P. Budgell, B.D. Cornuelle, E. Curchitser, E. Di Lorenzo, K. 

Fennel, W.R. Geyer, A.J. Hermann, L. Lanerolle, J. Levin, J.C. McWilliams, A.J. Miller, 

A.M. Moore, T.M. Powell, A.F. Shchepetkin, C.R. Sherwood, R.P. Signell, J.C. Warner, J. 

Wilkin, 2008: Ocean forecasting in terrain-following coordinates: Formulation and skill 

assessment of the Regional Ocean Modeling System, J. Comp. Phys., 227(7), 3595–3624, 

doi:10.1016/j.jcp.2007.06.016.

Halper, F.B. and D.W. McGrail, 1988: Long-term measurements of near-bottom currents and 

suspended sediment concentration on the outer Texas-Louisiana continental shelf, Cont. Shelf 

Res., 8(1), 23–36, doi:10.1016/0278-4343(88)90022-2.

Hunt, J. N., 1979: Direct solution of wave dispersion equation. J. of the Waterways, Port, 

Coastal, and Ocean Division, 105(4), ASCE, 457–459.

Hurlburt, H.E. and J.D. Thompson, 1980: A Numerical study of loop current intrusions and eddy 

shedding, J. Phys. Oceanogr., 10(10), 1611–1651.

Hutton, E., and J. Syvitski, 2004: Advances in the numerical modeling of sediment failure during 

the development of a continental margin, Mar. Geol., 203(3-4), 367–380.

Jackett, D. R. and T. J. McDougall, 1995: Minimal Adjustment of Hydrographic Profiles to 

Achieve Static Stability, J. Atmos. Oceanic Technol., 12(2), 381–389, doi:10.1175/1520-

0426(1995)012<0381:MAOHPT>2.0.CO;2.

Jackson, D.R. and M.D. Richardson, 2007: High-Frequency Seafloor Acoustics, The Underwater 

Acoustics Series, Springer Science & Business Media, New York, NY, 616 pages.
108

Jaffe, B.E., J.H. List, A.H. Sallenger Jr., 1997: Massive sediment bypassing on the lower 

shoreface offshore of a wide tidal inlet–Cat Island Pass Louisiana, Mar. Geol., 136(3-4), 

131–149, doi:10.1016/S0025-3227(96)00050-3.

Jaramillo, S., A. Sheremet, M. Allison, and E. Rogers, 2008: Mud induced wave-dissipation in 

the Atchafalaya Shelf, Louisiana, USA, Paper presented at the 2008 Ocean Sciences 

Meeting, Orlando, FL.

Jeng, D.S., B. Seymour, and J. Li, 2006: A new approximation for pore water pressure 

accumulation in marine sediment due to water waves, University of Sydney, Department of 

Civil Engineering, Research Report No R868.

Jeng, D.S., B. Seymour, and J. Li, 2007: A new approximation for pore pressure accumulation in 

marine sediment due to water waves, Int. J. Numer. Anal. Methods Geomech., 31(1), 53–69, 

doi:10.1002/nag.547.

Jenkins, C.J., 2015: Sediment Drainage Streams in Benthic Seafloor Classification, ESRI Ocean 

Solutions, Earth Solutions, 2nd EDITION, (submitted).

Julien, P.Y., 1998: Erosion and Sedimentation. Cambridge University Press, 280 pages.

Kindinger, J.C.L., 1988: Seismic stratigraphy of the Mississippi-Alabama shelf and upper 

continental slope, Mar. Geol., 83(1-4), 79–94, doi:10.1016/0025-3227(88)90053-9. 

Kineke, G.C., E.E. Higgins, K. Hart, and D. Velasco, 2006: Fine-sediment transport associated 

with cold-front passages on the shallow shelf, Gulf of Mexico, Cont. Shelf Res., 26(17-18), 

2073–2091, doi:10.1016/j.csr.2006.07.023.

Kuo, M.Y-H., M.D. Bolton, A.J. Hill, M. Rattley, 2010: New evidence for the origin and 

behavior of deep ocean ‘crusts’, Proc. 2nd int. Symposium on Frontiers in Offshore 

Geotechnics, ISFOG, Perth.

Landsea, C.W., 1993: A climatology of intense (or major) Atlantic hurricanes, Mon. Wea. Rev., 

121(6), 1703–1713, doi:10.1175/1520-0493(1993)121<1703:ACOIMA>2.0.CO;2.

Liu, J.Y. and W.R. Bryant, 2000: Sea floor morphology and sediment paths of the northern Gulf 

of Mexico deepwater, in Fine-Grained Turbidite Systems, Bouma, A.H. and C.G. Stone, eds.,

AAPG Memoir 72/SEPM Special Publication 68, 33–45.

Maaß, C., and U. Schumann, 1996: Direct numerical simulation of separated turbulent flow over 

a wavy boundary. In: Hirschel EH, editor. Flow simulation with high performance 

computers. Notes on numerical fluid mechanics, 52, p. 227–241.

Margolin, L. and P.K. Smolarkiewicz, 1998: Antidiffusive velocities for multipass donor cell 

advection, SIAM J. Sci. Comput., 20(3), 907–929, doi:10.1137/S106482759324700X.

Maxey M.R., and J.J. Riley, 1983: Equation of motion for a small rigid sphere in a nonuniform 

flow, Phys. Fluids, 26, 883–889.
109

McKee, B.A., R.C. Aller, M.A. Allison, T.S. Bianchi, and G.C. Kineke, 2004: Transport and 

transformation of dissolved and particulate materials on continental margins influenced by 

major rivers: benthic boundary layer and seabed processes, Cont. Shelf Res., 24(7-8), 899–

926, doi:10.1016/j.csr.2004.02.009.

Meade, R.H. and J.A. Moody, 2010: Causes for the decline of suspended-sediment discharge in 

the Mississippi River system, 1940-2007, Hydrol. Process., 24(1), 35–49, 

doi:10.1002/hyp.7477.

Mei, C.C. and M.A. Foda, 1981: Wave-induced responses in a fluid-filled poro-elastic solid with

a free surface-a boundary layer theory, Geophys. J. Roy. Astron. Soc., 66(3), 597–631, 

doi:10.1111/j.1365-246X.1981.tb04892.x.

Mohrig, D., and J. G. Marr, 2003: Constraining the efficiency of turbidity current generation 

from submarine debris flows and slides using laboratory experiments, Mar. Pet. Geol., 20(6-

8), 883–899, doi:10.1016/j.marpetgeo.2003.03.002.

Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, and J. Zavala-Garay, 2011a: 

The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation 

systems, Part I: System overview and formulation, Prog. Oceanogr., 91(1), 34–49, 

doi:10.1016/j.pocean.2011.05.004.

Moore, A.M., H.G. Arango, G. Broquet, C. Edwards, M. Veneziani, B. Powell, D. Foley, J.D. 

Doyle, D. Costa, and P. Robinson, 2011b: The Regional Ocean Modeling System (ROMS) 4-

dimensional variational data assimilation systems, Part II: Performance and Application to 

the California Current System, Prog. Oceanogr., 91(1), 50–73, 

doi:10.1016/j.pocean.2011.05.003.

Moore, A.M., H.G. Arango, G. Broquet, C. Edwards, M. Veneziani, B. Powell, D. Foley, J.D. 

Doyle, D. Costa, and P. Robinson, 2011c: The Regional Ocean Modeling System (ROMS) 4-

dimensional variational data assimilation systems, Part III: Observation Impact and 

Observation Sensitivity in the California Current System, Prog. Oceanogr., 91(1), 74–94, 

doi:10.1016/j.pocean.2011.05.005.

Morehead, M.D., J.P. Syvitski, E.W.H. Hutton, and S.D. Peckham, 2003: Modeling the temporal 

variability in the flux of sediment from ungauged river basins, Global Planet. Change, 39(1-

2), 95–110, doi:10.1016/S0921-8181(03)00019-5.

Morinishi, Y., T.S. Lund, O.V. Vasilyev, and P. Moin, 1998: Fully conservative higher order 

finite difference schemes for incompressible flow, J. Comp. Phys., 143(1), 90–124, 

doi:10.1006/jcph.1998.5962.

Nasr-Azadani, M.M., and E. Meiburg, 2011: TURBINS: An immersed boundary, Navier–Stokes 

code for the simulation of gravity and turbidity currents interacting with complex 

topographies, Comput. Fluids, 45(1), 14–28, doi:. 

Necker, F., C. Härtel, L. Kleiser, and E. Meiburg, E., 2002: High-resolution simulations of 

particle-driven gravity currents, Int. J. Multiphase Flow, 28(2), 279–300, doi:10.1016/S0301-

9322(01)00065-9.
110

Necker, F., C. Härtel, L. Kleiser, and E. Meiburg, 2005: Mixing and dissipation in particledriven gravity currents, J. Fluid Mech., 545, 339–372, doi:10.1017/S0022112005006932.

Neill, C.F. and M.A. Allison, 2005: Subaqueous deltaic formation on the Atchafalaya Shelf, 

Louisiana, Mar. Geol., 214(4), 411–430, doi:10.1016/j.margeo.2004.11.002.

Nielsen, P. 1986: Suspended sediment concentration under waves, Coast. Eng. 10(1), 23–31, 

doi:10.1016/0378-3839(86)90037-2. 

Nielsen, P. 1992: Coastal bottom boundary layers and sediment transport, World Scientific, 

Advanced Series on Ocean Engineering, 4, 340 pages.

Nittrouer, J. A., M. A. Allison, and R. Campanella, 2008: Bedform transport rates for the 

lowermost Mississippi River, J. Geophys. Res., 113(F3), F03004, 

doi:10.1029/2007JF000795.

Nodine, M.C., R.B. Gilbert, S.G. Wright, J.Y. Cheon, M. Wrzyszczynski, M. Coyne, and E.G. 

Ward, 2007: Impact of Hurricane-induced mudslides on pipelines, Proc. Offshore 

Technology Conference, Houston, Texas, OTC paper No. 18983.

Nodine, M.C., S.G. Wright, and R.B. Gilbert, 2006: Mudslides during Hurricane Ivan and an 

Assessment of the Potential for Future Mudslides in the Gulf of Mexico, Phase I Project 

Report Prepared for the Minerals Management Service, MMS Project Number 552, 37 pages.

Oey, L.-Y., H.-C. Lee, and W. J. Schmitz Jr., 2003: Effects of winds and Caribbean eddies on the 

frequency of Loop Current eddy shedding: A numerical model study, J. Geophys. Res., 

108(C10), 3324, doi:10.1029/2002JC001698.

Pantin, H. M., and M. C. Franklin, 2009: Predicting autosuspension in steady turbidity flow: 

Ignition point and their relation to Richardson numbers, J. Sediment. Res., 79(12), 862–871, 

doi:10.2110/jsr.2009.095.

Parchure, T. and A. Mehta, 1985: Erosion of Soft Cohesive Sediment Deposits, J. Hydraul. Eng., 

111(10), 1308–1326, doi:10.1061/(ASCE)0733-9429(1985)111:10(1308).

Parker, G., 1982: Conditions for the ignition of catastrophically erosive turbidity currents, Mar.

Geol., 46(3-4), 307–327, doi:10.1016/0025-3227(82)90086-X. 

Parker, G., Y. Fukushima, and H.M. Pantin, 1986: Self-accelerating turbidity currents, J. Fluid 

Mech., 171, 145–181, doi:10.1017/S0022112086001404.

Paulson, C.A. and J.J. Simpson, 1977: Irradiance measurement in the upper ocean, J. Phys. 

Oceanogr., 7(6), 952–956, doi:10.1175/1520-0485(1977)007<0952:IMITUO>2.0.CO;2.

Pepper, D.A. and G.W. Stone, 2002: Atmospheric forcing of fine-sand transport on a low-energy 

inner shelf: south-central Louisiana USA, Geo-Mar. Lett., 22(1), 33–41, doi:10.1007/s00367-

002-0094-x.
111

Prior, D. B., J.M. Coleman, and L.E. Garrison, 1979: Digitally acquired undistorted side scan 

sonar images of submarine landslides, Mississippi Delta, Geol., 7(9), 423–425, 

doi:10.1130/0091-7613(1979)7<423:DAUSSI>2.0.CO;2.

Rhodes, R.C., A.J. Wallcraft, and J.D. Thompson, 1985: Navy corrected geostrophic wind set for 

the Gulf of Mexico, NORDA Technical Note 310, Stennis Space Center, Mississippi, 103 

pages.

Rinehimer, J.P., C.K. Harris, C.R. Sherwood, and L.P. Sanford, 2008: Estimating cohesive 

sediment erosion and consolidation in a muddy, tidally-dominated environment: model 

behavior and sensitivity, Estuarine and Coastal Modeling, Proceedings of the Tenth 

Conference, November 5-7, Newport RI. M.L. Spaulding, ed. 819–838, 

doi:10.1061/40990(324)44.

Roberts, H.H., 1998: Delta Switching: Early Responses to the Atchafalaya River Diversion, J. 

Coastal Res., 14(3), 882–899.

Ross, C.B., W.D. Gardner, M.J. Richardson, and V.L. Asper, 2009: Currents and sediment 

transport in the Mississippi Canyon and effects of Hurricane Georges, Cont. Shelf Res., 

29(11-12), 1384–1396, doi:10.1016/j.csr.2009.03.002.

Sadler, P. M., 1981: Sediment accumulation rates and the completeness of stratigraphic sections, 

J. Geol., 89(5), 569–584, doi:10.1086/628623.

Safak, I., A. Sheremet, M.A. Allison, and T.-J. Hsu, 2010: Bottom turbulence on the muddy 

Atchafalaya Shelf, Louisiana, USA, J. Geophys. Res., 115(C12), C12019, 

doi:10.1029/2010JC006157

Sandwell, D. T., R.D. Müller, W.H.F. Smith, E. Garcia, and R. Francis, 2014: New global 

marine gravity model from CryoSat-2 and Jason-1 reveals buried tectonic structure, Science, 

346(6205), 65–67, doi:10.1126/science.1258213.

Sanford, L.P., 2008: Modeling a dynamically varying mixed sediment bed with erosion, 

deposition, bioturbation, consolidation, and armoring, Comput. Geosci., 34(10), 1263–1283, 

doi:10.1016/j.cageo.2008.02.011.

Senneset, K. and N. Janbu,1985: Shear Strength Parameters obtained from Static Cone

Penetration Tests, Strength Testing of Marine Sediments: Laboratory and In-Situ 

Measurements, ASTM STP 883, Chaney, R.C. and K.R. Demars, Eds., American Society for 

Testing and Materials, Philadelphia, 41–54. 

Shay, L.K., and R.L. Elsberry, 1987: Near-inertial ocean current response to Hurricane Frederic,

J. Phys. Oceanogr., 17(8), 1249–1269, doi:10.1175/1520-0485(1987)017<1249:NIOCRT>

2.0.CO;2.

Shchepetkin, A.F., and J.C. McWilliams, 2003: A method for computing horizontal pressuregradient force in an oceanic model with a nonaligned vertical coordinate, J. Geophys. Res., 

108(C3), 3090, doi:10.1029/2001JC001047.
112

Shchepetkin, A.F., and J.C. McWilliams, 2005: The Regional Oceanic Modeling System

(ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model, 

Ocean Modelling, 9(4), 347–404, doi:10.1016/j.ocemod.2004.08.002.

Shchepetkin, A.F., and J.C. McWilliams, 2009: Computational kernel algorithms for fine-scale, 

multi-process, long-term oceanic simulations, In: Handbook of Numerical Analysis: 

Computational Method for the Ocean and the Atmosphere, eds. T. Temam & J. Tribbia, 

Elsevier Science, 14, 121–183, doi:10.1016/S1570-8659(08)01202-0.

Sheremet, A. and G.W. Stone, 2003: Observations of nearshore wave dissipation over muddy sea 

beds, J. Geophys. Res., 108(C11), 3357, doi:10.1029/2003JC001885.

Sheremet, A., A. J. Mehta, B. Liu, G. W. Stone, 2005: Wave-sediment interaction on a muddy 

inner shelf during Hurricane Claudette, Estuar. Coast. Shelf Sci., 63(1), 225–233, 

doi:10.1016/j.ecss.2004.11.017.

Spalart, P.R., R.D. Moser, and M.M. Rogers, 1991: Spectral methods for the Navier-Stokes 

equations with one infinite and two periodic directions, J. Comp. Phys., 96(2), 297–324. 

Spierenburg, S. E. J., 1987: Seabed response to water waves. Masters thesis, Delft Univ. of

Technology, Delft, The Netherlands.

Spierenburg, S. E. J., 1989: Wave-induced stresses in a layered seabed in relation to the 

liquefaction potential, In Numerical Models in Geomechanics, NUMOG III, S. Pietruszczak 

and G. N. Pande, Eds., Elsevier Applied Science, 407–414.

Sterling, G.H., and G.E. Strohbeck, 1975: The failure of South Pass 70B platform in Hurricane 

Camille, J. Petrol. Technol., 27(03), 263–268, doi:10.2118/4475-PA.

Stone, G.W. and R.A McBride, 1998: Louisiana Barrier Islands and their Importance in Wetland 

Protection: Forecasting Shoreline Change and Subsequent Response of Wave Climate, J. 

Coast. Res., 14(3), 900–915.

Sturges, W. and R. Leben, 2000: Frequency of Ring Separations from the Loop Current in the 

Gulf of Mexico: A Revised Estimate, J. Phys, Oceanogr., 30(7), 1814–1819.

Styles, R. and S.M. Glenn, 2000: Modeling stratified wave and current boundary layers in the 

continental shelf, J. Geophys. Res., 105(C10), 24119–24139, doi:10.1029/2000JC900115.

Sumer, B.M., and J. Fredsøe, 2002: The Mechanics of Scour in the Marine Environment, 

Advanced Series on Ocean Engineering, Vol. 17, 552 pages.

Syvitski, J.P.M., and J.D. Milliman, 2007: Geology, geography, and humans battle for 

dominance over the delivery of fluvial sediment to the coastal ocean, J. Geol., 115(1), 1–19, 

doi:10.1086/509246.

Tarboton, D. G., 2003: Terrain Analysis Using Digital Elevation Models in Hydrology, 23rd 

ESRI International Users Conference, San Diego, California, July 7-11.
113

Tarboton, D. G. and M.E. Baker, 2008: Towards an Algebra for Terrain-Based Flow Analysis,

In: Representing, Modeling and Visualizing the Natural Environment: Innovations in GIS 13, 

Eds: N. J. Mount, G., Harvey, L., Aplin, P. & Priestnall, G., CRC Press, Florida, 167–194.

Teague, W.J., E. Jarosz, T.R. Keen, D.W. Wang, and M.S. Hulbert, 2006: Bottom scour 

observed under Hurricane Ivan, Geophys. Res. Lett., 33(7), L07607, 

doi:10.1029/2005GL025281.

Thompson, D.M., E.E. Wohl, and R.D. Jarrett, 1999: Velocity reversals and sediment sorting in 

pools and riffles controlled by channel constrictions, Geomorphology, 27(3-4), 229–241, 

doi:10.1016/S0169-555X(98)00082-8.

Thomsen, L. and G. Gust, 2000: Sediment erosion thresholds and characteristics of resuspended 

aggregates on the western European continental margin, Deep-Sea Res. I, 47(10), 1881–

1897, doi:10.1016/S0967-0637(00)00003-0.

Tolman, H.L. and D.V. Chalikov, 1994: Development of a third-generation ocean wave model at 

NOAA-NMC, Proc. Waves - Physical and numerical modelling, M. Isaacson and M.C. 

Quick Eds., Department of Civil Engineering, University of British Columbia, Vancouver,

Canada, 724–733.

Toniolo, H., M. Lamb, G. Parker, 2006: Depositional turbidity currents in diapiric minibasins on 

the continental slope: Formulation and theory, J. Sediment. Res., 76(5), 783–797.

Traykovski, P., W.R. Geyer, J.D. Irish, J.F. Lynch, 2000: The role of wave-induced densitydriven fluid mud flows for cross-shelf transport on the Eel River continental shelf, Cont. 

Shelf Res., 20(16), 2113–2140, doi:10.1016/S0278-4343(00)00071-6.

Umlauf, L. and H. Burchard, 2003: A generic length-scale equation for geophysical turbulence 

models, J. Mar. Res., 61(2), 235–265, doi:10.1357/002224003322005087.

van Ledden , M., W.G.M. van Kesteren, J.C. Winterwerp, 2004: A conceptual framework for the 

erosion behavior of sand-mud mixtures, Cont. Shelf Res., 24(1), 1–11, 

doi:10.1016/j.csr.2003.09.002.

Verruijt, A., Computational Geomechanics, Kluwer Academic Publishers, Dordrecht, 1995.

Viesca, R. C., and J. R. Rice, 2012: Nucleation of slip-weakening rupture instability in landslides 

by localized increase of pore pressure, J. Geophys. Res., 117(B3), B03104, 

doi:10.1029/2011JB008866.

Wajsowicz, R.C, 1993: A consistent formulation of the anisotropic stress tensor for use in 

models of the large-scale ocean circulation, J. Comp. Phys., 105(2), 333-338, 

doi:10.1006/jcph.1993.1079.

Walsh, J.P., R. Corbett, D. Mallinson, M. Goñi, M. Dail, C. Loewy, K. Marciniak, K. Ryan, C. 

Smith, A. Stevens, B. Sumners, and T. Tesi, 2006: Mississippi Delta Mudflow Activity and 

2005 Gulf Hurricanes, Eos Trans. AGU, 87(44), 477–478, doi:10.1029/2006EO440002.
114

Walsh, J.P. and C.A. Nittrouer, 2009: Understanding fine-grained river-sediment dispersal on 

continental margins, Mar. Geol., 263(1-4), 34–45, doi:10.1016/j.margeo.2009.03.016. 

Warner, J.C, C.R. Sherwood, H.G. Arango, and R.P. Signell, 2005: Performance of four 

Turbulence Closure Methods Implemented using a Generic Length Scale Method, Ocean 

Modelling, 8(1-2), 81––113, doi:10.1016/j.ocemod.2003.12.003.

Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. Arango, 2008: Development 

of a three-dimensional, regiona1, coupled wave, current, and sediment-transport model, 

Comput. Geosci., 34(10), 1284–1306, doi:10.1016/j.cageo.2008.02.012.

Warner, J.C., W.R. Geyer, and H.G. Arango, 2010: Using a composite grid approach in a 

complex coastal domain to estimate estuarine residence time, Comput. Geosci., 36(7), 921–

935, doi:10.1016/j.cageo.2009.11.008.

Weaver, D.C., G.D. Dennis, and K.L. Sulak, 2001: Northeastern Gulf of Mexico Coastal and 

Marine Ecosystem Program: Community Structure and Trophic Ecology of Demersal Fishes 

on the Pinnacles Reef Tract; Final Synthesis Report, U.S. Department of the Interior, 

Geological Survey, USGS BSR-2001-0008 and Minerals Management Service Gulf of 

Mexico OCS Region, New Orleans, LA, OCS Study MMS 2002-034, 92 pp.

Wiberg, P.L. and C.R. Sherwood, 2008: Calculating wave-generated bottom orbital velocities 

from surface-wave parameters, Comput. Geosci., 34(10), 1243–1262, 

doi:10.1016/j.cageo.2008.02. 010.

Wiegel, R.L., 1964. Oceanographical Engineering, Prentice Hall, Inc, Englewood Cliffs, NJ, 

USA.

Wijesekera H.W., D.W. Wang, W.J. Teague, and E. Jarosz, 2010: High sea-floor stress induced 

by extreme hurricane waves, Geophys. Res. Lett., 37(11), L11604, 

doi:10.1029/2010GL043124.

Williams, S.J., M.A. Arsenault, B.J. Buczkowski, J.A. Reid, J.G. Flocks, M.A. Kulp, S. Penland, 

and C.J. Jenkins, 2006: Surficial sediment character of the Louisiana offshore Continental 

Shelf region: a GIS Compilation, U.S. Geological Survey Open-File Report 2006-1195, 

online at http://pubs.usgs.gov/of/2006/1195/index.htm.

Winter, W. R., R. J. Jahnert, and A. B. França, 2007: Campos Basin, Petrobras Geosciences

Bulletin, 15(2), 511–529.

Winterwerp, J.C., W. G. M. van Kesteren, B. van Prooijen, and W. Jacobs, 2012: A conceptual 

framework for shear flow-induced erosion of soft cohesive sediment beds. J. Geophys. Res., 

117(C10), C10020, doi:10.1029/2012JC008072.

Wisser, D., B.M. Fekete, C.J. Vörösmarty, and A.H. Schumann, 2010: Reconstructing 20th 

century global hydrography: a contribution to the Global Terrestrial Network-Hydrology 

(GTN-H), Hydrol. Earth Syst. Sci., 14(1), 1–24, doi:10.5194/hess-14-1-2010.
115

Wollheim, W. M., C. J. Vörösmarty, B. J. Peterson, P. A. Green, S. Seitzinger, J. Harrison, A. F. 

Bouwman, and J. P. M. Syvitski, 2008: Global N removal by freshwater aquatic systems: a 

spatially distributed within basin approach, Global Biogeochem. Cycles, 22(2), GB2026, 

doi:10.1029/2007GB002963.

Wright, L.D. and Coleman, J.M., 1971: Effluent expansion and interfacial mixing in the presence 

of a salt wedge, Mississippi River delta, J. Geophys. Res., 76(36), 8649–8661, 

doi:10.1029/JC076i036p08649.

Wright, L.D. and J.M. Coleman, 1974: Mississippi River Mouth Processes - Effluent Dynamics 

and Morphologic Development, J. Geol., 82(6), 751–778. 

Wright, L.D. and C.A. Nittrouer, 1995: Dispersal of river sediments in coastal seas, Estuaries, 

18(3), 494–508, doi:10.2307/1352367.

Wright, L.D., C.R. Sherwood, and R.W. Sternberg, 1997: Field measurements of fairweather 

bottom boundary layer processes and sediment suspension on the Louisiana inner continental 

shelf, Mar. Geol., 140(3-4), 329–345, doi:10.1016/S0025-3227(97)00032-7.

Xu, K., C.K. Harris, R.D. Hetland, J.M. Kaihatu, 2011: Dispersal of Mississippi and Atchafalaya 

sediment on the Texas–Louisiana shelf: Model estimates for the year 1993, Cont. Shelf Res., 

31(15), 1558–1575, doi:10.1016/j.csr.2011.05.008.

Xu, K., D.R. Corbett, J.P. Walsh, D. Young, K.B. Briggs, G.M. Cartwright, C.T. Friedrichs, C.K. 

Harris, R.C. Mickey, and S. Mitra. 2014: Seabed erodibility variations on the Louisiana 

continental shelf before and after the 2011 Mississippi River flood, Estuar. Coast. Shelf Sci., 

149, 283–293, doi:10.1016/ j.ecss.2014.09.002.

Yamamoto, T., H. L. Koning, H. Sellmeijer, and E. Van Hijum, 1978: On the response of a poroelastic bed to water waves, J. Fluid Mech., 87(01), 193–206.

Youd T.L., 1973: Factors controlling maximum and minimum densities of sands, Evaluation of 

Relative Density and Its Role in Geotechnical Projects Involving Cohesionless Soils, ASTM, 

STP 523, 98–112, doi:10.1520/STP37866S.
116

The Department of the Interior Mission

As the Nation's principal conservation agency, the Department of the Interior has 

responsibility for most of our nationally owned public lands and natural 

resources. This includes fostering sound use of our land and water resources; 

protecting our fish, wildlife, and biological diversity; preserving the 

environmental and cultural values of our national parks and historical places; and 

providing for the enjoyment of life through outdoor recreation. The Department 

assesses our energy and mineral resources and works to ensure that their 

development is in the best interests of all our people by encouraging stewardship 

and citizen participation in their care. The Department also has a major 

responsibility for American Indian reservation communities and for people who 

live in island territories under US administration.

The Bureau of Ocean Energy Management Mission

The Bureau of Ocean Energy Management (BOEM) promotes energy 

independence, environmental protection and economic development through 

responsible, science-based management of offshore conventional and renewable 

energy resources.

The BOEM Environmental Studies Program Mission

The mission of the Environmental Studies Program (ESP) is to provide the 

information needed to predict, assess, and manage impacts from offshore energy 

and marine mineral exploration, development, and production activities on 

human, marine, and coastal environments.