Trends and decadal oscillations of oxygen and nutrients at 50 to 300 m depth in the equatorial and North Pacific

A strong oxygen-deficient layer is located in the upper layers of the tropical Pacific Ocean and deeper in the North Pacific. Processes related to climate change (upper-ocean warming, reduced ventilation) are expected to change ocean oxygen and nutrient inventories. In most ocean basins, a decrease in oxygen (“deoxygenation”) and an increase in nutrients have been observed in subsurface layers. Deoxygenation trends are not linear and there could be multiple influences on oxygen and nutrient trends and variability. Here oxygen and nutrient time series since 1950 in the Pacific Ocean were investigated at 50 to 300 m depth, as this layer provides critical pelagic habitat for biological communities. In addition to trends related to ocean warming the oxygen and nutrient trends show a strong influence of the Pacific Decadal Oscillation (PDO) in the tropical and the eastern Pacific, and the North Pacific Gyre Oscillation (NPGO) in particular in the North Pacific. In the Oyashio Region the PDO, the NPGO, the North Pacific Index (NPI) and an 18.6-year nodal tidal cycle overlay the long-term trend. In most eastern Pacific regions oxygen increases and nutrients decrease in the 50 to 300 m layer during the negative PDO phase, with opposite trends during the positive PDO phase. The PDO index encapsulates the major mode of sea surface temperature variability in the Pacific, and oxygen and nutrients trends throughout the basin can be described in the context of the PDO phases. El Niño and La Niña years often influence the oxygen and nutrient distribution during the event in the eastern tropical Pacific but do not have a multi-year influence on the trends.

1 Introduction

Oxygen and nutrient distribution are key parameters controlling marine ecosystems. How oxygen and nutrient concentrations vary and co-vary ultimately controls biogeochemical cycles. Globally, oxygen has been estimated to have decreased in the ocean by 2 % during the past 5 decades, likely caused by climate-change-related temperature increases, with the largest oxygen decrease in the North and equatorial Pacific (Schmidtko et al., 2017). Increasing sea surface temperatures reduce the solubility of oxygen in sea water and increase stratification, which leads to less convection of oxygen-rich water to subsurface layers. The global mean surface temperature (GISTemp) anomaly (e.g., Hansen et al., 2010) shows an increase of about 0.3 ^∘C from early 1900 to 1950, stagnant values between 1950 and 1976, and an increase of about 0.8 ^∘C from 1976 to 2015 (Fig. 1). Results from the UVic ECSM (University of Victoria Earth System Climate Model) model indicate that ocean oxygenation varies inversely with low-latitude surface wind stress (Ridder and England, 2014).

