Drivers and uncertainties of future global marine primary production in marine ecosystem models

. Past model studies have projected a global decrease in marine net primary production (NPP) over the 21st century, but these studies focused on the multi-model mean rather than on the large inter-model differences. Here, we analyze model-simulated changes in NPP for the 21st century under IPCC’s high-emission scenario RCP8.5. We use a suite of nine coupled carbon–climate Earth system models with embedded marine ecosystem models and focus on the spread between the different models and the underlying reasons. Globally, NPP decreases in ﬁve out of the nine models over the course of the 21st century, while three show no signiﬁcant trend and one even simulates an increase. The largest model spread occurs in the low latitudes (between 30 ◦ S and 30 ◦ N), with individual models simulating relative changes between − 25 and + 40 %. Of the seven models diagnosing a net decrease in NPP in the low latitudes, only three simulate this to be a consequence of the classical interpretation, i.e., a stronger nutrient limitation due to increased stratiﬁcation leading to reduced phytoplankton growth. In the other four, warming-induced increases in phytoplankton growth outbal-ance the stronger nutrient limitation. However, temperature-driven increases in grazing and other loss processes cause a net decrease in phytoplankton biomass and reduce NPP despite higher growth


Introduction
By producing organic matter, marine phytoplankton form the base of the marine food web, control the amount of food available for higher trophic levels, and drive the majority of the ocean's biogeochemical cycles, particularly that of carbon.The net formation rate of organic carbon by phytoplankton, i.e., net primary production (NPP), is a key determinant for the export of organic carbon from the surface ocean, thereby governing how ocean biology impacts the ocean-atmosphere exchange of CO 2 (Falkowski et al., 2003;Sarmiento and Gruber, 2006).Accurate projections of future patterns of NPP may be crucial not only to estimate the potential impacts of climate change on marine ecosystems and fishery yields but also to properly assess the evolution of the ocean carbon sink under anthropogenic climate change.
Several authors have analyzed trends in future NPP and the underlying drivers, using models of strongly varying complexity and spatial resolution with regard to both the physical and the ecosystem components and also investigating different climate change scenarios.In the majority of these studies, global marine NPP was projected to decrease in response to future climate change (Bopp et al., 2001;Boyd and Doney, 2002;Steinacher et al., 2010;Bopp et al., 2013;Marinov et al., 2013;Cabré et al., 2014).The main mechanism suggested to explain this decrease in NPP was a decrease in the upward supply of nutrients in the low latitudes because of increased vertical stratification (Bopp et al., 2001;Steinacher et al., 2010) and reduced upwelling.Lower nutrient availability then resulted in a decrease in phytoplankton growth and therefore reduced NPP.
However, a few studies produced contradicting results, i.e., they reported increases in global NPP as climate change progresses over the 21st century (Sarmiento et al., 2004;Schmittner et al., 2008).Taucher and Oschlies (2011) showed that in the case of the model used by Schmittner et al. (2008), the simulated increase in NPP is caused by the warmer temperatures enhancing phytoplankton growth and overcoming the suppression of their growth owing to stronger nutrient stress.However, this result cannot be easily generalized, since some models used in Steinacher et al. (2010) still project a decrease in NPP even though they have a stronger temperature dependence of the growth rate than that used in the model by Schmittner et al. (2008).
The past century provides very little experimental constraint on the impact of long-term climate change on marine productivity, largely because of the lack of long-term (> 50 years) observations.Using a combination of in situ observations of chlorophyll and of ocean transparency, Boyce et al. (2010) suggested a substantial decrease in phytoplankton biomass over the last 50 years, implying a very strong response of phytoplankton to ocean warming.This result has been met with a lot of scepticism (e.g., Rykaczewski and Dunne, 2011), especially because an independent assessment of long-term trends in ocean color by Wernand et al. (2013) implied no overall global trend.Smaller decreases in NPP ( 6 % over 50 years) were suggested by a hindcast simulation, where a marine ecosystem model coupled to an ocean general circulation model was forced with observed atmospheric variability and changes over the last 50 years (Laufkötter et al., 2013).The satellite observations since late 1997 suggest a negative correlation between sea surface temperature and NPP (Behrenfeld et al., 2006), but the observation period is clearly too short to distinguish natural fluctuations from an anthropogenically driven trend in global marine NPP (Henson et al., 2011;Antoine et al., 2005;Gregg, 2003).
Far less work has been done regarding future trends in the biomass of specific plankton functional types (PFTs), despite their importance in shaping ecosystem structure and function (Le Quéré et al., 2005).Experiments have revealed a negative relationship between warmer waters and phytoplankton cell size, suggesting that future warming may tend to favor small phytoplankton (Morán et al., 2010).Moreover, using year-to-year variability associated with the North Atlantic Oscillation and the Southern Annular Mode, Alvain et al. (2013) found that more stagnant conditions and warmer temperatures tend to disfavor diatoms, suggesting that diatoms will become less prevalent in the future.The few modeling studies available support this view, i.e., they reported global decreases in the diatom fraction and a shift towards smaller size classes (Bopp et al., 2005;Marinov et al., 2010Marinov et al., , 2013;;Dutkiewicz et al., 2013).In these models, this shift was driven by increased nutrient limitation that affected diatoms more strongly than small phytoplankton.
While published studies emphasized the role of changes in bottom-up factors in explaining the changes in NPP, topdown control by zooplankton grazing may also drive future changes in total NPP or phytoplankton composition.This mechanism is intriguing, since top-down control was recently identified as one of the main drivers of phytoplankton competition during blooms in several ecosystem models (Hashioka et al., 2013;Prowe et al., 2011).Further, topdown control affects the onset of the spring bloom (Behrenfeld, 2010;Behrenfeld et al., 2013), influences primary production in a trait-based ecosystem model (Prowe et al., 2012) and affects NPP and export production changes on regional scales (Bopp et al., 2001).
Previous efforts in comparing different models with regard to future trends in NPP have analyzed the multi-model mean response and focused on identifying regions of consistent changes and mechanisms among models (Steinacher et al., 2010;Bopp et al., 2013;Cabré et al., 2014).By largely disregarding the regions of large inconsistencies, this focus may have underestimated the uncertainty associated with current projections of future marine NPP changes.This is well illustrated by the most recent model comparison study by Bopp et al. (2013), where the spread in the global NPP change between the 10 investigated global models for a given climate change scenario was larger ( 20 and +2 %) than the NPP difference between the different scenarios for the multimodel mean ( 9 to 2 %), demonstrating that the model uncertainty is larger than the scenario uncertainty.
Reasons why models differ are seldom investigated in model comparison studies.In particular, it is often not readily clear whether the large spread in model projections is mainly caused by differences in the underlying ocean circulation model, by differences in the complexity of the ecosystem models or by differences in the parameterizations leading to differing sensitivities to, e.g., changes in temperature, nutrients and light.Such information is needed, however, in order to improve the existing models and eventually obtain reliable future projections.
In this work we go beyond the basic analysis of the multimodel mean and the identification of regions of model consistency.Our aim is to identify where models differ and by how much and then determine why they do so, i.e., to identify the underlying drivers of change.To this end, we use results from a set of eight global marine ecosystem models coupled to or forced with nine coupled carbon-climate Earth system models, which have simulated the future evolution of marine NPP under the Intergovernmental Panel on Climate Change (IPCC) Representative Concentration Pathways (RCP) 8.5 (van Vuuren et al., 2011).We decompose the long-term changes in NPP into the contributions of the different phytoplankton functional types and then identify the relative importance and uncertainty of the main drivers.We demonstrate that (i) current marine ecosystem models are revealing more spread with regard to future changes in NPP than shown previously, and (ii) even where the models simulate consistent changes, the underlying drivers are quite different.In particular, we highlight the critical but not wellquantified role of temperature change in determining the future changes in NPP.

