Articles | Volume 19, issue 4
Research article
06 Jul 2023
Research article |  | 06 Jul 2023

Global variability of high-nutrient low-chlorophyll regions using neural networks and wavelet coherence analysis

Gotzon Basterretxea, Joan S. Font-Muñoz, Ismael Hernández-Carrasco, and Sergio A. Sañudo-Wilhelmy

We examine 20 years of monthly global ocean color data and modeling outputs of nutrients using self-organizing map (SOM) analysis to identify characteristic spatial and temporal patterns of high-nutrient low-chlorophyll (HNLC) regions and their association with different climate modes. The global nitrate-to-chlorophyll ratio threshold of NO3 : Chl > 17 (mmol NO3 mg Chl−1) is estimated to be a good indicator of the distribution limit of this unproductive biome that, on average, covers 92 × 106 km2 ( 25 % of the ocean). The trends in satellite-derived surface chlorophyll (0.6 ± 0.4 % yr−1 to 2 ± 0.4 % yr−1) suggest that HNLC regions in polar and subpolar areas have experienced an increase in phytoplankton biomass over the last decades, but much of this variation, particularly in the Southern Ocean, is produced by a climate-driven transition in 2009–2010. Indeed, since 2010, the extent of the HNLC zones has decreased at the poles (up to 8 %) and slightly increased at the Equator (< 0.5 %). Our study finds that chlorophyll variations in HNLC regions respond to major climate variability signals such as the El Niño–Southern Oscillation (ENSO) and Meridional Overturning Circulation (MOC) at both short (2–4 years) and long (decadal) timescales. These results suggest global coupling in the functioning of distant biogeochemical regions.

1 Introduction

High-nutrient low-chlorophyll (HNLC) areas are ocean regions where primary production should be potentially high but phytoplankton biomass remains relatively low and constant despite the perennial nutrient availability for growth (Martin and Fitzwater, 1988; Chisholm and Morel, 1991). They are interesting regions because they challenge the accepted paradigm of a positive relation between macronutrient concentrations and phytoplankton biomass in open waters but, most importantly, because they represent an important fraction of the global ocean carbon inventories, and therefore their extent influences the potential withdrawal of atmospheric CO2 to the deep ocean (Martin et al., 1990; De Baar et al., 1995; Boyd et al., 2005). It is estimated that HNLC biomes cover roughly between 20 % and 30 % of the world's oceans, comprising three major ocean areas: the subarctic North Pacific (SNP), the eastern equatorial Pacific (EEP), and most of the Southern Ocean (SO) (Martin, 1990; Coale et al., 1996; Parekh et al., 2005).

Because nitrogen is the mineral nutrient needed in greatest abundance by phytoplankton and owing to its generalized depletion in surface waters over much of the oceans, it is considered a key limiting nutrient for ocean production. In HNLC regions, where nitrogen is in excess, other non-exclusive factors such as rapid top-down control by zooplankton grazing, low irradiance, limitations by silicic acid availability, and/or iron (Fe) limitation have been hypothesized to explain the persistently low chlorophyll (Chl). While these factors may contribute to different degrees to the observed low Chl and determine the phytoplankton dynamics in HNLC regions (see Chavez et al., 1991; Cullen, 1995; Coale et al., 1996; Dugdale and Wilkerson, 1998; Landry et al., 2011), it is generally acknowledged that Fe availability is central to the productivity of HNLC regions (Boyd et al., 2007). All HNLC regions share a chronic Fe depletion in surface waters, and experimental results show highly positive productivity responses to Fe addition (Martin et al., 1994; Boyd et al., 2000, 2004; Tsuda et al., 2003; Coale et al., 2004). Indeed, Fe requirements are the largest among the trace metals for several metabolic processes, and not surprisingly, it has been considered the ultimate limiting nutrient (Moore and Doney, 2007). This has led to the proposal of a conceptual model of phytoplankton nutrient limitation in the modern ocean based on two functioning regimes: one in which the supply of nutrients is relatively slow with nitrogen availability limiting productivity and a complementary regime, with enhanced macronutrient supply, wherein Fe limits productivity (Moore et al., 2013).

Iron limitation influences the uptake of nitrogen, thereby explaining the unused nitrate concentrations in HNLC regions. Indeed, it has been proposed that a delicate balance between nitrogen and Fe availability modulates phytoplankton growth and that co-limitation is rather ubiquitous in the sea (Bryant, 2003; Browning et al., 2017). Other elements and compounds such as B vitamins, which are also scarce in Fe-limited areas, can also be co-limiting factors for phytoplankton growth in these regions (e.g., Koch et al., 2011; Bertrand et al., 2012). For example, it has been experimentally shown that the addition of Fe and B12 to Antarctic phytoplankton assemblages can synergistically increase phytoplankton growth (Bertrand et al., 2011; Cohen et al., 2017).

Despite their relevance for global ocean productivity and carbon fluxes, HNLC regions remain loosely defined, and knowledge of their temporal and spatial variability and trends is limited. Moreover, their response in a global warming scenario is uncertain. Only general aspects such as expected shifts in phytoplankton community composition and changes in Fe-cycling rates have been addressed to date (Fu et al., 2016; Lauderdale et al., 2020). The original description of HNLC systems by Minas et al. (1986) referred to a slowly growing phytoplankton standing stock despite the presence of high nutrient concentrations. However, there are no rigid criteria accurately defining the functioning of these ecosystems. Several ecosystem characteristics such as species composition, ecosystem structure, carbon utilization pathways, and response to climate change also differ between the HNLC and other ecosystems, reflecting differences in the limiting factor (Falkowski et al., 1998; Ono et al., 2008).

Of particular interest are the aspects related to the reduced variability and high permanence (i.e., temporal persistence) typically characterizing large HNLC regions. These features are distinctive from those of highly variable systems, which may temporarily present HNLC conditions. For example, some light-limited regions at high latitudes may present low productivity and enhanced nutrients during winter, but they respond to a transient situation that does not correspond to the generally accepted HNLC paradigm. Similarly, high nutrients and low Chl have been observed at the end of the spring bloom in some productive systems (Nielsdóttir et al., 2009; Birchill et al., 2019) and some areas located in coastal upwelling regions (Hutchins et al., 1998, 2002; Firme et al., 2003; Eldridge et al., 2004). While complying with the necessary conditions of high nutrients and low Chl, it is uncertain whether these ephemeral systems share structural and functional similarities with the large HNLC regions.

At a time when understanding biogeochemical responses to large-scale forcings, including climate change, has become a scientific priority, it seems appropriate to revisit some concepts of the functioning of HNLC regions. Their extent and variability are indicative of the dynamic changes in the bidirectional interrelationships of phytoplankton with the environment and with other organisms at large scales. Most of the information on the long-term variations of HNLC regions is drawn from global studies suggesting that their productivity is varying as a consequence of global warming and that they experience prominent interannual to decadal fluctuations superimposed on these long-term trends (i.e., Boyce et al., 2014; Martinez et al., 2020). Available evidence suggests that some HNLC regions may be decreasing in size as a result of increased ocean stratification (Ono et al., 2008). More recently, Yasunaka et al. (2016) determined that surface trends of phosphate and silicate in the North Pacific are associated with the shoaling of the mixed layer, reporting that surface nutrient concentration was correlated with the North Pacific Gyre Oscillation (NPGO). Some studies have shown that oligotrophic areas in the Northern Hemisphere are expanding (e.g., Polovina et al., 2008); however, with some exceptions (i.e., Radenac et al., 2012; Yasunaka et al., 2014), specific long-term studies on HNLC regions are scarce, and knowledge of their variability on the global ocean scale and their responses to climate change remains uncertain.

The objective of the present study is to provide a quantitative assessment of the large-scale patterns of variability of the three major HNLC regions (SNP, EEP, and SO) and their relationship with the main modes of climate variability. Systematically determining the boundaries of these HNLC regions has remained elusive since it requires coherent information on nutrients and Chl fields. The present study is based on the analysis of 20-year time series of monthly global ocean color data and nutrient concentrations from a biogeochemical model using machine learning techniques and wavelet analysis. First, based on the statistical analysis of global NO3 : Chl ratios, we determine a robust quantitative criterion to objectively define HNLC regions. Then we characterize the temporal variability patterns of HNLC regions based on their NO3 and Chl concentrations by using the self-organizing map (SOM) technique. We use the herein-established statistical criterion to assess the spatial variations of HNLC regions over the study period unveiled from the SOM analysis in the spatial domain of NO3 : Chl ratios. Finally, through a combined SOM–wavelet coherence analysis (WCA), we quantify the spectral power and the dynamic relationship between the observed Chl variability and two main global-scale forcings: the El Niño–Southern Oscillation (ENSO) and Meridional Overturning Circulation (MOC). We show that the combination of WCA with SOM-derived characteristic time series is an especially suitable tool for the analysis of driver–response relationships in the ocean.

2 Materials and methods

2.1 Ocean color data

We employ 20 years of monthly global composites of satellite Chl Level-3 products, derived from merging SeaWiFS, MERIS, MODIS AQUA, and VIIRS sensors using a Garver–Siegel–Maritorena (GSM) algorithm (Maritorena and Siegel, 2005), obtained from the GlobColour dataset (, last access: 12 June 2023). The chlorophyll product is spatially gridded, and the weighted average of the different merged Level-2 products is then calculated. The composite consists of a rectangular regular map product in degrees with a spatial resolution of 0.25 (i.e., around 28 km at the Equator) and covers the period from January 1998 to December 2017. We excluded results in the Arctic Ocean and the coastal Southern Ocean due to the interference of ice cover and prolonged gaps in the data. A total of 654 395 pixels were considered in the analysis. We are aware that the consistency of merged multi-mission ocean color satellite series may suffer from some limitations, influencing long trend analysis (Mélin et al., 2017). However, no significant increase or decrease is observed in the first-order trends of GlobColour data in more recent studies (e.g., Moradi, 2021). Therefore, while recognizing that some differences in regional and seasonal biases may occur in unified data products and acknowledging that discontinuities and trends of the median with time should be interpreted carefully according to the sensors used (Garnesson et al., 2019), merged Chl can be generally considered a good indicator of the magnitude of the overall phytoplankton trends.

2.2 Nitrate data

