Area of the extent of highly oligotrophic gyre (a threshold defined as chl a ≤ 0.07 mg m⁻³; Polovina et al., 2008, hereafter oligotrophic threshold chl a) was calculated from monthly SeaWiFS (R2017, completed January 2018) and Aqua-MODIS (R2017, completed December 2017) satellite ocean color derived chl a (Figure 1A). The SeaWiFS data record extends from autumn of 1997–2010 albeit with episodic gaps during 2008–2010 due to periodic sensor malfunction. Reprocessing (Franz et al., 2007; Franz, 2009) has removed much of the drift in the SeaWiFS and Aqua-MODIS instruments. Data were merged by comparing the per pixel seasonal cycle of SeaWiFS to that of MODIS-Aqua from 2003 to 2007, and then applying a seasonal correction to the historical SeaWiFS record. Merged data include seasonally adjusted SeaWiFS data until July of 2002, after which MODIS-Aqua data were used to December 2016. We also compared the trends between the OCI algorithm (Hu et al., 2012) to that of the previously standard OC3/OC4 (OCX) algorithms for SeaWiFS and MODIS-Aqua. The OCI algorithm is more accurate for [chl a] < 0.2 mg m⁻³, converges with the OCX values above 0.2 mg m⁻³, and has been the standard NASA chl a algorithm since 2014.

Near daily, gridded sea level anomalies (SLA) and absolute dynamic topography (ADT) fields derived from satellite altimetry were downloaded from AVISO (https://www.aviso.altimetry.fr/en/data/products.html). Eddy kinetic energy (EKE) was calculated from the near daily fields:

where u and v are the zonal and meridional velocities calculated from sea surface height anomalies. Monthly means of each field were calculated.

Because the determination of extent of features estimated from geostrophic flow lines can be confounded by the effect of steric sea level rise on ADT (Gille, 2014), we downloaded steric sea level anomalies (Levitus et al., 2012) for the upper 2,000 m from NOAA's National Center for Environmental Information (https://www.nodc.noaa.gov/OC5/3M_HEAT_CONTENT/index.html), and calculated the trend in steric sea-level rise from 1993 to 2015 on a per pixel basis. Dynamic height was adjusted by subtracting the interpolated steric sea-level trend from the monthly satellite dynamic height field. The gyre was then classified as the region contained within closing contours of original ADT and adjusted ADT (ADT_adj) at 1.05 m, the minimum height at which the contour lines recirculated (closed) completely to the western boundary (Figure 1B).

Seascapes were classified from Level 3 9-km merged chl a (described above), PAR, and SST from the Advanced Very High Resolution Radiometer (AVHRR) and MODIS-Aqua. Similar to chl a, monthly averages of PAR and SST were computed from the combined SeaWiFS and MODIS product; all three variables were binned to 1 degree to match model output. Model output of SST, PAR, and surface chl a from 1985 to 2005 were taken from the Biogeochemical Element Cycling Model embedded in the Community Earth System Model (CESM-BEC, Doney et al., 2009; Moore et al., 2013). Surface chl a concentration is sensitive to mixed layer depth in the subtropics (e.g., Letelier et al., 1996; Siegel et al., 2013); thus we chose the CESM-BEC because of its inclusion of photoacclimation. While this does not isolate the effect of photoacclimation, it assures that modeled chl a and satellite chl a are responding similarly to physical forcing, reflecting the phytoplankton physiological responses observed at Station ALOHA. All data sets were standardized prior to classification. Satellite and model seascapes were classified separately (Figures 1C,D).

Climatological seasonal means of PAR, SST, and chl a were defined for both data sets from 2003 to 2007. We used a probabilistic self-organizing mapping algorithm (PrSOM, Anouar et al., 1998) to reduce the 3-variable (SSTx,y,m, PARx,y, chl a x,y,m) spatiotemporal data set onto a 15 × 15 neuronal map resulting in 225 classes, each with its own 3-D weight based on the maximum likelihood estimation (3-D MLEs). The neural net size (15 × 15 nodes) was chosen to maximize sensitivity to mesoscale processes while preventing underpopulated nodes (defined as < 500 pixels). The 3-D MLEs were then further reduced using a hierarchical agglomerative clustering (HAC) with Ward linkages (Ward, 1963). This linkage method uses combinatorial, Euclidian distances that conserve the original data space with sequential linkages (McCune et al., 2002). Euclidian distances here are equivalent to within-group and total sum of squares.

Previously, we found that eight seascapes classified from seasonal climatologies of satellite-derived SST, chl a, and PAR represented seasonal shifts in biogeochemical patterns and planktonic assemblages (Kavanaugh et al., 2014a,b). Here, we followed the methodology of Kavanaugh et al. (2014b), with three exceptions. First, the seasonal variability of PAR was removed on a per-pixel basis to minimize discontinuities associated with the spring and fall transitions, with interannual and spatial variability remaining. Second, stepwise agglomerations were conducted using the Ward method until four seascapes for each of the model and satellite data sets were defined: a subtropical seascape, oligotrophic transition, a mesotrophic transition, and a subarctic seascape. These four regions bear the spatial signature of analogous seasonally evolving regions classified previously (Kavanaugh et al., 2014b) and represent a trade-off between variance explained in both the satellite and model records, and spatial match-up between satellite and modeled seascapes (Figure S1). Finally, once the seasonal and spatial vectors were classified, the means, variances and covariances within seascapes informed a multivariate Gaussian mixture model (GMM). Class assignments of individual months were then determined by their maximum posterior probabilities. New classes were then predicted from the CESM and satellite GMMs using the monthly means of chl a, SST, and PAR of the respective model and satellite time series.

We focus on the dynamics of the subtropical seascape and the boundary between the subtropical and oligotrophic transition seascape. Total area of oligotrophy, recirculation, and multivariate subtropical seascape was assigned based on the sum of the area of all pixels within the appropriate contour or identified seascape in the North Pacific. Surface area of each pixel was calculated by correcting for the spheroid-effect on distance between lines of longitude. For each time step, the distance from Station ALOHA to the subtropical seascape boundary was also calculated using the mean of the 5th quantile of all boundary distances in a quadrant 45° on either side of due NE from Station ALOHA Figure 1D); encroachment of less oligotrophic waters from the NE was most common. When the oligotrophic boundary encroached beyond (south and or east of) Station ALOHA, the distance values were negative. All area-based analyses were truncated at 10°N to minimize the effect of the equatorial upwelling region.