Model descriptions
We use projections for the 2012-2100 period of nine model simulations for the IPCC's RCP8.5 scenario from the "MARine Ecosystem Model Intercomparison Project" (MAREMIP, http://pft.ees.hokudai.ac.jp/maremip/ index.shtml,Vogt et al., 2013;Sailley et al., 2013;Hashioka et al., 2013) and/or the "Coupled Model Intercomparison Project 5" (CMIP5, Taylor et al., 2012).As we perform an analysis of the effect of PFT composition on NPP changes, we include only data from those models that possess at least two phytoplankton PFTs and at least one zooplankton PFT.For the models taken from the CMIP5 archive, only the first ensemble member (r1i1p1) was used.
These criteria led us to use data from eight different marine ecosystem models: diat-HadOCC, BEC, TOPAZ, PISCES, MEM, PELAGOS, REcoM2 and PlankTOM5.3(Table 1 lists the model acronyms, their main references, and further information, e.g., on spin-up times).Since the same ocean ecosystem model PISCES was used in two different Earth system models, we analyze a total of nine different simulations.In most simulations, the ecosystem model was embedded into a coupled climate model and integrated over thousands of years in order to spin up the model under preindustrial conditions (see Table 1).In two simulations (REcoM2 and Plank-TOM5.3), the ecosystem model was used within a forced ocean model and was initialized with observed climatologies.In these simulations, a control run showed considerably smaller drift than the climate change response.We do not correct the small drift in these models to keep the internal mechanisms in the models consistent.
We describe the most important features of the ecosystem models in the following and give the full equations and parameters for the offline calculations shown in this work in the Appendix.The ocean ecosystem models used in this study are structurally similar, but they differ substantially in their details (see Table for an overview of the model structures).Within our selection, all models simulate at least two phytoplankton PFTs, usually representing diatoms and a nanophytoplankton type, and one zooplankton PFT.BEC and TOPAZ have an additional diazotrophic phytoplankton PFT.Moreover, TOPAZ differentiates between diatoms and other large phytoplankton depending on the availability of silicic acid.In PELAGOS, the nanophytoplankton type is further divided into flagellates and picophytoplankton.Plank-TOM5.3 includes an explicit coccolithophore type, while in most other models coccolithophores are modeled implicitly as a fraction of nanophytoplankton.Regarding zooplankton PFTs, TOPAZ only has implicit zooplankton activity, diat-HadOCC, BEC, and REcoM2 have one zooplankton type, while PISCES and PlankTOM5.3differentiate between micro-and mesozooplankton.MEM and PELAGOS have three zooplankton types, i.e., in addition to the microand mesozooplankton, they include predatory zooplankton in MEM and heterotrophic flagellates in PELAGOS.Finally, PELAGOS is the only model that includes heterotrophic bacteria explicitly.

Analysis of NPP and its drivers
A change in NPP can be driven by (i) a change in the biomass-specific rate of photosynthesis, (ii) changes in autotrophic respiration, or (iii) changes in phytoplankton biomass through, e.g., zooplankton grazing, sinking and other loss processes of phytoplankton.However, only PELA-GOS and REcoM2 model photosynthesis (gross primary production, GPP) and autotrophic respiration separately.Rather, most models calculate NPP directly as the product of the growth rate µ and biomass of phytoplankton, P .In these latter models, changes in marine NPP can thus result only from (i) changes in the phytoplankton growth rate and (ii) changes in phytoplankton biomass.In order to disentangle these two main classes of drivers, it is helpful to consider the full mass balance equation for any phytoplankton type P i : where 0 is the sum of the time rate of change and the physical processes of advection, convection, and diffusion, and where the first term on the right-hand side is NPP.We consider any driver that alters the growth rate µ i as a bottomup driver, while those that alter P , i.e., grazing, sinking, and other losses, we consider as top-down drivers, even though only grazing is strictly speaking a top-down process.
In all models, the growth rate of phytoplankton is parameterized using a multiplicative function of a maximum growth rate µ max , the temperature limitation T P f and the nutrient and light limitation factors N lim ,L lim , i.e., In all eight models except for diat-HadOCC, the temperature dependence of phytoplankton growth, i.e., T P f , is described using an exponential function based on Eppley (1972), albeit with rather different temperature sensitivities (i.e., Q 10 values; see also Table 3).In diat-HadOCC, phytoplankton growth is independent of temperature.While in most models the same Q 10 value is used for all phytoplankton PFTs, mesozooplankton has a higher Q 10 in PISCES and PELA-GOS and each PFT and process has its own Q 10 value in PlankTOM5.3,derived from observations.In REcoM2 an Arrhenius function is used which results in a Q 10 that decreases with temperature.
The nutrient and light limitation factors have dimensionless values between 0 and 1, with higher values promoting higher growth.All models consider limitation by multiple nutrients, with six out of the eight models applying Liebig's law of the minimum (Liebig, 1840), such that the value of the strongest limiting nutrient sets the total nutrient limitation.Thus, these models do not consider nutrient co-limitation.Exceptions to this are PELAGOS and diat-HadOCC, where nutrient limitation is multiplicative.In all models, nanophytoplankton growth is limited by nitrate and iron, while diatoms are additionally limited by silicic acid.In several models, limitation with respect to phosphate and ammonia is additionally considered (see Table ).The limitation regarding a specific nutrient is calculated either with Michaelis-Menten functions (Michaelis and Menten, 1913), following optimal uptake kinetics (Smith et al., 2009), or using a cell quota representation of nutrient deficiency, often with strong differences in half-saturation constants.The values of the halfsaturation constants and the equations are given in the Appendix; Table lists the type of nutrient limitation for the different models.
For diat-HadOCC, the full model equations are not available; therefore we cannot describe the light limitation.In all other models light limitation is parameterized based on the work of Geider et al. (1998), Webb et al. (1974) and Platt et al. (1980).Most models (except for MEM) use the following equation: where the constant parameter ↵ denotes the initial slope of the photosynthesis-irradiance curve, ✓ chl : c is the chlorophyll-to-carbon ratio, PAR is the photosynthetically available radiation and µ T,N is the maximum growth rate multiplied with the temperature effect and nutrient limitation.PISCES models an additional strengthening in light limitation when the mixed layer depth is deeper than the euphotic zone.In PELAGOS, µ T,N is replaced by a constant p r for the maximum specific photosynthetic rate.TOPAZ replaces the instantaneous chlorophyll-to-carbon ratio with a variable ratio that depends additionally on the memory of irradiance over the scale of 24 h (see the Appendix).
MEM uses the light limitation function from Platt et al. (1980): where is a photo-inhibition index and ↵,p r , and PAR are as above.
Note that in most models, temperature and nutrient status influence also the light limitation, such that in addition to the direct effects of temperature and nutrients on the growth rate, there is an additional indirect effect through light limitation (Geider et al., 1998).
Since PELAGOS does not compute NPP directly and also uses a different formulation for the growth limitation terms, it requires a separate analysis: in this model, NPP is calculated for each phytoplankton type by subtracting autotrophic respiration and other loss processes from its GPP, i.e., NPP i = GPP i exudation i respiration i lysis i .GPP is calculated in PELAGOS analogously to how NPP is calculated in the other models, i.e., using the product of biomass, maximum growth rate, temperature, light limitation and iron and silicic acid limitation.Nitrate and phosphate limitation are accounted for in the phytoplankton exudation and lysis terms.The reason for this differentiation between the various limiting nutrients is to account for internal storage capabilities of the phytoplankton cells (Vichi et al., 2007).To be able to compare PELAGOS to all other ecosystem models within a common framework, we estimated a multiplicative nutrient limitation factor on the basis of temperature, light limitation and the growth rate that was given in the PELAGOS output: Regarding the loss terms for phytoplankton biomass, grazing is considered in all models.However, given the large diversity in the complexity and parameterizations associated with the modeling of zooplankton, the role of grazing may differ substantially among the considered models (Sailley et al., 2013).
Grazing of zooplankton Z on phytoplankton P is calculated as grazing(Z, P ) in all models except TOPAZ, where g P Z is the maximum grazing rate of zooplankton Z on phytoplankton P and T Z f is the temperature limitation of zooplankton feeding.TOPAZ simulates the effects of zooplankton implicitly, and the representation of grazing is based on Dunne et al. (2005).Most models employ the same temperature sensitivity for zooplankton as they use for phytoplankton, i.e., T Z f = T P f , with the exception of PISCES and PELAGOS, where the mesozooplankton has a higher temperature sensitivity, and PlankTOM5.3,where each PFT has a different Q 10 value.The food depen-  dence is modeled differently in each model and is shown in Table 5.