Since nutrient observations are still too scarce to allow obtaining time-resolved global-scale fields, we used global NO3 obtained from the biogeochemical hindcast model provided by Mercator-Ocean (, last access: 12 June 2023, see Fig. S1 in the Supplement). It consists of monthly mean fields of several biogeochemical variables at 0.25 horizontal resolution over the global ocean obtained using the PISCES model (Aumont et al., 2015). The model is forced by daily mean fields of ocean, sea ice, and atmospheric conditions. Ocean and sea ice forcings are obtained from the numerical simulation FREEGLORYS2V4 produced at Mercator-Ocean and the source of atmospheric forcings is the ERA-Interim reanalysis produced at ECMWF. Initial conditions are set from the World Ocean Atlas 2013 climatology. A complete model description can be found at We compared available observational nutrient data (NO3) from the sea surface (mean 0–20 m), obtained by merging bottle-cast data from the World Ocean Database (WOD18, Boyer et al., 2018;, last access: 12 June 2023) with model results. Generally, we found good agreement between nitrate in situ data and model results (r=0.98). Main deviations occur in the Southern Ocean where NO3 concentrations are overestimated (up to 7.2 mmol m−3).

2.3 Climatological data

Data on climate indices were obtained from available databases. The Bi-monthly Multivariate El Niño–Southern Oscillation Index (MEI.v2), hereafter the ENSO index, was obtained from the National Oceanic and Atmospheric Administration National Center for Environmental Prediction website (, last access: 12 June 2023). MOC data (Moat et al., 2022; Smeed et al., 2019) for the period (2004–2018) were obtained from the RAPID-WATCH MOC monitoring project (, last access: 12 June 2023).

2.4 Identification of HNLC regions

Presently, the best approximation to define the global distribution of HNLC regions in the world ocean is the use of National Oceanographic Data Center (NODC) maps of surface nutrients (, last access: 12 June 2023). However, excess nutrient availability by itself does not necessarily reflect HNLC conditions. In situ experiments are capable of discerning Fe limitation conditions, but a more manageable metric to assess the limits on the spatial extent of HNLC regions is required, in particular for remote sensing applications, as well as for allowing objective comparison between different environmental scenarios and studies.

To obtain a quantitative criterion for the definition of HNLC regions, we analyze the values of NO3 : Chl ratios (mmol mg−1) obtained from the SOM analysis on the time domain over the global ocean throughout the 20 years of data to identify a common statistical behavior representing HNLC conditions. In particular, we use the probability density function (PDF) of the extracted SOM NO3 : Chl temporal patterns to identify a threshold for defining HNLC conditions (PHNLC). Changes in the trend of the standard deviation calculated for each bin of the PDF function are employed to set a threshold ratio. To calculate the total extent of each region (km2), the spatial area of each pixel was calculated by considering its latitude.

2.5 Time and space domain SOM analyses

We use SOM (Kohonen, 1982) to elucidate spatial and temporal patterns in the complex relationship between nutrients and phytoplankton. SOM is a subtype of artificial neural network that uses an unsupervised machine learning algorithm to process and extract hidden structures in large datasets. The SOM algorithm is mainly based on a training process through which an initial neural network is transformed by iteratively presenting the input data. In this study, the architecture of the neural network is set in a sheet hexagonal map lattice of neurons, or units, to have equidistant neurons and to avoid anisotropy artifacts. Each neuron is represented by a weight vector with a number of components that is equal to the dimension of the input data vector, i.e., the number of rows or columns in the Chl and NO3 matrices, depending on whether the analysis is performed in either the temporal or spatial domain. We use an initial network composed of units of random values (random initialization). In each successive iteration during the training process, the neuron with the greatest similarity (excited neuron), called the best-matching unit (BMU), is updated by replacing its values with the Chl and NO3 values of the input sample data. The similarity is estimated by computing the Euclidean distance between the components of the input sample and the components of the weight vector of the unit. The unit most similar to the input sample is the one with the minimum distance. In the learning process, Chl and NO3 values of the topological neighboring neurons of the excited neuron (BMU) are also updated by replacing their values with values determined by a Gaussian neighborhood function. In these computations, we use the imputation batch training algorithm (Vatanen et al., 2015), whereby the SOM assumes that a single sample of data (input vector) contributes to the creation of more than one pattern, as the whole neighborhood around the best-matching pattern is also updated in each step of training. This yields a more detailed assimilation of particular features appearing on neighboring patterns. A final neural network with the NO3 : Chl patterns is obtained after repeating the training process until a stable convergence of the map is obtained.

For typical satellite datasets, the SOM can be applied to both space and time domains. By applying the SOM in the spatial domain, one can extract characteristic spatial patterns of the input data. If transposing the input data matrix and applying the SOM in the time domain, one can extract characteristic temporal patterns, i.e., the characteristic time series. Since each of these time series represents the temporal variability of a particular region, this method can be used to identify regions of differentiated variability on a map. The SOM, when applied to both space and time domains of the same data (called “dual SOM” analysis by Liu et al., 2016), provides a powerful tool for diagnosing ocean processes from these different perspectives. In this study, we focus on the second type. We have addressed the analysis separately in the time and space domains of the log-transformed NO3 and Chl datasets. In the time domain, we implement a [4×3] joint-SOM analysis of NO3 and Chl using as input weight vectors concatenating the time series of NO3 and Chl at each pixel, so each neuron corresponds to a characteristic joint NO3 and Chl temporal pattern over the total period of data. Since each pixel has an associated characteristic time series, we can obtain the location of a particular temporal pattern by computing the BMU for each pixel, providing a map of regions of differentiated NO3 : Chl temporal variability. For the analysis presented herein only the regions with NO3 : Chl >PHNLC are considered (regions R1 to R5).

An obstacle to the temporal domain analysis on a global scale is the opposed seasonality in both of Earth's hemispheres. The algorithm classifies the time series at each grid point corresponding to the period of the signal but does not consider time lags between the time series. Hence, pixels located in either the Northern or Southern Hemisphere displaying a similar significant period in the NO3 and Chl temporal variability are classified in the same regional pattern even if they are in antiphase when the signals are seasonally lagged (6-month delay). Regionalization is spatially coherent, but the seasonal variation in the characteristic pattern that represents the neuron mixes the phenological patterns of both hemispheres. Therefore, to properly analyze the properties and trends of each of the classified regions, we have calculated the mean features of the regions by segregating the grid points corresponding to each pattern obtained from the SOM analysis into the Northern Hemisphere, Equator, and Southern Hemisphere (see scheme in Fig. 1). Linear trends of NO3 and Chl concentrations in each region are estimated by decomposing the NO3 : Chl time series in a seasonal signal plus a residual component and applying Theil–Sen slope adjustment (Sen, 1968) of the residuals of the deseasonalized series. Correlation analyses were performed using the Pearson product moment correlation, computing best-fit linear trends using regression analysis.

Figure 1Scheme of a characteristic pattern and its distribution map obtained from SOM time domain analysis, as well as global and regional mean series calculated for each of the HNLC regions.

The SOM analysis in the spatial domain [3×3] array is addressed by using as input data weighted vectors consisting of spatial distributions over the global ocean of NO3 : Chl ratios at a particular time. The selection of the number of neurons depends on the complexity of the data, the features to be examined in the dataset, and the minimization of the errors. In this case, the resulting neurons after the training loop unveil the characteristic patterns describing the spatial variability of the HNLC regions on a global scale. Then, when computing the BMU for each time we designate the extracted characteristic spatial pattern that better describes the spatial distribution of NO3 : Chl ratios (P1 to P9) at each time, obtaining the time evolution of the characteristic spatial patterns throughout the considered period.

Because the SOM is based on the similarity computed from the Euclidean distance between samples, the input vectors of the different variables are normalized to the same range before initializing the SOM computations. This guarantees a consistent comparison of the weights of the components when computing the distance of two vectors.

The size of the neural network (number of neurons) depends on the number of samples and the complexity of the patterns. An optimal choice is important to maximize the quality of the SOM. In the present study, the map size is set to [4×3] with 12 neurons for the time domain analysis, and a [3×3] neural network is used in the spatial domain. Using larger map sizes, the patterns are slightly more detailed, and more regions of a particular variability emerge, but the occurrence of the probability of the patterns decreases, without affecting the results noticeably (Basterretxea et al., 2018; Hernández-Carrasco and Orfila, 2018). If a reduced neural map, such as [2×2], is used, patterns are concentrated together with the occurrence probability in a few rough patterns, but in this case increasing the topological error.

SOM computations have been performed using the MATLAB© toolbox of SOM v.2.0 (Vesanto et al., 1999) provided by the Helsinki University of Technology (, last access: 12 June 2023). Further information on SOM analysis is provided in the Supplement.

2.6 Combined SOM–wavelet coherence analysis

Joint SOM–wavelet power spectral analysis was demonstrated by Liu et al. (2016) in the study of characteristic time series of sea level variations in different regions of the Gulf of Mexico. Here, we expand it further to combined SOM–wavelet coherence analysis to assess the response of HNLC regions to global forcings. We use an approach based on the wavelet coherence analysis (WCA) between two time series (Grinsted et al., 2004; see the Supplement for further details). WCA characterizes cross-correlations by identifying the main frequencies, phase differences, and time intervals over which the relationship between the variability of HNLC regions and the main global forcings considered in this study, ENSO and MOC indexes, is strong. To do so, we first analyze the variability in both frequency and time of the characteristic time series of NO3 : Chl in the different HNLC regions extracted by the time domain SOM computations and the time series of the global forcings using the continuous wavelet transform (CWT).

Cross-wavelet transform (XWT) characterizes the association between the CWT of two signals, providing information on the common power and relative phase in the frequency–time domain of two time series. By applying the XWT to the NO3 : Chl ratios and climate forcings, we determine the cyclic changes in each of the HNLC regions and their relationship with the global forcings mentioned above. Finally, we quantify the correlation between the continuous wavelet transform of two signals using the wavelet coherence analysis (WCA). In the time–frequency space the wavelet coherence coefficient R2 is calculated as the squared absolute value of the smoothed cross-wavelet spectrum normalized by the product of the smoothed wavelet individual spectra for each scale (Torrence and Compo, 1998; Torrence and Webster, 1999; Grinsted et al., 2004). R2 is interpreted as a localized correlation coefficient in the frequency–time domain, and it takes values between 0 (no correlation) and 1 (perfect correlation). The statistical significance level of the wavelet coherence is estimated using Monte Carlo methods as described in Grinsted et al. (2004). We use the MATLAB software package (Grinsted et al., 2004) for wavelet coherence analysis. It should be noted that cross-wavelet analysis does not establish causative relationships but only allows identifying possible linkages between variables through the synchrony of their time series.

3 Results

3.1 Global characterization of HNLC regions

