OCS Study

BOEM 2017-030

U.S. Department of the Interior

Bureau of Ocean Energy Management

Gulf of Mexico OCS Region

Remote Sensing Assessment of 

Surface Oil Transport and Fate 

During Spills in the Gulf of Mexico

Synthesis Report
OCS Study 

BOEM 2017-030

U.S. Department of the Interior

Bureau of Ocean Energy Management

Gulf of Mexico OCS Region

Remote Sensing Assessment of 

Surface Oil Transport and Fate 

During Spills in the Gulf of Mexico 

Synthesis Report 

Authors 

Ian MacDonald, PhD 

Dmitry Dukhovskoy, PhD 

Mark Bourassa, PhD 

Steven Morey, PhD 

Oscar Garcia-Pindea, PhD 

Samira Daneshgar Asl 

Chuanmin Hu, PhD 

Mark Reed, PhD 

Jørgen Skancke, PhD 

Prepared under BOEM Award 

Contract # M12PC00003 

by

Florida State University, Center for Ocean-Atmospheric Prediction Studies (COAPS) 

2000 Levy Avenue, Suite 292, Tallahassee, FL 32310 

University of South Florida, College of Marine Science 

830 1st St SE, St. Petersburg, FL 33701 

Stiftelsen for Industriell og Teknisk Forskning (SINTEF)

Strindvegen 4, 7034 Trondheim, Norway
iii

DISCLAIMER

This report was prepared under contract between the Bureau of Ocean Energy Management 

(BOEM) and the Florida State University. This report has been technically reviewed by BOEM, 

and it has been approved for publication. Approval does not necessarily signify that the contents 

reflect the views and policies of BOEM, nor does mention of trade names or commercial 

products constitute endorsement or recommendation for use.

REPORT AVAILABILITY

To download a PDF file of this Gulf of Mexico OCS Region 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 2017-030.

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

MacDonald, I.R. D. Dukhovskoy, M. Bourassa, S. Morey, O. Garcia-Pineda, Daneshgar Asl, S. 

C. Hu, M. Reed, and J. Skancke. 2017. Remote Sensing Assessment of Surface Oil Transport 

and Fate during Spills in the Gulf of Mexico. US Dept. of the Interior, Bureau of Ocean Energy 

Management, Gulf of Mexico Region, New Orleans, LA OCS Study BOEM 2017-030. 129 pp.

ABOUT THE COVER

Report cover shows an oil slick in the Gulf of Mexico. (Photo credit: I.R. MacDonald.)
iv

ACKNOWLEDGEMENTS

Editorial assistance of T. Ippolito in preparing reporting products of this study is gratefully 

acknowledged. 
v

CONTENTS

1. Project Overview, Tasks, and Investigators.............................................................................. 1

1.1. Investigators at FSU and Roles...................................................................................... 2

1.2. Consulting Investigators and Roles .............................................................................. 3

2. Remote Sensing of Surface Oil Distributions ........................................................................... 4

2.1. SAR Texture Classifying Neural Network Algorithm (TCNNA).................................. 5

2.2. SAR Oil Emulsion Detection Algorithm (OEDA)......................................................... 6

2.3. Optical Remote Sensing: MODIS and AVIRIS ............................................................. 7

2.4. Remote Sensing of Natural Seeps................................................................................ 12

2.5. Chronic, Small-scale Oil Discharges............................................................................ 14

2.6. Remote Sensing Results for Deepwater Horizon Oil Spill .......................................... 18

3. Surface Oil Transport and Weathering Model and Forcing Fields ......................................... 21

3.1. Surface Oil Drift Model with Weathering in the GoM ................................................ 21

3.1.1. Setup of Deepwater Horizon 2010 oil spill simulation in OSCAR................ 21

3.1.2. Adaptation of the OSCAR model for the WAMOST Project......................... 22

3.2. Forcing Fields and Validation and Adjustment of the Oil Drift Model ....................... 23

4. Effects of Oceanographic and Wind Forcing on Surface Oil.................................................. 63

4.1. Wind Forcing................................................................................................................ 63

4.2. Ocean Surface Currents................................................................................................ 64

4.3. Surface Waves.............................................................................................................. 65

4.4. Parameterization of Oil Surface Roughness................................................................. 85

5. Surface Oil Transport Constraints and Barriers ...................................................................... 87

5.1. Boundaries inferred from Instantaneous SSH and Surface Velocity Fields ................ 87

5.2. Material Boundaries Based on Lagrangian Coherent Structures ................................. 92

5.3. Transport of Surface Oil in the Presence of Material Boundaries ............................... 93

6. Future Research Needs.......................................................................................................... 101

6.1. Continued publication ................................................................................................ 101

6.2. Remote sensing of surface oil .................................................................................... 101

6.3. Forcing fields and validation...................................................................................... 102

6.4. Surface transport constraints...................................................................................... 103

7. Literature Cited...................................................................................................................... 104
vi

LIST OF FIGURES

Figure 1-1. Flow chart of tasks and deliverables............................................................................ 2

Figure 2-1. TCNNA processing of SAR images ............................................................................ 6

Figure 2-2. Areas of thick oil and oil emulsion delineated by the Oil Emulsion Detection 

Algorithm (OEDA) in a SAR image of the DwH surface oil June 3, 1010 ............... 7

Figure 2-3. Examples of AVIRIS products for oil thicknesses ...................................................... 8

Figure 2-4. Statistics of sun glint strengths (LGN) where oil slicks can or cannot 

be detected from MODIS images ............................................................................... 9

Figure 2-5. Oil slicks from the Ixtoc oil spill in 1979 .................................................................... 9

Figure 2-6. Examples of the three different MODIS oil distribution maps.................................. 11

Figure 2-7. Natural seep zones in the GoM.................................................................................. 13

Figure 2-8. Examples of the oil slicks extracted from ENVISAT ASAR images 

using TCNNA........................................................................................................... 15

Figure 2-9. Composite of TCNNA and manually extracted oil slicks associated with the 

Taylor site in BOEM lease block MC020: August 2010 to January 2012 ............... 16

Figure 2-10. Location of GC600 on the bathymetric map of the northern GoM.......................... 17

Figure 2-11. Locations of floating oil slicks extracted by TCNNA from 41 SAR images........... 18

Figure 2-12. Map of surface oil from DwH visible to SAR: Distribution and 

average volume of surface oil (m3 km-2) from DwH discharge.............................. 19

Figure 2-13. Time series of surface oil from the DwH discharge ................................................ 20

Figure 3-1. OSCAR simulation of surface oil during DwH discharge ......................................... 23

Figure 3-2. Coordinate surfaces in the HYCOM and NCOM in the upper 15 m......................... 26

Figure 3-3. Demeaned SSH fields (m) from CCAR, NCOM and HYCOM 

time-averaged over May 12–19, 2010 ...................................................................... 27

Figure 3-4. Maximum northern extent of the LC estimated from the analysis 

and CCAR SSH fields .............................................................................................. 28

Figure 3-5. Comparison of the LC fronts from the NCOM and HYCOM vs CCAR................... 29

Figure 3-6. Near-surface currents from HYCOM GoM analyses 20.1

(a) and 31.0 (b) and reanalysis 50.1 (c) on May 8, 2010, 00Z ................................ 30

Figure 3-7. SSS ship observations used for NCOM and HYCOM validation ............................. 31

Figure 3-8. Scatter plot of SSS data from the models vs ship observations................................. 32

Figure 3-9. Scatter plot of SST data from the models vs ship observations (WTER data only).. 33

Figure 3-10. Vertical T and S profiles from ARGO buoy measurements and 

NCOM and HYCOM................................................................................................ 35

Figure 3-11. CCAR SSH fields on selected dates......................................................................... 36
vii

Figure 3-12. Ocean Color Index (left column) and SSS fields from NCOM 

(middle) and HYCOM.............................................................................................. 38

Figure 3-13. Mean MHD for OCI–HYCOM and OCI-NCOM.................................................... 39

Figure 3-14. RMS difference between CCMP and COAMPS ..................................................... 41

Figure 3-15. Wind speed bias (a) and wind vector difference (b) between 

CCMP and COAMPS............................................................................................... 42

Figure 3-16. Spreading of surface oil during passage of atmospheric front, 24-28 April 2010 ... 43

Figure 3-17. Spreading of surface oil during passage of atmospheric front, 

30 April–5 May 2010................................................................................................ 44

Figure 3-18. Spreading of surface oil during passage of atmospheric front, 

6 June–10 June 2010................................................................................................. 45

Figure 3-19. Locations of the NDBC buoys and closest CCMP and COAMPS grid points........ 46

Figure 3-20. CCMP, COAMPS, and NDBC buoy wind speeds................................................... 47

Figure 3-21. NDBC buoy 42040 surface winds and SSM/I wind speed retrieval 

before (top) and after (bottom) the DwH oil spill..................................................... 47

Figure 3-22. Schematic diagram of the point-to-point matching algorithm ................................. 50

Figure 3-23. Weighted Mean Displacements measured for the data sets A (a) and B (b)............ 51

Figure 3-24. Sensitivity tests with three randomly generated shapes (a–c).................................. 53

Figure 3-25. Metric scores against scale for three shapes (a–c) shown in Figure 3-24................ 54

Figure 3-26. Metric scores against clockwise rotation angle (degrees) for three shapes (a–c) .... 54

Figure 3-27. Metric scores for the linear translation test that the three shapes (in the rows) 

are subjected to ......................................................................................................... 55

Figure 3-28. Metric scores (vertical axis) from the noise test for the three shapes (in columns). 56

Figure 3-29. Random shapes subsampled from a continuous field .............................................. 57

Figure 3-30. Metric scores against scale for three concentration fields (a–c) 

shown in Figure 3-29 a–c ......................................................................................... 58

Figure 3-31. Metric scores against clockwise rotation angle (degrees) for 

three concentration fields (a–c)................................................................................. 58

Figure 3-32. Metric scores for the linear translation test (Figure 3-29g) that 

three shapes (in the rows) are subjected to ............................................................... 59

Figure 3-33. Metric scores (vertical axis) from the noise test for the three 

concentration fields inside the shapes (in columns) ................................................. 60

Figure 3-34. Metric scores (vertical axis) from the sensitivity experiments 

listed in Table 3-5 ..................................................................................................... 58

Figure 3-35. Time-integrated SAR-TCNNA surface oil volume (m3/km2), 

May 23, 2010, 00Z.................................................................................................... 59
viii

Figure 3-36. Contours of the time-integrated surface oil volume (0.09 m3/km2), 

May 23, 2010, 00Z.................................................................................................... 59

Figure 3-37. Skill metric scores for time-integrated surface oil volume from the surface oil 

drift model simulations with varying half-life parameter and SAT/TCNNA data ... 60

Figure 3-38. Contours of the time-integrated surface oil volume = (0.09 m3/km2,), 

May 23, 2010, 00Z.................................................................................................... 61

Figure 3-39. Skill metric scores for time-integrated surface oil volume from the 

surface oil drift model simulations with varying wind coefficient parameter 

and SAT-TCNNA data ............................................................................................. 62

Figure 4-1. Magnitude of vector difference between velocity vector of topmost 

grid cell and upper 10 m averaged velocity from HYCOM and NCOM ................. 64

Figure 4-2. Near-surface velocity profiles of the wind-driven current for three 

conceptual models .................................................................................................... 65

Figure 4-3. Example of NCEP peak wave period data (s), 1 May 2010 00Z with the 

problematic gap circled............................................................................................. 68

Figure 4-4. Maps showing extent of the model domains used for Wavewatch III....................... 69

Figure 4-5. Example of WW3 modeled peak wave period data (s), 1 May 2010 00Z................. 70

Figure 4-6. Probability density function of Stokes drift magnitude distributions 

in the GoM, April–July 2010.................................................................................... 71

Figure 4-7. Directional rose plot of Stokes drift distributions in the GoM, April–July 2010 ...... 72

Figure 4-8. Probability density function of Stokes drift magnitude distributions 

for deep water (depth> 100 m) in the Gulf Mexico, April–July 2010...................... 74

Figure 4-9. Probability density function of Stokes drift magnitude distributions 

for shelf water (depth  100 m) in the Gulf Mexico, April–July 2010..................... 75

Figure 4-10. Probability density function of 24-hour Stokes drift trajectory displacements 

for deep water (depth > 100 m) in the GoM, April–July 2010................................. 76

Figure 4-11. Probability density function of 24-hour Stokes drift trajectory displacements 

for shelf water (depth  100 m) in the Gulf Mexico, April–July 2010..................... 76

Figure 4-12. Probability density function of Stokes drift to surface current ratios for 

deep water (depth > 100 m) in the GoM, April–July 2010. ..................................... 77

Figure 4-13. Probability density function of Stokes drift to surface current ratios for 

shelf water (depth  100 m) in the GoM, April–July 2010 ...................................... 78

Figure 4-14. Probability density function of Stokes drift to wind drift ratios for 

deep water (depth > 100 m) in the Gulf Mexico, April–July 2010 .......................... 79

Figure 4-15. Probability density function of Stokes drift to wind drift ratios for 

shelf water (depth  100 m) in the Gulf Mexico, April–July 2010 .......................... 79

Figure 4-16. Map of Hurricane Alex storm track ......................................................................... 81
ix

Figure 4-17. Probability density functions of Stokes drift magnitudes in the 

area of the GoM directly affected by Hurricane Alex .............................................. 83

Figure 4-18. Probability density functions of Stokes drift magnitudes in the 

area of the GoM affected by the oil spill .................................................................. 83

Figure 5-1. Okubo-Weiss parameter (x109) estimated from the HYCOM 

surface layer velocities ............................................................................................. 88

Figure 5-2. Relative vorticity (x105 s-1) estimated from the HYCOM 

surface layer velocities ............................................................................................. 89

Figure 5-3. Gradient of the SSH from HYCOM........................................................................... 91

Figure 5-4. SAR TCNNA compared to simulation results........................................................... 94

Figure 5-5. Example of the LCS ridge extraction procedure........................................................ 95

Figure 5.6. Surface salinity gradient (PSU/m, color map) from the HYCOM 

GoM Reanalysis (50.1)........................................................................................ 96-97

Figure 5-7. SSH gradient from the HYCOM GoM Reanalysis (50.1) .................................... 98-99

Figure A1. LC identification algorithm based on the Kalman Filtering..................................... 117

LIST OF TABLES

Table 2-1. Summary statistics for inventory of oil-covered water detected in 

176 independent SAR image samples of the GoM................................................... 14

Table 3-1. Key setup information for simulation of the DwH oil spill 

with the OSCAR model............................................................................................ 22

Table 3-2. Ocean analysis data sets .............................................................................................. 25

Table 3-3. Time and coordinates of ARGO float measurements used in 

the validation analysis............................................................................................... 34

Table 3-4. Wind data sets.............................................................................................................. 40

Table 3-5. Oil model parameters .................................................................................................. 63

Table 3-6. Oil model parameters for the half-life time estimation ............................................... 59

Table 4-1. Stokes drift magnitudes, entire GoM, April–July 2010 .............................................. 73

Table 4-2. Stokes drift magnitudes, deep and shelf water, GoM, April–July 2010...................... 74

Table 4-3. Dates of weeks before, during and after Hurricane Alex for study consideration ...... 80

Table 4-4. Stokes drift magnitudes (km/day) of the weeks surrounding Hurricane Alex

in the area directly impacted by the storm................................................................ 82

Table 4-5. Stokes drift magnitudes (km/day) of the weeks surrounding Hurricane Alex 

in area impacted by oil spill...................................................................................... 82

Table 5-1. Oil model parameters for the half-life time estimation ............................................... 93
x

ABBREVIATIONS AND ACRONYMS

AVIRIS Airborne Visible/Infrared Imaging Spectrometer

CCAR Colorado Center for Astrodynamics Research

CFSR Climate Forecast System Reanalysis

COAMPS Coupled Ocean/Atmosphere Mesoscale Prediction System

COAPS Center for Ocean-Atmospheric Prediction Studies

DwH Deepwater Horizon

EOAS Earth Ocean and Atmospheric Sciences Department 

(at Florida State University)

FSU Florida State University

GoM Gulf of Mexico

LC Loop Current

LCE Loop Current eddies

LGN normalized glint reflectance

MHD modified Hausdorff distance

NAVOCEANO Naval Oceanographic Office

NCEP National Center for Environmental Prediction

NCODA Navy Coupled Ocean Data Assimilation

NOGAPS Navy Operational Global Atmospheric Prediction System

NRDA National Resource Damage Assessment

OCI ocean color index

OEDA oil emulsion detection algorithm

OSCAR Oil Spill Contingency and Response (model)

OSO oil slick origin

SAR synthetic aperture radar

SINTEF Stiftelsen for Industriell og Teknisk Forskning 

Norwegian Technical and Industrial Research Institute)

SIIRS visible infrared imaging radiometer suite

TCNNA Texture Classifying Neural Network Algorithm 

(for surface oil detection)

USF University of South Florida
1

1. PROJECT OVERVIEW, TASKS, AND INVESTIGATORS

The project entitled Wind and Advection Modeling for Oil Spill Transport (WAMOST) was 

a multi-disciplinary effort conducted by researchers at the Florida State University (FSU), 

University of South Florida (USF), and the Norwegian Technical and Industrial Research 

Institute - Stiftelsen for Industriell og Teknisk Forskning (SINTEF). The team completed a 

comprehensive and technically innovative analysis of synthetic aperture radar (SAR) and optical 

remote sensing images of the Deepwater Horizon (DwH) oil spill and of lesser discharges from 

natural and anthropogenic sources. Results of these analyses provided validation for 

development of weathering and advection models for oil spill tracking and for investigation of 

background natural seeps in the Gulf of Mexico (GoM). BOEM is charged with improving oil 

spill risk analysis and adding to understanding the mechanisms that control surface oil transport 

by investigation of persistent natural and anthropogenic oil releases and of the DwH oil spill.

During the DwH emergency, the scientific community produced a number of impromptu 

forecasts of surface oil drift to complement those performed by the federal agencies. These 

results were used to guide the oil cleanup operations in the Gulf and, after the spill emergency 

had passed, to inform investigations of oil fate and effects for assessment of damages. Major 

sources of uncertainty in the oil spill forecasts stemmed from the absence of reliable initialization 

data and uncertainties affecting surface oil drift models. Many of these predictions used

Lagrangian oil drift models without weathering processes, but SINTEF has developed a wellregarded weathering model, known as OSCAR (Oil Spill Contingency and Response model). Oil 

particle advection outcomes from the SINTEF OSCAR model were compared to results from 

models developed by FSU’s Center for Ocean Atmospheric Prediction Studies (COAPS) and a 

significant adaptation of the OSCAR model was implemented to capture the chemical state and 

relative age of oil particles. 

Most surface oil drift models combine wind effect and ocean surface currents to estimate 

drift of oil (typically modeled as discrete Lagrangian particles representative of fixed quantities 

of oil) in the ocean. Generally, calculation of oil trajectories employs ocean surface currents 

derived from hydrodynamic models. Though oil surface drift occurs at the very surface of the 

ocean, “surface currents” simulated in a hydrodynamic model more precisely represent the 

velocity integrated (or averaged) throughout the topmost vertical layer, which may be several 

meters thick, depending on the discrete vertical resolution. Thus, the so-called “surface current” 

from a model may deviate substantially from the true surface current (that is, the velocity of the 

water parcels directly on the ocean’s surface). The WAMOST results were able to implement 

and validate alternative approaches for estimation of the surface ocean currents that are crucial to 

accurate modeling of surface oil transport from hydrodynamic model velocity fields and nearsurface winds. The approach was to downscale the layer mean velocity from the ocean model to 

the surface current. We will assess the importance of this modification to wind and surface 

current fields for oil spill modeling.

Many existing methods for oil transport estimates ignore the influence of wave-driven 

transport. The WAMOST team investigated this to better understand limitations and advantages 

of using surface wave drift in oil transport modeling. The implemented model considered

contributions to surface roughness from three types of surface features that include an 

aerodynamically smooth surface (Nikuradse, 1933) (Kondo, 1975), capillary waves (Bourassa, et 

al., 1999), and gravity waves (Smith, et al., 1992). Results informed the following objectives: (1) 
2

Identification of material boundaries in the GoM; (2) Analysis of oil mixing across the 

boundaries; (3) Quantification of critical strength values to define permeability of the material 

boundaries to surface oil. 

The components of this project were closely integrated and compared model output to remote 

sensing observations for respective validation. The flowchart below (Figure 1-1) shows how

project components and tasks were integrated.

Figure 1-1. Flow chart of tasks and deliverables

1.1. INVESTIGATORS AT FSU AND ROLES

Dr. MacDonald: Professor of Oceanography, Department of Earth Ocean and Atmospheric 

Science (EOAS) at FSU; served as the Principal Investigator and Project Manager. MacDonald 

has previous experience managing multi-disciplinary, multi-institutional projects in related areas 

of marine science for the Minerals Management Service, the Department of Energy, and NOAA, 

and service for the National Resource Damage Assessment process.

Dr. Dmitry Dukhovskoy: Associate Research Scientist, COAPS; led the effort to validate 

forcing fields for the surface oil drift model, to develop and test objective methods for skill 

assessment of the oil model, develop and validate the surface oil drift model, employ the skill 

assessment algorithm to evaluate the performance of the oil model and adjust its parameters. 

Dukhovskoy also used the surface oil model to analyze the role of wind forcing on oil spreading. 

He investigated several approaches to define ocean fronts and developed and tested in realistic 

application the Kalman Filter-based technique for ocean front tracking. 

Dr. Mark Bourassa: Associate Professor of Meteorology, EOAS Department at FSU, and 

Associate Director, COAPS; incorporated and validated methodology for parameterization of the 

surface oil roughness and implemented it into the surface-oil drift model. 

Dr. Steve Morey: Associate Research Scientist, COAPS; worked with the FSU and SINTEF 

team members to integrate the oil weathering and transport algorithms into a modeling system 
3

for application to oil spill events and evaluation with remotely sensed data. He performed

analyses of model and observational data to identify material boundaries and their importance to 

horizontal mixing processes, and contributed to writing scientific manuscripts and project 

reports.

Dr. Oscar Garcia-Pineda: Assistant Scholar Scientist, EOAS Department at FSU; developed

and implemented methodologies for tracking oil spill and oil characteristics from Synthetic 

Aperture Radar satellite images.

Ms. Samira Daneshgar Asl: Research Assistant, EOAS Department at FSU; assisted with data 

coding and analysis and contributed research on chronic anthropogenic oil spills, natural seeps, 

and the DwH surface oil.

1.2. CONSULTING INVESTIGATORS AND ROLES

Dr. Chuanmin Hu: Professor of Oceanography, USF; an optical oceanographer who uses 

optical techniques (laboratory, field, and remote sensing) to study marine algal blooms and 

pollution including oil spills. He is currently leading the optical oceanography lab at the USF

College of Marine Science. His expertise is in atmospheric correction and optical algorithms to 

differentiate various ocean materials and organisms, including oil slicks, Trichodesmium blooms, 

Sargassum blooms, Ulva prolifera dark water slicks, Karenia brevis red tides, and turbidity 

events.

Dr. Mark Reed: Senior Researcher, SINTEF; oversaw design of simulation experiments, 

execution of simulations, reporting of results, and participated in preparing reports and 

publications. 

Dr. Jørgen Skancke: Researcher, SINTEF; oversaw design of simulation experiments, 

execution of simulations, reporting of results, and participated in preparation of reports and 

publications.

The Scientific Review Group consisted of: Dr. Villy Kourafalou, Research Professor, 

University of Miami Rosentiel School of Marine and Atmospheric Science; and Mr. William 

Pichel, Physical Scientist, NOAA-NESDIS, Office of Satellite Applications and Research.
4

2. REMOTE SENSING OF SURFACE OIL DISTRIBUTIONS

BOEM is responsible for oil spill risk assessment and understanding the fate and transport of 

oil to inform environmental documents. However, when a potential spill is detected, ascertaining 

the rate and source of the discharge can be problematic. Even in cases of a very large discharge,

such as the DwH oil spill, uncertainties can be quite profound, particularly in the immediate 

aftermath of the incident, when initial estimates of the rate of release turn out to be very different 

from the actual numbers (McNutt, et al., 2011). Remote sensing has a crucial role in addressing 

this problem. Results of work carried out during the WAMOST project help advance the state of 

the art for remote sensing of oil by improving the methods for detection of floating oil and in 

understanding the influences of environmental factors, such as wind and currents that spread and 

dissipate surface slicks. This section of the report describes several parallel efforts carried out 

under Task 2 of the project; these comprise the following topics, which also contribute to results 

that helped the project team accomplish Tasks 3, 4, and 5. 

 New and enhanced methods for semi-automated processing of SAR imagery;

 Analysis of optical remote sensing data for better understanding of oil layer thicknesses, 

and oil slick morphology;

 Background oil slicks from natural seeps quantified through analysis of the SAR archive;

 Surveillance of chronic and episodic small-scale releases; and

 A comprehensive time-series of surface oil released during the DwH incident.

Remote sensing techniques for oil slick detection and characterization include optical (both 

passive and active), microwave, and radar sensors that can be deployed on aircrafts and satellites. 

Each of these techniques has its own advantages and disadvantages. All these techniques were

used during the DwH oil spill, with varying success. The remote sensing approach in this project 

is based on multi-source and multi-platform techniques, with specific focus on radar and passive 

optical remote sensing. Indeed, SAR and optical remote sensing were two primary means used 

during the DwH spill to monitor the spill’s spatial distributions on a near-daily basis. 

SAR has been used to detect the presence of oil on the ocean surface for many years (Espedal 

& Wahl, 1999; De Beukelaer, et al., 2003; Brekke & Solberg, 2005). During the DwH oil spill 

response, the NOAA National Environmental Satellite, Data, and Information Service (NESDIS)

used SAR to outline the spill on a near-daily basis. SAR has several advantages over other 

satellite-based sensors. SAR spatial and temporal coverage exceeds that of other sensors because 

many satellites have SAR sensors. Additionally, radar does not require clear skies or daylight.

The post-spill assessment relied on data from Radarsat-1 and -2; TerraSAR-X; CosmoSkyMed-1, 

-2, -3, and -4; ENVISAT; ALOS-1; ERS-2; and PALSAR (Leifer, et al., 2012; Garcia-Pineda, et 

al., 2013a, 2013b, 2015; MacDonald, et al., 2015b)

In contrast to single-band SAR measurements, optical remote sensing provides multi-bands 

and occasionally hyperspectral information on oil slick locations, size, type, and thickness. The 

multi-band data collected by the Moderate Resolution Imaging Spectroradiometer (MODIS) 

instruments onboard NASA’s Terra and Aqua satellites provide a compromise between 

resolution (250–500m) and coverage (2300-km swath), thus can be used to map surface oil over 

broad scales. The proof-of-concept has been demonstrated in the case of the 2003 Lake of 

Maracaibo (Venezuela) oil spill, and further enhanced by another study in 2009, which showed
5

that under sun glint conditions MODIS could have capacity similar to SAR for detection of 

surface oil films. Such a capability enabled researchers during the DwH oil spill to continuously 

monitor the spill’s surface expression, which was further used by numerical modelers to trace the 

oil transport. However, no attempt was made to estimate the oil slick thickness or volume before

the WAMOST project, because it posed significant technical challenges.

Hyperspectral measurements such as those from the NASA Airborne Visible/Infrared 

Imaging Spectrometer (AVIRIS) have demonstrated potential for estimating surface oil thickness 

(or volume), although limited in both space and time. During the DwH oil spill, AVIRIS 

collected hyperspectral data from 380 to 2500 nm in 224 spectral bands at a ground resolution of 

about ~7.6–10 m per image pixel. On 17 May 2010, several flight lines were conducted over 

thick, thin, and oil-free water around the DwH platform. In addition, the U.S. Geological Survey 

(USGS) has developed preliminary surface oil volume from each AVIRIS pixel based on 

laboratory-measured relationships between reflectance and oil-water mixing ratio and oil 

volume. However, how to extend the narrow AVIRIS flight lines to cover the entire DwH spilled 

area remained a technical challenge.

During the WAMOST project, optical remote sensing was used to maximize the advantages

of both AVIRIS and MODIS to: 1) characterize oil slick morphology, 2) understand sun glint 

impact on oil-water contrast, 3) improve oil slick delineation, and 4) quantify oil thickness and 

volume at synoptic scales. The analysis was combined with SAR observations for an

independent consistency check. 

2.1. SAR TEXTURE CLASSIFYING NEURAL NETWORK ALGORITHM (TCNNA)

SAR is an effective tool for detecting floating oil and SAR results were crucial for 

verification of WAMOST modeling efforts. Layers of floating oil dampen capillary waves and 

yield a stronger radar return than can be produced by energy scattered by surface roughness on 

an unoiled sea (Brekke & Solberg, 2005; Leifer, et al., 2012). SAR resolution depends on many 

different factors, including the age, band, and type of the sensor; beam mode; and incidence 

angle. Typically the SAR images that were available in 2010 had resolutions ranging from 12.5 

m to 75 m per pixel. However, there are impediments to using SAR images directly for 

delineating oil slicks with sufficient resolution to be most useful for validation. Because SAR 

images are generated from progressive scans that sweep the targeted area as the instrument 

passes overhead, the image intensity varies strongly across images in a direction perpendicular to 

the orbital path. Returns at near-nadir angles are stronger than those from the edges of the image. 

This precludes use of global contrast enhancement to improve the delineation of oil slicks. The 

WAMOST researchers solved this by developing or adapting algorithms based on training sets 

that classified SAR image pixels as including oiled or unoiled ocean surface.

After preliminary processing, SAR images were analyzed for the presence or absence of

floating oil with use of the Texture Classifying Neural Network Algorithm (TCNNA), which has 

