PISCES-v 2 : an ocean biogeochemical model for carbon and ecosystem studies

PISCES-v2 (Pelagic Interactions Scheme for Carbon and Ecosystem Studies volume 2) is a biogeochemical model which simulates the lower trophic levels of marine ecosystems (phytoplankton, microzooplankton and mesozooplankton) and the biogeochemical cycles of carbon and of the main nutrients (P, N, Fe, and Si). The model is intended to be used for both regional and global configurations at high or low spatial resolutions as well as for short-term (seasonal, interannual) and long-term (climate change, paleoceanography) analyses. There are 24 prognostic variables (tracers) including two phytoplankton compartments (diatoms and nanophytoplankton), two zooplankton size classes (microzooplankton and mesozooplankton) and a description of the carbonate chemistry. Formulations in PISCES-v2 are based on a mixed Monod–quota formalism. On the one hand, stoichiometry of C /N /P is fixed and growth rate of phytoplankton is limited by the external availability in N, P and Si. On the other hand, the iron and silicon quotas are variable and the growth rate of phytoplankton is limited by the internal availability in Fe. Various parameterizations can be activated in PISCES-v2, setting, for instance, the complexity of iron chemistry or the description of particulate organic materials. So far, PISCES-v2 has been coupled to the Nucleus for European Modelling of the Ocean (NEMO) and Regional Ocean Modeling System (ROMS) systems. A full description of PISCES-v2 and of its optional functionalities is provided here. The results of a quasi-steady-state simulation are presented and evaluated against diverse observational and satellite-derived data. Finally, some of the new functionalities of PISCES-v2 are tested in a series of sensitivity experiments.


Introduction
Human activities have released large amounts of carbon into the atmosphere since the beginning of the industrial era leading to an increase in atmospheric CO 2 by more than 100 ppmv.The oceans play a major role in the carbon cycle and in its adjustment.Sabine et al. (2004) have estimated that the oceans have absorbed about one-third of the anthropogenic emissions.This role is tightly controlled by the physical and biogeochemical states of the marine system, i.e., by the characteristics of the solubility and biological pumps.Yet, the role played by the ocean in the carbon cycle is likely to be modified in response to climate and chemical changes induced by the anthropogenic carbon emissions (e.g., Orr et al., 2005;Steinacher et al., 2010;Bopp et al., 2013).Global ocean biogeochemical models represent powerful tools to study the carbon cycle and to predict its response to future and past climate and chemical changes.Since the pioneering work by Bacastow and Maier-Reimer (1990) based on a very simple description of the carbon cycle, the number and the complexity of models have rapidly increased (e.g., Six and Maier-Reimer, 1996;Moore et al., 2004;Quéré et al., 2005;Aumont and Bopp, 2006;Yool et al., 2011).However, Published by Copernicus Publications on behalf of the European Geosciences Union.

O. Aumont et al.: A description of PISCES-v2
a greater complexity of the models raises difficulties related to the lack of data for validation and to the theoretical justification of the parameterizations (e.g., Anderson, 2005Anderson, , 2010)).
PISCES (Pelagic Interactions Scheme for Carbon and Ecosystem Studies) is a biogeochemical model which simulates marine biological productivity and describes the biogeochemical cycles of carbon and of the main nutrients (P, N, Si, Fe).This model can be seen as one of the many Monod models (Monod, 1942) as opposed to the quota models (Mc-Carthy, 1980;Droop, 1983) which are alternative types of ocean biogeochemical models.Thus, it assumes a constant Redfield ratio, and phytoplankton growth depends on the external concentration in nutrients.This choice was dictated by the computing cost whereby the internal pools of the different elements (necessary for a quota model) requires many more prognostic variables.Ultimately, PISCES was assumed to be suited for a wide range of spatial and temporal scales, including, typically, several thousand year-long simulations on the global scale.
In contrast to the Monod approach, when modeling silicate, iron and/or chlorophyll, assuming constant ratios is not justified anymore as these ratios can vary substantially.For instance, the Fe / C ratio can vary by at least an order of magnitude, in particular as a result of luxury uptake, (e.g., Sunda andHuntsman, 1995, 1997) compared to the N / C ratio which varies by "only" 2 to 3 times.Equally, the Si / C ratio can vary significantly in response to the degree of iron stress (Hutchins and Bruland, 1998;Takeda, 1998).Thus, in PISCES, a compromise between the two classical types of ocean models was chosen.The Fe / C, Si / C and Chl / C internal ratios are prognostically predicted based on the external concentrations of the limiting nutrients as in the quota approach.Phytoplankton growth rates are predicted simultaneously using the Monod approach for N, P and Si and the quota approach for Fe.As a consequence, PISCES should be considered to be a mixed Monod-quota model.
Historically, the development of PISCES started in 1997 with the release of the P3ZD model which was a simple Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model with semi-labile dissolved organic matter (DOM) (Aumont, 1998;Aumont et al., 2002).Phytoplankton growth rate was only limited by one nutrient (effectively phosphate) and many shortcomings were apparent in this model, especially in the high nutrient-low chlorophyll (HNLC) regions.This served to justify the development, beginning in 1999, of a more complex model that includes three limiting nutrients (Fe, Si, P), two phytoplankton and two zooplankton size classes.This model was called HAMOCC5 (Aumont et al., 2003), as it was based on HAMOCC3.1 (Six and Maier-Reimer, 1996) and used in the LSG model (Maier-Reimer et al., 1993).When this code was embedded in the ocean model Ocean PArallélisé (OPA) (Madec et al., 1998), it required some major changes and improvements, partly because of the much finer vertical resolution.In addition to the numerical schemes, these changes were mostly an im-proved treatment of the optics and the separation of the particulate organic matter into two different size classes.All these changes and the major recodings it required led us to adopt a new name for the model: PISCES.This name can be translated as fishes from Latin.
PISCES has been used so far to address a wide range of scientific questions.Unfortunately, a complete list of the studies which have been based on or made use of PISCES is not available, but more than about hundred referenced studies explicitly rely directly or indirectly on this model.These range from process studies (Aumont and Bopp, 2006;Gehlen et al., 2006;Tagliabue et al., 2009a;Tagliabue and Völker, 2011) to operational oceanography (Brasseur et al., 2009).PISCES has been used to analyze intraseasonal (Gorgues et al., 2005;Resplandy et al., 2009) to interannual and decadal timescales (Raynaud et al., 2006;Rodgers et al., 2008).PISCES is part of the Institut Pierre et Simon Laplace (IPSL) and CEntre National de Recherche en Météorologie (CNRM) Earth system models which contribute to the different Intergovernmental Panel on Climate Change (IPCC)related activities including the Climate Model Intercomparison Project (CMIP5) modeling component (Séférian et al., 2013).Several studies have been conducted that consider the potential impact of climate change on ocean biogeochemistry (Dufresne et al., 2002;Bopp et al., 2005;Steinacher et al., 2010).Modeling studies focusing on paleoceanography have been based on PISCES (Bopp et al., 2003;Tagliabue et al., 2009b).Finally, PISCES is also used in regional configurations to study specific regions such as the Peru upwelling (Echevin et al., 2008;Albert et al., 2010) or the Indian Ocean (Resplandy et al., 2012).
However, PISCES-v2 is currently available only in the NEMO modeling system.The implementation of this updated version of PISCES in the ROMS_AGRIF modeling system is currently underway and should be finished and available by the end of 2015.Since 2001, PISCES has undergone active developments.In 2004, a stable release of the model was made available to the community on the OPA web site.Soon after, an earlier documentation of the model was published as a Supplement to the study by Aumont and Bopp (2006).Since then, the model has significantly evolved without any update of the documentation and this has effectively rendered the earlier documentation obsolete.After 6 years of intense developments, it is more than appropriate at this point to provide the current or future users of the model with an updated and accurate description of the current state of PISCES, called PISCES-v2.This paper describes the main aspects of the model.At its end, a description of a climatological simulation is proposed using the standard set of parameters available when the model is downloaded.Finally, the impact of several new parameterizations is evaluated through the performance of a set of sensitivity experiments.

Changes from previous release
As already mentioned, PISCES as a research tool is in perpetual evolution.Numerous changes have been made relative to the previously documented version, PISCES-v1.A brief list of the main changes is made below, with these changes organized thematically.These changes are detailed in the following sections.
-Changes made to the code structure and design: a. Transition to full native Fortran 90 coding.The model has also undergone a reorganization of its architecture and coding conventions following the evolution of NEMO.
b. I/O interface should now be set by default to IOM (the new input-output manager of the NEMO modeling system) to benefit from the major improvements this interface offers.
c. Memory and performance improvements have been made.This version should run slightly faster and take much less memory than v1.
d.The namelist now includes many more parameters that may thus be changed without recompiling the code.
-Changes made to the nutrients: a. iron chemistry can be described according to two different parameterizations: the simple old chemistry scheme based on one ligand and one inorganic species, and a new complex chemistry module based on five iron species and two ligands.
b. Scavenging of inorganic iron and coagulation of iron colloids have been redesigned.
-Changes made to the phytoplankton compartments: a. Nutrients limitation terms now include a simple description of the impact of cell size.
b. Iron content and growth rate limitation by iron is modeled following the quota formalism.Luxury uptake of iron can be represented by this new formulation.
c. Silicification, calcification as well as nitrogen fixation are redesigned by diazotrophs.
d.The relationship between growth rate (primary production) and light can be chosen between two different formulations.
-Changes made to the zooplankton compartments: a.The microzooplankton grazing formulation is now identical to that of mesozooplankton.
b. Thresholds can be selected for both total food or individual prey types.
c. Food quality affects the gross growth efficiency of both zooplankton compartments.
-Changes made to dissolved organic matter and particulate materials: a. Two different schemes for the description of particulate organic matter can be chosen: the traditional two-compartment model or the Kriest model.
b. Bacterial implicit description has been redesigned.
c. Dissolution of biogenic silica assumes two different fractions.
d.The dust distribution in the water column is modeled using a very crude parameterization.
e.The numerics of vertical sedimentation have been improved (time splitting scheme).
-Changes made to the external sources of nutrients and to the treatment of the bottom of the water column: a. Spatially variable solubility of iron in dust can be specified from a file.
b. River discharge of nutrients has been improved.c.Denitrification in sediments is now parameterized as well as variable preservation of calcite.
As a consequence of these changes, the user should be warned that results produced with PISCES-v1 cannot be reproduced by PISCES-v2.Furthermore, in the rest of this work, PISCES will designate PISCES-v2.

Model description
PISCES currently has 24 compartments (see Fig. 1).There are five modeled limiting nutrients for phytoplankton growth: nitrate and ammonium, phosphate, silicate and iron.It should be mentioned that phosphate and nitrate + ammonium are not really independent nutrients in PISCES.They are linked by a constant and identical Redfield ratio in all the modeled organic compartments, but the nitrogen pool undergoes nitrogen fixation and denitrification in the open ocean and the upper sediments.Furthermore, their external sources (rivers, dust deposition) are not linked by a constant ratio.This means that if the latter three processes (nitrogen fixation, denitrification, and external sources) are deactivated and if the initial distributions of nitrate + ammonium and phosphate are identical, the simulated fields of both nutrients should remain identical.Four living compartments are represented: two phytoplankton size classes/groups corresponding to nanophytoplankton and diatoms, and two zooplankton size classes which are microzooplankton and mesozooplankton.For phytoplankton, the prognostic variables are the carbon, iron, chlorophyll and silicon biomasses (the latter only for diatoms).This means that the Fe / C and Chl / C ratios of both phytoplankton groups as well as the Si / C ratio of diatoms are prognostically predicted by the model.For zooplankton, only the total biomass is modeled.For all species, the C / N / P / O 2 ratios are assumed constant and are not allowed to vary.In PISCES, the Redfield ratios C / N / P are set to 122/16/1 (Takahashi et al., 1985) and the −O / C ratio is set to 1.34 (Kortzinger et al., 2001).In addition, the Fe / C ratio of both zooplankton groups is kept constant.No silicified zooplankton is assumed.The bacterial pool is not yet explicitly modeled.
There are three non-living compartments: semi-labile dissolved organic matter, small sinking particles and large sinking particles.As for the living compartments, the C, N and P pools are not distinctly modeled.Thus, constant Redfield ratios are imposed for C / N / P. On the other hand, the iron, silicon and calcite pools of the particles are explicitly modeled.As a consequence, their ratios are allowed to vary.The sinking speed of the particles is not altered by their content in calcite and biogenic silicate ("the ballast effect", Honjo, 1996;Armstrong et al., 2002).The latter particles are assumed to sink at the same speed as the large organic matter particles.An earlier version of PISCES had included a simple description of this ballast effect (Gehlen et al., 2006) but it has been abandoned since as observations do not suggest a clear relationship between sinking speeds and mineral composition of particles (Lee et al., 2009).All the non-living compartments experience aggregation due to turbulence and differential settling as well as Brownian coagulation for DOM.
In addition to the ecosystem model, PISCES also simulates dissolved inorganic carbon, total alkalinity and dissolved oxygen.The latter tracer is also used to define the regions where oxic or anoxic degradation processes take place.

Model equations
The reader should be aware that in the following equations, the conversion ratios between the different elements (Redfield ratios) have been generally omitted except when particular parameterizations are defined.All phytoplankton and zooplankton biomasses are in carbon units (mol C L −1 ) except for the silicon, chlorophyll and iron content of phytoplankton, which are respectively in Si, Chl and Fe units (mol Si L −1 , g Chl L −1 , and mol Fe L −1 , respectively).Finally, all parameters and their standard values in PISCES are listed in Tables 1a-e at the end of this section.