The mean pattern of global ocean color data for the 20 years analyzed reveals the well-known contrast in phytoplankton biomass between the highly productive areas located at high latitudes as well as in coastal upwelling regions and the low-latitude oceanic waters, where mean values are <0.1 mg m−3 (Fig. 2). Low-Chl regions generally correspond to low surface NO3 concentrations, whereas the opposite relationship (high nitrate and high chlorophyll) is more exceptional. Indeed, nutrient-rich productive waters are mainly restricted to shelf regions (coastal upwelling regions and shelf seas) or to the vicinity of islands (i.e., Falkland Islands) and other topographical features where multiple and overlapping sources of other elements, such as trace metals, are abundant (e.g., Pollard et al., 2009; Boyd and Ellwood, 2010). As shown in Fig. 2, only in the North Atlantic, the Bering Sea, and the eastern region of the Antarctic Peninsula is Chl enhanced. Conversely, a large part of surface ocean waters, particularly in the Southern Ocean and in the equatorial Pacific, corresponds to regions of relatively low Chl concentrations but with excess nitrate (i.e., >4 mmol m−3).

Figure 2(a) Global mean satellite Chl concentrations (mg m−3) for the period 1998–2018 and superimposed mean surface NO3 contour lines from modeling data. (b) Normalized probability density function of the values of the NO3 : Chl ratio obtained from the SOM temporal patterns. Green and blue lines show the fit to a normal distribution for the first and second PDF modes, respectively. (c) Standard deviation of the probability density function (PDF) of the NO3 : Chl (mmol mg−1) monthly ratios obtained for the 20 years analyzed. Note that the y-axis scale is logarithmic. The critical point ratio PHNLC=17 mmol NO3 mg Chl−1 delimits HNLC regions from macronutrient-limited regions.

The analysis of the normalized PDF of the NO3 : Chl extracted from the temporal SOM analysis provides a good discrimination criterion to define HNLC regions (shown in Fig. S2). As shown in Fig. 2b, the normalized PDF of the NO3 : Chl ratio displays a marked bimodal distribution with the main mode centered at low NO3 : Chl ( 5 mmol mg−1). The second mode, which corresponds to high-nutrient low-chlorophyll regions, is characterized by mean and standard deviation values of μ=24.1 and σ=6.7 mmol mg−1, respectively. A critical NO3 : Chl ratio bounds the lower limit of this distribution and can be estimated as μ-σ=17.4 mmol mg−1. Consistently, the PDF bulk analysis of its associated standard deviation (SD) function also reveals a clear critical value located where the value of the slope varies (Fig. 2c). Both analyses allow establishing a solid statistical criterion to infer a minimum value of PHNLC=17 mmol mg−1 for delimiting HNLC regions from other ocean regions. It is worth mentioning that while the PDF of the NO3 : Chl values obtained from the SOM analysis shows a bimodal mode, the bulk PDF of the raw NO3 : Chl values (i.e., without performing an SOM analysis) is unimodal.

From the 12 characteristic time patterns of NO3 : Chl variability obtained in the [4×3] SOM analysis, five display NO3 : Chl exceeding PHNLC at all times (not partially) throughout the entire study period (Fig. 3a). These associated subregions (R1 to R5) match the three traditionally reported HNLC regions (Fig. 3b). In these regions surface chlorophyll rarely exceeds 0.8 mg m−3, and the mean values range between 0.21 and 0.5 mg m−3 (Table 1). The global extent of these five SOM-identified HNLC subregions encompasses 25 % of the ocean, with the SO being by far the broadest region (18 % of the ocean), whereas SNP and EEP respectively occupy some 4 % and 3 % of the ocean. Besides the obvious absence of HNLC regions in the northern and central Atlantic, some latitudinal asymmetries are observed in the distribution of these regions. For example, the SO region extends to lower latitudes than the SNP (i.e.,  40 S), loosely coinciding with the South Antarctic Zone limit (SAZ; Orsi et al., 1995). Likewise, consistent with previous studies of this region (Radenac et al., 2012), the EEP displays a larger extent in the Southern Hemisphere (Fig. 3b).

Figure 3(a) Characteristic temporal variability patterns of NO3 : Chl ratios (N1–N12) as unveiled by the 4×3 SOM analysis in the time domain. Red dashed lines indicate the PHNLC value. (b) Coherent regions of HNLC variability (SNP, EEP, and SO) and corresponding subregions (R1 to R5) associated with the SOM temporal patterns exhibiting only NO3 : Chl values larger than PHNLC at all times throughout the entire analyzed period (i.e., identified with N9, N6, N2, N5, and N1). Patterns corresponding to a subregion in the Northern and Southern Hemisphere present a similar pattern although seasonally lagged (6-month delay). Insets show the time series of the averaged Chl over the corresponding subregion (complete map of regions of NO3 : Chl variability and corresponding temporal variability patterns are shown in Fig. S1). The red line represents the 24-month filtered series, and the blue line indicates the trend (values shown in Table 1).

Table 1Mean ± SD characteristics of each of the SOM-defined HNLC subregions (R1 to R5). NO3 (mmol m−3) and Chl (mg m−3) values are respectively from model and satellite data. Decadal chlorophyll trends (ΔChl, mgChl m−3 per decade) are calculated from the mean time series of monthly deseasonalized chlorophyll.

Download Print Version | Download XLSX

The global pattern obtained from the SOM analysis reflects a clear latitudinal zonation, which is mainly due to latitudinal variations in nutrient availability; while the chlorophyll concentration doubles along the latitudinal gradient (R1 to R5), NO3 increases up to 7-fold (see Table 1). It is noteworthy that nutrient concentrations are generally lower in the SNP (i.e., <17 mmol m−3) than in the SO, while Chl is comparatively higher (see Table 1). Indeed, R1 in SNP only achieves the NO3 : Chl criterion for HNLC regions during some periods. This region exhibits distinctive eastern and western provinces, which are consistent with previous studies describing the western region as more productive and variable (Imai et al., 2002; Harrison et al., 2004).

Major differences among the characteristic NO3 : Chl patterns in the defined subregions are due not only to variations in mean values, but also to the intensified seasonal variability at higher latitudes. For example, seasonality in Chl is particularly evident in R5 and less so in other polar subregions (Fig. 3b). Conversely, the seasonal component of variability in the EEP is masked by the intense short-term variability.

An interesting feature depicted from the temporal SOM analysis is the positive trend in Chl experienced in the HNLC regions located in polar areas, suggesting an increase in their productivity. Decadal tendencies are in the range of 0.04 to 0.06 mg m−3 per decade in the most productive subregions (R2 to R5 in SO and R3 in SNP) but become negligible at the Equator (Table 1). A regional average indicates a Chl increase of 0.6 % yr−1 in the SNP and 1.9 % yr−1 in the SO. Nevertheless, in the case of the SO positive trends are highly influenced by a positive Chl shift occurring at the end of 2010 when mean Chl increased by 24 % (see next section).

3.2 Spatial variability of HNLC regions

The set of nine coherent spatial patterns resulting from the SOM analysis in the space domain and their respective probabilities of occurrence are shown in Fig. 4. The organization of the maps in the figure reveals a hierarchical classification of the maps or spatial patterns. The most differentiated patterns, also displaying the highest probability of occurrence (the probability of finding a pattern similar to the input data), are located in the corners of the neural network, and transitional stages connecting these patterns fill the center. For example, along the top and left side patterns (P1, P2, P3, and P4), generally occurring during winter (see Figs. 4 and 5b), the SNP extends over a larger region compared to P7, P8, and P9, which display a 3 % decrease from the mean extent. Conversely, Fe limitation in the SO, as inferred from high NO3 : Chl ratios, is markedly enhanced towards the top and right side of the figure (P4, P5, P7, P8, and P9). The extent of the EEP region displays little variation. It should be noted that HNLC spatial extent and NO3 concentrations are not necessarily coupled since the boundaries also depend on Chl concentration values. In addition, patterns in the proximity of the Antarctic continent are, in some cases, not well-defined during winter due to ice cover in this region.

Figure 4Characteristic spatial patterns (P1 to P9) of HNLC regions as defined by NO3 : Chl ratios >PHNLC obtained from monthly data. The value on top of each pattern indicates its probability of occurrence over the 20-year period analyzed. To preserve the topology, the SOM algorithm introduces some patterns with zero probability of occurrence, such as P5. The color bar indicates the different NO3 : Chl ranges represented.

Figure 5 displays the time series of the BMUs and the monthly frequency of occurrence for each pattern. The main feature observed is the marked seasonality in the patterns shown in Fig. 4. The patterns with the highest probability of occurrence, P3 and P9 (100 % in April and 70 % in July, respectively), represent spring and summer situations in the Northern Hemisphere. P4 and P8 characterize transitions toward these patterns. Other patterns such as P6 and P2 (mostly occurring in winter and summer) are rarer but become more frequent after 2010 (Fig. 5a). As discussed below, this variation in HNLC regional patterns (i.e., P1 is no longer observed) suggests an abrupt and major transition towards more productive HNLC regions (higher Chl is observed).

Figure 5(a) Time evolution of the spatial patterns as defined by the best-matching units (BMUs) for the period of 1998–2018. The blue dashed line indicates the regime shift occurring after 2010. (b) Monthly frequency of occurrence of the spatial patterns identified in Fig. 4.


From the nine spatial patterns shown in Fig. 4, we estimated the seasonal and interannual variation in the extent of the HNLC regions (Fig. 6). Note that this regional partitioning is made on a global scale with global criteria and therefore leads to a large-scale smoothing, which could impact the values of the variation of the areas. However, as this signal smoothing is common to all the areas, this should not have any effect on the regional comparison of the area variation. The magnitude of these variations contrasts remarkably between the equatorial and polar regions. While the extent of the EEP only varies by 8.9 % seasonally, changes in SNP extent can exceed 100 % (Fig. 6). The peak in extent for the SNP corresponds to the boreal spring (63 % of the mean value in April). In the case of the SO, seasonality is mainly driven by changes related to the ice limit at high latitudes. Indeed, the extent of the HNLC region in the austral winter is <20 % of the mean annual extent.

Figure 6(a) Seasonal and (b) interannual variations in the spatial extent of the three HNLC regions, represented as a percentage of variation from the mean extent of each region. Variations refer to the mean extent of each region. Dark- and light-colored bars indicated positive and negative values, respectively. The blue dashed lines indicate the regime shift occurring after 2010.