been described in previous publications (Garcia-Pineda, et al., 2009; Garcia-Pineda, et al., 2010). 

This semi-automated routine filters a grayscale image with a 25x25-pixel kernel (Figure 2-1), 

testing pixel-by-pixel for edge and shape detections based upon the Leung-Malik filter bank 

(Leung & Malik, 2001). The neural network algorithm interpolates these detections within a 

training set previously compiled by an expert operator through classification of several thousand 

pixels in images known to contain floating oil, which present distinctive linear features, broadest 

near the origin and tapering away down-wind. Training sets were also compiled for the floating 
6

oil targets resulting from the DwH discharge, which were much larger in area and contained a 

greater variety of regular and irregular shapes. Detection of floating oil with the TCNNA was 

repeatedly validated by comparison with sea-level observations. Output of the TCNNA routine 

was a pixel-by-pixel classification with georectification. Using this product, it was possible to 

segment the area of oil-covered water for the entire image (Figure 2-1). However, the backscatter

by oil slicks is sensitive to oil of any thickness and training sets did not discriminate between 

thin or thick oil; therefore, we classified oil detected by the TCNNA as “thin” unless other lines 

of evidence indicated otherwise. A second image processing routine was used to detect “thick” 

oil in SAR images acquired during the DwH spill.

 

Figure 2-1. TCNNA processing of SAR images

Left: Pixel neighborhood classified by standard image processing filters. Right: Inputs to the neural 

network routine that were compared with training set results to classify each pixel as oiled or unoiled 

water.

2.2. SAR OIL EMULSION DETECTION ALGORITHM (OEDA)

Orbital surveillance of the DwH oil spill resulted in SAR images covering areas near the 

Macondo well where multiple lines of evidence verified the presence of thick oil layers. Analysis 

of these data found persistent anomalies within the oil slicks; these anomalies appeared to 

correspond to regions of thick and emulsified oil. Garcia-Pineda and others (Garcia-Pineda, et 

al., 2013a, 2013b) proposed that the regions’ oil emulsions would yield radar returns with values 

intermediate between unoiled water and thinner sheens due to absorption of radar energy by the 

water present in the emulsions. They developed the Oil Emulsion Detection Algorithm (OEDA)

to segment small regions of SAR where these anomalies were present. The method was found to 

be effective for portions of SAR images with adequate signal-to-noise properties; these were 

typically the regions in images where the radar scan angle was in the range of about 20°–50°.

Results allowed scaling of SAR images to two separate thickness classes: Thin oil was 

nominally 1 μm based on literature values (NOAA Hazmat, 2012), whereas thick oil was 

nominally 70 μm based on AVIRIS acquisitions (described below) and on oil volumes reported 

by operations that collected and burned areas of actionable oil during the DwH spill (Garcia-
7

Pineda, et al., 2015; Mabile, 2013). Figure 2-2 shows the relatively small areas of thick oil 

identified by the OEDA routine in a SAR image collected by ENVISAT on 3 June 2010 at 

11:36UTC.

Figure 2-2. Areas of thick oil and oil emulsion delineated by the Oil Emulsion Detection Algorithm 

(OEDA) in a SAR image of the DwH surface oil June 3, 1010

(Garcia-Pineda, et al., 2013a, 2013b)

2.3. OPTICAL REMOTE SENSING: MODIS AND AVIRIS

Before comparing MODIS and AVIRIS to derived synoptic oil distribution maps, several 

operations were performed to fully understand oil slicks and oil-water contrast under sun glint 

and no-glint conditions. These included the following:

Characterizing oil slick morphology (length, width, length-to-width ratio) using AVIRISderived oil slick thickness (volume) maps

These maps were furnished by the USGS and this work provided the basis for determining

why a pixel-wise correlation between MODIS and AVIRIS would not work. This is because oil 

slicks are very small when compared to MODIS pixel size (250–500 m); the median lengths, 
8

widths, and length-width ratios of different oil classes ranged between 22–38 m, 7–11 m, and 

2.5–3.3 m, respectively. Even for the 30 m resolution Landsat, most of the pixels still offer only 

partial oil coverage. The sub-pixel partial coverage for the MODIS pixels is also much smaller. 

These results are published in (Sun, et al., 2016). Figure 2-3 shows an example of the slick 

characteristics from AVIRIS estimates collected during DwH.

Figure 2-3. Examples of AVIRIS products for oil thicknesses

a) sub-region oil thickness products generated by the USGS based on AVIRIS measurements on 17 May 

2010 (the oil thickness classes for this study were ≤50, 50–200, 200–1000 and >1000μm); and b) the 

corresponding AVIRIS RGB true color image (R: 638.2 nm; G: 550.3 nm; B: 462.8 nm). Symbols (Å) mark 

locations where spectra were extracted and analyzed. (Figure adapted from Sun, et al., 2016)

Understanding oil-water contrast under sun glint 

The goal of this part of our work was to understand when the oil-water contrast reverse signs 

(from negative to positive) using radiative transfer and MODIS observations. The results were 

summarized in Lu, et al. (2016). Another goal was to quantify sun glint intensity threshold under 

which oil films could not be observed. This is critical to avoiding false negative detection. The 

work was based on statistics of MODIS-Terra, MODIS-Aqua, and Visible Infrared Imaging 

Radiometer Suite (VIIRS) measurements of thin oil slicks from oil seeps in the western GoM, 

with results summarized in (Sun & Hu, 2016). Figure 2-4 is an example of these results in which 

sun glint intensity thresholds have been identified.
9

Figure 2-4. Statistics of sun glint strengths (LGN) where oil slicks can or cannot be detected from 

MODIS images

(Source: Sun & Hu, 2016)

Extracting oil slicks from MODIS and VIIRS images

Though methods have existed for semi-automatic delineation of oil slicks in SAR imagery, 

until recently no effective method has been available for MODIS or other optical imagery 

because of the interference of clouds, cloud shadow, and sun glint gradient. Based on a series of 

filters and morphological operators, (Wang & Hu, 2015) developed a technique to extract oil 

slicks in noisy MODIS and VIIRS images in a semi-automatic fashion. The technique was then 

used by Sun, et al. (2015) to analyze the oil footprint and trajectory from the Ixtoc-I oil spill in 

the southern GoM on the Campeche shelf. This is the first time that such satellite-based statistics 

have been obtained and provided to the community (Figure 2-5).

Figure 2-5. Oil slicks from the Ixtoc oil spill in 1979

The slicks were delineated from CZCS and Landsat-MSS imagery after eliminating natural slicks. The 

colors represent different days, and the red triangles are the locations of field observations. (Source: Sun 

et al., 2015)
10

The fundamental difficulty in applying multi-band sensors to estimate oil thickness is their 

lack of enough spectral bands and a lack of field measurement of oil thickness. This latter issue 

does not appear to be solvable in the near future. Therefore, AVIRIS-based oil volume estimates 

derived from laboratory-based look-up-tables between hyperspectral reflectance and oil-water 

mixing ratio and oil volume were used to develop MODIS algorithms to estimate oil volume. 

After expending a significant amount of effort to test various methods of band combinations, a 

histogram-matching method has been developed and applied that allows us to scale limited 

AVIRIS observations to more synoptic and frequent MODIS observations. A manuscript 

describing this technique has been submitted (Hu, et al., Submitted 2016) resulting from a 

coordinated team effort. This effort takes into account cloud masking, masking of other artifacts, 

oil homogeneity (patchiness), sun glint impact on oil and water reflectance, coarse pixel 

resolution, and wide spectral bands. The outcome is a series of MODIS oil maps showing both 

the oil footprint and relative thickness (oil volume) for the selected 19 MODIS images collected 

on 18 days between late April and mid-July 2010. In addition to the surface oil footprint, there 

are three types of maps:

 Continuous surface oil volume,

 Classified oil thickness (three classes: sheen, thin, and thick), and

 Probability of different oil thickness

The MODIS data were preprocessed before analysis. This included: 1) A customized 

atmospheric correction that removed the gaseous absorption effect and molecular scattering 

(Rayleigh) effect to derive the Rayleigh-corrected reflectance (Rrc()). The aerosol effect was 

removed using a nearest-neighboring approach; 2) Mapping of all data to a cylindricalequidistance (rectangular) projection; 3) Delineation of the oil-water boundary using gradient 

method to derive the surface oil footprint. Figure 2-6 shows an example of these three different 

oil maps for the case of 17 May 2010. 
11

Figure 2-6. Examples of the three different MODIS oil distribution maps

MODIS oil distribution maps on 17 May 2010 (16:40 GMT), after scaling up AVIRIS observations. Top left: 

RGB image showing oil footprint; Top right: Surface oil volume (liters per MODIS 250 m pixel); Bottom 

left: Surface oil classes (sheen, thin, and thick); Bottom right: Surface oil probability for each class. No 

coverage: MODIS data are not useable; Oil-free water: ocean surface is free of oil; 1 m, 1–50 m, and 

50 m: ocean surface is covered by oil with thickness of < 1 m, 1–50 m, and > 50 m, respectively. 

The probability of each of the five classes is color coded, for example 0.10 represents that 10% of the 

surface is covered by the corresponding class). Results for the other 18 MODIS images are also 

available. (Adapted from Hu, et al., Submitted 2016)
12

2.4. REMOTE SENSING OF NATURAL SEEPS

Natural seeps are broadly distributed across the GoM (De Beukelaer, et al., 2003; 

MacDonald, 2011), with oil discharge rates from individual natural seeps thought to be generally 

<1 m3d

-1

(National Resource Council, 2003). Consequently, understanding of natural seeps is an 

important component of the cumulative impact analyses in NEPA environmental assessment 

documents. These analyses account for the Non-OCS Oil and Gas-related impact-producing 

factors, including from natural events and processes such as natural seeps.

Satellite SAR can readily detect these oil slicks under favorable wind conditions, i.e., speed 

2–10 m s-1(Espedal & Wahl, 1999). In the 1990s and early 2000s, a provisional remote sensing 

inventory indicated that natural oil slicks cover an average aggregate footprint of 850 km2in the 

Gulf at any given time (Mitchell, et al., 1999), from which the annual discharge rate for the 

entire Gulf was estimated to be 0.38–1.0 × 105 m

3

y

-1

 or 2.5–6.9 × 105bbl y-1and natural seeps 

were cited as the predominant source of oil entering the ocean (Mitchell, et al., 1999; National 

Resource Council, 2003). Following the DwH discharge, several authors have noted the 

prevalence of natural seeps in the vicinity of the DwH wreckage (Crooke, et al., 2015), the 

challenge of distinguishing natural oil slicks from anthropogenic releases (Aeppli, et al., 2013; 

Garcia-Pineda, et al., 2015) and the possibility that oil from natural seeps can have a deleterious 

impact on deep sea fauna (Boehm & Carragher, 2012) comparable to that attributed to oil from 

the DwH discharge (White, et al., 2012). 

To quantify oil discharged from natural seeps and in partial fulfillment of BOEM interests, 

(MacDonald, et al., 2015a) reviewed 1,158 prospective images obtained through the Alaska 

Satellite Facility (ASF) and analyzed 254 individual SAR images collected by the Radarsat-1, 

JERS-1, and ERS1 satellites between May 1997 and November 2007. Note that SAR images, as 

they are cataloged in the ASF archive, can overlap with the data acquired during an orbital pass 

of a satellite. Combining frames that overlapped along individual orbital paths yielded a total of 

176 distinct SAR acquisitions, which were then processed using the TCNNA routines described 

in Section 2.1 to determine the oil slick origin (OSO) points for each slick and the area of oilcovered water that it comprised.

The analysis identified 914 distinct seep zones (Figure 2-7), most of which are concentrated 

along the continental slope in the northern GoM. Seeps also exist in other regions of the Gulf 

including the Sigsbee Knolls and the Campeche Knolls (Bryant, et al., 1991) and a cluster of 

sources that occur within the super-giant Cantarell Complex oil field on the Campeche Bank of 

the Yucatan Peninsula (de Miranda, et al., 2004) . Oil slicks or similar features were eliminated 

from the inventory when a no-seep zone was observed within its 10 km grid cell and a one-cell 

buffer surrounding that cell. This filter removed anomalies that were detected in solitary SAR. 

Integration of multiple images predicts the seafloor locations of persistent natural oil seeps. 

Geographically, the largest numbers of seeps are associated with the salt-tectonic stratigraphy of 

the GoM, which is most pronounced in the central and western region of the continental slope 

(McBride, et al., 1998; Watkins, et al., 1978). A progressive curtailment of natural discharges 

was detected by SAR eastward along the Mississippi-Alabama Slope in the direction of the 

Florida Platform, a massive carbonate structure where sediment accumulation has been limited 

during basin history and where salt bodies are absent (Hine, et al., 2003; Pindell & Kennan, 

2007). 
13

Figure 2-7. Natural seep zones in the GoM

Each point is the geographic average of oil slick origins discharged from a discrete geologic seep 

approximately 2 km in diameter. Estimated total active seep zones total 914 locations. Figure from 

MacDonald, et al. (2015a).

For estimating the amount of oil released from natural seeps, the targets identified in SAR 

images were accumulated in a 10x10 km grid and the bootstrapped average oil cover per cell was 

compiled for the entire GoM. Determining the volume of oil required to form a slick depends 

upon approximations of slick thickness that are inherently difficult to verify. The maximum 

average area of floating oil in a 10×10 km grid cell was 15 km2, but the oil-covered area per cell 

was typically much smaller than that. Applying methods previously used to estimate the total 

flux of natural oil to the surface of the GoM (MacDonald, et al., 1993)National Research Council 

Committee on Oil in the Sea 2003), the average thickness of natural oil slicks was scaled at ~0.1 

μm. The maximum mean volume of surface oil from natural oil seeps in a grid cell was 0.015 m3

km-2(0.09 bbl km-2). The combined footprint of all natural oil slicks, based on average oil cover 

per grid cell, is 779 km2(SD 86.24). In research described below, the residence time for natural 

oil slicks was found to be strongly dependent on wind speed, but was conservatively 8–24 hr. 

Taking these time constants as the formation rate for the composite area of natural slicks, the

annual discharge from seeps for the entire Gulf equates to a flux at the sea surface of 2.51–9.43 × 

104m

3

y

-1

(1.58–5.94 105bbl y-1). Similar calculations can be performed for each of the four 
14

Gulf sub-regions (Table 2-1). This result suggests that the rates of natural seepage were 

substantially concentrated in the northwestern region of the Gulf rather than in the northeastern

region. Previous estimates for natural seepage have included only the aggregate northern GoM.

Table 2-1 provides summary statistics for an inventory of oil-covered water detected in 176 

independent SAR images of the GoM. The Gulf is gridded in 10x10 km cells, which were 

divided into four subregions along the lines of 25° north latitude and 90° west longitude. Each 

cell was sampled in multiple SAR images. Standard deviations were calculated from 100,000 

bootstrap replications. Annual discharge volume estimates assume that natural oil slicks of area 

shown and 0.1 μm thickness persist for 8 to 24 hours, with one standard deviation added to or 

subtracted from the upper and lower bounds, respectively.

Table 2-1.

Summary statistics for inventory of oil-covered water

detected in 176 independent SAR image samples of the GoM

GoM Bootstrap Results Estimated Annual Discharge of Natural Oil

Subregion Area

Mean oil 

cover

Fraction 

of total Lower Bound Upper bound

(102km 

grid) (km2) (km2) (%) SD (m3) (barrels) (m3) (barrels)

Northwest 321,000 530 68% 68.45 1.69×1041.06×1056.55×1044.12×105

Northeast 421,000 52.4 7% 28.11 8.87×1025.58×1038.82×1035.54×104

Southwest 501,000 189 25% 43.00 5.33×1033.40×1042.54×1041.65×105

Southeast 230,000 3.35 <1% 10.96 0 0 1.57×1039.86×103

Total Gulf 1,473,000 775 100% 86.24 2.51×1041.59×1059.43×1045.96×105

2.5. CHRONIC, SMALL-SCALE OIL DISCHARGES

Natural seeps are a common source of surface oil in offshore waters. In addition to mapping 

and quantifying natural seep sources on a regional and basin-wide scale, WAMOST investigators 

examined the occurrence of incidental and chronic anthropogenic releases in coastal waters

(Daneshgar Asl, et al., 2016). Results confirm the efficacy of SAR for accurately detecting 

small-scale spills and validate many of the reports compiled by the EPA’s National Response 

Center (NRC). Detailed analysis of slicks produced by an active natural seep in the Green 

Canyon (GC) 600 lease block illustrated how chronic sources of oil differ under varying 

environmental conditions. The remote sensing results for anthropogenic releases and the GC600 

seep are presented in this section. Modeling results from GC600 are detailed under the Task 3 

topics.

Satellite-borne SAR with TCNNA processing was used to obtain more precise estimates of 

the magnitude of the chronic hydrocarbon discharges described in qualitative pollution reports

associated with the production and transportation network of the U.S. coast of the GoM. The 

NRC oil pollution reports were collected and filtered for the period from 2001 to 2012 to 

determine which of the reports coincided with archived SAR images (Figures 2-8 and 2-9). Some 

of the images covered multiple reports and some of the oil discharges described in one report 

could be observed in more than one image. In all, 177 reports could be investigated from 137 

SAR images collected on the corresponding report dates. Further analysis found that slicks 

observed on 66 of these SAR images could be attributed to 67 of the reported incidents. 
15

Objective measurements indicated that the area of these non-point source oil slicks visible in 

SAR images was, on average, significantly larger than what was reported to the NRC.

 

Figure 2-8. Examples of the oil slicks extracted from ENVISAT ASAR images using TCNNA

TNNA indicated in yellow. (a) Slick discharged at the Taylor site, (b) Slick discharged from undetermined 

source, (c) Slick discharged from a platform. Dots show platform locations. (Source: Daneshgar Asl, et al., 

2015)

The sole chronic point source was observed at the former site of the Taylor offshore 

platform, which produced chronic oil slicks that were consistently much larger than other 

anthropogenic discharges. The evidence in SAR images of floating oil discharged from the 

Taylor site was verified by visual inspection from a boat and aerial photography. For some of the 

slicks discharged from the Taylor site, the accuracy of SAR images for detecting oil slick areas 

was validated by comparing SAR results to Landsat 7 (ETM+) and MODIS images. Aerial 

photography, MODIS, and L7 (ETM+) images could complement the use of SAR images. 

However, coarser pixel sizes tend to indicate larger values of area for oil slicks. Surveillance by 

SAR would improve accuracy for estimates of chronic anthropogenic oil pollution particularly 

where continuous discharges are ongoing.
16

Figure 2-9. Composite of TCNNA and manually extracted oil slicks associated with the Taylor site 

in BOEM lease block MC020: August 2010 to January 2012

Gridded outlines show BOEM lease blocks in region. Energy production platforms are visible as radarbright targets. (Source: Daneshgar Asl, et al., 2015)

In contrast to the Taylor site, the vast majority of anthropogenic oil releases appear to be 

short duration discharges of relatively small amounts of oil that occur at or near the ocean 

surface. Tracking such spills and assessing when potential releases might pose a significant 

threat to wildlife and other resources will be aided by a better understanding of how long small 

quantities of oil persist under ambient marine conditions. The WAMOST team used remote 

sensing methods to compile a dataset that could be used in modeling studies for this purpose

(Daneshgar Asl, et al., 2017). This study was conducted at seep zones located in BOEM GC600 

lease block and three adjacent lease blocks (GC557, GC601, and GC644) in the GoM. GC600 

(Figure 2-10) is one of the most prolific natural hydrocarbon seeps discovered and a perennial 

source of surface oil slicks in the GoM (Garcia-Pineda, et al., 2010).
17

Figure 2-10. Location of GC600 on the bathymetric map of the northern GoM

Grid indicates lease block boundaries and designations. Location of all OSO points detected in 41 SAR 

images (green triangles). The solid lines indicate the topography of the research area. (Source:

Daneshgar Asl, et al., 2017)

The location of GC600 was matched against the coverage of the 41 satellite images collected 

from January 2003 to November 2011 to determine whether the oil slicks were visible in 

archived SAR images. After a preliminary screening, TCNNA was used to delineate any 

observed oil slicks in the SAR images. After oil slicks were extracted, the OSO points were 

selected as single points within detected oil slicks. The ArcToolbox Conversion and Spatial 

Statistics tool was used to measure the width at the OSO points, length, and area of the extracted 

oil slicks. Oil slicks were consistently imaged across GC600 and its surrounding region. Oil 

slicks were detected as elongated dark patches in 41 SAR images collected from January 2003 to 

November 2011 (Figure 2-11). The complex variability of slick shape and direction is clearly a 

result of changing winds and currents in the surface layer (MacDonald, et al., 2002). Each oil 

slick represents a snapshot of the oil that reached the surface OSO points and drifted with wind 

and surface currents for an unknown time before the SAR image collection. A wind and current 

history for the hours immediately following each oil slick was compiled for modeling research 

described under Task 3.
18

Figure 2-11. Locations of floating oil slicks extracted by TCNNA from 41 SAR images

In BOEM GC600 lease block and its surrounding region, January 2003 to November 2011. Each color 

represents a different image. (Source: Daneshagar Asl, et al., 2017)

2.6. REMOTE SENSING RESULTS FOR DEEPWATER HORIZON OIL SPILL

Remote sensing of surface oil from the DwH discharge provides a seminal scenario for the 

effective use of technology and for the importance of understanding of wind and advection 

processes during an oil spill of national significance. The geographic source of the DwH 

discharge was essentially constant during the 87 days of flow, but physical details of the release 

points underwent substantial changes as responders gradually regained control of the well. The 

critical shift was amputation of the fallen risers on 2–3 June 2010. Prior to this action, discharges 

were dispersed among several points of failure along the fallen pipes; after, the entire discharge 

escaped from a single point atop the dysfunctional blowout preventer. Although the gross flow 

rate then increased, recapture of oil and treatment with dispersants reduced the net discharge 

until installation of the riser stack on 15 July ended all releases (Federal Interagency Solutions 

Group, 2010; McNutt, et al., 2012). Therefore, the two periods, 20 April–1 June and 2 June–15 

July, offer significantly different conditions, which potentially affect the subsequent distribution 

and fate of the oil. Remote sensing data provided a means for tracking a critical component of 

this discharge: movement of oil across the ocean surface. It is this component of the oil spill

(Figure 2-12) that generated contaminated marine snow (Passow, 2014), injured mesophotic 
19

corals (Etnoyer, et al., 2015; Silva, et al., 2015), and oiled over 2100 km of the Gulf coast

(Nixon, et al., 2016). 

Figure 2-12. Map of surface oil from DwH visible to SAR: Distribution and average volume of 

surface oil (m3km

-2

) from DwH discharge

Gridded at 5x5 km scale, across a cumulative footprint of 149,000 km2, 24 April–3 August 2010. (Source: 

MacDonald, et al. 2015a)

SAR imaging of surface oil commenced on 24 April and continued at high capacity through 

3 August, after which floating oil was no longer detected. Working under the joint auspices of 

the WAMOST project and the National Resource Damage Assessment effort, MacDonald, et al., 

(2015a) analyzed 169 SAR images collected during this period; they used TCNNA routines 

(Garcia-Pineda, et al., 2009) to delineate areas of water covered by thin (~1 um) oil and OEDA 

routines (Garcia-Pineda et al. 2013) to detect much smaller areas of thick (~70 um) oil. 

Interpolation among the images produced a continuous time series of gridded values for 

floating oil and oil emulsion (m3km-2) in 5x5 km cells across the impacted region (MacDonald, 

Garcia Pineda et al. 2015). The surface oil covered a large and dynamically amorphous region 

that was focused over the release point, but was continuously driven into different distribution 

patterns over a 149,000 km2area of the northeastern Gulf under changing wind and current 

effects. Figure 2-12 shows the average values in these cells for the 24 April–3 August interval. 

Analysis of the daily aggregated values shows two prominent features of the surface oil. First, 

the magnitude of oil was highly sensitive to wind speeds; throughout the emergency, surface oil 

that was visible to SAR decreased sharply when winds exceeded about 5 m s-1and then gradually 

increased when winds subsided (Figure 2-13). Second, there was a state change in the geographic 

concentration and distribution of surface oil when the pre- and post-riser removal periods are 

compared. In summary, the total detected volume of oil decreased by 21% after riser removal. 

However, probably due to increased treatments with Corexit, the ocean area over which the 
20

remaining oil was dispersed increased by 49% (Figure 2-13). At face value, this result is 

consistent with the efficacy of response efforts to reduce surface oil by recapture and burning 

operations (Federal Interagency Solutions Group, 2010) and with the subsea application of 

dispersant. This benefit has to be weighed against increased exposure of planktonic larvae and 

pelagic organisms to oil, which can produce deleterious effects to developing fish at very low 

concentrations (Incardona, et al., 2014). 

Figure 2-13. Time series of surface oil from the DwH discharge

Time series of DwH discharge plotted with surface oil and average wind speeds. Release magnitudes 

show best daily estimates of oil escaping from the damaged well. Discharge time-series subtracts the oil 

recovered from the gross release, while treatment further subtracts burned and dispersed by aerial and 

subsea applications of Corexit at maximum efficacy. Response events (a-l) that potentially affected 

spread of surface oil included: a) Macondo Well blowout occurs, b) DwH drillship sinks and release 

begins, c) aerial dispersant application begins, d) containment dome attempt fails; burning surface oil 

begins, e) subsea dispersant campaign begins (5-May), f) flaring of recovered oil begins, g) top-kill 

attempt, h) riser cut from blow-out preventer; direct injection of subsea dispersant begins, i) Hurricane 

Alex makes landfall, j) capping stack closure stops release, k) Tropical Storm Bonnie makes landfall, l) 

well killed by static backfill. (Source: MacDonald, et al., 2015b.)
21

3. SURFACE OIL TRANSPORT AND WEATHERING MODEL AND 

FORCING FIELDS

3.1. SURFACE OIL DRIFT MODEL WITH WEATHERING IN THE GOM

OSCAR is an oil spill model which calculates and records the distribution (as mass and 

concentrations) of contaminants on the water surface, on shorelines, in the water column, and in 

sediments. This data is recorded in three physical dimensions plus time. The model databases 

supply values for water depth, sediment type, ecological habitat, and shoreline type. OSCAR is 

connected to an oil database that supplies physical and chemical parameters required by the 

model.

Processes covered by OSCAR include surface spreading, slick transport, and entrainment

into the water column, evaporation, emulsification and shoreline interactions; together, these 

processes determine trajectory and fate of spilled oil. In the water column, horizontal and vertical 

transport by currents, dissolution, adsorption, settling, and degradation are simulated. The 

varying solubility, volatility, and aquatic toxicity of oil components are accounted for, such that 

concentrations are relevant to quantification of impacts.

3.1.1. Setup of Deepwater Horizon 2010 oil spill simulation in OSCAR

Conditions at the release site of the DwH oil spill changed several times throughout its 

duration, cutting of the broken riser, direct application of dispersants, and collection of oil at the 

wellhead. To account for this, a setup with multiple release sites was used to simulate the oil spill 

in the OSCAR model, all of which were located at -88.36°longitude and 27.738°latitude at a 

depth of 1500m. From 22 April to 3 June, when the riser was cut, we simulated the release of a 

total of 9000 m3of oil through two release sites: one from the broken riser end and the other 

representing kinks in the pipe above the blowout preventer. The number of 9000 m3represents 

an average between the estimated 63k barrels per day (BPD) at the beginning of the release and 

53k BPD when the riser was cut (McNutt, et al., 2011). The 9000 m3 were divided between 2340 

m

3

from the kinks and 6660 m3from the broken riser end, where the fractionation was obtained 

from estimates of the flow rates from these two sites by in situ measurements from Camilli et al. 

(2011). In reality, multiple kinks were present during this time, each with a small, likely noncircular outlet geometry (Camilla et al. 2011). Because the exact geometries are unknown, we 

used 0.1 m as an effective release diameter for the combined kinks, which has the effect of 

generating small droplets as may be expected from a release from several small orifices. For the 

broken riser, we used a release diameter of 0.53 m as estimated by Camilli et al. (2011) from in 

situ measurements (Camilli, et al., 2011). After the broken riser was cut on 3 June, a single 

release site was formed and large quantities of oil started being collected at the source. It is 

estimated that approximately 20k BPD were collected at the wellhead from 3 June onward

(Federal Interagency Solutions Group, 2010). Subsurface dispersants had already been applied, 

but could now be focused to a single source. To accommodate the collection of oil, we deducted

200 BPD from the previous total and used a single discharge point releasing 5820 m3per day 

from 3 June onward until the well was sealed on 15 July. To accommodate the application of 

dispersants, we made use of a Weber-scaled model for droplet generation (Johansen, et al., 

2013), where the effect of dispersants on the droplet size distribution is estimated by reducing the 

oil-water interfacial tension by a factor of 100 from the estimated baseline value of 0.03 N/m. 

The composition and weathering characteristics of the oil used in the simulation is based on 
22

weathering studies of the Macondo crude (Daling, et al., 2014). The oil spill was simulated 

through 20 August, to follow the fate of the oil after the well had been capped.

The currents used in the simulation were from the HYCOM GoM 1/25° (GOMl0.04) model