Figure 1 Global mean surface temperature anomaly from the 1951 to 1980 mean (GISTemp; peak of the annual bars; left scale) with La Niña years (blue bars), El Niño years (orange-brown bars), neutral years (grey bars) and the PDO index (solid line; right scale), with PDO phases marked (compiled by Miriam O'Brien for 1900 to 2015). The data sources used GISS NASA (temperature), Australian Bureau of Meteorology (ENSO years based on the Southern Oscillation Index), Japan Meteorological Society (PDO index) and Trenberth 2015 (PDO phases).

In the tropical eastern Pacific, two subsurface low-oxygen zones exist north and south of the Equator, with a pronounced minimum in oxygen at ∼100 to 500 m depth, and these are referred to as oxygen minimum zones (OMZs) or oxygen-deficient zones (ODZs). These OMZs are suboxic (oxygen concentrations below ∼4.5–10.0 µmol kg⁻¹; e.g., Karstensen et al., 2008; Stramma et al., 2008). In suboxic regions, nitrate and nitrite become involved in respiration processes such as denitrification or anammox (e.g., Kalvelage et al., 2013). Decreasing and, in a few areas of the Pacific, increasing oxygen content in the OMZ layer over the last 50 years has been described (e.g., Stramma et al., 2010). In the subarctic North Pacific, surface nutrient concentration decreased from 1975 to 2005 and is strongly correlated with a multidecadal increasing trend of sea surface temperature (SST) (Ono et al., 2008). Below the surface, however, oxygen decreased and nutrients increased in the subarctic Pacific pycnocline from the mid-1980s to around 2010 (Whitney et al., 2013). Nutrients would be expected to vary inversely with oxygen, if the dominant process was the remineralization of marine detritus (Whitney et al., 2013).

Climate modes influence oxygen and nutrient distributions. Because of the influence of SST on the solubility of oxygen and changes in convection and thermocline depth, the most prominent control on oxygen changes in the Pacific might be exerted by the Pacific Decadal Oscillation (PDO). A shoaling thermocline, such as occurs in the eastern Pacific during La Niña or a cool (negative) PDO state, enhances nutrient supply and organic matter export in the eastern Pacific while simultaneously increasing the fraction of that organic matter that is respired in the low-oxygen water of the uplifted thermocline. The opposite occurs during El Niño or a warm (positive) PDO state; a deeper thermocline reduces both export and respiration in low-oxygen water in the eastern Pacific, allowing the hypoxic water volume to shrink (Deutsch et al., 2011; Fig. S7). Previous syntheses of tropical and North Pacific physical, biological and chemical conditions during warm and cold PDO regimes have shown the far-reaching influence of the PDO on the Pacific Ocean, for example in controlling the out-of-phase relationship between sardine and anchovy populations in the eastern Pacific (Chavez et al., 2003). In a pattern similar to the PDO SST during the warm regime, lower nutrients are shown for the eastern equatorial and near-coastal tropical Pacific and higher nutrients in the northern Pacific, and vice versa for the cold regime (Chavez et al., 2003; their Fig. 3). A similar relationship to the PDO was also observed for surface (5–10 m) nutrient concentrations in the North Pacific (north of 10^∘ N; Yasunaka et al., 2016). Model simulations for the eastern Pacific Ocean for typical PDO positive conditions show a volume expansion of the suboxic regions by 7 % in 50 years due to a slowdown of the large-scale circulation related to the decrease in the intensity of the trade winds (Duteil et al., 2018). Other climate modes that could influence the oxygen and nutrient trends are the North Pacific Gyre Oscillation (NPGO; Di Lorenzo et al., 2008), North Pacific Index (NPI) and El Niño–Southern Oscillation (ENSO) events.

Here we use publicly available oxygen and nutrient data augmented with recent ship data to investigate the influence of decadal climate oscillations on the oxygen and nutrient variability in the tropical, eastern and northern Pacific Ocean where long time series are available. For the negative PDO phase we use data between 1950 and 1976. As there are large data gaps in the 1990s and early 2000s, we use data since 1977 for the warm PDO despite the variable PDO conditions after 1998. In addition, possible influence of the NPGO, of the NPI, of ENSO and of an 18.6-year oscillation (Royer, 1993) in the North Pacific is investigated.

2 Climate signals and data sets

2.1 Climate signals

Several climate signals are investigated here with regard to a possible influence on the oxygen and nutrient distribution and some basic details are listed here.

According to Deser et al. (2010) the leading empirical orthogonal function (EOF) of monthly SST anomalies over the North Pacific (after removing the global mean SST anomaly) and its associated principal component (PC) time series are termed the Pacific Decadal Oscillation (PDO) after Mantua et al. (1997).

The temperature pattern and the PDO time series (Fig. S1 in the Supplement) show that the PDO was negative from 1944 to 1976 with a stagnant temperature period, and they show warming temperatures since 1977 with PDO positive from 1977 to 1998 and PDO variable after 1998. At the time of transition from negative to positive PDO, a climate shift in the eastern and central North Pacific Ocean occurred in 1976–1977, which was caused by unique atmospheric anomalies acting over several months before the 1976–1977 winter (Miller et al., 1994). Despite a positive PDO index from 2002 to 2006 (Fig. 1), the period 1998 to 2013 is dominated by negative seasonal mean PDO indices and is typically considered to be a cool (negative PDO) phase (Trenberth, 2015).

The North Pacific Gyre Oscillation (NPGO) index tracks the changes in strength of the central and eastern branches of the North Pacific gyres and of the Kuroshio–Oyashio Extension. Like the PDO, it is a mode of decadal climate variability and it is defined as the second dominant mode of variability in sea surface height anomaly in the northeast Pacific over the region 25–62^∘ N, 180–110^∘ W (Di Lorenzo et al., 2008). Subsurface nutrient variability in the Gulf of Alaska, Line P in the eastern North Pacific at about 50^∘ N and the California Current System have been shown to be correlated with the NPGO (Di Lorenzo et al., 2009).

Figure 2 Distribution of the areas used and mean 50 to 300 m silicate, phosphate, nitrate and oxygen (all in micromoles per kilogram), in the Pacific Ocean.

In the North Pacific, the NPI (North Pacific Index; anomaly of the sea surface level pressure in the wintertime of the North Pacific (30–65^∘ N, 160^∘ E–140^∘ W; Minobe, 2000)) was introduced. In observations in the northwest Pacific an oscillation every 2 decades related to minima in 1962 and 1983 and maxima in 1971 and 1991 (Watanabe et al., 2003, their Fig. 2) has been described.

The ENSO cycle of alternating warm El Niño and cold La Niña events is the climate system's dominant year-to-year signal. ENSO originates in the tropical Pacific through interaction between the ocean and the atmosphere, but its environmental and socioeconomic impacts are felt worldwide (McPhaden et al., 2006). The 3-month running mean SST anomalies (ERSST.v5 SST anomalies) in the Niño 3.4 region (equatorial Pacific: 5^∘ N to 5^∘ S, 120 to 170^∘ W) of at least +0.5 ^∘C and lasting for at least five consecutive 3-month periods are defined as El Niño events and five consecutive 3-month periods of at least −0.5 ^∘C are defined as La Niña events (http://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php, last access: 13 May 2019). In the eastern Pacific, ENSO variability is most pronounced along the Equator and coastal Ecuador, Peru and Baja California (Wang and Fiedler, 2006). Coastal warming during El Niño is caused by downwelling Kelvin waves generated by mid-Pacific westerly wind anomalies that deepen the eastern thermocline, nutricline and oxycline and allow warming to occur (Kessler, 2006). Consequently, during El Niño events the upwelled water off the coast is warmer, more oxygen-replete and less nutrient-rich, while during La Niña events the water is colder, oxygen-poor and nutrient-rich (e.g., Graco et al., 2017).

2.2 Data sources by regions

The main hydrographic data set is similar to the one used and described in Schmidtko et al. (2017), relying on Hydrobase and World Ocean Database bottle data for nutrient data. Quality control and handling is described in Schmidtko et al. (2017) for oxygen and used here similarly for nutrients. The only divergence to the described procedure was that bottle data with missing temperature and/or salinity were assigned the temporal and spatial interpolated temperature and salinity derived from MIMOC (Schmidtko et al., 2013). This was done to ensure all data were in micromoles per kilogram and not requiring the discarding of already sparse data. In Schmidtko et al. (2017) this was not performed, since the error introduced near or in boundary currents and fronts can be significant. In contrast the areas here were chosen to represent homogeneous patches with significant amounts of data in the open ocean; thus in the areas analyzed here, this may only lead to minor errors in density, resulting in an error of less than 0.05 %, negligibly small in micromoles per kilogram, compared to the oxygen or nutrient data accuracy. Similarly, a nutrient file was compiled (status 25 February 2019). This compilation of hydrographic, oxygen and nutrient data is referred to in the following as “hydrodata”. The final hydrodata set was used to extract data for our regions of interest.

As nutrient data are sparse in many regions of the Pacific, areas were selected due to their better temporal data coverage (most regular sampling over the longest time period). Three regions in the equatorial Pacific were selected which can be compared with earlier observations. In the equatorial region mainly hydrodata were used. The region at 5^∘ N–5^∘ S, 165–175^∘ W (area E in Stramma et al., 2008; shown in Fig. 2a) which had hydrodata until 2009 was supplemented with data from a RV Investigator cruise at 170^∘ W from June 2016. The region 5^∘ N–5^∘ S, 105–115^∘ W (area D in Stramma et al. 2008; shown in Fig. 2a), which had hydrodata until 2008, was supplemented with data from a RV Ron Brown cruise at 110^∘ W in December 2016. The area 2–5^∘ S, 84–87^∘ W (near the Galápagos Islands, marked in Fig. 2a as area G) had been used by Czeschel et al. (2015) to investigate nutrient trends and is supplemented here with recent cruises until 2017.

Off California, intense repeated measurements have been made since 1949 within the California Cooperative Oceanic Fisheries Investigations (CalCOFI) project and data from the bottle data set (downloaded from http://www.calcofi.org/ccdata.html, last access: 13 August 2018, status 13 August 2018; data period March 1949 to November 2017, however with sparse nutrient data prior to 1984) were used to investigate the oxygen and nutrient changes within the California Current. We use the CalCOFI 1^∘ × 1^∘ data subset at 34–35^∘ N, 121–122^∘ W without additional hydrodata and no additional historic nutrient data. This region is located near the center of the CalCOFI station grid and we referred to it as CalCOFIc (shown in Fig. 2a).

Station P, located at 50^∘ N, 145^∘ W in the North Pacific, was established as a weather observation site with a weather ship in 1949 which was manned continuously until 1981, and routine hydrographic measurements were started in the 1950s. After the termination of the weather ship program, shipboard measurements have been made on average three times a year since. The data set from the Institute of Ocean Sciences, Sidney, BC, Canada, covers (status September 2018) the time period May 1956 to August 2017 with data collected on research cruises. The Station P data were supplemented here by hydrodata in the surrounding region 48–52^∘ N, 143–147^∘ W, which we refer to as area P (shown in Fig. 2a) in the following.

The Hawaii Ocean Time-series for the region north of Hawaii was taken from hydrodata in the region 22–25^∘ N, 156–159^∘ W, and downloaded bottle data for the Aloha station were taken at 22^∘45^′ N, 158^∘ W (http://hahana.soest.hawaii.edu/hot/hot-dogs/bextraction.html, last access: 15 January 2019, status 15 January 2019; time period covered October 1988 to December 2017). This data set is called the Aloha region (shown in Fig. 2a) in the following. Different to the other regions, this station is located in the subtropics with high oxygen content and low nutrient concentrations in the subsurface layer (Fig. 2). As there were continuous measurements since the 1980s we included this region in our investigation.

The measurements in the Oyashio region east of Japan (39–42^∘ N, 144–149^∘ E; shown in Fig. 2a) are from hydrodata, augmented with updated data collections used in Whitney et al. (2013). Hydrographic and biogeochemical measurements have typically been made in this region in every season since 1954 (used in Sasano et al., 2018) and were included here up to 2017.

For 50 years a repeat section was covered in the western Pacific along 137^∘ E (Oka et al., 2018). The section was measured biannually in winter and summer. To cover the western tropical Pacific, the area 20–26^∘ N, 134–140^∘ E was selected from the hydrodata and added with nutrient data for the period 2008 to 2018 from the 137^∘ E section as the hydrodata do not cover the period after 2008.

In the eastern tropical Pacific (2–5^∘ S, 84–87^∘ W; area G) and off Peru (7–12^∘ S, 78–83^∘ W; called the Peru region in the following) the hydrodata could be extended with some RV Meteor cruises carried out across the Equator and near the Peruvian coast, with two cruise legs from December 2008 to February 2009 (M77/3 and M77/4; Czeschel et al., 2011), in November 2012 (M90; Stramma et al., 2013) at 85^∘50^′ W and in June 2017 (M138). Sections across the Peruvian shelf between 9 and 16^∘ S were made during the RV Meteor cruise M91 in December 2012 (Czeschel et al., 2015) and M135 in March 2017. In October 2015 an RV Sonne cruise from Guayaquil, Ecuador, to Antofagasta, Chile, was carried out with sections at the Equator and off the shelf of Peru (Stramma et al., 2016). Although these measurements were made during one of the strongest El Niño events since the 1950s, the October 2015 measurements are used for the trend computations to estimate the influence of the El Niño event. In the figures the El Niño events of 1982/83, 1997/98 and 2015/16 are defined as very strong El Niño events, and the strong El Niño events (1957/58, 1965/66, 1972/73, 1987/88 and 1991/92) are marked by circles and the strong La Niña events (1973/74, 1975/76, 1988/89, 1998/99, 1999/2000, 2007/08 and 2010/11) are marked by squares in these years.

2.3 Data processing

On recent cruises the conductivity–temperature–depth (CTD) oxygen sensors were calibrated with oxygen measurements obtained from discrete samples from the rosette applying the classical Winkler titration method, using a non-electronic titration stand (Winkler, 1888; Hansen, 1999). The root-mean-square uncertainty of the CTD oxygen sensor calibration was on the order of ±1.0 µmol kg⁻¹.

The nutrients nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$), nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), phosphate (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$) and silicic acid (Si(OH))₄ referred to as silicate hereafter) on the recent cruises were measured on board with a QuAAtro autoanalyzer (Seal Analytical). For recent autoanalyzer measurements precisions are 0.01 µmol kg⁻¹ for phosphate, 0.1 µmol kg⁻¹ for nitrate, 0.5 µmol kg⁻¹ for silicate and 0.02 mL L⁻¹ (∼0.9 µmol kg⁻¹) for oxygen from Winkler titration (Bograd et al., 2015). For older uncorrected nutrient data, offsets are estimated to be 3.5 % for nitrate, 6.2 % for silicate and 5.1 % for phosphate (Tanhua et al., 2010).

