Research ArticleECOLOGY

Invisible oil beyond the Deepwater Horizon satellite footprint

See allHide authors and affiliations

Science Advances  12 Feb 2020:
Vol. 6, no. 7, eaaw8863
DOI: 10.1126/sciadv.aaw8863


Major oil spills are catastrophic events that immensely affect the environment and society, yet determining their spatial extent is a highly complex task. During the Deepwater Horizon (DWH) blowout, ~149,000 km2 of the Gulf of Mexico (GoM) was covered by oil slicks and vast areas of the Gulf were closed for fishing. Yet, the satellite footprint does not necessarily capture the entire oil spill extent. Here, we use in situ observations and oil spill transport modeling to examine the full extent of the DWH spill, focusing on toxic-to-biota (i.e., marine organisms) oil concentration ranges. We demonstrate that large areas of the GoM were exposed to invisible and toxic oil that extended beyond the boundaries of the satellite footprint and the fishery closures. With a global increase in petroleum production–related activities, a careful assessment of oil spills’ full extent is necessary to maximize environmental and public safety.


The Deepwater Horizon (DWH) blowout was a mega oil spill, introducing ~795 million liters of live oil into the Gulf of Mexico (GoM) with oil slicks covering an estimated area of 149,000 km2 (1). As a result, vast areas of the GoM were closed for fishing, at one point reflecting a maximum of more than a third of the U.S. national exclusive economic zone (2). The application of the fishery closures was based on satellite and areal imagery by the National Environmental Satellite, Data, and Information Service (NESDIS), using mainly synthetic aperture radar (SAR), and tracking of the NESDIS footprints forward in time, predicting their approximated locations for up to 72 hours (3). The closures were reopened on the basis of a systematic seafood sampling, performing chemical and sensory testing in various seafood specimens (2).

The cumulative satellite oil slick footprint was largely accepted as the DWH oil spill extent from scientific, public, and management perspectives (2, 4, 5). Yet, accumulating field data support a much wider extent of the DWH spill beyond the satellite footprint, reaching the West Florida Shelf (WFS), the Texas Shores (TXS), the Loop Current (LC) system, and the Florida Keys (FK) (see regional map in fig. S1). Specifically, in the WFS, studies indicate high concentrations of oil, including toxic and mutagenic levels, in the water (6), in sediments (7), in sand patties (8), and on the coast (9). Furthermore, high levels of polycyclic aromatic hydrocarbons (PAHs) from DWH oil were found in red snappers’ livers, which co-occurred with a high frequency of skin lesions in bottom-dwelling fish (10). Last, satellite imagery and particle tracking showed that an oil slick present east of Pensacola (north WFS) was transported southeast along the WFS, reaching Tampa and the Dry Tortugas (FK) within a few weeks (11). In TXS, high and toxic levels of oil were found in the water (12) and in sediments (12, 13). In the LC system, the European Space Agency reported the presence of DWH oil during mid-May (14). Later, during early June, NESDIS satellite imagery revealed the presence of 12 oil slicks in the LC system, stretching from FK to the GoM interior. High concentrations of oil were reported in the LC between late June and mid-July west of the fishery closures (15). Furthermore, a deep intrusion (~1000 to 1300 m), which was documented (16) and represented in oil transport simulations (17), included toxic PAH concentrations within 13 km from the well and above-background concentrations extending southwest ~300 km beyond the satellite footprint and closure areas. Overall, these observations indicate that DWH oil extended beyond the satellite footprint and the fishery closures (18).

Oil transport models are effective tools for predicting and reconstructing marine oil spills. Numerous modeling efforts were previously carried out in an attempt to reconstruct spatiotemporal dynamics of the DWH oil spill. The surface oil movement was modeled using Lagrangian Coherent Structure (LCS) core analysis, successfully reconstructing the oil footprint from the satellite imagery (19). Similarly, the three-dimensional (3D) oil application of the Connectivity Modeling System (oil-CMS) (20) and the General NOAA (National Oceanic and Atmospheric Administration) Operating Modeling Environment (GNOME) (21) simulated the oiled coastline areas, both successfully reconstructing the observed beaching patterns (9). The GNOME was also used to conduct probabilistic Monte Carlo simulations, running hundreds of oil transport scenarios based on historical wind and currents data (22). While some of these models provided possible scenarios of the DWH spill, a robust spatiotemporal examination of the spill’s extent, resolving toxic and nontoxic concentrations, does not exist for the DWH. Similarly, the relationship between satellite detection and in situ oil concentrations was rarely examined. This is a key point because satellite detection is the main source of information used for determining the oil spill extent and the corresponding management actions, e.g., fishery closures (2). Having the capacity to quantify the toxic-to-biota extent, which is not necessarily equal to the satellite-detected extent, would improve the understanding and, therefore, the management of oil spills.

