In this study, we used the “Nucleus for European Modeling of the Ocean” (NEMO) modeling framework (Madec, 2008). The NEMO configuration used displays a coarse resolution with 31 vertical levels and a 2° horizontal grid with a refined 0.5° resolution in the equatorial band. The model includes a biogeochemical component, the Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES; Aumont et al., 2015). PISCES is a model of intermediate complexity designed for global ocean applications (Aumont and Bopp, 2006), which uses 24 prognostic variables and simulates biogeochemical cycles of oxygen, carbon and the main nutrients controlling phytoplankton growth (nitrate, ammonium, phosphate, silicic acid, and iron). It simulates the lower trophic levels of marine ecosystems distinguishing four plankton functional types based on size: two phytoplankton groups (small = nanophytoplankton and large = diatoms) and two zooplankton groups (small = microzooplankton and large = mesozooplankton). Chl from PISCES (hereafter referred to as Chl_PISCES) is defined as the sum of the simulated diatoms and nanophytoplankton Chl content.

The NEMO-PISCES simulation is forced with atmospheric fields from the interannual Drakkar Forcing Set 5 (DFS5.2, Dussin et al., 2014) for wind, air temperature and humidity, precipitation, shortwave and longwave radiations. It is initialized with the World Ocean Atlas 2005 (WOA05) climatology for temperature, salinity, phosphate, nitrate and silicate (Garcia et al., 2006), while iron initial state is similar to the model climatology employed by Aumont and Bopp (2006). The model simulation was spun up using 3 repetitions of the 30 years’ DFS5.2 forcing set, and finally ran over 1979–2010.

Although successfully used in a variety of biogeochemical studies (e.g., Bopp et al., 2005; Gehlen et al., 2006; Lengaigne et al., 2007; Schneider et al., 2008; Steinacher et al., 2010; Tagliabue et al., 2010; Séférian et al., 2013; Aumont et al., 2015; Keerthi et al., 2017; Parvathi et al., 2017 and references therein), the ability of the PISCES model to reproduce satellite surface Chl is briefly illustrated in section “Evaluation of Chl_PISCES at global scale.”

Satellite surface Chl for Case I waters is provided by the Ocean Color – Climate Change Initiative (OC-CCI, hereafter referred to as Chl_OC–CCI) from the European Space Agency¹. This product combines multi-sensor, global, ocean-color products while attempting to reduce inter-sensor biases for climate research (Storm et al., 2013). OC-CCI extends the time series beyond that provided by single satellite sensors and is performant in terms of long-term consistency than other products from multi-mission initiatives (Belo Couto et al., 2016).

Only deep oceanic areas (depth > 200 m) are considered to avoid coastal waters where specific non-case-1 waters products are required. The Chl Level-3 product is binned on a regular 1° grid with a monthly resolution over January 1998–December 2010. This time period does not extend beyond 2010 to be consistent with the NEMO-PISCES simulation. Chl_OC–CCI is used to evaluate the PISCES model performances in Section “Evaluation of Chl_PISCES at global scale,” and to train the statistical method in Section “Application to satellite radiometric observations.”

The variability of phytoplankton biomass is driven in many regions of the world ocean and at many timescales by physical processes (e.g., Wilson and Adamec, 2002; Wilson and Coles, 2005; Kahru et al., 2010; Feng et al., 2015; Messié and Chavez, 2015). Our statistical architecture relates to 12 predictors and one biological variable (Chl). A sample thus refers to 13 variables. The 12 predictors (7 physical variables from NEMO-DFS5.2, 2 temporal and 3 spatial parameters) are detailed in Table 1, including their influence on Chl variations and the references supporting this influence.

We purposely limited the predictors to surface variables because our objectives are (1) to reconstruct Chl from physical observations, which are mainly available through remotely sensed surface data (oceanic observations below the surface are indeed usually not accessible at large spatial-scales or interannual time-scales); (2) to build a statistical scheme that can complement more complex numerical models (here, NEMO-PISCES) which simulate complex three-dimensional processes and are costly to run.

A first SVR is trained on physical predictors from NEMO and DFS5.2 vs. Chl_PISCES. The reconstructed Chl time-series is referred to as Chl_Svr–PISCES. A second SVR is trained using the same physical predictors but vs. satellite Chl observations (Chl_OC–CCI). The reconstructed Chl time-series is referred to as Chl_Svr–CCI.