To investigate oxygen and nutrient trends, we use the hydrodata set and recent ship sections to construct time series of annual mean oxygen and nutrient profiles. Oxygen and nutrient data are presented in micromoles per kilogram, with data obtained in different units converted to micromoles per kilogram. Linear trends and their statistical significance were computed using binned annual data of the linearly standard depth interpolated profiles (all measurements from one year were attributed to that year, as in earlier investigations; e.g., Czeschel et al., 2015). This was done since intra-annual variability is small compared to the inter-annual variability and the modes analyzed are of greater time spans; the reduction of data clusters greatly reduces the possible systematic bias in the trends. We do not address the errors in the vertical interpolation on standard depths, because the layer analyzed is thin and characterized in general by a linear decrease or increase in observed values. Due to heterogeneous vertical sampling the error is assumed to be small compared to the observed inter-annual variability. For the trend computation, successive years of the time series are not necessarily statistically independent of each other. Therefore we determined the effective number of degrees of freedom for the computation of the confidence interval. For each time series the lag-dependent temporal half folding range of the autocorrelation function was used to compute the degrees of freedom. The length of the time series was divided by the length of the time derived from the autocorrelation function, giving a statistical measure of the degrees of freedom for each analyzed time series independently.

Table 1 Linear trends of solutes in micromoles per kilogram per year with 95 % confidence intervals (p values) where data are available for the entire period since 1950, for negative (1950–1976; PDO-) and positive (after 1976; PDO+) PDO periods in the 50 to 300 m depth layer for selected ocean areas shown in Figs. 3, 4, 5, 7, S3, S4, S5 and S6. Trends whose 95 % confidence interval includes zero are shown in italics. Areas are named CalCOFIc (34–35^∘ N, 121–122^∘ W), area P (48–52^∘ N, 143–147^∘ W), Aloha region (22–25^∘ N, 156–159^∘ W), Oyashio region (39–42^∘ N, 144–149^∘ E), 137^∘ E (20–26^∘ N, 134–140^∘ E) and Peru region (7–12^∘ S, 78–83^∘ W).

As historical measurements are focused on the upper ocean and as oxygen and nutrient changes will have the largest impact on the biology in the upper ocean, the trends have been computed for the subsurface layer 50 to 300 m as presented in Czeschel et al. (2015). The upper boundary at 50 m was selected to avoid influence from atmosphere–ocean interaction in the mixed layer. In the North Pacific seasonal variability will reach below 50 m depth. The two areas we consider in the North Pacific contain a lot of data, and a test omitting the winter month measurements (January to March: Supplement Table S1) shows similar trends as for the entire year (Table 1). Hence the seasonal cycle did not have a larger impact on the results. The oxygen and nutrient distribution at 50 to 300 m depth (Fig. 2) shows the large variation in the parameters across the Pacific. The subtropical gyres are clearly visible by enhanced oxygen and low nutrient content. Higher nutrient contents are seen in the equatorial and eastern Pacific in the areas of equatorial and coastal upwelling, respectively. As gradients are low in regions of low nutrient content and nutrient data are often sparse, it makes sense to investigate changes in regions with enhanced nutrient content and a sufficiently long time series.

The data used for the oxygen and nutrient time series were interpolated with an objective mapping scheme (Bretherton et al., 1976) with Gaussian weighting using a temporal half folding range of 0.5 years and a vertical one of 50 m, a maximum temporal range of 1 year, and spatial range of 100 m. The covariance matrix was computed from the nearest 100 local data points and 50 randomly distributed data points within the maximum range, for the diagonal of the covariance matrix a signal-to-noise ratio of 0.7 was set (see Schmidtko et al., 2013, for details). Due to the random data points, the computed trends and confidence intervals vary slightly in each computation run of the interpolation; however the variation is very small compared to the confidence interval.