This is particularly important in light of photoinduced toxicity studies that have demonstrated that the combined effect of PAH and ultraviolet (UV) radiation can be two orders of magnitude more potent than PAH alone (4). This marked increase in hazard to early life stage was never considered in oil transport simulations thus far, and this is the first study to do so. The aim of this study was to examine the full extent of the DWH oil spill based on satellite/areal imagery, in situ measurements, and oil transport model, focusing on the oil fraction that is invisible to satellites and toxic to biota. To compute this fraction, we quantify the relationship between PAHs, which are associated with oil toxicity, and total petroleum hydrocarbons (TPHs), which closely represent the oil modeled in oil transport simulations (23). We then quantify the satellite detection threshold in terms of oil concentrations. Last, we apply these elements in an oil transport model, compute the total and toxic-to-biota extents of the DWH oil spill, and discuss their implications on the environment and public health.


Toxicity and visibility thresholds from in situ and satellite/areal observations

Linear regression indicated a significant relationship (P < 0.001, R2 = 0.9, n = 34) between TPH and PAH from the BP Gulf Science Data (GSD) for surface waters such that (Fig. 1, A and B)log10(TPH+1)=1.733+1.007×log10(PAH+1)(1)

Fig. 1 Toxicity and visibility thresholds from in situ and areal/satellite observations.

(A and B) The relationship between TPH and PAH concentrations from the GSD for surface waters (depth, 0 to 1 m; linear regression P < 0.001, R2 = 0.9, n = 34; black line). The blue line represents the regression of the same dataset after removing the five influential points with the highest concentrations (linear regression P < 0.01, R2 = 0.2, n = 29). Green points represent background and naturally slick-dispersed water samples from (24) (linear regression P < 0.01, R2 = 0.21, n = 33, green line). (C) The relationship between TPH and PAH concentrations from the GSD for deeper waters (depth, >1 m; linear regression P < 0.001, R2 = 0.3, n = 641; gray line). (D) The area of the NESDIS anomaly footprint across time (days from blowout). The green line represents the NESDIS anomaly decay (linear regression P < 0.001, R2 = 0.63, n = 22) from day 87 when the spill was stopped (blue vertical line) and onward, until NESDIS signal disappeared (orange vertical line). (E) PAH concentrations at the surface (depth, 0 to 1 m) from the GSD. In (A) to (C) and (E), the red lines represent the toxic-to-biota threshold of PAH = 0.5 ppb for the surface (depth, 0 to 1 m), and PAH = 1 ppb for the water column (depth, >1 m). Dashed magenta lines in (A, B, and E) represent the NESDIS detection threshold in terms of PAH and TPH concentrations. (F) Toxic (red/dark), above-background (green/bright), and below-background (blue/smaller markers) GSD in situ PAH water samples across depth and time. PAH and TPH data were log10(x + 1)-transformed in (A) to (C), (E), and (F).

The relationship remained significant (P < 0.01, R2 = 0.2, n = 29) with similar coefficients when removing the five influential high-concentration samples, such that (Fig. 1, A and B)log10(TPH+1)=1.514+1.021×log10(PAH+1)(2)

Background and naturally dispersed water samples collected near slick-dispersed areas in a separate study (24) were added to Fig. 1 (A and B) (green points), demonstrating a significant linear relationship (Fig. 1, A and B) (P < 0.01, R2 = 0.21, n = 33) in agreement with the GSD regression lines, such that (Fig. 1, A and B)log10(TPH+1)=1.7+0.91×log10(PAH+1)(3)

Linear regression for deeper waters (depth, >1 m) for TPH and PAH concentrations from GSD indicated a significant linear relationship (P < 0.001, R2 = 0.3, n = 641), such that (Fig. 1C)log10(TPH+1)=1.584+0.853×log10(PAH+1)(4)

Regression Eqs. 1 to 4 demonstrate a significant and robust relationship between TPH and PAH for surface and deep waters.

Overall, these regression equations indicate that TPH consist of ~1.5 to 3% of PAHs and falls within the range of previous studies, e.g., 1.5% (12) and 0.2 to 7% (25).

After the spill was contained (15 July 2010), the decay of the NESDIS anomaly signal was rapid (~502 km2 day−1), disappearing 22 days later (Fig. 1D; linear regression, area = 11,054.1 − 501.6 × days, R2 = 0.63, P < 0.001, n = 22). In contrast, the dissolved oil concentrations in the water showed a much slower decay, with toxic-to-biota concentrations persisting 68 days later (Fig. 1E). From the time the NESDIS signal disappeared, the mean (±SD) of the five maximal sampled concentrations of PAH and TPH were 17.02 ± 8.9 and 964 ± 392 parts per billion (ppb), respectively, representing the NESDIS anomaly sensitivity threshold in terms of PAH and TPH concentrations (Fig. 1, A and B). The range between the toxicity (PAH = 0.5 ppb) (4) and the NESDIS sensitivity (PAH = 17 ppb) thresholds represents toxic to biota oil concentrations, which are invisible to satellite/areal detection methods (Fig. 1, A, B, and E). Last, in deeper waters (depth, >1 m), toxic-to-biota levels of PAH (PAH > 1 ppb) were present across the water column, especially in the deep plume at a depth of 1000 to 1300 m (Fig. 1F) (16), and were beyond satellite/areal detection capabilities.

Spatial DWH cumulative extents