A remarkably good inverse correlation (r=-0.97, n=20) is observed between the interannual variations in the extent of EEP and the SO, and a weaker though significant relationship exists between the SNP and the EEP (r=-0.50, n=20). Therefore, as the extent of HNLC in polar regions contracts, the equatorial region expands and vice versa. All three regions exhibit a shift in their annual extent after 2010 (Fig. 6). Both the SNP and the SO decrease after this year (5 % and 2.6 %), whereas the extent of the EEP slightly increases (0.4 %).

3.3 Climate drivers of HNLC region temporal variations

During the period considered here, 1997/1998, 2009/2010, and 2015/2016 showed particularly strong El Niño events, whereas weaker episodes occurred in 2002/2003, 2004/2005, and 2006/2007 (E1, E2, and E3 in Fig. 7a). The 1997/1998 and 2009/2010 El Niño events were followed by fast transitions to La Niña phase (L1 and L2 in Fig. 7a).

Figure 7(a) ENSO (MEI v2) index. E1–E3 indicate intense El Niño episodes, and L1 and L2 mark strong La Niña periods. Panels (b), (c), and (d) display the cross-wavelet coherence between NO3 : Chl ratio and ENSO for each of the HNLC regions. (e) Meridional overturning (AMOC) volume transport for the period 2004–2017 measured at 26.5 N (Smeed et al., 2018). The red line shows that the interannual component is obtained by filtering the data with a 540 d low-pass filter after the removal of the mean seasonal cycle. Panels (f), (g), and (h) display the cross-wavelet coherence between NO3 : Chl ratio and AMOC for each of the HNLC regions. The thick black contours in the cross-wavelet coherence figures designate the 95 % confidence levels, and the cone of influence where edge effects are not negligible is shown as a lighter shade. The arrows indicate the phase relationship between the signals, with the horizontal component indicating in phase (rightward) or out of phase (leftward) and the vertical component indicating a 90 phase difference lagging (upward) or leading (downward). Period units are months.


The WCA values between NO3 : Chl ratio in each HNLC region and ENSO are shown in Fig. 7b to d. Generally, small coherence structures are observed at semi-annual periods, in particular in the SO. However, the main coherence pattern corresponds to a band extending in the 2 to 4 years in the SNP and >2 years in the EEP. This coherence between NO3 : Chl and ENSO in the 2–4-year period is particularly clear after 2005 (Fig. 7b). In the EEP, the coherence between the two series expands to periods >4 years, but, unlike in the SNP region where the NO3 : Chl ratio is in phase with the ENSO signal, the signals are strongly anticorrelated in this case (antiphase: relative phase of 180 between the two signals).

Figure 7e shows the MOC transport index (hereafter MOI) measured at 26.5 N (Smeed et al., 2019). MOI displays intense interannual variability; however, a clear long-term variation is observed. Transport exceeds 17 Sv until 2009, but it weakens during 2010, stabilizing thereafter. As shown in Fig. 7f to h coherence with NO3 : Chl ratios is strongest at interannual timescales (1–1.5 years) when MOC is debilitated.

4 Discussion

4.1 Global characterization of HNLC regions

In the present study, we have addressed the extent of the HNLC regions, their long-term variability, and the potential drivers of these variations. Despite the relevance of precise characterization of the extent of this biome for the understanding of physical ocean processes (i.e., upwelling) and the estimation of the amount of carbon drawn into the ocean by phytoplankton, objectively determining the boundaries of HNLC regions has remained elusive as it requires coherent information on both nutrients and Chl. We demonstrate that a statistical approach, based on a threshold in the distribution of the global NO3 : Chl ratios, can robustly delineate these regions (PHNLC). This is also evidenced in Table 1 where the characteristics of each cluster obtained with SOM analysis are remarkably constant (i.e., low SD).

As in preceding studies (e.g., Moore et al., 2013), we assume that excess NO3 in surface oceanic areas is indicative of Fe limitation. This avoids relying on the scarce information available on Fe stress or more complex ecosystem modeling approaches. Inference of phytoplankton Fe stress from satellite ocean color data has been attempted but it is a methodology that still presents large uncertainties (Browning et al., 2014). Furthermore, while bioavailable Fe is known to be the primary limiting factor in this biome, the establishment of HNLC conditions is influenced by various other factors such as light availability, grazing pressure, rate of Fe remineralization, and community structure, highlighting the diverse interrelations among these factors. Despite these drawbacks, the method developed herein provides results consistent with previous descriptions of the large-scale spatial patterns of HNLC regions, mostly based on NO3 fields (i.e., Archer and Johnson, 2000; Ono et al., 2008; Fu and Wang, 2022). Also, the proposed method for biome definition may introduce a bias in that the resulting spatial fields are smoother compared to those based on Fe limitation, which is due to the greater variability of Fe concentrations compared to NO3 fields.

The PHNLC obtained from the PDF distribution of the NO3 : Chl ratios represents a statistical threshold that integrates complex biological processes, including competition for resources, grazing, changes in species composition, nutrient uptake rates, and Fe-regulated algal photochemistry. Unlike Redfield or C : Chl ratios, which respond to physiological factors within phytoplankton cells, PHNLC can be considered an environmental indicator of changes in the structure and functioning of marine phytoplankton.