The correlation coefficient and significance between the oxygen time series and the PDO or NPGO were computed for all years available, using a chi-squared test with the degrees of freedom as computed above. Often the correlation was considered to be significant for the 95 % confidence interval different from null; however, recently it was stated that it seldom makes sense to calibrate evidence as a function of p values (McShane et al., 2019), and it was suggested to retire the term statistical significance and only state direct p values (Amrhein et al., 2019). In the following we list the p value of the null hypotheses in parentheses. For some regions with a longer time series the correlation was also computed for temperature or nitrate. In the North Pacific the correlation of oxygen with an 18.6-year oscillating signal was computed for the entire period with existing oxygen data. To investigate the lag the measurement years were shifted from −15 to +15 by 1-year steps. Since the PDO or NPGO time series are continuous, data point reduction due to lag shift is small and only occurring for individual data collected in the last 15 years. The impact is assumed to be smaller than the here-given uncertainties.

3 Trends and influence of climate signals

3.1 Long-term trends

In the global ocean the long-term SST trend in 1901–2012 was positive everywhere except for a region in the North Atlantic (IPCC, 2013, Fig. 2.21). However, for 1951 to 1980 decreasing SSTs were reported across the North Pacific Ocean (IPCC, 2013, Fig. 2.22). For 1981 to 2012, while the western Pacific showed a warming trend, a large region with decreasing SSTs was seen in the eastern Pacific Ocean. However, a different pattern emerges when the analysis is applied to the subsurface layer (50 to 300 m). Subsurface temperature trends computed for all data since 1950 in each of the regions discussed in Sect. 2 showed weak temperature increases. Exceptions are (1) the area P with a temperature increase in the 95 % confidence interval of 0.0083±0.0073 ^∘C yr⁻¹ for the period 1954 to 2017, with the highest temperatures in the positive PDO periods from 1977 to 1999 and since 2013, and (2) the Oyashio region, where a strong temperature decrease in the 95 % confidence interval of $-\mathrm{0.0273}\pm \mathrm{0.0188}$ ^∘C yr⁻¹ was derived for the period 1952 to 2017 with the lowest temperatures in the period 1977 to 2010 (Fig. S2). This agrees with the surface layer (0 to 50 m), where all areas showed increasing temperature for the entire measurement period, except for the Oyashio region where temperature decreased (Fig. S2). The very strong temperature increase at area P after 2013 was impacted by the strong surface temperature anomaly in the Northeast Pacific Ocean from 2013 to 2015 during the marine heat wave nicknamed “The Blob” (Bond et al., 2015; Di Lorenzo and Mantua, 2016). For the subsurface layer discussed here, however, it appears that the marine heat wave peaked in 2016–2017 (Jackson et al., 2018). Except for the Aloha and the CalCOFIc regions, all 0 to 50 m trends were not within the 95 % confidence interval, probably due to interdecadal, seasonal and regional variations in the temperature measurements. The temperature trends are in agreement with the SST expression of the PDO (Fig. S1), which is positive in the area P region and strongly negative in the Oyashio region.

Figure 3 Annual mean oxygen concentration (µmol kg⁻¹, red crosses) for years available and trends for the layer from 50 to 300 m plotted for the entire time period (dashed red lines) and for the periods 1950 to 1976 for the negative PDO phase and after 1976 for the positive PDO phase (solid red lines) for (a) 5^∘ S–5^∘ N, 165–175^∘ W (area E) and (b) 5^∘ S–5^∘ N, 105–115^∘ W (area D). El Niño years defined as strong or very strong are marked by an additional magenta circle, strong La Niña years are marked by an additional blue square. The changes of the PDO status in 1977, 1999 and 2013 are marked by vertical dotted lines. The annual mean PDO index is shown as a grey curve.

In the open ocean a decline in mean oxygen solubility of ∼5 µM associated with a hypothetical warming of 1 ^∘C throughout the upper ocean would expand the reach of hypoxic conditions by 10 % while suboxic zones would nearly triple in volume (Deutsch et al., 2011). However the sensitivity of hypoxic zones to variations in the depth of the thermocline introduces a mechanism for counteracting this expansion (Deutsch et al., 2011). In the areas of the equatorial Pacific (Figs. 3 and 4), the CalCOFIc region off California (Fig. 5), the area P in the North Pacific (Fig. S3), the Aloha region north of Hawaii (Fig. S4), the Oyashio region in the western North Pacific (Fig. 7) and a region off Peru (Fig. S6), the linear trend of the oxygen content of the layer from 50 to 300 m is negative for the entire time period of available data since 1950, except for the 137^∘ E data between 134 and 140^∘ E in the western Pacific (Fig. S5) and the eastern equatorial region 2 to 5^∘ S, 84 to 87^∘ W (just northwest of the Peru region) which has a positive trend caused by low and variable oxygen content in the 1955 to 1965 period. However, except for the western equatorial areas and area P, the oxygen trends are not within the 95 % confidence interval with regard to the entire time series (Table 1).

The long-term trends of the western equatorial regions between 5^∘ N and 5^∘ S at 165–175^∘ W and 105–115^∘ W show a continuous oxygen decrease since 1950 for the 50 to 300 m layer of $-\mathrm{0.36}\pm \mathrm{0.22}$ and $-\mathrm{0.65}\pm \mathrm{0.37}$ µmol kg⁻¹ yr⁻¹ (Fig. 3; Table 1). These trends are larger than the trends for the same regions since 1960 for the 300 to 700 m layer of $-\mathrm{0.19}\pm \mathrm{0.20}$ and $-\mathrm{0.13}\pm \mathrm{0.32}$ µmol kg⁻¹ yr⁻¹ (Stramma et al., 2008), as the deeper layer is located at the low-oxygen core of the OMZ where no large changes are possible.

For the period 1956 to 2006 at area P in the depth range 100 to 400 m an ocean warming of 0.005–0.012 ^∘C yr⁻¹ has been described (Whitney et al., 2007). Our longer temperature time series 1954 to 2017 at 50 to 300 m of 0.0083 ^∘C yr⁻¹ confirms this trend (Fig. S2). The oxygen trend of $-\mathrm{0.24}\pm \mathrm{0.23}$ µmol kg⁻¹ yr⁻¹ in area P (48–52^∘ N, 143–147^∘ W) for the 50 to 300 m layer for 1954 to 2017 is smaller than that previously reported for area P at different depth layers between 100 and 400 m for the shorter time period 1956 to 2006 (0.39–0.70 µmol kg⁻¹ yr⁻¹; Whitney et al., 2007). For yet another time period (1987 to 2011) and the layer 100–500 m in area P, trends have been reported (Whitney et al., 2013): −0.9 µmol L⁻¹ yr⁻¹ for oxygen (density × µ mol kg⁻¹ yr⁻¹; ∼0.88 µmol kg⁻¹ yr⁻¹), 0.085 µmol L⁻¹ yr⁻¹ for nitrate, 0.30 µmol L⁻¹ yr⁻¹ for silicate and 0.0033 µmol L⁻¹ yr⁻¹ for phosphate. Our trends for 50 to 300 m for the longer time period 1977 to 2017 (Table 1) are much smaller for oxygen (−0.18 µmol kg⁻¹ yr⁻¹) and slightly smaller for the nitrate (0.093 µmol kg⁻¹ yr⁻¹), silicate (0.193 µmol kg⁻¹ yr⁻¹) and phosphate (0.001 µmol kg⁻¹ yr⁻¹).

