Role of zooplankton dynamics for Southern Ocean phytoplankton biomass and global biogeochemical cycles

. Global ocean biogeochemistry models currently employed in climate change projections use highly simpliﬁed representations of pelagic food webs. These food webs do not necessarily include critical pathways by which ecosystems interact with ocean biogeochemistry and climate. Here we present a global biogeochemical model which incorporates ecosystem dynamics based on the representation of ten plankton functional types (PFTs): six types of phytoplankton, three types of zooplankton, and heterotrophic pro-caryotes. We improved the representation of zooplankton dynamics in our model through (a) the explicit inclusion of large, slow-growing macrozooplankton (e.g. krill), and (b) the introduction of trophic cascades among the three zooplankton types. We use the model to quantitatively assess the relative roles of iron vs. grazing in determining phytoplankton biomass in the Southern Ocean high-nutrient low-chlorophyll (HNLC) region during summer. When model simulations do not include macrozooplankton grazing explicitly, they systematically overestimate Southern Ocean chlorophyll biomass during the summer, even when there is no iron deposition from dust. When model simulations include a slow-growing macrozooplankton and trophic cascades among three zooplankton types, the high-chlorophyll summer bias in the Southern Ocean HNLC region largely disappears. Our model results suggest that the observed low phytoplankton biomass in the Southern Ocean during summer is primarily explained by the dynamics of the Southern Ocean zooplankton community, despite iron limitation of phytoplankton community growth rates. This result has implications for the representation of global biogeochemical cycles in models as zooplankton faecal pellets sink rapidly and partly control the carbon export to the intermediate and deep ocean.

the relative roles of iron vs. grazing in determining phytoplankton biomass in the Southern Ocean high-nutrient lowchlorophyll (HNLC) region during summer. When model simulations do not include macrozooplankton grazing explicitly, they systematically overestimate Southern Ocean chlorophyll biomass during the summer, even when there is no iron deposition from dust. When model simulations include a slow-growing macrozooplankton and trophic cascades among three zooplankton types, the high-chlorophyll summer bias in the Southern Ocean HNLC region largely disappears. Our model results suggest that the observed low phytoplankton biomass in the Southern Ocean during summer is primarily explained by the dynamics of the Southern Ocean zooplankton community, despite iron limitation of phytoplankton community growth rates. This result has implications for the representation of global biogeochemical cycles in models as zooplankton faecal pellets sink rapidly and partly control the carbon export to the intermediate and deep ocean.

Introduction
Phytoplankton, zooplankton and heterotrophic bacteria (including both Bacteria and Archaea, herein called "bacteria") in the oceans control important ecosystem processes and services (Ducklow, 2008), including primary, secondary and export production. Primary production, i.e. the production of organic matter by photoautotrophs using inorganic nutrients, can be either particulate and serve as food for heterotrophs, from protists to fish larvae, or dissolved and used by bacteria. Secondary production, the fraction produced by zooplankton grazing on phytoplankton, other zooplankton, or organic detritus, serves as food for larger organisms in the ocean, including fish and mammals. Export production, the fraction of primary production that sinks below the surface mixed layer, exerts an influence on marine biogeochemistry and climate as sinking organic matter remineralised to inorganic matter at depths becomes isolated from the atmosphere for decades to centuries. Export production responds primarily to the activity of large plankton, particularly the production and sinking of faecal pellets of zooplankton (e.g. copepods and euphausids) as well as the aggregation of diatoms, for example, during intense blooms. Export production reduces the surface concentration of inorganic carbon and maintains atmospheric CO 2 about 200 ppm lower than it would be in the absence of biological activity (Maier-Reimer et al., 1996). In contrast, bacteria and small zooplankton (e.g. heterotrophic flagellates and ciliates) remineralise and recycle organic matter in the upper ocean, thus reducing the quantity of organic matter that is exported. These ecosystem processes are controlled by the state of the environment (e.g. temperature, light, available nutrients, vertical mixing), and are modulated by the ecosystem structure of the planktonic community.
Dynamic green ocean models have been developed and used in global biogeochemical studies to understand and quantify the interactions between marine ecosystems and the environment. In these models, phytoplankton and zooplankton are grouped by taxa into plankton functional types (PFTs) according to their specific and unique roles in marine biogeochemical cycles (Hood et al., 2006;Le Quéré et al., 2005). Although generally only a small number of PFTs are treated explicitly, their inclusion has been shown to improve the realism of model simulations. For example, the explicit inclusion of diatoms in marine ecosystem models is required to reproduce the observed response to natural or purposeful iron fertilisation in the ocean (Aumont and Bopp, 2006), and observed changes in export production during glacial cycles (Bopp et al., 2002). The representation of diazotrophs (i.e. N 2 -fixing organisms) is necessary to simulate the feedbacks between iron and the nitrogen inventories of the ocean Moore and Doney, 2007) and to reproduce observed N : P ratios Deutsch, 2010, 2012), of coccolithophores to simulate large blooms of phytoplankton (i.e. chlorophyll) biomass (Gregg and Casey, 2007) and phytoplankton succession (Gregg et al., 2003), and of Phaeocystis to reproduce the ecosystem structure in the Southern Ocean (Wang and Moore, 2011).
Fewer studies have examined the role of different zooplankton PFTs in global ocean biogeochemistry, even though there are zooplankton physiological data sets (e.g. Hirst and Bunker, 2003;Straile, 1997). The simulation of phytoplankton biomass was improved in published studies when more mechanistic parameterisations of zooplankton dynamics constrained by observations were included in a global model (Buitenhuis et al., 2006. Similarly, the seasonal cycle of phytoplankton (Aita et al., 2003) and the open-ocean oxygen depletion (Bianchi et al., 2013) were improved when the influence of zooplankton vertical migration was included in global biogeochemical models. The choice of the grazing formulation in particular was found to influence phytoplankton diversity (Prowe et al., 2012;Vallina et al., 2014b) and the resulting food-web dynamics Vallina et al., 2014a), and to have implications for energy flow to higher trophic levels (Stock et al., 2014).
Zooplankton can influence the fate of exported materials through several processes, including grazing, repackaging of organic matter in faecal pellets, and the vertical migrations and transport of carbon and nutrients into the mesopelagic zone (e.g. Stemmann et al., 2000;Steinberg et al., 2008). Furthermore, there are important interactions among grazing, nutrient cycles, and environmental conditions, as was shown in studies based on regional models and observations in the equatorial Pacific (Landry et al., 1997;Price et al., 1994), North Pacific (Frost, 1991), the Atlantic (Daewel et al., 2014;Steinberg et al., 2012) and the Southern Ocean (Banse, 1995;Bishop and Wood, 2009). The importance of grazing was also highlighted during iron enrichment experiments (Henjes et al., 2007;Latasa et al., 2014), in part explaining why Biogeosciences, 13, 4111-4133, 2016 www.biogeosciences.net/13/4111/2016/ some experiments led to increased carbon export and others did not (Martin et al., 2013). Thus, a more explicit representation of different zooplankton PFTs in global models could provide important clues for the functioning of marine biogeochemistry.
Here, we present a new dynamic green ocean model with ten PFTs. The parameterisation of vital rates associated with these PFTs is based on an extensive synthesis of published information on growth rates and other relevant parameters. We use the model to examine a long-standing paradox in biological oceanography: the low phytoplankton biomass in the Southern Ocean despite the high concentrations of macronutrients. This has been attributed to lack of iron (Fe) because of the distance to continental dust sources (Geider and La Roche, 1994;Martin, 1990). Increases in phytoplankton biomass have been produced in more than a dozen open ocean iron fertilisation experiments (Boyd et al., 2007;Smetacek et al., 2012). The influx of Fe has been proposed as a driver for the drawdown of atmospheric CO 2 during glaciations Watson et al., 2000), and intentional Fe-fertilisation has been considered as a means to both geo-engineer climate (Rickels et al., 2012) and to sell carbon credits (Tollefson, 2012). However, ocean biogeochemistry models that explicitly include the effect of Fe-limitation on phytoplankton growth fail to reproduce the low Chl biomass observed during summer in the Southern Ocean (Aumont and Bopp, 2006;Dutkiewicz et al., 2005;Le Quéré et al., 2005;Moore et al., 2004). This raises the question of the relative control exerted by Fe-limitation on biomass vs. that exerted by the grazing pressure of zooplankton (Banse, 1996;Price et al., 1994) and more generally on the suitability of the current generation of models to explore ecosystem-climate interactions. Our study addresses this question directly.