According to our analysis, some 25 % of the ocean (18 % of the Earth's surface) corresponds to unproductive HNLC waters. With 83 % of the global HNLC biome extent, the SO is the largest region and it is the only one presenting clear latitudinal variation in the characteristic Chl patterns (Table 2 and Fig. 3). This is consistent with available descriptions of the physical and chemical properties of the SO, which tend to be across latitude due to the meridional structure of the MOC and because of the rapid zonal redistribution imposed by circumpolar currents (Orsi et al., 1995). The SNP and the EEP each constitute  8 % of the total HNLC. However, while the EEP remains relatively stable (cv= 5; Table 2), the SNP can change up to 2-fold (Fig. 4 and Table  2).

Table 2Basic statistics of the annual extent of each of the SOM-defined HNLC subregions during the analyzed period (1998–2017). Maximum variation (max. variation) is calculated as the difference between the maximum and the minimum extent.

Download Print Version | Download XLSX

Our analysis reveals marked decadal tendencies in Chl in the most productive subregions, ranging between 0.04 and 0.06 mgChl m−3 per decade, and negligible trends in the EEP (Table 1). Positive Chl trends at high latitudes and no significant tendencies at the Equator at the 95 % interval have been previously reported by Hammond et al. (2017). Indeed, climate change projections for the 21st century predict declines in global marine net primary production but increasing Southern Ocean productivity (Hauck et al., 2015; Moore et al., 2018). Nevertheless, it is noteworthy that our analysis shows that trends in some regions of the Southern Hemisphere are influenced by a Chl shift occurring in 2009–2010 (particularly R4 and R5). The possible source of this shift in phytoplankton biomass is more extensively discussed below. However, it is unlikely to be related to satellite data merging. SeaWiFS ended operations in 2010 and the MERIS sensor ceased in April 2012; however, a decrease in biomass would be expected from these changes (Garnesson et al., 2019).

4.2 Spatial variability of HNLC regions

SOM analysis allowed the characterization of the seasonal and interannual variability of HNLC regions (Figs. 5 and 6). The SO and the SNP are regions of seasonal extremes in productivity as a consequence of the large fluctuations in the environment that they experience. Important seasonal variability is observed in the SO, which is attributed to light limitation during winter and to variations in iron stocks occurring in surface waters (Tagliabue et al., 2014). Therefore, while presenting HNLC characteristics for most of the year, the SO exhibits distinct Chl variability patterns that are well-captured in the SOM regionalization and the characteristic patterns of each subregion. In the case of the SO, the changes in the extent are reduced (20 % from mean values), which we attribute to the large-scale nature of the physical processes regulating the productivity of this region. Mid-depth and deep-ocean waters communicate with the ocean surface after following a circuit route driven by ocean overturning circulation (Lumpkin and Speer, 2007). As reported by Smeed et al. (2018), the variability of the MOC flow system has an important decadal component associated with thermohaline forcing. This long-term component of variability could dominate over higher-frequency variability, but the length of the observational record of the AMOC is still insufficient to resolve variations at this scale.

Variability patterns in the SNP are attributed to both marked seasonality in the local forcings and fluctuations in the regional circulation patterns. Nishioka et al. (2021) discuss the sources of phytoplankton variability in this region in detail, distinguishing three cases of biological response to Fe in this region. The annual bloom in open waters would be controlled by the sedimentary Fe supplied from the continental marginal regions and circulated laterally through the intermediate layer (e.g., Lam et al., 2006; Cummins and Freeland, 2007; Takeda, 2011). These nutrients are upwelled to the surface by several mixing processes (winter mixing and eddy diffusive mixing), including interactions of tidal currents with the rough topography (Nishioka et al., 2007, 2020). Massive spring phytoplankton blooms that occur around the coastal boundary current would be the result of the large amounts of Fe originating from coastal shelf areas (Cullen et al., 2009). Finally, intense, sporadic, and patchy phytoplankton production is occasionally observed in both spring and summer in response to the atmospheric dust Fe supply.

Marked differences in variability are observed between the subpolar regions and the EEP, where seasonality is weak. Nevertheless, it should be noted that the EEP is a peculiar region that integrates subregions with 6-month out-phased seasonal variations (north and south of the Equator). On average, the southern equatorial region contributes more (62 % of the total extent) to the mean extent of the EEP, whereas the northern subregion determines most of the observed variability (Fig. S3).

At interannual scales, variations in the extent of HNLC regions are more modest (up to 5 %; Fig. 6) but all three regions are notably correlated, suggesting that their interannual variations are driven by global-scale processes. In particular, the SO and the EEP are highly correlated (r=0.97). A coupling between Southern Ocean productivity and equatorial productivity was suggested by Dugdale et al. (2002) based on the observed nutrient ratios. They proposed that both regions were connected, out of phase, by the meridian Sub-Arctic Mode Water. Nevertheless, the seaways in the Pacific circulation pathways are complex, and determining the overturning signature is challenging. Some studies have suggested that temperature anomalies subducted into the pycnocline at subtropical latitudes may not reach the Equator with any appreciable amplitude (Schneider et al., 1999). However, mass water balances in the equatorial Pacific reveal that the strength of the equatorial upwelling is related to variations in the Pacific overturning (PMOC; McPhaden and Zhang, 2002), and therefore, an influence on EEP extent could be presumed. Nevertheless, it should be remembered that the extent of the HNLC regions is not necessarily solely dictated by upwelling intensity; it is also defined by the interaction with adjacent circulation patterns (i.e., the subtropical gyres in this case).

The pathways of PMOC in the North Pacific are structured in a multicell configuration. The SNP does not ventilate the deep ocean at significant rates, and the PMOC cell at this latitude rather corresponds to an independently functioning intermediate water cell (Warren, 1983). PMOC in this region is reportedly weak (1–4 Sv) and extends no further than 50 N (Ishizaki, 1994; Yaremchuk et al., 2001). However, there is evidence showing the response of the SNP to changes in PMOC (Burls et al., 2017; Holzer et al., 2021). For example, it has been observed in TOPEX altimeter data that meridional overturning transport influences the basin-scale baroclinic circulation in the SNP. This would explain the lower, yet significant, correlation with the variations exhibited by SO and EEP at interannual scales (Fig. 6). Indeed, the shift observed in 2009–2010 is common to all three regions (although with a different sign).

The causes of the drastic 2009–2010 variation in the extent of HNLC regions are uncertain. Several ocean-scale changes have been reported in the period 2009–2011. For example, rapid warming, salinification, and a concurrent dissolved oxygen decline were observed for the Bermuda Atlantic Time-series Study (BATS) during the 2010s (Bates and Johnson, 2020). There is also evidence indicating that a decadal intensification of Pacific trade winds weakened in 2011 (Bordbar et al., 2019). Trade wind intensity variability in the equatorial Pacific region is associated with sea surface temperature (SST) anomalies, weakening of the equatorial divergence, and changes in upper-ocean thermal structure (England et al., 2014; Bordbar et al., 2017). The relationship between the equatorial wind intensity and the equatorial undercurrent strength is well-established (McPhaden, 1993). These atmospheric changes, affecting the upwelling of Fe in the EEP (and indirectly to other oceanic regions), agree with the hypothesis of Winckler et al. (2016), suggesting that ocean dynamics, not dust deposition, control the equatorial Pacific productivity. In the EEP the period of strongest seasonal upwelling is in the boreal summer (e.g., Philander and Chao, 1991), which is associated with an expansion of the HNLC in Fig. 6a. We also observe an expansion of the extent of the HNLC region in the EEP after 2021 (Fig. 6). Changes in these forcings could explain variations in the seasonal variability, but the long-term changes in HNLC extent would require more a permanent variation in the stoichiometric equilibria or in the functioning of the marine food webs.

4.3 Drivers of HNLC region variability

ENSO is the primary source of the interannual variability in the EEP, and its occurrence is related to the decline in nutrient supply. For example, the deepening of the thermocline during ENSO is associated with the depression of the Equatorial Undercurrent (EUC) that transports Fe across the basin from the western Pacific (Gordon et al., 1997), also reducing nutrient supply and phytoplankton productivity (Chavez et al., 1999; Strutton and Chavez, 2000). Sub-decadal fluctuations in Chl in the EEP region displaying a good correlation with the ENSO index have been reported before (Oliver and Irwin, 2008; Boyce et al., 2010). By contrast, the SO only shows weak evidence of this relationship, which suggests that ENSO is not a major forcing driving the variability of NO3 : Chl in this region. This is consistent with reports from Ayers and Strutton (2013), who did not find a significant relationship between nutrients in this region and ENSO.

It can be argued that differences in the response to ENSO are due to the different nature of the forcings driving nutrient supply in each region. While the EEP and the SNP seem to respond to ENSO-related changes, NO3 : Chl ratios in the SO are more stable and respond to annual and semi-annual variations. The coupling between the EEP and the SNP dynamics has been reported before. Qiu (2002) observed progressive shoaling of the Alaska gyre caused by a strengthening of the cyclonic circulation. The interannual variability of this gyre was connected to ENSO-related sea surface height anomalies originating in the tropics. Several large-scale climate pattern indexes are invoked to explain physical and biological fluctuations in the SNP. For example, Di Lorenzo et al. (2008) defined the North Pacific Gyre Oscillation (NPGO), which explains strongly correlated fluctuations of salinity, nutrients, and chlorophyll related to the circulation in the North Pacific Gyre. It is beyond the scope of the present work to assess the relationships of all these indexes with the variability in the extent of HNLC regions. Nevertheless, the proposed climate indexes for the Pacific present a strong relationship among them, which highlights the strong dynamical linkages between tropical and extratropical modes of climate variability in the Pacific basin and the important role played by ENSO (Di Lorenzo et al., 2013).

The interpretation of the influence of MOC on global ocean productivity is more challenging. A stronger MOC should result in faster upwelling of macronutrients and Fe at high latitudes as well as in increased Ekman transport of nutrients equatorward and subsequent subduction (Ayers and Strutton, 2013). Consequently, a weakened MOC should plausibly favor Fe recycling rates, as reported by Rafter et al. (2017).

At high latitudes, the weakening of the AMOC is coherent with a decrease in the extent of the SNP and the SO (Fig. 5). This anomaly, starting in 2009–2010, is a global feature also reflected in the Intertropical Convergence Zone (ITCZ) time series (see Green et al., 2017; Ibánhez et al., 2017), suggesting a strong atmosphere–ocean coupling with an impact on ocean productivity. It is not clear how a reduced flow would favor the increase in biomass at high latitudes. It has been proposed that a reduced AMOC from increased precipitation and melting sea ice could contribute to reduced vertical mixing, which may increase productivity in polar regions (Riebesell et al., 2009). Other studies (Martínez-Garcia et al., 2009) showed a relationship between AMOC and Chl variations, mainly due to the interaction of the main pycnocline and the upper-ocean seasonal mixed layer. In addition, some paleoclimatic studies have demonstrated that AMOC weakening can increase the productivity north of the Polar Front, but only if an increase in the atmospheric soluble Fe flux occurs (Muglia et al., 2018). Paleoceanographical records showing a strong correlation between proxies of aeolian Fe flux and productivity have been reported in this region (Kumar et al., 1995; Martínez-Garcia et al., 2009) but, in present times, dust deposition in this area is notably different (i.e., McConnell et al., 2007), and this effect is unlikely to be important at the timescales considered here. Complex ecosystem processes including competition for Fe with bacteria, Fe remineralization rates, and organic complexation processes could determine the phytoplankton response under future scenarios. Further, biomass build-up is driven not only by nutrient availability. Changes in biomass can be produced by variations in the thermocline depth, affecting the vertical distribution of phytoplankton. Nevertheless, variations in phytoplankton composition, physiological adjustments in cellular pigmentation, and grazing could also modulate Chl variability. Indeed, the prevailing food web structure may play an important role in Fe fertilization (Schmidt et al., 2016). At larger scales, there are still unresolved questions about the couplings occurring at different temporal scales. For example, MOC variations are known to interact with ENSO, influencing its amplitude and variance (Dong et al., 2006; Dong and Sutton, 2007; Timmermann et al., 2007).

5 Conclusions

Variations in the boundaries of the HNLC regions can provide an integrative view of how climate-scale ocean variations influence ocean productivity. We established a statistical criterion to identify HNLC regions from global Chl and NO3 data that sets the basis for systematic analyses of HNLC regions and their response to climate variations.

Our results suggest that while local-scale processes can determine certain aspects of the productivity of HNLC regions, their long-term patterns are strongly influenced by variations in global atmospheric and oceanic circulation. We observed significant interannual variations in the extent of HNLC regions (up to 5 % in Fig. 6), which are associated with anomalies in global ocean circulation patterns (i.e., Fig. 7e). Accordingly, our findings suggest a scenario in which HNLC regions are vulnerable to interbasin teleconnections rather than local forcings. These general patterns may be modulated by feedback between different forcing mechanisms. Furthermore, our analysis reveals a shift in phytoplankton biomass and HNLC variation patterns occurring at the end of 2010, which evidences the occurrence of fast transitions in ocean biogeochemistry. The underlying drivers of these regime shifts and the resulting biological responses to these ocean-scale changes require further investigation since they are a fundamental aspect of long-term variations in marine ecosystem functioning.

Finally, the present study highlights the importance of maintaining long and coherent datasets beyond satellite-borne information to be able to disentangle the different components of variability, particularly at long timescales, and to evaluate the impact of climate change on marine ecosystems. Most of the geochemical information at this scale (i.e., nutrient and Fe fields) will probably require further global sampling programs and refined modeling.

Data availability

All data included in the present study are accessible from the following publicly available repositories: CARINA (, Key et al., 2010), E.U. Copernicus Marine Service (, Mercator-Ocean, 2023), GEOTRACES (, GEOTRACES, Intermediate Data Product Group, 2021), GlobColour (, Maritorena et al., 2010), NERC (, Smeed et al., 2019), and NOAA (, Boyer et al., 2018).


The supplement related to this article is available online at:

Author contributions

Data were processed and analyzed mainly by GB, JSFM, and IHC. Writing was conducted by GB, JSFM, IHC, and SASW.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We acknowledge support for the publication fee from the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI). The present research was carried out within the framework of the activities of the Spanish Government through the “María de Maeztu Centre of Excellence” accreditation to IMEDEA (CSIC-UIB) (CEX2021-001198-M).

Financial support

This work was partially supported by a SIFOMED grant (CTM2017-83774-P) from the Ministerio de Ciencia, Innovación y Universidades, the Agencia Estatal de Investigación (AEI), and the Fondo Europeo de Desarrollo Regional (FEDER, UE). Gotzon Basterretxea was supported by the Salvador de Madariaga PRX18/00056 scholarship. Joan S. Font-Muñoz received funding from an individual postdoctoral fellowship, “Margalida Comas” (PD/018/2020), from Govern de les Illes Balears and Fondo Social Europeo. Ismael Hernandez-Carrasco received financial support from the project TRITOP (grant no. UIB2021-PD06) funded by the Universitat de les Illes Balears and FEDER (EU).

We acknowledge support for the publication fee from the CSIC Open Access Publication Support Initiative through its Unit of Information Resources for Research (URICI).

Review statement

This paper was edited by Aida Alvera-Azcárate and reviewed by two anonymous referees.


Archer, D. E. and Johnson, K.: A model of the iron cycle in the ocean, Global Biogeochem. Cy., 14, 269–279, 2000. 

Aumont, O., Ethé, C., Tagliabue, A., Bopp, L., and Gehlen, M.: PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies, Geosci. Model Dev., 8, 2465–2513,, 2015. 

Ayers, J. M. and Strutton, P. G.: Nutrient variability in subantarctic mode waters forced by the southern annular mode and ENSO, Geophys. Res. Lett., 40, 3419–3423, 2013. 