Figure 4 Annual mean concentration (µmol kg⁻¹, red crosses) for years available (see Table 1) and trends for the layer from 50 to 300 m plotted for the entire time period (dashed red lines) and for the periods 1950 to 1976 for the negative PDO phase and after 1976 for the positive PDO phase (solid red lines) between 2 and 5^∘ S, 84 and 87^∘ W (area G) µmol kg⁻¹ yr⁻¹ for oxygen, nitrate, silicate and phosphate. El Niño years defined as strong or very strong are marked by an additional magenta circle and strong La Niña years by an additional blue square. For oxygen, measurements in 1982 were removed as the 50–300 m mean was much too high (104.8 µmol kg⁻¹), and for nitrate, measurements in 1967 which were too high (36.8 µmol kg⁻¹) were removed. The changes of the PDO status in 1977, 1999 and 2013 are marked by vertical dotted lines. The annual mean PDO index is shown in the oxygen time series as a grey curve and the NPGO index is shown in the nitrate, silicate and phosphate time series as a grey curve.

Figure 5 Annual mean concentration (µmol kg⁻¹, red crosses) for years available and trends for the layer from 50 to 300 m plotted for the entire time period (dashed red lines) and for the periods 1950 to 1976 for the negative PDO phase and after 1976 for the positive PDO phase (solid red lines) between 34 and 35^∘ N, 121 and 122^∘ W from the CalCOFIc bottle data in micromoles per kilogram per year for (a) oxygen, (b) nitrate, (c) silicate and (d) phosphate. El Niño years defined as strong and very strong are marked by an additional magenta circle and strong La Niña years by an additional blue square. The changes of the PDO status in 1977, 1999 and 2013 are marked by vertical dotted lines. The annual mean PDO index is shown in the oxygen time series as a grey curve and the NPGO index is shown in the nitrate, silicate and phosphate time series as a grey curve.

Figure 6 Correlation between oxygen at 50–300 m in the region 34–35^∘ N, 121–122^∘ W (CalCOFIc) and (a) the PDO and (b) the NPGO (black lines) shifted between −15 and +15 years for the years after 1976 and the p values (x) which often are used to declare non-significant correlation for p∼ > 0.05 (blue lines for 0 and 0.05). For positive years PDO–NPGO lags; for negative years they lead.

The Aloha region is located in the southern part of the North Pacific subtropical gyre where oxygen is high and the nutrient inventory low. The oxygen trend since 1951 is negative, although not in the 95 % confidence interval. Nitrate and silicate trends (available since 1984/1985) are positive, while phosphate decreases. However, only the trends in silicate (including one measurement in 1970) and phosphate are within the 95 % confidence interval (Table 1).

For a region similar to our Oyashio region, an oxygen decrease of 0.73 µmol kg⁻¹ yr⁻¹ between 1968 and 1998 has been reported for the density layer 26.8–27.4 kg m⁻³ (∼260–1030 m) (Watanabe et al., 2003). Our oxygen trend for the layer 50 to 300 m is slightly positive (+0.15 µmol kg⁻¹ yr⁻¹) for the period 1977–2017 and negative (−0.23 µmol kg⁻¹ yr⁻¹) for the longer time period of 1952 to 2017. The phosphate increase in the deep layer 26.8–27.4 kg m⁻³ was reported as 0.004 µmol kg⁻¹ yr⁻¹ in Watanabe et al. (2003) while in our 50 to 300 m layer the phosphate trends for the positive and negative PDO periods and the entire period are larger, up to 0.010 µmol kg⁻¹ yr⁻¹ (Table 1). We observe a large silicate increase in the Oyashio region, +0.667 µmol kg⁻¹ yr⁻¹ since 1981 (Table 1). This exceptional silicate enrichment in the Oyashio region was noticed before and it was speculated that warming and a reduction of dense water formation causes this enrichment, because less silicate is transported into deep (> 300 m) waters (Whitney et al., 2013). Despite the different depth layers, the comparison with earlier investigations indicates that changes are related to time periods and not just a linear trend.

Figure 7 Annual mean concentration (µmol kg⁻¹, red crosses) for years available and trends for the layer 50 to 300 m plotted for the periods 1950 to 1976 for the negative PDO phase and after 1976 for the positive PDO phase (solid red lines) and for the entire time period (dashed red lines) for the Oyashio region (39–42^∘ N, 144–149^∘ E) from hydrodata CTD and bottle data and a data collection used in Whitney et al. (2013) and Sasano et al. (2018) in micromoles per kilogram per year for (a) oxygen, (b) nitrate, (c) silicate and (d) phosphate. For nitrate, measurements in 1963 were removed as the 50–300 m mean was much too high (1.27 µmol kg⁻¹). El Niño years defined as strong are marked by an additional magenta circle and strong La Niña years by an additional blue square. The changes of the PDO status in 1977, 1999 and 2013 are marked by vertical dotted lines. In the oxygen time series the 18.6-year sinusoidal nodal cycle is included (green curve). The annual mean PDO index is shown in the oxygen time series as a grey curve, and the annual mean NPGO index is shown in the nitrate, silicate and phosphate time series as a grey curve.

In the western Pacific at 137^∘ E (20–26^∘ N) an opposite long-term trend compared to the Oyashio region was observed (Table 1). In the 50–300 m layer oxygen increased while nitrate, silicate and phosphate decreased (Fig. S5). According to a map of oxygen changes in the upper 1200 m since 1960 (Oschlies et al., 2018; Fig. 3a) this is a region with oxygen trends in the range −5 to +5 µmol kg⁻¹ per decade; hence the weak positive oxygen trend is possible in this region. Also, the weak positive temperature trend (not shown) at 137^∘ E is opposite to the negative temperature trend in the Oyashio region.

In the Peru region the oxygen, silicate and phosphate in the 50 to 300 m layer decrease while nitrate increases (Fig. S6). However, these long-term trends are not within the 95 % confidence interval (Table 1), possibly due to the paucity of data and the reversal of trends related to the PDO phases as described below.

The oxygen time series from most Pacific regions show decreasing oxygen, although with varying magnitudes. These variations could be related to different climate signals as investigated below. The nutrient time series often show an increase over time: as expected, this is the opposite trend to oxygen.

Table 2 Correlation coefficient and p value for available data since 1950 between annual layer from 50 to 300 m concentration and PDO and NPGO with PDO–NPGO lags (negative PDO–NPGO leads). Areas with large p values > 0.05 are considered not significantly correlated.

3.2 The influence of the Pacific Decadal Oscillation (PDO)

An influence of the PDO has been seen in oxygen measurements (e.g., Czeschel et al., 2012) and modeling studies (e.g., Frölicher et al., 2009). According to the description of the PDO influence on the thermocline depth (e.g., Deutsch et al., 2011; Chavez et al., 2003) it is expected that during cold PDO phases the oxygen will decrease and the nutrients will increase in the eastern equatorial and tropical Pacific, while during warm PDO periods the oxygen should increase and the nutrients decrease. However, visual inspection of the oxygen time series in the equatorial Pacific (areas E and D; see Fig. 2 in Stramma et al., 2008) indicates stagnant oxygen concentrations before 1976 during the cold PDO and enhanced oxygen depletion since the 1980s in the OMZ in the subsurface layer. The annual mean oxygen concentration for the layer 50 to 300 m for 5^∘ S–5^∘ N, 165–175^∘ W (area E) and 5^∘ S–5^∘ N, 105–115^∘ W (area D) shows a strong oxygen decrease after 1976 (Fig. 3a and b). As nutrient data in the open Pacific are sparse, no reliable nutrient trends could be derived for areas E and D. The correlation of the PDO with the 50 to 300 m oxygen is high, with p values of ∼0.001 for all areas along the Equator with the highest correlation coefficient of +0.83 at 2 to 5^∘ S, 84 to 87^∘ W (area G) and slightly decreasing towards the west (Table 2).