(similar model configuration is described in (Dukhovskoy, et al., 2015b; Zamudio & Hogan, 

2008) obtained from the HYCOM model (HYCOM GoM experiment 31.0, n.d.). The wind data 

was obtained with the COAMPS model (Chen, et al., 2003). The model setup is summarized 

above in Table 3-1.

Table 3-1.

Key setup information for simulation of the DwH oil spill with the OSCAR model

Input Value & Description Source

Well coordinate 88° 15.6’W, 27° 44.28’N

Release depth 1500 m

Release rate total

22.04.2010–03.06.2010

9000 m3/day McNutt et al. (2011)

Release rate kinks

22.04.2010–03.06.2010

2340 m3/day Camilli et al. (2011)

Release rate broken riser

22.04.2010–03.06.2010

6660 m3/day Camilli et al. (2011)

Wellhead release diameter 0.54 m Camilli et al. (2011)

Kink release diameter 0.1 m Estimate

Release rate cut riser

04.06.2010–15.07.2010

5820 m3/day Federal Interagency Solutions Group

(2010)

Ocean currents HYCOM GoM 1/24oexp. 31.0 HYCOM (GOMl0.04)

Winds COAMPS Model U.S. Naval Research Laboratory, Marine 

Meteorology Division, Monterey, CA

3.1.2. Adaptation of the OSCAR model for the WAMOST Project

The WAMOST project required that we present gridded data from OSCAR in the netCDF 

format. This necessitated rewriting OSCAR's code to expose internally calculated arrays and 

write them to disk. This work was performed in two steps: (1) we changed variables shared with 

so-called common blocks to be represented in modules; and (2) we adapted the code to use the 

netCDF C library and wrote transformations from the internal grid structure to the grid structure 

needed for the netCDF format. This allowed for data sharing from OSCAR's internal calculations 

that was previously not possible. From this, we were able to output surface grids showing the 

oil's mass, volume, viscosity, and lifetime on the surface. Part of this work has been incorporated 

into the latest version of OSCAR: 8.0. OSCAR is free for academic use, which means that the 

above mentioned changes will enhance the use of OSCAR in research settings.

Results of the simulation indicate the changes in output and dispersion of oil from the 

Macondo well that took place as the responders increased direct recovery of oil, the application 

of dispersants, and burning of surface oil. Figure 3.1 focuses on the simulated effects of 

dispersant application. Consistent with SAR observations (Figure 2-12), the simulation indicates 

that dispersants were correlated with a decrease in volume and mass of surface oil, but an 

increase in the total footprint area. The gridded dataset produced by the modified OSCAR code 

can support further investigation of the dispersion and fate of the surface oil during DwH.
23

Figure 3-1. OSCAR simulation of surface oil during DwH discharge

Left: Oil mass, area and volume for oil discharge without subsea dispersant treatment. Right: Oil mass, 

area and volume for oil discharge with subsea dispersant treatment.

3.2. FORCING FIELDS AND VALIDATION AND ADJUSTMENT OF THE OIL DRIFT 

MODEL

Surface oil drift models are forced by surface ocean currents and near-surface winds. In most 

oil hindcasts, oceanic and atmospheric forcing fields are derived from reanalysis datasets. These 

datasets combine observational data and numerical models, providing dynamically consistent 

gap-free fields of oceanic and atmospheric variables. Requirements for the ocean models that 

would provide forcing fields for the DwH hindcast were outlined in the call for proposals. 

According to the RFP, the hydrodynamic model should be an existing and previously validated 

data-assimilating ocean model for the GoM, with demonstrated veracity to provide accurate 

analyses of surface oil transport. At the time of the proposal award (2012), there were two 

candidate models that satisfied the specified requirements: the HYCOM (Bleck, 2002; 

Chassignet, et al., 2003; Chassignet, et al., 2006) GoM 1/25° (GOMl0.04) analysis, and the 

Naval Oceanographic Office (NAVOCEANO) operational ocean prediction system for the GoM 

and Caribbean (NCOM AmSeas) (Barron, et al., 2004, 2006)

3.2.1. Ocean Forcing Fields

Two data-assimilative models the HYCOM GOMl0.04 and NCOM AmSeas were compared for 

the following criteria:

 Time coverage, horizontal and vertical grids and resolution: Special attention was given 

to the near-surface vertical resolution that is important for accurate representation of the 

near-surface velocity shear. 

 Sea Surface Height fields: Both models demonstrated general agreement with the 

altimeter-based data. Some discrepancies in locations (the northwestern extent of the 

Loop Current (LC), positioning of the Loop Current eddies (LCE) and frontal cyclones) 

and shapes of the mesoscale features. These discrepancies, however, may play a crucial 

role in oil transport resulting in distinctly different oil forecasts.
24

 Loop Current contours: Representation of the mesoscale features in the HYCOM 

GOMl0.04 and NCOM AmSeas was assessed by comparing the LC front in the 

simulations to the CCAR altimeter-based data. In general, both models have similar skills 

in predicting the LC front.

 Near-surface velocity fields: Discrepancies in the near-surface velocity fields in the 

HYCOM analyses were demonstrated. The discrepancies are related to differences in 

positioning of the LC and cyclonic eddies. 

 SST and SSS were compared to the ship observations collected in the area during the 

DwH spill event. The models demonstrated poor agreement with the SSS from ship 

observations. Both HYCOM GOMl0.04 and NCOM AmSeas SST fields matched 

reasonably well with the Shipboard Automated Meteorological and Oceanographic 

System (SAMOS) ship data (SAMOS, n.d.). Vertical profiles were compared to ARGO T 

and S profiles. There was no obvious evidence for either model to demonstrate a better 

skill. The largest deviations of the model T and S profiles from ARGO profiles were 

related to cases when the drifter was close to a mesoscale feature (LC, LCE, or a 

cyclone).

 Representation of the river plume in the models: The relationship between SSS and 

concentrations of optically active substances near the Mississippi River Delta was 

utilized to evaluate representation of the Mississippi river plume and near-coastal SSS in 

the analysis fields. A satellite ocean color index (OCI) was used to validate SSS field in 

the NCOM AmSeas and HYCOM GOMl0.04 near the Mississippi River Delta during the 

DwH oil spill event of 2010. The analysis indicated that, compared to NCOM, HYCOM 

tends to have lower OCI (higher SSS in GoM) for 30 ≤SSS ≤35 except for 30.5, 31, and 

34.5. This result demonstrated that HYCOM has a more dispersed river plume and its low 

salinity water spreads farther offshore than in the NCOM forecast.

3.2.1.1. HYCOM GOMl0.04 and NCOM AmSeas

HYCOM analyses. The GoM model has 1/25° equatorial resolution and latitudinal 

resolution of 1/25° cos(lat) or ~3.5 km for each variable at mid-latitudes. It has 20 coordinate 

surfaces in the vertical (for more detail see (Dukhovskoy, et al., 2015b; Zamudio & Hogan, 

2008). The data assimilation is performed using the Navy Coupled Ocean Data Assimilation 

(NCODA) (Cummings, 2005) system with a model forecast as the first guess. NCODA 

assimilates available satellite altimeter observations (along tracks obtained via the 

NAVOCEANO Altimeter Data Fusion Center), satellite and in situ sea surface temperature 

(SST) as well as available in situ vertical temperature and salinity profiles from XBTs, ARGO 

floats and moored buoys using a 3D-VAR technique (Metzger, et al., 2014). The system is run in 

real time at the NAVOCEANO Major Shared Resource Center. The HYCOM NCODA GoM

predictions were the major source of oceanographic information for the oil drift forecasts 

performed by several research groups during the DwH oil spill.

HYCOM GOMl0.04 (hereinafter, referenced as HYCOM) has several analysis datasets 

(experiments) that partially or completely cover the time period of the DwH oil spill (Table 3-2). 

Experiment 20.1 was the only one available at the time of the DwH spill. Experiment 31.0 covers 

the entire DwH event, and thus was chosen as the primary HYCOM analysis dataset for 

providing oceanographic fields used in this study. HYCOM GoM Reanalysis (50.1) became 

available near the end of the WAMOST project; therefore, the reanalysis fields have not been 
25

thoroughly validated. It should be mentioned, however, that the reanalysis is nested in global 

HYCOM reanalysis (in contrast to the other HYCOM experiments that used climatology open 

boundary conditions), which presumably has better representation of the high-frequency 

variability propagating from the open boundaries into the interior of the GoM domain. This 

potentially can improve representation of the LC behavior and representation of the LC fronts.

The NAVOCEANO operational ocean prediction system for the GoM and Caribbean

(AmSeas) is based on the US Navy’s Coastal Ocean Model (NCOM) configured for the GoM 

and Caribbean Sea (Barron, et al., 2004, 2006; Ko, et al., 2003). The model’s computational grid 

has horizontal spacing of 1/36° (~3km) and 40 levels in the vertical. The ocean prediction system 

assimilates all quality controlled observations in the region including satellite SST and altimetry, 

surface and profile temperature and salinity data using the NCODA MVOI system (for 

experiments 20.1), MODAS MVOI (30.1), MODAS 3D-VAR (31.0 and 50.1). Boundary

conditions are applied from the NAVOCEANO operational 1/8°global NCOM. It should be 

noted that NCOM AmSeas (hereinafter, referenced as NCOM) starts on 8 May 2010 and thus 

does not cover completely the time period of the DwH event. 

Table 3-2.

Ocean analysis data sets

GoM Analysis, 

Reanalysis, 

Experiment

Horizontalvertical 

resolution

Period 

Covered

Archive 

Frequency, 

hours

Atmospheric 

Forcing

OB Forcing

HYCOM 

GOMl0.04, 20.1

Equatorial 

1/25° and 

latitudinal of 

1/25°cos(lat), 

~3.5 km

20 hybrid 

layers

01/2003–

07/2010

3 NOGAPS Climatology from 

1/12° HYCOM 

Atlantic Ocean

HYCOM 

GOMl0.04, 30.1

06/2010–

03/2012

3 NOGAPS Climatology from 

1/12° HYCOM 

Atlantic Ocean

HYCOM 

GOMl0.04, 31.0

04/2009–

07/2014

1 NOGAPS 1/12° global 

HYCOM 

reanalysis

HYCOM 

GOMl0.04 50.1

1/25° 

(~3.5km) 

36 hybrid 

layers

01/1993–

12/2012

3 NCEP CFSR 1/12° global 

HYCOM 

reanalysis

NCOM AmSeas 1/36° (~3km)

40 σ/z levels

05/08/2010–

present

3 COAMPS NAVOCEANO 

operational 1/8° 

global NCOM

3.2.1b. Vertical Grids in HYCOM and NCOM

Though the models have comparable horizontal resolutions, similar assimilation techniques 

and assimilated data, the vertical grids are different and lead to a different accuracy in 

representation of the near-surface layer in the models. This is of particular importance for the 

surface oil drift simulation because oil is driven by the ocean surface currents. When ocean 

simulation is used as a forcing field, ocean surface currents are approximated by the ocean 

currents in the upper model layer (see Section 3.2). The coordinate surfaces in HYCOM and 

NCOM are shown in the upper 15 m of the water column (Figure 3-2). HYCOM model output 

fields used in this study as forcing fields for oil models have been interpolated into the 

geopotential depths (shown with the blue circles in the figure). Thus, surface currents from 

HYCOM GOMl0.04 experiments 31.0 and 20.1 are approximated by velocities in the topmost 
26

layer that is 3 m thick. For ease of analysis and visualization by academic community and other 

users, HYCOM format archive fields are converted to NetCDF files that contain diagnostic fields 

at standard Levitus depth levels. Thus the uppermost layer in the archived fields has coarser 

vertical resolution than the native HYCOM grid as demonstrated in Figure 3-2.

Figure 3-2. Coordinate surfaces in the HYCOM and NCOM in the upper 15 m

The blue circles indicate interpolated depths for model archive. The red crosses are locations where 

velocities and tracers are simulated.

3.2.1c. Sea Surface Height fields in HYCOM and NCOM

Representation of the mesoscale ocean circulation in the hydrodynamic models controls oil 

transport by ocean currents, especially away from coastal regions. In the considered models, the 

accuracy of the mesoscale circulation is largely constrained by altimeter-based observations of 

the sea surface height (SSH). During 2010, SSH fields were used from the Jason-1-Jason-2 

tandem mission and Envisat. Representation of the SSH fields in the models is compared against 

the Colorado Center for Astrodynamics Research (CCAR) GoM gridded altimeter-based 

products for the time period of DwH event. 
27

Figure 3-3. Demeaned SSH fields (m) from CCAR, NCOM and HYCOM time-averaged over May 12–

19, 2010

The black lines are 0.17-m SSH contours used for identification of the LC and LCE fronts (Leben, 2005). 

The numbers in CCAR map identify mesoscale features discussed in the text.

Both models demonstrate general agreement with the altimeter-based data. There is some 

disagreement (~1 day) between the models on timing of the Loop Current Eddy (LCE) 

detachment-reattachment. Also, there are discrepancies in locations (the northwestern extent of 

the LC, positioning of the LCE and frontal cyclones) and shapes of the mesoscale features. These 

discrepancies, however, may play a crucial role in oil transport resulting in distinctly different oil 

forecasts. As an illustration, time-averaged SSH fields are shown in Figure 3-3. Both NCOM and 

HYCOM reproduce major features of the SSH field during this time, such as location and 

northwestern extent of the LC, remnants of a LCE in the western GoM (feature “5”), frontal 

cyclone formed in the “neck” of the LC (feature “2”). However, the following disagreements

exist: the shape and position of the forming eddy (feature “1”) is different both in HYCOM and 

NCOM; HYCOM does not simulate a cyclone to the north of the LCE (feature “3”) that is 

observed in the CCAR field; and simulated anticyclonic and cyclonic eddies in the western GoM 

(“4” and “5”) in HYCOM and NCOM are different from CCAR. Although both models 

reproduce the frontal cyclone “2”, its location and strength in the NCOM is different from the 

cyclone in CCAR. This results in a different shape of the LC and the width of the LC “neck” 

which in its turn controls LC eddy shedding process. In NCOM, the LC is in a stable formation 

whereas in HYCOM the LC is in the process of eddy shedding. Compared to the models, the LC 

in CCAR is somewhere between the two simulated LC states. 
28

3.2.1d. Comparison of the Loop Current contours

Representation of the mesoscale features in the HYCOM and NCOM is assessed by 

comparing the LC front in the simulations to the CCAR altimeter-based data. Following Leben 

(2005), the LC and LCE fronts are tracked using the 0.17 m contour in demeaned SSH fields. 

Demeaned fields are calculated by subtracting the spatial mean from each daily SSH field; this is 

necessary to remove bias in the surface elevation fields associated with different reference 

surfaces and seasonal height variations resulting from upper-ocean warming and cooling. The 

maximum northern extent of the LC is an important characteristic of the LC system for oil drift 

modeling during the DwH accident because the proximity of the strong currents associated with 

the LC or LCE to the oil spill determines oil transport and its spatial distribution. Figure 3-4 

shows the maximum northern extent of the LC system. This demonstrates the northernmost 

position of the front of the LC or LCE west of the 91 W longitude. The diagram demonstrates 

good agreement between the simulated and altimetry-derived northern extent of the LC. 

However, the LC front in the simulations tend to be located farther north when the LC is in the 

extended position. By contrast, the LC front stays in the more southward position when the LC 

retracts. Dukhovskoy, et al. (2015b) concluded that northern extents of the LC fronts derived 

from altimeter-based SSH fields have slight southward biases and is likely due to the along track 

detrending and smoothing applied to the altimeter data during the objective analysis. The study 

also found the tendency of the HYCOM LC to retract farther south than the altimetry-derived LC 

fronts. It is worth noting that time series of the CCAR LC northern extent in Figure 3-4 is 

smoother compared to the simulations, which is related to the coarser resolution of the altimetrybased SSH product.

Figure 3-4. Maximum northern extent of the LC estimated from the analysis and CCAR SSH fields

Data are shown for the time period from May 9 through September 15, 2010.

A Modified Hausdorff Distance (MHD, see section 3-4a for the description of the 

methodology) is calculated between the LC front derived from the model SSH fields and the LC 

front derived from CCAR data (Figure 3-5b). In general, both models have similar skills in 

predicting the LC front. Individual peaks (corresponding to high MHD scores) are related to 

episodic cases when simulated LC front deviates from the altimeter-based LC front as illustrated 

in Figure 3-5a. During this time instance, the NCOM predicted LC position shifted eastward 
29

relative both to CCAR and HYCOM. The differences (with negative values indicating better 

HYCOM skill and positive for better NCOM skill) between HYCOM and NCOM MHD scores 

are only slightly negative (-0.03) demonstrating fairly similar skills for both models.

Figure 3-5. Comparison of the LC fronts from the NCOM and HYCOM vs CCAR

(a) SSH fields on August 22, 2010. (b) MHD for the LC fronts from HYCOM GOMl0.04 and CCAR (blue) 

and NCOM and CCAR (red). (c) Histogram of differences between MHDs for HYCOM and NCOM. 

Negative values indicate better scores for HYCOM, positive for NCOM.

Importance of accurate representation of the mesoscale features in the ocean models is 

illustrated in Figure 3-6. In general, the near-surface circulation from three HYCOM analysisreanalysis data sets on May 8, 2010 00Z looks similar. However, discrepancies in smaller details 

are evident. In particular, the direction and strength of the surface currents to the southeast of the 

DwH location (marked by the yellow dashed circle) are different among the models. The flow 

structure in this area is primarily controlled by the LC (“3” in Figure 3-6a) and two anticyclonic 

eddies (“1” and “2” in Figure 3-6a) whose location and strength is different in the simulations. 

The correct flow structure in this particular area on this day is important for oil simulation. 

Figure 3-6d illustrates the formation of the “tiger tail” when some fraction of the surface oil 

became entrained into the LC after being advected southeastward from the DwH location. While 

the flow structure in the analysis 20.1 (Figure 3-6a) would advect surface oil southeast allowing 

oil to be entrained into the LC, the other two forecast fields would advect the surface oil 

primarily south or southwest. Note also the difference in the flow speeds in this area. The 

currents are noticeably stronger in the analysis 20.1 compared to the other two velocity fields. 
30

This would explain rapid development of the “tiger tail” observed in the satellite-based surface 

oil concentration.

Figure 3-6. Near-surface currents from HYCOM GoM analyses 20.1 (a) and 31.0 (b) and reanalysis 

50.1 (c) on May 8, 2010, 00Z

The green bullet designates location of the DwH. Red “1” and “2” mark anticyclonic eddies in HYCOM 

20.1 (a). The yellow dashed circle highlights the area where the “tiger tail” starts forming before oil gets

entrained into the LC as shown in the SAR-TCNNA image (d). (d): SAR-TCNNA surface oil volume 

(m3/km2) on May 8, 2010, 00Z.

3.2.1e. T-S fields

Sea surface temperature (SST) and sea surface salinity (SSS) fields from the analysis data 

sets are compared to the ship observations collected in the area during the DwH event. The data 

were provided by SAMOS. Quality-controlled data from six ships were used for the comparison 

(Figure 3-7a). The analysis was performed for the locations within the region limited by the 

94W and 82W longitudes and 24N and 31N latitudes. For comparison analysis, every ship 

observation was first collocated in time with NCOM and HYCOM fields. Then at every time 

level, ship observations were grouped by the model grid cells based on the shortest-distance 

criterion. Within the grid cells, ship observations were averaged and compared to the 

corresponding model SSS and SST.

Scatter plots of SSS and SST from the models and two cruises (WTER and WTEO) are 

shown in Figures 3-8 and 3-9. SSS fields from HYCOM and NCOM reveal poor agreement with 

ship observations especially for the low-salinity water (Figure 3-8). The correspondence between 

the models and observed SSS is better for higher salinities (>32). There is barely any correlation 

between the models and the SAMOS ship SSS for salinities less than 30. Both models tend to 
31

have strong positive bias for the low-salinity regions. The relationship between simulated SSS 

and WTER observations is better than between the models and SSS measurements collected 

along WTEO track. The reason is that WTEO track was predominantly in the regions influenced 

by the Mississippi River, whereas WTER track was most of the time at more distant locations 

from the Mississippi Delta and majority of observed SSS was in high-salinity regions unaffected 

by land runoff. A detailed analysis of discrepancies in SSS between HYCOM simulations and 

SAMOS data is given in Androulidakis & Kourafalou (2013) and Kourafalou & Androulidakis 

(2013).

For SST fields HYCOM and NCOM reveal different agreement with the SAMOS ship data 

(Figure 3-9). Whereas NCOM shows similar correlation with WTER and WTEO SST, HYCOM 

has surprisingly weak agreement with WTER SST and nearly perfect correlation with WTEO 

observations. 

Figure 3-7. SSS ship observations used for NCOM and HYCOM validation

(a) SSS measurements along six ship tracks (Gordon Gunter, WTEO; Henry B Bigelow, WTDF; Oregon 

II, WTDO; PISCES, WTDL; Oceanus, WXAQ; Nancy Foster, WTER) during the oil spill event in the GoM

in 2010. Models’ comparison against SSS measurements from WTER (b) and WTEO (c) are presented in 

Figure 3-8. 
32

Figure 3-8. Scatter plot of SSS data from the models vs ship observations

NCOM vs ship WTER (a) and WTEO (c) observations. HYCOM vs ship WTER (b) and (d) observations.

The red line is a linear regression between the model data (Sm) and observations (So). The estimates of 

the regression coefficients and statistics (the 0.95 Confidence Interval of the slope and the coefficient of 

determination R

2

) are listed in the figure. 
33

Figure 3-9. Scatter plot of SST data from the models vs ship observations (WTER data only)

(a) NCOM vs ship observations. (b) HYCOM vs ship observations. The red line is a linear regression 

between the model data (Tm) and observations (To). The estimates of the regression coefficients and 

statistics (the 0.95 Confidence Interval of the slope and the coefficient of determination R

2

) are listed in 

the figure.

Vertical profiles of T and S from the model analysis data sets are compared with the T-S 

profiles derived from ARGO floats that drifted in the study area during the DwH event (Table 3-

3).
34

Table 3-3.

Time and coordinates of ARGO float measurements used in the validation analysis

Date 

[MMDD]

Time 

[hhmm]

Lat (N) Lon (W)

05/24 0817 28.931 88.465

05/26 1653 26.407 86.223

05/30 0720 26.912 87.163

06/07 1720 27.294 85.492

06/17 0818 27.483 86.725

06/20 0412 27.496 87.927

06/29 0855 28.493 88.755

07/04 0242 27.989 89.813

07/23 0807 28.947 88.735

08/16 1024 28.639 88.296

08/26 0854 28.506 87.257

09/07 1808 26.58 89.124

ARGO floats data are available from the National Oceanographic Data Center (NODC) for 

the DwH oil spill event (NODC, 2013). Overall, the vertical T and S profiles in both model 

analyses match the ARGO floats. There is no obvious evidence for either model to demonstrate a 

better skill. The largest deviations of the model T and S profiles from ARGO profiles are related 

to cases when the drifter was close to a mesoscale feature (LC, LCE or a cyclone). Then even a 

slight deviation in the forecasted eddy location from its actual positioning would result in 

substantially different thermohaline structure of the water column at this point as exemplified in 

Figure 3-10 (b and d). Both T and S profiles measured by ARGO floats reveal notable 

differences from the NCOM and HYCOM analyses. ARGO salinity profile reveals a typical 

structure of the LC water with somewhat fresher near-surface water mass that reflects the 

influence of local diluation in the Caribbean Sea and remote influence of the Amazon River (

(Gibbs, 1970) Hu et al., 2004). Below the fresher layer is the subtropical underwater (labeled 

“STUW” and also referenced as Tropical Water) derived from the subtropical surface waters 

north and south of the equator (Vidal, et al., 1994). This suggests that the floats in cases shown in 

Figure 3-10 (b and d) were drifting inside the LC. Indeed, altimeter-derived SSH demonstrates 

that ARGO floats were drifting along the LC front on these days (Figure 3-11). However, in the 

models the drifter location was close yet still outside the LC. 
35

Figure 3-10. Vertical T and S profiles from ARGO buoy measurements and NCOM and HYCOM

For four dates (a–d). Locations of the buoys and corresponding SSH fields is shown in Figure 3-11. 

“STUW” is the abbreviation for subtropical underwater.
36

Figure 3-11. CCAR SSH fields on selected dates

Black bullets designate locations of ARGO floats whose T-S profiles are presented in Figure 3-10.

3.2.1f. Representation of the river plume 

Coastal waters in this area are characterized by high concentrations of optically active 

substances such as chlorophyll-a, colored dissolved organic matter (CDOM) and sediments. 

These substances are useful for detecting and tracking river plumes because of strong 

relationship between the concentrations of optically active substances and river runoff (Hu, et al., 

2003; Kudela & Palacios, 2007; Miller & McKee, 2004; Son & Gardner, 2011). Satellite 

observations still cannot provide sea surface salinity measurements of sufficient resolution and 

precision for accurate detecting river plumes (Lagerloef, 2012). Alternatively, the unique optical 

properties of CDOM, sediments and chlorophyll-a are applicable as a proxy for low-salinity 

water near shore. The optical properties associated with CDOM and sediments loading act as 

quasi-conservative tracers. Their depletion is primarily caused by mixing and to a lesser extent 

by gravity and biologically mediated processes and photo-oxidation. Ocean optical properties 

related to the presence of chlorophyll vary spatially and temporally depending on nutrient 

availability, light and temperature. Thus, CDOM and sediment can be used as tracers of the river 

plume, chlorophyll concentration is a practical marker of a distinct water mass. In the Mississippi 

River Delta region, there exists a relationship between salinity and both chlorophyll-a 

concentration and CDOM absorption. Specifically, surface chlorophyll-a concentration (0.1 to 17 
37

mg m

-3

) and CDOM absorption (0.07 to 1.1 m-1 at 412 nm) have been found to generally 

decrease with increasing salinity waters (D'Sa & Miller, 2003). Wetzel (2001) and Lohrenz, et al.

(1997) have also closely related salinity variability to variability in chlorophyll-a pigment in the 

northern GoM such that as salinity increases Chlorophyll-a concentration decreases. In addition, 

CDOM has been found to be a valuable parameter for tracking advective freshwater movements 

along the Louisiana shelf using remote sensing (e.g., Hu, et al., 2004; D'Sa & Korobkin, 2008; 

Green, et al., 2008). 

The relationship between SSS and concentrations of optically active substances near the 

Mississippi River Delta is used to evaluate representation of the Mississippi river plume and 

near-coastal SSS in the analysis fields. A satellite ocean color index (OCI) is used to validate 

SSS field in the NCOM and HYCOM near the Mississippi River Delta during the DwH oil spill 

event of 2010 (Figure 3-12). OCI is derived from chlorophyll-a concentration (Chla) that 

removes contamination from sun glint and thin clouds, and also naturally includes the effect of 

all other colored material in the coastal zone (Hu, et al., 2011). Visual inspection of the sevenday average fields reveals a better agreement between NCOM and OCI. Specifically, there is a 

notably broader region of low salinities in HYCOM compared both to NCOM and OCI (see 

more detail on river parameterization in HYCOM in Kourafalou & Androulidakis (2013).

Modified Hausdorff Distance (MHD) is employed for quantitative comparison of spatial 

distribution of the forecasted salinity fields (Figure 3-13). MHD is calculated between OCI 

contour of a randomly selected value of 0.5 and contours of simulated SSS drawn at 0.5 interval 

between 30 and 35. Both curves have a broad minimum (indicating the best correspondence 

between SSS and OCI of 0.5) within 32 and 33. Within this range, the MHD scores for both 

analyses are similar demonstrating same skills of the models in forecasting SSS over this 

interval. For the SSS values outside this interval, NCOM has notably smaller MHD than 

HYCOM. It should be noted that the true relationship between OCI and SSS for this region is not 

known. Smaller MHD indicates better similarity in shape between the SSS and OCI contours, 

however the value of the SSS contour, that has the best MHD score, may be biased from the true 

SSS value. To better highlight discrepancies in the simulations, the SSS - OCI correspondence 

for NCOM and HYCOM is determined (Figure 3-13). The relationship is established by seeking 

OCI that provides the minimum MHD for a given SSS. The analysis indicates that compared to 

NCOM, HYCOM tends to have lower OCI (higher SSS in GoM) for 30 ≤SSS ≤35 except for 

30.5, 31, and 34.5. This result demonstrates that HYCOM has a more dispersed river plume and 

its low salinity water spreads farther offshore than in the NCOM forecast (more details and 

examples of MHD application for comparison of river plumes in hydrodynamic models and OCI 

are in Hiester, et al. (2016).
38

Figure 3-12. Ocean Color Index (left column) and SSS fields from NCOM (middle) and HYCOM 

Right: For three seven-day time periods (rows). Middle and right: The black contours delineate OCI=0.5 

(left) and SSS=32.
39

Figure 3-13. Mean MHD for OCI–HYCOM and OCI-NCOM

The bullets (right vertical axis) indicate OCI value that provides the best (smallest) MHD score for the 

model salinity contours. The MHD scores have been averaged for the time period from May 24 through 

September 9, 2010.

3.2.2. Wind Forcing Fields

The surface oil slick is highly variable in space and time due to the forcing from near surface 

ocean currents and forcing due to the surface winds. In order to predict the surface oil transport, 

it is essential to have an accurate representation of the surface wind field. A number of reanalysis 

fields exist that could be used as a source for wind forcing fields (Reanalysis.org, 2016), yet 

ideally the near-surface ocean currents should be dynamically consistent with the wind field. 

Otherwise wind forcing may counteract the ocean forcing resulting in errors of the surface oil 

transport forecasts. If ocean surface currents are from a numerical simulation, such consistency is 

provided when atmospheric fields forcing the ocean model are used to derive the near-surface 

winds. Thus, in this study only wind fields that forced HYCOM and NCOM analyses were 

considered (Table 3-4). The CCMP wind data were not used as forcing fields in the ocean 

forecasts, however the data set provides accurate information about ocean vector surface winds

combining scatterometer and radiometer winds. Hence, CCMP fields can be employed for 

validation or comparison for other wind data. Besides in HYCOM simulations, NOGAPS and 

Climate Forecast System Reanalysis (CFSR) winds were bias-corrected towards scatterometer 

winds. Thus, CCMP winds can be used as forcing fields with HYCOM near-surface ocean 

currents in the discussed oil drift numerical experimentations. NOGAPS, CFSR, and CCMP 

wind fields have been validated in the previous studies (Atlas, et al., 2011; Sharp, et al., 2015; 

Stopa & Cheung, 2014; Wallcraft, et al., 2009). COAMPS winds have been extensively used in 

other studies (e.g., Hodur, et al., 2002; Meir, et al., 2013; Pullen, et al., 2011). A detailed 

validation of COAMPS winds has been performed by Fitzpatrick, et al. (2010). COAMPS and 

CCMP winds were the primary forcing fields for oil simulations performed within this study. 

Thus, a more detailed comparison of these wind fields is provided below with particular focus on 

the accuracy of representation of wind vectors and atmospheric fronts in the study region. 
40

Table 3-4.

Wind data sets

Wind Data Temporal 

Resolution

Spatial 

Resolution

Spatial Coverage

Coupled Ocean-Atmosphere Mesoscale 

Prediction System (COAMPS)

3 hr 0.2° 120°W–60°W, 0°N–

32°N

NCEP Climate Forecast System Reanalysis 

(CFSR)

1 hr 0.25° Global

Navy Operational Global Atmospheric Prediction 

System (NOGAPS)

3 hr 0.5° Global

Cross-Calibrated Multi-Platform Ocean Surface 

Wind Vectors (CCMP)

6 hr 0.25° 0°W–360°W, 

78.7°S–78.7°N

3.2.2a. COAMPS and CCMP: comparison of wind vectors and wind speed 

Here, the COAMPS (Coupled Ocean/Atmosphere Mesoscale Prediction System) wind data 

are compared to the CCMP (Cross-Calibrated Multi-Platform) wind data. To focus on the region 

of interest, the grids for both data sets were trimmed to 100°W to 80°W longitude and 18°N to 

31°N latitude. Also, for both data sets to match temporally, only the 00Z, 06Z, 12Z and 18Z 

points were used. Finally, a bilinear interpolation of the COAMPS data was implemented in 

order to have the grid points of each data set collocated.

The Root Mean Square (RMS) difference, wind speed bias, and vector difference plots are 

averaged over the time period from April 20, 2010 00Z through July 16, 2010 18Z. This time 

interval produced 352 individual fields for the u- and v-components of the wind, wind speed, and 

wind vectors. In addition, three frontal passage events were plotted to determine how well the 

fronts were represented in both the CCMP and COAMPS data sets.
41

Figure 3-14. RMS difference between CCMP and COAMPS

(a) u-wind component; (b) v-wind component; (c) wind speed. Note a positive bias in v-component and 

speed over the DwH spill area.

There are fairly large differences between CCMP and COAMPS seen in the RMS difference 

fields (Figure 3-14). The u-component RMSD field has values between 1–2 m/s throughout most 

of the GoM. However, higher values near 3 m/s are located along the northwestern coast of Cuba 

and also within the Bay of Campeche. Generally, the v-component RMSD field has larger 

RMSD values than the u-component, ranging from 1.5–2.5 m/s. In the v-component and the 

wind speed RMSD fields there is a peculiar peak located over the oil spill region, which may 

suggest that the spill itself is influencing the magnitude of the winds in at least one of the 

products being compared. The larger u and v component RMS differences compared to the speed 

RMS differences strongly indicate biases in direction.
42

Figure 3-15. Wind speed bias (a) and wind vector difference (b) between CCMP and COAMPS

Positive (negative) values in the wind speed bias plot (Figure 3-15a) indicate that the average 

CCMP winds are larger (smaller) than the COAMPS winds. Note a positive peak (~1.5 m/s) over 

the oil slick location showing that the CCMP winds are stronger than the COAMPS winds in this 

area. The black box designates the region where the wind vector difference is shown (Figure 3-

15b). The wind vector difference plot (Figure 3-15b) highlights the fact that there are significant 

differences in the v-component of the wind between the two datasets. This is the main 

contributor to the wind speed bias. (Background contours show the absolute values of the wind 

vector difference.)
43

Figure 3-16. Spreading of surface oil during passage of atmospheric front, 24-28 April 2010

(a) Time-integrated surface oil volume (m3/km2) derived from TCNAA SAR data on 04/24 and 04/28/2010. 

Blue lines are SSH contours from HYCOM (experiment 31.0). Solid lines are positive SSH drawn every 

0.1 m. Dashed lines are negative SSH every 0.05 m. (b) Surface weather map for April 25, 2010, 12Z. (c) 

CCMP and (d) COAMPS wind vectors 4/25/2010 12Z

3.2.2b. COAMPS vs CCMP: comparison of frontal structure 

Usually the GoM is characterized by a low wind speed regime during the summer months. 

However, three frontal passages occur in the area within the first month after the oil spill. These 

are large wind forcing events that also have the widest spread between CCMP and COAMPS. A 

frontal passage through the GoM occurred from April 25–26, 2010. The NOAA surface weather 

map (Figure 3-16b) shows the location of the front on April 25th at 12Z. Strong southerly flow 

ahead of the front forced the oil slick northward into coastal regions along Louisiana, 

Mississippi, and Alabama (Figure 3-16a). This example demonstrates the importance of the 

timing and structure of fronts in the wind data forcing oil transport models. As can be seen in the 

CCMP (Figure 3-16c) and COAMPS (Figure 3-16d) wind vector plots from 25 April 2010 12Z, 

the frontal structure is similar, but not an exact match. In COAMPS, the front extends further 

south than in CCMP, but overall the direction and magnitude of the wind vectors is similar. An 

important distinction between CCMP and COAMPS is that the low wind speed areas behind the 

front are stronger in magnitude in CCMP than in COAMPS. This likely contributes to the results 
44

of the RMSD and wind speed bias plots discussed earlier. It is worth mentioning that under the 

strong wind conditions oil spread across SSH contours shown in Figure 3-16a. 

Figure 3-17. Spreading of surface oil during passage of atmospheric front, 30 April–5 May 2010 

(a) Time-integrated surface oil volume (m3/km2) derived from TCNAA SAR data on 04/30 and 05/05/2010. 

Blue lines are SSH contours from HYCOM (experiment 31.0). Solid lines are positive SSH drawn every 

0.1 m. Dashed lines are negative SSH every 0.05 m. (b) Surface weather map for April 25, 2010, 12Z. (c) 

CCMP and (d) COAMPS wind vectors 05/03/2010 12Z.

A second frontal passage occurs about a week later on May 3, 2010 12Z. The NOAA surface 

weather map shows the front located directly over the oil spill region (Figure 3-17b). The oil 

slick map for April 30–May 4 is growing in area and beginning to impact coastal Louisiana 

(Figure 3-17a). The frontal structure is similar between CCMP (Figure 3-17c) and COAMPS 

(Figure 3-17d); however, COAMPS has much larger wind speeds ahead of the front. In contrast, 

CCMP has stronger winds behind the front. Interestingly, CCMP also has a wind speed 

maximum located over the oil slick area that is not seen in the COAMPS wind vector plot. 
45

Figure 3-18. Spreading of surface oil during passage of atmospheric front, 6 June–10 June 2010

(a) Time-integrated surface oil volume (m3/km2) derived from TCNAA SAR data on 05/06 and 05/10/2010. 

Blue lines are SSH contours from HYCOM (experiment 31.0). Solid lines are positive SSH drawn every 

0.1 m. Dashed lines are negative SSH every 0.05 m. (b) Surface weather map for May 09, 2010, 12Z. (c) 

CCMP and (d) COAMPS wind vectors 05/09/2010 12Z.

The final front passed through the oil spill region on 8–9 May 2010 (Figure 3-18). The 

orientation of the trough axis was directed east-west (Figure 3-18b), unlike the other two 

examples. This resulted in strong northerly winds (Figure 3-18c) that would help push the oil 

slick away from the coastline. Once again, the CCMP and COAMPS wind fields are very similar 

in nature. They share the same wind maxima south of the Florida peninsula and other mesoscale 

features within the gulf are comparable.
46

3.2.2c. CCMP and COAMPS winds and buoy data over the oil spill

The previous analysis of CCMP and COAMPS data has revealed that on average CCMP 

winds are stronger than COAMPS winds over the oil slick location. Also the wind vector 

difference highlights substantial changes in the v-component of the wind. It is conjectured that 

the presence of oil biased wind speed retrievals in CCMP and COAMPS. An oil and sea water 

mixture changes the physical properties of the ocean surface and therefore affects radiometer 

wind speed retrievals. Presence of oil slick on the sea surface damps capillary waves altering 

back-scattering of the impulse leading to biases in retrievals of the wind vectors. 

To further investigate this, CCMP and COAMPS winds are compared to NDBC buoys

located within the oil slick (Figure 3-19). Four buoys were located within the oil spill: 42040, 

42364, 42374, and 42376. One buoy, 42039, was located outside of the oil slick. The comparison 

reveals that CCMP winds are consistently large than COAMPS winds over the oil slick 

especially during frontal passages (Figure 3-20, top). CCMP winds are in better agreement with 

NDBC buoy measurements than COAMPS winds. Analysis of the Special Sensor 

Microwave/Imager (SSM/I) radiometer wind speed shows no significant bias relative to the 

NDBC buoy wind speed in oil free regions (Figure 3-21, top), whereas a small negative bias is 

evident in SSM/I radiometer wind retrievals in oil slick regions (Figure 3-21, bottom). 

Figure 3-19. Locations of the NDBC buoys and closest CCMP and COAMPS grid points
47

Figure 3-20. CCMP, COAMPS, and NDBC buoy wind speeds. 

Top: buoys within the DwH oil slick. Bottom: buoy 42309 outside the oil slick. Note substantial negative 

bias of the COAMPS winds compared to the buoy winds within the oil slick during strong wind events in 

April and May. CCMP wind speeds are in better agreement with NDBC buoy wind speeds than COAMPS 

within the oil slick.

Figure 3-21. NDBC buoy 42040 surface winds and SSM/I wind speed retrieval before (top) and after 

(bottom) the DwH oil spill
48

This analysis leads to the following conclusions regarding wind forcing based on different 

wind data sets:

 The CCMP and COAMPS wind data sets exhibit considerable differences in 

the RMSD plots. Over a majority of the GoM, the RMSD values range from 

1-3 m/s for the u- and v-components of the wind and also the wind speed. A 

peak in RMSD values collocated with the oil slick in the northern GoM

suggests that the oil plays a role either in modifying the winds or in effecting 

wind retrieval from scatterometers and radiometers. On average, the CCMP 

wind speeds are larger than COAMPS in the oil slick area. It is also 

interesting that the v-component RMSD values are larger than those for the ucomponent. 

 Three distinct frontal passages were analyzed during the early stages of the 

DwH oil spill. After mid-May a more typical summertime weather pattern 

setup, devoid of frontal passages. Overall, the timing and structure of the 

fronts depicted by CCMP and COAMPS were very similar and matched well 

with the NOAA surface weather maps. However, some differences in wind 

direction and magnitude were present during these events. Commonly, the 

CCMP data exhibited stronger wind speed magnitudes than COAMPS in areas 

with light winds. Because the GoM during summertime typically has light 

winds, this fact could be the reason for the high RMSD values.

 The oil slick may introduce a wind speed bias into gridded products such as 

CCMP and COAMPS through compromised radiometer wind speed retrievals 

and reduced surface roughness

 CCMP wind speeds are in better agreement with NDBC buoy wind speeds 

than COAMPS within the oil slick. 

3.2.3. Validation methodologies for surface oil drift models and quantitative 

comparison techniques for spatial patterns

The DwH oil spill clearly demonstrated the necessity for techniques of oil model comparison 

and validation. One of the major hurdles in oil inter-model comparison is the lack of an

objective, generally-accepted evaluation system or technique for model skill assessment and 

comparison, especially for spatial patterns. Such techniques are particularly useful for model 

adjustment during iterative process when convergence of the numerical solution to the control 

field (e.g., observation) has to be estimated. Currently, oil model evaluation is predominantly 

qualitative based on visual comparison of the simulated and observed oil spills. During the DwH 

forecasting, several validation metrics were employed to assess reliability of the oil forecasts 

performed by the COAPS group, such as oil particles matching between the model and 

observations, locations of the centers of mass, and dispersion of oil particles around the center of 

mass. These metrics provide a first insight on the accuracy of the model forecasts picking up 

obvious mismatch between the simulated and the control field (observations) mostly related to 

apparent divergence in oil spill locations. However, these methods were not able to distinguish 

more subtle, yet important, discrepancies in shape orientation and geometry. Below, several 

possible methods for objective quantitative evaluation and skill assessment of surface oil models 
49

based on spatial distribution of scalar fields (e.g., oil concentration or volume) are investigated. 

A complete analysis of considered methodologies can be found in Dukhovskoy et al. (2015a), 

where these techniques were demonstrated for the sea ice application and in Hiester, et al. (2016) 

for river plume application in the GoM. 

3.2.3a. Overview of Validation Methodologies

The following validation metrics will be considered for skill assessment of modeled surface oil 

fields: absolute deviation (AD), root mean square deviation (RMSD), mean displacement (MD), 

Hausdorff distance (HD), and modified Hausdorff distance (MHD). The metrics can be applied 

in both 2D, using contours of oil edge (e.g., some small concentration value), and 3D, using both 

the contours and scalar field within. 

Absolute Deviation

The absolute deviation (AD) metric used here is defined as 

(3-1)

where, in 2D, Ci and C0 are the areas inside the oil edge contour  derived from the tested and 

control data sets, respectively

(3-2)

AD measures an absolute difference in the oil area. For a 3D application (the metric is referenced 

as “AD3D”), Ci and C0 are defined as 

, (3-3)

, (3-4)

where gi and g0 are the tested and control scalar fields (e.g., oil thickness). If these functions 

describe oil concentrations, equations (3-3) and (3-4) are oil areas in the two data sets bounded 

by contour .

Root Mean Square Deviation

Root mean square deviation (RMSD) or root mean square error (RMSE) is a frequently used 

quantitative estimate of difference between model output and the control data. For 3D 

applications, when spatial distribution of oil characteristics is assessed, RMSD-based techniques 

are usually calculated through consideration of point-by-point differences over the domain of 

interest.

For 2D application, RMSD can be employed to analyze spatial patterns through comparison 

of contours such as oil spill edge or any other contours inside the oil spill. For two given sets of 

points A and B, the RMSD score is defined as





C  dA   



Ci gi x,ydA   



C0 g0 x,ydA
50

, (3-5)

where d(ai, bi) is an error (here, distance) between the i

th pair of points ai  A and bi  B. In the 

case of two contours, RMSD is estimated for n data points along the curves (Figure 3-22). Note 

that the definition of RMSD assumes a correspondence between the points in sets A and B (

). This information is generally not known in realistic applications. There are 

several possible techniques to determine the pair-wise correspondence of points between the two 

shapes (e.g., Agarwal & Phillips, 2006). Here, for two given sets of points, the two closest points 

from the sets are considered as the matching pairs (see Figure 3-22 for more detail). Another 

hurdle related to the RMSD application is when the number of points is different in the data sets 

A and B. Then extra steps are required before RMSD can be calculated in order to define a pointto-point correspondence between the contours (Figure 3-22b).

For 3D applications, when spatial distribution of oil characteristics inside a contour is assessed, 

RMSD-based techniques are usually calculated point-by-point over the domain of interest. In this 

analysis, the following methodology is employed for 3D applications 

, (3-6)

where h(i,j) and g(i,j) are tested and control scalar fields, respectively, of the grid cell (i,j). Note 

that this technique also requires a point-to-point correspondence between the scalar fields (as do 

most other RMSD-based methods). This means that the methodology does not account for 

information about position of the scalar fields. This can potentially debilitate sensitivity and 

robustness of this method to differences in oil patterns (see more detail in Dukhovskoy, et al., 

2015a). 

Figure 3-22. Schematic diagram of the point-to-point matching algorithm

(a) d(ai, bi) is the error (“distance”) used in calculation of the RMSD score between two contours. The 

pair-wise correspondence of points is determined as follows. Given two sets of points A and B containing 

the same number of points, a point a0 on set A is chosen at random as the initial point, and then the 

closest point b0 on set B is selected as its corresponding point (if more than one such point exists, one is 

selected at random). Point a1 is chosen to be the point in A\{a0}(point a0 is excluded from the set A) 



DRMSD A,B

d ai  ,bi

2

i1

n



n



ai AbiB

    

g g

n m

h i j g i j

D

m

j

n

i

R M S D D

m a x m i n

, ,

1 1

2

3 





 
51

closest to a0. Point b1 is chosen to be the point in B\{b0} closest to a1. The process iterates in this fashion 

until all pairs have been assigned. (b) Sets A and B have a different number of points. Suppose set A has 

m points and set B has n points, with m < n. Let k = n – m. The k longest edges on A’s contour are 

chosen and at each midpoint a point is added. Sets A and B now have the same number of points, and 

the algorithm from (a) is performed.

Mean Displacement

The mean displacement (MD) metric originates in spatial statistics as a measure of spread or 

dispersion of a data set around a reference location (P0) and has been successfully applied to 

shape recognition (Chang, et al., 1991; Zhang & Lu, 2004). This technique offers a natural 

metric for surface oil model evaluation. The data set is defined as every point bounded by a 

given contour (Figure 3-23). In the 2D cases, which consider contour comparison (oil spill edge), 

the data set is defined as points on the contour. P0 is taken as the centroid of the data set. The 

skill is assessed by considering the difference between dispersions of the data points measured 

relative to P0 in the data sets. MD can be generalized for a 3D application (referenced as 

“MD3D”) where each location is “weighted” by functions h and g that describe distributions of 

some property (e.g., oil thickness) within the data sets A and B, respectively. The difference in 

the weighted dispersions gives the MD3D score (DMD3D): 

(3-7)

(3-8)

(3-9)

where is the distance between a point ai and P0 and similarly for (Figure 3-23). 

In this case, both the field shapes and distribution of the analyzed property within the fields are 

compared. Note, when h = g = 1, (Equations 3-7 through 3-9) give MD. 

Figure 3-23. Weighted Mean Displacements measured for the data sets A (a) and B (b)

The metric considers the overall dispersion of data points around a reference location (P0). The arrow 

shows the distance between P0 and a data point.



D A

1

n

hid ai  ,P0

i1

n





D B

1

m

gi d bi  ,P0 

i1

m





d ai  ,P0



d bi  ,P0
52

Topological Methods: Hausdorff Distance and Modified Hausdorff Distance 

Topological shape matching descriptors use a metric distance between the objects as the measure 

of shape similarity (Zhang & Lu, 2004). The topological approach is a natural means of 

evaluation when spatial distributions of simulated properties are compared across datasets. The 

Hausdorff distance (HD) is a classical, widely used method for both object and shape matching 

(Daoudi, et al., 1999; Dubuisson & Jain, 1994a; Huttenlocher, et al., 1993; Rucklidge, 1997).

Here, HD is defined as

, (3-10)

where , , and d(a, b) is the distance between the point a

on A and the point b on B. The distance d(a, b) can be Euclidean distance or another appropriate 

distance depending on the application (e.g., great circle distance). Generally, a smaller HD 

indicates better resemblance between the shapes. However, HD is sensitive to outliers. This is an 

undesirable property for a model validation metric, because a key requirement is the ability to 

quantify the overall resemblance of compared data fields. 

A modified version of the classic Hausdorff distance that is resistant to outliers is also considered

(Dubuisson & Jain, 1994b). In the present study, this modified Hausdorff distance (MHD) 

between two finite point sets A and B in the plane is defined as

, (3-11)

where |A| and |B| are the cardinality of sets A and B (here, the number of points on the contours), 

d(a,B) and d(A,b) are defined by and . As with HD, the 

distance d(a, b) can be Euclidian distance or another appropriate distance (depending on the 

application). Both metrics can be easily extended to a 3D application (referred as “HD3D” and 

“MHD3D”, respectively). 

3.2.3b. Sensitivity and robustness testing of the skill metrics with a shape (2D 

application)

The considered metrics undergo sensitivity and robustness tests in order to evaluate the 

ability of each metric to quantify similarity in shape between two objects. When quantifying 

similarity in shape, special consideration should be given to the ability of a skill metric to 

penalize important (with respect to each particular application) dissimilarities and ignore 

unimportant ones. A validation metric is sought which can appropriately penalize differences in 

1) scale, 2) translation, 3) rotation, and 4) noise (Figure 3-24). During each test, scores from the 

validation metrics are computed between the original control contour and a modified shape 

obtained from the original shape being subjected to changes in the course of the sensitivity 

testing. The tests have been performed for 150 randomly generated shapes. Results of the tests 

for three randomly selected shapes are summarized below. A detailed description of the tests and 

discussion of their results are given in Dukhovskoy, et al. (2015a).

 

DH D A,Bm a x s u p

aA

da,B,s u p

bB

dA,b      

     



da,Bi n f

bB

da,b



dA,bi n f

aA

da,b

     

















   

aA bB

M H D d Ab

B

d a B

A

D A B D AB ,

1

, ,

1

( , ) , m ax



da,Bm i n

bB

da,b



dA,bm i n

aA

da,b
53

Figure 3-24. Sensitivity tests with three randomly generated shapes (a–c)

The black dot is the reference point P0 (centroid) for calculating MD. Illustration of sensitivity and 

robustness tests are shown for Shape 1: (d) The shape is linearly scaled up to double the original size 

(grey contour). (e) The contour is rotated clockwise. The grey contour is rotated 210clockwise relative to 

the control black contour. (f) The contour is linearly translated in all directions along a unit vector u. Only 4 

directions are shown. (g) Normally distributed noise is added to the coordinates of the contour. The grey 

contour is obtained by adding the normally distributed random noise with 

2

= 9x10-3.

Scale Test: For all shapes considered, each metric exhibits the desired behavior of being 

strictly increasing (Figure 3-25). However, the graph of RMSD has sharp discontinuities that 

exhibit an unstable response of the metric to scaling. This result indicates that although RMSD is 

generally responsive to the scaling, it is not responsive in a robust way. This may result in biases 

in skill assessment of model performance when using RMSD. It is noteworthy that the scale test 

is the simplest of the tests that the metrics are subjected to. A metric that is in any way a 

measurement of shape is expected to pass this test. All but RMSD have successfully passed the 

test.
54

Figure 3-25. Metric scores against scale for three shapes (a–c) shown in Figure 3-24

Different lines represent different skill metrics (see legend). AD score has been normalized by the 

maximum value for demonstration purposes.

Rotation Test: The above analysis shows that none of the tested metrics are completely 

responsive to rotation (Figure 3-26). However, all the metrics except AD can penalize

disorientation of two contours when one contour is rotated over a small angle relative to the 

other. This is particularly relevant to the practical application considered here because 

anticipated biases in spatial orientation of contours of geophysical fields (including oil 

concentration contours) have small range; therefore, it is most important that a skill metric is able 

to penalize rotation over angles much less than /2.

Figure 3-26. Metric scores against clockwise rotation angle (degrees) for three shapes (a–c)

The AD score is not shown, as it remains 0 for all rotation angles. All scores have been normalized by 

their maximum values.

Linear Translation Test: The resulting metric scores (Figure 3-27) are represented by a 

surface with a minimum at (0,0) where the two contours are collocated. With RMSD (Figure 3-

27a, e, i), the metric score surface has obvious discontinuities, demonstrating limited sensitivity 

of this metric to translation. The other metrics all illustrate a good response to the translation test. 
55

Figure 3-27. Metric scores for the linear translation test that the three shapes (in the rows) are 

subjected to

Each row corresponds to a different shape and different tests are shown in the columns. The horizontal 

and vertical axes are distances between the centroids of the shifted and control contours. The zero score 

corresponds to the initial location of the contour (0,0). (a, e, i) RMSD. Note discontinuities in the metric 

field. (b, f, j) MD. (c, g, k) HD. (d, h, l) MHD. The AD metric takes the zero value everywhere and is not 

shown. The black contours are drawn at intervals of 0.1.

Random Noise Test: The preferred behavior of the metric score function is an overall 

increasing trend with points representing individual test outcomes distributed around the trend. 

The smaller the uncertainty range, the more robust the metric to noise. Both the AD and MD 

approaches (Figures 3-28a, f, k, and 3-27c, h, m) completely fail the test. The uncertainty range 

increases rapidly with noise amplitude. The very large spread of AD and MD scores indicates no 

robustness of these methods to noise. For RMSD there is a positive trend in the scores with 

growing noise amplitude (Figures 3-28b, g, l). However, there is a clear discontinuity and a very 

large spread of the metric scores for noise with an amplitude greater than ~3×10-

-2 RMSD is, 

therefore, partially responsive to noise and demonstrates good performance for noise with small 

amplitude. For practical applications, this means that RMSD scores for the oil contours with 

larger perturbations can be inflated and hence are not robust. Both HD and MHD methods 

demonstrate a very robust response to added noise (Figure 3-28, two bottom rows). A slightly 

larger spread of HD scores for the noise with larger amplitude is due to the sensitivity of HD to 

outliers. 
56

Figure 3-28. Metric scores (vertical axis) from the noise test for the three shapes (in columns)

(a, f, k) AD, (b, g l) RMSD, (c, h, m) MD, (d, i, n) HD and (e, j, o) MHD. The horizontal axis is the noise 

amplitude. Note the differences in the uncertainty range (spread) of the scores for different metrics.

3.2.3c. Sensitivity and robustness testing of the skill metrics with a shape and 

distribution function (3D application)

The sensitivity tests described in section 3.2.3b compared 2D shapes (contours) of two 

objects. In surface oil modeling applications, differences of the distributions of oil spill

characteristics inside the shape (oil thickness, age, chemical composition, etc.) are also important 

and may be greater than the dissimilarity in the shapes of the objects. Subjected to the same tests, 

the metrics are used in the context of pattern recognition in order to assess their ability to 

quantify the resemblance of spatial patterns inside the shapes. In this test, field concentration is 

considered as the third dimension (Figure 3-29) and 3D Euclidean distance is used as a the 

distance d(a, b) in RMSD3D, MD3D, HD3D and “weights” in MD3D (g and h in Equations 3-7 

and 3-8) and AD3D. It should be noted that for the 3D application RMSD is calculated using 

Equations 3-6. 
57

Without normalization, Equation 3-6 is a 3D generalization of Equation 3-5 under the 

assumption that the grid points on two fields are collocated leading to zero horizontal distances. 

A detailed analysis of the test results for the 3D application is given in Dukhovskoy et al. (2015). 

Figure 3-29 summarizes the major findings of this effort.

Figure 3-29. Random shapes subsampled from a continuous field

(a) A randomly generated continuous field. Concentration fields inside random shapes (b–d) are obtained 

by “cutting” the random field with the shapes. For illustration, shape 1 is shown with the black contour. (b–

d) Random shapes (same as in Figure 3-24 a–c) with concentration fields inside. (e–h) Illustration of 

sensitivity and robustness tests are shown for Shape 1 similar to Figure 3-24 d–g. For reference, the 

empty contour indicates the initial position of the shape.

Scale Test: For all shapes, all metric scores are strictly increasing as scale increases (Figure

3-30). However, the graph of RMSD3D has a highly non-linear change in the response to 

linearly increasing scale. This suggests an unstable behavior of the metric to scaling with over 

sensitivity to small changes in scale but weak response to scaling for the remainder of the range 

of tested values. The other metrics maintain a nearly constant rate of the score change suggesting 

equal sensitivity to the tested scales.
58

Figure 3-30. Metric scores against scale for three concentration fields (a–c) shown in Figure 3-29

a–c

Different lines represent different skill metrics (see legend). AD3D score has been normalized by the 

maximum value for demonstration purposes.

Rotation Test: Adding a third dimension has improved performance of the metrics but only 

MHD3D is close to being completely responsive to rotation (Figure 3-31). 

Figure 3-31. Metric scores against clockwise rotation angle (degrees) for three concentration 

fields (a–c)

The AD3D score is not shown, as it remains 0 for all rotation angles. All scores have been normalized by 

their maximum values.
59

Liner Translation Test: All metrics except RMSD3D illustrate a good response to the linear 

translation test (Figure 3-32). RMSD3D (Figs. 3-32a, e, i) shows an initial large response, with a 

large score for small translation. 

Figure 3-32. Metric scores for the linear translation test (Figure 3-29g) that three shapes (in the 

rows) are subjected to

Each row corresponds to a different shape and different tests are shown in the columns. The horizontal 

and vertical axes are distances between the centroids of the shifted and control contours. The zero score 

corresponds to the initial location of the contour (0,0). (a, e, i) RMSD3D. (b, f, j) MD3D. (c, g, k) HD3D. (d, 

h, l) MHD3D. The AD3D metric takes the zero value everywhere and is not shown. The white contours 

are drawn at intervals of 0.1.
60

Random Noise Test: Similar to the 2D case, both the AD3D and MD3D metrics (Figures 3-

33a, f, k and 3.2.31c, h, and m) fail the 3D noise test. The RMSD3D method demonstrates a very 

robust response to added noise (Figures 3-33b, g, and i). HD3D and MHD3D again show a 

robust response to noise with monotonically increasing scores and a bounded spread through the 

whole range of noise amplitudes in the test. 

Figure 3-33. Metric scores (vertical axis) from the noise test for the three concentration fields

inside the shapes (in columns)

(a, f, k) AD3D, (b, g, l) RMSD3D, (c, h, m) MD3D, (d, i, n) HD3D and (e, j, o) MHD3D. The horizontal axis 

is noise amplitude normalized by the maximum score. 
61

Conclusions

Sensitivity tests with 2D contours and spatially distributed random scalar fields inside the 

contours demonstrate (3D) a better performance of the topological metrics compared to the other 

metrics. The AD metric is the simplest and least reliable skill metric, passing only the scaling

test. It is followed by RMSD, which has a better performance in 3D cases. MD exhibits good 

responsiveness to the three continuous deformation tests but fails the noise test. Both topological 

metrics, HD and MHD, have passed all four tests and although they are not fully responsive to 

rotation, they are sufficiently responsive for realistic applications. The study has determined that 

MHD is a mathematically tractable yet efficient method for model skill assessment and 

evaluation that has a particular focus on spatial patterns and distribution. It can be effectively 

applied to objectively evaluate and compare both 2D and 3D oil spill characteristics in the 

model-to-model and model-to-observations applications. MHD was shown to perform better than 

the other considered metrics, such as RMSD and AD, and can be easily applied to surface oil

model evaluation and assessment. MHD can also be applied to validation of simulated 

oceanographic fields where both shape and distribution are of importance, for example for river 

plumes, and sea ice models (Dukhovskoy, et al., 2015a; Hiester, et al., 2016). Thus, MHD is 

suggested as a validation metric for oil drift model for this study. 

3.2.4. Objective Oil Model Evaluation and Adjustment

Realistic application of the MHD to the surface oil model evaluation and adjustment is 

demonstrated for the sensitivity runs with a surface oil drift model described in section 3.2.4a. 

3.2.4a. Surface oil drift model

A two-dimensional surface oil drift model was developed for several applications in this 

project. This model is based on a Lagrangian particle-tracking trajectory algorithm commonly 

used in oil spill models (e.g., Abascal, et al., 2009; Garcia-Martinez & Flores-Tovar, 1999), with 

several user-modifiable options facilitating parameter sensitivity and other studies. The model 

simulates the evolution of an oil spill on the surface of the ocean by computing the twodimensional trajectories of large numbers of simulated discrete “particles”, or Lagrangian 

elements (LEs), each representing a volume of oil or related constituents. The initial positions of 

the LEs are prescribed such that the model may either be initialized from an observed map of 

surface oil distribution, or added to the surface of the ocean randomly distributed near the source

within the seeding radius (Rs). The trajectories of the LEs are computed within a time-varying 

velocity field that is some combination of ocean surface currents and 10 m winds from analysis, 

reanalysis, or forecast model products. LEs are removed from the system randomly based on a 

prescribed half-life (Th), or total life (Tmax), to simulate removal of oil from the surface through 

processes such as evaporation, dissipation into the water column, biological activity, or other oil 

removal mechanisms.
62

The oil particle trajectories are computed as a superposition of advective processes and 

turbulent diffusion

 t  t

dt

d

a d

u x, u x,

x

 

(3-12)

where ua is the advective (oil drift) velocity and ud is the diffusive velocity. The advective

velocity is calculated as a linear combination of the surface ocean currents, 10 m wind vector and 

waves

uaucCwu10 Csus

(3-13)

where uc is the ocean surface current velocity of seawater estimated from the topmost grid layer 

of an ocean model, Cw is the wind drag coefficient, u10 is the 10 m wind velocity over the sea 

surface and  is a unit vector being directed at angle from the wind (wind deflection angle), 

Cs is the wave coefficient and us is the wave-induced Stokes drift. Typically, the wind drag

coefficient varies from 0.025 to 0.044 (American Society of Civil Engineers Committee on 

Modeling Oil Spills, 1996). The role of waves in oil transport can be important in the nearshore 

areas and in the presence of swell waves. The inclusion of wave-induced transport provides a 

mechanism for beaching of surface oil (Sobey & Barker, 1997). Values of Cs have been 

estimated within the range around 0.05–1.5% of the wave-induced Stokes drift (e.g., Castanedo, 

et al., 2006). The Stokes drift can be estimated from the winds (under the assumption of fully 

developed wave field) or from the wave spectral field (more detail are in Section 4.3). In many 

applications, wind and wave effects are combined and both represented by the wind drag 

coefficient. Under light wind conditions without breaking waves, Reed, et al. (1994) found that 

Cw = 0.035 provides realistic simulation of oil slick drift in offshore areas. The wind deflection 

angle is typically 20 degrees clockwise from the wind direction (in the Northern Hemisphere). 

The model has an option to use wind-dependent wind deflection angle with negative values 

corresponding to the clockwise rotation from the wind direction (Samuels, et al., 1982)



  25oexp 108

 u10  k

g

(3-14)

where k is the kinemtatic viscosity of the seawater (1.0510-6 m

2

/s) and g is the gravitational 

acceleration. 

The last term in Equation 3-14 approximates oil spreading under the influence of turbulent 

diffusion and other unresolved lateral mixing processes. The turbulent diffusive velocities ud are 

approximated as random fluctuations defined for each time step based on “random walk”

(Garcia-Martinez & Flores-Tovar, 1999; Lonin, 1999; Abascal, et al., 2009); 



ud ud expi2

(3-15)

where  is a random variable from the standard Gaussian distribution and 



ud

6

t

(3-16)

where  is the diffusion coefficient (Table 3-5) and t is the time scale of the process. This 

diffusion can be spatially variable, for example it may be enhanced within a specified region 

about the source to simulate additional spreading by gravitational effects. 
63

Lagrangian oil particles represent a volume of oil that can be estimated as follows. Number 

of oil particles (Np) released at every time step is specified in the model. In the presented 

experiments, the number of oil particles released at every time step (Np) varies from 10 to 50 

particles per 1 hour of the model time. Then the volume of oil represented by each particle is

expressed as , where Qoil is surfacing rate of oil. For the time period between April 

23 and May 23, when the highest surface oil volume was observed (Figure 2-12), estimated 

surfacing rate of oil is Qoil ≈ 67 m3/hr. For reference, the total oil spill flow on May 29 was ~83–

125 m3/hr.

3.2.4b. Oil model evaluation and adjustment 

The MHD metric described in section 3.2.3 is employed as an objective measure for the 

proximity of the numerical solution to observations or other control field. The low MHD score 

indicates good resemblance of the simulated oil spill to the control field. This idea can be 

developed further by employing MHD to adjust model parameters by minimizing MHD score. 

For example, the shapes of the oil spill approximated by some contour of the oil concentration or 

thickness from the simulation is compared to the observed shape based on MHD. The model 

parameters were varied until MHD gives the lowest score (meaning the closest match). In order 

to test the performance of the suggested skill metric for this purpose, the following experiment 

was conducted. Using the surface oil drift model (section 3.2.4a) with a set of parameters a 

control simulation was performed (“Control” in Table 3-5). Five sensitivity experiments were 

designed where one of the parameters was not fixed. Each experiment had 10 simulations with 

one parameter varying within some range that includes the value of the control experiment. Then 

simulated oil distribution from each individual run was compared to the control experiment using 

MHD. The purpose was to see if the MHD technique would be able to correctly rank the 

sensitivity runs giving the best score to the simulations with parameters close to the control run. 

Table 3-5.

Oil model parameters

Parameter Sensitivity Experiments

Control Half-Life Wind 

Coefficient

Wind 

deflection 

angle

Eddy 

diffusivity

Seeding 

Radius

Np

(1) 80 80 80 80 80 80

Th

(2) 120 50.4:24:266.4 120 120 120 120

Tmax

(3) 360 360 360 360 360 360

Rs

(4) 4 4 4 4 4 1.5:1:10.5

dt(5) 24 24 24 24 24 24



(6) 5 5 5 5 0.9:0.9:9 5

Cw

(7) 0.035 0.035 0:0.006:0.054 0.035 0.035 0.035



(8)

-20 -20 -20 0:-4.5:-40.5 -20 -20

1 Number of particles seeded per time step

2 Half-life (hrs)

3 Maximum life (hrs)

4

Seeding radius (the radius within which the particles are seeded 

around the DwH location, km)

5 Model time step (hrs)

6 Eddy diffusivity (diffusion coefficient for a 2D random walk 

approximating eddy diffusivity, m2/sec2)

7 Wind coefficient

8 Wind deflection angle (degrees)



Qp Qoil Np
57

The MHD-based ranking of the simulations correctly identifies the case with the parameter 

value close to the control run (Figure 3-34). For example, in the half-life experiment (Figure 3-

34a) the simulation #4 identified as the best match to the control run has the half-life parameter 

equal to 122.4 hours that is the closest value to 120 hours in the control experiment. The 

sensitivity experiment demonstrates that MHD score function does not have a single minimum 

and for the wind deflection angle and eddy diffusivity parameters, the minimum is not well 

defined. This means that minimization of a MHD-based cost function for multiple parameters in 

the surface oil drift model may need a special approach. Minimization for eddy diffusivity 

parameter does not look promising. 
58

Figure 3-34. Metric scores (vertical axis) from the sensitivity experiments listed in Table 3-5

The horizontal axis is individual runs within each experiment. (a) Half-life; (b) Wind coefficient; (c) Wind 

deflection angle; (d) Eddy diffusivity coefficient; (e) Oil seeding radius. The minimum MHD score indicates 

the closest match to the control run. The vertical red dashed line indicates the true closest simulation 

based on the parameters of the control run.
59

3.2.4c. Estimation of half-life and wind coefficient from SAR observations and oil 

drift model by minimizing MHD score

The MHD technique is employed to estimate half-life and wind coefficient of surface oil 

particles from SAR observations by running sensitivity experiments with the oil drift model 

(section 3.2.4a). The choice of the model parameters to be adjusted is dictated by a relatively 

well-defined minimum of the MHD score function in the sensitivity experiments discussed above 

(Figure 3-34). MHD scores are calculated for the contours of 0.09 m3/km2derived from the 

simulations and SAR-TCNNA field. For the experiment, the model is integrated from April 22 

through May 23. The choice of the time period is based on the fact that until May 23 the surface 

oil volume was predominantly increasing except for a time period between 8 May and 18 May 

(Figure 2-12). After 23 May, the surface oil volume was decreasing. Because the surface oil rate 

in the model is constant, the surface oil discharge rate is estimated for the time period when the 

oil volume was increasing (section 3.2.4a). Another consideration was the length of the model 

integration given the number of sensitivity experiments. The considered time interval was long 

enough for oil spill to spread around the DwH and develop a very distinct and prominent shape 

with a “tiger tail” (Figure 3-35) that could be used for model evaluation and adjustment.

Table 3-6.

Oil model parameters for the half-life time estimation

Parameter Estimation 

of Half-life

Estimation of Wind 

Coefficient

Np(1) 150 150

Th

(2) 48–288 168

Tmax

(3) 840 840

Rs

(4) 20 20

dt(5) 6 6



(6) 5 5

Cw

(7) 0.035 0–0.0875



(8) Speed dependent Speed dependent

HYCOM analysis 20.1 20.1

Surface Winds CCMP CCMP

1 Number of particles seeded per time step

2 Half-life (hrs)

3 Maximum life (hrs)

4

Seeding radius (the radius within which the particles 

are seeded around the DwH location, km)

5 Model time step (hrs)

6 Eddy diffusivity (diffusion coefficient for a 2D 

random walk approximating eddy diffusivity, 

m

2

/sec2)

7 Wind coefficient

8 Wind deflection angle (degrees)
59

Figure 3-35. Time-integrated SAR-TCNNA surface oil volume (m3/km2), May 23, 2010, 00Z

The magenta line contours 0.09 m3/km2. The inset shows the control contour (blue) derived from the 0.09 

contour that is used for adjusting the surface oil drift model parameters by minimizing the MHD score.

Figure 3-36. Contours of the time-integrated surface oil volume (0.09 m3/km2), May 23, 2010, 00Z

Each panel presents 0.09 contours from the sensitivity model experiment (blue) and the control contour 

from SAR-TCNNA data (magenta) replicated from Figure 3-34. Half-life used in each experiment is listed.
60

Figure 3-37. Skill metric scores for time-integrated surface oil volume from the surface oil drift 

model simulations with varying half-life parameter and SAT/TCNNA data

The scores are calculated for the contours shown in Figure 3-35. The numbers on top of the bars are halflives (in hours) used in the simulations.

The control contour is derived from the SAR-TCNNA time-integrated surface oil volume 

fields (section 2-6). Only the longest contour that includes the DwH location is considered 

(Figure 3-35). Eleven simulations with parameters listed in Table 3-6 (column “Estimation of 

Half-life”) is initialized on April 22, 2010 and integrated through 23 May 2010. Each simulation 

has varying half-life parameter ranging from 48 to 288 hours. Contours of the 0.09 m3/km2are 

derived from the simulated surface oil volume in same way as for the SAR-TCNNA (Figure 3-

36). The MHD metric is employed to rank the model experiments with varying half-lives relative 

to the SAR-TCNNA in order to determine half-life value providing the best match with the 

observations. According to the MHD ranking (Figure 3-37), model simulations with half-lives 

144 and 168 hours have the closest agreement with the SAR-TCNNA for the considered time 

interval. Visual inspection of Figure 3-36 suggests that too short half-life excessively constrains 

spreading of oil from the DwH location. In contrast, too long half-lives result in a wide spreading 

of oil in the model compared to the SAR-TCNNA data. The “tiger tail” is fully developed in the 

simulations with half-lives greater than 168 hours. However, these simulations forecast excessive 

oil spreading in the south direction in contradiction to the SAR-TCNAA observations, which 

show small amount of oil in that direction (Figure 3-35). The persistence in forming the south 

“tentacle” in the surface oil drift model simulations is explained by the presence of a cyclone to 

the southwest of the DwH in the HYCOM analysis fields (cyclone “1” in HYCOM in Figure 3-5) 

during this time period. The southward currents to the east of the cyclone center transported oil 

from the DwH area. 

Next, the wind coefficient is evaluated in similar fashion. In this set of numerical 

simulations, the half-life is kept constant at 168 hours, based on the results from the previous 

sensitivity experiment. The wind coefficient varies from 0 (meaning no wind forcing) to 0.0875 

(Table 3-6, column “Wind Coefficient”). Results from the model experiments demonstrate high 

sensitivity of the surface oil transport to the wind forcing (Figure 3-38). The simulation with no 

winds (Cw = 0) provides a reasonably good forecast of the oil distribution because ocean surface 
61

currents that transport surface oil are partly driven by the wind stress. Thus, wind effect is 

already in the ocean currents. Nevertheless, the model experiments indicate that additional wind 

effect needs to be taken into account to correct oil spreading driven by ocean only. For example, 

the simulation with zero wind effect simulates general shape similar to the SAR/TCNNA contour 

with the “tiger tail” and orientation of the shape. However, oil has not spread as far north 

towards the coast as it is in the SAR-TCNNA. The MHD ranking shows that the solution 

improves as the wind effect increases up to Cw=0.0175 and degrades afterwards (Figure 3-39). 

More discussion of the results of this sensitivity experiment is provided in section 4-1. Too high 

values of Cw result in excessive spreading of oil in the southwest-northeast directions, most 

probably related to the predominant wind directions during front passages (discussed in section 

3.2.2b). It is worth noting that low Cw prevent development of the south “tentacle” in the 

simulations in agreement with the observations. It is conjectured that in the simulations with 

Cw>8.75e-3, southward winds pull oil to the south where it gets entrained into the southward 

flow driven by the cyclone producing too much oil in this area. It is reasonable to assume that 

location of the cyclone is slightly shifted northward in the HYCOM analysis 20.1. 

Figure 3-38. Contours of the time-integrated surface oil volume = (0.09 m3/km2,), May 23, 2010, 00Z

Each panel presents 0.09 contours from the sensitivity model experiment (blue) and the control contour 

from SAR-TCNNA data (magenta) replicated from Figure 3-34. Wind coefficient used in each experiment 

is listed.
62

Figure 3-39. Skill metric scores for time-integrated surface oil volume from the surface oil drift 

model simulations with varying wind coefficient parameter and SAT-TCNNA data

The scores are calculated for the contours shown in Figure 3-36. The numbers on top of the bars are 

wind coefficients used in the simulations.

3.2.4d. Hindcast modeling of natural oil slicks 

Parameters of the surface oil drift model were evaluated and adjusted in the course of the 

hindcast modeling of natural oil slicks. The trajectories of oil slicks observed in 41 SAR images 

were simulated using the two-dimensional model (section 3.2.4a). The hindcast experiments 

were forced by the ocean surface current fields obtained from the HYCOM experiment 31.0 and 

CCMP surface winds. The hindcast simulations showed a good agreement with satellite 

observations. Wind drift factor and wind deflection angle were adjusted to obtain the best 

agreement between simulated and observed oil slick shape and location. Simulations with the 

wind drift factor of 3.5% of the wind speed and a wind deflection angle of 20° to the right of the 

wind direction provided the best agreement between the model and SAR observations. 
63

4. EFFECTS OF OCEANOGRAPHIC AND WIND FORCING ON 

SURFACE OIL

4.1. WIND FORCING

Most surface oil drift models combine the effects of wind and ocean surface currents to 

estimate drift of oil particles in the ocean (section 3.2.4a). Primitive, although still often used, 

models use a constant wind factor of 0.03 or 0.035 of the wind speed and a constant deflection 

angle of 20 degrees to the right (for the Northern HS). It is uncertain, however, how this wind 

effect should be included when the ocean surface currents from modern high-resolution 

hydrodynamic models are used in the oil drift forecasts. The values of the wind coefficient were 

empirically determined from tracking drifters and surface floats in the ocean (see Reed et al., 

1994). The drifters are transported by the ocean currents integrated over the depth of the buoys’ 

draft. These currents approximate ocean “surface currents” that transport surface oil (are not true 

surface current that is), the velocity of the water parcels directly on the ocean’s surface. 

Even more substantial difference between the true surface current and its approximation can 

occur in numerical models. While oil surface drift occurs at the very surface of the ocean, what 

are termed “surface currents” simulated a hydrodynamic model more precisely represent the 

velocity integrated (or averaged) throughout the topmost vertical layer, which may be several 

meters thick depending on the vertical discrete resolution (see example in Figure 3-1). Thus, the 

so called “surface current” from a model may deviate substantially from the true surface current. 

This idea is demonstrated in Figure 4-1 showing a histogram of the magnitudes of vector 

difference between velocity vector of the topmost grid layer and upper 10 meter averaged 

velocity from HYCOM and NCOM analysis data sets. According to these results, a typical 

magnitude of the vector difference is 5–10 cm/s and the direction of the near-surface and 10 m 

averaged vectors (not shown) can be significantly different. Supposedly, in the case of 10 m 

averaged ocean velocities, a stronger correction to the surface current is needed in order to 

simulate the surface oil drift. Such a correction would be done by stronger wind coefficient. This 

difference highlights potential errors that may arise from using wind drift parameterizations 

developed from coarser vertical resolution ocean models (where surface currents may actually be 

upper 10 m averaged currents). However, if a model has a very well-resolved surface layer with 

correct physics in the turbulent boundary layer, then the model’s surface currents would be more 

representative of reality. It is likely, though, that even the sophisticated regional and basin scale 

real-time models used today still do not resolve well the near-surface current shear.

Sensitivity model experiments with varying wind coefficient presented in section 3.2.4.c 

have demonstrated high sensitivity of the surface oil drift to wind factor. The best match of the 

simulated surface oil distribution to the SAR-based estimate was in the simulation with wind 

coefficient of 0.0175. The result demonstrates that the best solution was derived with 

significantly lower wind effect on the oil transport, compared to a typically used value of 0.035. 

Further investigation is needed to provide a physical explanation of why the smaller wind factor 

provided the best solution. This could be due to the fact that in the considered HYCOM 

reanalysis the vertical grid provides a sufficiently high resolution of the upper ocean (Figure 3-1) 

allowing better representation of momentum transfer from winds into the ocean and the upper 

ocean currents are more responsive to the surface wind stress in this experiment (20.1). Thus, a

smaller wind factor was needed to adjust the ocean velocity from the topmost model layer to the 

true surface current.
64

Figure 4-1. Magnitude of vector difference between velocity vector of topmost grid cell and upper 

10 m averaged velocity from HYCOM and NCOM

The differences are calculated for daily time series of velocity fields from HYCOM (red) and NCOM (blue) 

analyses near DwH for the time period 25 May–30 September 2010. The horizontal axis is the magnitude 

of vector difference in m/s. Differences in HYCOM and NCOM shear are due to vertical grid, wind forcing, 

and turbulent closure schemes.

4.2. OCEAN SURFACE CURRENTS

The relative importance of surface currents have been examined in comparison to wave- and 

wind-induced transport. Separate studies have been done for shelf regions (where currents are 

weaker on average) and for the deep ocean. Currents are found to be the most important 

contributor to surface transport under many conditions, but waves are often also important. 

Probability distributions of the impact of currents will be provided.

We had planned to adjust from average currents in the model’s upper layer to a true surface 

current. However, investigation into the problem has shown that researchers working on this 

topic are much more divided than originally thought. Separately funded activities have indicated 

several ways forward on this problem, but in only one case are the data needed to address this 

question readily available. Differences in models are illustrated in Figure 4-2.
65

Figure 4-2. Near-surface velocity profiles of the wind-driven current for three conceptual models

Classical Ekman spiral for constant viscosity (blue; Ekman 1905), logarithmic boundary layer or K-Profile 

Parameterization (KPP) scheme (green; Large, et al., 1994) and momentum injection by breaking waves 

(red; Kudryavtsev, et al., 2008). Corresponding arrows show surface velocity vectors. Black arrow 

indicates wind stress and purple arrow is the depth-integrated Ekman transport (same to all models). The 

data to overcome this huge disparity do not exist at this time. (Source: Nikolai Maximenko)

For this study, surface currents are compared to Stokes drift and wind drift. Surface currents 

were provided by HYCOM (Chassignet, et al., 2007) and wind drift is approximated at 2% of the 

CFSR 10 m wind speed. The wind drift is reduced because it is found that Stokes drift accounts

for a large fraction of the wind drift (see section 4.3). Both products are models, which does lend 

some uncertainty to their accuracy. However, both HYCOM and CFSR are data-assimilative, 

meaning that they should be reasonably robust for this purpose (see caveats in section 4.1). The 

comparison is described in section 4.3. 

4.3. SURFACE WAVES

Stokes drift is the lateral displacement of a particle in the direction of wave motion due to the 

orbital motion of a passing wave train in a body of water. This effect was described by Stokes 

(1847) for finite-amplitude gravity waves. The displacement is due to the orbital motions not 

forming closed loops, which itself is due to the diminishing effects of horizontal displacement 

with increasing depth. The forward motion of a particle on the surface at the crest of a wave (top 

of the orbit) is larger than the counter-motion backward at the trough of the wave (bottom of the 

orbit). The net displacement over a single wave period is the result of a second-order term in the 

overall wave motion equation, and so it is sometimes disregarded in models.

There have been many past studies relating to the theoretical importance of Stokes drift to the 

overall transport of surface oil in the ocean. For example, Sobey and Barker (1997) used an 

idealized model of a near-shore region with along-shore current and onshore waves to determine 
66

the relative importance of Stokes drift. They found that Stokes drift was responsible for onshore 

beaching of surface oil in their model, in part due to the natural refracting of waves towards the 

shore by shoaling. LeHenaff, et al. (2012) used DwH-observed oil locations to verify a model in 

which Stokes drift was included via a more complex derivation of the wind. They also found that 

including wind-induced drift (which includes Stokes drift) was beneficial to accurately modeling 

the movement of surface oil.

Stokes' wave theory is predicated on the assumption of a homogeneous, incompressible fluid 

with uniform depth through which the wave is passing, and with the wave itself being of constant

velocity and form throughout. For the purpose of determining Stokes drift velocities for a single 

ocean wave over a single wave period, these assumptions can be considered valid. 

Stokes drift (Us) averaged over a single wave period is given by

2sinh ( )

cosh(2 ( ))

2

2

k H

a k k H z

Us







(4-1) 

where a is wave amplitude, k is the wave number, σ is the wave frequency, z is the depth being 

considered, and H is the depth of the water column on which the wave is occurring. This form is 

an approximation that is valid when wind waves and swell are not mixing. In the GoM, waves 

are almost always wind-driven, with very little swell entering from the Atlantic Ocean. When 

swell is caused by local hurricanes and fronts, the above equation works well because the Stokes 

drift from the swell dominates the calculation. We do not recommend this approach outside of 

semi-enclosed seas. Variables commonly found in wave models, however, include wave height 

(h), wave period (T), and wavelength (L) (Monismith & Fong, 2004). Using the equations

h  2a , (4-2) 



2

T  , (4-3) 

and

k

L

2

 , (4-4) 

and assuming z = 0 because for this study only Stokes drift at the surface is being considered, 

Equation 4-1 becomes

 

2sinh ( )

cosh 2

2

2 2

k H

k H

LT

h

US





(4-5) 

Stokes drift is often included in models via an approximation based on the wind speed and 

direction. This approximation includes surface motion due to Ekman transport, and the result is 

called “wind drift” (Weber, et al., 2014). The common range of values for this combined 

transport used in models is 2–5% of the 10 m wind speed at an angle 20° to the right of the wind 

direction. This approximation arose out of necessity, at a time when ocean and wave models 

were of insufficient resolution to allow for direct accounting of either Stokes drift or Ekman 

transport. However, this parameterization requires several assumptions, including wind wave 
67

equilibrium (that is, a steady sea state) and sufficient time for a full Ekman balance to develop, 

neither of which are reasonable in real-world conditions. Modern ocean models no longer require 

this approximation to be made, because Ekman transport can be more directly modeled. 

Additionally, modern wave models allow for Stokes drift to be calculated directly, which makes 

using a wind speed approximation for wind drift unnecessary and unreasonable (particularly if it 

results in double the Ekman motion).

One reason for considering Stokes drift separately from using the wind speed-derived 

approximation is that Stokes drift can be present even without the presence of wind. Swell, by 

definition, is a wave that has propagated away from its area of formation. Stokes drift will be 

present for any wave, even swell, so therefore it will be present when swell is the only present 

portion of the wave spectrum. This is a situation that often occurs when there is calm, or with a 

very weak wind or recent wind such that the local wind wave field is not fully developed. Swell 

can also be present for higher local wind speeds if the swell waves are very large. A tropical 

cyclone is one example of a storm system that can create swell in the GoM. 

Similarly, when considering oil as the material being transported, it is important to note the 

effect its presence on the ocean surface has on local waves. Oil is well-known for its tendency to 

inhibit the development of small, capillary waves on a water surface due to its viscosity. It also 

can reduce air-sea friction, further inhibiting local wave development. However, this impact is 

minimized on swell originating away from an oil slick, which results in the presence of the 

Stokes drift factor in regions affected by an oil spill. Additionally, when high wind speeds occur 

over an oil slick, the slick tends to break apart, reducing this effect, both through turbulence in 

the upper ocean layer and wave breaking, mixing oil into the water column. That said, wind 

speeds of sufficient magnitude to do this are not commonplace in the GoM except in tropical 

cyclones and winter cold frontal passages.

Basing Stokes drift calculation on wind speed at a given level also ignores the different 

manners in which the wind can interact with the ocean surface, and thus how the sea state will 

develop given a particular (e.g. 10 m) wind speed. A wind profile in a stable atmospheric surface 

layer (which occurs frequently at night) will result in weaker wind stress at the atmosphere-water 

boundary, leading to smaller waves. Conversely, an unstable atmospheric surface layer, even 

with the same 10 m wind speed, will result in greater wind stress at the atmosphere-water 

boundary, leading to larger wave heights, and thus larger Stokes drift magnitudes for the same 

wind speed at the given height. Additionally, drift is affected by white capping and wave 

breaking, which violate the assumptions listed above. In the GoM, winds are typically light 

enough that Stokes drift is a useful approximation. 

4.3.2. Data and Methods

4.3.2.1 Wavewatch III

In order to undertake a quantitative analysis of wave transport in the GoM, a full, continuous 

gridded wave dataset is needed. This means that relying on observational data, such as wave 

reports from buoys and ships, is insufficient, as those sources are often temporally discontinuous 

and too sparsely located to provide a meaningful representation of the entire Gulf. Additionally, 

ship-based wave reports in particular are usually estimated rather than measured. Therefore, for 

this study, it is necessary to use data from a wave model to calculate Stokes drift.
68

4.3.2.1.1 Existing Data

An existing dataset was considered. The U.S. National Center for Environmental Prediction 

(NCEP) offers model hindcasts globally from 1997 to the present time, at three-hour intervals 

[NOAA 2009]. These wave hindcasts are forced by Global Forecast System (GFS) wind input. 

However, this data was incomplete in the GoM, as the gridded peak wave period output 

contained spatial gaps for undetermined reasons (Figure 4-3). In addition, the highest-resolution 

output available (1/15° grid) is only available for regions within approximately 100 km of coasts. 

Coarser data (1/6°) was available covering the entire GoM, but that also contained the peak 

period gaps. So in order to acquire a full, complete gridded dataset to cover the study period of 

April–July 2010, it became necessary to use a wave model to create high-resolution continuous 

wave data for the entire GoM specifically for use in this study. Accordingly, the Wavewatch III 

model was used.

Figure 4-3. Example of NCEP peak wave period data (s), 1 May 2010 00Z with the problematic gap 

circled

4.3.2.1.2 Running Wavewatch III

Wavewatch III (Tolman, 2009) is a spectral wind wave model that can simulate windgenerated local wave fields and swell propagating from non-local areas. It was developed by the 

Marine Modeling and Analysis branch of the Environmental Modeling Center within the NCEP. 

Wavewatch works by separating wave spectrum at each grid point into partitions by energy 

density peaks, as well as calculating peak and mean wave variables for the entire spectrum. The 

model has available a number of parameterization and other options. Wavewatch can also accept 

several input parameters as wave forcing and limiting mechanisms, including near-surface 
69

atmospheric winds, sea ice concentrations, as well as air and sea surface temperatures. For this 

work, ice is not included (because the area of interest is the GoM, where ice is not present). 

Atmospheric wind and temperature and SST are included, in order to allow for surface stress 

adjustment due to stability.

The model was set up with a 1/15° grid covering the GoM (18–31° N, 80–100° W) with a 

time step of 450 seconds, nested within a coarser grid (1/2° spacing) covering all of the north 

Atlantic Ocean (5° S-55° N, 5–100° W) with a time step of 900 seconds (Figure 4-4). Boundary 

conditions for the coarse grid (except for the western boundary, which is entirely land) and the 

initialization of both grids were done using the idealized Joint North Sea Wave Observation 

Project (JONSWAP) spectrum (Hasselmann, et al., 1973). This initialization was done to provide 

a starting point, after which the model was run for a two-week period (prior to the beginning of 

the study period) using model data (see below) for forcing, ensuring that any wave energy in the 

model not driven by real conditions would be dissipated before analysis data was generated, 

leaving only waves driven by actual wind. Similarly, by placing the coarse grid boundaries a 

considerable distance from the fine grid boundaries (which received boundary conditions from 

the coarse grid), JONSWAP-influenced wave energy was dissipated before it propagated into the 

fine grid.

Figure 4-4. Maps showing extent of the model domains used for Wavewatch III

(a) 0.5° Atlantic Ocean grid with box indicating the location of the GoM grid. (b) 1/15° GoM grid with the 

green region indicating depths  100 m and the orange region indicating depths > 100 m

Although the period of interest for this work is April–July 2010, Wavewatch was initialized 

at 15 March 2010 00Z and run through 10 August 2010 00Z. The early start was intended to 

provide a "spin-up" of the model to reduce non-real world-driven wave energy from the 

JONSWAP initialization. The extended run time allowed for sufficient additional data to be 

generated to calculate trajectories initialized as late as 31 July.

Wavewatch was forced using both atmospheric wind and air temperature and SST. Using 

both temperatures is important, as doing so provides the model with the ability to approximate 

the wind profile between the height of the "measurement" wind and the surface. Forcing data for 

this work was obtained from NCEP's Climate Forecast System Reanalysis (CFSR) (Saha, et al., 

2010). This product originally existed only for the period from 1979–2009, but was recently 