De Baar, H. J. W., de Jong, J. T. M., Bakker, D. C. E., Löscher, B. M., Veth, C., Bathmann, U., and Smetacek, V.: Importance of iron for plankton blooms and carbon dioxide drawdown in the Southern Ocean, Nature, 373, 412–415, 1995. 

Basterretxea, G., Font-Muñoz, J. S., Salgado-Hernanz, P. M., Arrieta, J., and Hernández-Carrasco, I.: Patterns of chlorophyll interannual variability in Mediterranean biogeographical regions, Remote Sens. Environ., 215, 7–17,, 2018. 

Bates, N. R. and Johnson, R. J.: Acceleration of ocean warming, salinification, deoxygenation and acidification in the surface subtropical North Atlantic Ocean, Commun. Earth Environ., 1, 33,, 2020. 

Bertrand, E., Saito, M., Lee, P., Dunbar, R., Sedwick, P., and DiTullio, G.: Iron Limitation of a Springtime Bacterial and Phytoplankton Community in the Ross Sea: Implications for Vitamin B12 Nutrition, Front. Microbiol., 2, 160,, 2011. 

Bertrand, E. M., Allen, A. E., Dupont, C. L., Norden-Krichmar, T. M., Bai, J., Valas, R. E., and Saito, M. A.: Influence of cobalamin scarcity on diatom molecular physiology and identification of a cobalamin acquisition protein, P. Natl. Acad. Sci. USA, 109, E1762–E1771, 2012. 

Birchill, A. J., Hartner, N. T., Kunde, K., Siemering, B., Daniels, C., Gonzalez-Santana, D., Milne, A., Ussher, S. J., Worsfold, P. J., and Leopold, K.: The eastern extent of seasonal iron limitation in the high latitude North Atlantic Ocean, Sci. Rep.-UK, 9, 1–12, 2019. 

Bordbar, M. H., Martin, T., Latif, M., and Park, W.: Role of internal variability in recent decadal to multidecadal tropical Pacific climate changes, Geophys. Res. Lett., 44, 4246–4255, 2017. 

Bordbar, M. H., England, M. H., Sen Gupta, A., Santoso, A., Taschetto, A. S., Martin, T., Park, W., and Latif, M.: Uncertainty in near-term global surface warming linked to tropical Pacific climate variability, Nat. Commun., 10, 1990,, 2019. 

Boyce, D. G., Lewis, M. R., and Worm, B.: Global phytoplankton decline over the past century, Nature, 466, 591–596, 2010. 

Boyce, D. G., Dowd, M., Lewis, M. R., and Worm, B.: Estimating global chlorophyll changes over the past century, Prog. Oceanogr., 122, 163–173, 2014. 

Boyd, P. W. and Ellwood, M. J.: The biogeochemical cycle of iron in the ocean, Nat. Geosci., 3, 675–682, 2010. 

Boyd, P. W., Watson, A. J., Law, C. S., Abraham, E. R., Trull, T., Murdoch, R., Bakker, D. C. E., Bowie, A. R., Buesseler, K. O., and Chang, H.: A mesoscale phytoplankton bloom in the polar Southern Ocean stimulated by iron fertilization, Nature, 407, 695–702, 2000. 

Boyd, P. W., Law, C. S., Wong, C. S., Nojiri, Y., Tsuda, A., Levasseur, M., Takeda, S., Rivkin, R., Harrison, P. J., Strzepek, R., Gower, J., McKay, M., Abraham, E., Arychuk, M., Barwell-Clarke, J., Crawford, W., Crawford, D., Hale, M., Harada, K., Johnson, K., Kiyosawa, H., Kudo, I., Marchetti, A., Miller, W., Needoba, J., Nishioka, J., Ogawa, H., Page, J., Robert, M., Saito, H., Sastri, A., Sherry, N., Soutar, T., Sutherland, N., Taira, Y., Whitney, F., Wong, S.-K. E., and Yoshimura, T.: The decline and fate of an iron-induced subarctic phytoplankton bloom, Nature, 428, 549–553,, 2004. 

Boyd, P. W., Strzepek, R., Takeda, S., Jackson, G., Wong, C. S., McKay, R. M., Law, C., Kiyosawa, H., Saito, H., and Sherry, N.: The evolution and termination of an iron-induced mesoscale bloom in the northeast subarctic Pacific, Limnol. Oceanogr., 50, 1872–1886, 2005. 

Boyd, P. W., Jickells, T., Law, C. S., Blain, S., Boyle, E. A., Buesseler, K. O., Coale, K. H., Cullen, J. J., de Baar, H. J. W., Follows, M., Harvey, M., Lancelot, C., Levasseur, M., Owens, N. P. J., Pollard, R., Rivkin, R. B., Sarmiento, J., Schoemann, V., Smetacek, V., Takeda, S., Tsuda, A., Turner, S., and Watson, A. J.: Mesoscale Iron Enrichment Experiments 1993–2005: Synthesis and Future Directions, Science, 315, 612–617,, 2007. 

Boyer, T. P., Baranova, O. K., Coleman, C., Garcia, H. E., Grodsky, A., Locarnini, R. A., Mishonov, A. V., Paver, C. R., Reagan, J. R., Seidov, D., Smolyar, I. V., Weathers, K. W., and Zweng, M. M.: World Ocean Database 2018, A. V. Mishonov, Technical Editor, NOAA Atlas NESDIS 87, (last access: 12 June 2023), 2018. 

Browning, T. J., Bouman, H. A., and Moore, C. M.: Satellite-detected fluorescence: Decoupling nonphotochemical quenching from iron stress signals in the South Atlantic and Southern Ocean, Global Biogeochem. Cy., 28, 510–524, 2014. 

Browning, T. J., Achterberg, E. P., Rapp, I., Engel, A., Bertrand, E. M., Tagliabue, A., and Moore, C. M.: Nutrient co-limitation at the boundary of an oceanic gyre, Nature, 551, 242–246, 2017. 

Bryant, D. A.: The beauty in small things revealed, P. Natl. Acad. Sci. USA, 100, 9647–9649, 2003. 

Burls, N. J., Fedorov, A. V, Sigman, D. M., Jaccard, S. L., Tiedemann, R., and Haug, G. H.: Active Pacific meridional overturning circulation (PMOC) during the warm Pliocene, Sci. Adv., 3, e1700156,, 2017. 

Chavez, F. P., Buck, K. R., Coale, K. H., Martin, J. H., DiTullio, G. R., Welschmeyer, N. A., Jacobson, A. C., and Barber, R. T.: Growth rates, grazing, sinking, and iron limitation of equatorial Pacific phytoplankton, Limnol. Oceanogr., 36, 1816–1833, 1991. 

Chavez, F. P., Strutton, P. G., Friederich, G. E., Feely, R. A., Feldman, G. C., Foley, D. G., and McPhaden, M. J.: Biological and chemical response of the equatorial Pacific Ocean to the 1997–98 El Niño, Science, 286, 2126–2131, 1999. 

Chisholm, S. W. and Morel, F. M. M.: What controls phytoplankton production in nutrient-rich areas of the open sea?, Am. Soc. Limnol. Oceanogr. Symp., 36, U1507–U1511, 1991. 

Coale, K. H., Johnson, K. S., Fitzwater, S. E., Gordon, R. M., Tanner, S., Chavez, F. P., Ferioli, L., Sakamoto, C., Rogers, P., and Millero, F.: A massive phytoplankton bloom induced by an ecosystem-scale iron fertilization experiment in the equatorial Pacific Ocean, Nature, 383, 495–501, 1996. 

Coale, K. H., Johnson, K. S., Chavez, F. P., Buesseler, K. O., Barber, R. T., Brzezinski, M. A., Cochlan, W. P., Millero, F. J., Falkowski, P. G., and Bauer, J. E.: Southern Ocean iron enrichment experiment: carbon cycling in high-and low-Si waters, Science, 304, 408–414, 2004. 

Cohen, N. R., A. Ellis, K., Burns, W. G., Lampe, R. H., Schuback, N., Johnson, Z., Sañudo-Wilhelmy, S., and Marchetti, A.: Iron and vitamin interactions in marine diatom isolates and natural assemblages of the Northeast Pacific Ocean, Limnol. Oceanogr., 62, 2076–2096, 2017. 

Cullen, J. J.: Status of the iron hypothesis after the open ocean enrichment experiment, Limnol. Oceanogr., 40, 1336–1343,, 1995. 

Cullen, J. T., Chong, M., and Ianson, D.: British Columbian continental shelf as a source of dissolved iron to the subarctic northeast Pacific Ocean, Global Biogeochem. Cy., 23, GB4012,, 2009. 

Cummins, P. F. and Freeland, H. J.: Variability of the North Pacific Current and its bifurcation, Prog. Oceanogr., 75, 253–265, 2007. 

Di Lorenzo, E., Schneider, N., Cobb, K. M., Franks, P. J. S., Chhak, K., Miller, A. J., McWilliams, J. C., Bograd, S. J., Arango, H., and Curchitser, E.: North Pacific Gyre Oscillation links ocean climate and ecosystem change, Geophys. Res. Lett., 35, L08607,, 2008. 

Di Lorenzo, E., Combes, V., Keister, J. E., Strub, P. T., Thomas, A. C., Franks, P. J. S., Ohman, M. D., Furtado, J. C., Bracco, A., and Bograd, S. J.: Synthesis of Pacific Ocean climate and ecosystem dynamics, Oceanography, 26, 68–81, 2013. 

Dong, B. and Sutton, R. T.: Enhancement of ENSO variability by a weakened Atlantic thermohaline circulation in a coupled GCM, J. Climate, 20, 4920–4939, 2007. 

Dong, B., Sutton, R. T., and Scaife, A. A.: Multidecadal modulation of El Niño–Southern Oscillation (ENSO) variance by Atlantic Ocean sea surface temperatures, Geophys. Res. Lett., 33, L08705,, 2006. 

Dugdale, R. C. and Wilkerson, F. P.: Silicate regulation of new production in the equatorial Pacific upwelling, Nature, 391, 270–273, 1998. 

Dugdale, R. C., Wischmeyer, A. G., Wilkerson, F. P., Barber, R. T., Chai, F., Jiang, M.-S., and Peng, T.-H.: Meridional asymmetry of source nutrients to the equatorial Pacific upwelling ecosystem and its potential impact on ocean–atmosphere CO2 flux; a data and modeling approach, Deep Sea Res. Pt. II, 49, 2513–2531, 2002. 