Czeschel et al. (2015) showed decreasing oxygen and increasing nutrients in the 50 to 300 m layer of the area 2–5^∘ S, 84–87^∘ W (area G) since 1976. If PDO influence acts as expected, an even stronger gradient should exist prior to 1976. For the layer 50 to 300 m between 2 and 5^∘ S, 84 and 87^∘ W, for the post-1976 positive PDO period, the computed trend is slightly modified compared to the time period investigated in Czeschel et al. (2015). Here we see decreasing oxygen concentrations but increasing nitrate and phosphate and decreasing silicate concentrations since 1977 (Fig. 4). The main differences between these time series are improvements to the objective analysis, less smoothing being applied, and two more cruises added for 2015 and 2017. The resulting trends since 1977 for oxygen and nitrate are smaller, for phosphate larger and for silicate reversed compared to Czeschel et al. (2015). Oxygen for the negative PDO phase (1955 to 1976) showed a strong positive trend of 1.63±1.18 µmol kg⁻¹ yr⁻¹ for the 50 to 300 m depth layer in the area G which led to a small positive oxygen trend for the entire time period since 1955.

Near California the CalCOFIc bottle data in the region 34–35^∘ N, 121–122^∘ W (Fig. 5, Table 1) show trends similar to those for the equatorial area (2–5^∘ S, 84–87^∘ W; area G) for the period since 1977, with decreasing oxygen and increasing nitrate, phosphate and silicate. The long-term trend over all measurements since 1950 shows increasing nitrate and phosphate and decreasing oxygen and silicate for the CalCOFIc region similar to the long-term trends for the equatorial region. The correlation between the 50–300 m annual means and the PDO annual mean is 0.42 (p=0.002) for oxygen (Fig. 6a; Table 2) and −0.37 (p=0.019) for nitrate (Table 2) for the period after 1976 with the PDO leading by 11 years for oxygen.

For the area P (48–52^∘ N, 143–147^∘ W) the correlation with PDO is −0.30 (p=0.019) (Table 2) for oxygen, and the correlation coefficients are −0.27 (p=0.109) for nitrate and +0.46 (p < 0.001) for temperature (with the PDO lagging by 10 years for oxygen). It is remarkable that the trends in oxygen are unchanged, whether fitting to the entire record or to the positive and negative PDO time periods separately. Hence in this northern Pacific region the PDO has a weak influence on the 50 to 300 m biogeochemistry, probably caused by water masses propagating by 5 to 15 years from the Oyashio region into this part of the North Pacific (Ueno and Yasuda, 2003) while the correlation of the PDO and oxygen is larger in the other regions of the North Pacific.

The Aloha region is located at the transition of the warm and cold area of the PDO; thus one might expect the influence of the PDO to be weak in this region. However, the oxygen was observed to increase during cold PDO and the trend was close to zero in the warm PDO phase (Fig. S4). During the PDO warm phase, nitrate and silicate increased while phosphate decreased. There are too few nutrient data in the period of the PDO cold phase and only phosphate is available with a minor increase (Table 1).

In the global PDO distribution (Fig. S1, top) the Oyashio region in the western North Pacific is located in a reversed temperature anomaly pattern compared to the eastern Pacific, so one might expect to see the opposite trends. The oxygen decreases during the cold PDO phase and slightly increases during the warm PDO phase. The nutrients increased during both PDO phases except for silicate, which decreased during the cold PDO phase (Fig. 7, Table 1). Different to the other regions, the temperature of the 50 to 300 m layer shows a decrease in this region (Fig. S2). The oxygen and nitrate of the Oyashio region show strong correlations with the PDO (Table 2) with the PDO lagging by 4 to 5 years, indicating a delay between the surface signal and the changes in the subsurface layer. In a previous study where oxygen concentrations were investigated from 1954 to 2014, a decrease in oxygen was attributed to a reduction of ventilation in winter due to warming and freshening and reduction of dense water formation in the Sea of Okhotsk (Sasano et al., 2018).

In the 137^∘ E region, except for nitrate with a negative trend during the cold PDO phase and a positive trend during the warm PDO phase, oxygen increases and silicate and phosphate decrease during both phases (Fig. S5). This indicates that there should be no strong correlation to the PDO phase; however the correlation reaches 0.45 (p < 0.001). According to the PDO expression on SST (Fig. S1) the PDO signal should be weak in the southwestern North Pacific. For the period 1985 to 2010 Takatani et al. (2012) described an oxygen decrease in the 137^∘ E section; however during this time period our oxygen measurements are also quite similar and the positive oxygen trend is mainly caused by low oxygen values before 1970 and high oxygen values since 2010.

In the Peru region there is an increase in oxygen in the 50 to 300 m layer during the cold PDO phase and a decrease during the warm PDO phase. Nitrate, silicate and phosphate show trends opposite to oxygen. Despite the paucity of data in this region, these PDO-related trends are within the 95 % confidence limit for oxygen in the cold PDO phase and for nitrate and phosphate since 1977 (Table 1). The correlation of the PDO with the 50 to 350 m oxygen is +0.45 (p=0.012). This indicates that the observed changes in oxygen and nutrients off Peru are associated with the PDO.

Although the time period since the shift from negative to positive PDO in 2013 is short, most areas examined here show higher 50 to 300 m oxygen concentrations than the trend line for the period since 1978 and lower nitrate and silicate concentrations than the trend line since 1977, except for the Aloha and 137^∘ E regions, the oxygen in 2017 in the Peru region, the El Niño year 2015/16 and The Blob in area P.

Figure 1 shows a global mean surface temperature increase before 1945, a stagnant global mean surface temperature trend during the PDO cool phase between 1945 and 1976, and a global mean surface temperature increase after 1976 despite this period encompassing a PDO warm and cold phase. As the influence of the warm PDO phase from 1977 to 1999 and the cold phase from 1999 to 2014 is not related to major oxygen and nutrient trend changes, the increasing temperature seems to be a major component of setting the long-term oxygen trend in the Pacific Ocean.

3.3 The influence of the North Pacific Gyre Oscillation (NPGO)

As the NPGO is defined for a smaller region in the northeast Pacific than the PDO, its largest influence on oxygen and nutrients is expected to be in the North Pacific. The NPGO index shows higher variability than the PDO index (e. g. Fig. 4). Strong NPGO minima were present in the years 1967, 1980, 1994, 2006 and 2015 and maxima were present in the years 1961, 1976, 1988, 2000 and 2010.

In the area P, the correlation with NPGO is −0.28 (p=0.029) for oxygen in the 50 to 300 m layer for the period since 1977 with NPGO leading by 4 years and +0.58 (p < 0.001) for nitrate with NPGO lagging by 1 year (Table 2). The oxygen data are highly variable, which might have led to the low correlation with the NPGO; however the nitrate correlation as well as a correlation of NPGO with temperature of −0.40 (p=0.001) show a strong relationship with the NPGO in this region.

