ENSO drives near-surface oxygen and vertical habitat variability in the tropical Pacific

To accurately represent the oxygen conditions experienced by marine animals such as tuna, we report oxygen content in terms of partial pressures (pO₂) because it is pO₂, rather than dissolved oxygen concentration (O₂), that drives oxygen transfer and provision to animal tissue (Hofmann et al 2011, Seibel 2011). To compute pO₂ from input variables of temperature, salinity, and O₂, we first convert O₂ to percent oxygen saturation (Garcia and Gordon 1992). We then divide percent oxygen saturation by the fractional atmospheric concentration of oxygen (21%) and as a final step, correct for hydrostatic pressure at depth (Enns et al 1965). An additional advantage of using pO₂ is that because it is a function of percent oxygen saturation, it automatically accounts for changes in surface oxygen solubility owing to changes in sea surface temperatures (SST); in other words, pO₂ is not affected by variations in solubility, allowing us to rule out this potential confounding factor driving variations in oxygen content. To describe changes in the availability of oxygenated vertical habitat in the upper-ocean to tuna and other pelagic predators with similar physiologies, we use THD, defined as the shallowest depth at which pO₂ first falls below 15 kPa. This 15 kPa pO₂ threshold is approximately equivalent to the 3.5 ml l⁻¹ dissolved concentration threshold (assuming a salinity of 32.5 psu, a temperature of 23.5 °C, and depth of 0 m, as in the tank experiment by Gooding et al 1981) typically used to describe hypoxic conditions considered dangerous for skipjack and yellowfin tuna (Ingham et al 1977, Barkley et al 1978, Evans et al 1981, Gooding et al 1981, Bushnell et al 1990, Prince and Goodyear 2006, Arrizabalaga et al 2015), which make up the vast majority of tuna caught in the western and central tropical Pacific (Pons et al 2017) (see supplementary figure 1 available online at stacks.iop.org/ERL/14/064020/mmedia, which shows more precisely how this 15 kPa pO₂ threshold translates into dissolved oxygen concentrations in units of both ml l⁻¹ and µmol kg⁻¹ as a function of temperature, salinity, and depth across values of these properties that are typically encountered by skipjack and yellowfin tuna).

In situ profiles of near-surface (0–700 m) temperature, salinity, and O₂ collected between January 1955 and May 2017 were downloaded from the World Ocean Database (WOD) at https://nodc.noaa.gov/OC5/SELECT/dbsearch/dbsearch.html on 25 June, 2017. WOD compiles uniform and quality-controlled data throughout the global ocean from buoys, ships, gliders, and floats (Boyer et al 2013). All available datasets within WOD were collected and used including data from ocean stations, high resolution CTDs, expendable bathythermographs (XBT), mechanical bathythermographs (MBT), profiling floats (including Argo), drifting buoys, moored buoys, autonomous pinniped bathythermographs, undulating oceanographic recorders, buckets, thermosalinographs, and gliders. XBT/MBT temperatures were corrected using Levitus et al (2009). Any flagged data were excluded.

