Long-term trends in ocean plankton production and particle export between 1960 – 2006

We analyse long-term trends in marine primary and particle export production and their link to marine phytoplankton community composition for the period 1960–2006 using a hindcast simulation of the Biogeochemical Elemental Cycling Model coupled to the ocean component of the Community Climate System Model. In our simulation, global primary and export production decrease significantly over the last 50 yr, by 6.5% and 8% respectively. These changes are associated with an 8.5% decrease in small phytoplankton biomass and 5% decrease in zooplankton biomass. Diatom biomass decreases globally by 3%, but with strong temporal and spatial variability. The strongest decreases in primary and export production occur in the western Pacific, where enhanced stratification leads to stronger nutrient limitation and a decrease in total phytoplankton. The concurrent decrease in diatom fraction and in zooplankton biomass causes a lower export efficiency in this region. Substantial phytoplankton composition changes also occur in the Southern Ocean and North Atlantic, although these are masked in part by a high degree of interannual variability. In these regions, stronger wind stress enhances mixing, reducing the biomass of small phytoplankton, while diatoms profit from higher nutrient inputs and lower grazing pressure. The relative fraction of diatoms correlates positively with the export efficiency ( r = 0.8, p r = –0.5, p < 0.05). However, the long-term trends in global export efficiency are ultimately driven by the reduction in small phytoplankton and particularly decreases in coccolithophore biomass. The diagnosed trends point toward a substantial sensitivity of marine primary production and export to climatic variations and trends.


Introduction
In the recent decades, evidence for impacts of climate variability and change on the global oceans has been accumulating (Denman et al., 2007).The ocean surface has warmed by 0.2 • C per decade in the last 30 yr (Levitus et al., 2009), ocean pH has decreased globally by 0.1 units compared to preindustrial times (Feely et al., 2009) and salinity changes have occurred associated with the amplification of the hydrological cycle (Durack et al., 2012).Warming has tended to increase ocean stratification globally, despite the observation of a significant increase in salinity in the low latitudes over the same period (Durack et al., 2012), which acts to decrease stratification there.At the same time, the trend toward lower salinities observed in the high latitudes actually enhances the warming-induced stratification there.This overall strengthening of the vertical density gradient has acted to reduce the mixing of surface waters with nutrient-rich deeper waters, resulting in stronger nutrient limitations for phytoplankton growth.This is expected to have resulted in a longterm spatial expansion of the oligotrophic areas of the ocean, similar to what was observed over the period 1998 to 2006 in the Atlantic and Pacific oceans (Polovina et al., 2008).These physical and chemical changes might have already affected marine ecosystems and marine productivity in a substantial manner (Boyd, 2011;Gruber, 2011), but very little is known so far about how marine plankton have changed in the last 50 years.
Marine plankton are an important component of the global carbon cycle.The autotrophic component, i.e, phytoplankton, take up inorganic carbon in the surface ocean and convert it to biologically fixed carbon during photosynthesis.Some of the fixed carbon is lost to dissolved organic carbon (DOC)

Laufkötter et al.: PFTs and export production 1960-2006
by respiration, by exudation or upon cell lysis.Zooplankton feeding upon phytoplankton, and the subsequent excretion and aggregation, form most of the particles that sink into the ocean interior.Together with the downward advection and mixing of DOC, this process is known as carbon export production or the soft-tissue component of the biological pump (Volk and Hoffert, 1985).Phytoplankton are also the main driver for the hard-tissue (carbonate) component of the biological pump, consisting of the formation of CaCO 3 by coccolithophores, and the sinking of this mineral phase into the ocean's interior upon death of the phytoplankton (Sarmiento and Gruber, 2006).
Plankton are taxonomically very diverse (Falkowski et al., 2004, e.g.,), but it is neither really conceivable nor effective to fully represent this diversity when investigating how marine plankton may have responded to past climate changes, and how they might respond to the changes looming ahead.In several current ocean models, the diversity is aggregated based on some criteria, such as the function different plankton types have in the oceanic cycling of biogeochemically important elements (Iglesias-Rodríguez, 2002).The resulting concept of plankton functional groups (PFTs) has been adopted widely in the modelling of marine plankton (Le Quéré et al., 2005) as it permits the representation of different pathways of carbon and other biogeochemically important elements through the marine food web and into the ocean's interior.In particular, such PFT-based models permit the investigation of how changes in plankton community composition might have changed the export fraction, i.e., the fraction of primary production that is exported into the deep ocean.This is critical when one wants to understand and predict the impact of changes in ocean biology on the oceanic uptake of CO 2 from the atmosphere, as it is only the export fraction that causes a net uptake of CO 2 from the atmosphere.
Estimating the response of the biological pump to future anthropogenic disturbances, such as global warming or ocean acidification, has proven to be challenging (Doney, 1999;Bopp et al., 2013).Nevertheless, a consensus has emerged that in nutrient-limited areas with little vertical mixing and shallow mixed layers, warming-induced stratification will lead to stronger nutrient limitation.Therefore, productivity and consequently export production is projected to decrease until 2100 (Bopp et al., 2001;Steinacher et al., 2009;Bopp et al., 2013).Satellite-based observations over the 1997-2006 period support this negative relationship between productivity and stratification (Behrenfeld et al., 2006) in the low latitudes.There is less agreement for the high latitudes, although most model simulations suggest that in light-limited areas with strong mixing, increased stratification results in a reduction in light limitation and thereby enhanced productivity and export production (Bopp et al., 2001;Doney et al., 2006;Bopp et al., 2013).It would be highly desirable to assess these projected changes for the 21st century with knowledge about the response of primary and export production to climatic changes over the past few decades, but observations prior to 1997 are scarce and very few model-based assessments have been conducted so far.
The few studies published so far on observed changes in marine ecosystems include the expansion of warmwater species into intermediate waters in the North Atlantic (Barnard et al., 2004;Beaugrand et al., 2002) and the alteration of plankton community structure in some regions, such as the Humboldt Current, the North Sea and the northeast Atlantic (Alheit and Niquen, 2004;Beaugrand, 2004;Richardson and Schoeman, 2004).The only global long-term studies conducted so far are those of Boyce et al. (2010) and Wernand et al. (2013), who estimated trends in global chlorophyll concentrations for the last century.Boyce et al. (2010) combined in situ chlorophyll measurements with indirect measurements of the ocean's transparency to extend the record in time and space.They suggested an overall substantial decline in marine phytoplankton biomass of about 1 % of the global median per year over the last 50 yr, albeit with substantial regional differences.These results have been controversially discussed (Rykaczewski and Dunne, 2011;Mackas, 2011;McQuatters-Gollop et al., 2011), with problems raised including the very large uncertainties associated with the merger of the different data sets and the contradiction with existing plankton biomass time series.Wernand et al. (2013) analysed a database of ocean colour, the Forel-Ule scale record.They report increases and decreases in different ocean basins over the last century, but no overall global trend.Thus, it remains inconclusive whether and how the observed changes in the chemical and physical environments of the ocean over the past few decades have altered global marine plankton production.
The current knowledge regarding changes in plankton community structure is even more limited.Researchers have attempted to estimate the present distribution of PFTs (Alvain et al., 2008;Hirata et al., 2011;Buitenhuis et al., 2013) and to reproduce their distribution in climate models (Aumont, 2003;Moore et al., 2002;Kishi et al., 2007;Gregg and Casey, 2007;Buitenhuis et al., 2010).But no systematic investigation of trends in community composition and its impact on production and export has been conducted so far.Modelling studies have shown that stronger stratification will favor small phytoplankton species over larger phytoplankton such as diatoms and that the resulting decrease in diatom fraction might lead to a less efficient export of organic matter (Marinov et al., 2010;Bopp et al., 2005), but this has not been observationally confirmed.
Here, we analyse the simulated trends in export production, primary production and plankton community structure and their drivers on regional-and global scales.We focus less on the absolute magnitude of the changes, but rather on the relationship between changes in community structure, the efficiency of export production and the changes in how carbon is routed from PFT biomass to sinking particles.To this end, we use a hindcast simulation with a comprehensive three-dimensional coupled physical-biogeochemical-ecological ocean model for the 1960-2006 period.The hindcast simulation was forced with prescribed atmospheric boundary conditions for wind and fluxes of heat and freshwater, resulting in an oceanic environment that varies through time in a manner very similar to the real ocean over these decades.Our study differs in this respect from those that employed fully coupled oceanatmosphere climate models, as the simulated temporal evolution of their oceanic environment could be drastically different from the real one due to the stochastic nature of the internal climate variability.
Our aim is to determine and analyze trends for the 1960 through 2006 period, but we do not make any attempt to attribute these changes to anthropogenic forcing.The described trend over these nearly five decades could be due to anthropogenic climate change, but they equally could be part of a natural multi-decadal oscillation of the climate system.Regardless of whether the trends are due to natural or anthropogenic processes, they are indicative of how marine plankton responds to perturbations and hence help us to ultimately better understand and predict the future.