Since 1980, the maxima and minima of the CalCOFIc time series of nitrate, phosphate and silicate (Fig. 5) often agree with the NPGO maxima and minima. The correlation between the 50–300 m annual mean and the NPGO annual mean is −0.35 (p=0.006) for oxygen (Fig. 6b) and +0.35 (p=0.025) for nitrate with the NPGO lagging by 1 year. The correlation with the NPGO is slightly weaker than with the PDO at CalCOFIc. This confirms the described decadal variations linked to the NPGO in the North Pacific and the California Current System (Di Lorenzo et al., 2009).

In the Oyashio region, oxygen shows a correlation of +0.35 (p=0.005) with the NPGO, and the correlation of nitrate and temperature with the NPGO is of similar strength. In the western part of the North Pacific at 137^∘ E the correlation with the NPGO is weak (−0.19; p=0.163). At the southern part of the North Pacific subtropical gyre, in the Aloha area, the NPGO correlation coefficient of +0.39 (p=0.008) with oxygen is similar to the PDO correlation. In the central equatorial Pacific and off Peru the correlation of oxygen with NPGO is lower than the correlation with the PDO.

3.4 The influence of the North Pacific Index (NPI) and El Niño–Southern Oscillation (ENSO)

The oscillation every 2 decades related to NPI with minima in 1962 and 1983 and maxima in 1971 and 1991 (Watanabe et al., 2003, their Fig. 2) is difficult to see in our analysis of the Oyashio region due to large year-to-year variability (Fig. 7). The correlation of the NPI (November to March anomaly) with the 50 to 300 m oxygen time series in the North Pacific leads to correlations of 0.38 (p=0.002) in the Oyashio region, 0.37 (p=0.002) in area P, −0.42 (p=0.002) at CalCOFIc, 0.33 (p=0.025) at Aloha and −0.29 (p=0.029) at 137^∘ E.

An oscillation of 16.4–19.6 years in oxygen every 2 decades, possibly driven by nodal tidal cycles of 18.6 years, was described recently for the Oyashio region with maxima at about 1971, 1989 and 2008 and minima at about 1962, 1980 and 1998 (Fig. 3 in Sasano et al., 2018). The 50 to 300 m oxygen time series does not show a strong visual correlation for the period 1961 to 2008 between the trend-corrected oxygen and an 18.6-year oscillation. In area P, oxygen trends include periods of increased ventilation of deeper isopycnals on a ∼18-year cycle (Whitney et al., 2007). The correlation of the oxygen content of the 50 to 300 m layer with an 18.6-year cycle for the Oyashio region was 0.14 (p=0.244). The correlation with the 18.6 year oscillation is 0.084 (p=0.51) in area P, −0.139 (p=0.28) at CalCOFIc, 0.146 (p=0.339) at Aloha and −0.142 (p=0.287) at 137^∘ E. Although the 18.6-year tidal cycle could play a role as a basic forcing for the bi-decadal (every 2 decades) ocean variations (Yasuda et al., 2006), the correlation in the North Pacific is much stronger with the NPI than with the 18.6-year oscillation, likely due to short-term fluctuations in the NPI and the observed data or the sensitivity to phase shifts in the oscillation.

Parameter distributions similar to those for cold PDO periods exist for La Niña events and for warm PDO periods for El Niño events (Deutsch et al., 2011). Surprisingly, in the equatorial regions (Figs. 3 and 4) the subsurface oxygen concentration at 50 to 300 m depth shows no clear anomalies in years of ENSO events. The ENSO signal seems to be restricted to the near-surface layer in the equatorial Pacific. In the eastern Pacific during El Niño periods, oxygen in the upper ocean is higher and nutrients are lower in the upwelling regions (CalCOFIc and Peru region) due to either reduced upwelling or upwelling of oxygen-richer and nutrient-poorer water masses. In the CalCOFIc region, the measurements during the very strong El Niño events in 1997/98 and 2015/16 show higher oxygen concentration and very low nutrient concentration in the 50 to 300 m layer when compared to the trend line and the neighboring years (Fig. 5). The deviations are very strong for the 1997/98 El Niño while moderate for the 2015/2016 El Niño. This signal is also visible for the 2015/2016 El Niño in the Peru region (Fig. S6) and for the 1997/98 El Niño at a shelf station off Lima (Graco et al., 2017). Not all strong El Niño events are associated with similar anomalies. Offshore from the upwelling region, in area P, the anomalies for the very strong El Niños 1997/98 and 2015/2016 are opposite, with low oxygen and high nutrient concentrations during the earlier event (Fig. S3). For La Niña events a reversed trend is visible in the eastern Pacific for some events, e.g., with low oxygen and high nutrient concentrations for the 1988/1989, the 1998/1999 and the 1999/2000 La Niña events in the CalCOFIc region (Fig. 5). These ENSO-related signals disappear in the following year, and hence the ENSO-related changes in oxygen and nutrients do not show a multi-year signal.

4 Discussion

One might wonder if the observed changes in trends in oxygen and nutrients are influenced by evolving methods rather than changes in climate. While the Winkler titration method to measure oxygen has remained the same since the early 1900s, the methods to determine nutrients have varied over time. However, except for the very low nutrient data in the Aloha region in the center of the North Pacific, all areas analyzed here show a similar and relatively large range for nitrate, silicate and phosphate. Precision of 5 % for nitrate measurements of more than 10 µmol L⁻¹, of about 6 % or less for silicate, and 5 % for phosphate at 0.9 µmol L⁻¹ and larger are reported for early nutrient measurements (Hansen and Koroleff, 1999). Similar offsets for measurements after the 1990s were derived with 3.5 % for nitrate, 6.2 % for silicate and 5.1 % for phosphate (Tanhua et al., 2010); accordingly the offsets for the eastern equatorial region (Fig. 4) could be as high as ∼1 µmol kg⁻¹ for nitrate, ∼0.1 µmol kg⁻¹ for phosphate and ∼1.0 µmol kg⁻¹ for silicate and hence smaller than the observed long-term trends.

To put the regions examined here in context, we compare with previously published trends for the entire North Pacific as well as changes with impact on the ocean circulation. Yasunaka et al. (2016) reported trends of surface phosphate and silicate averaged over the North Pacific from 1961 to 2012 as $-\mathrm{0.012}\pm \mathrm{0.005}$ and $-\mathrm{0.38}\pm \mathrm{0.13}$ µmol L⁻¹ per decade, respectively, whereas the nitrate trend averaged over the North Pacific was 0.01±0.13 µmol L⁻¹ per decade. This is in contrast to the subsurface layer examined here, where nitrate tended to increase over a similar time period. In particular, high nitrate increase was observed for the Oyashio region with +0.143 µmol kg⁻¹ yr⁻¹ for the period 1977 to 2017 (Table 1). An increase in anthropogenic nitrogen emissions from northeastern Asia and subsequent deposition over the North Pacific resulted in a detectable increase in nitrate concentrations in the near-surface layer since the 1970s (Kim et al., 2014). The observed nitrate increase at 50 to 300 m might be the response of decreased water subduction. For example, the North Pacific Intermediate Water is a dominant pathway to enter the mid-depth waters in the North Pacific and was freshening in the period 1960 to 1990 (Wong et al., 2001). Since the overturning in the North Pacific originates from the Sea of Okhotsk through dense shelf water, the observed freshening to a depth of ∼500 m during the past 4 decades could possibly weaken the shallow overturning of the North Pacific (Ohshima et al., 2014). Thus overturning ocean circulation in the North Pacific is likely slowing because less dense water is forming in the Sea of Okhotsk (Sasano et al., 2015, 2018) and hence decreasing the supply of oxygen.