Climate indices are provided by the National Oceanic and Atmospheric Administration (NOAA) website² : the AMO, the Multivariate El Niño Southern Oscillation (ENSO) Index (MEI) and the Interdecadal Pacific Oscillation (IPO).

The statistical reconstruction technique is based on a SVR. This method belongs to kernel methods in Statistical Learning Theory and relates to the Support Vector Machine (SVM, Vapnik, 1998). SVM is a kernel-based supervised learning method (Vapnik, 2000) developed for classification purpose in the early 1990s and then extended for regression by Vapnik (1995). The basic idea behind SVR is to map the variables into a new non-linear space using the kernel function, so that the regression task becomes linear in this space. The learning step estimates the parameters of the regression model according to a linear quadratic optimization problem, which can be solved efficiently. SVR also uses a robust error norm based on the principle of structural risk minimization, where both the error rates and the model complexity should be minimized simultaneously. Because SVR can efficiently capture complex non-linear relationships, it has been used in a variety of fields, and more specifically for oceanographic, meteorological and climate impact studies (Aguilar-Martinez and Hsieh, 2009; Descloux et al., 2012; Elbisy, 2015; Neetu et al., 2020), as well as in marine bio-optics (Kim et al., 2014; Hu et al., 2018; Tang et al., 2019).

Predictors and Chl are normalized by removing their respective average and dividing them by their standard deviations. Two SVR are trained over 1998-2010: one on Chl_PISCES and one on Chl_OC–CCI (Step A in Figure 1). This time period has been chosen as 1998 is the first complete year of the satellite Chl_OC–CCI time-series, and 2010 is the last year available of the modeled Chl_PISCES. The two resulting SVR schemes are applied on the NEMO-DFS5.2 physical predictors over 1979–2010. Finally, the annual means and standard deviations initially removed are applied to perform the back transformation and reconstruct either Chl_Svr–PISCES or Chl_Svr–CCI (Step B in Figure 1).

Considering a Gaussian kernel, SVR only involves the selection of two hyperparameters: the penalty parameter C of the error term and the kernel coefficient gamma, driving the reduction of the cost function. C and gamma values are 1 and 0.1, respectively when the SVR is trained on Chl_PISCES, and 2 and 0.3 when trained on Chl_OC–CCI (see details in the Supplementary Material and Supplementary Figure 1A). Sensitivity tests to an increasing portion of the sample total number (from 0.2 to 9% of the full dataset) used in the training process are performed (see Supplementary Material and Supplementary Figure 1B). The mean absolute error stabilizes for a sample number higher than 6%, suggesting that the SVR skills don’t improve much afterward. This observation combined with computational limitations lead us to present the 9% experiment hereafter.

The SVR skills to reconstruct Chl interannual to decadal variations are investigated performing Empirical Orthogonal Function analysis on Chl_PISCES, Chl_OC–CCI, Chl_Svr–PISCES and Chl_Svr–CCI. First, Chl data are centered and reduced (i.e., the monthly climatology is removed and the induced anomalies are divided by their standard deviations) to avoid an overly dominant contribution of high values on the analysis (Emery and Thomson, 1997) over the periods of interest (i.e., 1998–2010 or 1979–2010). A 5-month running mean is applied to focus on the interannual/decadal signal. The analysis is separately performed for the Atlantic, Pacific and Indian Oceans north of 40°S until 60°N, and for the 40°S–60°S region hereafter referred to as the Austral Ocean. Indeed, the large area covered by the Pacific Ocean and its dominant modes in climate variability (i.e., ENSO/IPO), could regionally dampen other modes of variability. Basin-scale spatial maps are then gathered to a global one, referred to as EOF. The associated time-series refer to as the Principal Components (PCs).

The ability of the NEMO-PISCES model to reproduce the satellite Chl over 1998–2010 is briefly presented here. Boreal winter and summer climatology from Chl_PISCES compare reasonably well with those of Chl_OC–CCI (Figure 2A vs. 2B and 2C vs. 2D). The model correctly represents the main spatial patterns with, for instance, higher Chl and a stronger seasonal cycle at high latitudes, despite an overestimated biomass in the Southern Ocean (Launois et al., 2015). The model also captures low Chl in the subtropical gyres, with some underestimation. This discrepancy may be explained by the lack of acclimation dynamics to oligotrophic conditions or by the assumption of constant stoichiometry either in phytoplankton or in organic matter in the model (Ayata et al., 2013; Aumont et al., 2015). The model underestimates Chl values in the equatorial Atlantic and Arabian Sea. In this latter region, mesoscale and submesoscale processes unresolved by the model have been shown to be of critical importance (Hood et al., 2003; Resplandy et al., 2011). Finally, the parameterization of nitrogen-fixing organisms not explicitly modeled in that PISCES version could explain the Chl_PISCES underestimation in the western Pacific in austral summer (Dutheil et al., 2018).