Data processing
Our analysis is based on monthly mean output for all surface ocean variables for the period 2012-2100.In order to facilitate direct comparisons, we regridded the model to a common 1 ⇥ 1 grid using the Earth System Modeling Framework (ESMF) regridding routines included in the NCAR Command Language (NCL) version 6.1.2,with the interpolation method set to bilinear.All models provided vertically (0-100 m) integrated NPP and biomass (in carbon units) of all PFTs.Primary production by diatoms and small phytoplankton was not available for PlankTOM5.3,MEM and PELAGOS and was estimated offline using the product of biomass and growth rate.The temperature limitations and growth rates were recalculated for all models except for PELAGOS and TOPAZ, where the growth rates were given in the model output.The nutrient and light limitation factors were included in the output of BEC, REcoM2 and TOPAZ, while they were recalculated from the monthly mean data for all other models using the original (not interpolated) data.The equations used for the recalculations are given in the Appendix.A comparison of recalculated and true values in the BEC model showed that the error in the recalculation is on the order of less than 10 %.
Changes for all properties are computed by first averaging the data for two 20-year periods, i.e., 2012-2031 and 2081-  Westberry et al. (2008) for NPP.We compare nutrients for the 1990-1999 period, while chlorophyll and NPP data are from 1997 to 2006.The angular coordinate shows the correlation coefficient, the distance from the origin denotes the normalized standard deviation and the distance from point [1.1] describes the root mean squared error.2100, and then taking the difference.For the growth limitation factors, we show the ratio changes, i.e., for any limitation factor x, we show the ratio <x>(t=2081 2100) <x>(t=2012 2031) , where the chevrons denote temporal averages.This is because the product of the relative changes in the temperature, light and nutrient limitation results approximately in the relative change in growth rate and the factor with the strongest change also has the strongest effect on the change in growth rate.

