Primary Production, an Index of Climate Change in the Ocean: Satellite-Based Estimates over Two Decades

: Primary production by marine phytoplankton is one of the largest ﬂuxes of carbon on our planet. In the past few decades, considerable progress has been made in estimating global primary production at high spatial and temporal scales by combining in situ measurements of primary production with remote-sensing observations of phytoplankton biomass. One of the major challenges in this approach lies in the assignment of the appropriate model parameters that deﬁne the photosynthetic response of phytoplankton to the light ﬁeld. In the present study, a global database of in situ measurements of photosynthesis versus irradiance (P-I) parameters and a 20-year record of climate quality satellite observations were used to assess global primary production and its variability with seasons and locations as well as between years. In addition, the sensitivity of the computed primary production to potential changes in the photosynthetic response of phytoplankton cells under changing environmental conditions was investigated. Global annual primary production varied from 38.8 to 42.1 Gt C yr − 1 over the period of 1998–2018. Inter-annual changes in global primary production did not follow a linear trend, and regional differences in the magnitude and direction of change in primary production were observed. Trends in primary production followed directly from changes in chlorophyll- a and were related to changes in the physico-chemical conditions of the water column due to inter-annual and multidecadal climate oscillations. Moreover, the sensitivity analysis in which P-I parameters were adjusted by ± 1 standard deviation showed the importance of accurately assigning photosynthetic parameters in global and regional calculations of primary production. The assimilation number of the P-I curve showed strong relationships with environmental variables such as temperature and had a practically one-to-one relationship with the magnitude of change in primary production. In the future, such empirical relationships could potentially be used for a more dynamic assignment of photosynthetic rates in the estimation of global primary production. Relationships between the initial slope of the P-I curve and environmental variables were more elusive. T.Y.; writing—original draft preparation, G.K.; writing—review and editing, S.S., T.P., M.B., R.J.W.B., M.E., H.G.G., K.G., T.I., Ž.K., E.M., W.H.v.d.P., G.H.T. and V.v.D.-V.; visualisation, G.K. and B.F.J.; supervision, S.S. and T.P.; project administration, G.K. S.S.; acquisition, G.K., S.S. T.P. authors agreed to the published version of the manuscript.


