Potential influence of the Deepwater Horizon oil spill on phytoplankton primary productivity in the northern Gulf of Mexico

Nine years after the Deepwater Horizon (DwH) oil spill (20 April–15 July 2010), the recovery of primary productivity at the ocean surface remains to be investigated. Here, we used the normalized fluorescence line height (nFLH) from the Moderate Resolution Imaging Spectroradiometer as an indicator of chlorophyll a concentration (Chl a). First, from the spatiotemporal variations of nFLH between 2001 and 2017, a reduction of nFLH after the DwH oil spill was observed (for a relatively long period, from 2011 to 2014). Second, a stepwise multiple regression model was used to examine which of the following environmental factors could explain the annual variations in nFLH: river discharge, total nitrogen load, total phosphorus load, photosynthetically available radiation, sea surface temperature and wind speed. Results show that river discharge, sea surface temperature and wind speed are the primary factors that regulated the annual nFLH variations in the DwH area during the pre-spill years. In contrast, this same model could not explain the reduction of nFLH for the four years after the DwH oil spill. After 2015, nFLH appears to have resumed to the pre-spill concentrations. Here we suggest that the nFLH reduction between 2011 and 2014 could have originated from the DwH oil spill, although the exact mechanism is yet to be determined.


Introduction
The Deepwater Horizon (DwH) 6 Girard and Fisher (2018) demonstrated-using high-definition imagery datathat deep-sea corals were heavily impacted and had not recovered by 2017. Using state-structured models, Girard et al (2018) suggested that the complete recovery of corals will take up to three decades (depending on the initial level of impact). As for the food web, the oil spill may have led to a significant decrease of the fishery yield which could take more than 30 years (especially for some slowly-growing populations) to recover fully (Ainsworth et al 2018). Tatariw et al (2018) found that the salt marsh Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 6 A full list of the acronyms used in study is summarized in table S1 is available online at stacks.iop.org/ERL/14/094018/mmedia. denitrification capacity (over the moderate and heavy oiling areas) had not yet recovered six years after the DwH oil spill. Long term studies are required-especially given the long life cycles of many of the organisms of interest, and the interannual variability in the meteorological and environmental conditions of the Gulf of Mexico.
The impacts of the DwH spill on phytoplankton population dynamics in the northern Gulf of Mexico have been under debate since the accident occurred . Hu et al (2011) suggested that phytoplankton productivity may have been stimulated by the DwH oil spill in the short term (weeks-months). A more recent study, however, found that the Mississippi River discharge may have also accelerated the phytoplankton growth after the spill (O'Connor et al 2016), which is consistent with previous observations of primary productivity in the region (Quigg et al 2011). Bretherton et al (2018 mimicked the chemical conditions after the DwH oil spill and found that the biological and physiological responses of phytoplankton to crude oil vary with species. Although some members of the phytoplankton community respond positively to the addition of oil, this was found to not be the case when both oil and dispersant were present (Bretherton et al 2019). Moreover, the chemical structure of the crude oil, the concentration level, and various environmental factors (e.g. temperature, light, and nutrient level) also influence the phytoplankton response to oil and/or dispersant (Vargo et al 1982, Østgaard et al 1984, González et al 2009, Kamalanathan et al 2018. In general, most studies have focused on short-term responses of chlorophyll a (Chl a) concentrations to the DwH oil spill (daysweeks-months). An evaluation of the long-term response (multi-year) is still lacking.
Chl a, which is used as a proxy for the phytoplankton biomass and primary productivity in the ocean, is generally estimated globally using satellite data (Behrenfeld and Falkowski 1997). Approximately 40% of all carbon fixation on Earth comes from marine phytoplankton (Falkowski 1994). Compared to the traditional ship surveys and laboratory analysis-both of which are time-consuming, and limited to small spatial and temporal coverage-satellite remote sensing provides large-scale observations of the Chl a at near real-time (Behrenfeld and Falkowski 1997). To date, many satellite sensors have been used for monitoring ocean color, such as the Sea-viewing Wide Field-of-View Sensor The goal of this study was to evaluate whether the Chl a concentration level-a proxy for primary productivity-was affected by the DwH oil spill in the northern Gulf of Mexico over a multi-year time period (and, if so, has it recovered to pre-spill levels-and when)? First, the spatiotemporal variations of the MODIS nFLH were evaluated from 2001 to 2017. Then, a multiple regression model was developed by correlating the nFLH anomaly with the anomalies of six environmental variables. Finally, by comparing the MODIS nFLH anomaly with the estimated nFLH anomaly (using the statistical model), the potential influence of DwH on long-term productivity in the northern Gulf of Mexico was explored.

Study area
The northern Gulf of Mexico region where the ocean surface was once covered by the DwH oil spill was selected as the study area. Figure 1 shows the extent of the domain and the spatial distribution of 'oiling days'-that is, the total number of days that a surface location was covered by oil (during the period from 23 April to 11 August 2010). This was based on data downloaded from the Gulf of Mexico Environmental Response Management Application (https://erma. noaa.gov/gulfofmexico/erma.html). Specifically, oiling days were calculated based on satellite Synthetic Aperture Radar (SAR) image classifications (Garcia-Pineda et al 2013). Here, the footprint of the surface oiling is referred to as the 'overall-DwH area', while the portion that suffered the most severe impacts (i.e. which had more than 30 oiling days) is referred to as the 'central-DwH area'. The overall-DwH area includes the entire study domain of 96 278 km 2 , while the central-DwH area accounts for 7.4% (7121 km 2 ) of this. By comparing the Chl a in these two areas, the oil spill effects were evaluated across spatial scales. . The nFLH is defined as the difference between the water-leaving radiance at 678 nm and a linearly interpolated waterleaving radiance at two surrounding bands (667 and 748 nm) (Behrenfeld et al 2009). Monthly MODIS nFLH data with a spatial resolution of 4 km (representing the period from January 2001 to December 2017) were obtained from the NASA Ocean Color website (https://oceancolor.gsfc.nasa.gov/) (data accessed in spring 2018), and then delineated for the study area (Gao and Li 2019).

Environmental factors
In this study, we selected six physical and chemical environmental variables that may influence Chl a in the DwH area: river discharge (Q), sea surface temperature (SST), photosynthetically active radiation (PAR), wind speed (WS), total phosphorus load (TP), and total nitrogen load (TN). The daily river discharge data for the Mississippi River at Tarbert Landing (Gage ID: 01100Q) and the Atchafalaya River at Simmesport (Gage ID: 03045Q) from 2001 to 2017 were obtained from the US Army Corps of Engineers (http:// rivergages.mvr.usace.army.mil). Monthly discharge values were calculated for these two stations, and the values from the two stations were summed to represent the total monthly river discharge into the Gulf of Mexico. The 4 km monthly SST and PAR data were obtained from the NASA Ocean Color website (https://oceancolor.gsfc.nasa.gov/), and were further processed into monthly time series for the overall-DwH area (Gao and Li 2019). Wind speed data were collected from the NOAA National Data Buoy Center (https:// ndbc.noaa.gov/). Since the Deepwater Horizon station (Station ID: 42872) was destroyed on 20 April, 2010, the wind data were acquired from the nearest station-the Luke offshore test platform (Station ID: 42040). The hourly wind speeds were aggregated to monthly mean values, except in cases where all data were missing for a certain month. In these cases, the climatological mean values for that month were used. The monthly TP and TN loads were obtained from the USGS stations located at the Atchafalaya and Mississippi Rivers (https:// nrtwq.usgs.gov/mississippi_loads/#/GULF) (Lee et al 2017). The monthly load amount is the product of constituent concentration and discharge integrated over time. The time series of these environmental factors are shown in figure S1.

Approach
To evaluate the variations of the nFLH (and its anomaly) before and after the DwH spill at different spatial and temporal scales, the monthly data were further processed into monthly climatology mean values (e.g. the climatology in January is the mean of the nFLH values of all Januaries from 2001 to 2017) and annual mean values. The monthly nFLH anomaly time series values were calculated by subtracting the monthly climatology mean values from the monthly values. To further explore whether the DwH oil spill had impacted the nFLH, a stepwise multiple regression model was developed after Ho (2006). The anomalies of the nFLH and the environmental variables were used in the procedure so that the correlation between MODIS observations and model estimations would not be affected by the seasonal cycles. First, the monthly climatological values for these six environmental factors were derived, and then the anomaly time series were generated. In the stepwise multiple regression model, the monthly nFLH anomaly was set as the dependent variable, while the anomalies of the environmental forcings were the independent variables. For each of the independent variables in the multiple regression model, only those terms that passed the significance test (with p-values<0.05) were adopted. The nFLH anomaly values estimated by the model were then compared with the MODIS observations. Specifically, the correlation (R), bias, and standard error (SE) were used to evaluate the performance of the regression model.  coast, but not in the central-DwH area. These anomalies may be due to the high discharge (Q) in these two years ( figure S2). Regardless, the locations of these anomalies indicate that primary productivity along the coast, which is close to the outlet and is rich in nutrients from freshwater inflows, is more sensitive than that in the central-DwH area.
Moreover, the nFLH anomaly patterns in August were examined from 2001 to 2017 ( figure 5). An apparent positive anomaly can be observed surrounding the central-DwH area in August 2010. This also  agrees with the findings by Hu et al (2011) which suggests that Chl a was stimulated by the DwH oil spill shortly after it occurred. In addition, distinct positive anomalies also occurred within the central-DwH area in August 2008, but there was no apparent positive anomaly that emerged in this area for August in the years after 2010. Since SST is negatively correlated to nFLH, the positive anomalies in 2008 can be explained by the fact that the SST in that year was the lowest for the study period ( figure S3). However, the SST in 2010 was above average, which supports the above point about the oil spill stimulating Chl a.
In addition, obvious negative anomalies were present from 2011 to 2013 in the central-DwH area (either with regard to the annual anomaly or the August anomaly), indicating that the oil spill may not have stimulated the phytoplankton in the relatively long term. The US drought of 2012 may also have exacerbated the negative anomaly by reducing the nutrient loading into the Gulf of Mexico (Wetz et al 2011), but this effect might not have lasted for a very long time (as the annual mean Q in 2013 recovered to the climatological value of the period from 2001 to 2017; see figure S2).

The relationship between nFLH and environmental forcing factors
To examine the potential effects of the oil spill on the nFLH over the DwH area, the MODIS observed nFLH anomaly data were compared with their counterpart from the stepwise multiple regression model. Through the regression procedure, three variables-the anomalies of Q, SST, and WS-were identified as the significant predictors for nFLH anomaly. The stepwise regression results are summarized in table S2, with the coefficients listed in table S3. The nFLH anomaly is most related to the anomaly of Q, which shows a significant correlation (R=0.31, p=0.000). By incorporating the anomalies of SST and WS in a stepwise manner, the R values increased to 0.42 and 0.46, while the standard errors decreased to 0.0180 and 0.0177 mW cm −2 μm −1 sr −1 , respectively. This selected model (with three independent variables) is thus referred to as Model_3 hereafter.
The biases of the Model_3 estimated monthly nFLH anomaly values from the MODIS nFLH anomaly during the study period are shown in figure S4. The nFLH anomaly is clearly overestimated from 2011 to 2014 (positive bias), which may have been caused by the DwH oil spill. To evaluate the performance of Model_3 on an annual basis, we calculated the annual mean bias (

Discussion
The effects of environmental forcing factors on phytoplankton vary across coastal areas (Green and Gould 2008, Hu et al 2011, Quigg et al 2011, Le et al 2014. Based on the stepwise multiple regression analysis of the anomaly data, Q, SST, and WS were identified as the most significant factors regulating the anomaly of Chl a (as indicated by nFLH) in the DwH area. According to the model results, Q is the most important factor-which agrees with the the findings by Le et al (2014). However, previous studies have indicated that SST and WS could impact primary productivity in the Gulf of Mexico (Agawin et al 2000, Green and Gould 2008, Bianchi et al 2010. It has been reported that nitrogen and phosphorous are the two key limiting nutrients for primary productivity in the northern Gulf of Mexico (Dagg and Breed 2003, Lohrenz et al 1997, Lohrenz et al 2008, Quigg et al 2011. Although total nitrogen and total phosphorous were not included explicitly in the model, their impacts were implicitly included because of their high correlations with Q ( figure S5). The relatively low correlation between the model estimates and the MODIS observations (R=0.46) can be attributed to two sources of uncertainties. The first are the uncertainties of the input data. The data associated with the model were collected via different approaches (e.g. remote sensing, in situ) and are of various spatiotemporal resolutions. These inconsistencies unavoidably have impacted the correlation negatively. The second source is related to the environmental variables which were not included in the model. Even though our model included the principal environmental drivers that regulate the nFLH in the DwH area, there are still other factors that influence primary productivity given the complexity of the marine ecosystem (e.g. loop currents, tropical storms, etc).
The biases between the modeled and observed nFLH anomalies have allowed us to explore the effects of environmental factors beyond those from the regression model (with the assumption that the input data uncertainties only add noise to the biases). Figure 6 clearly shows that the model experienced a significant change point in 2010. A mild underestimation from 2001 to 2010 is contrasted with a strong overestimation (2.9 times larger than the magnitude of the underestimation) from 2011 to 2014. Since the switch from under-to over-estimation began in 2010 and lasted for several years-and the most common environmental drivers have been represented in the model-we propose that this suppression effect might be attributed to the DwH oil spill. This argument also agrees with the study by Parsons et al (2015) which found that chlorophyll biomass was 85% lower in 2010 compared to the baseline of the previous 20 years (primarily due to the lower quantity of phytoflagellate, which decreased by 95%). Given that the DwH oil spill may have stimulated some phytoplankton while inhibiting others , there was a short-term increase in Chl a concentrations observed after the incident (Hu et al 2011). While some reports have found that the surface water of the Gulf of Mexico has recovered well from the DwH oil spill (Ferris 2017 and Franz 2017), D'Sa et al (2016) found fluorescence intensities and dissolved organic carbon concentrations three years after the DwH oil spill suggestive of a potential longterm persistence of oil in the dissolved organic carbon pool in the northern Gulf of Mexico. After 2015, primary productivity appears to have resumed to the prespill concentration levels. Thus, collectively the results from these studies offer insights about the responses of the primary producers.
While the surface oil was transported away from the DwH area by surface ocean currents and cleanup efforts, the multi-year reduction of primary productivity may be attributed to spill residue in the underlying sediments. For instance, Duan et al (2018) investigated the persistence of residual oil in the sediment at Bay Jimmy five years after the spill. They found that the concentrations of total petroleum hydrocarbons, n-alkanes, and polycyclic aromatic hydrocarbons increased significantly within the sediment after the oil spill. Previous studies have shown that polycyclic aromatic hydrocarbons had depressive effects on the phytoplankton biomass (Marwood et al 1999, Sargian et al 2005, Pelletier et al 2006, Kamalanathan et al 2018. Five years after the spill, most of the n-alkanes and polycyclic aromatic hydrocarbons had degraded and recovered to pre-spill levels. Even though the concentrations of total petroleum hydrocarbons were still relatively high, they had decreased by 97% compared to the level at 1.5 years after the spill. These, and the in situ observations reported in Franz (2017) and Ferris (2017), support our finding that primary productivity started to recover around 2015.
It is also worth noting that this study is focused on surface water (i.e. the first few meters) because of the strong light attenuation in the wavelengths used to calculate MODIS nFLH (i.e. 667, 678, and 748 nm). Moreover, our findings are based on a statistical model, which does not represent physical mechanisms. Additionally, although earlier analysis by NASA showed sensor degradation after 2011, the degradation was corrected through improved time-dependent sensor calibration and vicarious calibration in NASA's data reprocessing version 2014.0 (Hu et al  2015). Therefore, the temporal patterns of MODIS nFLH are unlikely affected by sensor calibration changes.
To date, continuous long-term field observations before and after the spill are still lacking, which hamper the understanding of the variability within this large area. In the absence of such long-term field data, this study shows an alternative way to evaluate the long-term changes of primary productivity (via Chl a) in the DwH area, which may be adapted for other oil spill events.

Conclusions
The most significant finding from this study is thatalthough annual changes in MODIS nFLH (used as a proxy for Chl a, an indicator for primary productivity) before 2010 can be explained by environmental forcing factors (river discharge, sea surface temperature, and wind speed)-the same explanation does not hold between 2011 and 2014. The behavior between 2011 and 2014 is speculated to be a result of the long-term, chronic impact of the DwH oil spill. Although it is impossible to verify this hypothesis due to a lack of continuous field measurements, this study represents the first attempt to use long-term satellite data to evaluate the potential chronical effects of the DwH oil spill on primary productivity. This also suggests the importance of continuous field measurements to help pinpoint the reasons behind the changes of phytoplankton productivity in future oil spills.