The spatial extent of the NESDIS anomaly footprint was limited to the north central area of the GoM and was nearly fully contained by the fishery closures (Fig. 2A). In contrast, the simulated cumulative oil spill extent, integrated across time and depth, spread well beyond the areas denoted by the NESDIS anomaly footprint and the fishery closures (Fig. 2, A and B). Specifically, the oil extended beyond the NESDIS anomaly footprint to TXS, WFS, LC, FK, and along the gulf stream toward the East Florida Shelf (EFS) (for visualization, see movie S1; cumulative mass across time and space is in fig. S2). This spatial distribution is in agreement with the higher-than-background PAH samples (Fig. 2B and fig. S3), as well as with the spatially explicit information reported in previous studies for Apalachee Bay (AB) (9), WFS (7), TXS (12, 13), LC (15), and the Deep Plume (DP) (16) (Fig. 2B). NESDIS anomaly also confirmed the presence of multiple oil patches in the LC and the FK (Fig. 2A). Notably, the spatial shape and location of the NESDIS anomaly (Fig. 2A) correspond to the regions containing the highest simulated cumulative oil concentrations (Fig. 2B) and to the area computed as visible and toxic to biota (Fig. 2C). Yet, large areas with oil concentrations below ~1000 ppb were not captured by the NESDIS anomaly and by the fishery closures (Fig. 2C) but were in the toxic ranges for the surface (depth, 0 to 1 m; PAH > 0.5 ppb) and the water column (depth, >1 m; PAH > 1 ppb). The duration of the presence of toxic-to-biota concentrations (Fig. 2D) corresponds to areas that exhibited high oil concentrations (Fig. 2, B and C), with a maximum of 84 days of toxic-to-biota concentrations in a small area west of the wellhead. Photoinduced toxicity studies demonstrate that the PAH concentrations that cause 50% mortality (LC50) or other deleterious effects (EC50) to the test organisms are largely in the range of invisible and toxic oil. Oil at these concentrations for surface water (depth, 0 to 1 m) extended beyond the satellite footprint and the fishery closures, potentially exterminating a vast amount of planktonic marine organisms across the domain (Fig. 2F).

Fig. 2 Spatial DWH cumulative extents.

(A) Cumulative NESDIS anomaly daily composites integrated from 20 April 2010 to 21 July 2010. Daily fishing closures are marked with gray lines; the cumulative fishing closure area is marked with a thick dashed yellow line. The black star represents the location of the DWH blowout [adapted by permission from Springer-Nature: Scenarios and Responses to Future Deep Oil Spills: Fighting the Next War by S. Murawski, C. Ainsworth, S. Gilbert, D. Hollander, C. Paris, M. Schlüter, D. Wetzel, Eds., 2019 (18)]. (B) Cumulative value of daily average oil concentrations (ppb), integrated across the same time span as (A) and across water depths. Vertical depth layers are 0 to 1 m, 1 to 20 m, and in 20-m increments down to 2500 m. Sediment and water samples with higher-than-background concentration are marked in bright green and dark blue circles, respectively [adapted by permission from Springer-Nature: Scenarios and Responses to Future Deep Oil Spills: Fighting the Next War by S. Murawski, C. Ainsworth, S. Gilbert, D. Hollander, C. Paris, M. Schlüter, D. Wetzel, Eds., 2019 (18)]. Red crosses in (B) represent approximate locations of DWH-related oil detections reported in previous studies: (i) (9), (ii) (7), (iii) (12), (iv) (13), and (v) (15). Daily fishery closures are marked with black polygons; the cumulative fishery closure area is marked with a dashed thick polygon. AB, Apalachee Bay; DP, Deep Plume; EFS, East Florida Shelf; FK, Florida Keys; LC, Loop Current System; TXS, Texas Shores; WFS, West Florida Shelf. (C) Categorization of the modeled oil spill are as follows: (i) nontoxic, PAH concentrations above background level and smaller than 0.5 and 1 ppb at the surface (depth, 0 to 1 m) and in the water column (depth, >1 m), respectively; (ii) toxic-to-biota and invisible, PAH concentrations 0.5 to 17 ppb at the surface and above 1 ppb in the water column; and (iii) toxic and visible, PAH concentrations above 17 ppb. In (C), categories were computed according to maximal concentrations across time. (D) Duration of toxic concentrations across the domain. (E) LC50 of 12 experiments examining the photoinduced toxicity to blue crab (31), fiddler crab (33), mahi mahi (29, 30), red drum (32), and speckled sea trout (32) (for more details, see table S2). (F) The spatial extent of the toxic concentrations from (E); color codes in (F) are according to bar colors in (E), representing concentrations exceeding LC50. In (F), toxic PAH of 0.5 ppb was concentrations were considered for surface waters only (depth, 0 to 1 m).

Spatial dynamics of the spill across time

The oil-CMS model indicated that the fishery closures captured 82 ± 19% (SD) of the total oil mass in the domain (Fig. 3A), 54 ± 20% (SD) of the total area covered by oil (Fig. 3B), and 70 ± 20% (SD) of the area with toxic-to-biota oil concentrations (Fig. 3C). In contrast, the fishery closures captured 94 ± 9% (SD) of the NESDIS anomaly (Fig. 3D). The closures match the NESDIS anomaly well, but an important component of the toxic oil is missed by the NESDIS anomaly (Fig. 2C) and by the fishery closures (Fig. 3C). The areas where the closures failed to capture the oil spill extent were primarily the TXS, LC, and WFS.