extended through 2010. Atmospheric wind was taken from the 10 m wind velocities. 

Temperatures were used from the water surface and from a height of 2 m. Water depths were 

provided by NOAA's World Geophysical Data Center two-minute Gridded Global Relief Data 
70

(ETOPO2v2) (National Geophysical Data Center/NESDIS/NOAA/U.S. Department of 

Commerce, 2001). It should be cautioned that this implementation of Wavewatch III does not 

include the effects of currents, which would slightly affect the wind stress levels at the air-sea 

boundary (by changing the wind velocity relative to a particular point on the surface, which 

would then be initially moving with the current rather than only with the wave motion).

The model was set to save output data at each hour. Three output variables were used for this 

work: significant wave height (h), peak wave period (T), and peak wave direction (θ). Notably, 

the peak period output did not contain the gaps present in the already-existing data (Figure 4-5). 

Figure 4-5. Example of WW3 modeled peak wave period data (s), 1 May 2010 00Z

The lack of a gap in this data is circled. 
71

Figure 4-6. Probability density function of Stokes drift magnitude distributions in the GoM, April–

July 2010

All of the input data were chosen as they are commonly available in both observations and 

models, and can be used to calculate Stokes drift. Wavelength (L) was then calculated from the 

peak period (T), using the deepwater wavelength