Model description
We use hindcast simulations from the Biogeochemical Elemental Cycling (BEC) model (Moore et al., 2002(Moore et al., , 2004;;Doney et al., 2006) embedded in the ocean component (Parallel Ocean Program) of the global climate model CCSM3 (Collins et al., 2006).The spatial resolution is 3.6 degrees in longitude and 0.8 to 1.8 degrees in latitude, with a finer resolution around the equator (Yeager et al., 2006).The ocean component of the CCSM3 has 25 vertical levels with increasing thickness from approximately 12 to 450 m.
The ecosystem model BEC models one zooplankton functional type and three phytoplankton types: diatoms, diazotrophs and a small phytoplankton class, which represents nano-and picoplankton and includes an implicit representation of calcifiers.Nutrient limitation with regard to PO 3− 4 (phosphate), Fe (iron), Si(OH) 4 (silicic acid), NH + 4 (ammonia) and NO − 3 (nitrate) is calculated according to a Michaelis-Menten nutrient uptake kinetics.The total nutrient limitation is then set as the minimum limitation factor of all nutrients.Nutrients and light are co-limiting, with light limitation being calculated according to Geider et al. (1998).Temperature sensitivity is calculated with a Q 10 temperature function (Eppley, 1972) and is identical for all plankton types.
Diazotrophs constitute only a very small fraction of total phytoplankton biomass and net primary production in the BEC (Moore et al., 2002) and also do not contribute to the export of particulate organic matter.Since they do not constitute an important contributor to neither trends in net pri-mary production nor to trends in export production, we do not discuss this PFT any further.
The zooplankton type is parameterized to represent either micro-or mesozooplankton, i.e., its functionality switches, depending on the type of prey ingested.In particular, the zooplankton changes its growth rate and the routing of grazed organic matter to the detrital pools is weighted differently (Moore et al., 2004).
Two detrital pools are modeled, one representing dissolved organic matter (DOM) and the other particulate organic matter (POM).The latter is mainly produced through zooplankton grazing, but also by the aggregation of diatoms and small phytoplankton.Aggregation losses of phytoplankton are calculated using a quadratic function of biomass.Particles formed as a result of grazing are calculated as a fraction of the grazed matter, with diatoms contributing more strongly to particle production than small phytoplankton.The fraction of grazed matter is influenced by temperature, while aggregation depends only on biomass.
Sinking and remineralization of particulate organic matter is modeled according to the mineral ballast model by Armstrong et al. (2002) with an implicit scheme, i.e., the particles are assumed to sink instantaneously and are remineralized within the same water column as they had been formed.Phytoplankton types influence the sinking behaviour of organic matter through the production of SiO 2 (diatoms) and CaCO 3 (small phytoplankton).Those minerals lead to ballasted particles with different dissolution length scales.
Because Moore et al. (2004) and Doney et al. (2006) provide a full description of the equations and parameters of the model, we only list the equations describing phytoplankton growth, zooplankton grazing, and the production and sinking of particles in the Appendix.This permits us to better connect our results with the particular implementation of the respective processes in the model.Moreover, some parameters of the BEC were modified relative to those used by Doney et al. (2006), with a list with the differences given in the Appendix (Table 1).These changes were incorporated to improve the model results relative to observations, but did not lead to major changes in ecosystem structure.

Forcing and spinup
A two-step procedure was used to generate our hindcast simulation.First a 3000 yr preindustrial spin-up simulation forced at the surface with CORE CNYF v2 (Common Ocean-Ice Reference Experiments Corrected Normal-year Forcing; Large and Yeager, 2004) was conducted.This long spin-up resulted in a model with a negligibly small drift in surface nutrient concentrations and primary and export production.It also resulted in stable deep ocean radiocarbon values and a minimal net air-sea flux of CO 2 (Graven et al., 2012).Second, in 1950, the  and we ran the model forward in time through the end of 2006.This transient forcing data set was calculated using the NCEP reanalysis data set (Kalnay et al., 1996) and satellitebased estimates of radiation, sea surface temperature, sea-ice concentration and precipitation (Large and Yeager, 2009) and hence represents the full suite of physical forcings affecting the ocean.
Although the two CORE forcings were constructed in order to minimize the model perturbation when switching the forcing from the normal year forcing to the interannually varying one, we nevertheless do not consider the first 10 yr of the transient simulation.This results in 47 yr of model data for analysis, i.e., from January 1960 through December 2006.To ensure the absence of drifts over this period, we conducted also a control simulation forced with the normal year forcing (CORE CNYF v2) over the same period.As was the case with the spinup, plankton biomass and export exhibited a negligibly small drift, with changes of less than 0.1 % over these 47 yr.
The CORE forcings include annual mean river runoff, but we considered this flux in the freshwater forcing only, but not in the input of nutrients.There is also no atmospheric deposition of macro-nutrients in this version of BEC.We use a constant climatology to prescribe atmospheric iron fluxes on the basis of the data from Mahowald et al. (2009).

Biological fluxes
We determine the export production of particles (mol POC m −2 yr −1 ) as the average POC flux through the 100 m depth level for each year.NPP (mol POC m −2 yr −1 ) is given as the annual mean vertically integrated net primary production between the surface and 100 m depth.Phytoplankton and zooplankton biomass are calculated as the annual mean average concentration between surface and 100 m depth (mol C m −2 ).Export efficiency is the fraction of NPP that is exported through the 100 m depth level, i.e., particle export production divided by NPP.Diatom fraction is the fraction of diatom biomass compared to total phytoplankton biomass, i.e. diatom fraction = diatoms/(diatoms + small phytoplankton + diazotrophs).

Model evaluation
An extensive evaluation of the BEC coupled to the CCSM was done by Doney et al. (2009b).The evaluation focused on time-mean spatial patterns, seasonal cycle and interannual variability of several ecological variables, among others chlorophyll, NPP and surface nutrients.The data sets used for comparison have been taken from Conkright et al. (2001) for surface nutrients and from McClain et al. (2004) for Sea-WIFS surface chlorophyll.Integrated primary production has been calculated using SeaWIFS data and the VGPM model (Behrenfeld and Falkowski, 1997).We give a brief summary of the evaluation results and refer to Doney et al. (2009b) for more details.In contrast to our simulation they used the CORE-IAF version 1 forcing (Large and Yeager, 2004) and also employed a slightly different set of parameters (see Appendix, Table 1).However, these differences caused minor changes in ecosystem structure and do not alter the conclusion about the major strengths and weaknesses of this model simulation.Doney et al. (2009b) showed that while model-simulated NPP and chlorophyll have a very small bias in the global mean, the model tends to overestimate surface chlorophyll in the subtropical gyres and to underestimate it in the subpolar gyres.The spatial correlation between the model and the data of the long-term mean chlorophyll is relatively low (r < 0.4), while the correlations for the seasonal anomalies varies between r = 0.3 and r = 0.7.In the vertical, the correlation of NPP with annual mean vertical profiles is very high (r > 0.95), but substantially lower for chlorophyll, mostly because the model simulates the deep chlorophyll maximum at a too shallow depth.Interannual variability is well captured in the low latitudes, but not as well in the mid-and high latitude regions.The model skill in reproducing sea surface temperature and surface nutrient fields is consistently higher across most of the metrics than for the simulated ecological fields, and is intermediate for surface pCO 2 , and the air-sea fluxes of CO 2 and O 2 .
The model simulates a global mean annual particle export production of 5.9-6.7 Pg C yr −1 across the 100 m depth level.The flux across 75 m, which is the depth most studies reported their export fluxes for, is between 6.5 and 7.3 Pg C yr −1 .This is lower than previous estimates that reported particle export production of 10 Pg C yr −1 or above (Najjar et al., 2007;Gnanadesikan, 2004;Schlitzer, 2002;Dunne et al., 2007), but higher than another group of recent estimates that suggested a particle export production of between 4 and 5 Pg C yr −1 only (Henson et al., 2011;Lutz et al., 2007).It is beyond the scope of this article to resolve these discrepancies, but we note that the model estimate is nearly exactly in between these two groups.
While the focus in this work is on particle export production, we recognize the importance of DOC to total export production.In our simulation, between 0.9 and 1.1 Pg DOC yr −1 is exported across the 80 m level, i.e., about 12 % of the total organic matter export of about 8 Pg C yr −1 , lower, but not inconsistent with the model results summarized by Najjar et al. (2007), who suggested a DOC export fraction of about 20 %.
Recently, the representation of PFTs in the CCSM-BEC was compared to satellite-based estimates of PFT Biogeosciences, 10, 7373-7393, 2013 distribution (Alvain et al., 2008;Hirata et al., 2011), direct observations of PFT composition (Buitenhuis et al., 2006(Buitenhuis et al., , 2010) ) and to other ecosystem models by Hashioka et al. (2012), Sailley et al. (2013) and Vogt et al. (2013).The diatom fraction at the peak timing of blooms varies between about 20 % (Alvain et al., 2008) and 70 % (Hirata et al., 2011) in satellite-based estimates.The CCSM-BEC results are closer to the Hirata et al. (2011) estimate with high diatom fraction (> 80 %) at peak timing of bloom (Hashioka et al., 2012).In the annual mean, the location of regions in which diatoms dominate biomass are captured by the BEC.However, the extent of these regions is overestimated in the BEC compared to the satellite-based results (Vogt et al., 2013).In terms of seasonal succession, CCSM-BEC simulates a high diatom fraction at the peak timing of blooms (> 80 %), which is closer to the estimates of Hirata et al. (2011) (∼ 70 %) than those of Alvain et al. (2008) (∼ 20 %).Since the BEC simulates a generic zooplankton type that switches in its functionality, a comparison to data or other models is difficult.However, it turns out that the modelled generic zooplankton mostly resembles the microzooplankton class of other models (Sailley et al., 2013).

Calculation of trends
All trends presented in this work were computed using a linear regression on annual mean model output.Changes were obtained by multiplying the slope of the linear regression with the simulation length (47 yr).We tested for the significance of all trends with a two-sided Student t test (requiring a level of significance α = 0.05 or alternatively requiring a p value of the regression of less than 0.05).To account for autocorrelation in the time series, we reduced the degrees of freedom when performing the t test as described in e.g.Zwiers and von Storch (1995).Typically, this reduced the degrees of freedom by 36 %.We report primarily percent changes, which we obtained by normalising all results to the decadal mean of the first ten years (1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969).For the maps, the changes were calculated for each grid cell and non-significant changes are shown in white.
In this work, we describe the significant trends we found in our simulation and analyse the drivers of these changes.We do not attempt to attribute the trends to anthropogenic climate change.

Results
Our model simulation shows a significant decrease in global particle export production (EP) by −0.8 Gt POC yr −1 (−8 %) from 1960-2006, but with strong temporal and regional variability (Figs. 1 and 2d).This trend is highly significant, with a probability of < 1 % being of random nature.The strongest decreases occur in the western subtropical Pacific (−2 mol POC m −2 yr −1 resp.−40 %, Fig. 2d).In contrast, EP increases in the North Atlantic (+0.6 mol POC m −2 yr −1 resp.+30 %) and in parts of the Southern Ocean.However, these increases in the Southern Ocean are largely compensated by areas of decreases so that the Southern Ocean as a whole exhibits no significant trend in EP.All other areas have changes of less than 5 %.
The significant trends persist also for the period from 1979 onward, where satellite observations were used to build the forcing fields (see Sect. 2.2).The trends are observable in all seasons and are most pronounced in autumn (September-November, not shown).
Changes in simulated NPP are strongly correlated with the changes in EP (spatial Spearman correlation S c = 0.96).However, in relative terms, changes in NPP are between 2-20 % weaker than changes in EP.The amplification of trends in EP is due to changes in export efficiency (Fig. 1e, f).The changes in export efficiency show a pattern similar to the changes in EP (S c = 0.88) and also NPP (Figs. 1 and 2).As a result, the trends in EP are amplified strongest in areas where the largest changes in NPP occur, i.e., the western subtropical Pacific, the North and equatorial Atlantic and the Southern Ocean.
NPP (and EP) are negatively correlated with SST in the low latitudes, but show weak correlation in the high latitudes (Fig. 3).The global average NPP is negatively correlated with SST, reflecting mostly the changes in the low latitudes.
The simulated changes in NPP and EP are associated with extensive shifts in both the NPP of diatoms and of small phytoplankton (Fig. 4).While diatom NPP decreases globally only by −3 %, the regional changes are of much greater magnitude.In the North Atlantic and the Southern Ocean, diatom NPP increases by up to +30 %, except in the area west of the Antarctic Peninsula (120 to 80 • West, 50 to 70 • South), where diatom NPP decreases (−30 %).In the equatorial and western Pacific, and in particular in the subtropical gyres, diatom NPP decreases by up to 60 % over the 1960 to 2006 period.The changes in the distribution of diatom biomass correlate positively with changes in diatom NPP (S c = 0.87), with increases in the high latitudes and decreases in the low latitudes.
Small phytoplankton NPP shows strongly contrasting trends to diatom NPP in the high latitudes, but similar changes in the low latitudes (Fig. 4a, b).The global average decrease of small phytoplankton NPP is −8.5 %.The model simulates strong decreases in small phytoplankton NPP in the Southern Ocean (up to −45 %), except for the area west of the Antarctic Peninsula where small phytoplankton NPP increases.In almost all other areas small phytoplankton biomass decreases between −2 % and −40 %, with strongest decreases in the subtropical gyres, and in the North and equatorial Atlantic.Trends in the biomass of small phytoplankton follow the trends in the NPP of small phytoplankton (S c = 0.84), with the exception of the equatorial Pacific where the NPP of small phytoplankton decreases whereas its biomass increases.As a consequence of the changes in the The small panels to the right of each map show the zonal mean of the respective map using the same unit.Changes have been calculated using a linear regression on annual mean model output.Areas with insignificant changes (p value of the regression is less than 0.05) are shown in white color.More details on calculation of trends are given in Sect.2.5.Both total changes and changes in percent are given over the 47 yr of the simulation (not % yr −1 ).Export efficiency is the fraction of primary production that is exported through 100 m, calculated by dividing the particle export production at 100 m by integrated net primary production above 100 m.
NPP of diatoms and small phytoplankton, the diatom fraction substantially increases in the North Atlantic and the Southern Ocean by up to 40 %, except the area west of the Antarctic Peninsula.The diatom fraction decreases in the equatorial and subtropical Pacific, in particular in the western part (−20 %, −50 %).Zooplankton biomass follows changes in the NPP of small phytoplankton almost everywhere (Fig. 4), but zooplankton changes are smaller than the changes in the NPP of small phytoplankton.
The production of SiO 2 follows largely the trends in diatom fraction (S c = 0.77), as can be expected given that only diatoms produce SiO 2 in our model.The small deviations in the relative trends are caused by changes in the in Si : C ratio of diatoms.In total, the changes in the ratio account for only about 5 % of the total changes in SiO 2 production.Trends in CaCO 3 production largely follow trends in small phytoplankton biomass, which is a direct consequence of coccolithophores being modeled as a fraction of small phytoplankton.The trends are actually amplified in regions where nutrient limitation increases, e.g., the tropical Pacific, largely as a result of the model's calcification rates being dependent on nutrient limitation.In contrast, the temperature changes have little influence on the simulated CaCO 3 production.1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990 2000 2010 1960 1970 1980 1990   As CaCO 3 is only produced by small phytoplankton, whose biomass tends to be anticorrelated to that of diatoms, the rain ratio is strongly anticorrelated to the diatom fraction (S c = −0.81).Consequently, the rain ratio shows significant changes where the phytoplankton community composition changes strongly, i.e., in the Southern Ocean, North Atlantic and western Pacific.

What mechanisms drive PFT distribution, export
and export efficiency?

Physical changes
The model simulates substantial changes in the physical environment of the upper ocean over the analysis period, thereby driving many of the biological changes described in the preceding section.The surface ocean is simulated to have warmed by 0.5 to 1.  tropical Pacific (+1 psu), and decreases in the northern part of the tropical Pacific (−1 psu).These changes enhance stratification in the low latitudes, while stratification is simulated to decrease in the Southern Ocean and in the North Atlantic.
The forcing fields used for this simulation prescribe increases in wind stress, particularly in the Southern Ocean and North Atlantic, and weaker winds in the tropical Pacific.As a consequence of the combined changes in wind and stratification, the simulated mixed layer deepens in the North Atlantic and Southern Ocean by up to 40 %.In the equatorial Pacific, the mixed layer shoaled substantially (−40 %).The changes in stratification and in the depth of the mixed layer affect phytoplankton directly by changing the light availability and have strong effects on nutrient supply, with nutrient concentrations generally increasing where stronger mixing occurs and vice versa.

Growth limitation and the role of bottom-up control
Phytoplankton biomass, production, and PFT distribution is controlled by a complex interplay between bottom-up processes, i.e., changes in temperature, nutrient and light availability, and top-down processes, i.e., changes induced by grazing pressure from zooplankton.Trends in the biological properties can be caused by changes in nutrients or light availability (bottom-up controlled trend) or it can be caused by changes in grazing pressure (top-down controlled trend).We will first discuss the importance of bottom-up controls, and then discuss the role of top-down controls in explaining the model-simulated trends in biomass, production, and PFT composition.
The consequence of the altered physical conditions are substantial changes in both light and nutrient limitation, which are shown in Fig. 5 (zonal mean of changes) and Fig. 2b (time series for selected regions).
The general pattern for both phytoplankton types is one of increased light availability, but higher nutrient stress in the low latitudes, and decreased light availability, but lower nutrient stress in the high latitudes.The net effect is a substantial decrease of NPP in both regimes.
Temperature changes have only a weak effect on phytoplankton growth compared to changes in nutrients and light.Warming in the surface ocean slightly increases NPP, but this effect is compensated by decreases in NPP caused by decreased nutrient and light availability.The mechanisms we find are in accordance with several model studies projecting the impact of future climate change on marine productivity (Steinacher et al., 2009;Bopp et al., 2001;Boyd and Doney, 2002;Sarmiento et al., 2004).
While the general pattern in nutrient and light changes is similar for diatoms and small phytoplankton, the magnitude of growth limitation is different for the two PFTs, leading to differential responses to the changes in the physical environments, and hence changes in the relative PFT compo- The limitation factor is a unitless value between zero and one, with zero representing maximal limitation and one representing unlimited growth.An increase in limitation factor leads therefore to an increase in growth.For this plot, changes in the limitation factor have been calculated for each cell as described in Sect.2.5, and the zonal means of the resulting changes are shown.
sition.The growth of both diatoms and small phytoplankton is calculated with a multiplicative growth function where the growth rate of phytoplankton group x is calculated by , where µ max denotes the maximum growth rate for diatoms and small phytoplankton (given in [d −1 ]) and T f , N x and L x are dimensionless temperature, nutrient and light limitation terms (see Appendix and Table 1).Small phytoplankton are parameterized to have lower nutrient requirements than diatoms.They are also not limited by silicic acid, but they have higher light requirements compared to diatoms.In contrast, the diatoms are better adapted to low light, but thrive better in high nutrient regimes.Temperature affects both phytoplankton PFTs (p-PFTs) equally and therefore cannot be the cause for any difference in growth rates of the two PFTs.An illustration for the (resulting) difference in growth rate of diatoms and small phytoplankton at different iron, nitrate and light levels is shown in Fig. 6.These differences in growth rates, induced by the different limitation factors (Fig. 5), can explain many of the simulated trends as discussed in the following.
In the tropical Pacific the warming of the surface ocean, stronger stratification and lower wind stress led to higher nutrient stress for both p-PFTs and reduced their growth rates.As this is a regime where small phytoplankon have a significant growth advantage compared to diatoms (Fig. 6), the growth rate of small phytoplankton decreased less than that of the diatoms.Thus, the biomass and productivity by small phytoplankton increased relative to that of diatoms, explaining the shift toward a lower diatom fraction.
Under the conditions typical for the Southern Ocean and North Atlantic (high nutrient concentration but limiting light availability), small phytoplankton have only a small growth advantage compared to diatoms.During the last 50 yr, both regions experienced an increase in wind stress, which enhanced light limitation and reduced nutrient limitation for both pPFTs.
In the Southern Ocean, the growth rate of both pPFTs decreased over the course of the analysis period, but less so for diatoms than for small phytoplankton, leading to growth rates of more similar magnitude for both pPFTs at the end of the simulation.This implies an overall reduction in biomass for both groups, but with a shift toward a higher diatom fraction.
Also the changes in the North Atlantic can be explained by bottom-up processes.Here the observed shift toward a higher diatom fraction is consistent with the growth limitation of diatoms having become smaller over the course of the 1960-2006 period, while that of small phytoplankton has increased.
These bottom-up limitation factors are quite successful in explaining the reported changes in pPFT biomass and NPP.However, it is puzzling that diatoms can benefit at the expense of small phytoplankton, since small phytoplankton always have a higher growth rate than diatoms under the nutrient and light regimes that phytoplankton experience in the ocean (Fig. 6).Thus, if only bottom-up factors were considered, one would expect small phytoplankton to always dominate biomass.In addition, not all trends can be explained exclusively by the trends in the growth limitation factors.For example, diatom biomass increased in the Southern Ocean although its growth limitation stayed rather constant.This suggests that also top-down controls, i.e., those involving the impact of grazing pressure by zooplankton, are of importance in understanding changes in phytoplankton biomass.

Zooplankton grazing and the role of top-down control
Grazing on phytoplankton type x is calculated in the model according to a Holling type III functional response (Holling, 1965).Grazing has the same temperature dependence as phytoplankton (see Appendix).The maximal zooplankton growth rate is higher for small phytoplankton than for diatoms.Therefore, zooplankton biomass is strongly coupled to small phytoplankton biomass and closely follows trends in small phytoplankton NPP (S c = 0.79, Figs. 4 and 8).
In the real ocean changes in zooplankton biomass are influenced by a variety of factors such as the presence of larger predators or direct effects of climate change on the zooplankton life cycle and resulting phenological mismatches (Edwards and Richardson, 2004), while zooplankton biomass in the model depends only on temperature, phytoplankton biomass and phytoplankton composition.Of these variables, temperature is the only one that can be directly affected by climate change, while the others influence zooplankton and their grazing pressure via bottom-up controls.It turns out that the model-simulated temperature trends of about 0.3 • C on global average can only drive about 5 % of the observed trends in PFT biomass, implying that biomass trends in our simulation are mostly bottom-up controlled.However, zooplankton are very important in setting the overall phytoplankton biomass level and its PFT composition and production.For example, the presence of grazing pressure strongly contributes to the dominance of diatoms (see below and Hashioka et al., 2012;Sailley et al., 2013).But since phytoplankton are primarily bottom-up controlled in our simulation, zooplankton tend to respond to changes in phytoplankton and not vice versa.Trends in phytoplankton PFT biomass as well as PFT composition are initiated by changes in bottom-up factors and are then being modified by changes in zooplankton grazing.
To illustrate the response of grazing pressure to changing phytoplankton biomass, we analyse the specific grazing rate SG(P x ), which is the grazing rate G(P x ) on phytoplankton P x , normalized by zooplankton concentration Z (Hashioka et al., 2012): where u max P x is the maximum grazing rate, T f the temperature dependence and g x a zooplankton grazing coefficient; see Table 1 and Appendix.Figure 6c shows the difference in specific grazing rate in dependence of the prey biomass ratio.At equal biomass, small phytoplankton experience a higher grazing pressure than diatoms, with this difference in grazing pressure becoming larger as total biomass increases.
The higher grazing pressure of zooplankton on small phytoplankton partly compensates for the advantage in growth rate of small phytoplankton and makes it possible for diatoms to dominate biomass in regions where their growth rate is close to the growth rate of small phytoplankton, i.e., in high-nutrient, low-light regimes such as the Southern Ocean and North Atlantic.This top-down control of phytoplankton composition in the BEC has also been analysed in Hashioka et al. (2012).
In the tropical Pacific, the light-and nutrient-limitationinduced reductions in the biomass of small phytoplankton and diatoms caused a decrease in zooplankton biomass.This results in less grazing of both phytoplankton-PFTs, but little changes in the relative grazing pressure (Fig. 6, yellow circles).Trends in the tropical Pacific are therefore only weakly modified by top-down processes.
In the North Atlantic, the bottom-up-induced increase in diatom biomass (while small phytoplankton biomass remained relative constant) results in a stronger relative grazing on diatoms.Trends in the North Atlantic are therefore weakened by top-down processes.
In the Southern Ocean, zooplankton grazing turns out to be key for understanding the relative shifts in phytoplankton biomass.In this region, small phytoplankton biomass decreased over the course of the simulation because of changes in light and nutrient limitation.For diatoms a lower light limitation was compensated by a stronger nutrient limitation.Zooplankton does not follow the increase in total phytoplankton in the Southern Ocean, as the increase in total biomass is associated with a strong shift towards more diatoms.Zooplankton growth rates are lower when feeding on diatoms, thus, zooplankton biomass decreased in line with the decrease in small phytoplankton.The specific grazing pressure shifted towards a stronger diatom grazing.But as the zooplankton biomass decreased, the total grazing pressure decreased as well, permitting diatoms to thrive.
In summary, trends in PFT composition can mostly be explained by changes in light and nutrient availability.Temperature-induced changes in zooplankton grazing explain only a small fraction (< 5 %) of the simulated phytoplankton biomass changes.
Changes in top-down control, which have been caused by changes in phytoplankton biomass, weaken phytoplankton trends in the high latitudes, but top-down control has little influence on biomass trends in the low latitudes.Therefore we identify bottom-up processes as the primary driver of the changes in biomass.
In order to understand how the changes in PFT composition are connected to the amplification of EP compared to NPP, we analyse the pathways along which carbon is routed from PFT biomass to sinking particles.

Carbon pathways from PFT biomass to sinking particles
In the BEC, particles are formed through four different pathways.Two of them are grazing-based, i.e., particles are formed by zooplankton grazing on either diatoms or small phytoplankton, and two are aggregation-based, i.e., particles form by aggregation of either small phytoplankton or diatoms.
Figure 7 shows that in the long-term mean, the production of particles in the low and mid-latitudes is dominated by grazing on small phytoplankton (55 % of total POC), with only small exceptions.The reason is that small phytoplankton dominate biomass in the low and mid-latitudes, and grazing pressure is amplified by high temperatures.Aggregation in the low and mid-latitudes is less relevant, as biomass is kept low due to grazing and aggregation is a quadratic function of biomass (see Appendix).
In the high latitudes, diatom aggregation tends to dominate with 54 % of total POC in Southern Ocean originating from diatom aggregation.In these regions, the diatom fraction is high and grazing pressure is low because of the cold temperatures.Biomass reaches high levels, which results in significant aggregation fluxes.Small phytoplankton aggregation dominates particle formation in areas where small phytoplankton has a high biomass and a high fraction of calcifiers.Grazing on diatoms dominates POC production only in few localized areas (Fig. 7).
This time-mean dominance pattern varies little over time.Over the 50 yr of simulation, regions where POC production is dominated by small phytoplankton aggregation decreased slightly (−4 %) in favor of regions where diatom aggregation dominates POC production (not shown).
However, the composition of POC reflects the changes in PFT structure (Figs. 8 and 9).Between the beginning

Laufkötter et al.: PFTs and export production 1960-2006
Diatom Aggregation Grazing of Diatoms Small Phyto Aggregation Grazing of Small Phyto Fig. 7. Map showing which particle production mechanism is quantitatively dominant over the study period.The particle production mechanisms in the BEC are as follows: aggregation by diatoms (blue), aggregation by small phytoplankton (yellow), zooplankton grazing of diatoms (green) and zooplankton grazing of small phytoplankton (red).We calculated the study period mean  of all fluxes to determine which mechanism contributed strongest to the sinking particle pool.(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970) and the end (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) of our analysis period, the relative fraction of POC that can be attributed to small phytoplankton grazing decreased from 45 to 38 % in the North Atlantic, in favor of diatom aggregation.Note that total small phytoplankton aggregation and grazing shows little changes, but total POC production increased, driven by higher diatom aggregation.In the Southern Ocean, the fraction of POC produced by diatom aggregation increased from 29 to 33 % at the expense of the fraction produced by small phytoplankton grazing (−3 %) and aggregation (−3 %).This reflects increases in diatom biomass and decreases in small phytoplankton and zooplankton.In the tropical Pacific, the relative fraction produced by small phytoplankton grazing increased from 73 to 79 % at the expense of diatom grazing (−4 %) and small changes in aggregation.These changes are associated with decreases in all PFTs and decreases in diatom fraction.

Relation between particle production mechanism and export efficiency
An increase in the diatom fraction leads to a higher fraction of particles ballasted with SiO 2 and a stronger particle production during grazing.Consequently, the diatom fraction is strongly positively correlated with the export efficiency almost everywhere in our simulation (S c > 0.8, Fig. 10).However, there are two notable exceptions: at about 50 • South and around 50 • North in the Pacific, the correlation is significantly negative (S c between −0.8 and −1).
In those areas, two factors influence the export efficiency of small phytoplankton: first, small phytoplankton include a high coccolithophore fraction, which leads to a strong production of CaCO 3 .CaCO 3 is parameterized to be a more effective ballasting material than SiO 2 (following Klaas and Archer 2002), therefore causing a higher export efficiency.Second, aggregation is parameterized to depend quadratically on biomass in our model.Thus a decrease in small phytoplankton at a high biomass level causes a stronger reduction in particle production than an increase in diatoms at a low biomass level can compensate for.
A third factor influencing export efficiency in our model is the biomass level of small phytoplankton.At low biomass levels, small phytoplankton is assumed to be composed of nano/picophytoplankton species.The fraction of grazed matter that is routed to POC is parameterized to be smaller to reflect a stronger microbial loop.Hence, a decrease in small phytoplankton biomass leads to a lower particle export efficiency.In our simulation this effect dominates in the low latitudes, where diatom biomass is very low and grazing on small phytoplankton constitutes the most important particle export flux.
In summary, export efficiency (Fig. 1e and f) decreases strongly in the low latitudes, because of both a lower diatom fraction and decreases in small phytoplankton biomass.Export efficiency increases in the North Atlantic driven by a higher diatom fraction.In the area south of 30 • S, the diatom fraction increases, but small phytoplankton biomass and coccolithophore fraction decreases strongly, causing an overall decrease in export efficiency.On the global scale, diatom fraction increases but export efficiency decreases (Fig. 2e), largely because of the low latitude trends in diatom fraction and export efficiency.

Comparison with observed estimates of long-term changes in Chl, NPP and EP
Observations of long-term changes in Chl, NPP, and EP are scarce, particularly when requiring that they extend beyond 20 yr.In addition, the vast majority of reported trends concern changes in one region only, with only two studies so far (Boyce et al., 2010;Wernand et al., 2013) having made an attempt to determine long-term chlorophyll trends over the global ocean.However, the Boyce et al. ( 2010) study was very controversially received (see introduction), and Wernand et al. ( 2013) do not report results for the Southern Ocean (south of 30 • S).Thus, to our knowledge, there are no observations that would permit us to assess our results on a long-term and global scale.Nevertheless, there are several available data sets that are useful to evaluate the model at least on the local to regional level, and to assess the model performance over shorter time periods.
and BATS are located in regions where rather weak or insignificant trends occur over the full simulation period.For the 1989-2007 period our simulation shows an increase in chlorophyll and NPP at HOT of 15 %, but a decrease of about 5 % at BATS.Moving to larger scales, we only consider the studies making use of post-1997 satellite data (e.g., GlobColour product), since we consider the comparison of the satellite-based estimates of chlorophyll for 1978-1986 (Coastal Zone Color Scanner) with those from 1997 onward as too problematic despite several efforts (Antoine, 2005;Gregg, 2003).In con-Fig.10.Spearman correlation (S c ) of diatom fraction with export efficiency.Trends were removed from the time series before calculation.Significance has been tested at significance level α = 0.1.trast to Behrenfeld et al. (2006), our results do not show significant trends over the 1997 through 2006 period, but we do capture the correlation between NPP and SST (Spearman correlation = −0.81)over this period quite well.In fact, this correlation exists throughout the 1960-2006 period (Spearman correlation −0.87, see also Fig. 3).
Regarding long-term trends on the global scale, Wernand et al. (2013) report decreases in chlorophyll for the North and equatorial Pacific and Indian Ocean but increases for the North and equatorial Atlantic Ocean that sum up to no significant global trend in chlorophyll in the last century.They do not report results for the southern Atlantic and Pacific (south of 30 • S) because of too few observations.We capture the direction and magnitude of trends reported in Wernand et al. (2013) in all major ocean basins except the equatorial Atlantic, where we find a decrease in chl, and the North Atlantic, where our chlorophyll increase is significantly weaker.In addition, our results show a significant global decrease in chlorophyll of 3 %.Considering that a trend of such small magnitude may not be detectable in sparse data sets, our results are not inconsistent with the findings of Wernand et al. (2013).However, this contrasts with the strong global trend reported by Boyce et al. (2010).We compare our data only to that part of the Boyce et al. (2010) record that is based on direct chlorophyll measurements, omitting observations that are based on transparency observations.This avoids some of the controversial aspects of the Boyce et al. ( 2010) study (Rykaczewski and Dunne, 2011;Mackas, 2011).Even with this selection, their global decline in chlorophyll is nearly an order of magnitude larger than what we infer from our model.On the regional scale, our trends follow the direction (not the magnitude) of in situ chlorophyll trends described in Boyce et al. (2010) in the equatorial Atlantic and Pacific, northern Indian Ocean and parts of the Southern Ocean, i.e., in about half of the regions.Our simulated trends differ in several other regions.In particular, in the North Atlantic our simulation shows a strong increase in chlorophyll.Furthermore, Boyce et al. ( 2010) report increases in chlorophyll in the North Pacific, South Indian and South Pacific, where our simulation shows decreases (North and South Pacific) or no clear trend (South Indian).However, one needs to take into consideration that the number of measurements on which Boyce et al. ( 2010) based their trends are sparse in general and extremely low in the Southern Hemisphere, where our trends differ most from their results.At the moment it is not possible to distinguish whether the model is seriously underestimating the trends, particularly in the Southern Hemisphere, or if the observational record is not reliable enough to derive robust trends.

Comparison with model studies
While there are many model studies examining primary and export production under future climate change (Steinacher et al., 2009;Bopp et al., 2001;Boyd and Doney, 2002;Sarmiento et 2004;Bopp et al., 2005Bopp et al., , 2013)), we provide the first detailed analysis of changes in PFTs, NPP and EP in a model hindcast over the 1960 through 2006 period driven with prescribed atmospheric forcing.Other published model studies covering the last decades have used fully coupled models, which produce their own particular evolution of weather and climate over the second half of the 20th century (Steinacher et al., 2009;Bopp et al., 2013).
Our projected changes for the 1960-2006 period are consistent in sign, but larger in magnitude than the trends reported by Steinacher et al. (2009) and Bopp et al. (2013) over the same period.These coupled models simulated a linear decrease of less than 3 % for both NPP and EP for the same period.Most current fully coupled ocean models have limited skill in reproducing observational-based estimates of global NPP (M.Kidston, personal communication, 2012, see also Anav et al., 2013).Thus it is conceivable that these models underestimate the changes in the last decades.