Model description and development
The PlankTOM10 dynamic green ocean model is a global ocean biogeochemistry model that includes plankton ecosystem processes based on the representation of 10 PFTs and their interactions with the environment. Plank-TOM10 incorporates six autotrophic and four heterotrophic PFTs: picophytoplankton (pico-eukaryotes and non N 2fixing cyanobacteria such as Synechococcus and Prochlorococcus), N 2 -fixers (Trichodesmium and N 2 -fixing unicellular cyanobacteria), coccolithophores, mixed phytoplankton (e.g. autotrophic dinoflagellates and chrysophytes), diatoms, colonial Phaeocystis, bacteria (here used to subsume both heterotrophic Bacteria and Archaea), protozooplankton (e.g. heterotrophic flagellates and ciliates), mesozooplankton (predominantly copepods), and crustacean macrozooplankton (euphausiids, amphipods, and others, called "macrozooplankton" for simplicity; Fig. 1  and PlankTOM6 (bottom) marine ecosystem models. The arrows show grazing fluxes by protozooplankton (purple), mesozooplankton (red), and macrozooplankton (green). Only fluxes with weighing factors above 0.1 are shown (Table 3). ton are not included in the model. Diversity within groups is not considered, and the physiological parameters for each PFT are the same everywhere in the ocean, although some are dependent on environmental conditions (i.e. nutrients, light, food, temperature).
The current version of the PlankTOM10 model was developed from the model of Buitenhuis et al. (2013a), using the strategy for regrouping PFTs described by Le Quéré et al. (2005). It does not include new equations for growth and loss terms compared with previous versions of the Plank-TOM model, but it includes an additional trophic level in the zooplankton PFTs (i.e. macrozooplankton). Parameterisations are based on more data related to the vital rates of individual PFTs, where new information was available. Previous studies have shown that model results are highly sensitive to PFT growth rates (Buitenhuis et al., 2006, and considerable effort was made to constrain these rates using observations from Breitbarth (2005), Bissinger et al. (2008), Buitenhuis et al. (2008, Sarthou et al. (2005), Schoemann et al. (2005, Rivkin and Legendre (2001), Hirst and Bunker (2003), and .
The complete set of model equations and parameter values are provided in the Supplement. Here, we describe the elements that are most important for the analysis of the Southern Ocean and the strategy used to determine parameter values for PFT growth and loss processes.
PlankTOM10 simulates the growth of ten PFTs in response to environmental conditions. The PFT biomasses are produced by the model for each grid box based on the growth and loss term equations presented in the Supplement. The model includes three detrital pools: large and small particulate organic matter, and semi-labile dissolved organic matter. The sinking rate of large particles is based on the mineral (ballast) content of particles following Buitenhuis et al. (2001), while the sinking rate of small particles is constant at 3 m d −1 . The model includes full cycles of carbon (C), oxygen (O 2 ), and phosphorus (P), which are assimilated and released by biological processes at a constant ratio of 122 : 172 : 1 (Anderson and Sarmiento, 1994). Phytoplankton and particulate organic matter have a variable Fe / C ratio, while zooplankton and bacteria have a fixed ratio of 2 × 10 −6 , which is lower than the minimum phytoplankton Fe / C ratio (Schmidt et al., 1999). Zooplankton and bacteria relelase excess iron. The model also includes a full cycle of silica (Si) and calcite (CaCO 3 ) as in Maier-Reimer (1993), and simplified cycles for Fe and nitrogen (N). CO 2 and O 2 are exchanged with the atmosphere using the gas exchange formulation of Wanninkhof (1992). The Fe cycle is represented as in Aumont and Bopp (2006). Iron is deposited with dust particles using the monthly fields of Jickells et al. (2005), the Fe content of dust is assumed to be 3.5 % everywhere. We assume an Fe solubility from dust of 1 % (Jickells et al., 2005). Iron is also delivered to the ocean via river fluxes following the outflow scheme of da Cunha et al. (2007) with 95 % sedimentation in estuaries. Dissolved inorganic nitrogen (DIN) is the sum of nitrate and ammonium. The N : P ratio of organic processes is set to the Redfield ratio of 16 : 1. N 2 -fixers can use N 2 and thus have access to unlimited N from the atmosphere.
The growth rate parameters for the ten PFTs in Plank-TOM10 are based on a compilation of growth rates as a function of temperature (Sect. 2.2). Phytoplankton PFT growth rates are also limited by light and inorganic nutrients (P, N, Si, and Fe) using a dynamic photosynthesis model that represents the two-way interaction between photosynthetic performance and Fe / C and Chl / C ratios (Buitenhuis and Geider, 2010). Light limitation is constrained by the slope of the photosynthesis-irradiance curve (α) and the maximum Chl / C ratio (θ max ). We could not distinguish PFT-specific values for α (Geider et al., 1997) and used a mean value of 1.0 mol C m 2 (g Chl mol photons) −1 for all PFTs. Observed θ max for diatoms are systematically higher than those of other PFTs (Geider et al., 1997). There are too few direct observations to parameterize θ max for other PFTs, so we fitted the observations (Geider et al., 1997) for θ max to the maximum growth rate (µ max ) presented in that paper. The fit showed θ max increasing with growth rate (n = 19, p = 0.02). We thus used a θ max higher than average for Phaeocystis and diatoms, and a lower than average θ max for N 2 -fixers.
We used a two-step approach to define the nutrient limitation parameters, which are not well constrained by observa-  tions. Firstly, we assigned initial PFT-specific half-saturation values to each phytoplankton PFT based on literature-derived values, using the value for a similar-sized PFT when PFTspecific information was not available. We then examined the covariations of surface Chl concentrations with the limiting nutrient concentrations as shown in Fig. 3, and adjusted the magnitude of the half-saturation parameters of phytoplankton PFT to approximately fit the observations, keeping the ratios of k half values between phytoplankton PFTs approximately the same as the initial ratios. With this approach, we use the observed k half values as an initial starting point but tune the model to match the emerging properties highlighted in Fig. 3.
Initial values for the half-saturation concentrations of P (k P ) and N (k N ) for phytoplankton growth rates were based on observations. For N 2 -fixers, coccolithophores and diatoms, the half-saturation values for growth were computed using the half-saturation values of uptake reported in Riegman et al. (1998, and Sarthou et al. (2005) multiplied by the minimum/maximum N : C ratio Table 1. Growth rates of PFTs at 0 and 20 • C (µ 0 and µ 20 ) and rate increase for a 10 • C increase in temperature (Q 10 ). The uncertainty in µ 0 and Q 10 represents ±1 standard deviation from an optimal parameter value in the parameter space. Full references for the phytoplankton growth rate data are provided in the Supplement. The zooplankton growth rate data are from the published data synthesis cited here. (0.33) to account for the acclimation of nutrient saturated vs. nutrient limited growth (Morel, 1987). For picophytoplankton, reported values for the half saturation extend over 3 orders of magnitude. We assigned low half-saturation values as these organisms grow even under very low nutrient conditions . For mixed phytoplankton, we assigned a value intermediate between picophytoplankton and diatoms. For Phaeocystis, we used half-saturation values that characterise colonies (Schoemann et al., 2005). The selected set of parameter values shown in Fig. 3 is reported in Table 2. Iron uptake was computed using a cell quota model (Buitenhuis and Geider, 2010;Geider et al., 1997), where the Fe uptake by phytoplankton PFTs is explicitly regulated by the light conditions. The three parameters needed are the minimum, the maximum and the optimal Fe quotas. The minimum and maximum quotas were set at the same value of 2.5 and 20 µmol Fe / mol C for all PFTs based on the analysis of Buitenhuis and Geider (2010). The optimal quota was set to the minimum quota plus 2 * µ max 20 based on Sunda and Huntsman (1995) for all PFTs. In addition, phytoplankton PFTs also respond to the concentration of Fe in water which is parameterised with a half-saturation constant. The half saturation of Fe uptake (k Fe ) is lower for picophytoplankton  than other phytoplankton, and higher for N 2 -fixers (LaRoche and Breitbarth, 2005) and diatoms (Sarthou et al., 2005). Intermediate values for k Fe have been reported for the other phytoplankton PFTs (Le Vu, 2005;Schoemann et al., 2005). The selected set of parameter values after adjustments produces no systematic covariation between Chl and Fe, as observed (Fig. 3, Table 2).
The half-saturation parameters of zooplankton grazing rate were initially based on the relationship between metabolic rates and body volume of Hansen et al. (1997). We used the same approach as for nutrient limitation of the phytoplankton PFTs, and adjusted the half-saturation parameters for grazing based on the observed covariations between surface Chl concentrations and zooplankton biomass (Fig. 3). The selected set of parameter values that approximately fit the observed covariations in Fig. 3 is reported in Table 2.
Zooplankton food preferences were assigned based on predator-prey size ratio (Table 3), as there were insufficient data to determine these parameters directly across the range of zooplankton and phytoplankton considered here. This approach assumes that protozooplankton generally have a high preference for bacteria and a low preference for diatoms, that mesozooplankton have a higher preference for protozooplankton and a low preference for N 2 -fixers and bacteria, and macrozooplankton have a lower preference for N 2 -fixers, picophytoplankton and bacteria than other groups. Although some data were available to characterise grazing on Phaeocystis spp. (Nejstgaard et al., 2007), it is not used specifically here because it required knowledge on the life forms of Phaeocystis in situ. We assume that all zooplankton graze on organic particles (Table 3) but prefer to graze on other PFTs. The weighing factors influenced primarily the biomass of the prey and predators, but had little influence on their geographic distribution. We thus used the model results on biomass (Table 4) to guide the size of the relative preferences among PFTs for each grazer.
The gross growth efficiency (the part of grazing that is incorporated into biomass) was defined based on the mean across available observations: 0.21 for bacteria (data from    Rivkin and Legendre, 2001), and 0.29, 0.25, and 0.30 for protozooplankton, mesozooplankton and macrozooplankton, respectively (data from Straile, 1997  (3-47 %) * The biomass ranges have been computed using the method described in Buitenhuis et al. (2013b).
(Sect. 2.2), except for the mortality of macrozooplankton and mesozooplankton. There are nine observations on macrozooplankton mortality and we tuned this term based on the resulting biomass. The fitted relationship for the mortality of mesozooplankton was reduced by a factor of ∼ 2 to account for the explicit mortality from macrozooplankton represented in the model. This correction preserves the temperaturedependence of mortality, but it recognises that explicit grazing by macrozooplankton already takes place in the model, which does not represent the grazing by other organisms (e.g. salps, fish larvae). In total, grazing accounts for 2/3 to 3/4 of the mortality of mesozooplankton (Hirst and Kiorboe, 2002).

Growth rates as a function of temperature
The most important trait that distinguishes the various PFTs is the rate at which they grow under different conditions (Buitenhuis et al., 2006. We compiled maximum growth rates as a function of temperature (Table 1). We fit an exponential growth relationship to the observations by optimising the relation µ T = µ 0 × Q T /10 10 where T and µ T are the observed temperature and associated growth rate, µ 0 is the growth at 0 • C, and Q 10 is the derived temperaturedependence of growth (Table 1). The parameter values for µ 0 and Q 10 were estimated by minimising the error, quantified as the least squares cost function ((µ T − µ T obs )/µ T obs ) 2 . Normalising to observations helps ensure a good fit of µ T in cold waters where growth rates are low. We used exponential growth, rather than a temperature-optimal growth, to avoid biases caused by the lack of observations for some PFTs at low or high temperatures. The p-value of a linear regression between observations and the exponential fit (Table 1) provides a measure of how well the relationship is constrained by the observations. The fit assigns equal weight for all the data, rather than following the 99 % quantile (e.g. Eppley, 1972;Bissinger et al., 2008) to provide a better representation of the mean community for each PFT.  (Garcia, 2006b). Fe data are from (Tagliabue et al., 2012). The protozooplankton biomass data are updated from , the mesozooplankton biomass data from Buitenhuis et al. (2006), and the macrozooplankton biomass data include all krill data from Atkinson et al. (2004) and other crustacean data from . All data are monthly averages except for the mesozooplankton, which are seasonal. All data are for the surface, generally corresponding to the mixed layer, except for observed Chl, which is seen by satellite over one optical depth, and observed mesozooplankton and macrozooplankton, which are from depth-integrated tows and may underestimate surface concentrations (by a factor 1.5-2; see text). The black lines are medians, and grey shadings the 25-75 % interquartile range for Chl concentration. The median from the PlankTOM10 model is shown in red.
Growth rate parameters estimated with this method are well constrained (p values < 0.05) for seven of the ten PFTs, including all of the heterotrophic PFTs (Table 1). There are insufficient data to provide significant constraints on the growth rates of N 2 -fixers (p = 0.76), and some uncertainty in the growth data for coccolithophores (p = 0.06) and Phaeocystis (p = 0.23; Table 1). However, the growth of N 2 -fixers is less than that of other phytoplankton PFTs (Fig. 2), and the fitted relationship produces µ T less than that of other PFTs despite these uncertainties. An exponential function may not be appropriate for growth rates of coccolithophores and Phaeocystis (Schoemann et al., 2005). The growth rate of coccolithophores was overestimated at low temperatures due to high growth rates at 20 • C and the absence of observations for temperatures below 5 • C. We reduced the fitted growth rate of coccolithophores linearly to 0 below 10 • C to match the observed reduced coccolithophore biomass in cold regions (O'Brien et al., 2013).

Covariation between Chl and nutrients or zooplankton
We used relationships between observed concentrations of Chl and both inorganic nutrients (e.g. NO 3 , PO 4 and Fe), and zooplankton biomasses (protozooplankton, mesozooplankton and macrozooplankton; Fig. 3) to provide additional constraints on model parameters. Specifically, we used observations for in situ NO 3 and PO 4 concentrations from the World Ocean Atlas 2009; in situ Fe concentration data from Tagliabue et al. (2012); protozooplankton biomass data from ; mesozooplankton biomass data from Buitenhuis et al. (2006); macrozooplankton biomass data from Atkinson et al. (2004) and .
All the data were binned into 1 × 1 • grid boxes. Most observations are for the surface ocean. Mesozooplankton and macrozooplankton data are from depth-integrated tows of typically 200 m depth and may underestimate surface concentrations (by a factor 1.5-2 based on our model simulations). All data are monthly except for mesozooplankton, which are seasonal. Chl concentration is from SeaW-iFS satellite averaged over 1998-2009 and interpolated to the same grid. The model output was averaged over the same time period, and sampled for the same month and on the same grid box as the observations. The data intervals were chosen to include approximately the same number of grid boxes, except for macrozooplankton where the lowest interval was set to 0-0.05 µmol C L −1 because of the large number of grid boxes with very low macrozooplankton concentration. Ten concentration intervals were used for the nutrients (Fig. 3). Chlorophyll concentrations covary with NO 3 concentrations at < 3 µmol L −1 , and with PO 4 in the range 0.3-0.5 µmol L −1 (Fig. 3; Spearman ranked correlations for data in the 25-75 % interquartile range give r = 0.72 for NO 3 and r = 0.73 for PO 4 ). These relationships are consistent with our understanding of the growth limitation of phytoplankton in the subtropics, where NO 3 and PO 4 concentrations are low. There is no observed covariation between Chl and Fe concentration (r = −0.16). The strongest covariations are between Chl and protozooplankton at concentrations < 0.6 µmol C L −1 (r = 0.83) and between Chl and mesozooplankton at concentrations < 0.3 µmol C L −1 (r = 0.77). There is no covariation between Chl concentration and macrozooplankton biomass (r = −0.19; Fig. 3). We use these relationships to tune the growth limitation parameters in the model, so that the functional relationships between Chl and nutrients or zooplankton are close to the observed relationships overall.

Simulations
PlankTOM10 is coupled to the Ocean General Circulation Model (OGCM) NEMO version 3.1 (NEMOv3.1). We used the global configuration (Madec and Imbard, 1996), which has a resolution of 2 • of longitude and a mean resolution of 1.5 • of latitude, with enhanced resolution up to 0.3 • in the tropics and at high latitudes. The model resolves 30 vertical levels, with 10 m depth resolution in the upper 100 m. NEMOv3.1 calculates vertical diffusion explicitly and represents eddy mixing using the parameterisation of Gent and McWilliams (1990). The model thus generates its own mixed-layer dynamics and associated mixing based on local buoyancy fluxes and winds. NEMOv3.1 is coupled to a dynamic-thermodynamic sea-ice model (Timmermann et al., 2005). PlankTOM10 is initialised from observations of dissolved inorganic carbon (DIC) and alkalinity from Key et al. (2004), O 2 and nutrients from Garcia et al. (2006a, b), and temperature and salinity from the World Ocean Atlas 2005 Locarnini et al., 2006). Fe is initialised with a constant concentration of 0.6 nmol Fe L −1 north of 30 • S and 0.2 nmol Fe L −1 in the Southern Ocean, consistent with observations Tagliabue et al., 2012). The PFTs equilibrated within 3 years and were not influenced by initialisation. The model is forced by daily winds and precipitation from the ECMWF interim reanalysis (Simmons et al., 2006) from 1989 to 2009. Results for standard simulations are averaged over 1998-2009. A series of sensitivity tests are presented for the model parameters that influence the key results the most.
To understand the interaction pathways among ecosystems, biogeochemistry and climate, we developed a simplified version of the model that included only six PFTs (Plank-TOM6) (Fig. 1). PlankTOM6 is identical to PlankTOM10, except that the growth rates of N 2 -fixers, mixed phytoplankton, Phaeocystis, and macrozooplankton are zero, and the mortality of the mesozooplankton is increased to account for the lack of macrozooplankton predation until the point when primary production is at its maximum. Given the otherwise similar model structure, parameters, initialisation and simulation protocol, comparison of results from PlankTOM6 and PlankTOM10 provides information on the specific roles of zooplankton dynamics in the model.

Temperature and size -dependence of PFT growth rates
The data show systematic patterns in growth rates that differ among PFTs. The growth rates of all PFTs increase with increasing temperature, but not to the same extent (Fig. 2). The growth rate of phytoplankton PFTs increases with PFT size, from 0.15 d −1 for N 2 -fixers to 1.87 d −1 for Phaeocystis, and the growth rate of heterotrophic PFTs decreases with size, from 1.22 d −1 for bacteria to 0.19 d −1 for macrozooplankton (Table 1). The sign of the relationship between growth rate and size between phytoplankton PFTs is the opposite of the sign of this relationship within specific PFTs, including diatoms (Sarthou et al., 2005), picophytoplankton (Chen and Liu, 2010) and coccolithophores (Buitenhuis et al., 2008).

Ecosystem properties in the PlankTOM10 model
PlankTOM10 reproduces the main characteristics of observed surface Chl, with high concentrations in the high latitudes and low concentrations in the subtropics, higher Chl concentration in the Northern compared to the Southern Hemisphere, and in the South Atlantic compared to the South Pacific Ocean (Fig. 4). The global biogeochemical fluxes simulated by PlankTOM10 are generally below or at the low end of the range of observed values (Table 4), with global primary production of 42.4 PgC yr −1 , export production of 7.6 PgC yr −1 , export of CaCO 3 and SiO 2 of 0.4 PgC yr −1 and 2.9 PgSi yr −1 , respectively, and N 2 fixation of 165 TgN yr −1 . PlankTOM10 produces distinctive geographical distributions of carbon biomasses among PFTs (Fig. 5). About a third of the phytoplankton biomass occurs as picophytoplankton, followed in descending abundance by diatoms and Phaeocystis, mixed phytoplankton, coccolithophores and N 2 -fixers (Table 4). This distribution is broadly consistent with observations (Buitenhuis et al., 2013b), but the simulated phytoplankton biomass is generally on the low side of the observational range, which is consistent with the results of the global biogeochemical fluxes. The simulated biomass of coccolithophores is overestimated (i.e. 0.077 compared with 0.001-0.032 PgC), although CaCO 3 export is underestimated, suggesting either that the model calcification or aggregation rates are too low or that zooplankton calcifiers contribute significantly to CaCO 3 export.
The model underestimates bacterial biomass by a factor of 10 compared with observations. This possibly reflects the fact that the model only represents highly active bacteria and a substantial fraction of observed biomass is from low activity and ghost cells. The model underestimates protozooplankton by a factor of 1.5-5 (in absolute value) or 2-3 (as a fraction of total biomass value) compared to observations (Table 4). This discrepancy could be caused by the underestimation of bacterial biomass, as bacteria are an important source of food for protozooplankton. The simplified representation of the range of protozooplankton grazers in a single PFT representing both heterotrophic nanoflagellates and microzooplankton could also play a role. Simulated mesozooplankton biomass is only slightly below the observed range, while simulated macrozooplankton biomass is within the observed range, although the uncertainty here is large (0.010-0.64 PgC). Overall the balance is slightly skewed towards relatively more biomass than observed in the larger zooplankton www.biogeosciences.net/13/4111/2016/ Biogeosciences, 13, 4111-4133, 2016 (53 % compared to 3-47 %) compared to the smaller zooplankton groups (13 % compared to 27-31 %; Table 4). The geographic distribution of each simulated PFT is also distinctive (Figs. 6-7). Satellite data products indicate that small phytoplankton (picophytoplankton and N 2 -fixers) are generally dominant in the tropics, haptophytes (coccolithophores and Phaeocystis) in mid to high latitudes, and diatoms in high latitudes (Alvain et al., 2005;Brewin et al., 2010). The simulated phytoplankton distribution generally matches the distribution inferred from satellite normalised radiance (Fig. 6), except in the temperate zones where observations suggest a balance between picophytoplankton and haptophytes and the model shows a dominance of haptophytes. PlankTOM10 also reproduces the locations of blooms of colonial Phaeocystis and coccolithophores (Fig. 7). The simulated geographic distributions of zooplankton PFTs are particularly distinctive, with protozooplankton abundant in the tropics and subtropics, mesozooplankton at high latitudes of both hemispheres, and macrozooplankton with high biomass in the North Pacific and South Atlantic and along the coasts (Fig. 5).
The marine ecosystem as a whole appears to function realistically: mesozooplankton grazing on phytoplankton is somewhat overestimated relative to the 5.5 Pg yr estimated by Calbet (2001), so they have taken over the role of principal herbivores. Possibly the faster turnover rates of small copepods are overrepresented in the observational data on mesozooplankton, leading to a trophic position of mesozooplankton somewhat too low in the food chain. Export production, phytoplankton biomass and metazoan zooplankton biomass are realistic in the model, leading to realistic seasonal cycles, but the regenerated part of primary production is underestimated, concomitant with low protozooplankton biomass, which impacts the model on shorter timescales of days.

Comparison of PlankTOM6 and PlankTOM10
PlankTOM10 and PlankTOM6 generally produce similar results in surface Chl concentration, nutrient distribution, primary and export production (Fig. 8), except that PlankTOM6 fails to reproduce the observed low Chl concentration in summer in the Southern Ocean ( Fig. 4; Sect. 3.4). The over-all differences between the two models, quantified statistically using a Taylor distribution (Taylor, 2001), are less than 0.1 in either correlation or normalised standard deviation (Fig. 8). PlankTOM10 does slightly better than PlankTOM6 for the distribution of Chl, primary and export production, but slightly worse for the distribution of silica and nitrate, with similar performance for phosphate (Fig. 8). These differences are small in part because of the short duration of the simulations presented here (20 years), which allow equilibration of the ocean surface only. The models are generally sim-  and PlankTOM6 (in white) with observations. Chl, biomass and nutrient observations are as in Fig. 3. Export production is from (Schlitzer, 2004) and represents annual mean flux at 100 m. Primary production is from Buitenhuis et al. (2013a) and includes monthly mean values for the surface 300 m. The black dot shows the location where the model results should be if it was perfect and there were no errors in the observations. The distance from the black dot quantifies the performance of the model (Taylor, 2001).
ilar also in their representations of the distribution of biomass among phytoplankton PFTs, with most of biomass being in picophytoplankton in both models ( Fig. 9 and Table 4). However, PlankTOM6 allocates more biomass to protozooplankton compared to PlankTOM10, though PlankTOM6 is still at the low end of observed concentrations (Table 4). The failure of PlankTOM6 to reproduce the observed low Chl concentration in the Southern Ocean during summer is further highlighted in Fig. 10, which shows the seasonal cycle of mean Chl for the Northern Hemisphere and the Southern Ocean, where it is most pronounced. In PlankTOM6, the seasonal cycles in the north and south are very similar, with the slightly lower concentrations in the Southern Ocean during summer caused by a slightly deeper summertime mixed-layer depth (29 m compared to 19 m). In contrast, in PlankTOM10, the seasonal cycle of Chl in the south is smaller and concentrations are always below those in the north, as is the case for observations. As PlankTOM6 and PlankTOM10 have identical physical environments (including mixed-layer depth), the north-south differences are due to ecosystem structure. In the following sections, we focus our analysis on the model parameters that influence the low Chl concentration in the Southern Ocean the most.

Role of zooplankton dynamics for HNLC regions
The observed phytoplankton biomass, including the low Chl concentrations in high-nutrient low-chlorophyll (HNLC) regions, reflects the balance between phytoplankton growth and loss. Phytoplankton growth rates vary with temperature, light, and nutrient supply, whereas losses result mainly from grazing by zooplankton, respiration, cell death, sinking to depth, and dilution by vertical mixing. Any process that reduces the net rate of increase in phytoplankton biomass (i.e. differences between growth and loss) may lead to low residual Chl concentration. For example, Platt et al. (2003a) showed that deep mixing by wind dilutes Chl in the surface layer and reduces the average irradiance experienced by the phytoplankton. This results in low growth rate and demand for nitrate, the conditions generally observed in HNLC regions. Here we further examine the consequences of high zooplankton-mediated grazing losses.
We use the north / south ratio in surface Chl concentration as a metric to quantify model performance, focusing on the Pacific Ocean where the contrast between the Northern Hemisphere and the Southern Ocean is most pronounced. This metric is simple and easy to quantify with data (geographic locations: boxes in Fig. 4). Satellite observations indicate a north / south Chl ratio of 2.16 ± 0.35 (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) mean ± 2 SD of annual values). To ensure that the ratio is not affected by potential biases in the SeaWiFS Southern Ocean data (Johnson et al., 2013), we also used in situ data from the World Ocean Atlas, which indicates a similar north / south Chl ratio of 2.0. This ratio is 1.72 ± 0.051 in the PlankTOM10 and 1.21 ± 0.074 in the PlankTOM6 simulations (Fig. 11). Controlling factors in this ratio are examined here through a set of sensitivity tests.

Role of trophic level and top zooplankton
We tested the specific effect of macrozooplankton on Chl by running four additional model experiments (Fig. 11): in the Z1 simulation, we added macrozooplankton to PlankTOM6, in Z2 we parameterised the top grazer in PlankTOM6 using the same growth and loss rate parameters as macrozooplankton, in Z3 we removed macrozooplankton from Plank-TOM10, and in Z4 we parameterised the top grazer in Plank-TOM10 using the same growth and loss rate parameters as mesozooplankton. These sensitivity studies were identical to the PlankTOM10 (or PlankTOM6) simulation in all other respects. Experiments Z1 and Z2 both include macrozooplankton, but in different food-web positions. These experiments maintain a high north / south Chl ratio of 1.64 and 1.46, respectively (Fig. 11). Experiments Z3 and Z4 did not include macrozooplankton, but had grazing structures as in the standard PlankTOM6 and PlankTOM10 models; the north / south Chl ratio was 1.26 and 1.11, respectively. These four experiments show that the presence in the model of slow-growing zooplankton, such as macrozooplankton, plays a pivotal role in determining the relative average concentrations of Chl in the Northern vs. Southern Hemisphere (difference between PlankTOM6 and both Z1 and Z2). More realistic patterns are achieved by including a third zooplankton food-web compartment (higher ratio in Z1 than in Z2) and three additional phytoplankton compartments (higher ratio in PlankTOM10 than in Z1).

Role of macrozooplankton growth rate
We examined the impact of macrozooplankton grazing in sensitivity tests in which the grazing rate of macrozooplankton was varied within the range of the observed growth rates ( Fig. 2; Table 1). These simulations show that macrozooplankton grazing rate has a strong influence on the Chl north / south ratio (Fig. 12). The PlankTOM10 simulation that uses the mean growth rate from observations (Sect. 2.2) produces results that are closest to the observed north / south Chl ratio. When the grazing rate is decreased (by up to 2σ ), the macrozooplankton biomass decreases by over 50 % and the north / south Chl ratio decreases from 1.72 to 1.05. When the grazing rate is increased, the macrozooplankton biomass decreases because of pressure on the food sources (Fig. 12), and the Chl north / south ratio also decreases. These simulations suggest that the observed Chl north / south distributions are a consequence of trophic balances among PFTs.

Role of atmospheric iron deposition
We tested the relative role of atmospheric iron deposition compared with grazing for the north / south Chl distribution by applying five different dust deposition scenarios, all (except one) with realistic but different regional distributions, to the PlankTOM10 and PlankTOM6 models: D0 is an extreme case with no atmospheric dust deposition (where phytoplankton use iron sources from deep waters), D1 dust deposition including the effect of dust particle size on iron solubility , and D2-D4 iron deposition using the three distinct dust fields (Ginoux et al., 2001;and Luo, 2003;Tegen et al., 2004) averaged by Jickells et al. (2005). The simulated north / south Chl ratios vary from 1.62 to 1.85 in these experiments (Fig. 11). These differences are smaller than the differences between the PlankTOM10-like (1.46-1.85) and PlankTOM6-like simulations (1.08-1.26) for all experiments. In PlankTOM6, even the simulation with no iron deposition from dust (D0) produces Southern Ocean Chl concentrations that are too high during summer. This result is consistent with the observation that although Fe is lower in the Southern Ocean than elsewhere, concentrations average around 0.3 nmol Fe L −1 (range of 0.15-0.6 nmol Fe L −1 ) in the summer (January and February, n = 79) in the Subantarctic region (Tagliabue et al., 2012), which is near the half saturation for growth of most phytoplankton as well as those used in the model (Le Quéré et al., 2005;Sarthou et al., 2005). Thus Fe concentrations may be limiting for phyto-  (Arora et al., 2011;Dufresne et al., 2013;Giorgetta et al., 2013;Jones et al., 2011). Results from Séférian et al. (2013) mainly differ through their representation of sub-gridscale processes, with improvements in the representation of summer mixed-layer depth from Models 1 to 3. All data are averaged between 30 and 55 • latitude in both hemispheres: 140-240 • E in the north and 140-290 • E in the south, as highlighted in Fig. 4. plankton growth, but nevertheless the observed very low Chl concentrations during summer months seem to reflect losses due to other processes, such as grazing mortality rather than reduced growth rates from low Fe supply. As a means of validating the model results, we also tested the response of PlankTOM10 to Fe-fertilisation to verify that the model reproduced the observed Chl blooms under Fe enrichment conditions (Boyd and al., 2007). This was done by saturating the surface layer of the ocean with Fe for 1 month (February). In this experiment, surface Chl south of 40 • S increased by 2.1 ± 2.2 mg Chl m −3 (mean ± 1 SD) with a maximum concentration of 14.2 mg Chl m −3 . This is similar to the responses observed at sea during Fe-fertilisation experiments (Boyd and al., 2007). Thus Planktom10 predicts that net phytoplankton growth can escape the constraint imposed by zooplankton grazing and bloom when superabundant Fe is provided, as is the case during the meso-scale Fe-fertilisation experiments. The response of the model to Fe enrichment lends further support to our hypothesis that grazing is responsible for the low Chl concentration in the Southern Ocean during summer under realistic Fe inputs.

Role of combined effects
Model simulations could be influenced by the model structure and parameters, the physical transport, meteorological data, or the choice of dust deposition fields. We assessed the combined effects of model choices by comparing our results with outputs from seven other models: a version of the PISCES model (Aumont and Bopp, 2006), the CCSM-BECs model (Doney et al., 2009), and the NEMURO model (Kishi et al., 2007), IPSL-CM5A-LR (Dufresne et al., 2013), GRDL-ESM2M (Jones et al., 2011), HadGEM2-ES (Giorgetta et al., 2013, and CanESM2 (Arora et al., 2011). All of these other models focus on the representation of phytoplankton groups and parameterise grazing pathways in a simpler fashion than PlankTOM10. They produce a north / south Chl ratio in the range from 0.60 to 1.36, lower than the value (1.72) obtained using PlankTOM10. Previous studies have suggested that the overestimation of Chl may result from a generalised model bias towards too shallow a mixing depth in the Southern Ocean in summer, but Séférian et al. (2013) have shown that while better representation of sub-grid-scale processes and mixed-layer depth improves the simulation of Chl overall, it does not lead to a more realistic north / south Chl ratio (Fig. 11). Thus, the comparison between Plank-TOM10 and other ocean biogeochemistry models supports our contention that it is important to simulate grazing pathways explicitly.

Discussion
The development of PlankTOM10 has benefited from the existence of the very extensive range of observations to develop realistic parameterisations of key processes, particularly PFT growth rates. Although the simulated global biogeochemical fluxes are generally below or at the low end of the range of observed values and several regional discrepancies exist between observed and modelled biomass and fluxes, the model reproduces both the relative importance of different PFTs and the geographic patterns in their abundance. Thus, while not perfect, the model is sufficient to explore the role of ecosystem dynamics in determining ocean biogeochemistry.
Our analyses suggest that Southern Ocean Chl during summer is primarily controlled by zooplankton grazing, particularly the presence of a slow-growing zooplankton, and the structure of the pelagic food web, rather than the low supply rate of iron. Trophic cascading appears to account for the differences between the results from PlankTOM10 and PlankTOM6 ( Fig. 13; Zollner et al., 2009  ple, protozooplankton graze on phytoplankton (and bacteria), which reduces their prey's biomass. However, mesozooplankton graze on phytoplankton and protozooplankton, and macrozooplankton graze on phytoplankton and both protozooplankton and mesozooplankton. Thus the grazing pressure of larger zooplankton on smaller zooplankton can indi-rectly reduce the overall grazing pressure on phytoplankton. In PlankTOM10, macrozooplankton concentration is higher in winter in the Northern Hemisphere Pacific sector where the surface layer is more stratified and food is abundant, compared with the Southern Ocean Pacific sector where the surface layer is more mixed and food is scarce. Thus when the spring bloom starts in the north, the biomass and grazing pressure exerted by macrozooplankton is high enough to reduce the biomass of smaller zooplankton, consequently reducing the grazing pressure on Chl and leading to an increase in Chl. However, in the south, macrozooplankton biomass is too low to cause significant losses of smaller zooplankton. Hence, the high proto-and meso-zooplankton biomasses prevent a phytoplankton bloom from developing in that region. Although PlankTOM6 simulates some degree of trophic cascade with the presence of two zooplankton PFTs, our sensitivity tests presented in Fig. 11 show that the difference in growth rates between the two zooplankton PFTs is too small to impact the phytoplankton significantly. The higher concentration of macrozooplankton biomass in the north compared to the south is consistent with the observations, where the mean biomasses of macrozooplankton were reported to be 3 times higher in the Northern Hemisphere compared to the Southern Hemisphere . A similar contrast is found between the Atlantic and Pacific sectors of the Southern Ocean, where the high macrozooplankton biomass observed in the Atlantic (Atkinson et al., 2004) would reduce the abundance of smaller zooplankton, resulting in higher Chl concentrations in the Atlantic sector, as simulated in PlankTOM10 (Fig. 4). Such trophic cascades have been observed in diverse ecosystems on land and in the ocean (Casini et al., 2009). Furthermore, many observational-based studies have highlighted the important role of zooplankton grazing for controlling phytoplankton biomass (Atkinson et al., 2001;Banse, 1996;Dubischar and Bathmann, 1997;Granĺi et al., 1993). Although some processes are missing from the model (e.g. vertical migration of zooplankton, which mostly contributes to downward export), the model suggests that the primary cascading effect of grazing is sufficient to account for a large part of the north / south Chl differences.
Our results indicate that zooplankton grazing exerts an important control on Southern Ocean Chl. This propagates through to influence phytoplankton biomass. Indeed, the north / south ratio of phytoplankton biomass at the surface is greater in PlankTOM10 (1.62) compared to PlankTOM6 (1.18), very close to the modelled north / south ratio of Chl. The difference between PlankTOM10 and PlankTOM6 also persists through depth until about 300 m. Because of these marked differences, it is clear that the representation of global biogeochemical cycles in ocean models is influenced by the ecosystem structure. In both PlankTOM6 and Plank-TOM10, the mesozooplankton and macrozooplankton faecal pellets aggregate into the same large, fast-sinking particle pool, thus limiting the effect of different size classes of zooplankton on carbon export. To distinguish the effects of different food-web structures on export production, a wider spectrum of particle size classes sinking at different speeds are needed (e.g. Kriest, 2002). In addition, an improved vertical dynamics of the mesopelagic zone, together with the enhanced representation of zooplankton dynamics in the present study would allow further exploration of the interactions between iron fertilisation, grazing, and mixed-layer dynamics, which have led to large differences among ocean iron fertilisation experiments (Smetacek and Naqvi 2008;. There are a number of limitations to the current version of PlankTOM10, including simplified overwintering strategies for zooplankton, the use of a coarse Fe model, and the lack of representation of semi-refractory organic matter. In addition, the model does not include some ecosystem pathways, such as viral lysis (Evans et al., 2009), and the zooplankton representation does not include salps, pteropods, and auto-and mixo-trophic dinoflagellates. The nano-and micro-zooplankton are also combined into a single compartment. The realism of the simulations may also be affected by the relatively coarse resolution of the physical ocean model. However, these biases affect both PlankTOM6 and Plank-TOM10, and thus the experiments still provide information on the processes that differ between the two models. Our work suggests that improved representation of the zooplankton components could help further constrain the processes that regulate Chl distribution in models. The effect of further ecosystem model developments will be explored in followup studies.

Conclusions
The development of global marine ecosystem models is hampered in particular because of our poor understanding of several critical ecosystem processes and food-web interactions (Smetacek et al., 2004), and the paucity of global-scale observation of physiological rates and biomass for parameterisation and validation (Le Quéré and Pesant, 2008;Barton et al., 2013). For example, the wide range in observed growth rates for the same temperature is an indication of the challenges met by marine ecosystem modellers, particularly in representing the within-PFT diversity, which is unaccounted for in our model. In addition, the lack in knowledge of trophic relationships means that semi-arbitrary choices have to be made to characterise the predator-prey relationships based on size. Much more work is needed to understand the specific pathways by which matter circulates within ecosystems, taking into account the regional distributions of zooplankton groups and interactions with the environment including seasonal mixed-layer dynamics.
The role of macrozooplankton highlighted here has implications for carbon export to depth because faecal pellets of some macrozooplankton have very fast sinking rates (Fortier et al., 1994;Turner 2002). Hence, a more explicit representation of the pelagic food web in global models is needed to capture the full range of interactions between marine ecosystems, marine biogeochemistry and climate. The synthesis and analysis of observations and model results by the MAREDAT and MAREMIP projects provide valuable insights into the processes that control marine ecosystems, including the contributions that different PFTs make to ocean biomass (Buitenhuis et al., 2013a;Hashioka et al., 2013;Sailley et al., 2013).
Our simulations examining the effects of grazing on phytoplankton biomass raise questions about the biological and biogeochemical bases for the current projections of the feedbacks between climate (and other environmental changes) and marine ecosystems. It also highlights potential complications for the large-scale proposed use of purposeful Fe-fertilisation to enhance the deep ocean storage of CO 2 (Ciais et al., 2013). Assessments of the impact of such geo-engineering techniques will be unreliable, at least until the full ecosystem response including the grazing pathways (Landry et al., 1997) and the relationship between ecosystem dynamics and deep water carbon export (Smetacek et al., 2012) can be reproduced with models, which could be used to make quantitative predictions of deliberate Fe-fertilisation over large areas.
Our results on the important role of grazing do not contradict the results on the importance of Fe-fertilisation as highlighted in Fe enrichment experiments (Boyd and al., 2007), because additional Fe would trigger further growth provided that Fe were initially below an optimal concentration . However, our results suggest that low Fe concentrations by themselves are insufficient to account for the very low Chl levels observed in the Southern Ocean HNLC region in summer, and that differences in zooplankton trophic and community structure, and concomitant grazing dynamics play an important role in controlling phytoplankton blooms and maintaining very low Chl levels in that region. Although previous studies emphasised the role of phytoplankton community structure (Arrigo et al., 1999) and mixed-layer dynamics for nutrient supply and demand (Platt et al., 2003a, b) in ocean biogeochemical cycles, our analysis makes it clear that it is important to consider the whole pelagic ecosystem, including the zooplankton, when studying and predicting ecosystem responses to Fe (or any essential nutrient) fertilisation. This complex interplay has received less attention than either the drivers of primary production or the representation of Fe cycling in global biogeochemical modelling. Our results suggest that representing zooplankton interactions more explicitly in models would improve the representation of biogeochemistry-climate interactions, and could bring new insights to understand changing global biogeochemical cycles.
The Supplement related to this article is available online at doi:10.5194/bg-13-4111-2016-supplement.