Model evaluation
Most of the models analyzed in this study have been evaluated individually in their respective documenting publications (see references in Table 1).Therefore, we restrict ourselves to an evaluation of the variables that are most relevant for this work, i.e., vertically integrated NPP, chlorophyll (chl), surface NO 3 , surface PO 4 and surface SiOH 4 (Fig. 1 and Tables 6, 7).We compare modeled NPP, using a 1998-2007 climatology for each model, with results from the updated Carbon-based Production Model-2 algorithm derived from Sea-viewing Wide Fieldof-view Sensor (SeaWiFS) satellite data (Westberry et al., 2008), downloaded from http://www.science.oregonstate.edu/ocean.productivity/index.php.For chlorophyll, we use chlorophyll a from the SeaWiFS Project generated by the NASA Goddard Space Flight Centre (ftp://oceans.gsfc.nasa.gov).We used monthly means computed from Level 3 binned daily products.For both NPP and chlorophyll data we removed coastal values (depth < 500 m) prior to the calculations.For the nutrients, we used the respective objectively analyzed climatologies from the World Ocean Atlas 2013 (Garcia et al., 2014) and compared it to model output for the 1990-1999 period.
On a global scale, the model-simulated nitrate fields correlate reasonably well with the observations, with all models showing correlations between 0.62 and 0.85 and normalized standard deviations (NSD) between 0.86 and 1.10.However, the bias is rather large, with values between 4.24 and +4.89 mmol N m 3 , corresponding to a bias of approximately ±70 % of the global average.For phosphate (not shown), the results are very similar to those of nitrate but for silicic acid the models perform less successfully.The correlations are lower and between 0.45 and 0.76, the normalized standard deviations scatter more, and the biases are larger (see Table 6).
The correlations for chlorophyll are mostly between 0.5 and 0.72; however, the normalized standard deviations are rather low (most models have NSD values < 0.5).The higher standard deviation in the observations stems mostly from the coastal ocean (standard deviation decreases from 1.8 Least well simulated is the distribution of NPP.The correlations are relatively low (0.18-0.69), the range of normalized standard deviation is as large as that of silicic acid (0.78 to 1.49), and in some of the models, the bias is very large ( 8.8 to +6.8 mol C m 2 yr 2 ).Global annual NPP ranges between 17 and 83 Pg C yr 1 (40.1 Pg C yr 1 in the multimodel mean), compared to on average 50.7 Pg C yr 1 in the satellite-based estimates (Carr et al., 2006) and 58 ± 7 based on 14 C NPP (Buitenhuis et al., 2013a).
However, global correlations in nutrients and NPP are strongly influenced by the globally dominant gradient between the Southern Ocean and the low latitudes.While this www.biogeosciences.net/12/6955/2015/Biogeosciences, 12, 6955-6984, 2015 gradient is generally well reproduced by the models, the model skill in reproducing the regional nutrient and NPP patterns is considerably lower (not shown).

Twenty-first-century changes in primary production
Starting from very different levels, the models simulate global NPP to change under the RCP8.5 scenario anywhere from 15 to +30 % ( 4.3 to +10 Pg C yr 1 ) over the 2012 to 2100 period (Fig. 2).One model shows an increase, five models show a decrease and three models project changes of less than 1 %, which are not significant (p value > 0.05) when compared to the level of interannual variability.The models suggest a median decrease of 7.2 % with an interquartile range (IQR) of 13.4 % ( 2 Pg C yr 1 with an IQR of 4.5 Pg C yr 1 ).This is comparable to the results reported by Bopp et al. (2013) using 10 Earth system models from the CMIP5 project under RCP8.5 ( 8.6 % ± 7.9 %) and also to another recent multi-model comparison conducted by Steinacher et al. (2010) under the IPCC's special report emission scenario A2 ( 10 ± 3 %, 2.9 ± 1.4 Gt C yr 1 ).However, the range of projections covered by our study with respect to NPP (45, 16 % without PlankTOM5.3)is higher than the 14 and 6 % reported by Bopp et al. (2013) and Steinacher et al. (2010), respectively.The regional pattern of the multi-model median change in NPP (Fig. 3b) shows distinct regional differences.The multimodel median suggests NPP increases in the Southern Ocean (south of 40 S, +10 %), in the Arctic Ocean (+40 %), in the southern Indian Ocean and in the southern subtropical Pacific, while decreases of 10.9 ± 23.5 % are projected for the low latitudes (30 S-30 N), with strongest decreases in the North Atlantic ( 30 %) and along the Equator in all basins.The range of NPP projections in different regions is given in Table 8.In most models as well as in the multi-model median, the decreases in the low latitudes are stronger than the increases in the high latitudes, resulting in the global decrease in NPP.This partial regional compensation was noted both by Bopp et al. (2013) and Steinacher et al. (2010).However, these changes are spatially heterogeneous and the multimodel mean masks differences between the individual models.
To illustrate these inter-model differences, we show the IQR (Fig. 3c) of the absolute change in NPP at each location.The IQR of NPP is around 1 mol C m 2 yr 1 in the high and intermediate latitudes, which is of the same magnitude as the trends in the multi-model median changes.In the low latitudes the IQR is significantly higher with values between 3 and 5 mol C m 2 yr 1 , exceeding the multi-model median substantially.Thus, the model projections lack consistency, making direct interpretation of the multi-model median response difficult.

Changes in bottom-up versus top-down control
The changes in NPP in the different models can be driven either by changes in the growth rates (bottom-up) or phytoplankton biomass (top-down control; see Sect.2.2 above).In order to obtain a first impression of the potential reasons underlying the NPP changes, we split the change in NPP into a component representing the change in the biomass of the whole phytoplankton community and a component representing the whole community growth rate.As the growth rates are only available at the surface in many models, we calculate the components for surface NPP changes.We computed these two components by first calculating a first-order Taylor decomposition of NPP into the changes in growth rate weighted with biomass and the changes in biomass weighted with growth rate within each model and for each phytoplankton PFT j : We then determine the median across all models (Fig. 4).We find that the multi-model median growth rates increase nearly everywhere, while the median biomass decreases in the low latitudes but increases in the Southern Ocean, mimicking the changes in NPP.As was the case for NPP, the model spread is large for both factors driving NPP, and particularly so in the low latitudes (not shown).
We focus next on the drivers affecting the growth rates, i.e., the bottom-up factors temperature, light, and nutrients, and afterwards discuss the factors affecting phytoplankton biomass, i.e., the top-down control, and do so from a global perspective.We then extend the analysis to the level of individual phytoplankton PFTs, which is best done on the regional scale, across which the responses are relatively homogeneous in contrast to the global scale.

Global analysis of bottom-up factors
Figure 5 shows the projected changes in sea surface temperature, photosynthetically active radiation (PAR) and surface concentrations of NO 3 and Fe as a zonal average for all models.Figure 6 shows the resulting relative changes in the growth rates and the limitation factors for temperature, light and nutrients for all models where the equations describing the limitation factors were available.Note that an increase in any limitation factor corresponds to an alleviation of this limitation, i.e., a positive impact on the growth rate.To simplify the plot, for each model only the values for the phytoplankton PFT with the strongest temperature (or light or nutrient lim-  In the low latitudes, sea surface temperature is projected to warm by about 2-3 C with some model variance (Fig. 5a).In the Southern Ocean, the warming is less pronounced and even more consistent among models (+1 ± 1 C), while in the Arctic Ocean, the warming is not only stronger but also differs strongly among the models (projections range between no change and +4 C).This surface ocean warming stimulates phytoplankton growth everywhere and in all models, although given the different temperature sensitivities and the different levels of warming, the spread is large (Fig. 6a).In the low latitudes, the temperature limitation factor is simulated to increase by +10 and +30 % (corresponding to weaker limitation).In the Southern Ocean the increase remains small (0-10 %), reflecting the small temperature changes, while in the northern high latitudes, the temperature limitation factor increases by up to 40 %.
In contrast to the large changes in temperature, the PAR at the surface changes little globally, with the important exception of the high latitudes (Fig. 5b), where light availability is affected by changes in sea ice.In the Arctic, PAR is modeled to increase (projections range between increases of 2 and 18 W m 2 ), while in the Southern Ocean, models disagree even on the direction of change, reflecting the divergent trends in sea ice (Mahlstein et al., 2013).Consequently, most models show little changes and also little spread in the surface light limitation term between 60 N and 60 S (Fig. 6b).In the high latitudes the spread is generally larger, with projections in light limitation factor ranging between 10 and +40 %.However, in all but one model, relative changes in light limitation are of a similar magnitude as the relative changes in temperature limitation in the high latitudes.
The iron concentrations are projected to change in a latitudinally relatively uniform manner with changes between 0.05 and +0.2 µmol Fe m 3 , with one exception (diat-HadOCC), where a strong increase is simulated (+0.5 µmol Fe m 3 ) in the Arctic.These generally small and uniform changes are reflecting the constant dust deposition in all models.Regionally, models differ most in the change in iron concentration in the equatorial Pacific (not shown), potentially related to differences in the transport of iron-rich water by the Equatorial Undercurrent (Vichi et al., 2011a;Ruggio et al., 2013).There is little agreement among the models with regard to the direction of changes in the surface concentration of nitrate, with decreases and increases of up to ±3 mmol N m 3 .Similar changes are modeled for phosphate (not shown).The large range of projected trends leads to very wide ranges for the relative changes in the nutrient limitation factor.In fact, with increases and decreases of up to ±90 % in the low latitudes, ±15 % in the Southern Ocean and 0 and 40 % in the region north of 30 N, the nutrient limitation factor is changing the most.
In nearly all models, the magnitude of the nutrient limitation term is determined solely by the most limiting nutrient (Liebig limitation, see Sect. 2).Except for PlankTOM5.3,the limitation patterns for different PFTs within the same model are rather similar, but the differences between models are large.Therefore, we show in Fig. 7 the limitation pattern only for diatoms.In the Southern Ocean, most models agree on iron-limiting phytoplankton growth in the annual mean, while PlankTOM5.3only simulates iron limitation in parts of the Southern Ocean and near the Antarctic continent in summer.In the low latitudes, models show substantial differences in the equatorial upwelling region in the Pacific.Only some models capture the iron limitation shown in data (C.M. Moore et al., 2013).There is substantial variation in the extent of the iron-limited region and also the direction of change in iron concentration.As this is a region with high NPP values in the annual mean (see Table 8), uncertainties in this region significantly affect the range in NPP projections.In the remaining low latitudes, models show either phosphate or nitrate limitation.
As half of the models use specified N : P Redfield ratios instead of modeling an explicit PO 4 tracer, nitrate and phosphate limitation cannot be distinguished in these models.However, as nitrate and phosphate are usually highly correlated, a differentiation between nitrate and phosphate limitation might not significantly increase the uncertainty in nutrient limitation projections.Most models agree on stronger nutrient limitation (a decrease in the nutrient limitation factor of between 0.01 and 0.05) in the low latitudes ex-  In summary, the changes in nutrients and temperature emerge as the most important determinants for the changes in the growth rates, with light generally playing a lesser role, except for the very high latitudes, particularly the Arctic.The changes in the bottom-up factors determine the changes in phytoplankton growth rate, which are shown in Fig. 6d; again, the PFT with the strongest changes is shown.In two models (CNRM-PISCES and IPSL-PISCES), the growth rate decreases in all latitudes except for the Southern Ocean, resulting from the high nutrient stress in these two models.However, all other models predominantly simulate increases in growth rates, owing to the temperature effect outweighing the decrease in nutrient availability in these models.The decreases in low-latitude NPP in the six models that show increases in growth rates are thus not bottom-up driven but caused by a loss of biomass.Future changes in phytoplankton biomass can be caused by several top-down factors, which we will discuss in the following.

Analysis of top-down control
Possible reasons for the simulated phytoplankton biomass decrease between 2000 and 2100 are (1) changes in circulation or mixing leading to a stronger lateral or vertical loss of biomass, (2) increased aggregation or mortality of phytoplankton if explicitly modeled or (3) a higher grazing pressure.Unfortunately, none of these fluxes have been stored by most models.Further, recalculated values are not precise enough to analyze the difference between NPP and loss processes.Therefore, we cannot quantitatively differentiate our analysis into the changes in grazing loss, aggregation and physical biomass loss across all models.We nevertheless try to shine some light on this critically important issue by using qualitative arguments and the partial information we have from those models that were able to provide the phytoplankton grazing loss.
We hypothesize that the loss of biomass caused by physical transport does not significantly increase, as all models show an increase in stratification over the next century.Furthermore, phytoplankton aggregation (and mortality) depends exponentially (linearly) on biomass but are temperature independent, so neither aggregation nor mortality losses can increase at lower biomass levels, eliminating this set of processes as well.This leaves us with increased grazing pres- sure as the most likely driver of the simulated biomass loss in the low and intermediate latitudes and the high northern latitudes.This hypothesis is supported by the fact that in all five models for which the grazing fluxes were available the fraction of grazed NPP increases throughout the 21st century north of 50 S (Fig. 8), i.e., the grazing pressure increases.In TOPAZ the increase is comparatively small (+0.1 %).However, grazing is the only loss process in this model and changes in the ratio between grazing and NPP have a direct and strong impact on phytoplankton biomass.In the models where aggregation and mortality are explicitly modeled, the increase in the grazed fraction of NPP is stronger (+5 ± 3 %).This larger change in the grazed fraction in these models can be understood when considering that as biomass decreases, the aggregation losses decrease as well.This automatically leads to a shift in the loss pathways toward grazing, even though the grazing pressure per se does not change.The fact, however, that biomass decreases in the first place, strongly indicates that increases in grazing pressure are the driver of the phytoplankton biomass losses diagnosed in the models.
To better understand the potential drivers of this increase in the grazing pressure, we analyze the fraction of NPP that is grazed by zooplankton, given by Here, g P Z is the grazing rate, T Z f and T P f are the temperature limitation for phytoplankton and zooplankton, respectively, P and Z denote phytoplankton and zooplankton biomass and µ max is the maximum phytoplankton growth rate, as introduced in Eqs. ( 2) and ( 5).P dependence is the dependence on prey concentration, as shown in Table 5.
Climate change affects the ratio between grazing and NPP via temperature and also via changes in nutrient and light limitation.Furthermore, the grazing : NPP ratio is affected by changes in zooplankton biomass, i.e., increases in total grazing and zooplankton mortality indirectly play a role.In the models where the same temperature function for both phytoplankton growth and zooplankton grazing is used (i.e., T P f = T Z f ; see Table 3), the temperature limitation cancels out.Still, with a higher temperature the total grazing intensifies due to an increase in zooplankton growth rate and thus a larger zooplankton biomass, which will intensify grazing (see Eq. 6).On the other hand, the grazing : NPP ratio can in-crease through a decrease in the phytoplankton growth rate µ because of stronger light or nutrient limitation, thus decreasing NPP in the equation above.
Lacking the model output from the three-dimensional models, we use a one-box model to explore the sensitivity of the grazing : NPP ratio to changes in temperature and nutrient limitation instead.To this end we consider only one phytoplankton and one zooplankton group, using a simplified form of the equations and parameterizations of the BEC model.We did not include further phytoplankton loss terms like aggregation or mortality and used a quadratic temperatureindependent mortality as loss process for zooplankton.We performed a spin-up until the model reached an equilibrium state under conditions representative of the low latitudes (temperature limitation of 0.8 corresponding to about 27 C, strong nutrient limitation of 0.1 corresponding to less than 0.5 mmol NO 3 m 3 and weak light limitation).As grazing is the only loss process of phytoplankton, 100 % of NPP is grazed in the equilibrium state.To test the sensitivity of grazing pressure to temperature changes, we increased the temperature from 27 to 30 C over a time period of 10 years but kept light and nutrient limitation constant.The experiment showed an 8 % decrease in phytoplankton biomass within the 10 simulation years even though the phytoplankton growth rate was increasing, caused by a temperature-driven increase in zooplankton biomass and thus grazing.On average, about 101 % of NPP was grazed per month during the 10-year period.
To test the sensitivity of grazing pressure to nutrient changes, we enhanced nutrient limitation by 30 % (nutrient limitation factor decreases from 0.1 to 0.07) over 10 years, while keeping temperature constant at 27 C.In this experiment, phytoplankton biomass decreased by 15 %.Besides the decrease in phytoplankton growth in this experiment compared to the equilibrium state and the first experiment, 102.5 % of NPP was grazed on average each month, indicating that the change in nutrient limitation has a similar effect on grazing to the temperature increase.These results indicate that the grazing pressure can be increased by both stronger nutrient limitation and higher temperatures.As the basic structure of the NPP and grazing equations is similar in most models, this mechanism might explain the observed biomass loss in the low and northern high latitudes.However, the specific grazing parameterizations and also the zooplankton mortality parameterizations differ substantially between models, such that the strength of the grazing response and the magnitude of the biomass loss is most likely different between models.

Regional changes in phytoplankton community structure
In the following we will refine our analysis to include differences in PFT responses and focus on two example regions: the low latitudes (30 S-30 N) and the Southern Ocean (50-90 S).The low latitudes have been chosen because they explain a large part of the global NPP change (Table 8).Moreover, they exhibit the largest interquartile range (Fig. 3c) and are therefore the main reason for the large range in global NPP projections.The Southern Ocean has been chosen to demonstrate the mechanisms underlying NPP changes for a region where NPP increases in the multi-model median.This is a also region where diatoms form a significant fraction of biomass and drive NPP changes in several models.The drivers of the NPP changes in the North Atlantic and North Pacific will be described briefly at the end of this section.

Low-latitude phytoplankton community changes
All models analyzed in this study except one agree that NPP decreases between the 2012-2031 and the 2081-2100 periods in the low latitudes (Fig. 9), albeit with different magnitudes (between 0.004 and 0.09 mol C m 3 yr 1 ).The exception is PlankTOM5.3that shows a strong increase of, on average, 0.1 mol C m 3 yr 1 .This will be discussed separately below.In three models the trend is caused by similar decreases in both diatom and nanophytoplankton NPP (BEC, TOPAZ and diat-HadOCC).Diatom changes contribute about a third of total NPP changes in both PISCES simulations and the decrease is mainly driven by a decrease in the NPP by nano-or picophytoplankton in PELAGOS and MEM, with little changes in diatom NPP.In REcoM2, diatoms and nanophytoplankton trends almost fully compensate for each other.Changes in diazotrophs (modeled in BEC and TOPAZ) and large non-diatom phytoplankton contribute less than 10 % to the total trend.Figure 10 shows the relative change in temperature, light and nutrient limitation, growth rate, biomass and NPP for diatoms, nano-or picophytoplankton and coccolithophores in the low latitudes.Diat-HadOCC could not be included in the figure as the equations for the limitation factors are not available.In three models (BEC, MEM, REcoM2) diatoms show a stronger response to nutrient limitation than nanophytoplankton, which translates into a smaller increase or even decrease in growth.However, in all models except MEM diatoms show larger relative biomass and NPP losses than nanophytoplankton, indicating that in TOPAZ, PELA-GOS and the PISCES simulations top-down control is the main reason for the decrease in diatom relative contribution to biomass.
The PlankTOM5.3 trend is caused by an increase in coccolithophore NPP (+0.14 mol C m 3 yr 1 ), partly compensated for by a decrease in nanophytoplankton NPP ( 0.04 mol C m 3 yr 1 ).We note that export production changes do not follow the increase in NPP but decrease strongly (not shown), indicating a very large increase in the recycling efficiency in this model.This is caused by a strong increase in microzooplankton biomass and their grazing on phytoplankton, with rapid recycling of the nutrients back to their inorganic forms explaining the increase in nutrient .Relative change in temperature limitation factor (red), light limitation factor (yellow), nutrient limitation factor (orange), growth rate (green), biomass (light blue) and NPP (purple) for nanophytoplankton (full), diatoms (hatched) and coccolithophores (dotted) at the surface of the low latitudes for all models where the full equations were available.An increase in limitation factor denotes weaker limitation, which leads to stronger growth.The relative change in a variable is the ratio between the 2081-2100 average and the 2012-2031 average.A value of 1 means no change, 1.2 corresponds to a 20 % increase, 0.8 corresponds to a 20 % decrease.The product of the relative change in temperature, light and nutrient limitation results approximately in the relative change in growth rate.See main text for further details.
availability.This greatly enhances regenerated production, even as new production decreases.

Southern Ocean phytoplankton community changes
All models simulate an increase in surface NPP in the Southern Ocean south of 50 S, but the magnitude of the change varies by several orders (+0.006 and +0.11 mol C m 3 yr 1 ; Fig. 11).Furthermore, the contributions of the different phytoplankton PFTs to these NPP trends differ strongly between the different models.Four models show a stronger increase in the NPP by nanophytoplankton compared to that by diatoms (PlankTOM5.3,MEM, TOPAZ, CNRM-PISCES), three models show an exclusively diatom-driven NPP change (BEC, PELAGOS, REcoM2) and two models show similar changes in the NPP by diatoms and nanophytoplankton (diat-HadOCC, IPSL-PISCES).Only one model shows a significant decrease in diatom NPP (PlankTOM5.3).
In seven out of eight models, surface ocean warming is the most important driver of the increase in phytoplankton growth for both diatoms and nanophytoplankton.All but the CNRM-PISCES and PELAGOS model show a relief from nutrient stress for all phytoplankton types, i.e., an increase in the nutrient limitation factor (1-15 % increase), although these models remain iron limited throughout the 21st century.Diatoms respond more strongly to changes in nutrient concentrations than nanophytoplankton in all models except for PlankTOM5.3.In addition, in many models a stronger topdown control of nanophytoplankton than diatoms becomes apparent, indicated by differences in biomass changes despite similar growth rate changes.Only in MEM and Plank-TOM5.3 do diatoms seem to be more strongly top-down controlled.In PELAGOS the diatom fraction is almost 100 % south of 50 S and shows little changes.The final result is a stronger increase in diatom NPP compared to nanophytoplankton NPP in BEC, TOPAZ, IPSL, CNRM and REcoM2 and a weaker increase or even decrease in diatom NPP in MEM and PlankTOM5.3.

NPP changes and their drivers
Our finding of the key role of temperature in defining the response of NPP to future climate change contrasts with the conclusion of the majority of the past studies, which attributed the decrease in NPP to a decrease in nutrient availability, particularly in the low latitudes (Bopp et al., 2001;Moore et al., 2002;Steinacher et al., 2010;Marinov et al., 2010;Cabré et al., 2014).To explain this discrepancy, we focus on the temperature functions in the models used in the studies above.First, several of the earlier models had either no temperature dependence of phytoplankton growth at all (the model HAMOCC5.1) or the temperature sensitivity was rather weak, with a Q 10 value of 1.13 for temperatures higher than 15 C (in the models HAMOCC3 and NCAR CSM1.4carbon).It is thus not surprising that ocean warming did not significantly affect global productivity in these model simulations compared to the models analyzed in this study that have a Q 10 of at least 1.68.A further model analyzed in several of these studies is the IPSL model with PISCES as the ecological or biogeochemical component.A later version of this model is analyzed in our study.Like previous authors we find that changes in nutrient limitation are the main driver of NPP changes in PISCES, mostly because decreases in nutrient limitation are significantly stronger compared to other models.Finally, several authors attribute projected decreases in NPP to concurrent decreases in macronutrient availability (Steinacher et al., 2010;Marinov et al., 2013;Cabré et al., 2014).However, our analysis shows that in many models the global NPP decrease, and particularly that in the low latitudes, is not caused by decreasing growth rates, such as one would expect from increasing nutrient limitation.Rather the decrease in NPP is caused by biomass losses, presumably a result of a warming-induced increase in grazing pressure.We conclude that the temperature effect and top-down control might have been underestimated in several earlier studies.The importance of warming which we have identified for future NPP is more in line with another group of studies, where global NPP was projected to increase with climate change and a temperature-driven increase in metabolic rates was identified as the cause (Schmittner and Galbraith, 2008;Sarmiento et al., 2004;Taucher and Oschlies, 2011).However, Schmittner and Galbraith (2008) and Taucher and Oschlies (2011) considered only the temperature dependence of phytoplankton growth and remineralization, while the growth of zooplankton and hence the grazing pressure on phytoplankton were independent of temperature.Likewise, the algorithm used to estimate chlorophyll in Sarmiento et al. (2004) is based on the assumption that chlorophyll is purely bottom-up controlled.
Finally, Dutkiewicz et al. (2013) aimed to separate the direct temperature effect from the altered nutrient input and light availability caused by stratification.In their study, temperature, nutrient and light changes compensate for each other nearly perfectly, resulting in very little change in global NPP.Still, the importance of temperature for phytoplankton growth and zooplankton grazing shown by most models in our study indicates that temperature might play a major role in the response of NPP to climate change.
While we emphasize here the role of temperature in the models, our understanding of how temperature controls the most important ecological and biogeochemical processes in real marine ecosystems is not well established.Most models base their parameterizations of temperature effects on laboratory studies that show -within favorable thermal rangesan exponential increase in growth with increasing temperature (Eppley, 1972;Bissinger et al., 2008).However, there are major uncertainties in quantifying the temperature sensitivities of different physiological processes and of functional types (Ikeda et al., 2001;Lomas et al., 2002;Hirst and Bunker, 2003;Hancke and Glud, 2004;Sand-Jensen et al., 2007).Several authors suggest a stronger temperature response of heterotrophs than autotrophs (López-Urrutia et al., 2006;Rose and Caron, 2007), which would lead to major consequences for the metabolic balance of the oceans under rising temperatures (Duarte et al., 2013;Williams et al., 2013;Ducklow and Doney, 2013;García-Corral et al., 2014).Furthermore, in current implementations both phyto-and zooplankton grow faster with increasing temperatures without any upper thermal limit beyond which growth rates may come down.The underlying assumption is that if the temperature rises to values outside the optimal range of a certain species, the species will be replaced by another species with a higher temperature tolerance.However, particularly in the tropics, it is unclear if this assumption holds.Thomas et al. (2012) showed that warming might lead to a decrease in diversity in the tropics, which could potentially lower NPP due to the loss of highly productive species.Finally, due to the lack of measurements, synergistic effects of multiple stressors are barely considered in current models.Recently, temperature sensitivity has been shown to be reduced under nutrient limitation (Staehr and Sand-Jensen, 2006;Tadonléké, 2010;Marãnón et al., 2014), which would result in an overestimation of temperature sensitivity and therefore NPP in the oligotrophic regions of the ocean.Overall, the temperature assumptions on which current model projections are based are affected by high uncertainties.

Changes in phytoplankton community
Seven out of nine models in our study show a global decrease in the relative abundance of diatoms with decreases in low latitudes but increases in the Southern Ocean, confirming results reported by Bopp et al. (2005); Marinov et al. (2010); Dutkiewicz et al. (2013); Manizza et al. (2010) and Marinov et al. (2013).The difference between the diatom and nanophytoplankton nutrient response has been identified as the primary driver of the decrease in diatom fraction in Bopp et al. (2005); Marinov et al. (2010) and Marinov et al. (2013).Our results show that, while models currently agree on a global decrease in diatom fraction, there is no agreement on regional changes and models do not agree on whether a stronger nutrient response or a higher susceptibility to grazing pressure is the cause.As diatom biomass tends to be overestimated by several of the models (Vogt et al., 2013;Hashioka et al., 2013), the relative importance of changes in diatom biomass may constitute an upper bound for future global NPP changes.

Identifying and reducing uncertainties
The spread in globally integrated NPP projections in our study is 45 %, with the PlankTOM5.3model causing 25 % of it alone.Given this wide spread in NPP projections, we attempt to identify the different sources of uncertainty in the following and then investigate whether there is a way to narrow the uncertainty of the projections using emergent constraints.

Sources of model uncertainties
The biogeochemical and biological parameterizations that contribute the largest uncertainties are initial nutrient concentrations: models (except Plank-TOM5.3) agree on similar decreases in nutrient concentration in the low latitudes, despite disparities with regard to the identification of the most limiting nutrient.However, the differences in relative nutrient limitation change are very large (±90 %, see Fig. 6).Particularly the PISCES simulations show a strong relative decrease in nutrient limitation, which is caused by it having nutrient concentrations that are (too) low at the beginning of the simulation (see Sect. 3).On the other hand, a positive bias in nutrients as observed in other models might lead to too weak a response in nutrient limitation.
the relative importance of iron versus nitrate limitation and projections for iron concentrations: increases in iron availability allow the small global increase in nanophytoplankton NPP in REcoM2 and attenuate or even outbalance the low-latitude NPP decrease in BEC and TOPAZ.This is of particular relevance in the equatorial upwelling region in the Pacific (see Fig. 3), which is iron limited according to observations (C.M. Moore et al., 2013) and is responsible for 14-33 % of global NPP at present in the different models (Table 8).The differences in the projected changes in iron concentration in the equatorial upwelling region in the Pacific are potentially related to differences in circulation: according to Vichi et al. (2011b) and Ruggio et al. (2013), the Equatorial Undercurrent may intensify and shoal with climate change and this may bring more iron to the eastern equatorial upwelling, partly off-setting the reduced nutrient input due to the warming surface.Note that the dust deposition is held constant in current projections.Variable iron forcing in future simulations will lead to more realistic NPP projections but might further increase this uncertainty.
different Q 10 values (between 1.68 and 2.08) and different projections for SST (sea surface temperature) increase (+2, +3 C), which together result in a high uncertainty in the temperature response of both phytoplankton growth and zooplankton grazing.Further un-certainty is introduced by the stronger temperature response of zooplankton types parameterized in some models.
the relative importance of the response of top-down controls versus that of the microbial loop, potentially related to different Q 10 values and differences in the partitioning of the grazed material.
the fact that there is no agreement with regard to the direction of change in light limitation in the Southern Ocean, reflecting the wide range in projected sea-ice changes and other factors influencing surface light such as cloud cover.However, light limitation currently only introduces a minor uncertainty compared to the nutrient and temperature effects, at least for surface NPP.
In order to reduce the spread in NPP projections, we also need to understand how much of the uncertainty arises from the underlying physical forcing (transport, mixing, and temperature) and how much is caused by the different ecosystem parameterizations.Unfortunately, the design of our study does not allow for a clear distinction between uncertainty from physical forcing and from the use of different ecosystems, as this would require us to compare projections from different Earth system models using the same ecosystem model with projections of one Earth system model coupled to different ecosystem models.Results from Sinha et al. (2010), who compared two different circulation models coupled to the same biogeochemical model, indicate that differences in the underlying physics lead to substantial differences in PFT biogeography but only have small effects on total NPP.A more ambitious program is currently being undertaken, in which a larger group of ecosystem models are being coupled to the same circulation model (the iMarNet project, Kwiatkowski et al., 2014).The outcome of this project will help to better separate the ecosystem model uncertainty from the uncertainty introduced by different physical models.

Constraining NPP projections
The concept of emergent constraints (e.g., Allen and Ingram, 2002) has been used with success to reduce uncertainties for future projections.The basic premise is that models that provide a better fit to a specific set of current constraints provide a better estimate for the future changes.The emergent constraint is usually established by finding a good correlation between an observable parameter for the present and the future change in NPP.We have tested for correlations between the different models' skill to predict current NPP and its projected changes, using both the 2012-2031 average of globally integrated NPP and the slope between chlorophyll a and sea surface temperature as a measure for model skill.Although chlorophyll a is a poor indicator of biomass in the low latitudes (Siegel et al., 2013), it can be used as indicator of model skill and is comparatively well constrained by obser- Figure 12.Relative changes between 2012-2032 and 2080-2100 in annual mean temperature limitation factor (red), light limitation factor (yellow), nutrient limitation factor (orange), growth rate (green), biomass (light blue) and NPP (purple) for nanophytoplankton (full), diatoms (hatched) and coccolithophores (dotted) at the surface of the Southern Ocean (50-90 S).An increase in limitation factor denotes weaker limitation, which leads to stronger growth.PELA-GOS has a relative diatom contribution of more than 95 % of total biomass; therefore we show only results for diatoms.
vations.As the metric for the projected changes we used the change in NPP defined as the difference between the 2012-2031 average and the 2081-2100 average and the change in NPP weighted with the temperature increase.Moreover, as regions with positive and negative changes might cancel each other out, leading to little net NPP change despite strong local changes, we also tested for a relation between the mean of absolute NPP changes and model skill.
We did not find any significant correlation between model skill and NPP changes, neither on regional nor on global scales, and the relation is weak at best between globally integrated NPP and the absolute change in NPP (Fig. 13).We hypothesize that the cause for this lack of relationship is the uncertainty in the relative importance of the net effect of temperature on NPP and on nutrient limitation.This hypothesis is supported by results from Taucher and Oschlies (2011), who compared two simulations, one temperature dependent and one independent.Both simulations fitted equally well to observations, but the direction of NPP change was the opposite.It seems that matching the current observations is not sufficient in order to estimate which sign of future NPP change is more realistic.Thus, we need a better understanding of the mechanisms in order to reduce the uncertainty in projections.Efforts to extend the amount of data that is available for model parameterization and evaluation (Buitenhuis et al., 2013b) will hopefully help achieve that goal.show the observed NPP range as reported by Carr et al. (2006); for the low latitudes we give the observed NPP range spanned by the estimates of Behrenfeld and Falkowski (1997) and Westberry et al. (2008).

Conclusions
In this work we present a multi-model comparison of nine model simulations with regard to NPP and its underlying drivers.We show projected changes in global NPP between 15 and +30 % by the end of this century for the highemission scenario RCP8.5, with the largest inter-model discrepancies stemming from the low latitudes.All but one model simulate either a decrease in NPP or changes of less than 0.5 % in the low latitudes but for very different reasons.The main drivers are warming-induced enhancement of phytoplankton growth, increased nutrient limitation and decreases in phytoplankton biomass, which are most likely caused by temperature-enhanced grazing by zooplankton.Only three models show reduced phytoplankton growth rates due to increased nutrient limitation.Thus, in this set of models, temperature and nutrient concentrations are at least equally important drivers for changes in NPP, contradicting many prior studies that emphasized the sole importance of a stronger nutrient limitation.
One major difficulty faced in this study is the limited availability of model output variables related to ecosystem growth and loss rates, particularly limitation factors and grazing rates.The changes in growth rate, temperature limitation, and light and nutrient limitation reported in this work have been recalculated in six out of nine models using surface monthly mean fields.The obtained results are therefore an approximation of the original values.We have compared recalculated values with original values in the models where the limitation factors were given, and we estimate the error to be less than 10 %.We conclude that, while the absolute values reported might be inaccurate, the relative importance of nutrient vs. temperature limitation shown in this work is correct.Furthermore, we can discuss only surface NPP changes.For the models where three-dimensional limitation factors were available (BEC, REcoM2), we compared our results for the surface with the 100 m average, and we can confirm that the same mechanisms that govern the surface changes also hold for the 100 m average.In addition, the changes in surface NPP correlate with the changes in integrated NPP in all models, except for the Arctic Ocean.It therefore seems likely that our surface drivers also describe the changes in integrated NPP.To ease future studies of NPP changes, we recommend the inclusion of mixed layer averages of growth rate, light and nutrient limitation and grazing fluxes in the standard model output of future model intercomparison projects.The availability of changes in growth rates could prevent common misinterpretations of drivers by analyzing univariate correlations with only one of several possible drivers.
To reduce the uncertainty in NPP projections, the representation of present-day nutrient concentrations and resulting limitation patterns should be further improved.Particularly a bias in present-day nutrient concentration strongly affects relative changes in nutrient limitation and therefore NPP projections.Furthermore, given the importance of top-down control shown in this work, we need a better understanding of zooplankton mortality and further potential drivers of zooplankton biomass like phenological or trophic mismatches, diseases or changes in predation from higher trophic levels.Finally, a better understanding of the temperature dependency of all key ecological or biogeochemical processes is needed.In particular, this includes the determination of the different temperature response functions for the different PFTs and trophic levels.

A1 BEC
In the following, we give the equations and parameters governing NPP in all models except for diat-HadOCC, where the full equations are currently not available.We abbreviate nanophytoplankton with "nano", diatoms with "diat", zooplankton with "zoo", and meso-and microzooplankton with "meso" and "micro", respectively.

Growth rate of phytoplankton PFT i
Temperature function (for all PFTs) Total nutrient limitation

Fe
Phosphate limitation of PFT i Silicate limitation of diatoms Nitrate and ammonium limitation of PFT i ) Light limitation of PFT i Grazing Diatoms of large phytoplankton depends on the silicate concentration Iron limitation of PFT i Nitrate and ammonium limitation of PFT i Silicate limitation of diatoms I Mem is the memory of irradiance over the scale of 24 h and was provided in the model output.

A3 PISCES
Growth rate of phytoplankton PFT i

Total nutrient limitation
Nitrate and ammonium limitation of PFT i where MXL denotes the mixed layer depth and Heup the depth of the euphotic zone.

Microzooplankton grazing
The food options I for mesozooplankton are nanophytoplankton, diatoms and microzooplankton.

A4 MEM
Growth rate of phytoplankton PFT i Temperature function (for all PFTs) Silicate limitation of diatoms

Nitrate and ammonium limitation of PFT
Light limitation of PFT i G meso!nano ,G meso!diat ,G meso!micro ,G pred!diat are all calculated using the same equation but different parameters.
Nutrient limitation with respect to phosphate and nitrate is not included in the phytoplankton growth rate but acts through the exudation and lysis terms.The exudation and lysis terms have not been recalculated in this work, instead we estimated a multiplicative nutrient limitation factor (see Sect. 1).We refer to Vichi et al. (2007) for a full description of the nutrient limitation in PELAGOS.

Temperature function for PFT
T 10 10 Light limitation of PFT i

Grazing
Grazing of zooplankton type i on phytoplankton type j is calculated as where F denotes the total food available and is calculated as e i j denotes the capture efficiency of zooplankton i when grazing on phytoplankton j .e i j is set to 1.0 for mesozooplankton.For microzooplankton and heterotrophic flagellates, it depends on prey density:  The food sources F for microzooplankton are small phytoplankton, diatoms, coccolithophores and small particulate organic carbon.The food sources F for mesozooplankton are small phytoplankton, diatoms, coccolithophores and small particulate organic carbon.

A7 REcoM2
Growth rate of phytoplankton PFT i

Figure 3 .
Figure 3. Spatial patterns of multi-model annual mean integrated net primary production (NPP) for (a) the 2012-2031 average, (b) changes between 2081-2100 and 2012-2031 under RCP8.5 and (c) interquartile range of the changes in NPP projections.The unit is mol C m 2 yr 1 .The blue boxes in (a) mark the regions which are discussed in more detail in this work, namely the Southern Ocean south of 50 S, the low latitudes (30 S-30 N) and the equatorial upwelling region in the Pacific.

Figure 4 .
Figure 4. First-order Taylor decomposition of the surface NPP changes between 2012-2032 and 2080-2100 in (a) biomassweighted changes in growth and (b) growth-weighted changes in biomass.The unit is mol C m 3 yr 1 .

Figure 5 .
Figure 5. Zonal mean of projected sea surface temperature change, photosynthetically active radiation (PAR) change and change in surface Fe and NO 3 concentrations.We calculate the change as the difference between the 2012-2031 average and the 2081-2100 average.Different line colors denote different models, as in the legend of Fig. 2.

Figure 6 .
Figure 6.Zonal mean of the relative change in temperature, nutrient and light limitation and growth rate.For each model only the values for the phytoplankton PFT with the strongest temperature (or light or nutrient limitation factor or growth rate) response is shown.We calculate the relative change as 2081-2100 average 2012-2031 average .A value of 1 means no change and is indicated by the dotted line.Different line colors denote different models, as in the legend of Fig. 2.

Figure 7 .
Figure 7. Changes in relative diatom nutrient limitation (calculated as the 2081-2100 average divided by the 2012-2031 average) in all models that use Liebig limitation (smallest individual nutrient limitation term determines total nutrient limitation).The colors indicate changes in the nutrient limitation value, with positive values indicating an increase in nutrient limitation factor which is equivalent to lower nutrient limitation and an increase in growth.The hatching indicates the limiting nutrient.A change in limiting nutrient during the simulation period is shown with dots.REcoM2 does not simulate the Arctic; these missing values are shown in white.

Figure 8 .
Figure8.Fraction of NPP that is grazed (grazing / NPP) normalized to the 2012-2031 average at the surface of the low and high northern latitudes ( 50 S-90 N).This plot shows data from all models where total grazing on phytoplankton is available in the output.

Figure 9 .
Figure 9. Decomposition of annual mean area-averaged low-latitude surface NPP changes between 2012-2032 and 2080-2100 (red bar, in mol C m 2 yr 1 ) into change in nanophytoplankton (yellow) and diatom (orange) surface NPP.Changes in diazotrophs (green) and picophytoplankton (light blue) have been included in the bar indicating nanophytoplankton changes for the models that simulate these functional types.For TOPAZ, changes in large non-diatom phytoplankton (dark blue) are included in the bar indicating diatom changes.Changes in coccolithophore NPP are shown in purple.Note the change in scale between the first three plots (models with large surface NPP changes) and the remaining six plots.While for diat-HadOCC, BEC, IPSL-PISCES, CNRM-PISCES, REcoM2 and TOPAZ the surface NPP of the PFTs was included in the model output, we show recalculated values for PlankTOM5.3,MEM and PELAGOS.

Figure 11 .
Figure 11.Decomposition of Southern Ocean (50-90 S) surface NPP changes between 2012-2032 and 2080-2100 (red bar, in mol C m 2 yr 1 ) into change in nanophytoplankton (yellow) and diatom (orange) surface NPP.Changes in diazotrophs (green) and picophytoplankton (light blue) have been included in the bar indicating nanophytoplankton changes for the models that simulate these functional types.For TOPAZ, changes in large non-diatom phytoplankton (dark blue) are included in the bar indicating diatom changes.Changes in coccolithophore NPP are shown in purple.Note the change in scale between the first three plots (models with large surface NPP changes) and the remaining six plots.While for diat-HadOCC, BEC, IPSL-PISCES, CNRM-PISCES, REcoM2 and TOPAZ, the surface NPP of the PFTs was included in the model output, we show recalculated values for PlankTOM5.3,MEM and PELAGOS.

Figure 13 .
Figure 13.Relationship between the change in NPP and the 2012-2031 average NPP for all models.Change in NPP has been calculated as the sum of the differences between the 2012-2031 average and the 2081-2100 average for each grid cell (open dots).We additionally show the negative absolute differences of the changes (full dots), calculated by taking the sum of the negative absolute differences between the 2012-2031 average and the 2081-2100 average for each grid cell.Each color represents a model, (a) shows global values and (b) shows the low latitudes.The gray area marks the range of current observational NPP estimates.For global values we show the observed NPP range as reported byCarr et al. (2006); for the low latitudes we give the observed NPP range spanned by the estimates ofBehrenfeld and Falkowski (1997) andWestberry et al. (2008).

I
denotes the food options and consists of diatoms and nanophytoplankton for microzooplankton.Grazing on diatoms is calculated accordingly.C. Laufkötter et al.: Drivers of future marine primary production Mesozooplankton grazing µ micro,flagellates .

Table 1 .
Vichi et al. (2011a)3)lations used in this work.For differences between the two PISCES simulations, seeSéférian et al. (2013). 2 Land and ocean carbon pools have been adjusted to the atmospheric preindustrial CO 2 with an acceleration method described inVichi et al. (2011a).