∂P ∂t
= (1 − δ P )µ P P − m P P K m + P P − sh × w P P 2 In this equation, P is the nanophytoplankton biomass, and the five terms on the right-hand side represent growth, mortality, aggregation and grazing by micro-and mesozooplankton.The mortality term is modulated by a hyperbolic function of P to avoid extinction of nanophytoplankton at very low growth rates.
In PISCES, the growth rate of nanophytoplankton (µ P ) can be computed according to two different parameterizations:   terms in these equations are defined below.The choice between the two different formulations is made through a parameter in the namelist (ln_newprod).When ln_newprod is set to true, which is the default option of PISCES, Eq. ( 2a) is used.In the previous equations, L day is day length (∈ [0, 1]).f 1 (L day ) expresses the dependency of growth rate to the length of the day (Gilstad and Sakshaug, 1990;Thompson, 1999).z mxl is the depth of the mixed layer and f 2 (z mxl ) imposes an additional reduction of the growth rate when the mixed layer depth exceeds the euphotic depth: where z eu is the depth of the euphotic zone defined as the depth at which there is 1 ‰ of surface PAR.t P dark is set to 3 days for nanophytoplankton and 4 days for diatoms, as diatoms generally better cope with prolonged dark periods.t dark is an estimate of the mean residence time of the phytoplankton cells within the unlit part of the mixed layer, assuming a vertical diffusion coefficient of 1 m 2 s −1 ; 86 400 converts t dark from s −1 to day −1 .Figure 2 displays f 2 (z mxl ) as a function of z.
µ P is defined as follows (Eppley, 1972): In PISCES, vertical penetration of the photosynthetic available radiation (PAR) is based on a simplified version of the model by Morel (1988), which is described in Lengaigne et al. (2007).Visible light is split into three wavebands: blue (400-500 nm), green (500-600 nm) and red (600-700 nm).
For each waveband, the chlorophyll-dependent attenuation coefficients are fitted to the coefficients computed from the full spectral model of Morel (1988) (as modified in Morel and Maritorena, 2001) assuming the same power-law expression.At the sea surface, visible light is split equally between the three wavebands.PAR can be a constant or a variable fraction of the downwelling shortwave radiation, as specified in the namelist (ln_varpar).
Light absorption by phytoplankton depends on the waveband and on the species.The normalized coefficients β i have been computed for each phytoplankton group by averaging and normalizing, for each waveband, the absorption coefficients published in Bricaud et al. (1995).
In PISCES, the nutrient limitation terms are defined as follows: L P Fe = min 1, max 0, θ Fe,P − θ Fe,P min θ Fe,P opt .