Eldridge, M. L., Trick, C. G., Alm, M. B., DiTullio, G. R., Rue, E. L., Bruland, K. W., Hutchins, D. A., and Wilhelm, S. W.: Phytoplankton community response to a manipulation of bioavailable iron in HNLC waters of the subtropical Pacific Ocean, Aquat. Microb. Ecol., 35, 79–91, 2004. 

England, M. H., McGregor, S., Spence, P., Meehl, G. A., Timmermann, A., Cai, W., Gupta, A. S., McPhaden, M. J., Purich, A., and Santoso, A.: Recent intensification of wind-driven circulation in the Pacific and the ongoing warming hiatus, Nat. Clim. Change, 4, 222–227, 2014. 

Falkowski, P. G., Barber, R. T., and Smetacek, V.: Biogeochemical Controls and Feedbacks on Ocean Primary Production, Science, 281, 200–206,, 1998. 

Firme, G. F., Rue, E. L., Weeks, D. A., Bruland, K. W., and Hutchins, D. A.: Spatial and temporal variability in phytoplankton iron limitation along the California coast and consequences for Si, N, and C biogeochemistry, Global Biogeochem. Cy., 17, 1016,, 2003. 

Fu, W. and Wang, W.: Biogeochemical Equilibrium Responses to Maximal Productivity in High Nutrient Low Chlorophyll Regions, J. Geophys. Res.-Biogeo., 127, e2021JG006636,, 2022. 

Fu, W., Randerson, J. T., and Moore, J. K.: Climate change impacts on net primary production (NPP) and export production (EP) regulated by increasing stratification and phytoplankton community structure in the CMIP5 models, Biogeosciences, 13, 5151–5170,, 2016. 

Garnesson, P., Mangin, A., Fanton d'Andon, O., Demaria, J., and Bretagnon, M.: The CMEMS GlobColour chlorophyll a product based on satellite observation: multi-sensor merging and flagging strategies, Ocean Sci., 15, 819–830,, 2019. 

GEOTRACES, Intermediate Data Product Group: The GEOTRACES Intermediate Data Product 2021 (IDP2021), NERC EDS British Oceanographic Data Centre NOC [data set],, 2021. 

Gordon, R. M., Coale, K. H., and Johnson, K. S.: Iron distributions in the equatorial Pacific: Implications for new production, Limnol. Oceanogr., 42, 419–431, 1997. 

Green, B., Marshall, J., and Donohoe, A.: Twentieth century correlations between extratropical SST variability and ITCZ shifts, Geophys. Res. Lett., 44, 9039–9047, 2017. 

Grinsted, A., Moore, J. C., and Jevrejeva, S.: Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlin. Processes Geophys., 11, 561–566,, 2004. 

Hammond, M. L., Beaulieu, C., Sahu, S. K., and Henson, S. A.: Assessing trends and uncertainties in satellite-era ocean chlorophyll using space-time modeling, Global Biogeochem. Cy., 31, 1103–1117, 2017. 

Harrison, P. J., Whitney, F. A., Tsuda, A., Saito, H., and Tadokoro, K.: Nutrient and plankton dynamics in the NE and NW gyres of the subarctic Pacific Ocean, J. Oceanogr., 60, 93–117, 2004. 

Hauck, J., Völker, C., Wolf-Gladrow, D. A., Laufkötter, C., Vogt, M., Aumont, O., Bopp, L., Buitenhuis, E. T., Doney, S. C., and Dunne, J.: On the Southern Ocean CO2 uptake and the role of the biological carbon pump in the 21st century, Global Biogeochem. Cy., 29, 1451–1470, 2015. 

Hernández-Carrasco, I. and Orfila, A.: The Role of an Intense Front on the Connectivity of the Western Mediterranean Sea: The Cartagena-Tenes Front, J. Geophys. Res.-Ocean., 123, 4398–4422, 2018. 

Holzer, M., DeVries, T., and de Lavergne, C.: Diffusion controls the ventilation of a Pacific Shadow Zone above abyssal overturning, Nat. Commun., 12, 4348,, 2021. 

Hutchins, D. A., DiTullio, G. R., Zhang, Y., and Bruland, K. W.: An iron limitation mosaic in the California upwelling regime, Limnol. Oceanogr., 43, 1037–1054, 1998. 

Hutchins, D. A., Hare, C. E., Weaver, R. S., Zhang, Y., Firme, G. F., DiTullio, G. R., Alm, M. B., Riseman, S. F., Maucher, J. M., and Geesey, M. E.: Phytoplankton iron limitation in the Humboldt Current and Peru Upwelling, Limnol. Oceanogr., 47, 997–1011, 2002. 

Ibánhez, J. S. P., Flores, M., and Lefèvre, N.: Collapse of the tropical and subtropical North Atlantic CO2 sink in boreal spring of 2010, Sci. Rep., 7, 1–9, 2017. 

Imai, K., Nojiri, Y., Tsurushima, N., and Saino, T.: Time series of seasonal variation of primary productivity at station KNOT (44 N, 155 E) in the sub-arctic western North Pacific, Deep-Sea Res. Pt. II, 49, 5395–5408, 2002. 

Ishizaki, H.: A simulation of the abyssal circulation in the North Pacific Ocean. Part II: Theoretical rationale, J. Phys. Oceanogr., 24, 1941–1954, 1994. 

Key, R. M., Tanhua, T., Olsen, A., Hoppema, M., Jutterström, S., Schirnick, C., van Heuven, S., Kozyr, A., Lin, X., Velo, A., Wallace, D. W. R., and Mintrop, L.: The CARINA data synthesis project: introduction and overview, Earth Syst. Sci. Data, 2, 105–121,, 2010 (data available at:, last access: 12 June 2023). 

Koch, F., Marcoval, M. A., Panzeca, C., Bruland, K. W., Sañudo-Wilhelmy, S. A., and Gobler, C. J.: The effect of vitamin B12 on phytoplankton growth and community structure in the Gulf of Alaska, Limnol. Oceanogr., 56, 1023–1034, 2011. 

Kohonen, T.: Analysis of a simple self-organizing process, Biol. Cybern., 44, 135–140, 1982. 

Kumar, N., Anderson, R. F., Mortlock, R. A., Froelich, P. N., Kubik, P., Dittrich-Hannen, B., and Suter, M.: Increased biological productivity and export production in the glacial Southern Ocean, Nature, 378, 675–680, 1995. 

Lam, P. J., Bishop, J. K. B., Henning, C. C., Marcus, M. A., Waychunas, G. A., and Fung, I. Y.: Wintertime phytoplankton bloom in the subarctic Pacific supported by continental margin iron, Global Biogeochem. Cy., 20, GB1006,, 2006. 

Landry, M. R., Selph, K. E., Taylor, A. G., Décima, M., Balch, W. M., and Bidigare, R. R.: Phytoplankton growth, grazing and production balances in the HNLC equatorial Pacific, Deep-Sea Res. Pt. II, 58, 524–535, 2011. 

Lauderdale, J. M., Braakman, R., Forget, G., Dutkiewicz, S., and Follows, M. J.: Microbial feedbacks optimize ocean iron availability, P. Natl. Acad. Sci. USA, 117, 4842–4849, 2020. 

Liu, Y., Weisberg, R. H., Vignudelli, S., and Mitchum, G. T.: Patterns of the loop current system and regions of sea surface height variability in the eastern Gulf of Mexico revealed by the self-organizing maps, J. Geophys. Res.-Ocean., 121, 2347–2366, 2016. 

Lumpkin, R. and Speer, K.: Global ocean meridional overturning, J. Phys. Oceanogr., 37, 2550–2562, 2007. 

Maritorena, S. and Siegel, D. A.: Consistent merging of satellite ocean color data sets using a bio-optical model, Remote Sens. Environ., 94, 429–440, 2005. 

Maritorena, S., Hembise Fanton d'Andon, O., Mangin, A., and Siegel, D. A.: Merged satellite ocean color data products using a bio-optical model: Characteristics, benefits and issues, Remote Sens. Environ., 114, 1791–1804,, 2010 (data available at:, last access: 12 June 2023). 

Martin, J. H.: Glacial-interglacial CO2 change: The iron hypothesis, Paleoceanography, 5, 1–13, 1990. 

Martin, J. H. and Fitzwater, S. E.: Iron deficiency limits phytoplankton growth in the north-east Pacific subarctic, Nature, 331, 341–343,, 1988. 

Martin, J. H., Fitzwater, S. E., and Gordon, R. M.: Iron deficiency limits phytoplankton growth in Antarctic waters, Global Biogeochem. Cy., 4, 5–12, 1990. 

Martin, J. H., Coale, K. H., Johnson, K. S., Fitzwater, S. E., Gordon, R. M., Tanner, S. J., Hunter, C. N., Elrod, V. A., Nowicki, J. L., and Coley, T. L.: Testing the iron hypothesis in ecosystems of the equatorial Pacific Ocean, Nature, 371, 123–129, 1994. 

Martinez, E., Gorgues, T., Lengaigne, M., Fontana, C., Sauzède, R., Menkes, C., Uitz, J., Di Lorenzo, E., and Fablet, R.: Reconstructing global chlorophyll-a variations using a non-linear statistical approach, Front. Mar. Sci., 7, 464,, 2020. 

Martínez-Garcia, A., Rosell-Melé, A., Geibert, W., Gersonde, R., Masqué, P., Gaspari, V., and Barbante, C.: Links between iron supply, marine productivity, sea surface temperature, and CO2 over the last 1.1 Ma, Paleoceanography, 24, PA1207,, 2009. 

McConnell, J. R., Aristarain, A. J., Banta, J. R., Edwards, P. R., and Simões, J. C.: 20th-Century doubling in dust archived in an Antarctic Peninsula ice core parallels climate change and desertification in South America, P. Natl. Acad. Sci. USA, 104, 5743–5748, 2007. 

McPhaden, M. J.: Trade wind fetch-related variations in equatorial undercurrent depth, speed, and transport, J. Geophys. Res.-Ocean., 98, 2555–2559, 1993. 

McPhaden, M. J. and Zhang, D.: Slowdown of the meridional overturning circulation in the upper Pacific Ocean, Nature, 415, 603–608, 2002. 

Mélin, F., Vantrepotte, V., Chuprin, A., Grant, M., Jackson, T., and Sathyendranath, S.: Assessing the fitness-for-purpose of satellite multi-mission ocean color climate data records: A protocol applied to OC-CCI chlorophyll-a data, Remote Sens. Environ., 203, 139–151, 2017. 

Mercator-Ocean: Global Ocean Biogeochemistry Hindcast product: GLOBAL_MULTIYEAR_BGC_001_029, E.U. Copernicus Marine Service [data set], last access: 12 June 2023. 