High Chl are accurately simulated in the eastern boundary upwelling systems. In two of the three main High Nutrient Low Chlorophyll (HNLC) regions, i.e. the equatorial Pacific and the eastern subarctic Pacific, the model successfully reproduces the moderate Chl_OC–CCI. However, the model overestimates Chl_OC–CCI east of Japan because of an incorrect representation of the Kuroshio current trajectory. This common bias in coarse resolution models (i.e., Gnanadesikan et al., 2002; Dutkiewicz et al., 2005; Aumont and Bopp, 2006) is potentially related to too deep mixed layer simulated in winter inducing very strong spring blooms (Aumont et al., 2015). In the Southern Ocean, the third and largest main HNLC region, the model overestimates Chl_OC–CCI values, especially during summer. However, the standard satellite algorithms that deduce Chl from reflectance tend to underestimate in situ observations by a factor of about 2–2.5, especially for intermediate concentrations (e.g., Dierssen and Smith, 2000; Kahru and Mitchell, 2010). It is to note that Chl in physical-biogeochemical coupled models is commonly overestimated in the Southern Ocean, and systematically underestimated in the oligotrophic gyres (Séférian et al., 2013).

The 1^st mode of the EOF analysis performed on interannual Chl displays close percent of total variance for Chl_OC–CCI and Chl_PISCES (16.6% vs. 21.1%, respectively). Their PCs in the Pacific Ocean are well correlated with the MEI (r = 0.71 and 0.89 with p = 0.0015 and p < 0.001, respectively; Figure 3C). PCs show the greatest positive values in January 1998 during the peak of the strong 1997/1998 El Niño event and the greatest negative values during the following La Niña beginning of 1999. The associated EOFs display a Chl horseshoe pattern (Figures 3A,B), reminiscent of the ENSO pattern on SST (Supplementary Figure 2; Messié and Chavez, 2012). While the tropical Pacific experiences a Chl decrease during El Niño events, the North and South Pacific display a Chl increase, and inversely during La Niña. This typical ENSO pattern is also related to remote Chl anomalies outside the Pacific induced by atmospheric teleconnections, such as a Chl decrease in the tropical North Atlantic and in the South Indian Ocean during El Niño. Although the Atlantic and Indian Ocean’s PCs are not correlated with the MEI (0.14 and 0.05, respectively), their EOFs are similar to those obtained from analysis performed at global scale (vs. basin scale here) and which have been largely discussed in the past (e.g., Behrenfeld et al., 2001, 2006; Yoder and Kennelly, 2003; Chavez et al., 2011). Chl_PISCES reasonably well captures the first mode of Chl_OC–CCI interannual variability over 1998–2010 in the Pacific and Atlantic Oceans, with 0.89 and 0.77 (p < 0.001) correlations between their PCs, respectively, but not in the Indian Ocean, where the PCs correlation is far weaker (0.13) and insignificant (Figures 3C–E).

In this section, the SVR uses the same physical predictors from NEMO-DFS5.2 as in Section “Synthetic reconstruction from a physical-biogeochemical ocean model,” but it is trained on satellite radiometric observations (e.g., Chl_OC–CCI). The same procedure is followed (see Supplementary Figures 5A,B). A first validation is performed for 20% of 9% of the full data set and over the 1998–2010 training period showing a high determination coefficient of 0.87 and RMSE of 0.37 between Chl_OC–CCI and Chl_Svr–CCI (Supplementary Figure 5C).

As expected, the regression lines between the whole dataset of Chl_OC–CCI vs. Chl_Svr–CCI for each oceanic basin and at global scale are farther away from the 1:1 line than for the synthetic study over the training period, but still remain close (higher slope than 0.8, except in the Austral Ocean; Figure 7). The SVR trained on NEMO-DFS5.2 predictors vs. satellite Chl is expected to be less efficient than the SVR trained on the coherent NEMO-DFS5.2-PISCES physical-biogeochemical dataset. Some of the biological interactions/processes (such as the diversity of the prey-predator relationships, the complexity of photoacclimation phenomena) are not yet optimally formulated by model equations inducing that Chl derived from numerical modeling is oversimplified compared to the complexity of the real ocean. Not to mention that satellite Chl may itself be partially affected by other components that are not Chl, such as colored dissolved organic matter (CDOM; Morel and Gentili, 2009) and suspended particulate matter (SPM). Phytoplankton can also adjust their intracellular Chl according to light and nutrient availability (e.g., Laws and Bannister, 1980; Behrenfeld et al., 2015). The induced Chl changes are no longer ascribed to changes in biomass. All these signatures on satellite Chl could explain Chl_Svr–CCI underestimation. Nevertheless, determination coefficients between Chl_Svr–CCI and Chl_OC–CCI remain high over the training time period (0.85, 0.89, and 0.86 for the Indian, Pacific and Atlantic Oceans, respectively, Figure 7).