Implications for future change
We discuss two questions to assess the possible implications for future changes: first, are our trends in NPP and EP likely to continue in the next century?And second, which effects will changes in PFT distribution have for future export production?
The trend towards increased stratification in the low latitudes and related decreases in NPP seen in our simulation for the past five decades is projected to continue in nearly all studies that simulate anthropogenic climate change in the next 100 yr (Steinacher et al., 2009;Bopp et al., 2001;Boyd and Doney, 2002;Le Quéré et al., 2003;Bopp et al., 2013).Human-induced warming of the surface ocean, with important regional modifications by salinity is the driver for these increases in stratification (Capotondi et al., 2012).
In the North Atlantic, future climate change predictions show robust increases in stratification and declines in NPP driven by increases in greenhouse gases (Capotondi et al., 2012;Steinacher et al., 2009).This contrasts with the reduced stratification and higher NPP we simulated in these regions over the past 5 decades, driven by stronger wind stress caused by by an overall positive trend in the North Atlantic Oscillation (Hurrell et al., 2001).In most climate simulations, the NAO does not continue to become even more positive, therefore permitting the stratification impact to emerge from the potential masking effect of the NAO trend over the last five decades.
For the next century, increases in NPP and EP are predicted for the polar Southern Ocean, but there is little agreement among models for the region between 40-60 • S (Bopp et al., 2013).Over the past five decades, westerly winds increased in the Southern Ocean, most likely driven by ozone depletion above Antarctica and increases in greenhouse gases (e.g., Marshall 2003).However, the NCEP reanalysis (which is the basis of the CORE-CIAF v2 forcing used in our simulation) is known to overestimate winds in the Southern Ocean (Marshall, 2003).Therefore wind-driven changes in our simulation might be overestimated in the Southern Ocean.In addition, the strong variability in the Southern Ocean makes it difficult to detect a climate change signal (Boyd et al., 2008).We conclude that caution is required when interpreting our trends in the Southern Ocean and possible consequences for future climate.
Regarding the second question, the effects of PFT distribution on export efficiency, our results show a negative correlation between stratification and diatom fraction.This correlation has also been found in studies that analyse future climate change in different ocean models (Marinov et al., 2010;Bopp et al., 2005).Furthermore, the effects of changes in diatom fraction on export efficiency are complex in our simulation.We find a positive correlation between interannual variability Biogeosciences, 10, 7373-7393, 2013 of diatom fraction and export efficiency, except in those regions with high coccolithophore fraction (Fig. 10).While this is in accordance with previous assumptions about the relationship between diatom fraction and export efficiency (Le Quéré et al., 2005), diatom fraction is not the main driver of long-term global trends in export efficiency in our simulation.Other important factors influencing export efficiency are the coccolithophore fraction and the concentration of small phytoplankton biomass, which is associated with the efficiency of the microbial loop in our model.In contrast to our results, Bopp et al. (2005) showed that the change in diatom fraction resulted in a reduced efficiency of the biological pump on global scale under a future 4 × CO 2 scenario.Changes in export production (−25 %) have been significantly stronger than changes in NPP (−15 %) in their study.
Considering the contrasting mechanisms in our hindcast simulation compared to the mechanisms identified in the future study by Bopp et al. (2005), projections of export efficiency seem to be strongly dependent on the ecosystem representation in the model.Further model studies are needed to explore the range of possible interactions between PFT composition and export efficiency.