Table 3 .
Comparison of temperature limitations in ecosystem models.Meso stands for for mesozooplankton, micro for microzooplankton, cocco for coccolithohores, and nano for nanophytoplankton.

Table 4 .
Comparison of nutrient limitation of phytoplankton growth in ecosystem models.

Table 5 .
Comparison of prey dependence of grazing.For the full equations, see the Appendix.

Table 6 .
Westberry et al. (2008)ting global NPP, measured in Spearman's rank correlation, normalized standard deviation (NSD) and bias.The NPP data are fromWestberry et al. (2008); the average global NPP value is 12.6 mol C m 2 yr 1 .The chlorophyll data are from the SeaWiFS Project; the average global chlorophyll value is 0.28 mg Chl m 3 .
to 0.5 mg Chl m 3 when removing coastal areas with water depths < 500 m).Most models capture the lower open-ocean variability, however, in the two models that have a variability comparable to the observations (diat-HadOCC and Plank-TOM.3),the variability arises from the open ocean and is therefore significantly higher than the observed open-ocean variability.

Table 8 .
NPP changes (total and in percentage of the global value) in different regions.The Pacific upwelling region is shown in Fig. 3. Changes describe the difference between the 2012-2031 and the 2081-2100 average.

Table A1 .
Symbols used in the model equations.
P min G large = min(k graz max ,u max ⇥ T G nano = min(k graz max ,u max ⇥ T f ⇥ f ⇥ {N graz large }) ⇥ P large

Nitrate and ammonium limitation of PFT i
i ⇥ Z micro

Table A5 .
MEM parameters.Threshold value for meso.zoo.grazing on micro zoo.
i ⇥ Z meso

Table A7 .
PlankTOM5.3 parameters.Chl W d) 1 Initial slope of P I curve ↵ nano 0.83e 6 mol C m 2 (g Chl W d) 1 Initial slope of P I curve ↵ cocco 1.25e 6 mol C m 2 (g Chl W d) 1 Initial slope of P I curve