The NRMSE between Chl_OC–CCI vs. Chl_Svr–CCI is lower than 20% over most of the global ocean (Figure 8A). Correlations higher than 0.9 (p < 0.001) are evident over large subtropical areas in the Atlantic, Indian and Pacific Oceans as well as in the Equatorial Pacific (Figure 8B). Interestingly, the SVR generally does a better job at reconstructing the satellite Chl than the modeled one (Figures 5A,C vs. Figure 8). NRMSE are higher at high latitudes and along the oligotrophic area boundaries, although to a less extent than for Chl_PISCES. Because Chl_OC–CCI can only be retrieved under clear sky conditions, gaps in satellite observations (especially during wintertime) likely alters the SVR learning and could explain such a degradation of Chl_Svr–CCI as moving toward high latitudes.

The SVR ability to replicate Chl_OC–CCI interannual variability is now investigated over 1998–2010 (Figure 9). In the Pacific Ocean, Chl_OC–CCI and Chl_Svr–CCI 1^st EOFs are close (Figure 9A vs. 9B), their PCs are highly correlated (r = 0.89, p < 0.001; Figure 9C), and their percent variance are similar (Table 4). As presented in Section “Evaluation of Chl_PISCES at global scale,” this mode of Chl variability can be attributed to ENSO, given their EOFs pattern as well as their PCs highly correlated with the MEI (r_OC–CCI/_MEI = 0.71 and r_Svr–CCI/_MEI = 0.91, with p = 0.0015 and p < 0.001, respectively). Interestingly, Chl_Svr–CCI EOFs are closer to Chl_OC–CCI than Chl_PISCES in several areas such as in the north-western Pacific, the south-western Atlantic and the Indian Ocean from Madagascar to the western coast of Australia (Figures 9A,B vs. Figure 3B). Consistently, correlations between Chl_OC–CCI and Chl_Svr–CCI PCs in the three basins and for the 1^st two modes are higher than between Chl_OC–CCI and Chl_PISCES (Table 4).

Chl_OC–CCI linear trends over 1998–2010 exhibit large areas of increase or decrease (red and blue areas in Figure 10A, respectively). Productive regions at high latitudes and along the equatorial and upwelling areas generally exhibit positive Chl_OC–CCI trends, albeit many underlying regional nuances. Contrastingly, trends are generally negative in the center of the gyres. These regional trends are consistent with those extracted from the first 13 years of the SeaWiFS record and discussed by Siegel et al. (2013) (see their Figures 5B, 8B). The negative trends in the oligotrophic gyres were also reported by Signorini et al. (2015) who attributed this behavior to MLD shallowing trends. Surface water density variability induced by changes in temperature and salinity, combined with wind stirring, are effective drivers of vertical mixing, which in turn control the renewal of nutrients from the rich-deep layers toward the euphotic zone. Thus, shallower MLD would decrease nutrient uplift and phytoplankton growth in the oligotrophic areas.

Chl_Svr–CCI trends agree qualitatively well with those of Chl_OC–CCI at global scale (Figure 10B vs. 10A, respectively). Indeed, decline of Chl_Svr–CCI can be observed in the center of the gyres, while outside Chl_Svr–CCI generally increases in a similar way to Chl_OC–CCI. Chl_Svr–CCI accurately captures the largest Chl_OC–CCI increase observed in the Southern Ocean along the Antarctic Circumpolar Current. While Gregg and Casey (2004) reported a substantial negative bias in the SeaWiFS data for this region when compared to in situ observations, which could hamper the reliability of satellite trends discussed in this area (e.g., Siegel et al., 2013), the SVR remains able to reproduce the positive observed trend. Despite qualitative spatial agreements, it is noteworthy that the SVR underestimates by half the magnitude of the satellite trend (see scales in Figure 10A vs. 10B).