Minas, H. J., Minas, M., and Packard, T. T.: Productivity in upwelling areas deduced from hydrographic and chemical fields 1, Limnol. Oceanogr., 31, 1182–1206, 1986. 

Moat B. I., Frajka-Williams E., Smeed D. A., Rayner D., Johns, W. E., Baringer, M. O., Volkov, D., and Collins, J.: Atlantic meridional overturning circulation observed by the RAPID-MOCHA-WBTS (RAPID-Meridional Overturning Circulation and Heatflux Array-Western Boundary Time Series) array at 26N from 2004 to 2020 (v2020.2), British Oceanographic Data Centre – Natural Environment Research Council, UK [data set],, 2022. 

Moore, C. M., Mills, M. M., Arrigo, K. R., Berman-Frank, I., Bopp, L., Boyd, P. W., Galbraith, E. D., Geider, R. J., Guieu, C., and Jaccard, S. L.: Processes and patterns of oceanic nutrient limitation, Nat. Geosci., 6, 701–710, 2013. 

Moore, J. K. and Doney, S. C.: Iron availability limits the ocean nitrogen inventory stabilizing feedbacks between marine denitrification and nitrogen fixation, Global Biogeochem. Cy., 21, GB2001,, 2007. 

Moore, J. K., Fu, W., Primeau, F., Britten, G. L., Lindsay, K., Long, M., Doney, S. C., Mahowald, N., Hoffman, F., and Randerson, J. T.: Sustained climate warming drives declining marine biological productivity, Science, 359, 1139–1143, 2018. 

Moradi, M.: Evaluation of merged multi-sensor ocean-color chlorophyll products in the Northern Persian Gulf, Cont. Shelf Res., 221, 104415,, 2021. 

Muglia, J., Skinner, L. C., and Schmittner, A.: Weak overturning circulation and high Southern Ocean nutrient utilization maximized glacial ocean carbon, Earth Planet. Sc. Lett., 496, 47–56, 2018. 

Nielsdóttir, M. C., Moore, C. M., Sanders, R., Hinz, D. J., and Achterberg, E. P.: Iron limitation of the postbloom phytoplankton communities in the Iceland Basin, Global Biogeochem. Cy., 23, GB3001,, 2009. 

Nishioka, J., Ono, T., Saito, H., Nakatsuka, T., Takeda, S., Yoshimura, T., Suzuki, K., Kuma, K., Nakabayashi, S., and Tsumune, D.: Iron supply to the western subarctic Pacific: Importance of iron export from the Sea of Okhotsk, J. Geophys. Res.-Ocean., 112, C10012,, 2007. 

Nishioka, J., Obata, H., Ogawa, H., Ono, K., Yamashita, Y., Lee, K., Takeda, S., and Yasuda, I.: Subpolar marginal seas fuel the North Pacific through the intermediate water at the termination of the global ocean circulation, P. Natl. Acad. Sci. USA, 117, 12665–12673, 2020. 

Nishioka, J., Obata, H., Hirawake, T., Kondo, Y., Yamashita, Y., Misumi, K., and Yasuda, I.: A review: iron and nutrient supply in the subarctic Pacific and its impact on phytoplankton production, J. Oceanogr., 77, 561–587, 2021. 

Oliver, M. J. and Irwin, A. J.: Objective global ocean biogeographic provinces, Geophys. Res. Lett., 35, L15601,, 2008. 

Ono, T., Shiomoto, A., and Saino, T.: Recent decrease of summer nutrients concentrations and future possible shrinkage of the subarctic North Pacific high-nutrient low-chlorophyll region, Global Biogeochem. Cy., 22, GB3027,, 2008. 

Orsi, A. H., Whitworth III, T., and Nowlin Jr., W. D.: On the meridional extent and fronts of the Antarctic Circumpolar Current, Deep-Sea Res. Pt. I, 42, 641–673, 1995. 

Parekh, P., Follows, M. J., and Boyle, E. A.: Decoupling of iron and phosphate in the global ocean, Global Biogeochem. Cy., 19, GB2020,, 2005. 

Philander, S. G. H. and Chao, Y.: On the contrast between the seasonal cycles of the equatorial Atlantic and Pacific Oceans, J. Phys. Oceanogr., 21, 1399–1406, 1991. 

Pollard, R. T., Salter, I., Sanders, R. J., Lucas, M. I., Moore, C. M., Mills, R. A., Statham, P. J., Allen, J. T., Baker, A. R., and Bakker, D. C. E.: Southern Ocean deep-water carbon export enhanced by natural iron fertilization, Nature, 457, 577–580, 2009. 

Polovina, J. J., Howell, E. A., and Abecassis, M.: Ocean's least productive waters are expanding, Geophys. Res. Lett., 35, L03618,, 2008. 

Qiu, B.: Large-scale variability in the midlatitude subtropical and subpolar North Pacific Ocean: Observations and causes, J. Phys. Oceanogr., 32, 353–375, 2002. 

Radenac, M., Léger, F., Singh, A., and Delcroix, T.: Sea surface chlorophyll signature in the tropical Pacific during eastern and central Pacific ENSO events, J. Geophys. Res.-Ocean., 117, C04007,, 2012. 

Rafter, P. A., Sigman, D. M., and Mackey, K. R. M.: Recycled iron fuels new production in the eastern equatorial Pacific Ocean, Nat. Commun., 8, 1100,, 2017. 

Riebesell, U., Körtzinger, A., and Oschlies, A.: Sensitivities of marine carbon fluxes to ocean change, P. Natl. Acad. Sci. USA, 106, 20602–20609, 2009. 

Schmidt, K., Schlosser, C., Atkinson, A., Fielding, S., Venables, H. J., Waluda, C. M., and Achterberg, E. P.: Zooplankton gut passage mobilizes lithogenic iron for ocean productivity, Curr. Biol., 26, 2667–2673, 2016. 

Schneider, N., Miller, A. J., Alexander, M. A., and Deser, C.: Subduction of decadal North Pacific temperature anomalies: Observations and dynamics, J. Phys. Oceanogr., 29, 1056–1070, 1999. 

Sen, P. K.: Estimates of the regression coefficient based on Kendall's tau, J. Am. Stat. Assoc., 63, 1379–1389, 1968. 

Smeed, D. A., Josey, S. A., Beaulieu, C., Johns, W. E., Moat, B. I., Frajka-Williams, E., Rayner, D., Meinen, C. S., Baringer, M. O., and Bryden, H. L.: The North Atlantic Ocean is in a state of reduced overturning, Geophys. Res. Lett., 45, 1527–1533, 2018. 

Smeed, D., Moat, B. I., Rayner, D., Johns, W. E., Baringer, M. O., Volkov, D. L., and Frajka-Williams, E.: Atlantic meridional overturning circulation observed by the RAPID-MOCHA-WBTS (RAPID-Meridional Overturning Circulation and Heatflux Array-Western Boundary Time Series) array at 26N from 2004 to 2018, British Oceanographic Data Centre, National Oceanography Centre, NERC, UK [data set],, 2019. 

Strutton, P. G. and Chavez, F. P.: Primary productivity in the equatorial Pacific during the 1997–1998 El Niño, J. Geophys. Res.-Ocean., 105, 26089–26101, 2000. 

Tagliabue, A., Aumont, O., and Bopp, L.: The impact of different external sources of iron on the global carbon cycle, Geophys. Res. Lett., 41, 920–926, 2014. 

Takeda, S.: Iron and phytoplankton growth in the subarctic North Pacific, Terrapub, 4, 41–93, 2011.  

Timmermann, A., Okumura, Y., An, S.-I., Clement, A., Dong, B., Guilyardi, E., Hu, A., Jungclaus, J. H., Renold, M., and Stocker, T. F.: The influence of a weakening of the Atlantic meridional overturning circulation on ENSO, J. Climate, 20, 4899–4919, 2007. 

Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998. 

Torrence, C. and Webster, P. J.: Interdecadal changes in the ENSO–monsoon system, J. Climate, 12, 2679–2690, 1999. 

Tsuda, A., Takeda, S., Saito, H., Nishioka, J., Nojiri, Y., Kudo, I., Kiyosawa, H., Shiomoto, A., Imai, K., and Ono, T.: A mesoscale iron enrichment in the western subarctic Pacific induces a large centric diatom bloom, Science, 300, 958–961, 2003. 

Vatanen, T., Osmala, M., Raiko, T., Lagus, K., Sysi-Aho, M., Orešič, M., Honkela, T., and Lähdesmäki, H.: Self-organization and missing values in SOM and GTM, Neurocomputing, 147, 60–70, 2015. 

Vesanto, J., Himberg, J., Alhoniemi, E., and Parhankangas, J.: Self-Organizing Map in Matlab: the SOM Toolbox, in: Matlab DSP Conference 1999, Espoo, Finland, 16–17 November 1999, 35–40, (last access: 12 June 2023), 1999. 

Warren, B. A.: Why is no deep water formed in the North Pacific?, J. Mar. Res., 41, 327–347, 1983. 

Winckler, G., Anderson, R. F., Jaccard, S. L., and Marcantonio, F.: Ocean dynamics, not dust, have controlled equatorial Pacific productivity over the past 500,000 years, P. Natl. Acad. Sci. USA, 113, 6119–6124, 2016. 

Yaremchuk, M., Bindoff, N. L., Schröter, J., Nechaev, D., and Rintoul, S. R.: On the zonal and meridional circulation and ocean transports between Tasmania and Antarctica, J. Geophys. Res.-Ocean., 106, 2795–2814, 2001. 

Yasunaka, S., Nojiri, Y., Nakaoka, S., Ono, T., Whitney, F. A., and Telszewski, M.: Mapping of sea surface nutrients in the N orth P acific: Basin-wide distribution and seasonal to interannual variability, J. Geophys. Res.-Ocean., 119, 7756–7771, 2014. 

Yasunaka, S., Ono, T., Nojiri, Y., Whitney, F. A., Wada, C., Murata, A., Nakaoka, S., and Hosoda, S.: Long-term variability of surface nutrient concentrations in the North Pacific, Geophys. Res. Lett., 43, 3389–3397, 2016. 

Short summary
We examine global ocean color data and modeling outputs of nutrients using SOM analysis to identify characteristic spatial and temporal patterns of HNLC regions and their association with different climate modes. HNLC regions in polar and subpolar areas have experienced an increase in phytoplankton biomass over the last decades, particularly in the Southern Ocean. Our study finds that chlorophyll variations in HNLC regions respond to major climate variability signals.