Caveats and limitations
In this work we determine and discuss significant trends in ocean NPP and EP in the last 47 yr, but without making an attempt to attribute the trends to anthropogenic climate change.Any trend detection is very sensitive to the level of variability, which is substantial, and potentially underestimated in our model relative to observations.Thus, our model may detect significant trends too early (see also Henson et al., 2011); particularly in regions with weak trends but high variability (e.g.Southern Ocean), our trend estimates may be biased.
Our model-based results are potentially also quite sensitive to the details of the ecosystem model, such as the employed parameterizations of biological processes and the associated parameter choices.Important mechanisms that are not considered in our study include the effects of ocean acidification on PFT growth and on calcification rates, and a potential decoupling of the C:N:P ratio.These processes affect the PFT distribution and global biogeochemical cycles and have previously been found to be of relevance for the export efficiency (Weber and Deutsch, 2012;Doney et al., 2009a).The consideration of these processes in ecosystem models has been hampered by sparse availability of data on PFT distribution, traits and physiology.The recent development of new data sets describing PFT abundance (Buitenhuis et al., 2013) and nutrient utilization traits of PFTs (Edwards et al., 2012) is an important step forward and will help model development and improve performance of ecosystem models.Another limitation of our study related to ecosystem representation is the temperature dependency of plankton growth.Taucher and Oschlies (2011) argue that a better representation of the temperature sensitivity of primary pro-duction could even change the direction of projected NPP trends.Schmittner et al. (2008) and Taucher and Oschlies (2011) both find a temperature-driven increase in global NPP but a decrease in export under SRES A2.Furthermore, heterotrophic processes such as grazing and remineralization have been assumed to show stronger responses to elevated temperatures (Pomeroy and Wiebe, 2001;Riebesell et al., 2009) than NPP, but our model assumes the same temperature sensitivity for phytoplankton growth and for zooplankton grazing.We need further research to better understand the balance between NPP and grazing processes under increasing temperatures (Riebesell et al., 2009) and to temperature sensitivity of phytoplankton growth in general.In addition, the representation of the second trophic level in the BEC as one single zooplankton class mimics the behaviour of microzooplankton (Sailley et al., 2013).Effects of mesozooplankton grazing and other higher trophic levels could modify phytoplankton biomass and might be insufficiently considered in our simulation.Likewise, impacts of climate change on larger predators or on the zooplankton life cycle might exert additional changes in top-down control on phytoplankton, which are not represented here.
Further uncertainties are associated with our forcing.In particular, the NCEP reanalysis (which is the basis of the CORE-CIAF v2 forcing used in our simulation) is known to overestimate winds in the Southern Ocean (Marshall, 2003).Therefore wind-driven changes in our simulation might be overestimated in the Southern Ocean, requiring some caution when interpreting our trends in the Southern Ocean and possible consequences for future climate.
Finally, we focus our analysis exclusively on the particle export production, while DOC may contribute up to 20 % to total export production (Najjar et al., 2007).Changes in DOC export may thus cause substantial changes in total carbon export and needs to be included in further studies to complete our picture of future changes in marine export production.