2

2

0

gT L 

(4-6) 

and the dispersion relation















L

H

L L

2

0

tanh

(4-7) 

to reach a final wavelength for use in Equation 4-5. Stokes drift was then calculated at each point 

in the GoM at each hourly time step from 1 April 2010 00Z to 10 August 2010 00Z.

4.3.2.2.1 Stokes Drift

In order to provide a basic overview of Stokes drift in the GoM, Stokes drift was calculated 

(using Equations 4-1 through 4-5) at each point in the model domain, restricted to the Gulf (i.e., 

excluding the portions of the Caribbean Sea and Atlantic Ocean present in the model domain) at 

each hourly time step from 01 April 2010 00Z to 31 July 2010 23Z, so as to include the time 

period during and some time before and after the DwH oil spill. Any points for which data 

needed for calculating Stokes drift (wave height, period, length, or direction) were missing were 

excluded, which occurred 0.08% of the time. Additionally, grid points at which the water depth 

was less than 1 m were also excluded, because these locations were so near to land that waves 
72

are either breaking or likely directed onshore. An example of the daily drift is shown in Fig. 

4.11.

It is also useful to consider Stokes drift as a comparison between deep water and shelf water 

areas. Because waves in sufficiently shallow water (how shallow is dependent on wavelength) 

interact with the ocean bottom, which results in changes to wave parameters such as height and 