Introduction
The oceans play a key role in biogeochemical processes on Earth. Phytoplankton are responsible for almost half of the total global net primary production [1][2][3][4][5]. This does not only provide the basis for the marine food web, but also has a strong impact on carbon sequestration in the ocean's interior [6]. Marine primary production, estimated to be of the order of 50 Gt C per annum [2][3][4][5]7], is one of the largest fluxes of carbon on our planet. Because of its importance, phytoplankton primary production has received considerable attention from the scientific community. Studies based on in situ observations are now supplemented by satellite-based calculations to estimate global primary production patterns at high spatial and temporal resolutions. Yet, trends in biological fields estimated from remote-sensing observations have not been taken into account in recent studies on global carbon budgets and pools and fluxes of carbon in the ocean [8,9]. In recent years, considerable efforts have been made to correct inter-sensor biases and merge data from multiple ocean-colour satellite sensors to provide a long (over two decades) record of phytoplankton biomass in the global oceans through the Ocean Colour Climate Change Initiative of the European Space Agency [10]. This time series now offers the opportunity to undertake a systematic study of changes in phytoplankton primary production over the last 20 years.
Phytoplankton primary production is forced by physico-chemical conditions in the water column, including temperature, light and micro-and macronutrients. These drivers are influenced by seasonal, inter-annual and multidecadal variations in oceanic and atmospheric processes. For example, phytoplankton primary production in polar regions is strongly influenced by seasonal solar irradiance patterns and the formation of surface mixed layers due to ice melt in spring and summer [11][12][13][14]. In contrast, at lower latitudes where trade winds prevail, phytoplankton primary production can be nutrient-limited year-round and seasonal patterns are less obvious [15,16]. Superimposed on seasonal cycles are the variations associated with inter-annual and multidecadal ocean-atmospheric oscillations. These oscillations are associated with anomalies in Sea Surface Temperature (SST), precipitation and wind patterns, leading to changes in water column stability and nutrient loading into the euphotic zone [17][18][19]. The El Niño-Southern Oscillation (ENSO), North Pacific Gyre Oscillation (NPGO), Indian Ocean Dipole (IOD) and Atlantic Multidecadal Oscillation (AMO) have all been shown to affect phytoplankton primary production [17][18][19][20][21]. These natural variations in water column conditions can cause a 10-fold variation in primary production between different regions, with low-nutrient subtropical waters at the lower end of production and highly eutrophic coastal regions at the upper end [22,23].
Given these natural variations in physico-chemical conditions and in phytoplankton primary production, it is expected that the physical changes associated with climate change will redistribute phytoplankton primary production. Over the past decades, increases in SST and ocean heat content, along with enhanced precipitation relative to evaporation and sea ice melt, have caused significant variations in physico-chemical conditions of the water column [23,24]. Subsequent changes in temperature and density stratification, and nutrient loading into the euphotic zone, are expected to affect phytoplankton growth and primary production under global climate change [23,24]. Several studies based on in situ, satellite and/or modelling observations have shown that changes in global primary production associated with climate change ranged from a 0.57-13% decrease [25][26][27][28] to a 2% increase [29]. Discrepancies between these estimates may be based on differences in methodology or in variations in temporal and spatial scales. It therefore seems that care has to be taken in estimating global primary production, especially considering that some regions will experience additional local forcing under climate change, such as melting of sea ice in polar regions. Overall, it is expected that primary production will decrease in temperate to tropical oceanic regions and will increase at high latitudes, while there is uncertainty in the direction, magnitude and differences of changing primary production in shelf and coastal regions [30].
One of the major challenges in estimating primary production from remote-sensing observations lies in the assignment of the photosynthetic efficiency of phytoplankton cells [31][32][33]. Models based on ocean-colour remote-sensing observations typically use a relationship between phytoplankton biomass and Photosynthetic Active Radiation (PAR, 400-700 nm) to compute primary production [4,31,34,35]. One such relationship is the photosynthesis versus irradiance (P-I) curve, which can be represented by a variety of mathematical equations [36,37]. The initial slope (α B ) and the assimilation number (P B m ) of the P-I curve may vary with environmental conditions such as irradiance, temperature and nutrient concentrations, and the taxonomic composition and size structure of phytoplankton communities [33,[38][39][40][41][42][43]. One of the strategies for the assignment of photosynthetic parameters in the computation of primary production on a global scale is to assign parameters on the basis of ecological provinces of the ocean [2,16,31,44,45], allowing for variations in photosynthetic parameters with season and with province. This strategy was adopted in the present study, and an existing global database of P-I parameters [33] was extended to improve spatial and temporal coverage. The global P-I database was subsequently partitioned using Longhurst's geographical classification system of biomes and provinces [16]. The biogeographic classification is based on physical conditions that shape the structure and function of phytoplankton communities over large (basin) scales, and the supply of nutrients and the average irradiance within the surface mixed layer that impact the physiological capacity of phytoplankton cells [16,33]. Another challenge in the estimation of primary production from satellite observations lies in the requirement to specify the vertical structure in phytoplankton biomass, given that the satellite observations are confined to a finite surface layer that is much smaller than the euphotic depth and is not resolved with depth. To overcome this limitation, we have used a large database of in situ chlorophyll-a profiles to extrapolate ocean-colour remote-sensing observations of surface chlorophyll-a through the water column [2,44,45]. Seasonal means of P-I parameters and chlorophyll-a profile parameters were then used together with a 20-year time series of remotely-sensed chlorophyll-a concentrations and surface PAR to establish global primary production and its changes over the two decades. The results are discussed in the context of the sensitivity of computed primary production to potential changes in the photosynthetic response of phytoplankton cells under changing environmental conditions.

Surface Chlorophyll-a Data from Satellites
Surface chlorophyll-a concentrations at 9 km spatial resolution and monthly temporal resolution for the period of 1998-2018 were obtained from the European Space Agency (ESA) Ocean Colour Climate Change Initiative project (OC-CCIv4.1, https://esa-oceancolour-cci.org/). The dataset contains merged products of observations from the Sea-viewing Wide Field-of-View Sensor (SeaWiFS, 1997-2010), the Medium Resolution Imaging Spectrometer (MERIS, 2002(MERIS, -2012, the Moderate Resolution Imaging Spectroradiometer (MODIS, 2002-present) and the Visible Infrared Imaging Radiometer Suite (VIIRS, 2012-present) that are climate-quality controlled, bias-corrected and error-characterised (see details below) [10].

Primary Production Model
Several models have been described to estimate primary production based on ocean-colour remote-sensing observations [29,35,[46][47][48]. All models calculate daily water column production as a function of some measure of phytoplankton biomass and the photosynthetic response of phytoplankton to light. However, the different models can be categorised as linear or non-linear; spectral or non-spectral; vertically-uniform or vertically-non-uniform; or a combination of these [46,47]. They have also been categorised as depth-integrated or resolved and as wavelength-integrated or resolved [35]. Reducing models to a canonical form helps analyse similarities and differences between models and highlights the importance of model parameters [46,47,49]. The differences between spectral and non-spectral models are systematic and significant [47], but they can be corrected for [47,50,51]. In a study at the scale of the entire North Atlantic Ocean, Sathyendranath et al. [52] showed that ignoring the vertical structure in chlorophyll-a concentration reduced the computed primary production by about 9%, but in individual provinces, the difference could be higher (maximum reported value was about 16%). But the differences are systematic, and therefore, information on vertical structure should be taken into account when available. Furthermore, primary production within the deep chlorophyll maximum is likely fuelled by new production and would be important in calculations of new and export production [52].
In this study, we used a spectrally-resolved model that incorporates vertical structure in chlorophyll-a concentration [2,31,44,45,52], with recent updates [49]. This model simulates changes in photosynthesis as a function of irradiance using a two-parameter photosynthesis versus irradiance (P-I) function. The model has consistently performed well when compared with other models [5,35,53] and has been implemented on a global scale [2]. In the present study, considerable improvements have been made to the global coverage of the parameter database, while data provided by the Ocean Colour Climate Change Initiative (OC-CCI) project [10,54] allowed for the use of over 20 years of remote-sensing observations. The OC-CCI products [10] are multisensor products (reducing missing data), in which biases between sensors have been corrected (avoiding artificial trends in data arising from systematic differences between biases) and have been processed with a common protocol for calculation of chlorophyll-a concentration (minimising any systematic differences arising from differences between algorithms). Melin et al. [55] have shown that the chlorophyll-a trends calculated with OC-CCI time series are consistent with those calculated from single sensor products, demonstrating the fitness of the data for climate change studies. All of these, along with the length of the time-series data, are key considerations when studying the variability in ocean primary production in the context of climate change. The model used here is identical to the one described in Sathyendranath et al. [49] (see Appendix A for a brief description of the steps involved), but with a notable improvement to the P-I parameter assignment, based on an enhanced in situ database.

Photosynthesis Versus Irradiance Parameters
Parameters of the photosynthesis versus irradiance (P-I) curve were obtained from a global database [33,56] and additional literature sources . A quality check was performed on all data (9765 experiments) following Bouman et al. [33], using lower limits of the initial slope of the P-I curve α B (0.002 mg C mg Chl-a −1 h −1 (µmol photons m −2 s −1 ) −1 ) and the assimilation number of the P-I curve P B m (0.2 mg C mg Chl-a −1 h −1 ) and an upper limit for the maximum quantum yield of carbon fixation φ m (0.15 mol C mol photons −1 ). The value of φ m was calculated as α B /ā * B × 0.0231 [95] using α B and either simultaneous measurements of the mean specific chlorophyll-a absorption coefficientā * B (in m 2 mg Chl-a −1 ) or an estimate ofā * B based on Brewin et al. [96,97] (see Appendix A). In addition, major outliers in the dataset were identified using the outermost fences of the interquartile range. After quality control, 8676 experiments were used for further analysis. Note that this is a significant improvement over the P-I parameter database that was available at the time of Longhurst et al. [2]; they had access to 1862 P-I observations at that time, mostly from the North Atlantic Ocean. To estimate regional primary production, P-I data were assigned to biogeographic provinces according to Longhurst [16] (Appendix B). The P-I database covered 53 provinces, representing 96.6% of the world's oceans ( Figure 1). No in situ P-I experiments could be found for the coastal areas of Africa (EAFR) and India (INDE) and two regions in the North Pacific Ocean (NPPF, NPSE). The data were divided into seasons using 3-month intervals; i.e., March-May for spring/autumn, June-August for summer/winter, September-November for autumn/spring and December-February for winter/summer in the Northern/Southern Hemisphere. Mean and standard deviations of α B and P B m were calculated for each season and biogeographic province available in the P-I database ( Table 1). Temporal and spatial data gaps in α B and P B m were filled by statistical analysis of the relationships between seasons within each biogeographic province and the relationships between adjacent biogeographic provinces ( Figure 2). To this end, values of α B and P B m were log-transformed and significance (p < 0.05) was tested using ANOVA analysis followed by Tukey-Kramer post-hoc testing for unequal sample sizes (Past 3, Hammer et al. 2001). Results were used to assign mean and standard deviations of α B and P B m for missing seasons and/or biogeographic provinces ( Table 1) respecting boundaries of the ocean basins and biomes [16]. Linear-and log-scaled mean values of α B and P B m were highly similar (r 2 = 0.989, p < 0.001, with the majority of data normally distributed on regional and seasonal scales), and calculations of primary production were performed with the linear-scaled mean and standard deviation of each P-I parameter (Table 1) to support interpretation of (linear) trends in the sensitivity analyses (see below).  Relationships of photosynthesis versus irradiance (P-I) parameters between adjacent biogeographic provinces. Seasonal relationships are indicated by colour blocks, with significant differences (p < 0.05) denoted for the initial slope (α B ) and assimilation number (P B m ) of the P-I curve (red), α B only (dark blue) and P B m only (light blue). Comparisons were not available for the coastal areas of Africa (EAFR) and India (INDE), two regions in the North Pacific Ocean (NPPF, NPSE) and some seasons in other biogeographic regions due to lack of data (light grey). Biogeographic provinces are listed in Appendix B.

Analyses of Primary Production
The sensitivity of primary production to changes in photosynthetic parameters was estimated using three separate model runs for the period between 1998 and 2018, with three different P-I parameters assignments as follows: (a) with the mean of the P-I parameters for each of the provinces and seasons (the main run); (b) with the mean minus one standard deviation of the two P-I parameters (-1 SD); and (c) with the mean plus one standard deviation of the two P-I parameters (+1 SD) ( Table 1).
All other input variables (i.e., light, chlorophyll-a and chlorophyll-a profile parameters) were kept the same for each of these model runs. To assess the magnitude of change in primary production with change in individual P-I parameters, an additional sensitivity analysis was performed for a sample year (arbitrarily chosen to be 2003) in which either α B or P B m was adjusted by ±1 standard deviation, while the other parameter was maintained at its mean value. In the full model computations in which α B and P B m were adjusted simultaneously, the relationship between α B and P B m , given by the light adaptation parameter (I k in µmol photons m −2 s −1 ), remained largely unchanged. In the additional sensitivity analysis I k was allowed to vary and increased for −1 SD α B and +1 SD P B m and decreased for +1 SD α B and −1 SD P B m . For all model computations, both annual and seasonal primary production rates were estimated on global and regional scales. Regions were selected based on Longhurst's definitions of ocean basins (Antarctic, Atlantic, Indian and Pacific) and biomes (Coastal, Polar, Trades and Westerlies) [16]. Table 1. Sample size (n), mean and standard deviation (SD) of the initial slope (α B in mg C mg Chl-a −1 h −1 (µmol photons m −2 s −1 ) −1 ) and the assimilation number (P B m in mg C mg Chl-a −1 h −1 ) of the photosynthesis versus irradiance (P-I) curve for each biogeographic province and season. Values in blue are obtained from statistical comparisons between seasons and biogeographic provinces, while other values are directly obtained from the P-I database. Biogeographic provinces are listed in Appendix B. Statistical analyses were used to assess relationships between P-I parameters and other (environmental) variables available in the P-I database (latitude, depth, chlorophyll-a, PAR, temperature and nutrients; correlation and regression analysis), relationships between primary production estimates and input variables (chlorophyll-a, PAR and P-I parameters; correlation and regression analysis) and changes in primary production trends on various spatial and temporal scales (correlation analysis). To estimate the rate and direction of change in annual primary production between 1998 and 2018 for each grid point, monthly means were corrected for seasonality by subtracting monthly climatologies. The rate of change over time and its significance were calculated using linear regression and Student's t tests following Santer et al. [98]. Using the slope and intercept from the regression analysis, the percentage change per year in primary production (PP) was calculated as 100 · ((12 · slope) / (intercept + PP climatology )). Before all statistical procedures, data were tested for normality and homogeneity of variances and transformed for further statistical analysis when necessary. Differences were considered significant when p < 0.05.
The impacts of different climate indices on global and regional primary production were characterised based on annual mean anomalies that were corrected for seasonality following Racault et al. [19].

Global and Regional Annual Primary Production
Global annual primary production computed using mean photosynthesis versus irradiance (P-I) parameters (for each biogeographic province and for each season) varied from 38.8 to 42.1 Gt C y −1 in the period 1998-2018 (Table 2; Figures 3 and 4A). Summer (11.6-12.9 Gt C) was the most productive season in each of the years, followed by spring (11.0-11.8 Gt C), autumn (8.7-9.5 Gt C per season) and winter (7.5-8.0 Gt C per season) ( Figure 4B). On regional scales, annual primary production was highest in the Pacific Ocean (43.4-44.6%), followed by the Atlantic (27.7-28.8%), Indian (16.0-17.0%) and Antarctic Oceans (10.7-11.6%) (Table 2; Figure 3). In addition, the highest annual primary production rates were found at low latitudes in the Trades biome (39.6-41.1%), followed by the Westerlies (29.0-30.6%), Coastal (22.6-24.7%) and Polar biomes (5.9-6.8%) (Table 2; Figure 3). These regional differences in annual primary production were related to the surface areas of the specific ocean basins and biomes (r 2 = 0.643, p < 0.01), with the coastal regions being relatively more and polar regions relatively less productive than the other regions when computed as a rate per unit area (Table 2; Figure 3A).

Trends in Primary Production
Linear trends in annual primary production between 1998 and 2018 varied considerably on regional scales ( Figure 3B). At low and mid latitudes, trends in primary production were generally weak and negative (up to −3.0%), although large areas of positive trends were also observed in the South Atlantic Ocean and the South Pacific Ocean. In polar and coastal (upwelling) regions, stronger, positive trends in primary production were observed (up to +4.5%). Although significant linear trends were observed at individual pixels, the observed inter-annual changes in primary production on global and regional scales did not follow a linear pattern.
Inter-annual trends in global primary production showed an increase in rates between 1998 and 2003; relatively stable rates between 2003 and 2011; and a subsequent decrease in rates until 2015, after which rates showed a minor increase ( Figure 4A). Annual primary production in the Atlantic and Pacific Oceans showed similar inter-annual trends to global primary production (r 2 = 0.866, 0.926, p < 0.001) ( Figure 4C). Trends in annual primary production in the other ocean basins varied from the global trend, with relatively lower production between 2003 and 2011 in the Antarctic Ocean (r 2 = 0.675, p < 0.001) and a relatively early decrease in production in the Indian Ocean (r 2 = 0.828, p < 0.001) ( Figure 4C). Annual primary production in the Coastal, Trades and Westerlies biomes showed inter-annual trends comparable with that in global primary production (r 2 = 0.815-0.856, p < 0.001) with the highest rates observed between 1998 and 2000 in the Trades biome, and a relatively slow increase in production between 1998 and 2011 in the Westerlies biome ( Figure 4E). In the Polar biome, production decreased relatively early between 2004 and 2011 and was relatively high after 2015 compared with trends in global annual primary production (r 2 = 0.510, p < 0.001). Table 2. Climatological mean and standard deviation (n = 21) of annual primary production (in Gt C y −1 ) between 1998 and 2018 for each ocean basin and biome as defined by Longhurst (2007). Range in annual primary production between 1998 and 2018 is given in parenthesis. Results are given for primary production estimates based on mean, −1 standard deviation and +1 standard deviation photosynthesis versus irradiance (P-I) parameters. Surface areas (in km 2 ) for each ocean basin and biome are also provided. Trends in seasonal global primary production were highest in late spring to mid-summer, with the lowest rates observed in December for the Northern Hemisphere and in June for the Southern Hemisphere ( Figure 4B). Most regions showed similar seasonal trends in primary production with the peak occurring either earlier (Pacific Ocean and Westerlies and Coastal biomes) or later (Antarctic and Atlantic Oceans and Polar biome) in summer (r 2 = 0.782-0.962, p < 0.001) ( Figure 4D,F). Monthly primary production in the Trades biome was more variable from spring to autumn compared with the global trend (r 2 = 0.782, p < 0.01) ( Figure 4D). Trends in seasonal primary production in the Indian Ocean deviated most from the global trend, with two peaks in monthly primary production observed in spring and autumn and the lowest rates observed in summer ( Figure 4F).
Inter-annual and seasonal trends in global primary production were closely related to chlorophyll-a biomass (Spearman's rank correlation coefficient r s = 0.742-0.939, p < 0.05) ( Figure 3C). In the Antarctic and Indian Oceans and Trades biome, annual primary production was also related to Photosynthetic Active Radiation (PAR) (r s = 0.484-0.600, p < 0.05) ( Figure 3E Longhurst (2007). The dotted lines illustrate the relative global primary production per year (C,E) and month (D,F). Estimates of monthly primary production for the Southern Hemisphere were shifted to depict the summer season (December-February) along with that of the Northern Hemisphere (June-August) in months 6-8. Relative trends for each basin and biome were calculated by subtracting the minimum primary production from the annual (C,E) or monthly (D,F) primary production and dividing this by the difference between the minimum and maximum primary production between 1998 and 2018 or between January-December.

Sensitivity of Primary Primary Production to Changes in Photosynthetic Parameters
Global annual primary production varied from 20.4 to 22.2 Gt C y −1 between 1998 and 2018 when both P-I parameters were reduced simultaneously by one standard deviation (-1 SD), whereas the values ranged from 56.5 to 61.2 Gt C y −1 when the P-I parameters were increased by one standard deviation (+1 SD) (-46.5% and +44.9% compared with the results using the mean P-I estimates) (Table 2; Figures 3D,F and 5). The magnitude of the decrease in primary production when the P-I parameters were adjusted by −1 standard deviation was always greater than the increase in production when the P-I parameters were adjusted by +1 standard deviation ( Figure 5). The sensitivity of primary production to changes in P-I parameters was highest in the Atlantic Ocean, followed by the Pacific, Antarctic and Indian Oceans ( Figures 3D,F and 5; Table 2). The sensitivity was highest in the Trades biome and lowest in the Westerlies biome ( Figures 3D,F and 5; Table 2). Trends in global and regional annual primary production for the sensitivity analyses (data not shown) were similar to those observed for the main model run with mean P-I parameters ( Table 2; Figures 3B and 4) (r 2 = 0.978-0.999, p < 0.001).
On a seasonal basis, global primary production changed between −50.1 to −43.7% and +42.0 to +48.6% when the photosynthetic parameters were adjusted by −1 and +1 standard deviation, respectively ( Figure 5). The highest deviation from the mean P-I based primary production estimates was observed during spring and summer in the Atlantic Ocean and in the Trades biome, whereas the lowest deviation was observed during autumn and winter in the Antarctic Ocean and Westerlies biome. Trends in seasonal primary production were similar to those observed for the mean photosynthetic parameters estimates (Figure 4) when the photosynthetic parameters were adjusted by +1 standard deviation (data not shown). When the photosynthetic parameters were adjusted by −1 standard deviation, seasonal trends changed in the Indian Ocean and Coastal and Trades biomes. Primary production in these regions became relatively lower in spring and summer compared with other seasons (data not shown). No changes in seasonal primary production trends were observed in the Antarctic, Atlantic and Pacific Oceans and Polar and Westerlies biomes when photosynthetic parameters were adjusted by −1 standard deviation.   Figure 5. Percentage change in primary production (PP) for estimates based on mean photosynthesis versus irradiance (P-I) parameters ±1 standard deviation compared with estimates based on mean P-I parameters. Mean percentage differences in annual and seasonal primary production for each ocean basin and biome are given. Data were obtained from model computations in which both P-I parameters were adjusted simultaneously and the light adaptation parameter (I k ) was unchanged.

Relationship between Photosynthetic Parameters and Primary Production
It was expected that the changes in the magnitude of global and regional primary production were driven by variations in photosynthetic parameters, as all other input variables remained unchanged between the different model computations. When the relative change in primary production was compared with that of the P-I parameters for −1 SD and +1 SD estimates, variations were shown to be closely coupled (the light adaptation parameter I k was unchanged) (Figure 6). Both the initial slope of the P-I curve (α B ) (r 2 = 0.490 for −1 SD and r 2 = 0.508 for +1 SD estimates) and the assimilation number (P B m ) (r 2 = 0.750 for −1 SD and r 2 = 0.719 for +1 SD estimates) showed positive linear relationships with primary production for each season and biogeographic province. The weaker sensitivity of daily water column primary production to change in α B , relative to that of P B m , could be explained by the importance of α B under light-limited conditions, as opposed to P B m , whose effect is dominant in light-saturating conditions. It is important to note that the ratio of P B m to α B (i.e., I k ) remained unchanged between these different estimates of primary production. Independent variations in α B and P B m that modify I k could lead to higher sensitivity of primary production to the change [100][101][102][103]. The sensitivity analysis in which α B and P B m were independently adjusted by ±1 standard deviation (variable I k ) showed that changes in P B m caused greater variation in global annual primary production than changes in α B (Figure 7). Significant relationships between P-I parameters and primary production were also observed when α B and P B m were varied independently (-1 SD α B : y = 0.570 x, r 2 = 0.836; +1 SD α B : y = 0.322 x, r 2 = 0.440; −1 SD P B m : y = 0.741 x, r 2 = 0.908; +1 SD P B m : y = 0.492 x, r 2 = 0.733). When I k increased (-1 SD α B and +1 SD P B m ), primary production became more sensitive to changes in P B m compared with those in α B (see slope of relationships above).

Variation in Photosynthetic Parameters
In the global P-I parameter database α B ranged between 0.002-0.085 mg C mg Chl-a −1 h −1 (µmol photons m −2 s −1 ) −1 and P B m between 0.20-8.00 mg C mg Chl-a −1 h −1 . Mean values for each biogeographic province ranged between 0.005 and 0.054 mg C mg Chl-a −1 h −1 (µmol photons m −2 s −1 ) −1 for α B and between 1.01 and 6.25 mg C mg Chl-a −1 h −1 for P B m (Table 1). Lowest mean values of α B and P B m were observed in the Mediterranean (MEDI, summer) and Antarctic (ANTA, spring) provinces, whereas the highest values were observed in the Gulf Stream (GFST, winter) and Caribbean (CARB, summer) provinces, respectively (Table 1). Standard deviations varied between 0.2 and 99.1% (average of 43.8%) for α B and between 8.6 and 111.6% (average of 47.1%) for P B m ( Table 1). Similar to observations reported in Bouman et al. [33], spatial and temporal variations in photosynthetic parameters could be related to local environmental conditions. Relationships between α B and environmental conditions were variable between ocean basins and biomes resulting in relative weak relationships on a global scale (Figure 8). The initial slope α B increased with daily PAR in the Atlantic and Indian Oceans and with nitrate concentrations in the Antarctic, Atlantic and Pacific Oceans ( Figure 8). Positive relationships between α B and chlorophyll-a were observed at mid latitudes (Trades biome), but a negative relationship was observed in the Coastal biome. The standard deviation of α B increased at lower levels of PAR, at lower nitrate and silicate concentrations and at higher chlorophyll-a concentrations, but no other significant relationships with environmental parameters were observed ( Figure 8). The assimilation number P B m showed overall stronger relationships with environmental conditions compared with α B (Figure 8). Notably, P B m increased with temperature and PAR, possibly coinciding with latitudinal differences (Figure 8). The Pacific Ocean deviated from these results with an opposite trend observed between P B m and temperature at higher latitudes (data not shown). A negative relationship was observed between P B m and depth in all ocean basins and biomes, consistent with the known decline in P B m with decreasing temperature and light. P B m was generally lower at low phosphate concentrations, with the strongest relationships observed in the Antarctic, Atlantic and Indian Oceans at higher latitudes in the Polar and Coastal biomes. Variation in P B m as expressed by the standard deviation increased at higher temperatures and lower latitudes (Figure 8). The standard deviation of P B m also showed a positive relationship with depth and negative relationships with nutrient conditions.

Discussion
In the present study, a global database of photosynthesis versus irradiance (P-I) parameters, together with a 20-year time series of remote-sensing based chlorophyll-a concentrations, was used to study the magnitude and variability in marine primary production on a global scale. The estimate for global annual primary production of 38.8-42.1 Gt C y −1 between 1998 and 2018 in this study was within the range reported before (32.0-70.7 Gt C y −1 ) [5,104] and close to earlier reported values for depth-and wavelength-resolved primary production models (∼44 Gt C y −1 ) [2,4,5,7,22]. According to the model used in this study, primary production depends on phytoplankton biomass (in chlorophyll units), Photosynthetic Active Radiation (PAR, 400-700 nm; total value and its spectral and angular distribution) and on the assigned values of the photosynthetic and chlorophyll-a profile parameters. Although the model does not explicitly include the effects of environmental variables such as temperature and nutrients, or mixed-layer dynamics, these were implicitly accounted for through the photosynthetic and chlorophyll-a profile parameters which were assigned by season and biogeographic province [2,16]. Based on an inter-comparison of various primary production models, it has been reported that primary production generally increases at higher chlorophyll-a concentrations, higher PAR and shallower mixed-layer depths, whereas variability in temperature could either increase or decrease primary production [4]. In the present study, trends in global and regional annual primary production were best explained by variations in chlorophyll-a concentration, which in turn may vary with seasonal, inter-annual and multidecadal variations in physico-chemical conditions of the water column [17][18][19]. This study confirmed that global annual primary production varied with the ENSO and AMO [17][18][19]26], but not all variation in global annual primary production could be explained by large scale ocean-atmospheric oscillations. The previously reported negative (linear) trend in global annual primary production [25,27,28] was not observed in the present study. Instead, a more dynamic pattern of inter-annual trends in primary production was revealed at global and regional scales (also see [26,29]).
The assignment of photosynthetic parameters remains one of the major challenges in the assessment of global annual primary production using numerical models based on remote-sensing observations [31][32][33][34]. In this study, we have tackled this problem by assembling a database of around ten thousand observations that covered the majority of the biogeographic provinces of Longhurst [16]. The sensitivity of primary production to variations in the photosynthetic parameters was further studied by investigating the effect on primary production of changing P-I parameters from their mean values. P-I parameters may vary 2-10 fold among different biogeographic provinces (this study; [33,86,105]). This variation may reflect natural variability, but might also be affected to some extent by small differences in measurement protocols from author to author [33,86]. In the database used here, we tried to minimise the latter source of variability, for example by correcting values of the initial slope of the P-I curve (α B ) for the spectral quality of the light source used in the P-I experiment (also see [33]). A sensitivity analysis in which P-I parameters were adjusted by ±1 standard deviation revealed that the variation in photosynthetic rates may lead to a decrease or increase in the magnitude of global annual primary production by 45-47%. Global annual primary production remained within the range of earlier observations (32.0-70.7 Gt C y −1 ) [2,4,5,22] when both P-I parameters were adjusted by +1 standard deviation (+1 SD) (56.5-61.2 Gt C y −1 ), but adjustments by −1 standard deviation (-1 SD) resulted in considerably lower global annual primary production rates (20.4-22.2 Gt C y −1 ). Seasonal trends in global primary production were little affected, as the magnitude of change in P-I parameters was similar among seasons. The sensitivity analysis illustrated the importance of the parameters that describe the relationship between phytoplankton biomass and PAR in the calculations of primary production, but adjusting P-I parameters by ±1 standard deviation would represent the lower and upper limits of change in the photosynthetic response of phytoplankton cells. It would therefore be important to better understand the variability in P-I parameters and subsequent estimates of primary production under natural variations in environmental conditions and under global climate change.
Over the past three decades, considerable efforts have been made to establish a global database of P-I parameters ( [2,31,33,86]; this study) and to decipher their empirical relationships with physico-chemical and optical properties to enable prediction of photosynthetic parameters on regional and global scales [38,40,41,[105][106][107][108]. The observed relationships between physico-chemical conditions and P-I parameters in the present study confirmed earlier observations that temperature may be a good predictor of the assimilation number (P B m ), especially in coastal regions and temperate oceanic regions where temperature and associated water column stability dictates seasonal changes in the taxonomic and size structure of phytoplankton communities [40,86,109]. We note however, that the correlation coefficient between P B m and temperature is nowhere higher than 0.42, indicating that the importance of other factors (such as light and nutrient availability) in determining the variability in the assimilation number cannot be ruled out. The temperature dependence of P B m is of particular interest for assigning photosynthetic rates on regional and global scales, as Sea Surface Temperature (SST) can be obtained from remote-sensing observations on similar spatial and temporal scales to chlorophyll-a concentrations. Moreover, SST is a strong predictor of global climate change [24]. However, in regions with different underlying physical forcing that experience a smaller range in temperature, such as the Arabian Sea and open ocean gyres, the relationship between temperature and P B m is less obvious (this study; [40,41,105,107,110]). In such regions, chlorophyll-a concentration and the taxonomic and size structure of the phytoplankton community may be better indicators of variability in P B m [38,41,86,107,108]. The initial slope of the P-I curve seems to be more difficult to predict based on empirical relationships with physico-chemical conditions (this study, [41,86]), and it has been suggested that the simplest approach to estimate α B would be to relate α B to the assimilation number [33,110,111]. This approach may be supported by the strong dependence of I K (P B m /α B ) on latitude and depth, two spatial indicators that can be seen as general proxies of water column conditions [16,33,105].
The relationship between photosynthetic parameters and temperature is of particular interest in understanding the scope of change in primary production under global climate change. Over the past decades, SST has increased by 0.5 • C and is projected to increase a further 1.5-4.0 • C under different CO 2 emission scenarios [24]. The rise in SST and subsequent changes in stratification and nutrient loading into the euphotic zone are expected to affect phytoplankton growth and primary production [23,24]. One estimate of a potential change in annual primary production arising from variations in photosynthetic parameters under global climate change can be arrived at by using SST as the main driver of change in P B m . Assuming a simplified linear relationship between P B m and temperature in the Coastal biome (where temperature dependence of P B m was highest; P B m = 0.13 * T + 1.82, r 2 = 0.872 for T < 20 • C), P B m might be expected to increase by 8.3% under a rise of SST of +2 • C. Based on the relationships between P B m and primary production estimates presented in this study (Figure 7; assuming I K is unchanged), annual primary production in the Coastal biome could increase by +0.69 Gt C y −1 . Depending on the specific relationship with temperature, variations in P-I parameters and subsequent estimates of global primary production may vary on regional scales (for example, +13.4% in P B m in the Polar biome). The actual variation in P-I parameters and primary production under global climate change would be more complex and the interplay between different physico-chemical conditions will have a major effect on the direction of change.

Conclusions
It is the first time that highly quality-controlled, multisensor, inter-sensor-bias-corrected, ocean-colour observations extending over some two decades have been combined with increased spatial and temporal coverage of in situ observations of the photosynthetic parameters of phytoplankton, to compute the magnitude and variability of primary production on a global scale. This has led to a more accurate assessment of global annual primary production and its trends over the past 20 years. Variability in global annual primary production could be related to inter-annual and multidecadal oscillations, such that the present record of ocean-colour observations is not of sufficient length to detect trends associated with climate change [112]. Here, we report an inter-annual variability (standard deviation) of ±2.9% around a mean of 40.7 Gt C y −1 within the two decades studied. The importance of accurately assigning photosynthetic parameters in global and regional calculations of primary production has been illustrated by a sensitivity analysis. With the recent development of a global database of in situ measurements of P-I parameters [33] and the subsequent enhancement of the database (this study), photosynthetic parameters could be assigned to almost all biogeographic provinces (defined by Longhurst [16]). This considerably improved the confidence with which regional primary production could be estimated, especially in those regions that were previously known to be different from others, such as the Arabian Sea and the Antarctic Ocean [110]. Yet, the need to improve P-I data coverage in large areas of the global ocean still remains (this study, Figure 1; [33,49,86]). In particular, large areas of the Pacific and Indian Oceans remain poorly sampled. Methods designed to assign photosynthetic parameters based on their relationships to other variables amenable to remote-sensing [106,110], could, in the future, lead to a more dynamic assignment of these parameters. Sea surface temperature and phytoplankton community size structure (this study; [33,38,40,41,86,105,107,108]) could be suitable variables for further development of such methods for different ocean basins and biomes.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Model of Daily Water-Column Primary Production
Appendix A. 1

. Phytoplankton Biomass
In the model of Platt and Sathyendranath [31] and Sathyendranath et al. [49], ocean-colour remote-sensing products and a standard Gaussian function are used to calculate the distribution of phytoplankton biomass (B(z) in mg m −3 ) at depth. Depth profiles of chlorophyll-a were computed as shifted Gaussian functions: using the background biomass (B 0 in mg m −3 ), the total peak biomass (h in mg m −2 ), the depth of the chlorophyll maximum (z m in m), the standard deviation around the peak value (σ in mg m −3 ) and the ratio of the chlorophyll peak height to the total peak biomass at z m (ρ , dimensionless) [52]. Each profile parameter can vary independently, resulting in a versatile expression that can describe the biomass profile in a wide variety of oceanographic regimes. Profile parameters (h, σ and ρ ) were obtained for 57 biogeographic provinces [16] and four seasons from an archived global database of 26,232 in situ chlorophyll-a measurements [2,101]. At each pixel, the profile parameters were scaled such that the surface biomass matched the satellite chlorophyll-a value. Phytoplankton biomass profiles were then used for calculating the underwater light field and for estimating primary production (see below).

Appendix A.2. Irradiance Field
Spectrally-resolved irradiance at the sea surface was computed using a clear-sky model and expressed as the sum of a direct sunlight component and a diffuse skylight component. The surface direct and diffuse components were then scaled to match the daily Photosynthetic Active Radiation (PAR, 400-700 nm) products from the National Aeronautics and Space Administration (NASA) (https: //oceancolor.gsfc.nasa.gov/) and corrected for reflection and refraction at the sea surface assuming a flat ocean. The spectrally-resolved irradiance just below the surface was then used to construct the underwater light field (I(z, λ, θ) in µmol photons m −2 s −1 ), as the sum of a direct (d) and a diffuse (s) component of solar irradiance [113]: where z is the depth (in m), λ is the wavelength (in nm), θ d is the zenith angle of sun in water (in degrees), K d is the light attenuation coefficient (in m −1 ) for direct sunlight, K s is the light attenuation coefficient (in m −1 ) for diffuse skylight, a(z, λ) is the volume absorption coefficient (in m −1 ), b b (z, λ) is the backscattering coefficient (in m −1 ) and cos θ s is the mean cosine for the angular distribution of diffuse skylight after refraction at the sea surface [31,102]. The calculations make use of the chlorophyll-a profile to account for the influence of depth-dependent biomass on the light attenuation coefficient at each depth. The value of a(z, λ) is expressed as the sum of the contributions to absorption from pure seawater, phytoplankton, coloured dissolved organic matter and detritus. Absorption by phytoplankton depends on the concentration of chlorophyll-a (B(z)): where a * B (λ) is the chlorophyll-a specific absorption coefficient (in m 2 mg Chl-a −1 ) at wavelength λ. Concentrations of B(z) were obtained from ocean-colour remote-sensing observations (see main text) and expressed as the sum of chlorophyll-a concentrations contained in three size classes: pico-(p), nano-(n) and microphytoplankton (m). Phytoplankton absorption a B (z, λ) was estimated as the sum of the contributions of pico-, nano-and microphytoplankon to total phytoplankton absorption following Brewin et al. [96,97]: Here, the size-specific absorption coefficients a * p (λ), a * n (λ) and a * m (λ) (in m 2 mg Chl-a −1 ) are the values reported by Brewin et al. [96], and the fitted parameters B m p (0.13 mg Chl-a m −3 ) and B m p,n (0.77 mg Chl-a m −3 ) are the maximum concentrations attainable by picophytoplankton and combined pico-and nanophytoplankton, respectively. The parameters S p (6.15 m 3 mg Chl-a −1 ) and S p,n (1.26 m 3 mg Chl-a −1 ) determine the rate of change in the chlorophyll-a concentrations associated with picophytoplankton and the combined concentration of pico-and nanophytoplankton with changes in total chlorophyll-a concentration (model parameters are from Brewin et al. [97]). Similar to absorption, the total backscattering coefficient b b used in Equations (A4) and (A5) depends on back-scattering by pure seawater and by chlorophyll-a concentration, as in Sathyendranath et al. [114]: with b bw (z, λ) being the backscattering coefficient of water according to Morel [115] and b bB (z) the particle backscattering coefficient modelled as a function of chlorophyll-a concentration as in Sathyendranath et al. [114], following Ulloa et al. [116] and Loisel and Morel [117].

Appendix A.3. Daily Primary Production over the Water Column
The model of Platt and Sathyendranath [31], with the updates as in Sathyendranath et al. [49], uses a local algorithm based on surface biomass fields from ocean-colour remote-sensing, chlorophyll-a profile parameters, irradiance resolved with wavelength, angular distribution and depth, and photosynthesis versus irradiance (P-I) parameters to estimate water column primary production. Here, by the word "local", we imply that the model is implemented with parameters that are specific to the location and time. Primary production at depth z and time t (P B (z, t) in mg C mg Chl-a −1 h −1 ) is given by: where Π B (z, t) = where α B (z, t, λ) is the photosynthetic action spectrum (in mg C mg Chl-a −1 h −1 (µmol photons m −2 s −1 ) −1 ) and integrals are taken over the range of PAR (400-700 nm) [31,37]. In Equation (A10), the shape of α B (z, t, λ) is scaled such that the mean value is equal to the non-spectral value of α B for flat, white light [118] and the spectral shape of α B is taken to be the same as that of the phytoplankton absorption spectrum. Note that the P-I parameters do not change with depth in the present primary production model. Model calculations were performed at 9 km spatial resolution using a wavelength resolution of 5 nm, a depth interval of 0.5 m from the surface to the euphotic depth (depth at which light is reduced to 1% of its value at the surface) and at 12 time steps from dawn till local noon. The computed production at each depth and at each time step was summed over depth and time, and then multiplied by two to obtain daily water column primary production. In the event of any missing data in monthly OC-CCI chlorophyll-a fields, the computed primary production in each biogeographic province and in each month was scaled to full coverage using the mean primary production and the area of that province, with a weighting function accounting for variability in PAR in the specific biogeographic province. Mean monthly production (in mg C m −2 d −1 ) in each biogeographic province was then summed to obtain global annual primary production (in Gt C y −1 ) for each year between 1998 and 2018.