Model results indicate that more than 50 % of the total internal variability of oxygen is linked to the PDO in the North Pacific surface and subsurface waters (Frölicher et al., 2009). The long-term 50–300 m trends since the 1950s in the eastern equatorial Pacific and the CalCOFIc region (Figs. 4 and 5; Table 1) indicate a long-term increase in nitrate and phosphate and a decrease in silicate, but often with reversed trends in oxygen and nutrients when separated into cold and warm PDO phases. From the 1980s to the 2010s in the North Pacific, oxygen decreased while nitrate, phosphate and silicate increased (Whitney et al., 2013), similar to what was observed in the eastern equatorial Pacific and the CalCOFIc region for the period after 1976 (Figs. 4 and 5). For the California Current System the decadal oxygen changes seem to be primarily controlled by ocean circulation dynamics (Pozo Buil and Di Lorenzo, 2017). Pozo Buil and Di Lorenzo (2017) found that subsurface anomalies in the core of the North Pacific Current propagate the oxygen signal downstream within about 10 years to the coastal regions and predicted a strong decline in oxygen by 2020 in the California Current System. The recent measurements in the CalCOFIc region, shown here, support this prediction.

The nutrient increase in the CalCOFIc region since the 1980s could be related to upwelling variability in the California Current. A strong nitrate flux from 1980 to 2010 was driven almost entirely by enhanced equatorward winds, negating a weak negative trend associated with increased surface heat flux (Jacox et al., 2015). However, changes in the properties of source waters (primarily from the eastern tropical Pacific via the California Undercurrent) have likely driven most of the biogeochemical trends observed in the southern California Current (Meinvielle and Johnson, 2013; Bograd et al., 2015; Nam et al., 2015).

Despite the low data coverage, the measurements in the Peru region confirm the expected opposite trends for oxygen and nutrients related to the PDO phases. For a shallow shelf station (145 m depth) near Lima, measurements in the upper 100 m show increasing oxygen concentrations for the period 1999 to 2011 (Graco et al., 2017). This contrasts with the decreasing oxygen trend we observe from 1977 to 2017 in the Peru region (Fig. S6, Table 1) and indicates that different processes and trends might exist on the shelf compared to the open ocean.

In the eastern tropical and subtropical Pacific very strong El Niño and some strong La Niña events are apparent in the oxygen and nutrient distribution but do not result in a multi-year signal or trend.

5 Concluding remarks

In this study, we investigated the influence of well-documented atmosphere–ocean decadal oscillations on the trends in oxygen and nutrients in the upper subsurface layer of the Pacific Ocean. Due to the limited subsurface nutrient data, only select areas were investigated and the results may have larger uncertainties for the areas with low data coverage and in particular for the combination of the warm and cold PDO periods after 1976. Especially in the South Pacific data are sparse, in part because not all existing data have been made public. Statistically significant trends and correlations hold true for the period and for data analysis, though in data-sparse regions, these findings should be the subject of future scrutiny. Regional- and small-scale processes, related to or independent of large-scale PDO forcing, may alter the signals seen so far.

A test excluding the winter months (January to March) in the area P and the Oyashio region (Table S1) showed that the seasonal cycle had little influence on the trends derived for the 50 to 300 m layer. The depth layer of 50 to 300 m was selected as this is the major layer of biological subsurface activity influenced by oxygen and nutrient variability. For example, several warm-water mesopelagic species in the California Current, which are apparently adapted to the shallower, more intense OMZ off Baja California, were shown to be increasing despite declining midwater oxygen concentrations and becoming increasingly dominant (Koslow et al., 2019). Enhanced biological activity in coastal regions might lead to larger nutrient variability and obscure climate-related signals. One also has to keep in mind that the results might be influenced to a degree by changes in the propagation of water masses, as the ocean dynamics are not stationary and nor are trends in adjacent water masses.

Agreeing well with the regions with the largest SST signal of the PDO (Fig. S1a), the PDO seems to have the strongest influence in the 50 to 300 m layer in the equatorial and eastern Pacific (Table 2). During the cold PDO phase in the eastern Pacific and the stagnant global mean surface temperature signal, in the CalCOFIc region and the Peru region the 50 to 300 m oxygen increases and the nutrient concentrations decrease (Figs. 5, S6) and show the opposite trends during the global warming period since 1977, which we call the warm PDO phase despite a period of a PDO cool phase (Fig. 1). An increase in oxygen from 1950 to 1980 and a decrease after 1980 have also been described for some isopycnal surfaces in subsurface waters of the northeast Pacific (Crawford and Peña, 2016). In the western Pacific at 137^∘ E the influence of the PDO should be weak and even a long-term oxygen increase was observed despite a large correlation coefficient with the PDO.

With respect to other climate indices, the results are more mixed and statistical significance should be treated with caution. The NPGO has the largest impact on decadal variations in the North Pacific in area P (Table 2). The NPGO influence is also visible in the central and eastern North Pacific and the equatorial Pacific, although weaker than the correlation with the PDO. The NPI is well correlated to the oxygen changes in the North Pacific with the largest correlation at CalCOFIc. The 18.6-year nodal tidal cycle has the largest correlation in the Oyashio region, the area P and CalCOFIc; however the correlation is weak compared to the correlation with the NPI. In the Oyashio region a combination of the PDO, the NPGO, the NPI and the 18.6-year nodal cycle contributes to the trends in temperature, oxygen and nitrate. The oxygen decline in the Oyashio region has different controlling mechanisms for different depth layers. In the upper layer above σ_θ=26.7 kg m⁻³ (∼200 m) the oxygen decline is primarily attributed to the reduction of winter convection upstream and in part to the deepening of isopycnal surfaces due to warming and freshening in the upper layers, while below σ_θ=26.7 kg m⁻³ the oxygen decline is attributed to the reduction of dense shelf water formation in the Sea of Okhotsk associated with the reduction of sea ice production and freshening (Sasano et al., 2018).

The PDO and NPGO strongly influence the trends in oxygen and nutrient inventories; nevertheless the long-term global mean surface temperature trend (Fig. 1) seems to play a role in the oxygen trend as long-term trends indicate an oxygen decrease throughout the cold and warm PDO phases in all but two areas. Note that while the most likely basin-wide drivers of oxygen and nutrient variability were investigated here, other contributors might exist depending on the region, as shown in detail for the Oyashio region by Sasano et al. (2018).

As greenhouse gas concentrations rise further, a variable but positive trend of increasing global mean surface temperature is expected, which should lead to a continuing decrease in oxygen and increase in nutrients in the subsurface layer important for biological activity. These biogeochemical changes can be expected to have significant economic consequences through changes in the availability of living marine resources.