speed that are included in the calculation of Stokes drift, it is prudent to consider whether those 

interactions make a significant difference in the magnitudes of Stokes drift as compared to 

deeper water, where waves do not interact with the bottom. Additionally, shelf water is 

characteristic of having weaker overall surface currents compared to deep water (this will be 

important later). For this study, the boundary between shallow-shelf and deep water is set at 100 

meters. This provides for the approximate separation of the shallow continental shelf areas from 

the deep ocean.

Figure 4-7. Directional rose plot of Stokes drift distributions in the GoM, April–July 2010

Directions are in oceanographic convention.

4.3.2.2.2 Trajectories

Another method of comparing Stokes drift between deep and shelf water is with the use of 

trajectories. This allows for considering not how Stokes drift magnitudes change at a point, but 

instead considering what would happen to a theoretical particle (such as a patch of oil floating on 

the surface) over time due to Stokes drift. Here, with Stokes drift velocities calculated at each 

grid point every hour, it is possible to consider the net displacement of a particle over a period of 

time (here chosen to be 24 hours) due solely to Stokes drift. Particles are considered massless 

and infinitely small, which means they offer no resistance to their theoretical movement. In order 

to determine net displacements, at each hour in the data, a tracer grid was initialized at each 

velocity grid point in the GoM. For each of 24 successive hours, a new position was calculated 

for each position grid point based on the Stokes drift velocity field for that hour, and except for 
73

the initial advection (when all grid points were co-located with the initial velocity grid points), 

the velocity was applied according to a Runge-Kutta interpolation of that hour’s Stokes drift 

velocity field. If a position grid point was at any time advected off the velocity field, it was 

considered stopped at its last known point for the remainder of the time. Once the tracer points 

had been advected for 24 hours, the distance they had ended up from their initial locations was 

then calculated.

4.3.3. Stokes Drift in the GoM

4.3.3.1 Stokes Drift in the Full GoM

For the GoM as a whole, the average Stokes drift magnitude was 3.99 km/day, while the 

median was 3.40 km/day (see Table 4-1) during the study period. This indicates that the 

distribution of Stokes drift magnitudes skewed towards smaller values (in fact, this is the case for 

all GoM-wide Stokes drift-related distributions considered in this study) (e.g. Figure 4-6). In

addition, there was a wide variation in the distribution of wave (and therefore Stokes drift) 

directions, although the vast majority of waves during the period did have at least some 

westward component (Figure 4-7).

Table 4-1.

Stokes drift magnitudes, entire GoM, April–July 2010

Parameter Value (km/day)

Mean 3.99

9th percentile 0.76

25th percentile 1.79

Median (50th percentile) 3.40

75th percentile 5.63

91st percentile 8.13

Standard deviation 2.88

4.3.3.2 Stokes Drift Comparison: Deep and Shelf Water

For waves occurring over shelf water, Stokes drift magnitudes averaged 3.70 km/day with a 

median of 3.09 km/day, while for deep water Stokes drift magnitudes averaged 4.12 km/day with 

a median of 3.52 km/day (Table 4-2). In addition, magnitudes in shelf water exhibited a smaller 

variation, with a standard deviation of 2.78 vs. 2.91 in deep water (Figures 4-8 and 4-9). It is 

unclear whether the larger average Stokes drift magnitudes in deep water are simply the result of 

larger waves due to higher wind speeds, or to what extent, if any, the slowing of waves due to 

bottom interaction was responsible.
74

Table 4-2.

Stokes drift magnitudes, deep and shelf water, GoM, April–July 2010

Parameter

Deep Water Value 

(km/day)

Shelf Water Value 

(km/day)

Mean 4.12 3.70

9th percentile 0.86 0.59

25th percentile 1.95 1.46

Median (50th percentile) 3.52 3.09

75th percentile 5.71 5.43

91st percentile 8.29 7.79

Standard Deviation 2.91 2.78

Figure 4-8. Probability density function of Stokes drift magnitude distributions for deep water 

(depth> 100 m) in the Gulf Mexico, April–July 2010. 
75

Figure 4-9. Probability density function of Stokes drift magnitude distributions for shelf water 

(depth  100 m) in the Gulf Mexico, April–July 2010. 

The small difference between the average Stokes drift magnitudes for deep and shelf water 

was not unexpected. Though conventional wisdom does hold that “shallow water” waves 

generally have larger Stokes drift magnitudes, this is in reference to shallow water waves which 

are interacting with the ocean bottom. However, that is not the comparison being made here. The 

large majority of waves occurring over shelf water (again, defined for this study as depths of 100 

m or less) are in fact still deep water waves by that definition. Waves in the GoM are rarely large 

enough to become “shallow water” waves except in very shallow water (for example, depths of 

less than 10 m), which in this study constitutes a very small number of grid points.

4.3.3.3 Comparison of Stokes Drift Trajectories

For the April–July 2010 period, the 24-hour trajectories exhibited similar differences to the 

distribution of Stokes drift magnitudes when comparing those initialized in shelf water vs. deep 

water, with a mean displacement of 3.46 km for shelf water compared to a 3.99 km mean 

displacement for deep water (Figures 4-10 and 4-11). While this is likely primarily caused by the 

larger Stokes drift magnitudes found in deep water areas, this difference is also influenced to a 

small degree by the fact that shelf water trajectories are more likely to reach a stopping point via 

encountering land than deep water trajectories are by leaving the domain (because deep water 

grid points mostly if not completely originate too far from land to allow Stokes drift to advect 

them to a land point within the 24-hour limit except under the most extreme of circumstances).
76

Figure 4-10. Probability density function of 24-hour Stokes drift trajectory displacements for deep 

water (depth > 100 m) in the GoM, April–July 2010

Figure 4-11. Probability density function of 24-hour Stokes drift trajectory displacements for shelf 

water (depth  100 m) in the Gulf Mexico, April–July 2010
77

4.3.3.4 Stokes Drift Compared to Current and Wind Drift

Though Stokes drift showed a small difference when comparing deep and shelf water areas, 

this difference is much more pronounced when comparing Stokes drift magnitudes to surface 

current magnitudes. In deep water areas, Stokes drift magnitudes were an average of 20.8% of 

the collocated surface current (Figure 4-12). However, in shelf water areas, Stokes drift was 

36.0% of the surface current magnitude on average, with a much wider distribution of ratios 

(Figure 4-13). This means that Stokes drift is a more significant relative factor in surface 

transport in shelf water areas, and so accurate representation of Stokes drift in transport models is 

more important in these areas.

Figure 4-12. Probability density function of Stokes drift to surface current ratios for deep water 

(depth > 100 m) in the GoM, April–July 2010.
78

Figure 4-13. Probability density function of Stokes drift to surface current ratios for shelf water 

(depth  100 m) in the GoM, April–July 2010

Stokes drift is a larger contribution to the total surface transport of oil in shelf water primarily 

due to these areas having smaller surface current magnitudes. As seen previously, there is not a 

large difference in Stokes drift magnitudes themselves between deep and shelf water, while large 

magnitudes of surface current are primarily found in the loop current and loop current eddies, 

which are largely confined to deep water areas.