(6f)
As already stated in the introduction, PISCES is a mixed Monod-quota model.Thus, N and P limitations are based on a Monod parameterization where growth depends on the external nutrient concentrations, whereas Fe limitation is modeled according to a classical quota approach.It should be noted here that for iron, an optimal quota (θ Fe,P opt ) is used in the denominator which allows luxury uptake as in the model proposed by Buitenhuis and Geider (2010).
The choice of the half-saturation constants is rather difficult as observations show that they can vary by several or- ders of magnitude (e.g., Perry, 1976;Sommer, 1986;Donald et al., 1997).However, in general, these constants increase with the size of the phytoplankton cell as a consequence of a smaller surface-to-volume ratio (diffusive hypothesis) (Eppley et al., 1969).Thus, diatoms will tend to have larger half-saturation constants than nanophytoplankton.However, in PISCES, phytoplankton are modeled by only two compartments, each of them encompassing a large range.Experiments performed with the model have shown that results are sensitive to the choice of these half-saturation constants.
Following these remarks, it appeared not appropriate to keep the nutrient half-saturation constants constant.It was then decided to make them vary with the phytoplankton biomass of each compartment because the observations show that the increase in biomass is generally due to the addition of larger size classes of phytoplankton (e.g., Raimbault et al., 1988;Armstrong, 1994;Hurtt and Armstrong, 1996): where S P rat is the size ratio of the larger size class over the smaller size class.K P ,min i is the half-saturation constant of the smaller size class.This parameterization assumes that half-saturation constants increase linearly with size (Eppley et al., 1969).The size dependence of these constants with cell size is not necessarily linear but has been suggested to follow a power-law function with an exponent lower than 1 (Litchman et al., 2007).However, in a recent review, Edwards et al. (2012) found an exponent close to 1 for nitrogen (linear relationship) and larger than 1 for phosphorus.Thus, considering these uncertainties, we decided to keep a linear relationship, which remains within the estimated range and which can also be derived from simple volumetric considerations (surface-to-volume ratio).The three parameters in this equation (P max , K P ,min i , and S P rat ) can be independently specified for each phytoplankton group.Finally, observations also suggest that these half-saturation constants should vary with the mean nutrient concentrations, probably as an acclimation to the local environment (Collos et al., 1980;Smith et al., 2009).This acclimation mechanism is not included in PISCES, except for the case of silicate (see Sect. 4.1.2).
The distinction between new production based on nitrate and regenerated production based on ammonium is computed as follows (O'Neill et al., 1989): where µ P NO 3 and µ P NH 4 are the uptake rates of nitrate and ammonium, respectively.
The nanophytoplankton aggregation term w P depends on the shear rate (sh), as the main driver of aggregation is the local turbulence.This shear rate is set to 1 s −1 in the mixed layer and to 0.01 s −1 below.This means that the aggregation is reduced by a factor of 100 below the mixed layer.

Diatoms
In this equation, D is the nanophytoplankton biomass, and the five terms on the right-hand side represent growth, mortality, aggregation and grazing by micro-and mesozooplankton.
As for nanophytoplankton, the absorption coefficients of diatoms depend on the considered waveband: The production terms for diatoms are defined as for nanophytoplankton, except that the limitation terms also include Si: As for the other nutrients, the half-saturation factor of silicate can vary significantly over the ocean.In the tropical and temperate regions, this factor is around 1 µM, whereas values as high as 88.7 µM have been measured for Antarctic species (Sommer, 1986;Martin-Jézéquel et al., 2000).In that case, rather than an effect of the cell size, these variations are a consequence of an acclimation of the cells to their local environment.When plotted against maximum local yearly concentration of silicate, a crude relationship can be inferred (Pondaven et al., 1998): where Si here is the maximum Si concentration over a year (note that during the first year of a pluriannual simulation, Si is set to a constant).For the other nutrients, we use the same parameterization as for nanophytoplankton (see Eq. 7).
The diatoms aggregation term w D p is increased in case of nutrient limitation because it has been shown that diatoms cells tend to excrete a mucus (exocellular polysaccharides, EPS) which increases their stickiness.As a consequence, collisions between cells yield to a more efficient aggregation process (Smetacek, 1985;Decho, 1990): Furthermore, as for nanophytoplankton, the aggregation is multiplied by the shear rate.Enhanced aggregation rates when diatoms are stressed result in a rapid decline of the diatoms blooms when nutrients become exhausted and produce strong export events.

Chlorophyll in nanophytoplankton and diatoms
Chlorophyll biomass I Chl (where I denotes P or D, typical units are µg Chl L −1 or mg Chl m −3 ) for both phytoplankton groups is parameterized using the photo-adaptive model of Geider et al. (1997): where I is the phytoplankton group and θ Chl,I is the chlorophyll-to-carbon ratio of the considered phytoplankton class; 12 represents the molar mass of carbon; ρ I Chl the ratio of energy assimilated to energy absorbed as defined by Geider et al. (1996): In this equation, 144 is the square of the molar mass of C and is used to convert from mol to mg, as the standard unit for Chl is generally in mg Chl m −3 .It should be noted that for chlorophyll synthesis, the second parameterization of phytoplankton growth is used to compute μI (see Eq. 2b).This is necessary because of the expression for ρ I Chl .

Iron in nanophytoplankton and diatoms
The temporal evolution of the iron biomass of phytoplankton I Fe (model units are mol Fe L −1 ), where I denotes P or D, is driven by the following equation Iron in phytoplankton is modeled in PISCES according to a classical quota approach.However, to be consistent with chlorophyll and silica, we model the iron biomass of phytoplankton (I Fe ) rather than the iron quota (θ Fe,I ) directly.Growth rate of the iron biomass of phytoplankton is parameterized according to As in Flynn and Hipkin (1999), iron uptake is also downregulated via a feedback from θ Fe,I using a normalized inverse hyperbolic function with a small shape factor set to 0.05.
In the former equation, L I Fe lim,1 is the iron limitation term and is modeled as follows: where bFe is the concentration of bioavailable iron (see Sect. 4.5.3).The half-saturation constant for iron uptake is also increasing with phytoplankton biomass as for the other half-saturation constants (see Eq. 7).At low iron concentrations, observations suggest that iron uptake might be enhanced, at least for some species (Harrison and Morel, 1986;Doucette and Harrison, 1991), giving surge uptake.Morel (1987) proposed a parameterization of both this surge uptake and the downregulation of iron uptake at high iron quota (see above) which has been included in the recent model of Buitenhuis and Geider (2010).In PISCES, a different parameterization has been chosen since downregulation is already included in Eq. ( 17): L lim,2 equals 4 at very low iron concentrations and 1 at high iron concentration.Overall, the downregulation in Eq. ( 17) together with the surge uptake induced by the previous equation results in a behavior of the system that is qualitatively equivalent to what results from the parameterization of Buitenhuis and Geider (2010).

O. Aumont et al.: A description of PISCES-v2
The demands for iron in phytoplankton are for photosynthesis, respiration and nitrate/nitrite reduction.Following Flynn and Hipkin (1999), we assume that the rate of synthesis by the cell of new components requiring iron is given by the difference between the iron quota and the sum of the iron required by these three sources of demand, which we defined as the actual minimum iron quota: In this equation, the first right term corresponds to photosynthesis, the second term corresponds to respiration and the third term estimates nitrate and nitrite reduction.The parameters used in this equation are directly taken from Flynn and Hipkin (1999).The modeled iron quota in PISCES varies thus between this minimum quota θ Fe,I min and the maximum quota θ Fe,I max , i.e., between about 1 and 40 µmol Fe (mol C) −1 when using the standard set of parameters (see Table 1).

Silicon in diatoms
The elemental ratio Si / C (or Si / N) has been observed to vary by a factor of about 4 to 5 over the global ocean with a mean value around 0.14 ± 0.13 mol mol −1 (Sarthou et al., 2005).Light, N, P, or Fe stress has been demonstrated to lead to heavier silicification (e.g., Takeda, 1998;Franck et al., 2000;Martin-Jézéquel et al., 2000).It has been suggested that these elevated elemental ratios result from the physiological adaptation of the silicon uptake by the cell depending on the growth rate and on the G2 cycle phase during which Si is incorporated (Martin-Jézéquel et al., 2000;Claquin et al., 2002).Lighter silicification can only result from silicate limitation.
We model the variations of the Si / C ratio following the parameterization proposed by Bucciarelli et al. (2002, unpublished manuscript): Relative to the original parameterization, an additional limitation term by Si has been added (F D Si lim,2 ) to produce a lighter silicification in case of Si exhaustion.
The different terms in Eq. ( 22) are defined as follows: where ϕ is the latitude.In the Southern Ocean, observations show that diatoms are very heavily silicified.After correcting for the potential effects of iron limitation, silicification in the Southern Ocean is at least 3 times stronger than in the tropical regions, which can only be explained by the diatoms morphological types (Baines et al., 2010).To reproduce those high Si / C ratios, we have introduced the term L D Si lim,2 which increases the Si / C ratio by a factor of up to 3 when silicate concentrations are high, a specific characteristics of the Southern Ocean.This increase is restricted to the Southern Hemisphere and is controlled by the parameter K 2 Si .This parameter is set in the namelist and thus, if it is set to a very high value, then no increase of Si / C at high silicate concentrations is predicted by the model.

Microzooplankton
In this equation, Z is the microzooplankton biomass, and the four terms on the right-hand side represent growth, grazing by mesozooplankton, quadratic and linear mortalities.
The grazing rate depends on temperature according to a typical exponential relationship similar to what is used for phytoplankton: where g 0,Z max is the maximum grazing rate at 0 • C, b Z is the temperature dependence and T is the temperature.In their review, Buitenhuis et al. (2010) have found a Q 10 (Q 10 = b 10 Z ) between 1.7 and 2.2.Lower temperature dependences were found in laboratory experiments compared to what as been identified in the field.In PISCES, we have set Q 10 to 2.14 which is not only close to the value found in the field but also close to the value chosen for mesozooplankton (see below).All terms driving the temporal evolution of microzooplankton have been assigned the same temperature dependence.Mortality is enhanced when oxygen is depleted.In other words, microzooplankton (but also mesozooplankton, see below) are treated as being unable to cope with anoxic waters.This increased mortality also avoids respiration in waters devoid of oxygen.
Grazing on each species I is defined as where J denotes all the species microzooplankton can graze upon (P , D and POC) and p Z J is the preference microzooplankton has for each J .In PISCES, we have chosen a Michaelis-Menten parameterization with no switching and a threshold (F Z thresh ) (Gentleman et al., 2003).This choice is rather arbitrary.Another very popular formulation in models is the Michaelis-Menten parameterization with active switching introduced by Fasham et al. (1990).However, this parameterization exhibits anomalous dynamics such as suboptimal feeding (Gentleman et al., 2003).In our parameterization, a threshold for each individual resource (J Z thresh ) can be specified in addition to the global threshold (F Z thresh ).For low food abundance, this global threshold is allowed to slowly decrease to 0 as a function of the total food level to maintain some grazing pressure, in particular in the ocean interior.
Responses of zooplankton to quality of their preys have been termed stoichiometric modulation of predation (SMP) by Mitra and Flynn (2005).A complete review of the different expected responses has been presented by Mitra et al. (2007).For instance, when confronted with poor food quality, zooplankton can increase their ingestion rate (Plath and Boersma, 2001;Darchambeau and Thys, 2005), or decrease it as the food can become deleterious (Flynn and Davidson, 1993).Accounting for the complexities of these different types of behavior has not been implemented within PISCES as this would require a model with flexible stoichiometry.Additionally, it would require a correct parameterization of the different potential responses and the apparently contradictory nature of observed responses implies that this task will be very complicated.In PISCES, food quality is assumed to only affect gross growth efficiency (e Z ): When food quality becomes poor (either the Fe / C ratio θ Fe,I or the N / C ratio θ N,I of the preys decreases), e Z decreases: When the Fe / C ratio of the ingested preys becomes lower than the zooplankton Fe / C ratio, the excess carbon (and nutrients) is lost as dissolved inorganic and organic carbon (and nutrients).This is described in PISCES by a decrease in the carbon gross growth efficiency (Eq.27b).By construction in PISCES, the N / C quota is constant, so this quota is estimated by solving the classical Droop equation assuming that it is at steady state (see above the definition of θ N,I ).

Mesozooplankton
In this equation, M is the mesozooplankton biomass, and the three terms on the right-hand side represent growth, quadratic and linear mortalities.All terms in this equation have been assigned the same temperature dependence using a Q 10 of 2.14 (Buitenhuis et al., 2005).
Parameterization of mesozooplankton grazing is similar to microzooplankton.In addition to the "conventional" concentration-dependent grazing described by Eq. ( 26a), flux feeding is also accounted for in PISCES.This type of grazing has been shown to be potentially very important for the fate of particles in the water column below the euphotic zone (Dilling and Alldredge, 2000;Stemmann et al., 2004).Flux feeding depends on the flux and thus, on the product of the concentration by the sinking speed.In PISCES, both the small and the large particles experience this type of grazing: This importance of flux feeding has been analyzed in PISCES by Gehlen et al. (2006).They have shown that flux feeding is the most important process that controls the flux of particulate organic carbon below the surface mixed layer.
In Eq. ( 28), the term with a quadratic dependency to mesozooplankton does not depict aggregation but grazing by the higher, non-resolved trophic levels.Following Anderson et al. (2013), the upper trophic levels are modeled assuming an infinite chain of carnivores.This assumption permits one to easily compute the production of fecal pellets as well as the respiration and excretion by these non-resolved carnivores: where function f up (x) is It should be noted here that a similar quadratic term is also included in the equation for microzooplankton (see Eq. 26a) despite the fact that their predators are (at least partially) represented in PISCES.In that case, this term rather represents other density-dependent mortality factors such as viral diseases.As a consequence, the assumption of an infinite chain of carnivores is not used for microzooplankton and everything is routed to POC.

DOC
The temporal evolution of DOC is driven by the following equation where I includes P , D and POC for microzooplankton and P , D, Z and POC for mesozooplankton (see Eqs. 24 and 28, respectively).In the following, DOM and DOC will be used indifferently since the stoichiometric ratios in dissolved organic matter are assumed constant in PISCES.Marine DOM has traditionally been divided into several fractions characterized by their lability.DOM, which recycles over timescales of a few months to a few years, is called semi-labile DOM (Anderson and Williams, 1999).Transport of this pool of dissolved organic matter can make a significant part of the carbon pump (Carlson et al., 1994;Anderson and Williams, 1999).As a consequence, this important pool of DOM is modeled in PISCES.The labile and refractory pools of DOM are not explicitly modeled.
The degradation of semi-labile DOC is parameterized as follows: Remineralization of DOC can be either oxic (Remin) or anoxic (Denit) depending on the local oxygen concentration.The distinction between the two types of organic matter degradation is performed using a factor (O 2 ) that varies between 0 and 1 (see Sect. 4.5.1 for the formulation of this factor).It is assumed that the specific rates of degradation (λ DOC ) specified for respiration and denitrification are identical.
Depending on the quality of the organic matter, bacteria may take up nutrients from seawater (e.g., Goldman and Dennett, 1991;Thingstad and Lignell, 1997), and thus may be limited by their availability.Of course, bacterial production is also limited by the abundance of dissolved organic matter.Therefore, we parameterize the regulation of the degradation of DOM by bacterial activity (L bact ) according to The half-saturation constants of the P and N limitation terms (K bact i ) are set in the namelist.In PISCES, bacterial biomass is not explicitly modeled; Instead, we use the following formulation: Otherwise. (35b) In the previous equation, 0.7(Z + 2M) is a proxy for the bacterial concentration.This relationship has been constructed from an unpublished version of PISCES (already mentioned in Aumont and Bopp, 2006) that includes an explicit description of the bacterial biomass.Below a certain depth (z max ), this biomass decreases with depth via a powerlaw function (Aristegui et al., 2009).
In Eq. ( 32), the terms DOC denote aggregation processes and are described hereafter (see Sect. 4.4.1).For DOM, we consider turbulence-induced as well as Brownian aggregation processes.

Particulate organic matter
PISCES includes two different schemes for particulate organic matter: -A simple model based on two different size classes for particulate organic matter.In that case, particulate organic matter is modeled in PISCES using two tracers corresponding to the two size classes: POC for the smaller class (1-100 µm) and GOC for the larger class (100-5000 µm).
-A more complex model proposed by Kriest and Evans (1999) in which the size spectrum of the particulate organic matter can be represented by a power-law function.Here, particulate organic matter is represented by two variables: the first (POC) is the carbon concentration and the second (NUM) is the total number of aggregates by unit volume of water.
By default, the simplest parameterization is used.The Kriest model is activated by a cpp key key_kriest.

Two-compartment model of Particulate Organic Matter (POM)
The temporal evolution of POC is written: where w POC is the vertical sinking speed.For POC, it is set to a constant value, in general to a small value on the order of a few meters per day.The fate of mortality and aggregation of nanophytoplankton depends on the proportion of the calcifying organisms (R CaCO 3 ).We assume that 50 % of the organic matter of the calcifiers is associated with the shell.Since calcite is significantly denser than organic matter, 50 % of the biomass of the dying calcifiers is routed to the fast sinking particles.The same is assumed for the mortality of diatoms as a consequence of the denser density of biogenic silica.The specific degradation rate λ POC depends on temperature with a Q 10 of about 1.9, the same as for phytoplankton.Furthermore, observations generally tend to show slower degradation rates when waters are anoxic (Harvey et al., 1995;Mooy et al., 2002).In Mooy et al. (2002), the attenuation coefficient (b) for the flux was found to be about 0.4 instead of the standard value 0.86 (Martin et al., 1987).This corresponds to a 45 % decrease of the degradation rate in anoxic waters relative to oxic waters, which is implemented as POC experiences aggregation due to turbulence and differential settling: In this equation, the first two terms correspond to turbulent aggregation, and the two last terms to differential settling aggregation.The values of the parameters controlling these processes have been computed offline assuming a steadystate power-law size spectrum for particles with an exponent of 3.6.Subsequently, the different coagulation kernels (e.g., Jackson, 1990;Kriest and Evans, 1999) have been integrated over the size ranges corresponding to the different compartments.A constant stickiness of 0.1 has been chosen.
The temporal evolution of GOC is written The equation controlling the temporal evolution of GOC is similar to that of POC.However, some observations have shown that the mean sinking speed of particulate organic matter increases with depth (e.g., Berelson, 2002).Such an increase is consistent with the power-law formulation proposed by Martin et al. (1987).Such an increase in the settling speed is parameterized in PISCES for GOC as follows: The parameters in this equation have been adjusted using a model of aggregation/disaggregation with multiple size classes (Gehlen et al., 2006).The maximum sinking speed is set to 200 m d −1 and is reached at about 5000 m depth over most of the ocean since z max is generally less than 100 m.We have not included any ballasting effect due to the higher density of biogenic silica or calcite (Klaas and Archer, 2002;Armstrong et al., 2002).In fact, observations are rather contradictory on this ballast effect (Lee et al., 2009).In particular, the greater efficiency of the vertical sedimentation of organic matter when associated with calcite and biogenic silica may be due rather to the protection of an organic matter fraction by the inorganic matrix (Moriceau et al., 2009;Engel et al., 2009 Here we present a brief overview of the model of Kriest and Evans (1999).The reader is referred to the literature, where the method has been presented (e.g., Kriest andEvans, 1999, 2000;Kriest, 2002), for more detail.The model postulates that the carbon content (m(d i )), the sinking speed (w(d i )) and the abundance of the aggregates (n(d i )) can be described by power-law functions of their diameters (d i ): It is also assumed, as in Kriest (2002), that aggregates above a certain size L have a constant sinking speed w L .
The slope of the size spectrum can be computed from the total number of aggregates (NUM) and the total mass of particles (POC), which are the two state variables of the model: where m l is the mass of the smallest aggregate (of size l).
Having , the average sinking speed of numbers (w NUM ) and mass (w POC ) can be computed following Kriest (2002): The number of particles and the mass of particles change independently.For instance, sinking tends to remove larger particles.As a consequence, the relationship between the number of particles and their mass evolves with time and space and so does .As a result, the sinking speeds for both mass and number vary with space and time.
Aggregation (ξ ) depends on the particle abundance, their size distribution, rate of turbulent shear and the difference in particle sinking speeds as well as the stickiness (the probability that two particles stick together after contact).The approach implemented in PISCES follows that described in Kriest (2002); see Kriest (2002) for the term ξ and its computation.Currently it is assumed that turbulent shear rate is high in the mixed layer (1 m s −1 ), and low below (0.01 m s −1 ).Summing up the number of collisions due to turbulent shear and differential settlement, C sh and C ds , respectively, the decrease of the number of particles due to aggregation is then: In PISCES, the stickiness (the efficiency of the collisions) is set to a constant value in the namelist.
The temporal evolution of the mass of particles is given as This is exactly equal to the sum of the two equations used for the temporal evolution of POC and GOC in the twocompartment model of PISCES (see Eqs. 37 and 40).
In this equation, each process affecting the mass of the particles is divided by the mean mass ( m) of the compartment exerting this process to convert to numbers.

Iron in particles
In this subsection, the description corresponds to the twocompartment version of the model.To obtain the Kriest version, the equations for both SFe and BFe, the iron content of the small and big particles, respectively, should be simply summed.
θ Fe,P 0.5R CaCO 3 (m P P P + K m P + sh × w P P 2 ) where Fe is the free form of dissolved iron.Its determination is detailed in Sect.4.5.3.Bactfe is the amount of iron taken up by bacteria which is lost as particulate organic iron.Its computation is detailed in Sect.4.5.3.The free form of dissolved iron Fe is the only form of iron that is assumed to be susceptible to scavenging.The scavenging rate of iron is made dependent upon the particulate load of the seawater as follows (e.g., Honeyman et al., 1988;Parekh et al., 2004): Implicitly, in this equation, it is assumed that the affinity of iron for the different types of biogenic particles is the same.Iron is also scavenged by lithogenic particles originating from dust deposition as evidenced by mesocosm experiments (Wagener et al., 2010).The concentration of lithogenic particles is estimated as described in Eq. ( 84).Model estimates (Ye et al., 2011) suggest a different affinity for these particles compared to biogenic particles, which justifies the split between biogenic and lithogenic materials in Eq. ( 50).The amount of iron that is scavenged by POC (λ Fe POC Fe ) and GOC (λ Fe GOC Fe ) is then allocated to SFe and BFe, respectively.

PSi
The dissolution rate of PSi depends on in situ temperature and on silicic acid saturation following the parameterization proposed by Ridgwell et al. (2002): Si eq = 10 6.44− 968 The evolution of λ PSi as a function of Si and of temperature is shown on Fig. 4.
Laboratory experiments show that the diatom frustule is made of two biogenic silica phases which dissolve simultaneously, but at different rates (e.g., Kamatani et al., 1980;Van Capellen et al., 2002;Moriceau et al., 2009;Loucaides et al., 2012).The first phase dissolves significantly faster than the second phase.It is associated with membrane lipids and amino acids and represents about one-third of the frustule (Moriceau et al., 2009).However, the existence of these two phases is still a matter of debate as it has been hypothesized to be a result of the experimental design of the dissolution experiments (Loucaides et al., 2012).In PISCES, despite this uncertainty, we model silica dissolution using two phases.The proportion of the most "labile" phase is set to a constant (χ lab ) in the upper ocean and is computed in the rest of the ocean assuming steady state: Nitrification (Nitrif) corresponds to the conversion of ammonium to nitrate due to bacterial activity.It is assumed to be photo-inhibited (e.g., Horrigan et al., 1981;Yoshioka and Saijo, 1984) and reduced in suboxic waters: where PAR is the PAR averaged over the mixed layer and (O 2 ) varies between 0 (oxic conditions, O 2 > O 2 min,1 ) and 1 (anoxia) according to When waters become suboxic, nitrate instead of oxygen is consumed during the remineralization of organic matter, i.e., denitrification (Denit).The N / C stoichiometric ratio of denitrification R NO 3 can be computed from R −O 2 /NO 3 and is found to be 0.86 (Paulmier et al., 2009).Equation (57), implies that denitrification stops at oxygen concentration above 6 µM (Lipschultz et al., 1990).We further assume complete oxidation by nitrate of the ammonia released from organic matter during denitrification.This oxidation rate has been arbitrarily set to the same value as nitrification rate (λ NH 4 ).
Finally, nitrogen fixation is parameterized in PISCES as follows: This very crude parameterization is based on the following assumptions that have been inferred from studies of Trichodesmium (e.g., Mills et al., 2004;Masotti et al., 2007;Zehr, 2011): -Nitrogen fixation is restricted to warm waters above 20 • C (µ P > µ P (20) = 2.15).
-Nitrogen fixation is restricted to areas with insufficient nitrogen (L P N < 0.8).-Nitrogen fixation requires iron and phosphorus.
-Nitrogen fixation needs high light levels, i.e., E fix is high.
The scaling factor N m fix is set from the namelist and thus, may be chosen by the user.

Phosphate
All terms in this equation have been described previously.

Iron
Iron scavenging (Scav) has been described previously in Sect.4.4.3.Iron is present in seawater largely as colloids (e.g., Wu et al., 2001;Wu and Boyle, 2002;Boyd and Ellwood, 2010).These colloids may aggregate with dissolved organic matter as it forms gels.Thus, they may be transferred to the particulate pool, and settle to the ocean floor.Very few models have incorporated this potential important sink of dissolved iron (Ye et al., 2009(Ye et al., , 2011)).In PISCES, we model this process following the approach chosen for DOM (see Sect. 4.3): where Fe coll is computed from the iron chemistry model (see below).
When dissolved iron concentration exceeds the total ligand concentration L T , scavenging is enhanced as it is done in many other biogeochemical models (e.g., Moore et al., 2004;Dutkiewicz et al., 2005): This scavenging loss term is assumed to be definitive, i.e., iron is permanently removed from the ocean by this process.
Heterotrophic bacteria acquire iron from seawater using siderophore-based iron transport systems (Haygood et al., 1993;Martinez et al., 2000).Observations show that they have quite elevated Fe / C ratios and account for a significant fraction of the total biological uptake of iron (Tortell et al., 1996(Tortell et al., , 1999)).The bacterial uptake of iron is parameterized according to where θ Fe,Bact max denotes the maximum Fe / C ratio of bacteria.The different iron pools are computed using a chemistry model.Two different chemistry models are available in PISCES: -A simple chemistry model based on one ligand (L) and two dissolved iron forms: dissolved inorganic iron (Fe ) and dissolved complexed iron (FeL).
-The complex chemistry model of Tagliabue and Arrigo (2006) as modified by Tagliabue and Völker (2011).This model is based on two ligands (L W and L S ) and five iron forms: free Fe(II) (Fe(II) ) and Fe(III) (Fe(III) ), Fe(III) bound to the weak ligand (FeL W ), Fe(III) bound to the strong ligand (FeL S ) and solid iron (Fe P ).
The complex iron model is activated in PISCES setting the Boolean variable ln_fechem to true.Our main purpose is not to provide a fully detailed description of both chemistry models as they have been described fairly extensively elsewhere.For the simple chemistry model, the reader should refer to Aumont and Bopp (2006), whereas the complex model is detailed in Tagliabue and Völker (2011).For the complex model, all chemical constants have identical values to what was chosen in Tagliabue and Völker (2011) and are thus not listed in Table 1a-e.Only a very brief description of both models will be given here, especially for the complex model.Both models are based on the assumption that chemical reactions are fast enough relative to the other biogeochemical processes affecting iron (for instance phytoplankton uptake) that they can be considered at equilibrium.

Simple chemistry model
Dissolved iron is assumed to be in the form of free inorganic iron Fe and of "complexed" iron FeL.Both forms of iron are assumed to be equally susceptible to consumption by phytoplankton despite recent observations suggest that this may be not the case (Nishioka and Takeda, 2000;Chen and Wang, 2001;Chen et al., 2003).In other words, the total bioavailable concentration of iron is equal to the total dissolved iron concentration (Fe).The chemical speciation of iron is deduced from the three following equations The chemical equilibrium constant K Fe eq is computed from the formulation proposed by Liu and Millero (2002).Solving this set of equations is equivalent to solve a second-order polynomial equation in Fe , whose solution is Colloidal iron is assumed to represent 50 % of FeL: The total ligand concentration L T can be either constant over the ocean, using a value defined in the namelist or can be variable using the relationship proposed by Tagliabue and Völker (2011): where L T is in nmol L −1 and DOC in µmol L −1 .

Complex chemistry model
The iron chemical system is governed by the following set of four equations A supplementary reaction has been added relative to the original set of equations.In the Pacific Ocean, thermal (dark) reduction of Fe(III) organic complexes has been shown to produce the accumulation of a sizeable amount of Fe(II) in the mesopelagic zone (Hansard et al., 2009).Additional constraints are given by the conservation of total dissolved iron (Fe), L WT and L ST over the fast timescale: Solving this system of equations is equivalent to solving a third-order polynomial equation in Fe(III) (Eq.16 in Tagliabue and Völker, 2011).Because thermal aphotic reduction of FeL W has been added here, the definition of some coefficients in the original study has changed: where k th has been set to 0.0048 h −1 .Then, knowing Fe(III) , the other four iron species can be computed.
Observations suggest that the weak ligand (L W ) is ubiquitous in the water column and is probably produced by the degradation of organic matter sinking from the upper layers of the ocean.The strong ligand is present in the upper ocean and is most probably produced by autotrophic and heterotrophic bacteria (for instance siderophores) (e.g., Boyd and Ellwood, 2010).In PISCES, we assume that two-thirds of the total ligand concentration above 0.6 nmol L −1 is going to L S , the rest is attributed to L W : As in the simple chemistry model, the ligand concentration L T can be either constant over the ocean, using a value defined in the namelist or can be variable using the relationship proposed by Tagliabue and Völker (2011) (see Eq. 67).
The rate constants required by the model are identical to those described by Tagliabue et al. (2009a) as modified by Tagliabue and Völker (2011).Furthermore, we have slightly changed the formulation of the oxidation rate constant used in the original model: This avoids numerical problems in strongly anoxic areas where oxygen concentration is close to 0. Bioavailable iron can be defined either as Fe(II) + Fe(III) + FeL S or as Fe(II) + Fe(III) + FeL S + FeL W . k th has assigned the value computed from the observations by Hansard et al. (2009), consistent with the data of Pullin and Cabaniss (2003).Colloidal iron and dissolved inorganic iron are defined as Fe coll = 0.5(Fe p + FeL W + FeL S ), (73a) Fe = Fe(III) + Fe(II) . (73b) We assumed that 50 % of the iron bound to ligands and of the particulate inorganic iron is colloidal iron.

Si
All terms in this equation have been already defined previously.

Calcite
In PISCES, calcium carbonate is assumed to exist only in the form of calcite.Thus, aragonite is not considered, for instance, for the computation of chemical dissolution in the water column.
The biological production of sinking calcite is defined as The rain ratio R CaCO 3 is variable.We propose the following parameterization for this ratio: This parameterization is based on a set of very simple assumptions, mainly inferred from the review by Zondervan (2007): -Coccolithophores are not very abundant in very oligotrophic waters.
-Calcification tends to be maximum at intermediate light levels and decrease at either high and low light levels, around 30 and 4 W m −2 , respectively.
-Coccolithophores are not found when the temperature of sea water is below 0 • C.
-Coccolithophores are found in stratified waters.Their abundance decreases when the mixed layer depth (z mxl ) exceeds 50 m.
-Maximum levels of coccolithophores are found in the mid-latitudes, where temperature is around 10 • C.
We recognize that this parameterization is quite ad hoc and may seem arbitrary.But as it will be shown, it simulates reasonable calcification patterns and alkalinity distribution (yet we recognize that it could be for the wrong reasons).Furthermore, it avoids an explicit modeling of the coccolithophores which is far from being trivial.Only part (η I ) of the grazed shells are routed to sinking calcite.The rest is taken to dissolve in the acidic guts of zooplankton (Jansen and Wolf-Gladrow, 2001).This dissolution is still debated.However, observations tend to show that a significant proportion of the sinking shells is lost in the upper ocean, with this being associated with grazing as well as other mechanisms (Milliman et al., 1999).
The dissolution of calcite is modeled as in Gehlen et al. (2007): 4.7 The carbonate system All terms in the above equations have been described previously in this document.In addition to these biogeochemical fluxes, the ocean exchanges CO 2 with the atmosphere at the sea surface.The gas exchange coefficient is computed from the relationship proposed by Wanninkhof (1992).No exchange is allowed with the atmosphere across sea ice: where % ice is the concentration of sea ice which varies between 0 and 1.The carbonate chemistry follows the OCMIP protocols (more information at http://ocmip5.ipsl.jussieu.fr/OCMIP/) except that it has been simplified to reduce the computing cost: alkalinity only includes carbonate, borate and water (H + , OH − ).Atmospheric pCO 2 can be set as an external tunable parameter via a namelist parameter or read from a file.Its value is uniform over the global ocean (no spatial gradient) and is not allowed to vary in response to the air-sea fluxes.This means that PISCES does not include an interactive atmospheric (box or more complex) model (although this functionality can be added very easily).Finally, the impact of atmospheric pressure on pCO 2 can be accounted for by setting the Boolean ln_presatm to true in the namelist.In that case, the 2-D spatial distribution of atmospheric pressure should be read in a file.

Oxygen
In this equation, the stoichiometric ratio O ut 2 represents the change in oxygen relative to carbon when ammonium is converted to organic matter, whereas O nit 2 denotes the consumption of oxygen during nitrification.Their values have been set respectively to 131/122 and 32/122 so that the typical Redfield ratio for oxygen is equal to 1.34 as proposed by Kortzinger et al. (2001).
Oxygen is exchanged with the atmosphere using the parameterization of Wanninkhof (1992) to compute the gas exchange coefficient.The atmospheric concentration of oxygen is constant over time and space and cannot be specified by the user.As for CO 2 , no air-sea fluxes are allowed when the ocean is covered by sea ice (see Eq. 82).

External supply of nutrients
Nutrients are supplied to the ocean from five different sources: atmospheric dust deposition, rivers, sea ice, sediment mobilization and hydrothermal vents.

Atmospheric deposition
The model can include the atmospheric supply of Fe, Si, P and N. The former three sources (Fe, Si and P) are dependent on each other as they are computed from the same dust input file.They are activated in PISCES by setting the Boolean ln_dust to true.Otherwise, no atmospheric source of Fe, P and Si is prescribed.Furthermore, in that case, the dust concentration in the ocean (used for instance in Eq. 50) is set to 0. The iron content of dust is set to a constant value specified in the namelist.Its default value is 3.5 % which is the average content of crustal material (e.g., Taylor and McLennan, 1985;Jickells and Spokes, 2001;Jickells et al., 2005).The solubility of dust iron in sea water can be either set to a constant value in the namelist or can be read from a file if ln_solub is set to true.Once it has left the surface layer, particulate inorganic iron from dust is still assumed to experience dissolution.The dissolution rate is computed assuming that mineral particles sink at a constant speed specified in the namelist and that about 0.01 % of the particulate iron dissolves in a day (Bonnet and Guieu, 2004).This is equivalent to a remineralization length scale of 20 000 m if the sinking speed is set to a typical value of 2 m day −1 , on the same order as the length scale prescribed for the same process by Moore et al. (2004).Atmospheric deposition of Si is also considered following Moore et al. (2002b) and is restricted to the first layer of the model.Atmospheric deposition of P is computed from dust deposition assuming that the total phosphorus content of dust is 750 ppm (Mahowald et al., 2008) and that the solubility in surface sea water is 10 % (Ridame and Guieu, 2002; Mahowald et al., 2008).As for Si, deposition is restricted to the first level of the ocean model.Atmospheric deposition of N is treated separately from the deposition of the other nutrients and can be activated in the model by the Boolean ln_ndepo.All nitrogen deposited at the ocean surface is assumed to dissolve.We made the quite strong assumption for all nutrients that sea ice does not alter the deposition fluxes.
The dust (Dust) concentration in the ocean is modeled in a very simplistic way in PISCES.It is computed from dust deposition assuming a constant sinking speed (the same as the sinking speed used to compute iron dissolution from dust in the interior of the ocean).Furthermore, dust is not transported by the ocean currents.This assumption is made in PISCES to avoid adding another prognostic tracer in the model.As a consequence, the concentration of dust is computed as where D dust is dust deposition at the surface and w dust is the prescribed sinking speed of dust.

River discharge
River discharge is activated by setting the Boolean variable ln_river to true in the namelist.The river discharge of the different elements is then read from a file that must be provided in that case by the user.The river supply of DIN, DIP, DON, DOP, Si, DIC, alkalinity and DOC need to be provided.As DON, DOC and DOP are not separately modeled in PISCES (fixed stoichiometry), dissolved organic matter is assumed to remineralize instantaneously at the river mouth and thus, DON, DOP and DIC are added to DIN, DIP and DIC, respectively.As a default in PISCES, river supply of all elements but DIC and alkalinity is taken from the GLOBAL-NEWS2 data sets (Mayorga et al., 2010).For DIC and alkalinity, we use results from the Global Erosion Model (GEM) of Ludwig et al. (1996), neglecting the POC delivery as most of it is lost in the estuaries and in the coastal zone (Smith and Hollibaugh, 1993).All fields are interpolated onto the ORCA grid and co-localized with the river runoff prescribed in the physical model.Iron is also delivered to the ocean by rivers.The amount of supplied iron is computed from the river supply of inorganic carbon, assuming a constant Fe / DIC ratio.This ratio is determined so that the total Fe supply equals 1.45 Tg Fe yr −1 as estimated by Chester (1990).

Reductive mobilization of iron from marine sediments
Reductive mobilization of iron from marine sediments have been recognized as a significant source to the ocean (Johnson et al., 1999;de Baar and de Jong, 2001;Moore et al., 2004).Fe concentrations in the sediment pore waters are often several orders of magnitude larger than in the seawater.
A large part of the iron released to the ocean either by diffusion or by resuspension is likely to be oxidized in insoluble forms and trapped back to the sediments, at least in oxygenated waters (de Baar and de Jong, 2001).Yet, some of this iron should escape as observations clearly show increasing concentration gradients of particulate and dissolved iron toward the coastal zones.Unfortunately, almost no quantitative information is available to parameterize this potentially important source.Observations from benthic chambers indicate that this source may be controlled by the oxygen concentrations overlying the sediments (Raiswell and Anderson, 2005;Severmann et al., 2010) and perhaps the magnitude of the organic carbon export to the sediments (Elrod et al., 2004).Such potential relationships are not yet embedded in PISCES.
In a way similar to Moore et al. (2004), we apply a maximum constant iron source from the sediments.Since anoxic sediments are more likely to release iron to the seawater, we have modulated this source by a factor (Fsed) computed from the metamodel of Middelburg et al. (1996): ζ Fsed = −0.9543+ 0.7662 ln(z Fsed ) − 0.235(ln(z Fsed )) 2 , (85b) Fsed = min (1, exp(ζ Fsed )/0.5) . (85c) From this metamodel, it is possible to estimate the relative contribution of anaerobic processes to the total mineralization of organic matter in the sediments, and thus to have an indication on how well the sediment is oxygenated (Soetaert et al., 2000).Our modulation factor is sim-ply set equal to this relative contribution.The maximum iron flux from the sediments has been set by default to 2 µmol Fe m −2 d −1 by adjusting the modeled iron distribution to the few iron observations available over the continental margins.This value is identical to that used by Moore et al. (2004) in their model.The maximum iron flux constant can be specified in the namelist and thus, may be changed from the default value by the user.
Unfortunately, as a consequence of the relatively coarse resolution of ORCA2, the model bathymetry is not able to correctly represent the critical spatial scales of the ocean bathymetry.An example is the continental shelves, which typically have a width scale of 10-30 km, which can be approximately an order of magnitude less than the horizontal resolution of the model.In order to take sub-model grid scale bathymetric variations into account in the Fe source function, the model grid structure has been compared with the highresolution ETOPO5 data set.An algorithm was developed whereby for each and every horizontal grid cell, the corresponding region in the ETOPO5 data set is considered.For each vertical level in the model corresponding to a particular horizontal grid point, the corresponding ocean-bottom area from ETOPO5 (in fractional units) is saved, with the end result being a three-dimensional array containing an equivalent area for the bottom bathymetry of the ocean for the ETOPO5 data set.The iron flux computed as described above is then multiplied by this fractional area % sed (which varies between 0 and 1): This corresponds to a global flux of 34 Gmol Fe yr −1 .

Iron from hydrothermalism
Recent studies have shown that hydrothermalism may deliver to the deep ocean a significant amount of dissolved iron (e.g., Mackey et al., 2002;Boyle and Jenkins, 2008;Bennett et al., 2008;Toner et al., 2009).Despite very large uncertainties, this source has been estimated, based on discrete data and a model, to 3 to 9 × 10 8 mol Fe yr −1 globally (Bennett et al., 2008;Tagliabue et al., 2010).In PISCES, this source is included following the modeling study by Tagliabue et al. (2010) and may be activated by setting the Boolean ln_hydrofe to true.The hydrothermal flux of iron has been computed based on observed correlations between 3 He and dFe (Boyle et al., 2005;Boyle and Jenkins, 2008) and using a data compilation of dFe / 3 He (see the Supplement of Tagliabue et al., 2010).Then, the spatial distribution of this flux has been derived from previous modeling works on 3 He, which relate the 3 He flux to the ridge-spreading rates (Farley et al., 1995;Dutay et al., 2004); 0.2 % of the delivered iron is assumed to be soluble.

Iron from sea ice
The last external source of nutrients which is taken into account in PISCES is the exchange of iron between the ocean and the sea ice associated with formation and melting.This source is activated by setting the Boolean variable ln_ironice to true.The receding ice edge is often characterized by intense phytoplankton abundance which can be explained by ocean stratification promoted by the melting of sea ice (Smith and Nelson, 1985) as well as the releases of iron accumulated in sea ice during winter (Sedwick and Di Tullio, 1997;Tagliabue and Arrigo, 2006).Measurements in sea ice have found iron concentrations of more than 1 order of magnitude higher than in adjacent sea water (Lannuzel et al., 2007(Lannuzel et al., , 2008)).About 90 % of this iron has been shown to be of oceanic origin (Lannuzel et al., 2007).Thus, iron is taken up from sea water when ice forms and is released back to the ocean when it melts.Lancelot et al. (2009) have studied the impact of this source in the Southern Ocean and shown that it is of primary importance in the seasonal ice zone.Their approach relies on the modeling of iron concentration within sea ice.In PISCES, we have simplified this model by assuming that iron concentration in sea ice is constant.In that case, the iron fluxes between the ocean and the sea ice can be computed from the water fluxes between these two reservoirs: where EP oi is the water flux (in kg m −2 s −1 ) from the ice to the ocean and Fe ice is the iron concentration in sea ice which has been found to be on the order of 10 nmol L −1 .In this equation, F ice,− Fe is thus the loss of iron from the ocean when sea ice forms and F ice,+ when sea ice melts.It should be noted here that since we do not model iron in sea ice, the exchange of iron between both reservoirs is not conservative.In the model configuration presented here, ice represents a net source of iron of 0.024 Gmol Fe yr −1 .

Bottom boundary conditions
At the bottom of the ocean, the exchange between the sediments and the ocean can be represented either with or without a sediment model.The sediment model is activated by using the cpp key key_sed.This model will not be described in this document.It is basically identical to the model of Heinze et al. (1999) with some modifications as described by Gehlen et al. (2006).The main modification is the addition of denitrification to the set of early diagenetic reactions.Parameter values are identical to those in Heinze et al. (1999).
When the sediment model is not activated, very basic but different treatments are applied at the bottom of the ocean depending on the tracer considered.For biogenic silica, the amount of particulate material that is permanently buried in the sediments is assumed to exactly balance the external input from dust deposition and river discharge, described in the previous section.Then, we assume that the part of biogenic silica that is not permanently buried redissolves back to the water column instantaneously.
For particulate organic carbon, we first determine the proportion of organic matter reaching the seafloor that is permanently buried.The burial efficiency is computed using the algorithm proposed by Dunne et al. (2007): where E burial is the burial efficiency and F OC is the flux of organic carbon at the bottom (in mmol C m −2 d −1 ).We then use the metamodel by Middelburg et al. (1996) to determine the proportion of degradation of the remaining organic matter that is due to denitrification: where the tracer concentrations are in µmol L −1 and F OC is the flux of organic carbon at the bottom (in µmol cm −2 d −1 ).
In this equation, oxygen and nitrate concentrations are not allowed to be below 10 µmol L −1 and 1 µmol L −1 , respectively.Then, the fluxes of nitrate and oxygen to the sediment as a consequence of denitrification and oxic degradation, respectively, can be computed: Particulate organic carbon which has been degraded by denitrification and oxic processes is released in the bottom box as ammonium.
A specific treatment of calcite at the sediment interface is embedded in PISCES.The preservation of calcite in the sediments is represented as a function of the saturation level of the overlying waters: where is the calcite saturation level.This relationship has been deduced from the study by Archer (1996).The permanent burial of calcite is modulated by % CaCO 3 .The amount of calcite that is not buried, instantaneously dissolves back to the ocean.

Model parameters and their default values
Table 1a-e list model parameters, their respective units and default values as well as a brief description of each of them.Many of these parameters can be specified in the namelist_pisces file.As much as possible, the parameter values have been derived from the literature.However, many parameters, such as the mortality rates, are either not constrained at all, or only poorly constrained by the observations.Their values have been adjusted by successive simulations evaluated against the observational data sets presented below.
In addition to the parameters above, PISCES includes a number of control parameters defined as Boolean variables that appear in the namelist file namelist_pisces.These variables either allow one to switch between different functional forms or activate additional functionalities.These con-  trol parameters are listed in Table 2. Finally, some functionalities, such as the Kriest model of particulate organic matter, require a major reorganization of the code, for instance a change in the number of prognostic variables.In that case, these functionalities are activated through CPP keys which force the model to be recompiled.These CPP keys are listed in Table 3.

Model results
The objective of this section is not to present a full and exhaustive validation of the model results.This has already been presented in a wide range of publications using different configurations of the model (see the Introduction).
Here we present instead a brief comparison of PISCES with available observations, in its standard global configuration.This configuration is the default setup available when downloading the code from the NEMO web site (the standard ORCA2_OFF_PISCES configuration).All the necessary input files can be obtained from this web site.

Model setup
The dynamical state of the ocean has been simulated using the ocean physical model ORCA2-LIM in version 3.2 (Madec, 2008).This model is based on an ocean general circulation model OPA9, coupled with the sea ice model Louvain-la-Neuve Ice Model (LIM2) (Timmermann et al., 2005).The spatial resolution is about 2 • by 2 • cos (where is the latitude) with a focusing of the meridional resolution to 0.5 • in the equatorial domain.The model has 30 vertical layers, with an increased vertical thickness from 10 m at the surface to 500 m at 5000 m.Representation of the topography is based on the partial step thicknesses (Barnier et al., 2006;Penduff et al., 2007).Lateral mixing along isopycnal surfaces is performed both on tracers and momentum as in Lengaigne et al. (2003).The parameterization of Gent and McWilliams (1990) is applied poleward of 10 • to represent the effects of non-resolved mesoscale eddies.Vertical mixing is parameterized using the turbulent kinetic energy (TKE) scheme of Gaspar et al. (1990), as modified by Madec (2008).
The fields used to drive the ocean are identical to those used by Aumont and Bopp (2006).However, the resulting physical circulation state simulated by the ocean model is different as several new parameterizations and new algorithms have been included in ORCA2-LIM.Climatological atmospheric forcing fields have been constructed from various data sets consisting of daily NCEP/NCAR 2 m atmospheric temperature averaged over 1948-2003(Kalnay et al., 1996)), monthly relative humidity (Trenberth et al., 1989), monthly ISCCP total cloudiness averaged over 1983-2001(Rossow and Schiffer, 1999), monthly precipitation averaged over 1979-2001(Xin and Arkin, 1997) and weekly wind stress based on European Remote-Sensing Satellite (ERS) satellite product and TAO observations (Menkes et al., 1998).Surface heat fluxes and evaporation are computed using empirical bulk formulas as described by Goose (1997).To avoid any strong model drift, modeled sea surface salinity is restored to the monthly WOA01 data set (Conkright et al., 2002) with a nudging timescale of 40 days applied through local freshwater forcing (thereby conserving salt).The ocean dynamical model has been spun-up for 200 years, starting from rest and from the climatology of Conkright et al. (2002) for temperature and salinity.
Phosphate, oxygen, nitrate and silicic acid distributions have been initialized at uniform concentrations inferred from  observed climatologies (Garcia et al., 2010).Initial values for dissolved inorganic carbon and alkalinity are taken from the OCMIP guidelines (Orr, 1999).The ecological tracers are initialized uniformly to arbitrary low values.Iron concentrations are set everywhere to 0.6 nM.The model is then spun-up offline for 4000 years using the circulation state predicted by the dynamical model.Atmospheric pCO 2 is set to a pre-industrial value of 278 ppm.After this integration, primary productivity as well as CO 2 fluxes drift by less than 0.001 Gt C yr −1 .As the external sources and sinks of nutrients are not fully balanced (see the model description), the global inventories of phosphate, nitrate, alkalinity and silicate are restored toward the observed inventories, once a year on 1 January.In practice, this correction is done by scaling the 3-D concentrations with a constant uniform factor so that the simulated total inventories do not drift away from the observed inventories.Thus, we do not restore the simulated 3-D distributions to 3-D observed fields so that the predicted spatial and temporal patterns are not corrected in any way to better match the observations.However, the predicted global inventories of P, N, Si and alkalinity can not be used to evaluate the model skill since they are not prognostically predicted.Anyhow, this correction is very small and corresponds to a relative change in the concentration of the tracers on the order of 1-5 × 10 −5 yr −1 ; therefore, that no significant jump is introduced by this technique.The activation of this technique as well as the frequency at which it is applied are controlled by a Boolean parameter and a parameter respectively, in the namelist file namelist_pisces (see Table 2).

Global budget
Table 4 presents the global carbon budget as simulated by PISCES, when embedded in ORCA2-LIM.The annual net predicted primary production is 44 Gt C yr −1 .This value falls on the lower bound of the broad estimates given by satellite observations which give values between 37 and 67 Gt C yr −1 (Longhurst et al., 1995;Antoine et al., 1996;Behrenfeld and Falkowski, 1997;Behrenfeld et al., 2005).Using PISCES in a higher resolution model would certainly produce a significantly larger number as mesoscale and submesoscale processes have been shown to stimulate biological productivity (McGillicuddy et al., 1998;Oschlies and Garçon, 1998;Lévy et al., 2001), and coastal regions, characterized by a intense primary productivity, are not properly resolved by the coarse grid.About 17 % of the primary production is due to diatoms.Global estimates of the contribution of diatoms to total production are rather uncertain and broad.Nelson et al. (1995) have suggested that diatoms may be responsible for up to 40 % of the total primary production.However, as discussed by Aumont and Bopp (2006), this value is certainly overestimated.In recent years, algorithms, which attempt to retrieve the composition of phytoplankton from space, have been developed (e.g., Alvain et al., 2005;Uitz et al., 2006;Hirata et al., 2008;Brewin et al., 2010).Only a few of these methods give quantitative estimates of the contribution of the different species or size classes to total biomass or primary productivity (Brewin et al., 2011).The estimated global contribution of diatoms from these methods ranges from as low as 7 % to as high as 32 % of the total phytoplankton (Uitz et al., 2010;Hirata et al., 2011) (if one assumes crudely that microphytoplankton are effectively equivalent to diatoms).Finally, ocean biogeochemical models predict the contribution of diatoms to be between 15 and 30 % (e.g., Moore et al., 2002a;Aumont et al., 2003;Dutkiewicz et al., 2005;Yool et al., 2011).
Export production at 150 m is estimated to be 6.9 Gt C yr −1 ; 86 % of this export is related to settling particles (one-third by the small sinking particles and twothird by the fast sinking particles).The remainder is due to vertical advection and diffusion of dissolved organic carbon, which occurs mainly in the mid-ocean gyres (vertical advection) and in the high latitude regions during winter (vertical diffusion).Constraining export production is rather difficult, if not impossible, considering the very broad range given by estimates either based on models or observations and the different definitions of export production, in particular the depth horizon at which it is estimated (e.g., Eppley and Peterson, 1979;Schlitzer, 2000;Moore et al., 2002a;Yool et al., 2011).Mesozooplankton grazes about 9 % of total primary production.This value is close to other estimates either based on observations (Calbet, 2001) or models (Moore et al., 2002a;Buitenhuis and Geider, 2010).Total gazing by mesozooplankton is predicted to be 11.2Gt C yr −1   2010) when extrapolating observations.Despite estimates of grazing by microzooplankton are quite badly constrained, this might suggest that it is overestimated in the model.
Table 5 shows the calcite and silicon budgets for the upper 150 m of the ocean.Production of calcite and export at 150 m are simulated to be, respectively, about 1.6 and 0.8 Gt C yr −1 by PISCES.These numbers fall within the limits of the quite large range of 0.4-1.8Gt C yr −1 estimated either for global calcification or export of particulate inorganic carbon (PIC) (Murnane et al., 1999;Lee, 2001;Moore et al., 2002a;Balch et al., 2007;Berelson et al., 2007).For silicate, the model predicts a vertical export of biogenic silicate of 106 Tmol Si yr −1 .This value is within the 105 ± 17 Tmol Si yr −1 estimated for the global ocean (Tréguer and De La Rocha, 2012).Global production of biogenic silica by diatoms is 146 Tmol Si yr −1 in our model.This value is quite low compared to the 239 Tmol Si yr −1 given by Tréguer and De La Rocha (2012).About 27 % of biogenic silica dissolves in the top 150 m of the ocean, half the estimate of Nelson et al. (1995) and Tréguer and De La Rocha (2012).However, as already mentioned, because of its coarse resolution, the physical model configuration does not properly resolve the coastal zones.For the open ocean only (in a strict sense), Tréguer and De La Rocha (2012) estimated biogenic silica production to be about 103 Tmol Si yr −1 .Not surprisingly then, considering the limitations due to the spatial resolution, our modeled estimate is between the open ocean and global values.The mean Si / C for uptake of diatoms as predicted by PISCES is thus 0.23, which is high relative to the optimal Si / C of 0.13 (Brzezinski, 1985).This suggests thus that over most of the ocean, diatom cells are stressed, not a very surprising result.Furthermore, a large part of the biogenic silica production oc- curs within the Southern Ocean, a region where diatom cells are very heavily silicified (Baines et al., 2010).Table 6 presents the global nitrogen budget as simulated by PISCES.River discharge and atmospheric deposition of nitrogen are given by the prescribed input fields to PISCES.By definition, burial in the sediments is set exactly equal to river discharge.Nitrogen fixation is predicted to be 111.8Tg N yr −1 .This value is close to the mean value of about 140 Tg N yr −1 estimated from direct observations or nutrients analysis (Capone et al., 1997;Deutsch et al., 2007).Figure 6 shows a comparison between the spatial distribution of observed nitrogen fixation rates from the MA-Rine Ecosystem DATa (MAREDAT) project and that as simulated by PISCES.This indicates that, despite a quite simplistic formulation, the model is able to capture the main observed patterns, at least on an annual-mean basis.Modeled denitrification in the water column and in the sediments are about 78 and 93 Tg N yr −1 , respectively.Sediment denitrification estimates are significantly higher, in the range of 130-300 Tg N yr −1 (Codispoti et al., 2001;Galloway et al., 2004;Gruber, 2004).However, considering the coarse spatial resolution of the model, this is expected as most of benthic denitrification occurs over the continental margins.The sources and sinks of nitrogen are slightly unbalanced, with the sources exceeding the sinks by about 21 Tg N yr −1 .

Chlorophyll
The modeled chlorophyll distribution is compared to GLOB-COLOUR satellite observations for two seasons in Fig. 7.The seasons have been defined to roughly correspond to bloom periods in the high latitudes.The observed patterns are qualitatively reproduced by the model.Slightly too low chlorophyll concentrations are simulated in the subtropical gyres.This discrepancy may be explained by the lack of acclimation dynamics to oligotrophic conditions in the model or by the assumption of constant stoichiometry ei- ther in phytoplankton or in organic matter (Ayata et al., 2013).Chlorophyll concentrations are quite strongly underestimated in the equatorial Atlantic and in the Arabian Sea.In the latter region, mesoscale and submesoscale processes have been shown to be of critical importance (Lee et al., 2000;Kawamiya, 2001;Hood et al., 2003).A model study, using PISCES coupled to a higher resolution version of NEMO, has been shown to simulate chlorophyll distribution in much better agreement with the observations (Koné et al., 2009).Chlorophyll concentrations are high in the eastern boundary upwelling systems.The sedimentary source of iron plays a role in these systems.When this iron source is not included in models, modeled chlorophyll concentrations are much lower (Aumont and Bopp, 2006;Moore and Braucher, 2008).
In two of three main HNLC regions, i.e., the equatorial Pacific and the eastern subarctic Pacific, the model succeeds in reproducing the moderate chlorophyll concentrations.In spring, chlorophyll levels are strongly overestimated east of Japan.As in all coarse resolution models, the ocean circulation in this region is not correctly represented with an incorrect trajectory of the Kuroshio current (i.e., Gnanadesikan et al., 2002;Dutkiewicz et al., 2005;Aumont and Bopp, 2006).Simulated mixed layer depths are too deep in winter and as a consequence the spring bloom is very strong (similar features occur in the North Atlantic).In the equatorial Pacific Ocean, a minimum threshold value has been imposed on iron (0.01 nmol L −1 ) in the model.If not used, chlorophyll concentrations become much too low on both sides of the Equator, resulting in an accumulation of macronutrients and a poleward migration of the southern (northern) boundary of the northern (southern) subtropical gyre (see Fig. 5 in Tagliabue et al., 2009a).The existence of such a threshold suggests that either a minor but regionally important source of iron is missing in PISCES (for instance the dissolution of particulate inorganic iron) or that the standard iron chemistry is too simple (Tagliabue et al., 2009a;Tagliabue and Völker, 2011).
In the Southern Ocean, the third and largest of the principal HNLC regions, chlorophyll concentrations appear to be strongly overestimated by the model when evaluated against satellite-derived observational products, especially during summer.Furthermore, the increase in phytoplankton in late spring and early summer occurs too early.However, numerous studies comparing satellite chlorophyll to in situ data have shown that the standard algorithms used to deduce chlorophyll concentrations from reflectance tend to underestimate in situ observed values by a factor of about 2 to 2.5, especially for intermediate concentrations (e.g., Dierssen and Smith, 2000;Korb et al., 2004;Garcia et al., 2005;Kahru and Mitchell, 2010).Clearly, evaluating the model in the Southern Ocean is quite challenging and requires a more thorough systematic analysis of both the model and the available data sets.

Iron
Figure 8 shows the distribution of iron at three different depth ranges for the model and for the observations.The observational distributions come from the recently published database of Tagliabue et al. (2012) augmented with about 1000 recent observations.The data set can be downloaded from http://pcwww.liv.ac.uk/~atagliab.A complete and exhaustive validation of the model is made difficult by the relative sparsity of the data.
As expected, the highest concentrations of iron in the open ocean are found in the subtropical North Atlantic Ocean and in the Arabian Sea.Those high values are produced by the enhanced dust deposition, mainly emanating from the Sahara desert.The model tends to underestimate the maximum values found in both basins.Interestingly, the local minimum, which is observed west off Mauritania just below the maximum Saharan dust plume, is well captured by the model.Such a minimum is explained by the combination of very low solubilities of the iron contained in the Saharan dust particles when they are close to their source region (Bonnet and Guieu, 2004;Luo et al., 2005) with enhanced scavenging by the dust particles deposited at the ocean surface (Wagener et al., 2010).Very high iron concentrations, typically above 1 nmol L −1 are both observed and modeled along the coasts and over the continental margins as a result of sediment mobilization.As already mentioned in the previous section, this strong source of iron sustains the high productivity observed along the coasts (Johnson et al., 1999), in the eastern boundary upwelling systems (Bruland et al., 2005) but also downstream of the islands, especially in the Southern Ocean (Blain et al., 2007;Pollard et al., 2007;Korb et al., 2008).In the rest of the open ocean, iron concentrations are typically low, generally below 0.2 nmol L −1 , especially in the HNLC regions.PISCES tends to exaggerate these low concentrations.
Iron concentrations increase with depth due to the remineralization of organic particles settling from the surface waters (Johnson et al., 1997;Moore and Braucher, 2008).However, except near the coasts, concentrations rarely exceed 1 nmol L −1 .Again, PISCES captures the main observed patterns both at intermediate depths and in the deep ocean.In the Atlantic Ocean and in the Arabian Sea, iron concentrations remain relatively elevated at intermediate depth in the observations and in the model.In the model, these high values are due to the slow but significant release of iron by the dust particles which sink out from the surface.In the Pacific Ocean, the coastal signature extends far beyond the coastal domain.For instance, it has been proposed as a potential explanation for the episodic blooms observed at station P in the northeastern subarctic Pacific Ocean (Lam et al., 2006;Misumi et al., 2011).In the deepest waters of the Pacific and Indian oceans, iron concentra- tions tend to decrease to the bottom of the ocean and they often fall below 0.6 nmol L −1 .Despite the fact that ligands concentrations in seawater are highly variable, they are typically larger than this value which is the uniform ligand concentration chosen in the model experiment shown here (e.g., Wu and Luther, 1995;Boyé et al., 2001Boyé et al., , 2003;;Hunter and Boyd, 2007;Ibisanmi et al., 2011).The model explains this decrease by the aggregation of iron colloids which are transferred to the particulate pool and thus sink out of the ocean as hypothesized by several studies (Wu et al., 2001;Ye et al., 2009;Geldhill and Buck, 2012).The lowest iron concentrations in the intermediate and deep ocean are found in the Southern Ocean.Iron concentrations slowly increase with depth to reach about 0.4 nmol L −1 in the deep ocean.Higher values are found along Antarctica due to sediment mobilization.

Nutrients, oxygen, alkalinity and DIC
In this section, the simulated distributions of macronutrients, oxygen, alkalinity and DIC are evaluated against available observations.The observations comprise the World Ocean Atlas 2009 for nutrients and oxygen (Garcia et al., 2010), and the GLobal Ocean Data Analysis Project (GLODAP) database for DIC and alkalinity (Key et al., 2004).
Figures 9 and 10 show the surface distributions of nitrate and silicate and zonally averaged sections in the Atlantic and Pacific oceans.At the surface, the model compares quite well with the observations, especially for nitrate.Nitrate concentrations seem to be slightly overestimated along the Antarctic  (Garcia et al., 2010).Panels are the same as on Fig. 9. coast.However, as most of the data have been collected during the productive season in this region, the climatology is likely to be biased toward low values.The surface silicate distribution is less well represented by PISCES, in particular in the Southern Ocean.The silicate front (defined as the latitude at which silicate becomes exhausted) is located too far north in the model.At depth, both modeled nutrients exhibit the same deficiencies.In the Atlantic Ocean, concentrations in the deep ocean are strongly overestimated.Too shallow North Atlantic deep waters (NADW), with strongly underestimated transport simulated for lower NADW, accounts for this problem (Arsouze et al., 2008;Griffies et al., 2009;Smith et al., 2010).As a result, Antarctic bottom waters, characterized by high silicate and nitrate concentrations, tend to dominate over too large part of the deep Atlantic Ocean.In the Pacific Ocean, both nitrate and silicate concentrations are underestimated in the deep waters of the Northern Hemisphere.
In Fig. 11, the modeled oxygen distribution is evaluated against observations.Not surprisingly, the surface distribution compares quite well to the observations as oxygen is  (Garcia et al., 2010).Panels are the same as on Fig. 9.
close to its solubility value and is thus strongly constrained by sea surface temperature.At depth, the main deficiency is the overestimation of oxygen concentrations in the Pacific Ocean.Ventilation along Antarctica, mainly in the Ross and Weddell seas, is too strong in the physical model.Inspection of the simulated mixed layer depths shows that the mixed layer reaches the bottom of the ocean at several locations along Antarctica (not shown), which is not realistic (de Boyer-Montégut et al., 2004).The nearly homogeneous oxygen concentrations south of 60 • S are a consequence of this too intense winter mixing, which thus ventilates the deep ocean with too much oxygen.Figures 12 and 13 display the modeled and observed distributions of DIC and alkalinity at the surface and along zonally averaged sections in both the Atlantic and the Pacific.Modeled DIC does not include the anthropogenic perturbation since atmospheric CO 2 was set to its pre-industrial value.We have estimated the observed pre-industrial distribution of DIC as the difference between total DIC and anthropogenic carbon, which are both available in GLODAP Key et al. (2004).It should be also mentioned here that no observations were available north of 60 • N. Values north of this latitude have been extrapolated for plotting purpose.At the surface, several modeled features are not visible in the observations.Very low alkalinity and DIC concentrations are predicted in the Bay of Bengal, in the Gulf of Guinea, close to the Indonesian islands and generally at the mouths of the tropical rivers.The lack of observations in these regions may explain this difference, as the GLODAP database is based on a rather coarse sampling coverage.In the deep ocean, the main deficiencies noticed for the macronutrients are apparent in the simulated distributions.

Skill assessment
In this section, we quantitatively estimate the model performance using Taylor diagrams (Taylor, 2001).Taylor diagrams evaluate both the correlation normalized by the observed standard deviation (SD) (circumference axis) and the relative variability (radial axis) of the model and observations.The distance between the model points and the (1,1) coordinate point (defined as the reference point) is equal to the standard root mean error, normalized by the observed SD.The closer the model is to the observations, the closer the points should be to the reference point.Although a number of means and diagnostics exist (Allen et al., 2007;Doney et al., 2009;Vichi and Masina, 2009), Taylor diagrams have become quite popular as they synthesize, in a quite convenient way, several statistical diagnostics.
Figures 14 and 15 show Taylor diagrams for surface chlorophyll and mesozooplankton averaged over the top 150 m of the ocean.The agreement is rather modest for both variables, especially for mesozooplankton.For chlorophyll, the model performs slightly better for annual-mean distributions, which suggests biases in the representation of the seasonal cycle.The Southern Ocean exhibits the poorest agreement.In particular, the model tends to strongly underestimate the spatial variability since the SD is smaller for the annualmean distribution than for seasonally varying fields.In the other basins, the variability is overestimated, especially in the Atlantic Ocean where the spring blooms in the subarctic domain are too intense, at least relative to satellite observations (see Fig. 7).Mesozooplankton variability is strongly underestimated by PISCES in all basins.The use of a square closure scheme for mortality may partly explain this bias as this scheme tends to dampen extremes.Preliminary tests with PISCES coupled to the upper trophic layer model Apex Predators ECOSystem Model (APECOSM) (Maury et al., 2007) produce a much greater spatial and temporal variability for mesozooplankton, especially in the high latitudes and along the continental margins.
Figure 16 shows Taylor diagrams for nutrients, oxygen, alkalinity and DIC.Overall, except for the carbonate system and iron, the model performs quite well, as expected from the comparison made in the previous section.The poorest agreement is found for both alkalinity and iron.For iron, the model tends to strongly underestimate the spatial variability, both at the surface and in the interior of the ocean.Through a reinspection of Fig. 8, we can see that this weak bias is not surprising.In particular, the gradients from the coastal regions to the open ocean are generally too small.This suggests that the sediment source of iron is too small and should either be increased and/or made more variable.For the carbonate system, the predicted spatial variability is overestimated, in particular in the interior of the ocean.In fact, the data distribution which has been used to produce the observed climatology is rather coarse (Key et al., 2004).As a consequence, the interpolation procedure strongly smooths the DIC and alkalinity distribution.Thus, the GLODAP database probably underestimates the real variability of these tracers.To avoid this problem, we should have used a non-interpolated data product as for iron or mesozooplankton.To estimate the potential uncertainty associated with the use of GLODAP, we have used another alkalinity database only available at the surface (Lee et al., 2006).The agreement between the model and this database is much better (see Fig. 16), thus confirming that interpolation in GLODAP potentially leads to a strong underestimate of the real spatial variability.

Sensitivity tests with some new parameterizations
A number of new parameterizations has been introduced in the current version of PISCES.The objective of this section is to briefly document the impact of some of these.To do so, we have run a series of sensitivity experiments for a duration of 10 years in which specific parameterizations have been either changed or removed.Table 7 summarizes the different experiments performed.The objective of these tests is not to unequivocally demonstrate that the new formulations improve the model skills but is rather to show the consequences of their utilization on the model behavior.

Dependence of growth rate on light
In the first two experiments, PAR and LIGHT, the sensitivity of the model results to the dependence of growth rate to light has been tested.In the PAR experiment, PAR is set as a constant fraction of incident shortwave radiation, here 43 %, as usually done in ocean biogeochemical models.Chlorophyll distribution is almost identical to the standard simulation (not shown).Furthermore, global primary production and export production remain almost unchanged (see Table 7).Model results are thus almost insensitive to the variability of the fraction of shortwave radiation that is PAR.In the second experiment, we use an alternative formulation of light limitation which corresponds to the standard parameterization as proposed by Geider et al. (1997) (see Eq. 2b).In this formulation, the light saturation parameter E k directly depends on temperature and nutrient limitation.Thus, since the Q 10 of phytoplankton is close to 2, E k is then predicted to be 6 to 8 times smaller in the very high latitudes than in the tropical Black dot corresponds to NO 3 , brown dot to O 2 , red dot to PO 4 , green dot to SiO 3 , light-blue dot to DIC, purple dot to alkalinity and gray dot to iron.The additional purple dot labeled as Alk-Lee uses the database constructed by Lee et al. (2006) to compare with the model.domain.Furthermore, in the very oligotrophic regions, such as the central subtropical gyres, E k is close to 0 as a consequence of a very intense nutrient limitation.In the LIGHT experiment, the initial slopes of P -I curves have been prescribed so that the resulting E k are identical to those of the standard case at 15 • C for no nutrient limitation.
Figure 17a and b show the difference in chlorophyll between the LIGHT experiment and the standard case for two seasons.The alternative parameterization of light limitation produces changes in surface chlorophyll at both seasons.In the very high latitudes of both hemispheres, surface chlorophyll is strongly increased during the corresponding growing season.The temperature dependence in the alternative parameterization produces lower light saturation parameters and thus, a weaker light limitation.On the contrary, in the mid-to high latitudes of both hemispheres, surface chlorophyll is significantly lower, especially in the Southern Ocean and in the Pacific Ocean.The temperature dependence of the light saturation parameter results in a weaker light limitation during winter.As a consequence, chlorophyll concentrations and primary productivity are predicted to be higher during this season generating a significant consumption and export of nutrients.At the beginning of the growing season, the stock of nutrients in the upper ocean is then lower which leads to weaker and shorter spring blooms.In the very high latitudes, the absence of light during winter and the presence of sea ice explain the different modeled response.In the low latitudes, the differences are relatively small.Surface chlorophyll concentrations tend to be higher in HNLC and productive regions.The alternative formulation tends to produce a stronger light limitation in the subsurface and thus, reduces the nutrient uptake below the surface.More iron and macronutrients are advected into the surface layer (not shown) which results in higher chlorophyll concentrations and in some cases, in larger productive regions (for instance in the tropical Atlantic Ocean and in the Arabian Sea).
Figure 18 shows the day at which blooms reach their maximum intensity in the Sea-Viewing Wide Field-of-View Senso (SeaWiFS) data, in the standard case and in LIGHT.Over the low and mid-latitudes as well as in the North Atlantic Ocean, the timing of the bloom maximum predicted by the standard model is in broad agreement with the satellite data.However, in the central part of the subarctic gyre of the North Pacific, the model simulates a bloom maximum which occurs much too early in the growing season, in January compared to August in the satellite observations.A similar bias is also predicted in part of the Southern Ocean, especially in the eastern part of the three sectors of this ocean.When the alternative parameterization of light limitation is used, the bloom timing remains unchanged over most of the ocean, except in the high latitudes in areas where the winter mixed layer remains relatively shallow.Such a result is not surprising because the alternative formulation predicts a much lower light saturation parameter in cold waters which alleviates light limitation at the beginning of the growing season.As a consequence, the bloom occurs earlier in the growing season, which tends to worsen the model behavior in the high latitudes of both hemispheres.In the North Pacific, the strong bias is not modified by the alternative formulation which suggests that this bias is not related to an incorrect description of light limitation.In fact, the model predicts a very strong limitation of phytoplankton growth by iron during summer and thus, simulated chlorophyll concentrations are very low.In winter, the mixed layer deepens supplying the surface with iron.However, it remains relatively shallow preventing thus phytoplankton from being severely light limited.Chlorophyll concentrations are then maximum during winter and minimum during summer, which is identical to what is observed in the subtropical gyres, at BATS for instance (Lévy et al., 2005;Fernández I et al., 2005).Yet, it is completely out of phase relative to the observations, suggesting that in that region, the model either strongly overestimates iron limitation during summer or that iron-light co-limitations are incorrectly parameterized in PISCES.
The sensitivity experiment presented here shows that model results are very sensitive to how light limitation is parameterized.Primary production, export production as well as the magnitude of the bloom are strongly impacted by the choice of the formulation describing light limitation of phytoplankton growth.The parameterization proposed by Geider et al. (1997) shares some similarities with the Liebig's law of the minimum.When nutrients are very limiting, light limitation becomes negligible since E k tends to 0. When light is strongly limiting, nutrients limitation becomes unimportant and growth rate becomes linearly related to light and Chl / C. The parameterization used in the standard case is similar to the multiplicative description of the limiting factors.As a consequence, the standard parameterization predicts lower phytoplankton growth rates, smaller primary production and less intense blooms.On the other hand, the timing of the bloom maximum is much less sensitive to the formulation of light limitation, except in the strongly stratified areas of the high latitudes.At low latitudes, light limitation at the surface is of secondary importance, despite that light limitation in the subsurface appears to partly control the amount of nutrients supplied to the surface.In the mid-and high latitudes, in areas characterized by deep winter mixed layers, the timing of the bloom maximum (but not its magnitude) appears to be virtually insensitive to the description of light limitation.This means that other factors, such as the timing of stratification, drive the timing of the bloom maximum.

Simple parameterization of cell size
In PISCES, a very basic parameterization of phytoplankton cell size has been developed to compute the values of the half-saturation coefficients for the different nutrients (see Eq. 7).This parameterization is based on the classical hypothesis, supported by observations, that the mean cell size of a phytoplankton community increases as the biomass increases (e.g., Raimbault et al., 1988;Armstrong, 1994;Hurtt and Armstrong, 1996).In the SIZE experiment, this simple parameterization has been removed, i.e., the half-saturation constants are kept constant to their minimum values as specified in Table 1.
Figure 17c and d display the differences in surface chlorophyll between the SIZE experiment and the standard configuration of the model.The largest differences are simulated in the high latitudes of both hemispheres, during the growing season.A closer inspection of the model results show that the largest changes occur at the end of the spring or summer bloom, when the exhaustion in nutrients becomes a major limiting factor.In the standard experiment, the cell-size parameterization produces high half-saturation constants during the phytoplankton bloom since they directly depend on the biomass level.Thus, nutrient limitation occurs earlier and is more severe leading to a shorter and less intense bloom.In the eastern boundary upwelling systems, the biomass is also very high.However, unlike in the high latitudes, the phyto-plankton biomass is mainly controlled by grazing so that nutrient concentrations are generally much higher than the values of the high saturation constants.In the subtropical oligotrophic gyres, the impact is negligible since the mean cell size is predicted to be at its minimal value in the standard experiment, which is equivalent to what is imposed in the SIZE experiment.The impact of the cell-size parameterization on nutrients is small, except for silicate in the equatorial Pacific Ocean (not shown).In this region, nanophytoplankton become strongly favored in the SIZE experiment because in the standard case, their cell size is not predicted to be minimum, whereas this is the case for diatoms.When the cell-size parameterization is removed, nanophytoplankton biomass increases and completely out compete diatoms.As a consequence, silicate consumption in the equatorial Pacific Ocean is strongly reduced which explains the simulated higher values in the SIZE experiment.However, the total chlorophyll concentration is nearly identical because the decrease in diatoms compensates for the increase in nanophytoplankton.Furthermore, the total chlorophyll biomass is regulated by the total supply in iron, whereas the contribution of the different phytoplankton species is driven by their competitive abilities (here specified by the values of their half-saturation constants).

Food quality and grazing
Food quality may have profound impacts on the grazing activity by zooplankton as discussed by Mitra et al. (2007).When absorbing prey with poor nutritional value, zooplankton may have two different options: (1) increase the retention time of the prey to extract as many metabolites as they can (Plath and Boersma, 2001), or (2) decrease the retention time of the preys to maintain the highest possible metabolite concentration in the digestive apparatus and thus to increase the probability to absorb valuable compounds (Tirelli and Mayzaud, 2005;Dutz et al., 2008).In the first case, growth efficiency is increased whereas it is decreased in the second case.In PISCES, poor food quality is assumed to impair gross growth efficiency (e Z ) of both microzooplankton and mesozooplankton based on the stoichiometric ratios of their preys (Fe / C and N / C, see Eq. 27).In the FOOD sensitivity experiment, the effect of food quality on the gross growth efficiency has been removed, i.e., e Z N is set to 1. Surface chlorophyll concentrations are almost unaltered when the impact of food quality is removed (see Fig. 17e and  f).The only noticeable differences are simulated from the equatorial Pacific Ocean where very strong iron limitation causes very low Fe / C ratios in phytoplankton.In the FOOD experiment, these low Fe / C ratios do not reduce zooplankton growth efficiency.Grazing pressure on phytoplankton is then higher.The nutrients distributions are also very close to those predicted in the standard experiment.Thus, food quality appears to have minimal consequences on chlorophyll and nutrients, at least in terms of their absolute values.
Figure 19 shows the relative changes in phytoplankton, microzooplankton and mesozooplankton biomasses (in carbon).A significant reduction in the carbon biomass of phytoplankton is predicted in the FOOD experiment.This reduction is maximum in the subtropical gyres where it may exceed 40 % because of more intense grazing by zooplankton.These changes are not perceptible in chlorophyll concentrations (at least with the color scale chosen on Fig. 17) because of the extremely low Chl / C in the gyres.Both on microzooplankton and mesozooplankton, the differences between the FOOD and the standard experiments are even more pronounced.Both zooplankton biomasses increase by more than 100 % in the subtropical gyres of all oceans and this increase even exceeds 200 % in the subtropical gyre of the South Pacific Ocean.
Food quality may thus have very important impacts on zooplankton, especially in the very oligotrophic regions.Furthermore, the importance of food quality is predicted to be more critical in regions depleted in nitrogen, characterized by very low N / C ratios in phytoplankton, than in iron limited areas.Several points may explain this greater sensitivity.First, even in the most severely iron limited areas, the Fe / C ratio in phytoplankton drops very rarely below half the value of the Fe / C ratio in zooplankton.In the central part of the subtropical gyres, where nitrogen limitation is the most intense, N / C ratios in phytoplankton can reach 0.04, that is about 3 times less than the N / C ratio of zooplankton.Second, the available food in the intense oligotrophic areas is much lower than in the iron limited regions.Chlorophyll concentrations in the typical HNLC regions are generally around 0.2 to 0.3 mg chl a m −3 , whereas it is below 0.1 mg chl a m −3 in the subtropical gyres.As a consequence, zooplankton biomass is lower in the subtropical gyres which increases the magnitude of the relative changes.

Conclusions
In this paper, we have presented a full and thorough description of the current state of the ocean biogeochemical model PISCES, called PISCES-v2.Since the latest published version of the model (Aumont and Bopp, 2006), PISCES-v2 has undergone major changes both in terms of the modeled processes and of the model structure and performance.Relative to its previous version PISCES-v1, key changes are a ma-jor redesign of phytoplankton growth description, including a quota-based representation of iron limitation, an improvement of the zooplankton compartment, a better description of the benthic processes and a simple description of nitrogen fixation by diazotrophs.A complete list of the changes made in PISCES-v2 relative to its previously published version is detailed in Sect. 2. The performance of the model has been then evaluated using a climatological simulation run to quasi-steady state.The model produces reasonable surface distributions of chlorophyll, mesozooplankton and nutrients (including iron) and simulates consistent vertical distributions of the main biogeochemical tracers.Some of the main deficiencies of the model are the spatial distribution of the oxygen minimum zones, the silicic acid distribution in the Southern Ocean, too elevated nutrients concentrations in the deep Atlantic Ocean and an out-of-phase predicted seasonal cycle of chlorophyll in the subarctic Pacific Ocean.
PISCES includes several optional parameterizations that may be activated from the namelist.In this study, we have presented the impacts of some of these optional formulations evaluated in a set of sensitivity experiments.The choice of the light limitation scheme has the largest effect on the model solution, especially on chlorophyll.The amplitude of the seasonal cycle in the high latitudes is profoundly impacted whereas the timing of the bloom maximum is in general only very moderately altered.The effect of food quality on the growth efficiency of zooplankton has been shown to lead to important relative changes in the oligotrophic subtropical gyres.The model suggests that it is critical to maintain sufficiently high chlorophyll levels in these regions.It may also contribute to, at least partly, explaining the too low primary productivity simulated by other biogeochemical models in the subtropical gyres (Yool et al., 2013).
The description of PISCES presented here has been restricted to the core scheme which can be obtained online from different SVN repositories depending on the dynamical framework in which it is embedded (see the Introduction for a list of theses repositories).In addition to the description of the lower trophic levels of marine ecosystems, and the biogeochemical cycles of carbon and of the main nutrients (N, P, Si, Fe), as described in this manuscript, a few additional modules have been embedded into PISCES.These modules enable the model to compute the cycles of climate-relevant gases emitted by the ocean such as dimethylsulfide (DMS) (Bopp et al., 2008), and nitrous oxide (N 2 O) (Martinez-Rey et al., 2015).An explicit representation of paleo-proxies, such as δ 13 C (Tagliabue et al., 2009b), Pa / Th (Dutay et al., 2009), Nd (Arsouze et al., 2009), is also available.
PISCES is still in a phase of active development despite the fact that its development started more than 10 years ago.Avenues for future improvements are large and numerous and concern all aspects of the model.The challenges confronting marine biogeochemical modeling have been identified in many dedicated studies (e.g., Doney, 1999;Hood et al., 2006;Merico et al., 2009;Smith et al., 2011;Mitra et al., 2014).Setting priorities in a long list of potential necessary modifications is a rather difficult task which relies not only on the diagnostic of the major deficiencies of the current model but also on the future research scope envisioned for the model.In the coming years, PISCES will evolve along two main avenues.First, a more sophisticated treatment of phytoplankton physiology will replace the current relatively simple scheme.A main consequence is the representation of variable elemental ratios for all major elements (N, P, Fe, Si, C).Redfield-Monod models have been shown to exhibit serious deficiencies which advocate for their replacement by more detailed mechanistic schemes (Flynn, 2010;Smith et al., 2011).Second, almost all marine biogeochemical mod-els have been built on the classical distinction between phytoplanktonic autotrophic organisms and zooplanktonic heterotrophic organisms.However, this dichotomy has been increasingly challenged in the recent years as observations have shown that most protists, probably with the exception of diatoms, have to a lesser or greater degree a mixotrophic status (e.g., Stoecker, 1998;Flynn et al., 2013).The conceptual schemes on which biogeochemical models, including PISCES, should then be revised, in particular the distinction between phytoplankton and microzooplankton.

Figure 1 .
Figure 1.Architecture of PISCES.This figure only shows the ecosystem model omitting thus oxygen and the carbonate system.The elements which are explicitly modeled are indicated in the left corner of each box.
1. (a) Model parameters for phytoplankton with their default values in PISCES.(b) Model parameters for zooplankton with their default values in PISCES.(c) Model parameters for DOM with their default values in PISCES.(d) Model parameters for particulate organic and inorganic matter with their default values in PISCES.(e) Model parameters for various processes with their default values in PISCES.

Figure 2 .
Figure2.Reduction of growth rate when the mixed layer depth exceeds the euphotic depth for nanophytoplankton (continuous line) and diatoms (dashed line).Depth corresponds to Z.

Figure 3 .
Figure 3. θ Si,D opt as a function of Si concentration and F D Si lim,1 .The vertical axis corresponds to log(Si).

FeFigure 4 .
Figure 4. Dissolution rate of PSI (λ PSi ) normalized to its value at 0 • C with no silicate.Temperature is in • C.

a
Carbon fluxes are all in Gt C yr −1 .b The total vertical flux due to sinking POC is 7.3 Gt C yr −1 at 100 m depth.

a
Calcite fluxes are all in Gt C yr −1 .b Biogenic silica fluxes are all in Tmol Si yr −1 .byPISCES, quite similar to the value of 10.4 ± 3.7 Gt C yr −1 estimated byHernández-León and Ikeda (2005) for the global respiration of mesozooplankton in the upper 200 m of the ocean.About 80 % of total primary production, i.e., 35.8 Gt C yr −1 , is consumed up by microzooplankton above the upper bound of the 25-33 Gt C yr −1 given byBuitenhuis  et al. (

Figure 5 .
Figure 5. Sediment source of iron as a function of depth.This plot displays the vertical variation of Fsed (see Eq. 85 for the definition of this factor).

Figure 7 .
Figure 7. Surface seasonal mean chlorophyll concentrations (mg chl a m −3 ) in April-May-June (panels a and c) and November-December-January (panels b and d).Panels (a) and (b) display satellite observations from GLOBCOLOUR.Panels (c) and (d) are model results.

Figure 8 .
Figure 8. Spatial distribution of annual-mean iron concentrations (in nmol L −1 ) as observed (left column) and as simulated by PISCES (right column).On panels (a) and (b), iron has been averaged over the top 50 m of the ocean.On panels (b) and (c), iron has been averaged over 200-1000 m.The bottom two panels display the iron distributions average over the depth range 1000-5000 m.Model values have been sampled at the same location and month as the data.

Figure 9 .
Figure 9. Annual-mean NO 3 concentrations in µmol N L −1 .Observations are from the World Ocean Atlas 2009 (Garcia et al., 2010).(a) Observed surface.(b) Model run surface.(c) Observed transect zonally averaged over the Atlantic.(d) Same as (c) but for the model.(e) Observed transect zonally averaged over the Pacific.(f) Same as (e) but for the model.

Figure 12 .
Figure 12.Annual-mean natural DIC concentrations in µmol L −1 .Observations are from GLODAP.The pre-industrial distribution of DIC has been estimated in GLODAP as the difference between total DIC and anthropogenic carbon.Panels are the same as on Fig. 9.

Figure 13 .
Figure 13.Annual-mean alkalinity concentrations in µmol eq L −1 .Observations are from GLODAP.Panels are the same as on Fig. 9.

Figure 14 .
Figure 14.Taylor diagrams of model-observation comparisons for surface chlorophyll (log10-transformed) using monthly mean fields (a) and annual-mean fields (b).Black dot corresponds to global comparison; red dot to the Atlantic Ocean, green dot to the Pacific Ocean, brown dot to the Indian Ocean and gray dot to the Southern Ocean (south of 45 • S).

Figure 15 .
Figure 15.Taylor diagram of model-observation comparisons for mesozooplankton using monthly mean fields.Data come from the Green Ocean Project web site.Black dot corresponds to the global ocean; red dot to the Atlantic Ocean, green dot to the Pacific Ocean, brown dot to the Indian Ocean and gray dot to the Southern Ocean (south of 45 • S).

Figure 16 .
Figure 16.Taylor diagrams of model-observation comparisons for nutrients using monthly mean fields.The data are identical to those used in previous plots.Panel (a) corresponds to the global ocean.Panel (b) shows the comparison restricted to the top 100 m of the ocean.Black dot corresponds to NO 3 , brown dot to O 2 , red dot to PO 4 , green dot to SiO 3 , light-blue dot to DIC, purple dot to alkalinity and gray dot to iron.The additional purple dot labeled as Alk-Lee uses the database constructed byLee et al. (2006) to compare with the model.

Figure 17 .
Figure 17.Surface seasonal mean chlorophyll anomaly (mg chl a m −3 ) relative to the standard simulation in April-May-June (left column) and November-December-January (right column).Panels (a) and (b) correspond to the LIGHT test; panels (c) and (d) show to the SIZE test; panels (e) and (f) display the FOOD test.

Figure 18 .
Figure 18.Day of the year at which sea surface chlorophyll is maximum.Panel (a) corresponds to the observations; panel (b) displays the standard simulation.Panel (c) shows the difference between the LIGHT and the standard experiments.Only the regions where the amplitude of the seasonal cycle exceeds 0.1 mg chl a m −3 are shown.

Figure 19 .
Figure 19.Annual-mean relative change in the surface carbon biomass of total phytoplankton (panel a), microzooplankton (panel b), and mesozooplankton (panel c) in the FOOD experiment compared to the standard case.
Model parameters for phytoplankton with their default values in PISCES.(b) Model parameters for zooplankton with their default values in PISCES.(c) Model parameters for organic and inorganic matter with their default values in PISCES.(d) Model parameters for various processes with their default values in PISCES.

Table 2 .
Boolean variables in the namelist.These variables activate functionalities of PISCES.ln_check_mass Check mass conservation (T) * The frequency at which the restoring technique is applied is specified by the parameter nn_pisdmp. *

Table 3 .
Available CPP keys in PISCES.

Table 4 .
Global annual budget of C in the top 150 m of the ocean.

Table 5 .
Global annual budget of calcite and Si in the top 150 m of the ocean.

Table 6 .
Annual budget * of N over the global ocean.Net budget of nitrogen (Sources minus Sinks) * All nitrogen fluxes are in Tg N yr −1 .

Table 7 .
Sensitivity experiments performed with PISCES to evaluate the impact of specific parameterizations.Primary production (PP) and export production at 150 m (EP) are in Gt C yr −1 .