Conclusions
We present the first analysis of changes in plankton community structure, NPP and EP in a model hindcast from 1960 onwards driven with prescribed atmospheric forcing, i.e., winds, and fluxes of heat and freshwater.Our results suggest a significant global decline in NPP (−6.5 %) and EP (−8 %).Simulated changes in NPP and EP go along with global decreases in ocean chlorophyll of 3 % over the last five decades.Our downward trend supports the conclusions drawn by Boyce et al. ( 2010), but our simulation suggest a tenfold smaller trend.We could not resolve this large discrepancy, but tend to question the robustness of the observational trends.Simulated changes over the past 5 decades are more pronounced than those estimated from simulations with fully coupled atmosphere-ocean models forced with reconstructed CO 2 emissions.A possible reason could be that fully coupled We identified increased stratification in the low latitudes and increased wind stress and mixing in the North Atlantic and Southern Ocean as the main drivers for changes in PFTs and NPP.The trends in our simulation are mostly driven by bottom-up controls, with decreased nutrient concentrations being the main driver in the low latitudes and decreased light availability the main driver in the higher latitudes.However, the representation of top-down control is restricted to temperature effects on grazing of one generic zooplankton type in our model.Adding a representation of higher trophic levels and their sensitivity to physical changes might emphasize the role of top-down control.
Finally, our analysis reveals that the dynamics of particle production depend on a complex interplay of diatom, coccolithophore and zooplankton fraction, temperature and total biomass.This is in contrast to previous studies which show a linear relationship between diatom fraction and export efficiency (Bopp et al., 2005).The link between PFT distribution and export efficiency is currently not well understood.Further model studies should focus on determining the full range of possible responses of export efficiency to changes in PFT composition.Furthermore, an improvement of the representation of PFT distributions and their contribution to the particle flux is necessary.For this we need both better estimates of PFT distributions and measurements of the contribution of different PFTs to the particle flux.New data sets of PFT biomass and pigment concentrations (Buitenhuis et al., 2013) and more accurate satellite-based PFT estimates (Alvain et al., 2008;Hirata et al., 2011) provide a promising step towards improved model evaluation.Climate change is expected to accelerate in the upcoming century, and stronger perturbations of the marine ecosystems are likely to affect PFT distribution, primary and export production to a wider extent than shown in this work.Including improved representations of PFT distribution and behaviour in simulations of future climate change might reveal an important feedback of the biological pump to climate change.