Similar comparisons were made with Stokes drift and a percentage of the wind speed ("wind 

drift"). Here, 2% of the wind speed is considered, which is at the bottom of the range of wind 

drift parameterizations used in trajectory models. As shown in Figures 4-14 and 4-15, there is 

little difference in ratios when comparing between deep and shelf water areas (which is expected 

since both regions typically experience similar wind speeds). However, it can be seen that there 

is a wide distribution of ratios of Stokes drift to 2% of the wind speed in both figures. This 

indicates that there is poor correlation of Stokes drift and wind speed, implying that using wind 

speed as a proportional proxy for Stokes drift is not especially accurate. Because there is 

sometimes swell propagating into a region from elsewhere, and sea state does not 

instantaneously change in response to changing wind speeds, this is not an unexpected result.
79

Figure 4-14. Probability density function of Stokes drift to wind drift ratios for deep water (depth > 

100 m) in the Gulf Mexico, April–July 2010

Wind drift is considered to be 2% of the 10 m wind speed.

Figure 4-15. Probability density function of Stokes drift to wind drift ratios for shelf water (depth 

100 m) in the Gulf Mexico, April–July 2010

Wind drift is considered to be 2% of the wind speed.
80

4.3.4. Stokes Drift during a Hurricane

Though Stokes drift is induced by any wave, the largest magnitudes of Stokes drift are 

generally produced by the largest waves. Correspondingly, the largest waves in the GoM are 

produced by the strongest winds, which are almost always found in tropical cyclones. In 

addition, large waves produced by these storms propagate away to become significant swell in 

locations within the Gulf well away from their origins. This produces an extreme case in which 

Stokes drift as estimated from the local wind speed can be especially inadequate as a means of 

accounting for particle displacement.

During the four months of this study, two tropical cyclones passed through the GoM. 

Hurricane Alex occurred in late June, while Tropical Depression Bonnie (previously a tropical 

storm) occurred in late July. Bonnie was not considered for this study, due to being below 

tropical storm-force for its entire presence in the Gulf.

4.3.4.1 Hurricane Alex

Hurricane Alex formed in the Caribbean Sea on 24 June 2010 as a tropical depression, then 

strengthened into a tropical storm, crossing the Yucatan Peninsula and entering the southwestern 

GoM on 27 June with maximum sustained wind speeds of 35 kt. The storm then moved 

northwest across the western Gulf, strengthening into a category 2 hurricane with maximum 

sustained wind speeds of 95 kt, before making a second landfall on the northern Mexico coast on 

July 1 (Pasch, 2010). This resulted in a period of approximately 72 hours during which largeheight waves were being generated by increasingly strong winds across the southwestern GoM. 

These waves were of sufficient size and energy that they could propagate throughout the GoM as 

swell before dissipation.

To examine how Hurricane Alex affected Stokes drift magnitudes across the Gulf, two 

regions within the Gulf are compared before, during, and after the storm. The first of these 

regions was the southwestern Gulf, where the storm had a direct impact on wave heights, while 

the second region was the northeastern Gulf, where the oil spill was occurring and distant from 

the hurricane winds (Figure 4-16). For each region, Stokes drift magnitudes were compared 

during three seven-day periods: one before the storm entered the Gulf, one encompassing the 

entire time any part of the storm's circulation was over the Gulf, and one following the storm's 

landfall (see Table 4-3). Seven-day periods are intended to be short enough to prevent the effects 

of the hurricane on Stokes drift from being lost as noise, and are long enough to not be affected 

by daily and day-to-day random weather events.

Table 4-3.

Dates of weeks before, during and after Hurricane Alex for study consideration

Storm Location Dates Considered

Before entering GoM 6/19/2010–6/25/2010 

Impacting GoM 6/26/2010–7/2/2010 

After exiting GoM (Over land) 7/3/2010–7/9/2010 
81

Figure 4-16. Map of Hurricane Alex storm track

Lower left box is area considered to be directly affected by storm. Upper right box is area considered to 

be affected by oil spill. Adapted from the NOAA National Hurricane Center.
82

4.3.4.2 Stokes Drift During Hurricane Alex

During the week before Alex entered the GoM, Stokes drift magnitudes in the southwestern 

Gulf study region averaged 2.69 km/day, which was below the region's four-month average of 

3.03 km/day, with a standard deviation of 1.47 km/day (Figure 4-17 red line, and Table 4-4). 

Table 4-4.

Stokes drift magnitudes (km/day) of the weeks surrounding Hurricane Alex, in the area directly 

impacted by the storm

Parameter Before Alex During Alex After Alex

Mean 2.69 7.64 3.80

9th percentile 1.20 1.91 0.66

25th percentile 1.70 3.42 1.67

Median (50th percentile) 2.35 6.67 3.51

75th percentile 3.31 10.61 5.39

91st percentile 4.94 15.20 7.43

Standard deviation 1.48 5.16 2.58

Similarly, Stokes drift magnitudes in the northeastern Gulf averaged 1.89 km/day with a 

standard deviation of 1.82 km/day, which was also below that region's four-month average of 

3.22 km/day (Figure 4-17, red line, and Table 4-5). As Alex traversed the southwestern GoM, its 

winds generated large waves. This resulted in that week's average Stokes drift magnitude in the 

region rising sharply to 7.64 km/day, which lies at the 88th percentile of the four-month period 

(Figure 4-18, blue line, and Table 4-4). The distribution of magnitudes during that time period 

also increased drastically, with the standard deviation rising to 5.16. In the northeastern Gulf, the 

mean Stokes drift magnitude also rose, but only to 2.82 km/day, which was still below the area's 

four-month average. However, the variability of the magnitudes also increased to 2.41 (Figure 4-

18, blue line, and Table 4-5). It is likely that most (but not all) of the storm-produced wave 

energy during this period was still remaining in the southwestern GoM, with only a small amount 

having propagated away into the northeastern GoM to boost Stokes drift values there.

Table 4-5.

Stokes drift magnitudes (km/day) of the weeks surrounding Hurricane Alex in area impacted by oil 

spill

Parameter Before Alex During Alex After Alex

Mean 1.89 2.82 3.99

9th percentile 0.16 0.44 0.42

25th percentile 0.42 0.98 1.46

Median (50th percentile) 1.35 2.75 3.32

75th percentile 2.75 3.90 6.23

91st percentile 4.81 6.70 8.66

Standard deviation 1.82 2.41 3.01
83

Figure 4-17. Probability density functions of Stokes drift magnitudes in the area of the GoM

directly affected by Hurricane Alex

The red line with + marks is before the storm; the blue line with x marks is during the storm, and the green 

line with * marks is after the storm.

Figure 4-18. Probability density functions of Stokes drift magnitudes in the area of the GoM

affected by the oil spill

The red line with + marks is before the storm; the blue line with x marks is during the storm, and the green 

line with * marks is after the storm.
84

Once Alex made landfall in northeastern Mexico, its winds were no longer influencing wave 

development in the GoM. In the southwestern Gulf, this resulted in average Stokes drift 

magnitudes dropping to 3.80 km/day, or just under half of their during-storm values. The overall 

distribution of Stokes drift magnitudes in this region during this time was still more spread out 

compared to the week before the storm, with a standard deviation of 2.58, indicating a residual 

effect from the storm itself (Figure 4-17, green line, and Table 4-4). In the northeastern Gulf, 

Stokes drift magnitudes were even higher than they were during the storm, averaging 3.99 

km/day, with a similarly larger standard deviation of 3.01 (Figure 4-18, green line, and Table 4-

5). This is most likely due to propagating swell from the storm more readily influencing Stokes 

drift in this region.

Hurricane Alex is an extreme example of swell impacting Stokes drift magnitudes in an area 

distant from the storm. Large waves generated by the storm propagated northeastward as swell 

into the oil spill region, leading to larger Stokes drift magnitudes throughout the area even 

though there was no corresponding large increase in wind speeds. Estimating Stokes drift as a 

function of wind speed would have been especially inaccurate in this instance. Additionally, this 

case demonstrates the high variability of Stokes drift over a short time due to a single extreme 

weather event. 

4.3.5. Conclusions Regarding Currents, Stokes Drift, and Wind Drift

Stokes drift has been shown to be an important mechanism in the transport of surface oil, 

with an average magnitude of 3.99 km/day and an average 24-hour Lagrangian displacement of 

3.84 km for any particle during the study period of 1 April 2010 to 31 July 2010 in the GoM. 

Although this study was limited to a single body of water over a relatively short period of time, 

the basic principle, that Stokes drift is an important surface transport mechanism and should be 

accounted for in models in the most precise manner available, is logically applicable for any 

body of water on which wind waves form.

When comparing Stokes drift for waves occurring in water of relatively shallow depth (the 

continental shelf, delineated as water of depth less than 100 m) against that for waves occurring 

in deeper water (depth greater than or equal to 100 m), no physically significant difference in 

magnitude was found during the study period. This is most likely because waves are primarily 

wind-driven, and there was little if any difference in wind speeds over waters of differing depths. 

It is possible that for areas of sufficiently shallow depth where waves interact with the ocean 

bottom, there may have been a more notable difference, but the GoM has relatively small waves 

when compared with other basins, due to lighter winds and scarcity of swell, especially during 

the spring and summer months included in this study. This means that the number of grid points 

and times wherein waves would be impacted by bottom interaction was small.

However, when comparing Stokes drift magnitudes to surface current magnitudes, there is a 

notable difference between shelf and deep water. Stokes drift is a larger relative component in 

overall surface transport compared to surface current in shelf water, approximately double the 

percentage of the local current magnitude, when compared to deep water. This is primarily due 

to weaker surface currents in shallower water. Therefore, if calculating Stokes drift in an oil spill 

model with limited computing resources, it is more important to do so for shallow and shelf 

water areas, as Stokes drift is a larger part of total surface transport there. 
85

Though Stokes drift is often approximated in models as a percentage of the wind speed at a 

specific angle from the wind direction, it is not the best way to account for Stokes drift if the 

model includes a wave component. This is because Stokes drift has significant variation from the 

wind, due to the lag between changing wind speed and sea state response, as well as swell 

propagating into a region from elsewhere. This means that there are often waves (and with them, 

Stokes drift) occurring even when the local wind is calm. Thus, when designing a trajectory 

model which includes waves, it is preferable to calculate the Stokes drift component of transport 

directly from the wave parameters rather than using the wind speed approximation.

During weather events that involve high wind speeds over long time periods and large areas, 

such as hurricanes, waves (and with them, Stokes drift magnitudes) can grow very large. Swell 

propagating out from these areas into more distant areas can make a wind speed approximation 

of Stokes drift especially inaccurate. Additionally, Stokes drift can become a much more 

significant fraction of the total surface transport. This effect, at least in the GoM, can last for a 

few days after the storm has exited the basin.

The results of this section have been submitted for publication in JGRO-O in the special 

issue on oil spills. They have been accepted with minor revisions (Clark, et al., 2016).

4.4. PARAMETERIZATION OF OIL SURFACE ROUGHNESS

A modification of surface roughness was added as a new option in the Coupled-OceanAtmosphere-Wave-Sediment-Transport (COAWST) model. This is largely accomplished 

through a modification of the atmospheric model (WRF) to have a new option for calculating 

surface turbulent fluxes. One further complication is that the surface current data had to be 

passed out of the ocean model (ROMS) to be used in the flux calculations. 

The model considers contributions to surface roughness from three types of surface features 

that include an aerodynamically smooth surface (Kondo, 1975; Nikuradse, 1933) capillary waves

(Bourassa, et al., 1999), and gravity waves (Smith, et al., 1992). The roughness length 

parameterization equation for oil is

0.5 2 2 2

*

0 2

*

*

0.11

s c g

w

b a

z

u g

 

   



                            

u

u

, (4-8)

where the 



s

, 



c

, 



g

terms are binary weights for the roughness length associated with an 

aerodynamically smooth surface, capillary waves, and gravity waves; 





is the molecular 

viscosity of air, represents the oil damping effects on capillary waves and short gravity waves, 

and here is set to 0.25 when oil is present; b = 0.019 is a dimensionless constant that is 

determined from laboratory observations (Bourassa, et al., 1999); 





is surface tension, which is 

determined by the SST; 



 w

is water density, and g is gravitational acceleration. This is different 

than the anisotropic roughness length that was proposed, but is mathematically identical in the 

coordinate system we apply, and the algorithm is faster.

The following modifications have been made from Bourassa (2006). The value of is no 

longer restricted to zero (no capillary waves) and one, which was well suited for wave tanks, but 

not for much larger spatial scales. The new value of represents an average over enough space











c



c
86

and time for a smooth transition from an aerodynamically smooth surface to a rough surface.

Similarly, =0 applies only to wind waves.

For water surface, = 0 and = 0 for Ueff < Ulim; otherwise, (4-9)



ct anh0.4Ue f fUl i m

3  

        

, (4-10)



gt anh0.2Ue f fUl i m

3  

        

. (4-11)

For oil surface, = 0 and = 0 for Ueff < Ulim; otherwise, (4-12)



ct anh0.4Ue f fUl i m

3  

        

, (4-13)



gt anh0.3Ue f fUl i m

3  

        

, (41-14)

where



Ue f f u* l nz z  01z,z0   ,L/Kv , 



z,z0  ,L

is a stability term, and



Kv

and L are von 

Karman’s constant and the Monin-Obukhov stability length, respectively.

For a water surface Ulim = 1.0 m s

-1

, and for the Deepwater Horizon oil-covered region, the 

value of Ulim is 7.0 m s

-1

. This value of Ulim for oil was crudely chosen based on scatterometer 

observations of the surface roughness. The value of Ulim is not general, and will vary greatly 

depending on the oil conditions. Further work could be done to better estimate the values of Ulim

and ; however, that was beyond the scope of this work. 

We had hoped to incorporate this modification into COAPS’s modeling capability for the 

GoM; however, the code requires wave information, and that was not available in a coupled 

sense (which is important). We put a great deal of people-hours into solving that problem, but in 

the end found it was easier to modify the public version of the COAWST model. One advantage 

of using this model is that some of the modeled sediment in the water column can be given the 

characteristics of oil, allowing the simulated oil to track and interact with the surface roughness 

parameterization. 

This work is being carried on through separate funding. Under that funding we have found 

that the ocean currents substantially modify air and sea coupling, and expect oil and waves to 

have greater impacts over most of the ocean. That will be investigated within the next year, and 

will lead to a doctoral dissertation and presumably to peer-reviewed publications that will 

acknowledge BOEM support.



g



c



g



c



g
87

5. SURFACE OIL TRANSPORT CONSTRAINTS AND BARRIERS

Mesoscale and submesoscale features in the ocean are commonly associated with large 

velocity shear and anisotropic mixing along jets (e.g., Smith & McWilliams, 2003), forming near 

barriers to cross-boundary transport. In some places in the ocean, such as the northern GoM with 

its large coastal freshwater input and seasonal horizontal thermal gradients, such features may be 

evident in surface hydrographic properties. Knowledge of the roles of such features in 

constraining transport may serve to improve forecasting of oil drift, given that they may be 

observable from satellite imagery at certain times. The goals of this part of the research program 

are to: (1) identify features that may serve as boundaries or constraints to surface transport in the 

northern GoM, (2) analyze the transport of oil across these boundaries, and (3) identify and 

quantify critical strengths of such boundaries for serving as barriers to surface oil transport.

5.1. BOUNDARIES INFERRED FROM INSTANTANEOUS SSH AND SURFACE 

VELOCITY FIELDS

A material boundary separates different bodies of fluids, as the boundary always marks the 

same fluid material as the fluid evolves in time. Material boundaries are dependent on the time 

evolution of the flow field and may not be identifiable in instantaneous oceanic fields. For slowly 

evolving oceanographic features, such as mesoscale eddies and fronts, there may be some 

correspondence between the material boundaries and structures in other scalar fields, such as 

strong gradients (indicative of fronts) in surface salinity, temperature, height, and optical 

properties. The purpose of this exercise is to explore whether or not useful information on 

barriers to oil transport can be inferred from observations of these variables, which would aid in 

prediction of surface oil drift.

Of these variables, satellite observations of temperature, salinity and color may be of limited 

use for identifying potential boundaries for surface oil transport because the very presence of 

surface oil alters the remotely sensed observations of these oceanic properties (e.g., Leifer, et al., 

2012). Thus, SSH is at present the most likely observable field that can yield information about 

potential constraints on surface oil transport. SSH measurements do not provide fields of 

observations at a single time, however, as these measurements are made along satellite altimeter 

tracks with wide spatial (hundreds of kilometer) and temporal (O [10 days]) spacing. These 

observations from multiple satellites may be used to produce SSH maps through statistical 

methods (e.g., Leben, 2005) or by assimilation into a numerical model. The resulting maps may 

be an estimate of the true SSH field at a particular time, but they are really created from 

observations collected over a time window. For near-real-time products, this window is 

obviously lagged in time. The benefit of using model SSH analyses is that related velocity data 

are also available for concurrent analysis. In this section, boundaries inferred from SSH fields 

from a data assimilating Hybrid Coordinate Ocean Model (HYCOM, experiment, 31.0, are 

compared to SAR-derived observations of oil coverage. In particular, a range of dates in midMay 2010 during which oil was transported to the region of the Loop Current forming a “tiger 

tail”-like structure (e.g., Olascoaga & Haller, 2012) is analyzed. The focus of this section is on 

maps of scalar variables that are synoptic in time.
88

5.1.1. Techniques based on Instantaneous Velocity Fields

Many front definition techniques have been developed as methodologies for tracking 

mesoscale eddies, i.e., closed-contour features (e.g., closed contours of relative vorticity fields in

McWilliams, 1990; the curvature center method in de Leeuw & Post (1995); the winding-angle 

method in (Sadarjoen & Post, 2000); the threshold-free identification method in Chelton, et al. 

(2011). The oceanic fronts analyzed here may be closed- or open-contour features. One such 

front definition is based on the Okubo-Weiss Parameter (Okubo, 1970; Weiss, 1991), which

assesses deformation (shear and strain) and relative rotation of the flow field. This method is the 

most effective for identification of eddies because in eddies, rotation dominates the deformation 

terms making the parameter negative (Figure 5-1). Thus negative regions in the Okubo-Weiss 

parameter field are possible eddies (the parameter does not distinguish the sense of rotation in 

the eddies). Additional steps are needed to isolate eddies from other vorticity-dominated features. 

Usually a threshold value is chosen (Calil, et al., 2008). The major drawback of this 

methodology is the field is noisy. It may not be useful for identification of oceanic fronts other 

than those associated with eddies. 

Figure 5-1. Okubo-Weiss parameter (x109) estimated from the HYCOM surface layer velocities

The blue contours are SSH isolines at 0, 0.25 and 0.5 m. The black contour is the surface oil boundary 

from SAR-TCNNA data.
89

A method of identifying mesoscale eddies based on relative vorticity was proposed by 

McWilliams (1990). The eddies are identified by the last closed contour around the local 

vorticity extremum. The method has not been used for identification of open-contour oceanic 

fronts. Although, the field is smoother than the Okubo-Weiss parameter and fronts are more 

evident (Figure 5-2). Oceanic velocity fronts are associated with extrema of the relative vorticity 

gradient matching the maximum of the flow speed. Potentially, these fields can be employed for 

front tracking if velocity fields are available. Application of the potential vorticity to study GoM 

mesoscale processes is illustrated in (Le Hénaff, et al., 2012).

Figure 5-2. Relative vorticity (x10

5

s

-1

) estimated from the HYCOM surface layer velocities

The blue contours are SSH isolines at 0, 0.25 and 0.5 m. The black contour is the surface oil boundary 

from SAR-TCNNA data.

5.1.2. Techniques based on SSH

SSH-based techniques have been successfully used in the past to track the mesoscale features 

in the GoM. Because the LC and its associated LCEs are approximately in geostrophic balance, 

fixed SSH contour levels will very nearly follow streamlines in the flow. Past studies have used 

SSH contours for tracking of the LC front in SSH fields assuming that the LC front coincides 

with a single streamline. Other LC tracking methods have also been proposed (e.g., AndradeCante, et al., 2013; Lindo-Atichati, et al., 2013). The most basic and simple, yet reliable 

technique for LC and LCE tracking was suggested by Leben (2005) where the LC and LCE 
90

fronts are tracked using the 0.17-m contour in demeaned SSH fields. Demeaned fields are 

calculated by subtracting the spatial mean from each daily SSH field, which is necessary to 

remove bias in the surface elevation fields associated with different reference surfaces and 

seasonal height variations due to upper-ocean warming and cooling. The 0.17 m threshold for LC 

identification seems to be natural because under the assumption of geostrophic balance, the SSH 

contours are streamlines of the instantaneous geostrophic flow. At the same time, it is not 

obvious that one particular SSH contour can precisely follow the LC front, especially in light of 

the fact that the pathline of a fluid parcel can cross streamlines in a time dependent flow field. 

This implies a possibility of discrepancies in LC metrics derived from different definitions of the 

LC front. In order to test the sensitivity of LC statistics to alternate definitions of the LC front, 

the LC is also tracked using a discrete Kalman Filtering algorithm (Kalman, 1960) that identifies 

frontal positions using a combined analysis of SSH and SSH gradient fields (Figure 5-3). 

In the Kalman Filtering LC tracking algorithm, the supposition is made that the LC front 

closely follows the high-velocity core of the LC. Under a geostrophic assumption, the maximum 

gradient of the SSH closely follows the core of the LC and thus should be a more natural 

candidate to use as a criterion for eddy identification. In the following tracking methodology, the 

LC is defined employing the discrete Kalman Filtering algorithm (Kalman, 1960) to obtain the 

frontal position from two model fields (see details in Appendix A and its associated Figure A1). 

Two model fields provide information for an a priori estimate and correction (referenced as a 

“measurement” in the traditional application of correcting model prediction) to obtain the final 

frontal location (a posteriori estimate). In this application, the SSH field provides the first guess 

of the LC front location (a priori estimate). The second field is the SSH gradient. The gradient 

field renders the information about the dynamics of the upper ocean and is used as a reference 

field to correct the first-guess approximation of the LC front from the SSH field. In theory, 

following the maximum gradient would delineate the location of the core and frontal position of 

the LC. However, the SSH gradient field, as with many dynamic fields, has local extrema and 

cannot be objectively tracked to draw a single continuous contour from the Yucatan Channel to 

the Straits of Florida. This algorithm uses information about the LC location from two fields, 

“deciding” at every step whether to trust the first or the second field more. Any other 

oceanographic field capable of capturing mesoscale structures can be used as a “measurement.” 

For instance, the relative vorticity could be a good potential alternative as it highlights dynamical 

fronts. The SSH gradient fields have been chosen to demonstrate the utility of extra information 

derived from the original SSH field that may be easily derived from satellite observations. 

Validation of the technique is provided in Dukhovskoy et al. (2015). 
91

Figure 5-3. Gradient of the SSH from HYCOM

The grey contour is 0.17 m. The red contour is the LC front obtained by Kalman Filtering technique. The 

methodology corrects the 0.17m contour location to match it closer to the maximum SSH gradient. The 

yellow contour is the surface oil boundary from SAR-TCNNA data.

From inspection of the evolution of the oil coverage maps and analysis of velocity and SSH 

fields from a data assimilating model analysis product (which can be thought of as a dynamical 

mapping system for observational data, an alternative to a purely statistical mapping approach), 

it is clear that none of these analysis methods revealed a barrier to the oil drift during this critical 

stage of the oil spill during which the tiger tail formation occurred. It is during this time that oil 

first entered the Loop Current. Although the oil in the Loop Current was observed to dissipate 

rapidly, it is worth noting that the warm environment likely sped the evaporation or 

biodegradation of oil. Had this event happened during the winter months, this oil entrained in the 

Loop Current may have survived a much longer transit.

The surface oil not only moved toward the Loop Current region it crossed over extrema in 

the Okubu-Weiss parameter (Figure 5-1), crossed local extrema of vorticity of both positive and 

negative sign (Figure 5-2), and crossed the Loop Current frontal position defined from SSH 

fields (Figure 5-3). It is possible that the LC frontal region may indeed serve as a barrier to 

transport into the interior of the LC SSH maximum region, but there seemed to be a conduit for 
92

oil to enter the LC rather than any barriers to this transport. Such a conduit mechanism is 

explored in the following sections which analyze the time-dependent velocity field. 

5.2. MATERIAL BOUNDARIES BASED ON LAGRANGIAN COHERENT 

STRUCTURES

Lagrangian Coherent Structures (LCSs) are becoming a tool of much interest for the study of 

advection of fluids. LCSs are material surfaces in that fluid particles existing along an LCS will 

remain on the LCS as it evolves over time. LCSs are lines (in two-dimensional cases) with often 

strongly repelling or attracting nature, and thus have been shown to be useful in forecasting the 

evolution of patterns in movement of material (oil) on the ocean surface (Olascoaga & Haller, 

2012). The flux across an LCS is negligible (comparable to errors in the computation), so they 

will serve as barriers in the transport of fluid from one side to the other (Shadden, et al., 2005). 

However, because these, like all oceanic material boundaries, evolve in time, their identification 

at any one instant may be of limited use in establishing barriers to eventual oil transport. Because

LCSs move as the flow field evolves, LCSs should not be thought of as barriers in relation to 

particular geographic regions because these LCS may eventually migrate across regions of 

concern as the flow field evolves, such that a particular location initially on one side of an LCS 

may be on the other side at some later time. Additionally, identification of LCSs requires 

knowledge of the flow field evolution over some time interval, which cannot be known from 

observations alone, particularly for future times. 

LCSs are identified through Lagrangian analysis of the time-dependent flow field. A 

commonly used approach is based on computation of Finite Time Lyapunov Exponent (FTLE) 

fields, where local ridges approximately define hyperbolic (attracting or repelling) LCSs. Other 

analysis methods have been developed to additionally detect other coherent structures in the flow 

field, including elliptic LCSs (Lagrangian vortices) and parabolic LCSs (Lagrangian jet cores) 

(e.g., Haller, 2015). Because the attracting or repelling structures are of particular interest for 

analyzing potential barriers for movement of oil on the ocean surface, the FTLE method is used 

here. In essence, the FTLE is a scalar value σt

T

(x) that characterizes the rate of stretching of 

trajectories about a trajectory of a point x over the time interval [t, t+T]. Details and derivation of 

FTLEs are omitted here, but can readily be found in texts such as Shaden (2012). Practically, the 

FTLE scalar field is computed from a time-dependent velocity field (here, surface currents from 

an ocean model). A flow map is determined by computing trajectories for a large number of 

initial locations distributed (commonly along a grid) in space over the time interval [t, t+T]. 

From the gradient matrix of the flow map is computed the right Cauchy-Green strain tensor. The 

natural logarithm of its largest eigenvalue, divided by the time of integration, yields the FTLE 

(actually, each eigenvalue corresponds to an FTLE, but the largest FTLE associated with the 

largest eigenvalue is usually considered as the FTLE without distinction). The ridges of the 

FTLE field represent the LCSs. Note that FTLEs may be computed either forward in time or 

backward in time to determine regions of fastest separation or repelling (forward in time) or 

attraction (backward in time). These LCSs consisting of ridges in the forward-in-time FTLE 

field, indicating repelling regions, are analogous to stable manifolds in time invariant flow fields, 

and attracting LCSs are analogous to unstable manifolds.
93

5.3. TRANSPORT OF SURFACE OIL IN THE PRESENCE OF MATERIAL 

BOUNDARIES 

Oil transport pathways and their relationships to structures in various oceanographic scalar

fields are investigated using the output from a simplified oil drift model driven by a data 

assimilative hydrodynamic model hindcast. This approach is used so that the Lagrangian 

pathways of oil parcels can be analyzed in conjunction with a variety of oceanic fields that 

cannot be readily determined from observations alone (e.g., salinity and properties of the timedependent flow field). 

A simulation of the oil spill was conducted using the simple surface oil drift model described 

in Section 3.2.4a. This simulation used the HYCOM 1/25° GoM Reanalysis (50.1) as the driving 

hydrodynamic model. This HYCOM simulation has only recently been made available (late 

March 2016) so its inclusion in the different phases of this overall project is limited. 

Nevertheless, its updated parameterization of riverine input yields more realistic plumes and 

salinity gradients that made it a useful tool for this study of oceanic boundaries and surface 

transport so it is employed here. This HYCOM simulation was forced by CFSR-derived 

atmospheric variables. Several oil spill simulations were conducted and compared to SARderived time-composited oil coverage maps (Figure 5-4) and the best agreement was achieved by 

advecting the oil in the model surface currents with no additional wind drift term and an oil halflife of 168 hours with maximum oil life of three half-lives. Additional parameters of this 

simulation are shown in Table 5-1.

Table 5-1.

Oil model parameters for the half-life time estimation

Parameter Oil Transport Pathways

Np 50

Th 168

Tmax 504

Rs 20

dt 1

 10*

Cw 0

 not used

HYCOM analysis 50.1

Surface Winds not used

*

enhanced eddy diffusion within 50 km of the source
94

Figure 5-4. SAR TCNNA compared to simulation results

Left: Fraction of the time period 24 April 2010 1200UTC to 14 July 2010 0000UTC that each 1/20° x 1/20° 

bin had oil present from the numerical simulation (left) and SAR TCNNA analyses (right).

From the HYCOM GoM 50.1 surface velocity fields, FTLE maps were computed over the 

time intervals [T-3 days, T] during the time period T ε [28 April 2010, 9 July 2010]. The FTLE 

maps were constructed from flow maps computed both forward in time and backward in time. 

From these FTLE maps, LCS positions corresponding to the most prominent ridges in the FTLE 

maps were estimated. Ridge identification is a non-trivial task, and multiple algorithms have to 

be employed for this purpose (e.g., Shaden, 2012). For this study, a ridge extraction procedure

(as described below) using image analysis tools was implemented. First, an FTLE map is 

converted to an image of intensity values, mapping the range of FTLE values onto [0,1]. This 

image then undergoes adaptive contrast enhancement, where the FTLE values are adjusted to fit 

a uniform distribution separately within 24 geographic regions (Figure 5-5). Next, this adjusted 

image is interpolated from the native 1/25° resolution to a finer resolution (1/100°) grid. This 

map is then converted to a black-white image where values exceeding an 85th percentile 

threshold value are retained. A skeletonization procedure is then applied to this thresholded 

image to extract approximate ridge locations marking LCS boundaries (Figure 5-2 – bottom). 

LCSs computed from the forward-in-time FTLE maps indicate regions that are repelling in 

nature, and LCSs from the backward-in-time FTLE maps indicate regions that are attracting. For 

convenience, these will henceforth be referred to as “repelling” and “attracting” LCSs. Though 

the skeletonization occasionally produces undesirable artifacts in the boundaries (for example, 

multiple short connecting lines between proximate boundaries or disconnected line segments 

along a ridge), the procedure works well for the purposes here and has reasonably fast 

computational performance.
95

Figure 5-5. Example of the LCS ridge extraction procedure

Top Left: Backward-in-time FTLE map for the three-day time period ending on 21 June 2010. Top Right: 

The FTLE values remapped to the range [0,1] and processed with adaptive contrast enhancement. 

Bottom: Prominent ridges extracted from the adjusted FTLE map representing “attracting” LCSs, 

analogous to unstable manifolds.

Though LCS boundaries do not necessarily correspond to features seen in instantaneous 

oceanographic fields due to their time-dependence, certain scalar fields can at times be seen to 

align with attracting LCS boundaries. This may be due to convergence of oceanographic 

properties along these attracting LCSs. For example, regions of high magnitude salinity gradient 

are many times aligned along attracting LCS boundaries offshore of the Mississippi Delta 

(Figure 5-6). This gives the appearance that these high salinity gradient regions may be serving 

as barriers to surface oil transport. However, dynamically, it is not the salinity itself that serves 

as a barrier but rather features in the flow field, evident by LCSs. It is not clear whether these 

convergent zones in the flow field are formed by buoyancy driven circulation associated with the 

riverine water, or if low salinity water advected along these structures, which may be linked to 

mesoscale circulation features, give rise to the salinity gradients. Kourafalou & Androulidakis 

(2013) investigated this issue and suggested that the river induced salinity front and associated 

different regimes of the buoyancy-driven circulation restrained and facilitated the onshore 

transport east-west of the Mississippi Delta. The relative persistence of a high-magnitude salinity 

gradient in the vicinity of the Mississippi River Delta and associated LCS along which oil tends 

to congregate in the simulation supports this theory. Alignment of weaker salinity gradients with 
96

attracting LCSs and oil convergence further offshore appears to be more the result of advection 

of low salinity water and oil along deepwater mesoscale circulation features.

Figure 5.6. Surface salinity gradient (PSU/m, color map) from the HYCOM GoM Reanalysis (50.1)

With attracting (white) and repelling (gray) LCS positions (estimated from FTLE fields computed from flow 

maps over the previous three-day intervals) overlaid. Simulated surface oil positions are shown in black.
97

Figure 5-6. Continued

Mesoscale circulation features are more clearly seen from SSH fields. This is also the 

oceanographic variable that is most readily obtainable from satellite observations, having 

potential for forecasting future movement of surface oil. Strong currents, such as the Loop 

Current and along boundaries of eddies, are revealed by the gradient of SSH. Inspection of SSH 

gradient maps with LCS boundaries overlaid (Figure 5-7) show that while LCS boundaries are 

commonly found in the vicinity of large SSH gradient magnitudes, instantaneous flow patterns 

inferred from the SSH gradient do not align with the LCS boundaries. This is due to the nonsteady nature of the currents. It is clear that surface oil tends to converge and be advected along 

attracting LCSs. 
98

The evolution of a “tiger tail”-like feature in the oil distribution has been a topic of much 

discussion (e.g., Olascoaga & Haller, 2012). This feature formed in mid-May as oil was advected 

southeastward along the periphery of a small cyclone toward the vicinity of the anticyclonic 

circulation of the Loop Current (or its recently detached eddy). In this simulation, the oil 

interacts with this cyclone forming the tiger tail somewhat later, beginning at the end of May 

(Figure 5-7). This oil was advected along an attracting LCS structure (which the algorithm 

employed here represents as discontinuous lines shown in purple in Figure 5-7). Another such 

event was simulated maturing on 30 June. Here, oil can be seen converging along a clearly 

defined attracting LCS stretching from near the DwH site toward the LC region.

Figure 5-7. SSH gradient from the HYCOM GoM Reanalysis (50.1)

Attracting (white) and repelling (gray) LCS positions (estimated from FTLE fields computed from flow 

maps over the previous three-day intervals) overlaid. Simulated surface oil positions are indicated in 

black.
99

Figure 5-7. Continued

There is a prevalence of attracting LCSs in the region of (but not necessarily aligned with) 

strong SSH gradients This suggests that the flow field that forms these attracting structures 

actually favors entrainment of oil into the strong currents associated with the large magnitude 

SSH gradients. What prevented more oil from finding its way into the Loop Current or its large 

eddies during the Deepwater Horizon oil spill was the fact that the flow field within the DeSoto 

Canyon region favors retention. Eddies within the region are partially topographically 

constrained, thus remaining in the region over long periods entraining surface material in their 

localized circulation. Interaction of these eddies with eddies or the Loop Current to the south 

intermittently forms pathways by which the surface oil may be advected to the Loop Current. 

The formation of these pathways is clearly shown by attracting LCSs stretching southward, such 

as that on 30 June in the simulation (Figure 5-7). 
100

Though the attracting LCSs show potential conduits for transport, they also serve as barriers 

to transport across them. It is therefore reasonable to ask whether these boundaries can be 

inferred from readily obtainable observational data. However, inspection of Figure 5.7 reveals 

that these LCS, while prevalent along regions of large SSH gradient, also occur in regions of 

small SSH gradient with oil advected along them. It is necessary to have velocity data to estimate 

LCS locations. Beron-Vera, et al. (2008) used geostrophic velocity fields derived from satellite 

SSH to compute LCSs. This technique could possibly be used in the context of an oil spill. 

However, surface oil drift is significantly impacted by ageostrophic velocity components that 

cannot be deduced from altimetry. Further, altimetry-derived velocity fields are lagged in time. 

Thus, hydrodynamic forecast models remain the best tools for forecasting movement of surface 

oil, either through analysis of the kinematics of the flow field or through direct application in an 

oil spill forecast model

.
101

6. FUTURE RESEARCH NEEDS

A major concern of oil spill response is tracking actionable oil; meaning oil that is amenable 

to treatment or recovery. As the foregoing report has shown, WAMOST project has achieved a 

number of peer-reviewed results that can improve the quantitative observation of surface oil 

through remote sensing and has generated datasets that will facilitate continuing research on 

potential ecological impact and numerical assessment of oceanographic forcing factors that 

affect the advection of surface oil. This process has also identified promising directions for 

further research to improve our ability to quantify risk associated with oil spills including their 

fate and effects. These are briefly outlined below following the approximate order of 

presentation used in this report.

6.1. CONTINUED PUBLICATION

All of the task groups have additional results to publish in the peer-reviewed record. As 

academics, the PIs are committed to continuing this process beyond the period of performance 

for the BOEM project. All future papers that are developed from the work performed under 

WAMOST will acknowledge BOEM support and will be submitted to the WAMOST technical 

officer for review prior to publication.

6.2. REMOTE SENSING OF SURFACE OIL

The fundamental difficulty is accurate characterization and quantification of surface oil state 

(i.e., emulsified or not), oil volume, and oil thickness in the lab and in the field. To date, all 

validation of remotely sensed oil maps is through indirect measures (e.g., burning, visual 

interpretation). More effort needs to be dedicated to overcome such a fundamental difficulty.

For SAR applications, new technologies initiated during and after the DwH spill have the 

potential to greatly increase the resolution and sensitivity of SAR observations using C-band 

radar. Examples of new satellite platforms include the Canadian RADARSAT Constellation, 

which produces quad-polarized images at resolutions as fine as 1 x 3 m (Canadian Space 

Agency, 2015) and the European Sentinel-1, which is currently undergoing operational testing 

(ESA. Earth Online, 2016). Adaptation of TCNNA and OEDA algorithms to the new image 

formats will be required in preparation for any oil spill applications. Research in this direction 

could capitalize on the facilities at the National Oil Spill Response Test Facility – OHMSETT. 

Additionally, the ongoing release of oil from the MC020 lease block (Taylor Energy site) offers 

a real-world test-bed for remote sensing validation and calibration.

For optical remote sensing applications, two challenges still exist that require significant amount 

of future effort: 1) determine oil presence or absence, including differentiating oil from other 

look-alikes (Trichodesmium, Sargassum); 2) determine oil volume (or thickness) from multiband, coarse-resolution data. For the oil’s heterogeneity and for its lack of spectral signatures in 

the visible, these still require future effort to improve the oil quantification. 
102

6.3. FORCING FIELDS AND VALIDATION

The structure of the near-surface velocity profile and the related dynamics still remain a 

fundamental knowledge gap for simulating surface oil transport. Results from this project 

demonstrate the differences in how hydrodynamic models represent surface velocity given their 

numerical schemes and vertical grid discretization, as well as the impact Stokes drift can have in 

inducing current shear at the surface. There is a lack of velocity observations within the upper 

few meters of the ocean surface making this a difficult area of study. However, recent and 

upcoming advances in observing this surface layer shear provide great opportunities for further 

research aimed at improving understanding the near-surface physics, which should in turn allow 

improvements to be made in estimating the “true” surface current from upper layer model 

currents when applied to oil drift simulations.

During this project, advances were made in characterizing oil thickness. However, the 

motion of this thicker oil, including emulsions of different densities, likely differs from thin oil. 

To improve prediction capabilities and better direct response efforts toward the more harmful 

thick oil, the drift of thick oil emulsions should be addressed in future research tied with 

developing better understanding of near-surface ocean velocity structure and interfacial stresses 

between air, oil and water under different wind, wave and current regimes. 

The DwH oil spill provided valuable data for understanding the weathering of oil. However, 

this event occurred during a time of year in which surface waters and air temperatures of the 

Gulf were very warm and incoming solar radiation was at its maximum. During the event, there 

was much concern about oil becoming entrained in the Loop Current, with some groups showing 

simulations of possible transport pathways through the Florida Straits and into the Atlantic 

Ocean. Though oil was entrained in the Loop Current, the relatively rapid weathering of oil 

prevented significant amounts of actionable or observable oil from traveling these great 

distances. The DwH oil was also light compared to some well-studied oil spills in colder 

climates, such as the Exxon Valdez spill, in which oil remained for much longer periods of time. 

Thus, it is not clear how the fate of the oil released during the DwH spill would have been 

different had the event happened during the winter when northern Gulf water and air 

temperatures can be more than 15° C colder, and shortwave radiation is only about one third of 

that in the warm months. Further research should be conducted to understand potential fates of 

different oil types when released in cold climate conditions.

Results from this research and recently published research by others (e.g., Kourafalou &

Androulidakis, 2013) show that strong salinity gradients may serve as constraints to oil transport. 

These are poorly observed, however. Further research should be conducted to develop better 

techniques of characterizing these features from remotely sensed data, including recent surface 

salinity observations from the SMAP (Soil Moisture Active Passive) satellite, and in the use of 

other radiometric satellite sensors to detect small-scale oceanographic features in the presence of 

surface oil.
103

6.4. SURFACE TRANSPORT CONSTRAINTS

The actual surface current is likely substantially greater than the velocity in the top layer of 

the ocean model. Atmospheric models solve this problem by parameterizing the wind profile in 

the layer close to the surface. A similar approach can be used in an ocean model, but 

observations of surface currents (and ideally near surface shear) are needed to determine which 

parameterization is appropriate and the conditions for which it is appropriate. In the last year in 

situ and remotely sensed technologies have been developed to measure the surface current. These 

technologies should be used to gather the data needed to determine appropriate 

parameterizations, and a model of the current shear near the ocean surface should be added to the 

models.

Two-way coupled ocean-wave-atmosphere models now exist and could be used to examine

the importance of this coupling in assessing flow transport. Preliminary work suggests that this 

coupling has non-negligible impacts on transport.
104

7. LITERATURE CITED

Abascal, A. J. et al., 2009. Calibration of a Lagrangian transport model using drifting buoys 

deployed during the Prestige oil spill. J Coastal Res, 25(1), pp. 80-90.

Aeppli, C. et al., 2013. Recurrent oil sheens at the Deepwater Horizon disaster site fingerprinted 

with synthetic hydrocarbon drilling fluids. Environ Sci Technol, 47(15), pp. 8211-8219.

Agarwal, P. K. & Phillips, J. M., 2006. On bipartite matching under the RMS distance. ON, 

Queen's University, p. 143–146.

American Society of Civil Engineers Committee on Modeling Oil Spills, 1996. State-of the-art 

review of modeling transport and fate of oil spills. J Hydraul Eng, 122(11), p. 594–609.

Andrade-Cante, F., Sheinbaum, J. & Zavala Sanson, L., 2013. A Lagrangian approach to the 

Loop Current eddy separation. Nonlinear Processes in Geophysics, 20(1), pp. 85-96.

Androulidakis, Y. S. & Kourafalou, V. H., 2013. On the processes that influence the transport 

and fate of Mississippi waters under flooding outflow conditions. Ocean Dynamics, 63(2), p. 

143.

Atlas, R. et al., 2011. A cross-calibrated, multiplatform ocean surface wind velocity product for 

meteorological and oceanographic applications. Bull Amer Meteor Soc, 92(2), pp. 157-174.

Barron, C. N. et al., 2004. Sea surface height predictions from the Global Navy Coastal Ocean 

Model (NCOM) during 1998-2001. J. Atmos. Oceanic Technol, 21(12), p. 1876–1894.

Barron, C. N. et al., 2006. Formulation, implementation and examination of vertical coordinate 

choices in the global Navy Coastal Ocean Model (NCOM). Ocean Modeling, 11(3-4), p. 

347–375.

Beron-Vera, F. J., Olascoaga, M. J. & Goni, G. J., 2008. Oceanic mesoscale eddies as revealed 

by Lagrangian coherent structures. Geophys Res Lett, 35(L12603).

Bleck, R., 2002. An oceanic general circulation model framed in hybrid isopycnic Cartesian 

coordinates. Ocean Model, 4(1), p. 55–88.

Boehm, P. D. & Carragher, P. D., 2012. Location of natural oil seep and chemical fingerprinting 

suggest alternative explanation for deep sea coral observations. P Natl Acad Sci, 109(40), p. 

E2647.

Bourassa, M. A., 2006. Satellite-based observations of surface turbulent stress during severe 

weather. In: W. Perrie, ed. Atmosphere Ocean Interaction, Vol. 2. Southampton, UK: Wessex 

Institute of Technology Press, pp. 35-52.

Bourassa, M. A., Vincent, D. G. & Wood, W. I., 1999. A flux parameterization including the 

effects of capillary waves and sea state. J Atmos Sci, 56(9), p. 1123–1139.
105

Brekke, C. & Solberg, A. H., 2005. Oil spill detection by satellite remote sensing. Remote Sens 

Environ, 95(1), pp. 1-13.

Bryant, R. B., Lugo, J., Cordova, C. & Salvador, A., 1991. Physiography and bathymetry. In: A. 

Salvador, ed. The Geology of North America, The Gulf of Mexico Basin. Boulder: J Geol Soc 

Am, p. 561.

Calil, P. H. P., Richards, K., Jia, Y. & Bidigare, R., 2008. Eddy activity in the Lee of the 

Hawaiian Islands. Deep Sea Res, 55(10-13), p. 1179–1194.

Camilli, R. et al., 2011. Acoustic measurement of the Deepwater Horizon Macondo well flow 

rate. P Natl Acad Sci, 109(50), p. 20235–20239.

Canadian Space Agency, 2015. Satellite Characteristics. [Online] Available at: http://www.asccsa.gc.ca/eng/satellites/radarsat/radarsat-tableau.asp [Accessed 15 05 2016].

Castanedo, S. et al., 2006. The Prestige oil spill in Cantabria (Bay of Biscay). Part I: Operational 

forecasting system for quick response, risk assessment and protection of natural resources. J 

Coastal Res, 22(6), p. 1.

Chang, C. C., Hwang, S. M. & Buehrer, D. J., 1991. A shape recognition scheme based on 

relative distances of feature points from the centroid. Pattern Recognition, 24(11), p. 1053–

1063.

Chassignet, E. P. et al., 2006. Generalized vertical coordinates for eddy-resolving global and 

coastal ocean forecasts. Oceanography, 19(1), p. 20–31.

Chassignet, E. P. et al., 2007. The HYCOM (HYbrid Coordinate Ocean Model) data assimilative 

system. J Mar Syst, 65(1-4), p. 60–83.

Chassignet, E. P., Smith, L. T., Halliwell, G. R. & Bleck, R., 2003. North Atlantic simulation 

with the HYbrid Coordinate Ocean Model (HYCOM): Impact of the vertical coordinate 

choice, reference density, and thermobaricity. J Phys Oceanogr, 33(12), p. 2504–2526.

Chelton, D., Schlax, M. G. & Samelson, R. M., 2011. Global observations of nonlinear 

mesoscale eddies. Prog Oceanogr, 91(2), p. 167–216.

Chen, S. et al., 2003. COAMPSTM Version 3 Model Description - General Theory and Equations, 

Monterey: NRL.

Clark, M., Heath, N., Bourassa, M. A. & Chassignet, E., 2016. Quantification of Stokes drift as a 

mechanism for surface oil advection in the Gulf of Mexico. J Geophys Res.. J Geophys Res, 

Volume Forthcoming.

Crooke, E. et al., 2015. Determination of sea-floor seepage locations in the Mississippi Canyon. 

Mar Petrol Geol, 59(January 2015), p. 129–135.
106

Cummings, J. A., 2005. Operational multivariate ocean data assimilation. Quart J Royal Met 

Soc, 131(613), p. 3583–3604.

Daling, P. S. et al., 2014. Surface weathering and dispersibility of MC252 crude oil. Mar Pollut 

Bull, 87(1-2), p. 300–310.

Daneshgar Asl, S. et al., 2016. Chronic, anthropogenic hydrocarbon discharges in the Gulf of 

Mexico. Deep Sea Research Part II: Topical Studies in Oceanography. Deep Sea Research 

Part II: Topical Studies in Oceanography, Volume 129, p. 187–195.

Daneshgar Asl, S., Dukhovskoy, D., Bourassa, M. & MacDonald, I. R., 2017. Hindcast modeling 

of oil slick persistence from natural seeps. Remote Sens Environ, Volume 189, p. 96–107.

Daoudi, M. et al., 1999. Shape distances for contour tracking and motion estimation. Pattern 

Recogn, 32(7), p. 1297–1306.

De Beukelaer, S. M., MacDonald, I. R., Guinnasso, N. L. & Murray, J. A., 2003. Distinct sidescan sonar, RADARSAT SAR, and acoustic profiler signatures of gas and oil seeps on the 

Gulf of Mexico slope. Geo-Mar Lett, 23(3-4), p. 177–186.

de Leeuw, W. C. & Post, F. H., 1995. A statistical view on vector fields. In: M. Gobel & H. 

Muller, eds. Visualization in Scientific Computing: Rostock. Rostock, Germany: SpringerVerlag, Wien, p. 53–62.

de Miranda, F. P. et al., 2004. Analysis of RADARSAT-1 data for offshore monitoring activities 

in the Cantarell Complex, Gulf of Mexico, using the unsupervised semivariogram textural 

classifier (USTC). Can J Remote Sen, 30(3), pp. 424-436.

D'Sa, E. J. & Korobkin, M., 2008. Colored dissolved organic matter in the northern Gulf of 

Mexico using ocean color: seasonal trends in 2005. Cardiff, Wales, s.n.

D'Sa, E. J. & Miller, R. L., 2003. Bio-optical properties in waters influenced by the Mississippi 

River during low flow conditions. Remote Sens Environ, 84(4), p. 538–549.

Dubuisson, M. P. & Jain, A. K., 1994a. 2D matching of 3D moving objects in color outdoor 

scenes. In: Proceedings CVPR ’94. Proceedings of the IEEE Computer Society Conference 

on Computer Vision and Pattern Recognition. Seattle, WA, IEEE Computer Society Press 

and the Institute of Electrical and Electronic Engineers, pp. 887-891.

Dubuisson, M. P. & Jain, A. K., 1994b. A modified Hausdorff distance for object matching. 

Jerusalem, IEEE, p. 566–568.

Dukhovskoy, D. S. et al., 2015b. Characterization of the uncertainty of Loop Current metrics 

using a multidecadal numerical simulation and altimeter observations. Deep-Sea Res Pt I, 

Volume 100, p. 140–158.

Dukhovskoy, D. S. et al., 2015a. Skill metrics for evaluation and comparison of sea ice models. J 

Geophys Res Oceans, Volume 120, p. 5910–5931.
107

ESA. Earth Online, 2016. ESA. Earth Online. What is Sentinel 1?. [Online] Available at: 

https://earth.esa.int/web/guest/missions/esa-operational-eo-missions/sentinel-1 [Accessed 15 

05 2016].