Fig. 3 Spatiotemporal dynamics of the spill.

(A) Mass, (B) area of total oil spill, (C) area of oil with toxic concentrations, and (D) area of NESDIS anomaly footprint, captured by the fishery closures (magenta) with respect to the total (black) present in the domain across 93 days from blowout. The NESDIS anomaly footprint and oil-CMS oil concentrations partitioned to visible and toxic to biota, invisible and toxic to biota, and invisible and nontoxic for 15 May 2010 (E to G), 6 to 18 June 2010 (H to J), and 2 July 2010 (K to M). In (F), (I), and (L), the color bar represents cumulative oil concentrations across time and depth. Sediment and water samples from the GSD in (F), (I), and (L) are marked in light green and dark blue circles, respectively, with light red outlines representing samples with higher-than-background concentrations. Toxicity was considered for the surface (depth, 0 to 1 m; PAH > 0.5 ppb) and for the water column (depth, >1 m; PAH > 1 ppb). (A), (B), (D), (F), (I), and (L) are adapted by permission from Springer-Nature: Scenarios and Responses to Future Deep Oil Spills: Fighting the Next War by S. Murawski, C. Ainsworth, S. Gilbert, D. Hollander, C. Paris, M. Schlüter, D. Wetzel, Eds., 2019 (18).

Looking at a few specific dates as examples, the oil-CMS indicated expansions of the spill extent on 15 May toward the LC (Fig. 3E), which is in agreement with the European Space Agency report for that same date (14). Here, too, the shape and location of the NESDIS anomaly (Fig. 3E) agree with the areas with high oil concentrations computed by the oil-CMS (Fig. 3F) and with the visible category of the oil-CMS (Fig. 3G). As the spill evolved, the cumulative NESDIS anomaly expanded (Fig. 3H), but to a much smaller extent compared with the oil distribution computed by the oil-CMS (Fig. 3I), which indicated an expansion to the WFS, TXS, and LC regions, which agrees with the GSD in situ water and sediment samples (Fig. 3I). Here, too, the shape and location of the NESDIS anomaly (Fig. 3H) agree with the areas with high oil concentrations computed by the oil-CMS (Fig. 3I) and with the visible category of the oil-CMS (Fig. 3J). According to the oil-CMS, some of the DWH oil flowed northbound with the gulf stream (Fig. 3I) toward the EFS, which is also confirmed by multiple higher-than-background water and sediment field samples for that region (fig. S3). Between 28 June and 1 July 2010, Hurricane Alex swept through the GoM, with strong southeasterly winds and waves enhancing the mixing and beaching of the DWH oil, causing a steep decline of more than 50% in the NESDIS anomaly coverage area (Fig. 3D). Here, too, the shape of the NESDIS anomaly (Fig. 3K) aligned with the areas of high oil concentrations (Fig. 3L) and with the visible category (Fig. 3M) computed by the oil-CMS.


The spatial synthesis of previously published studies (716), combined with the GSD water and sediment samples (26, 27), clearly shows that DWH oil spill extended beyond the satellite detection and the fishery closure’s limits, reaching TXS, WFS, DP, and LC (Fig. 2B) (18). Our results from the oil-CMS agree with that extent, providing the oil transport context for these observations and resolving an apparent spatial discrepancy between satellite footprint and field studies. When taking into account modes of toxicity such as photoinduced toxicity and cardiac malperformance (4, 2833), a substantial portion of the spill is “toxic and invisible.” The significant relationship between PAH and TPH (Fig. 1, A to C), which is in agreement with previous studies (12, 25), allowed us to express the PAH toxicity threshold in TPH and to apply that to the model output. Subsequently, we quantified the degree to which the applied fishery closures captured different fractions of the DWH spill and show that “toxic-to-biota and invisible” oil extended beyond the satellite footprint and the fishery closures (Figs. 2 and 3).

Previous spatial modeling efforts of DWH provided similar spatial extents to that of the oil-CMS. NOAA’s GNOME produced a nearly identical spatial probability heat map to that of the oil-CMS (Fig. 2B) (22). Similarly, other modeling studies using lagrangian coherent structures (LCS) core analysis and particle tracking (34, 35) indicated that DWH oil reached Texas Shores (TXS), WFS, LC, and FK. To the best of our knowledge, there is no published case of cumulative DWH spatial extent from oil transport models, contradicting the extent presented here. To account for a debated topic among oil spill modelers—the droplet size distribution (DSD)—we performed simulations using three types of debated droplet size ranges: 1 to 500, 80 to 2210, and 1 to 2400 μm. The oil spill spatial extents from the three DSDs were similar, with larger DSDs producing slightly higher concentrations at the surface compared with the lower DSD range (see section S1 and fig. S4).