Model equations and parameters
We give the most important model equations in the following.Model parameters are written in lowercase and the values are given in Table 1.

A1 Changes in PFT biomass
The changes in biomass of a p-PFT P x are calculated as where µ x denotes the growth rate, G x the loss due to grazing, agg x the loss due to aggregation and loss x representing the non-grazing mortality.
The growth rate for phytoplankton type x is calculated as where µ max denotes the maximum growth rate, T f the temperature dependence and N x , L x represent nutrient and light limitation.The temperature sensitivity T f is calculated as where q 10 is a temperature reference factor, T temperature and t ref a reference temperature (s.Table 1).Nutrient limitation N x is the minimum of the specific nutrient limitations: For n = Fe, P and SiO 2 , specific nutrient limitation is calculated as where K n x are the Michaelis-Menten half saturation coefficients.For NO 3 and NH 4 specific nutrient limitation is calculated as The light limitation L x is calculated according to Geider et al. (1998) as where I par is irradiance, α is the initial slope of the photosynthesis-irradiance (P -I ) curve and ( Chl C ) x is the chlorophyll to carbon ratio of phytoplankton x.Note that light limitation depends on nutrient limitation, with a higher nutrient limitation value leading to a lower light limitation value.

A2 Grazing and changes in zooplankton biomass
The fraction of phytoplankton biomass P x that is grazed by zooplankton Z is calculated with Holling type III function (Holling, 1965): ).
Biogeosciences, 10, 7373-7393, 2013   is the maximum zooplankton growth rate on phytoplankton x, g is a zooplankton grazing coefficient and f x z a scaling factor for grazing (see Table 1).The grazed fraction is partly incorporated into new zooplankton biomass, a part is remineralized and a fraction is routed to the sinking detritus pool.