Raw profiles of temperature, salinity, and O₂ were binned onto 5-by-5 horizontal degree, monthly mean maps with standard WOD vertical levels (5 m depth resolution down to 100 m, 25 m depth resolution between 100 and 500 m, and 100 m depth resolution beneath 500 m). These monthly mean maps were then used to calculate corresponding pO₂ maps at grid points where all three input variables were available. The mean seasonal cycle was then computed at each grid point for all variables. No lower threshold was placed on the number of years required to compute a climatological monthly mean (e.g. if a given grid cell had only one monthly mean data point in January over the entire time series, then that data point became January's climatological mean at that grid cell). Anomalies were then calculated by subtracting the mean seasonal cycle at each corresponding grid point. The resulting 4D (x, y, z, time) anomaly maps were used for all subsequent map-based analyses (figures 1–4; supplementary figures 3–6, 8, 9), where missing data is denoted with empty grid cells.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For all other analyses (table 1; figure 5; supplementary figure 7), raw profiles of the variables were grouped by month and region of interest, including the Western Equatorial Pacific box (WEP), the Eastern Equatorial Pacific box (EEP), and exclusive economic zones (EEZs) (defined in figures 2(a) and (c)). Raw profiles of pO₂ were computed whenever corresponding raw profiles of all three input variables were available. The mean seasonal cycle was then computed within each region of interest, and anomaly profiles were subsequently calculated by subtracting the mean seasonal cycle within each corresponding region.

Zoom In Zoom Out Reset image size

All aforementioned mapped and regional anomalies were also computed by subtracting out World Ocean Atlas 2013 climatologies (Locarnini et al 2013, Zweng et al 2013, Garcia et al 2013), which did not yield any notable differences in our results. To maintain data source consistency, we thus show only anomalies computed by subtracting out the WOD-computed climatologies discussed above. See supplementary information—data coverage for a discussion of the effects of spatially and temporally variable amounts of data on our results.

One potential driver of near-surface oxygen variability is changes in the vertical structure of differentiated water masses. Subsurface temperature data is much more plentiful than that of O₂ (see supplementary information—data coverage), and thermo-and oxycline depths covary strongly throughout the tropical Pacific; we therefore use variations in TCD alone to represent the vertical movements of water masses with distinct temperature and oxygen signatures. TCD was calculated from depth-resolved monthly mean temperature maps (figures 1, 3–4; supplementary figures 4, 6(d), 9(b), (d), (f), (h)) and profiles (table 1) using the variable representative isotherm method, as recommended for tropical waters by Fiedler (2010), who compared several objective methods of computing TCD with empirical TCD estimates made by eye. Other objective methods of computing TCD, such as the depth of the 20 °C isotherm or the depth of the maximum vertical gradient in temperature, led to similar results here. In addition to its computational efficiency, the advantage of the variable representative isotherm method is that it accounts for the wide range of temperatures in which the thermocline occurs, since there is no single isotherm that well represents empirical TCDs (e.g. empirical thermocline temperatures range from 16 °C to 26 °C within the tropical Pacific alone) (Fiedler 2010). According to the variable representative isotherm method, the thermocline spans the water column from the base of the mixed layer to the depth at which temperature has dropped halfway toward the temperature at 400 m; TCD is then defined as the midpoint of this layer. Mathematically, the isotherm at which the TCD occurs is as follows:

$$\begin{eqnarray*}&&\begin{array}{l}{T}\left({\rm{thermocline}}\,{\rm{representative}}\,{\rm{isotherm}}\right)\\ \,=\,{T}\left({\rm{MLD}}\right)-0.25\left[{T}\left({\rm{MLD}}\right)-{T}\left(400\,{\rm{m}}\right)\right],\end{array}\end{eqnarray*}$$

where T(MLD) = SST − 0.8 °C.

All variables (either binned onto 5°-by-5° monthly mean maps or spatially averaged over specified regions) were examined for ENSO-related variability using the Oceanic Niño Index (ONI). ONI is calculated as the 3 month running mean of SST anomalies in the Niño 3.4 region (5°N–5°S, 120°–170°W) (Huang et al 2016). El Niño and La Niña periods occur when ONI is greater than 0.5 °C and less than −0.5 °C, respectively, for at least five consecutive 3 month running mean periods (supplementary figure 2). Time series of the ONI were downloaded from http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php on 4 June, 2018.

To avoid incorrect or overstated interpretations, we control our false discovery rate (FDR) whenever analyzing the results of multiple-hypothesis tests. The FDR is defined as the statistically expected fraction of local (i.e. at a single grid point in the case of map analyses or within a single region in the case of comparisons across EEZs or other specified areas) null hypothesis test rejections for which the respective null hypotheses are actually true. Controlling FDR in the context of multiple-hypothesis tests thus requires smaller p-values compared to single hypothesis testing in order to reject local null hypotheses. The procedure for determining significance under a controlled FDR is as follows (Wilks 2016): (1) sort the collection of p-values from N (i.e. the total number of grid points in the case of map analyses and the total number of regions in the case of comparisons across EEZs or other specified areas) local hypothesis tests _pi, with i = 1, ..., N, in ascending order. (2) Reject local null hypotheses if their respective p-values are no larger than a threshold level p_FDR, which is equal to the largest _pi that is no larger than the fraction of α_FDR specified by i/N. p_FDR is thus calculated as follows: ${{p}}_{{\rm{FDR}}}=\,\max \left[{{p}}_{{i}}:{{p}}_{{i}}\leqslant \left({i}/{N}\right){\alpha }_{{\rm{FDR}}}\right],$ where α_FDR is the chosen FDR (expressed as a fraction).

This method assumes that the multiple local tests are statistically independent, but is also valid even when the results of the multiple tests are strongly correlated (such as is the case with geospatial data). Indeed, when spatial correlations are high, the achieved FDR will be even smaller (stricter) than the chosen FDR (Ventura et al 2004, Wilks 2006), so that the chosen α_FDR should be approximately double the desired level when analyzing highly spatially correlated data grids (Wilks 2016).

The MATLAB code required to make all of the figures and tables generated here can be found at https://doi.org/10.5281/zenodo.2648131.

Data in the form of *.mat files required to the make all of the figures and tables generated here can be found at http://doi.org/10.5281/zenodo.2648124.

In the equatorial Pacific (between 7.5°N and 7.5°S), mean near-surface pO₂ differs substantially from east to west, with higher pO₂ values at any given depth in the west compared to the east (figure 1(a); supplementary figures 3(a), (b)). The mean oxycline, where pO₂ changes most quickly with depth, is located at around 100 m depth near the western boundary, descends to 150 m near the dateline, and shoals to around 50 m near the eastern boundary (figure 1(a)). Mean THD (figure 1 solid gray contour) lies close to the oxycline across the entire basin. The within-oxycline vertical pO₂ gradient is strongest in the east, coincident with the top of the shallow oxygen minimum zone there.

Monthly variability in equatorial Pacific pO₂ shows distinct east–west differences as well (figure 1(b); supplementary figures 3(c), (d)). Near-surface pO₂ variability is greater in the east due to the presence of sharp vertical gradients in pO₂, with maximum variability occurring within the relatively shallow oxycline between 140°W and 80°W (figure 1(b); supplementary figures 3(c), (d)). pO₂ variability in the west, between the western boundary and 160°W, is about half as large as that in the east and is again greatest surrounding the local oxycline (figure 1(b); supplementary figures 3(c), (d)).

Interannual variability in near-surface pO₂ is driven primarily by ENSO throughout the tropical Pacific. During El Niño, pO₂ decreases to the west of 160°W and increases to the east of 160°W throughout the upper water column (figure 1(c); supplementary figures 3(e), (f)). The opposite occurs during La Niña (figure 1(d); supplementary figures 3(g), (h)). In the east above the oxygen minimum zone, these results are in agreement with both field (Gutiérrez et al 2008, Czeschel et al 2012, Llanillo et al 2013, Graco et al 2017, Stramma et al 2016) and model-based studies (Yang et al 2017). In the west, the effects of ENSO on near-surface oxygen conditions have not been well-characterized before this study. In both the east and west, ENSO-driven pO₂ variability is greatest within the vicinity of the oxycline (around 100–150 m depth in the west and 50 m depth in the east) (figures 1(c), (d)), as is the case for total interannual pO₂ variability (figure 1(b)). The ENSO-associated pO₂ variability dominates the total variability, in agreement with model results from Resplandy et al (2015), who showed that ENSO is the major driver of interannual variations in tropical Pacific air–sea oxygen fluxes within three Coupled Model Intercomparison Project 5 (CMIP5) Earth System Models (CESM1-BGC, GFDL-ESM2G, and GFDL-ESM2M).

Interannual variability in tropical Pacific THD is also driven primarily by ENSO. Relative ENSO-driven THD variations are greatest from the western boundary (around 125°E) to 160°W and from 130°W to the eastern boundary (around 80°W) (figures 2(b), (c)). Within the western region (125°E to 160°W), THD is as much as 22 m shallower during El Niño versus mean conditions, and 34 m shallower during El Niño versus La Niña conditions (figures 1(c), (d) contours). Averaged over the WEP defined in figure 2(a), THD variations compress vertical habitat space by 10.3% (14.5 m) during El Niño (figure 2(b); table 1) and expand it by 9.5% (13.3 m) during La Niña (figure 2(c); table 1), compared to mean conditions. Vertical habitat space is thus compressed by 27.8 m during El Niño relative to La Niña within the WEP (p < 0.05, unequal variances t-test). In the east, ENSO acts in the opposite direction, shoaling THD during La Niña and deepening it during El Niño. Averaged over the EEP defined in figure 2(a) and compared to mean conditions, hypoxia-induced vertical habitat compression of 2.9% (1.1 m) occurs during La Niña (figure 2(c); table 1), while expansion of 12.9% (4.9 m) occurs during El Niño (figure 2(b); table 1). Vertical habitat space is thus compressed by about 6 m during La Niña relative to El Niño within the EEP (p < 0.05, unequal variances t-test).

Based on the observations analyzed here, ENSO-related variations in oxygenated tuna vertical habitat space appear to be driven primarily by the up-and-down motions of distinct water masses separated by well-defined thermo- and oxyclines. ENSO has different effects on the vertical positioning of these water masses and therefore TCDs depending on the region. In the WEP (as defined in figure 2(a)), El Niño shoals TCD by 12.2% (18.1 m) (supplementary figure 4(b); table 1), while La Niña deepens TCD by 9.0% (13.3 m) compared to mean conditions (supplementary figure 4(c); table 1). TCD is thus approximately 31.5 m shallower during El Niño relative to La Niña within the WEP (p < 0.05, unequal variances t-test). Within the EEP (as defined in figure 2(a)), El Niño deepens TCD by 20.4% (10.5 m) (supplementary figure 4(b); table 1), while La Niña shoals TCD by 15.6% (8.1 m) compared to mean conditions (supplementary figure 4(c); table 1). TCD is thus about 18.5 m shallower during La Niña relative to El Niño within the EEP (p < 0.05, unequal variances t-test).

The strongest temporal correlations between pO₂ and TCD (figure 3(a)) occur at depths surrounding the mean thermocline, where both pO₂ and TCD variability are largest (figures 1(b)–(d); 3(a)). The prevalence of these strong positive correlations between pO₂ and TCD, in addition to the strong positive correlations between THD and TCD throughout the tropical Pacific (figure 3(b); table 1), suggest the following: as the thermocline deepens, hypoxic waters beneath the thermocline are pushed downward, expanding oxygenated vertical habitat space and increasing pO₂ in the upper water column; conversely, as the thermocline shoals, hypoxic waters beneath the thermocline are pulled upward, compressing oxygenated vertical habitat space and decreasing pO₂ in the upper water column.

Though first-order variations in THD anomalies are well-explained by the up-and-down motions of the thermocline and accompanying water masses, concurrent changes in pO₂ along isopycnals within the thermocline—caused by changes in primary production and physical transport—can either reinforce or dampen this primary effect (isopycnals are roughly equivalent to isothermals in the tropical Pacific because density below the mixed layer is primarily determined by temperature there). During El Niño, decreases in equatorial primary and export production in the east should decrease microbial oxygen consumption and thus increase pO₂ along isopycnals; in the west, increases in local production and export should lead to a decrease in pO₂ along isopycnals. During La Niña, these conditions should reverse (Chavez et al 1999, Behrenfeld, et al 2001, Lehodey 2001). Contrary to these expectations, however, isopycnal pO₂ within the thermocline decreases in the east and increases in the west during El Niño (figure 4(b); supplementary figure 5(b)), with the opposite occurring during La Niña (figure 4(c); supplementary figure 5(c)). The isopycnal pO₂ variations observed here are therefore likely caused by ENSO-associated changes in lateral ventilation and transport within the thermocline rather than local changes in primary production. ENSO-associated lateral transport changes are, in turn, likely driven by variations in the strength of subsurface zonal currents, the most important being the Equatorial Undercurrent (EUC). The EUC carries oxygen-rich water from the subtropics and western Pacific to the central and eastern Pacific along the thermocline (Stramma et al 2010), and undergoes large variations in strength associated with ENSO (Johnson et al 2000, Izumo 2005). During the 1997/98 El Niño, for instance, EUC mass transport (measured along the equator at 140°W) dropped to nearly zero from a mean of around 34 Sv (Izumo 2005). Decreases in EUC strength of this magnitude greatly reduce within-thermocline transport of oxygen from the western to the central and eastern equatorial Pacific during El Niño. Increases in EUC strength during La Niña, on the other hand, bring anomalously large amounts of oxygen from the west to the east, raising and lowering within-thermocline isopycnal pO₂ values in the east and west, respectively. During both ENSO phases, these transport-driven isopycnal oxygen changes dampen and oppose the effects of TCD variations, but because they have much a weaker impact on pO₂ and THD, the vertical movement of water masses remains the dominant factor influencing ENSO-associated pO₂ and THD anomalies.

The effects of transport-driven changes on ENSO-related THD and pO₂ variations are particularly small in the western Pacific (figures 4(b), (c)). As a result, ENSO-induced TCD variations alone can accurately predict concomitant THD variations, leading to comparable all-profile and ENSO-driven THD versus TCD regression coefficients within the WEP (table 1). All-profile THD versus TCD regression coefficients are computed by relating THD and TCD anomalies from all simultaneous profiles, and thus represent the magnitude of THD changes driven primarily by changes in TCD. ENSO-driven THD versus TCD regression coefficients are computed as the difference between El Niño and La Niña THD values divided by the difference between El Niño and La Niña TCD values, and thus represent the magnitude of ENSO-related THD changes driven by both ENSO-associated TCD and isopycnal pO₂ changes. In contrast to the situation in the west, ENSO-driven variations in eastern Pacific, within-thermocline isopycnal pO₂ noticeably dampen the THD variations induced by vertical shifts in water masses and TCD (figures 4(b), (c)). Because of this, variations in TCD alone (that is, not accounting for within-thermocline isopycnal pO₂ changes) slightly overpredict differences in THD between ENSO phases within the EEP, as evidenced by a larger all-profile THD versus TCD regression coefficient compared to the ENSO-driven one here (table 1). Despite the noticeable dampening effect of physical transport changes on pO₂ and THD variations in the EEP, the up-and-down movements of water masses remain the primary factor driving ENSO-related pO₂ and THD variations, as can be seen from the net direction and sign of these variations. In sum, based on observations alone, the vertical motions of the thermocline appear to be the main driver of ENSO-induced variations in near-surface pO₂ and THD throughout the tropical Pacific, while ventilation-driven variations in isopycnal pO₂ within the thermocline play an important (especially in the east) but secondary role. This is again in agreement with model results from Resplandy et al (2015), who found that the main drivers of interannual variations in tropical Pacific air–sea oxygen flux were changes in physical circulation and vertical transport rather than primary production.

In both the WEP and EEP regions (as defined in figure 2(a)), El Niño induces greater deviations of THD from mean conditions compared to La Niña (table 1), which is consistent with observed asymmetrical effects of ENSO on other well-studied climate phenomena (Hoerling et al 1997, Kang and Kug 2002, An et al 2005). Averaged over the WEP, El Niño shoals THD by 14.5 m, while La Niña only deepens THD by 13.3 m (table 1). In the EEP, El Niño deepens THD by 4.9 m, while La Niña only shoals THD and thus compresses tuna vertical habitat space by 1.1 m (table 1). These asymmetric variations in THD are caused by corresponding asymmetric variations in TCD between opposing ENSO phases, such that ENSO-driven TCD changes are also larger during El Niño compared to La Niña within both the WEP and EEP regions (table 1).

We now narrow our focus to the western tropical Pacific both because it is the most tuna-rich region in the world and ENSO-driven oxygenated vertical habitat variations are most pronounced here. THD variations of the magnitude observed here can significantly affect catchability and CPUE using common industrial surface fishing gear such as longlines and purse seines, potentially resulting in overfishing in years when vertical habitat space is compressed (Green 1967, Barkley et al 1978, Evans et al 1981, Nakano et al 1997, Prince and Goodyear 2006, Bigelow and Maunder 2007, Prince et al 2010). Though tuna occur throughout the Western and Central Pacific Ocean, approximately 60% of the total tuna catch in this region is taken from within the EEZs of only eight Pacific island nations that jointly enacted the Nauru Agreement to better manage their fisheries (Havice 2013). Of the eight Parties to the Nauru Agreement (PNA), the EEZs belonging to Micronesia and the Marshall Islands undergo the greatest reductions in mean vertical habitat space during El Niño (14.0 m or 9.9% and 13.4 m or 7.9%, respectively), while the EEZs belonging to Palau, Papua New Guinea, and Nauru exhibit the greatest expansions in mean vertical habitat space during La Niña (14.3 m or 12.5%, 11.2 m or 8.1%, and 14.5 m or 11.1%, respectively) (figure 5). The EEZs of Palau, Micronesia, Nauru, and the Marshall Islands undergo the largest total mean ENSO-driven variations in vertical habitat space, with 19.5 m, 23.9 m, 19.5 and 19.3 m differences in THD between El Niño and La Niña, respectively (figure 5).