Comparing oil transport modeling with in situ observations is challenging, especially since the GoM is highly active in terms of oil drillings and natural oil seeps (1). Hence, it is possible that some in situ samples contained non-DWH oil. Yet, few methods aid in determining the pollution source. Explicit forensic analyses were applied in samples of sediments (13), sand patties (8), and animal tissues (10), demonstrating that DWH oil was the source of contamination. For water samples, because of quick degradation of many compounds, a forensic fingerprinting analysis was not possible, yet spatial kriging interpolation methods (36) demonstrated a clear concentration gradient from the Macondo well, further supporting the DWH as the pollution source. Moreover, because of the patchy nature of the oil distribution in both time and space, many observations from multiple studies found no evidence of oil. Even within a radius of 10 km of the Macondo well during the oil spill, ~52% of the water samples were below PAH background level (26). This is important when weighing “no indication of oil” versus “positive indication of oil.” In other words, when studies report that they found no evidence of contamination, it does not necessarily mean that oil was not present in the wider area. Therefore, reports of absence of oil, for example, in the LC (37) during July 2010, do not necessarily indicate that the area was free of oil, particularly since the methods used in that study were not sensitive to low oil concentrations.

Similarly, oil slicks detected by satellites positively indicate the presence of toxic oil concentrations. Yet, when these slicks disappear, toxic concentrations may still be present (Fig. 2C) (38) and may persist for several days and weeks after the disappearance of satellite/areal signal (Fig. 1E) (38). For example, an oil slick identified at the north part of WFS, east of Pensacola, which was tracked forward in time, showed that even after 11 days, plume core concentrations remained high (~150 ppb) yet undetectable by satellites (11), corresponding to toxic PAH concentrations of ~1.8 ppb (Fig. 1, A and B). Similarly, a small in situ oil spill experiment demonstrated high toxicity of dissolved and dispersed oil beyond the edges of an oil slick and after the slick was no longer visible (38).

Deep-sea oil spills are dynamic and complex, and their hindcasts are always limited due to simplifying assumptions, stochasticity, and uncertainty in parameterization. For example, the extent of the spill detected by remote sensing is dependent on oil slick thickness, which is affected by wave energy, turbulence, and oil viscosity (39). As a consequence, turbulent regions such as the loop current are likely to create more mixing and dispersion of the oil slicks, which will limit the satellite detection ability for these regions. Only fragmented slicks were identified by satellites in the LC region, compared with larger and more consistent areas predicted by the oil transport model (Fig. 2, A and C). Other sources of uncertainty in the model include the currents output from the GoM–HYCOM (Hybrid Coordinate Ocean Model) with a spatial resolution of 0.04° and a daily temporal resolution. Finer resolution of currents in time and space will improve the accuracy of the model and reduce the associated uncertainties. Similarly, eddy diffusivity coefficient, the decay rates of the pseudocomponents, and the evaporation rates may vary in time and space. Implementation of realistic ranges of these parameters from opining studies will provide important information regarding the uncertainty of the model’s output and could further improve the agreement between the satellite footprint and the visible portion of the oil predicted by the model (Fig. 2, A and C).

It is important to state that satellite detection is the immediate, primary, and most reliable source of information regarding the presence of high oil concentrations at the surface, particularly for purposes of oil spill response and clean up. The oil transport and fate models provide an important complementary spatiotemporal information about the 3D distribution of the full concentration range of the oil spill, including toxic and invisible oil. Fishery closures could be determined using the computed toxic concentrations, applying any toxicity threshold set by management. In addition, this information from oil transport and fate models is relevant for impact assessment and for natural resource management when considering risks of potential oil spills.

Until recently, the estimated satellite detection threshold was roughly equal to the estimated level of concern (LOC). For example, in (22), predicted probabilities of regions receiving oil concentrations greater than LOC were computed using the GNOME (22), based on slick thickness of ~1 μm, which corresponds to the visible state of “dull sheen” of the oil (22). The lower range of NESDIS-detectable oil slick thickness for the DWH corresponds to approximately TPH = ~1000 ppb (Fig. 1), which is in agreement with the upper concentration ranges for slick-dispersed oil (24). However, recent studies revealed that toxicity to early-life-stage organisms may occur at much lower concentrations because of the combined effect of PAH and UV, which increases the toxicity by 10- to 200-fold compared with PAH alone (4, 28), as well as due to cardiac toxicity and developmental abnormalities at similarly low PAH levels (4, 40).

Nearly all organisms tested for photoinduced toxicity in recent studies exhibited LC50 at the ranges of invisible and toxic oil (Fig. 2E). For example, the LC20 of speckled sea trout larvae was PAH =12 ppb across 72 hours of exposure without UV, compared with ~0.15 ppb, when exposed to PAH and UV for a total test duration of only 24 hours (4). Note that photochemical reciprocity means that UV exposure concentrations are required for a more complete comparison, which is outside the scope of this study [see (31) for phototoxic dose metrics]. Considering this and other toxic outcomes at low concentrations (4), a much lower toxicity threshold of PAH = 0.5 ppb for surface waters was concluded in the federal DWH environmental impact assessment (4) and applied in the current study, showing that toxic concentrations extend beyond the satellite footprint (Fig. 2, C, E, and F).