To understand the role of mesoscale variability on basin-scale geography, we examined the relationship between mean EKE and the difference in extent between the gyre defined by circulation and that defined by oligotrophy (see above). Pixels were categorized as belonging in (1) both the oligotrophic and physical gyres, (2) only in the physical gyre, or (3) only in the oligotrophic gyre. The mean EKE was calculated over the pixels for each of the three aforementioned groups.

To infer the effect of climatic forcing on the interannual dynamics of gyre geography, three different indices were used. The multivariate El Niño Southern Oscillation Index (MEI) is computed from a principal component analysis (PCA) of six variables including sea level pressure, zonal and meridional wind components, cloudiness, and sea surface and air temperatures (Wolter and Timlin, 1993). These data are available from NOAA (http://www.esrl.noaa.gov/psd/enso/mei/). The North Pacific Gyre Oscillation (NPGO, Di Lorenzo et al., 2008; http://www.o3d.org/npgo/) is the second dominant mode calculated from North Pacific sea surface height anomalies and is associated with accelerated North Pacific, Alaska, and California Currents (Chavez et al., 2011). The Pacific Decadal Oscillation (PDO, Mantua et al., 1997; Zhang et al., 1997) is associated with the dominant mode of North Pacific SST anomalies and the data are available from the University of Washington (http://jisao.washington.edu/pdo/).

In situ NPP, Total chl a (Tchl a, the sum of divinyl and monovinyl chl a), and in situ physical patterns in the subtropical North Pacific were assessed using archived data from Station ALOHA (http://hahana.soest.hawaii.edu/hot/hot-dogs/interface.html). Net PP was determined by daytime incubation in bottles spiked with ¹⁴C bicarbonate incubated in situ to maintain natural light and temperature (Letelier et al., 1996). Chl a and photosynthetic accessory pigments were measured by high performance liquid chromatography (HPLC) according to Wright et al. (1991). Mean pigment concentrations and NPP patterns were quantified separately for the light-saturated surface (0–45 m) and light-limited region of the euphotic zone that contained the deep chlorophyll maximum (DCM, 75–150 m). Potential temperature and bottle salinity were also recorded for the upper 45 m. Mixed layer depth was calculated using a potential density gradient threshold of 0.005 kg m⁻⁴ and a potential density offset of 0.125 kg m⁻³ from the surface (Karl and Lukas, 1996).

Cross correlation analyses were conducted to determine the correlations between climate and geographic metrics, the relationship between seascape dynamics and in situ conditions at Station ALOHA, and the time scales or lags at which the correlation was strongest. Cross correlation analysis was also conducted to determine the effect of EKE on spatial mismatch between geographic metrics. Prior to analysis, monthly climatological means were calculated by removing outliers (exceeding ± 3 standard deviations), then missing in situ data were interpolated using a 3-month LOESS filter to minimize the effect of within-year data density on interannual trends. Then, anomalies for each time series were calculated by subtracting the monthly climatological mean, calculated from the smoothed data, and then smoothing with a 12-month LOESS filter prior to correlation analysis. Statistical significance of correlation coefficients were determined using a critical value which adjusted degrees of freedom to account for autocorrelation in covariance (Glover et al., 2011). Correlations that exceeded this value are considered significant (p < 0.05).

Trends over time were adjusted to account for the effect of climate oscillations. We conducted a PCA across the NPGO, PDO, and MEI indices; multiple linear regression was conducted using the PCA scores as predictors and individual seasonally-detrended anomalies as responses. Trends were then a result of a simple linear fit to the residuals of the multiple linear regression model.

From 1985 to 2016, three to four cycles of the NPGO, MEI, and PDO were evident (Figure 2). In the early part of the record (1985–1998), the PDO and MEI were out of phase but moved into phase following the 1997–1998 El Niño and subsequent La Niña, remaining in phase for the rest of the data record (r = 0.6, p < 0.001). Throughout our study period, the NPGO was anticorrelated (r = −0.5, p < 0.001) with the other two indices, but also led the PDO and MEI by 5 and 8 months, respectively (Figure S2). Thus, any secular trend reported is that which remains after accounting for the effect of oscillations in the sign and magnitude of the climate indicators.

Patterns and trends among the different gyre geographic metrics varied and were associated with different climate forcing. The area of low-level chl a was dependent on algorithm, with the lower sensitivity of the OCX resulting in a larger expanse of oligotrophy (Figure 3A, Table 1). Large interannual oscillations were evident, with any linear trend being highly dependent on record length (Table S1). From 1998 to 2016, expansion of the area of low-level chl a occurred at a rate of 0.36% year⁻¹ (Table 1). Seasonal and interannual oscillations as well as weak increases (0.5% per year) in the physical extent of the gyre were evident from the altimetry-based time series (Figure 3B, Table 1). Trends were weaker in data adjusted for steric sea-level changes than those unadjusted for absolute dynamic height (ADT_adj and ADT, respectively, Table S1). As with the extent of oligotrophy, presence or magnitude of a trend depended on the length of the time series, with no trend evident from 1998 to 2006 (Figure 3B, Table S1), despite the strong ENSO signal during that span. Because modeled and satellite seascapes were highly correlated and in phase for both indices (Figures S1B,C; expansion: r = 0.64, p < 0.05; isolation r = 0.76, p < 0.05), the data records were combined and a single index of expansion or isolation was used that spanned the entire data record (1985–2016). There was evidence of slight expansion of the subtropical seascape from 1985 to 2016 (Table 1). Isolation distance also varied from year to year, although aseasonal, within-year variability was higher relative to the other metrics (Figure 3D). A linear trend was evident from 1985 to 2016 (but not over a shorter duration), with isolation distance increasing by 0.9% per year (Table 1).

Gyre extent, defined by oligotrophic threshold chl a, was positively correlated to the MEI, with positive MEI preceding expansion of low-level chl a area by 6–12 months (Figure 4A). The areal extent of recirculation was positively correlated with the NPGO and strongly negatively correlated with the PDO and MEI (Figure 4B). Similar to the area of low-level chl a, the subtropical seascape was strongly correlated to the MEI and was also related to the PDO (Figure 4C), whereas the isolation of Station ALOHA within the subtropical gyre was positively correlated with the NPGO and negatively correlated with the MEI and PDO (Figure 4D). While covering a much larger area, areal extent of multivariate seascapes was weakly correlated to that of threshold chl a patterns, (Figure S3A: r = 0.24, p < 0.05); however, patterns preceded the chl-only index by approximately 7 months. Multivariate seascape expansion was anti-correlated to both the dynamic topography extent described above (r = −0.48, p < 0.05), and the distance of Station ALOHA to the edge of the subtropical seascape (Figure S3B: r = −0.36, p < 0.05). Isolation distance was positively correlated to gyre extent as measured by geostrophy (r = 0.33, p < 0.05), with isolation shifts preceding that of gyre extent also by 7 months. Isolation was not correlated to the low-level chl a metric.

The NPSG is most contracted in winter and at its greatest extent in September for both definitions of the gyre (Figures 5A,B). The northern winter boundary of the physical gyre exceeds that of low-level chl a, and the summer eastern boundary extends less than the low-level chl-a. However, on interannual scales, the physical extent and low-level chl a-based extent were anti-correlated, due to different climate forcing. The physical extent is strongly positively correlated with the NPGO, whereas the extent of low-level chl a is positively correlated with the PDO and MEI. Increased EKE resulted in expansion of the areal extent of oligotrophic-threshold chl a. The result is that off of the Kuroshio, the spatial mismatch between the two definitions is minimized with increased EKE between their boundaries (Figure 5C), where instabilities associated with the Kuroshio extension appeared to blend the physical and chl a-based boundaries. On the eastern and southern edge, the spatial mismatch increases with increased EKE (Figure 5D), where expanded oligotrophy is associated with greater EKE, particularly in summer.

Local SST averaged over the upper 45 m at Station ALOHA exhibited strong oscillations in response to changing climate forcing over the data record, with warmer temperatures associated with positive ENSO and PDO and negative NPGO (Figure 6A, Figure S4). While there was a significant linear trend from 1989 to 2016 (Table 1: 0.010 ± 0.004°C year⁻¹, p < 0.01), it was very weak (R² = 0.04). Sea surface salinity exhibited large oscillations in the opposite direction of SST and also a relatively strong increasing trend over the data record (0.008 ± 0.001 psu year⁻¹, p < 0.001, R² = 0.18; Figure 6B). Mixed layer depth had large interannual variability over the data record, but also oscillated in response to changing surface salinity and temperature. A weak linear increase was apparent over the data record for mixed layer calculated via the surface offset method (0.29 ± 0.13 m year⁻¹, p = 0.02, R² = 0.02; Figure 6D), although none was apparent with the mixed layer depth calculated via the gradient method (Figure 6C, Table 1).

With the exception of isolated events (e.g., early 2001), seasonal patterns of T-chl a in the surface and at depth were out of phase (Figure 7A). Concentrations in the light-limited euphotic zone were approximately double that of the surface, and also increased over the data record (Table 1: 0.001 ± 0.0002 mg chl a m⁻³ year⁻¹, p < 0.001, R² = 0.08). No analogous trend was evident in surface chlorophyll.

NPP increased over the data record in both the surface and at depth (Figure 7B, Table 1). Mean rates in the surface were ~4x higher than that in the 75–150 m depth bin (Table 1: 5–7 mg C m⁻³ d⁻¹ vs. 1–2 mg C m⁻³ d⁻¹). However, the relative rate of increase in NPP over time was stronger in the light limited euphotic zone vs. the surface (Table 1: surface ~1% year⁻¹, trend = 0.053 ± 0.010 mg C m⁻³ d⁻¹year⁻¹, p < 0.001, R² = 0.09; At depth ~3% year⁻¹, trend = 0.045 ± 0.004 mg C m⁻³ d⁻¹year⁻¹, p < 0.001, R² = 0.34). Together, these trends contributed to a relatively strong linear increase in NPP integrated over 150 m (Figure 7C, 5.63 ± 0.85 mg C m⁻² d⁻¹year⁻¹, p < 0.001, R² = 0.18). Oscillations were also evident in both chl a and NPP time series.

Increased mixed layer depth followed isolation of Station ALOHA (Figure 8A, r = 0.32, p < 0.05) and preceded a contracted subtropical seascape (Figure 8B, r = −0.32, p < 0.05) by ~6 months. Salinity was positively associated with isolation but at longer lags (r = 0.28–0.32, lags > 20 months) and negatively associated with, but preceding expansion of the subtropical seascape (r = −0.35, p < 0.05, ~1year). SST was positively associated with subtropical seascape expansion (Figure 8B), but preceded expansion by ~1 year (r = 0.3, p < 0.05). These patterns generally followed what would be expected by the climate forcing of geography (Figure 4, Figure S4) with SST at Station ALOHA positively correlated with the PDO and MEI (and negatively with the NPGO), and surface salinity and mixed layer depth positively correlated with the NPGO (and negatively correlated with the PDO and MEI) (Figure S4). Time lags and wavelengths were much longer and broader when considering climate indices directly, rather than their geographic signatures.

Station ALOHA T chl a and NPP were also correlated to gyre geography following local physical forcing. Isolation was associated with increased chl a in the surface (Figure 8C, r = 0.32, p < 0.05), and increased surface and integrated NPP (Figure 8E, r = 0.32 to 0.4, p < 0.05) followed isolation by ~5 months. Decreased Tchl a at depth was strongly associated with expansion of the subtropical seascape (Figure 8D: r = −0.41, p < 0.05) as was decreased NPP throughout the water column (Figure 8F: r = −0.28 to −0.39, p < 0.05).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.