A3 Formation of POC
The particulate organic carbon originates from grazing and aggregation and is calculated as where G POC diat and G POC small are the fractions of grazed diatoms and small phytoplankton that are routed to the sinking particles and agg diat and agg small aggregation by diatoms and small phytoplankton, respectively.

A4 Routing of grazed matter to POC
The fraction of grazed matter of phytoplankton x that is routed to POC, G POC x , is calculated as small is a small phytoplankton grazing factor (see Table 1).

A5 Aggregation
Aggregation of phytoplankton x is calculated as agg x = min(a max x , p x × (P x ) 2 ), where a max x and p x are the maximum aggregation rate and a quadratic mortality rate of phytoplankton x (see Table 1).
In addition, at least 1 % of diatom biomass aggregate and sink.

Fig. 1 .
Fig. 1.Simulated changes in net primary production (a and b), particle export production (c and d) and export efficiency (e and f) from 1960 to 2006.The small panels to the right of each map show the zonal mean of the respective map using the same unit.Changes have been calculated using a linear regression on annual mean model output.Areas with insignificant changes (p value of the regression is less than 0.05) are shown in white color.More details on calculation of trends are given in Sect.2.5.Both total changes and changes in percent are given over the 47 yr of the simulation (not % yr −1 ).Export efficiency is the fraction of primary production that is exported through 100 m, calculated by dividing the particle export production at 100 m by integrated net primary production above 100 m.