Interestingly, trends in Chl_PISCES generally differ from Chl_OC–CCI (Figure 10C). This is striking for the North Pacific and Atlantic high latitudes, but also in the equatorial Atlantic and Arabian Sea with opposite trends when compared with Chl_OC–CCI and Chl_Svr–CCI, and in a more mitigated manner in the Austral Ocean.

Chl_Svr–CCI is also compared with the only historical consistent dataset built by Antoine et al. (2005) who reanalyzed ocean color time series from CZCS (1979–1983) and SeaWiFS (1998–2002). A 22% global mean increase of Chl between the two era was reported. It was mainly due to large increases in the intertropical areas and to a lesser extent in higher latitudes, while oligotrophic gyres displayed declining concentrations (Figure 11A). SST from the SODA reanalysis was used as a proxy of ocean stratification and opposite Chl and SST changes over 60% of the ocean between 50°S and 50°N was reported (light blue and yellow in Figure 11B, adapted from Martinez et al., 2009). This inverse relationship was used to hypothesized that multidecadal changes in global phytoplankton abundances were related to basin-scale oscillations of the ocean dynamics. Briefly, SST changes were related to a regime shift of the PDO (although the use of the basin-scale IPO would have been more appropriate) from a warm to a cold phase in the Pacific and Indian Oceans leading to an increase of Chl, and inversely in the Atlantic Ocean with a regime shift from a cold to a warm phase of the AMO leading to a Chl decrease.

Observed Chl changes over the last decades are accurately reproduced by Chl_Svr–CCI, including a Chl increase in the equatorial Pacific and the southern tropical Indian Oceans, as well as a Chl decline in both the Atlantic and Pacific oligotrophic gyres (Figure 11C). However, the magnitude of the SVR reconstructed Chl is underestimated (note that the Chl ratio is multiplied by 2 in Figure 11C to allow the comparison with Figure 11A). On average, the inverse relationship between Chl_Svr–CCI and SST_NEMO (Figure 11D) occurs over 69.4% of the global ocean between 50°S and 50°N in a similar way to that reported by Martinez et al. (2009), especially in the Pacific Ocean (see Figure 11D vs. 11B). In the Indian Ocean, although Chl mainly increases in both studies, it is here associated with a SST decrease. Interestingly, this inverse Chl-SST relationship in the Indian Ocean (yellow area in Figure 11D) was reported in Behrenfeld et al. (2006) over the SeaWiFS era, suggesting that the SST dataset used in Martinez et al. (2009) may have decadal discrepancies for this region.

In their study, Martinez et al. (2009) analyzed two 5-year time periods apart from each other by 15 years. They suggested that averaging observations separately over the two time-periods may have dampen the effect of interannual variability and reveal the decadal one. Most of the changes observed between the time periods covered by the two satellites are here confirmed based on the reconstructed Chl_Svr–CCI. However, the continuous 30-year time series of Chl_Svr–CCI provides new insights on the observed regime shifts (Figure 12). In the Pacific Ocean, the 1^st EOF of Chl_Svr–CCI (Figure 12A) is close to the Chl spatial patterns obtained from the CZCS to SeaWiFS era (Figure 11C) and the PC remains highly correlated with the IPO over 1979-2010 (r = 0.94 with p < 0.001, Figure 12B). The Chl increase in the Indian Ocean, north-east of Madagascar toward the west coast of Australia, between the 1980’s and the 2000’s also appears on the Chl_Svr–CCI EOF. These temporal changes might also be related to the IPO variability (correlation between the IPO index and Chl_Svr–CCI PC = 0.6, p < 0.001; Figure 12C).

In the Atlantic Ocean, CZCS-SeaWiFS Chl and Chl_Svr–CCI 1^st EOF also share some similarities, including a decrease of Chl in the subtropical gyres and an increase in the equatorial/tropical regions. The associated PC (Figure 12D), exhibits a shift between 1979–1983 and 1998–2002 consistently with Figure 2C of Martinez et al. (2009). In this latter study, this change was attributed to a regime shift of the AMO. However, the AMO index is not correlated with the 1^st Chl_Svr–CCI PC (r = 0.03, p = 0.43) but rather with the 2^nd mode (r = 0.43 with p = 0.003, Supplementary Figure 6), likely explaining the spatial discrepancies in Figure 11A vs. 11C. Although the detailed analysis of Chl decadal variability is beyond the scope of the present study, these initial findings underscore the importance of continuous time series at regional/global scales to combine spatial and temporal information’s and properly investigate Chl long-term variability.