Photoinduced PAH toxicity has the potential to present a hazard to a wide spectrum of taxa including fish (2830), invertebrates (31, 33, 41), and plants (42), which play key roles in the ecosystem, forming the base of the food web and primary productivity and consisting the base for the adult organisms pool. Damaging these important components at large scales can have deleterious effects on the ecosystem (43), especially considering that photoinduced toxicity may similarly affect many other taxa, which were not tested. Moreover, beyond the direct lethal effects of photoinduced toxicity, sublethal physiological effects such as reduced growth and impaired reproduction and health have also been observed in the field and in the laboratory (4). It is therefore important to map these concentrations in the domain to better understand which regions were potentially affected and to what possible extent. Since near-surface organisms are mainly affected by photoinduced toxicity, the deleterious effect can be expected throughout the pelagic realm both in LC and in TXS (Figs. 2 and 3). This information can be useful for impact assessment purposes and for restoration efforts (4).

The importance of our findings is not only limited to environmental effects but also could have implications for human health. While PAH toxicity for various marine species can be directly obtained from dose-response toxicity experiments, this direct approach cannot be applied for humans. In terms of seafood consumption, the DWH seafood safety protocol (2) considered the LOC of only 13 PAH compounds and did not consider the possible toxicity of photomodified compounds, i.e., the formation of new chemical compounds produced as a result of UV and PAH interaction, or PAH degradation products. Many of these compounds have been shown to be mutagenic and carcinogenic (4446), and their persistence across trophic levels is not well studied. In addition, the sample size of the DWH food sampling was simply too small to produce meaningful inference regarding the absence of contaminated seafood (section S2). Hence, the possible hazard associated with consumption of seafood containing oil and oil breakdown products, including photomodified compounds, ought to be carefully evaluated in the future.

In conclusion, our results change established perceptions about the consequences of the oil spill by showing that toxic and invisible oil extend beyond the satellite footprint at concentrations that present potentially lethal and sublethal hazards to a wide range of taxa under a range of circumstances found in the GoM. We did so by computing toxicity and visibility thresholds and applying these in an oil transport model. Our results are strongly supported by in situ and remote sensing observations. In future oil spills, toxic and invisible portions could be computed on the basis of the quantitative framework presented here. Similarly, our computed thresholds could be used as null-hypotheses benchmarks until more accurate data are obtained for a given scenario. We do not provide an exact oil concentration benchmark for fishery closures since that involves a trade-off analysis between ecosystem health risks, human health risks, and fishery revenue loss, which is beyond the scope of our current work. Instead, we quantify and highlight the concept of toxic-to-biota and invisible oil that can promote informed decision-making for future oil spills for managers, fishers, and seafood consumers. This is especially important given the global increase in deep-sea drilling efforts and petroleum-related activities.


Experimental design

To examine the extent and the spatiotemporal dynamics of the DWH oil spill, we compared data from three sources: (i) in situ water and sediment measurements (26, 27, 47) and other published studies (716), (ii) satellite/areal imagery footprint, and (iii) oil transport model (17). We quantified the relationship between PAH and TPH, since oil toxicity is mostly associated with PAH concentrations and TPH accounts for nearly all of the modeled oil (23). While previous studies indicate that oil contains 0.2 to 7% PAHs (25), it is important to quantify this relationship for the specific case of the Macondo oil similarly to the practice in (48) and (12). To quantify the satellite/areal imagery detection threshold of the oil, we examined the maximal PAH and TPH concentration values after the spill was no longer detected by satellite/areal imagery. We then applied these elements in the oil-CMS model to compute the oil spill extent including its toxic-to-biota and invisible fractions. To quantify the effectiveness of the fishery closures in capturing the different portions of the oil, we computed the daily oil mass, total area, and area with toxic-to-biota concentrations, which were captured by the fishery closures.

Quantitative and statistical analyses

Toxicity of oil is attributed to PAH concentrations, with thresholds of PAH > 0.5 ppb for the water surface (depth, 0 to 1 m) and PAH > 1 ppb for water column (depth, >1 m) (4). To quantitatively link between PAH and TPH, we used linear regression similarly to the practice in (15) and (12), expanding the regression analysis presented in (49). We further examined the robustness of the linear relationship for surface waters by computing linear regression after removing the five most influential data points, which had orders of magnitude higher concentrations compared with the rest of the dataset. Last, we used a different TPH-PAH dataset digitized from (24) to examine the robustness of our computed linear relationship.

We also used linear regression to quantify the decay rate in the NESDIS anomaly signal from the time the well was sealed (15 July 2010). Exponential decay regression was computed as well, but it did not provide a better fit compared with the linear model, possibly because oil droplets continued to rise to the surface even after the well was sealed (50). To determine the satellite detection threshold, we computed the mean of the maximal five PAH and TPH GSD samples since the last detection of DWH oil by NESDIS (12 August 2010). The obtained satellite detection threshold together with biota-toxic thresholds (4, 28) was used to define the toxic-to-biota and invisible ranges in the oil-CMS oil transport model. R software version 3.03 (51) was used for all data and statistical analyses.

In situ water and sediment samples