Fig. 2 .
Fig. 2. Time series in percent for selected regions from 1960 to 2006.(a) Indicates the location of the different regions.Column (b)presents time series for the total limitation, with a higher limitation value leading to more production.Column (c) shows time series for PFT biomass, column (d) particle export production (EP) and integrated net primary production (NPP), and column (e) shows export efficiency and diatom fraction.Trends have been calculated with linear regression, using annual mean model output.Significance of trends has been tested with a t test (p value of the regression is less than 0.05), for significant trends regression lines have been added.More details on calculation of trends are given in Sect.2.5.

Fig. 3 .
Fig. 3. Relationship between anomalies of annual mean depthintegrated NPP [mol C m −2 yr −1 ] and anomalies of annual mean SST [ • C] for three regions.Each dot represents one year, with the low latitude average (30 • S-30 • N) shown in dark blue, high latitude average shown in light blue and the global average shown in red.Anomalies are defined here as the deviation of annual mean SST (NPP) from the 1960-2006 mean of SST (NPP) in the respective regions.
5 • C in the low and mid-latitudes, but there are also large regions that experienced substantial cooling, i.e., in the North Pacific and parts of the Southern Ocean with a net cooling of up to −1 • C. The global surface ocean warmed by +0.3 • C, which is slightly below the estimate by Smith et al. 2008 (+0.4 • C).Salinity increases in the North Atlantic (+0.8 psu) and in the southern part of the

Fig. 4 .
Fig. 4. Simulated percentual changes between 1960 and 2006 in (a) small phytoplankton NPP and (b) diatom NPP, (c) zooplankton biomass and (d) resulting changes in diatom fraction.The small panels to the right of each map show the zonal mean of the respective map.Changes have been calculated using a linear regression on annual mean model output.Areas with insignificant changes (p value of the regression is less than 0.05) are shown in white.More details on calculation of trends are given in Sect.2.5.

Fig. 5 .
Fig. 5. Zonal mean changes in (a) light and (b) nutrient limitation of diatoms and (c) light and (d) nutrient limitation of small phytoplankton.The limitation factor is a unitless value between zero and one, with zero representing maximal limitation and one representing unlimited growth.An increase in limitation factor leads therefore to an increase in growth.For this plot, changes in the limitation factor have been calculated for each cell as described in Sect.2.5, and the zonal means of the resulting changes are shown.

Fig. 6 .
Fig. 6.(a) and (b) Difference between small phytoplankton growth rate and diatom growth rate at different nitrate, iron and light values.Typical values for regions discussed in the text are marked with coloured circles.The Southern Ocean is a mainly iron-limited region, while the subtropical gyres and the North Atlantic are mainly nitrate-limited.The growth rate [ 1 d ] is multiplied with the biomass [mmol C m −3 ] to obtain the total growth.(c) Small phytoplankton-specific grazing/diatom-specific grazing at different phytoplankton biomass levels and constant temperature.
report an average increase in NPP of 2 % per year for the 1989-2006 period for the tropical North Atlantic (BATS) and tropical East Pacific (HOT).According to our simulation, both HOT

Fig. 9 .
Fig. 9. Changes in POC composition in % in the Southern Ocean,tropical Pacific and North Atlantic between beginning (1960-1970)  and end (1996-2006)  of simulation period.We obtain the changes by calculating the percentage of each mechanism's contribution to total POC and then computing the difference between the percentaged fraction of POC at the beginning and the end of the study period.The changes are significant within the 95 % confidence interval.
own internal variability which does not correspond to the observed variability.
that are routed to POC, and e POC

Table A1 .
Doney et al. (2009b)osystem model BEC.Bold numbers in brackets denote values that have been used in the version ofDoney et al. (2009b).