Espedal, H. & Wahl, T., 1999. Satellite SAR oil spill detection using wind history information. 

Int J Remote Sens, 20(1), pp. 49-65.

Etnoyer, P. et al., 2015. Decline in condition of gorgonian octocorals on mesophotic reefs in the 

northern Gulf of Mexico: before and after the Deepwater Horizon oil spill. Coral Reefs, 5(1), 

p. 77–90.

Federal Interagency Solutions Group, 2010. Oil Budget Calculator Deepwater Horizon 

Technical Documentation. [Online] Available at: 

http://www.restorethegulf.gov/sites/default/files/documents/pdf/OilBudgetCalc_Full_HQPrint_111110.pdf

[Accessed 08 04 2013].

Fitzpatrick, P., Lau, Y., Hill, C. & Karan, H., 2010. COAMPS wind validation in the NCOM 

Intra-American Seas (AMSEAS) domain in the Gulf of Mexico during 20 June to 10 July 

2010, and the use of NCOM and COAMPS data to examine the impact of cyclones on the

Deepwater Horizon oil spill, s.l.: Geosystems Research Institute.

Garcia-Martinez, R. & Flores-Tovar, H., 1999. Computer modeling of oil spill trajectories with a 

high accuracy method. Spill Sci Technol B, 5(5/6), p. 323–333.

Garcia-Pineda, O. et al., 2013a. Detection of floating oil anomalies from the Deepwater Horizon 

oil spill with synthetic aperture radar. Oceanography, 26(2), p. 124–137.

Garcia-Pineda, O. et al., 2013b. Oil spill mapping and measurement in the Gulf of Mexico with 

Textural Classifier Neural Network Algorithm (TCNNA). IEEE J Sel Top Appl Earth Obs 

Remote Sens, PP(99), p. 2517–2525.

Garcia-Pineda, O. et al., 2010. Remote-sensing evaluation of geophysical anomaly sites in the 

outer continental slope, northern Gulf of Mexico. Deep-Sea Res Pt II, 57(21-23), p. 1859–

1869.

Garcia-Pineda, O. M. et al., 2015. Transience and persistence of natural hydrocarbon seepage in 

Mississippi Canyon, Gulf of Mexico. Deep-Sea Res Pt II, Volume 129, pp. 119-129.

Garcia-Pineda, O. et al., 2009. Using SAR images to delineate ocean oil slicks with a Texture 

Classifying Neural Network Algorithm (TCNNA). Can J Remote Sens, 35(5), pp. 411-421.

Gibbs, R. J., 1970. Circulation in the Amazon River estuary and adjacent Atlantic Ocean. 

Journal of Marine Research, Volume 28, pp. 113-121.

Green, R. E., Gould, Jr, R. W. & Ko, D. S., 2008. Statistical models for sediment/detritus and 

dissolved absorption coefficients in coastal waters of the northern Gulf of Mexico. Cont Shelf 

Res, 28(10-11), p. 1273–1285.
108

Haller, G., 2015. Lagrangian Coherent Structures. Annu. Rev. Fluid Mech, Volume 47, p. 137–

162.

Hasselmann, K. et al., 1973. Measurements of wind-wave growth and swell decay during the 

Joint North Sea Wave Project (JONSWAP), Ergänzungsheft 8-12: Deutches 

Hydrographisches Institut.

Hiester, H. et al., 2016. A topological approach for quantitative comparisons of ocean model 

fields to satellite ocean color data. Meth Oceanogr. Meth Oceanogr, Volume 17, pp. 232-

250.

Hine, A. C. et al., 2003. The west-central Florida inner shelf and coastal system: A geologic 

conceptual overview and introduction to the special issue. Mar Geol, 200(1-4), pp. 1-17.

Hodur, R. M. et al., 2002. The Coupled Ocean/Atmosphere Mesoscale Prediction System 

(COAMPS). Oceanography, Volume 15, pp. 88-98.

Hu, C. et al., Submitted 2016. Remote sensing estimation of surface oil volume during the 2010 

Deepwater Horizon oil blowout in the Gulf of Mexico: Scaling up AVIRIS observation with 

MODIS measurements. Intl J Applied Earth Observation and Geoinformation.

Hu, C. M. et al., 2003. Comparison of ship and satellite bio-optical measurements on the 

continental margin of the NE Gulf of Mexico. J Remote Sens, 24(13), p. 2597–2612.

Hu, C. et al., 2004. Linkages between coastal runoff and the Florida Keys ecosystem: A study of 

a dark plume event. Geophys. Res. Lett. 31, L15307,, Volume 31, p. L15307.

Hu, C. et al., 2011. Did the northeastern Gulf of Mexico become greener after the Deepwater 

Horizon oil spill?. Geophys Res Lett, Volume 38, p. L09601.

Huttenlocher, D. P., Klanderman, G. A. & Rucklidge, W. J., 1993. Comparing images using the 

Hausdorff distance. IEEE Trans Pattern Anal Mach Intel, 15(9), p. 850–863.

HYCOM GoM experiment 31.0, n.d. HYCOM + NCODA Gulf of Mexico 1/25° Analysis 

Experiment 31.0. [Online] Available at: https://hycom.org/dataserver/gom-analysis

Incardona, J. P. et al., 2014. Deepwater Horizon crude oil impacts the developing hearts of large 

predatory pelagic fish. P Natl Acad Sci, 111(15), p. E1510–E1518.

Johansen, O., Brandvik, P. J. & Farooq, U., 2013. Droplet breakup in subsea oil releases - Part 2: 

Predictions of droplet size distributions with and without injection of chemical dispersants. 

Mar Pollut Bull, 73(1), p. 327–335.

Kalman, R. E., 1960. A new approach to linear filtering and prediction problems. Trans ASME J 

Basic Eng, 82(Series D), p. 35–45.
109

Ko, D. S., Preller, R. & Martin, P., 2003. An experimental real-time Intra-Americas sea ocean 

nowcast/forecast system for coastal ocean prediction. Seattle, WA, American Meterology 

Society.

Kondo, J., 1975. Air-sea bulk transfer coefficients in diabatic conditions. Bound-Layer Meteor. , 

Volume 9, p. 91–112.

Kourafalou, V. H. & Androulidakis, Y. S., 2013. Influence of Mississippi River induced 

circulation on the Deepwater Horizon oil spill transport. J Geophys Res, 118(8), p. 3823–

3842.

Kourafalou, V. H. & Androulidakis, Y. S., 2013. Kourafalou, V. H., & Androulidakis, Y. S. 

(2013). Influence of Mississippi River induced circulation on the Deepwater Horizon oil spill 

transport. Journal of Geophysical Research: Oceans, 118(8), pp. 3823-3842.

Kudela, R. M. & Palacios, S. L., 2007. Optical characterization and age estimates of river 

plumes on the US west coast. Bellingham, WA, Society of Photo-optical Instrumentation 

Engineers, pp. 66800L-1– 66800L-10.

Kudryavtsev, V., Shrira, V., Dulov, V. & Malinovsky, V., 2008. On the vertical structure of 

wind-driven sea currents. J Phys Oceanogr, 38(10), p. 2121–2144.

Lagerloef, G., 2012. Satellite Mission Monitors Ocean Surface Salinity. Eos, 93(25), p. 233–240.

Large, W. G., McWilliams, J. C. & Doney, S. C., 1994. Oceanic vertical mixing: a review and a 

model with a nonlocal boundary layer parameterization. Rev Geophys, 32(4), p. 363–403.

Le Hénaff, M., Kourafalou, V. H., Morel, Y. & Srinivasan, A., 2012. Simulating the dynamics 

and intensification of cyclonic Loop Current frontal eddies in the Gulf of Mexico. J. 

Geophys. Res, Volume 117, p. C02034.

Leben, R. R., 2005. Altimeter-derived Loop Current metrics. Circulation in the Gulf of Mexico: 

Observations and Models, pp. 181-201.

LeHenaff, M. et al., 2012. Surface evolution of the Deepwater Horizon oil spill patch: Combined 

effects of circulation and wind-induced drift. Environmental science & technology, 46(13), 

pp. 7267-7273.

Leifer, I. et al., 2012. State of the art satellite and airborne marine oil spill remote sensing: 

Application to the BP Deepwater Horizon oil spill. Remote Sens Environ, 120(0), pp. 185-

209.

Leung, T. & Malik, J., 2001. Representing and recognizing the visual appearance of materials 

using three-dimensional textons. Int J Comput Vision, 43(1), pp. 29-44.

Lindo-Atichati, D., Bringas, F. & Goni, G., 2013. Loop Current excursions and ring detachments 

during 1993–2009. Int. J. Remote Sens, 34(14), p. 5042–5053.
110

Lohrenz, S. E. et al., 1997. Variations in primary production of northern Gulf of Mexico 

continental shelf waters linked to nutrient inputs from the Mississippi River. Mar Ecol Prog 

Ser, Volume 155, p. 45–54.

Lonin, S. A., 1999. Lagrangian model for oil spill diffusion at sea. Spill Sci Technol B, 5(5/6), 

pp. 331-336.

Lu, Y. et al., 2016. Lu Y, Sun S, Zhang M, Murch B, Hu C. 2016. Refinement of the critical 

angle calculation for the contrast reversal of oil slicks under sunglint. Refinement of the 

critical angle calculation for the contrast reversal of oil slicks under sunglint. J Geophys Res 

Oceans. 121: 148–161, Volume 121, pp. 148-161.

Mabile, N. J., 2013. Considerations for the application of controlled in-situ burning. Oil Gas 

Facil., 2(2).

MacDonald, I. et al., 1993. Natural oil-slicks in the Gulf of Mexico visible from space. J 

Geophys Res Oceans, 98(C9), p. 16351–16364.

MacDonald, I. R., 2011. Remote sensing and sea-truth measurements of methane flux to the 

atmosphere (HYFLUX project) Final Report, Sugar Land, TX: National Energy Laboratory, 

Department of Energy.

MacDonald, I. R. et al., 2015a. Natural and unnatural oil slicks in the Gulf of Mexico. J Geophys 

Res Oceans, 120(12), p. 8364–8380.

MacDonald, I. R. et al., 2015b. Dataset: Neural network analysis determination of oil slick 

distribution and thickness from satellite Synthetic Aperture Radar, April 24 - August 3, 2010. 

[Online] 

Available at: https://data.gulfresearchinitiative.org/data/R1.x132.137:0045/

[Accessed 21 01 2016].

MacDonald, I. R. et al., 2002. Transfer of hydrocarbons from natural seeps to the water column 

and atmosphere. Geofluids, 2(2), p. 95–107..

McBride, B. C., Rowan, M. G. & Weimer, P., 1998. The evolution of Allochthonous salt 

systems, northern Green Canyon and Ewing Bank (offshore Louisiana), northern Gulf of 

Mexico. AAPG Bull, 82(5B), pp. 1013-1036.

McNutt, M. K. et al., 2011. Review of flow rate estimates of the Deepwater Horizon oil spill. P 

Natl Acad Sci, 109(50), p. 20260–20267.

McNutt, M. K. et al., 2012. Applications of science and engineering to quantify and control the 

Deepwater Horizon oil spill. P Natl Acad Sci, 109(50), p. 20222–20228.

McWilliams, J. C., 1990. The vortices of geostrophic turbulence. J Fluid Mech, Volume 219, p. 

387–404.
111

Meir, T. et al., 2013. Forecasting the New York City Urban Heat Island and Sea Breeze during 

Extreme Heat Events. Wea. Forecasting, Volume 28, p. 1460–1477.

Metzger, E. J. et al., 2014. US Navy operational global ocean and Arctic ice prediction systems. 

Oceanogr, 27(3), pp. 32-43.

Miller, R. L. & McKee, B. A., 2004. Using MODIS Terra 250 m imagery to map concentrations 

of total suspended matter in coastal water. Remote Sens Environ, Volume 93, p. 259–266.

Mitchell, R., MacDonald, I. R. & Kvenvolden, K., 1999. Estimates of total hydrocarbon seepage 

into the Gulf of Mexico based on satellite remote sensing images. EOS Supp, 80(49), p. 

OS242.

Monismith, F. & Fong, D., 2004. A note on the potential transport of scalars and organisms by 

surface waves. Limnol Oceanogr, Volume 49, p. 1214–1217.

National Geophysical Data Center/NESDIS/NOAA/U.S. Department of Commerce, 2001. 

ETOPO2, Global 2 Arc-minute Ocean Depth and Land Elevation from the US National 

Geophysical Data Center (NGDC). [Online] Available at: 

http://dx.doi.org/10.5065/D6668B75 [Accessed 06 06 2013].

National Oceanographic Data Center, 2013. Glider and Float Data: Deepwater Horizon Support, 

Float Data, 17 Aug 10-29 Nov 10. [Online] Available at: 

http://www.nodc.noaa.gov/deepwaterhorizon/glider_float.html. [Accessed 02 11 2016].

National Resource Council, 2003. Oil in the Sea: Inputs, Fates, and Effects. Washington, 

D.C.(Washington, D.C.): The National Academies Press.

Nikuradse, J., 1933. Stromungsgesetze in rauben Rohren. 361 ed. Verein Deut: VDI-Verlag 

Edn..

Nixon, Z. et al., 2016. Shoreline oiling from the Deepwater Horizon oil spill. Mar Pollut Bull, 

107(1), pp. 170-178.

Office of Response and Restoration, Emergency Response Division (ERD or Hazmat), 2012. 

Open Water Oil Identification Job Aid for Aerial Observation. Seattle, WA: s.n.

Okubo, A., 1970. Horizontal dispersion of floatable particles in the vicinity of velocity 

singularity such as convergences. Deep Sea Res, Volume 17, p. 445–454.

Olascoaga, M. J. & Haller, G., 2012. Forecasting sudden changes in environmental pollution 

patterns. P Natl Acad Sci, 109(13), p. 4738–4743.

Pasch, R., 2010. Tropical Cyclone Report: Hurricane Alex (AL012010) 25 June–2 July 2010, 

Miami, FL: NOAA/National Weather Service.

Passow, U., 2014. Formation of rapidly-sinking, oil-associated marine snow. Deep-Sea Res Pt II.
112

Pindell, J. L. & Kennan, L., 2007. Rift models and the salt-cored marginal wedge in the northern 

Gulf of Mexico: Implications for deep water Paleogene Wilcox deposition and basin-wide 

maturation. Houston, TX, Gulf Coast Section, Society for Sedimentary Ecology 

(GCSSEPM), p. 146–186.

Pullen, J. et al., 2011. Atmospheric and oceanic processes in the vicinity of an island strait. 

Oceanography, 24(1), p. 113–121.

Reanalysis.org, 2016. Atmospheric Reanalyses Comparison. [Online] Available at: 

https://reanalyses.org/atmosphere/comparison-table [Accessed 15 02 2016].

Reed, M., Turner, C. & Odulo, A., 1994. The role of wind and emulsification in modelling oil 

spill and surface drifter trajectories. Spill Sci Technol B, 1(2), p. 143–157.

Rucklidge, W. J., 1997. Efficient locating objects using Hausdorff distance. Int J Comput Vision, 

24(3), p. 251–270.

Sadarjoen, I. A. & Post, F. H., 2000. Detection, quantification, and tracking of vortices using 

streamline geometry. Comput Graph, Volume 24, p. 333–341.

Saha, S. et al., 2010. NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly TimeSeries Products, January 1979 to December 2010. [Online] Available at: 

http://dx.doi.org/10.5065/D6513W89 [Accessed 20 09 2013].

SAMOS, n.d. Shipboard Automated Meteorological and Oceanographic System. [Online] 

Available at: https://samos.coaps.fsu.edu/html/related.php [Accessed 08 2014].

Samuels, W. B., Huang, N. E. & Amstutz, D. E., 1982. An oilspill trajectory analysis model with 

a variable wind deflection angle. Ocean Eng, 9(4), p. 347–360.

Shadden, S. C., Lekien, F. & Marsden, J. E., 2005. Definition and properties of Lagrangian 

coherent structures. Physica D, Volume 212, p. 271–304.

Shaden, S. C., 2012. Lagrangian coherent structures. In: Transport and Mixing in Laminar 

Flows: From Microfluidics to Oceanic Currents. s.l.:Wiley-VCH Verlag GmbH & Co. 

KGaA..

Sharp, E., Dodds, P., Barrett, M. & Spataru, C., 2015. Evaluating the accuracy of CFSR 

reanalysis hourly wind speed forecasts for the UK, using in situ measurements and 

geographical information. Renewable Energy, Volume 77, p. 527–538.

Silva, M., Etnoyer, P. J., Demopoulos, A. W. & MacDonald, I. R., 2015. Coral injuries observed 

at mesophotic reef communities following the Deepwater Horizon oil discharge. Deep-Sea 

Res Pt II..

Smith, R. D. & McWilliams, J. C., 2003. Anisotropic horizontal viscosity for ocean models. 

Ocean Model, 5(2), p. 129–156.
113

Smith, S. D. et al., 1992. Sea surface wind stress and drag coefficients: the HEXOS results. 

Bound-Layer Meteorol, Volume 60, p. 109–142.

Smith, S. D. et al., 1992. Sea surface wind stress and drag coefficients: the HEXOS results. 

Bound-Layer Meteorol, Volume 60, p. 109–142.

Sobey, R. J. & Barker, C. H., 1997. Wave-driven transport of surface oil. J Coastal Res, 13(2), p. 

490–496.

Son, Y. B. & Gardner, W. D., 2011. Climatological variability of surface particulate organic 

carbon (POC) and physical processes based on ocean color data in the Gulf of Mexico. 

Korean J Remote Sens, 27(3), p. 235–258.

Stopa, J. E. & Cheung, K. F., 2014. Stopa and Cheung. 2014. Intercomparison of wind and wave 

data from the ECMWF Reanalysis Interim and the NCEP Climate Forecast System 

Reanalysis. Ocean modelling, Volume 75, p. 65–83.

Sun, S. & Hu, C., 2016. Sun glint requirement for the remote detection of surface oil films. 

Geophys Res Lett, Volume 43, p. 309–316.

Sun, S. et al., 2016. Oil slick morphology derived from AVIRIS measurements of the Deepwater 

Horizon oil spill: Implications for spatial resolution requirements of remote sensors. Mar 

Pollut Bull, 103(102), p. 276–285.

Sun, S., Hu, C. & Tunnell Jr, J. W., 2015. Surface oil footprint and trajectory of the Ixtoc-I oil 

spill determined from Landsat/MSS and CZCS observations. Mar Pollut Bull, Volume 101, 

p. 632–641.

Tolman, H., 2009. User manual and system documentation of WAVEWATCH IIITM, version 

3.14. Contribution 276.. Camp Springs(MD): Environmental Modeling Center, Marine 

Modeling and Analysis Branch for the National Centers of Environmental Prediction.

Vidal, V. M. et al., 1994. Winter mass distributions in the western Gulf of Mexico affected by a 

colliding anticyclonic ring. J. Oceanography, 50(5), p. 559–588.

Wallcraft, A. J. et al., 2009. Comparisons of monthly mean 10 m wind speeds from satellites and 

NWP products over the global ocean. J Geophys Res, Volume 114, p. D16109.

Wang, M. & Hu, C., 2015. Extracting oil slick features from VIIRS nighttime imagery using a 

Gaussian filter and morphological constraints. IEEE Geosci Remote Sens Lett, 12(10), p. 

2051–2055.

Watkins, J. S. et al., 1978. Occurrence and evolution of salt in deep Gulf of Mexico. In: A. H. 

Bouma, G. T. Moore & J. M. Coleman, eds. Framework, Facies, and Oil-Trapping 

Characteristics of the Upper Continental Margin. Tulsa, OK: AAPG Studies in Geology, pp. 

3-20.
114

Weber, T. C. et al., 2014. Acoustic estimates of methane gas flux from the seabed in a 6000 km2 

region in the northern Gulf of Mexico. Geochem Geophys Geosyst., 15(5), p. 1911–1925.

Weiss, J., 1991. The dynamics of enstrophy transfer in two dimensional hydrodynamics. Physica 

D, Volume 48, p. 273–294.

Wetzel, R. G., 2001. Limnology, Lake and River Ecosystems. Third ed. San Diego: Academic 

Press.

White, H. K. et al., 2012. Impact of the Deepwater Horizon oil spill on a deep-water coral 

community in the Gulf of Mexico. P Natl Acad Sci, 109(50), p. 20303–20308.

Zamudio, L. & Hogan, P. J., 2008. Nesting the Gulf of Mexico in Atlantic HYCOM: 

Oceanographic processes generated by Hurricane Ivan. Ocean Model, Volume 21, p. 106–

125.

Zhang, D. & Lu, C., 2004. Review of shape representation and description techniques. Pattern 

Recognition, 37(1), pp. 1-19.
115

APPENDIX A.

Kalman Filter Algorithm for Identification of the Loop Current Front

In the following analysis, the assumption is made that the LC front closely follows the core 

of the LC. Under a geostrophic assumption, the core coincides with the maximum SSH gradient. 

Thus, the true LC front is determined by the maximum SSH gradient. Rationale of the suggested 

alternative approach for LC tracking is based on the idea of combining the knowledge about 

frontal position from two dynamic fields by employing the Kalman Filter technique (Kalman 

1960; Maybeck 1979). In the described application, GoM-HYCOM SSH field provides the first 

guess of the LC front location (a priori estimate). GoM-HYCOM SSH gradient field is used to 

obtain the final frontal location (a posteriori estimate). Any other oceanographic field capable of 

capturing mesoscale structures can be used as a “measurement.”

In following the LC front from the Yucatan Channel to the Straits of Florida, the task is to 

estimate the next position along the contour

 

T

k

kkx ,

( and  are coordinates) given the 

previous position xk-1. It is assumed that the process can be described by the linear stochastic 

difference equation

xk Axk1wk1

, (B1)

where xk is the kth discrete position along the LC contour moving from Yucatan Channel to 

Florida Straits, xk-1 is the previous position, and wk-1 is the process noise. The matrix A relates 

the state at the previous step (k-1) to the state at the current step k. The matrix A changes with 

each step. The prediction is corrected on the basis of a measurement, whose state is

zk Hx k νk

. (B2)

The matrix H relates the current state to the measurement zk. In this application, H does not 

change. k represents the measurement noise. 

Both random variables are assumed to be independent and normally distributed,

w ~ N0,Q

and 

ν ~ N0,R, where Q and R are process noise and measurement noise 

covariance matrices, respectively. They change at every step. Smaller Q or R suggests more 

confidence in predictions or measurements. 

At every step, the a priori state estimate at step k is obtained on the basis of knowledge of the 

process prior to step k. Then the a priori state is corrected on the basis of measurements, resulting 

in the a posteriori state estimate. Also defined are the a priori and a posteriori estimate errors

 

k

 k ke x x

ˆ

(B3)

k k k

e  x  x

ˆ . (B4)

Then a priori and a posteriori error covariance matrices are

 

T

k

E k k

  

P  e e , (B5)

 

T

k

E k kP  e e . (B6)
116

An a posteriori state estimate 

xk

ˆ

is defined as a linear combination of an a priori estimate 



k x

ˆ

and a weighted difference between a measurement zk and predicted measurement 



k Hx

ˆ

(also 

known as innovation or residual)

 

 

k

 k k kx

ˆ

x

ˆ K z Hxˆ . (B7)

The matrix K is the gain (blending) factor that minimizes the a posteriori error covariance (B6)

 

1

 

K  P H HP H R

T

k

T

k k

. (B8)

At every step, the LC tracking algorithm assesses the reliability of the a priori estimate and 

measurement depending on predefined criteria and changes Q and R accordingly. 

The algorithm for identification of the LC using the Kalman Filter technique is as 

follows.

1) Obtain the a priori state estimate. The a priori state estimate (red dot in Figure A1) is 

obtained by following in the direction perpendicular to the local SSH gradient from

1

ˆ

xk

(red 

arrow in Figure A1). Now matrix A, which is required for projecting the error covariance matrix 

to the next step and is unknown at k

th location, is estimated as

 

1

1

ˆ ˆ







A xkxk

. (B9)

2) Update Q. The process noise covariance Q and the measurement error covariance R are 

updated following a set of criteria to judge the trustworthiness of



k x

ˆ

and zk. The SSH at the a 

priori estimate of the frontal location (k=( )) and SSH averaged over the previous five 

locations (



) are used to update the process noise covariance Q,

 

   







  

  



   

   

k k

k k

if

if C

1000 ,

1000 ,

1.2

1.2

(B10)

,

where agrees with the uncorrelated process noise assumption.



x ˆ k_

Q Q0 C



















0 0.01

0.01 0

Q0
117

Figure A1. LC identification algorithm based on the Kalman Filtering

The inset is GoM-HYCOM SSH gradient (cm/km) field in the GoM. The red box is the region where the 

LC tracking technique is illustrated. In the zoomed region, the background field (gray shading) is the SSH 

gradient (cm/km). The black and red lines are the 0.05-m and 0.17-m SSH contours, respectively. The 

blue lines are contours of the SSH gradient drawn at 0.1 cm/km intervals. The red dot is the a priori state 

estimates, the green dot is measurement, and the purple circles are the a posteriori state estimates of 

frontal locations at discrete positions. See text for more detail.

3) Project the error covariance from state k-1 to k

P  AP A  Q

 T

k k 1

(B11)

4) Obtain measurement location. The measurement of the frontal location Zk (green dot Figure

A1) is given by the local maximum of the SSH gradient (||max). The local maximum is 

searched in the direction (green arrow in Figure A1) determined by the local second derivative of 

the SSH at the a priori location (



k x

ˆ

).

5) Update R. Local SSH is obtained for measurement location Zk (z = (Zk)) and R is updated 

as

k z

c0100 

(B12)

  

1.2 4

0

exp 10  DAO C c

(B13)
118 R  C R0

, (B14)

where DAO is distance (meters) between 



k x

ˆ

and Zk, and 



















0 0.01

0.01 0

R0

.

6) Compute the Kalman gain. Equation (B8) is used. The matrix H is 



















0 1

1 0

H . (B15)

This design of H means that in the absence of process error and measurement error, the process 

is identical to the measurement. 

7) Update the a priori estimate using (B7) to obtain the a posteriori estimate.

8) Update the error covariance 

 

Pk Pk KkHPk . (B16)

Steps 1–8 are then repeated to obtain the next discrete location along the LC front. 
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 the 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 communities.

The Bureau of Ocean Energy Management Mission

The Bureau of Ocean Energy Management (BOEM) works to manage the 

exploration and development of the nation's offshore resources in a way that 

appropriately balances economic development, energy independence, and 

environmental protection through oil and gas leases, renewable energy 

development and environmental reviews and studies.