To examine the DWH spatiotemporal extent of the oil spill from in situ measurements, we used the BP GSD (26, 27), which includes more than 25,000 water and sediment samples from at least 67 response and Natural Resource Damage Assessment studies collected by multiple response agencies, trustees, and BP (26, 27). We used TPH (C9 to C44) and PAH50 concentrations (ppb), representing the summed concentrations of the 50 main PAH components (table S1). Water sample data included only unfiltered samples; therefore, concentrations account for both the dissolved and the particulate phases of hydrocarbons. The temporal window for the current study ranged from 5 May 2010 (the earliest available sample) to 3 October 2010, spanning across 16 to 167 days after blowout. Background PAH thresholds for water and sediment samples were 0.056 ppb (50) and 300 ppb (52), respectively.

The spatial extent of the modeled oil spill was compared with the distribution of higher-than-background water and sediment concentrations from the GSD. Predicted spill extents that coincided with observed higher-than-background oil concentrations represent a positive correspondence between modeled extent and in situ measurements. In addition, we compared the spatial distribution of the GSD data to the NOAA DIVER dataset (47).

Satellite and areal imagery: NESDIS anomaly

NESDIS anomaly data represent SAR, high-resolution visible/multispectral imagery, and/or other types of ancillary data sources representing the estimated areas with oil slicks on the water surface during the DWH (, using mainly SAR from the Center of Southeast Tropical Advanced Remote Sensing (CSTARS). The NESDIS anomaly together with daily fishery closures data were downloaded from the Environmental Response and Management Application web page ( ArcMap (Esri) software was used for computing the area of the NESDIS anomaly daily composites and for computing the daily intersected regions between the NESDIS footprint and the fishery closures. State waters that covered less than 2% of the domain were omitted from the spatial analysis because of limited data availability regarding their closures; only federal closures were, thus, considered.

Oil transport model: Oil-CMS

We implemented the existing oil application (17) of the CMS (53) to compute the transport and fate of the live oil spilled during the DWH blowout. The oil-CMS performed Lagrangian particle tracking of oil droplets released at the trap height above the well location. Particle transport calculations considered ocean 3D currents, temperature, salinity, multifractional droplet buoyancy, biodegradation, and surface oil evaporation (17, 20, 5355). The fourth-order Runge-Kutta integration scheme formed the basis for particle advection in the model. Computations of the vertical terminal velocity of a droplet were based on the droplet’s density and size, its Reynolds number, as well as on water temperature, salinity, density, and kinematic viscosity (17, 20, 53, 54).

The model output was saved every 2 hours and included oil droplets’ effective density, size, location, and depth. The oil-CMS horizontal grid spacing was 0.04° and included 20 vertical layers. The oil-CMS applied a multifraction droplet approach in which each droplet included multiple hydrocarbon fractions (54). The biodegradation dynamics were based on high-pressure experiments and applied different decay rates for the different fractions (54). This allowed accounting of dissolution processes where the droplets shrink with the partitioning of the oil compounds in the water column (56). Postprocessing algorithms translated model output into oil concentrations (18, 49, 57, 58).

Oil-CMS experimental setup

Three thousand oil droplets were released above the blowout location (28.736°N, 88.365°W) every 2 hours for 87 days until 15 July 2010, the day of successful capping of the oil well. The release depth was 1222 m, i.e., 300 m above the Macondo well depth, which was the estimated height of oil and gas separation above the wellhead (17). Initial droplet diameters were drawn from a uniform distribution between 1 and 500 μm. Two other DSDs, 80 to 2210 and 1 to 2400 μm, were simulated for spatial extent comparison (section S1). The postprocessing algorithm modified the obtained concentrations to fit to a log-normal curve of the DSDs. Each droplet released by the oil-CMS model contains three pseudocomponents (fractions) accounting for the differential oil densities as follows: 10% of light oil with a density of 800 kg m−3, 75% of intermediate oil with a density of 840 kg m−3, and 15% of heavy oil with a density of 950 kg m−3. The biodegradation half-life rates for the light, intermediate, and heavy fractions were set to 30, 40, and 180 hours based on laboratory and observation studies (54). Evaporation half-life rate was set to 250 hours (17, 55). Horizontal diffusion was set to 10 m2 s−1.

Ocean hydrodynamic forcing for the present study used daily output from the HYCOM for the GoM region on a 0.04° horizontal grid, including 40 vertical levels spanning from the surface to 5500 m. HYCOM uses data assimilation using the Navy Coupled Ocean Data Assimilation, which assimilates available satellite altimeter and sea surface temperature observations, as well as available temperature and salinity profiles from moored buoys and drifting hydrographic floats. The HYCOM output variables used for oil-CMS simulations included horizontal and vertical velocity components, temperature, and salinity.

The simulation included parameterization of the surface wind drift effects (20). Wind stress components from the 0.5° Navy Operational Global Atmospheric Prediction System were interpolated into the HYCOM GoM 0.04° grid, and 3% of their values were added to the top-level ocean velocity horizontal components, taking into account the wind stress rotation. The corrected ocean velocity fields were then implemented in the oil-CMS.

Flow rate was modeled as a simplified case of a constant flow rate during the oil spill event. The estimated 7.3 × 105 tons of oil, which corresponds to ~5.0 million barrels of oil (59), were represented by a total of 3.132 × 106 droplets. These values translate into 233 kg of oil represented by a single oil droplet at each release time in the oil-CMS model. We further computed the scaling factor for each droplet as the ratio between the current and the initial masses. Mass was estimated from the droplet diameter and effective density. Oil mass at each output time was then scaled to obtain effective oil mass at a given location and then summed for all the droplets found in each postprocessing domain 3D grid cell and at a given time step. Last, since TPH accounts for ~97% of spilled oil (23), we multiplied the oil mass by a factor of 0.97 to obtain TPH. The 4D oil concentration model output of oil-CMS can be found online in (60). A similar experimental setup was applied in (18, 49, 57, 58).

Spatial toxicity and visibility computation

We used MATLAB 2017a (MathWorks Inc.) to compute the spatial analyses related to the oil-CMS output. The spatial toxicity-visibility thresholds were based on the GSD analysis (see the “Quantitative and statistical analyses” section above), i.e., invisible and nontoxic (background level < PAH < toxicity threshold), invisible and toxic (toxicity threshold < PAH < visibility threshold), and visible and toxic (PAH > visibility threshold). We used a toxicity threshold of PAH = 0.5 ppb for surface waters (depth, 0 to 1 m) and 1 ppb for the rest of the water column (depth, >1 m) (4), and background level of PAH concentration of 0.056 ppb (50). The surface toxicity threshold is based on photoinduced toxicity, i.e., the combined effect of PAH and UV, whereas the toxicity threshold for the water column is based both on the decreased effect of photoenhanced toxicity with depth [Fig. 4.3-10 in (4)], as well as on cardiac toxicity and developmental abnormalities (4, 40). The latter was reviewed in (4), showing that fish larvae developed edema, decreased heart rates, and developmental abnormalities at low oil concentrations, in some cases at PAH < 1 ppb [Table 4.3-5 in (4)]. Similarly, to compute the exposure duration, we summed the number of days in which a given grid cell was exposed to toxic concentrations for the surface and the water column.

Species-specific spatial expression of toxic concentrations according to LC50 values was obtained from photoinduced toxicity studies (2831, 33) for early life stages of mahi mahi (Coryphaena hippurus), blue crab (Callinectes sapidus), red drum (Sciaenops ocellatus), speckled sea trout (Cynoscion nebulosus), and fiddler crab (Uca longisignalis). For this analysis, we considered surface waters only (0 to 1 m), where photoinduced toxicity is the most potent, although this should not be taken as the extent of the depth range to which photoinduced toxicity can occur.

The visibility threshold was obtained from the GSD decay dynamics at the surface (see the “Quantitative and statistical analyses” section). We assumed that oil is undetectable by satellites in waters deeper than 1 m regardless of the oil concentration (1).

Spatiotemporal dynamics of the spill

To compute the spatial dynamics of the DWH spill across time, we examined the following variables: oil mass, the modeled area of the oil spill, the modeled area of the toxic extents, and the area of the NESDIS anomaly footprint. Specifically, we compared the daily amount captured by the fishery closures to the total daily amount across the domain. Here, we considered both the surface (0 to 1 m) and the water column (depth, >1 m) for toxicity computation. The analysis applied here expands the spatial analysis applied in (18).


Supplementary material for this article is available at

Section S1. Comparing oil-CMS outputs of three possible DSDs

Section S2. The effectiveness of the seafood sampling scheme in determining the absence of contaminated seafood

Fig. S1. Spatial representation of the satellite footprint and relevant regions in the domain.

Fig. S2. Cumulative mass (in kilograms) across time and space.

Fig. S3. In situ water and sediment PAH concentrations from BP GSD (26, 27) and NOAA DIVER (47) data bases.

Fig. S4. Oil spill extents under different DSDs.

Fig. S5. Contaminated area in which seafood was monitored but was not closed for fishing.

Table S1. The chemical components of PAH50 from the GSD.

Table S2. List of photoinduced experiments used for the spatial analysis for the following organisms: mahi mahi (C. hippurus), blue crab (C. sapidus), red drum (S. ocellatus), speckled sea trout (C. nebulosus), and fiddler crab (U. longisignalis).

Movie S1. Daily oil spill extents from the oil-CMS transport model.

Movie S2. Oil spill extents under different DSDs.

Data file S1. R code S1.

Data file S2. DataPAH_TPH_surface.csv.

Data file S3. DataFrom_Bejerano_et_al_2013.csv.

Data file S4. DataPAH_TPH_water_column.csv.

Data file S5. SatSignalDecay.csv.

Data file S6. DataPAH_167d.csv.

References (6164)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the Paris lab and R. Faillettaz for the useful comments and suggestions. Funding: Support was provided by the National Academies of Sciences–Gulf Research Program (NAS-GRP) award: Understanding oil spill impacts on fishing communities of the Gulf of Mexico: From Deepwater Horizon to future spill scenarios. Author contributions: All authors participated in the design and the writing of the paper. C.B.P. and N.P. developed the oil-CMS application. S.B.J. supported the chemical aspect of this paper. I.B. performed the quantitative and spatial analyses of this paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article