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

Abstract

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${}_{\mathrm{2}}$ by more than 100 $\mathrm{ppmv}$. The oceans play a major role in the carbon cycle and in its adjustment. 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.,. 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 based on a very simple description of the carbon cycle, the number and the complexity of models have rapidly increased e.g.,. However, 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.,.

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 as opposed to the quota models 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., 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 . 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) . 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 , as it was based on HAMOCC3.1 and used in the LSG model . When this code was embedded in the ocean model Ocean PArallélisé (OPA) , 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 improved 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 to operational oceanography . PISCES has been used to analyze intraseasonal to interannual and decadal timescales . 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 . Several studies have been conducted that consider the potential impact of climate change on ocean biogeochemistry . Modeling studies focusing on paleoceanography have been based on PISCES . Finally, PISCES is also used in regional configurations to study specific regions such as the Peru upwelling or the Indian Ocean .

PISCES is currently embedded into two modeling systems: NEMO and ROMS_AGRIF . It can be downloaded from their respective web sites:

• http://www.nemo-ocean.eu for the NEMO ocean modeling framework;

• http://www.romsagrif.org for the ROMS_AGRIF modeling framework.

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 . 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:

  • 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.

  • 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.

  • Memory and performance improvements have been made. This version should run slightly faster and take much less memory than v1.

  • The namelist now includes many more parameters that may thus be changed without recompiling the code.

• Changes made to the nutrients:

  • 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.

  • Scavenging of inorganic iron and coagulation of iron colloids have been redesigned.

• Changes made to the phytoplankton compartments:

  • Nutrients limitation terms now include a simple description of the impact of cell size.

  • 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.

  • Silicification, calcification as well as nitrogen fixation are redesigned by diazotrophs.

  • The relationship between growth rate (primary production) and light can be chosen between two different formulations.

• Changes made to the zooplankton compartments:

  • The microzooplankton grazing formulation is now identical to that of mesozooplankton.

  • Thresholds can be selected for both total food or individual prey types.

  • Food quality affects the gross growth efficiency of both zooplankton compartments.

• Changes made to dissolved organic matter and particulate materials:

  • Two different schemes for the description of particulate organic matter can be chosen: the traditional two-compartment model or the Kriest model.

  • Bacterial implicit description has been redesigned.

  • Dissolution of biogenic silica assumes two different fractions.

  • The dust distribution in the water column is modeled using a very crude parameterization.

  • 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:

  • Spatially variable solubility of iron in dust can be specified from a file.

  • River discharge of nutrients has been improved.

  • 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

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.

[Figure omitted. See PDF]

PISCES currently has 24 compartments (see Fig. ). 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 $/$ ${\mathrm{O}}_{\mathrm{2}}$ ratios are assumed constant and are not allowed to vary. In PISCES, the Redfield ratios C $/$ N $/$ P are set to $122/16/\mathrm{1}$ and the $-$O $/$ C ratio is set to 1.34 . 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”, ). 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 but it has been abandoned since as observations do not suggest a clear relationship between sinking speeds and mineral composition of particles . 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 ($\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$) except for the silicon, chlorophyll and iron content of phytoplankton, which are respectively in Si, Chl and Fe units ($\mathrm{mol}\mathrm{Si}{\mathrm{L}}^{-\mathrm{1}}$, $\mathrm{g}\mathrm{Chl}{\mathrm{L}}^{-\mathrm{1}}$, and $\mathrm{mol}\mathrm{Fe}{\mathrm{L}}^{-\mathrm{1}}$, respectively). Finally, all parameters and their standard values in PISCES are listed in Tables a–e at the end of this section.

(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.

(a)
Parameter	Units	Value	Description
${\mathit{\mu }}_{\text{max}}^{\mathrm{0}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.6	Growth rate at 0 ${}^{\circ }\mathrm{C}$
${\mathit{\mu }}_{\text{ref}}$	${\mathrm{d}}^{-\mathrm{1}}$	1.0	Growth rate reference for light limitation
$b_{\text{resp}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.033	Basal respiration rate
$b_P$	–	1.066	Temperature sensitivity of growth
${\mathit{\alpha }}^I$	$(\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	2; 2	Initial slope of $P$–$I$ curve
${\mathit{\delta }}^I$	–	0.05; 0.05	Exudation of DOC
${\mathit{\beta }}_{\mathrm{1}}^I$	–	2.1; 1.6	Absorption in the blue part of light
${\mathit{\beta }}_{\mathrm{2}}^I$	–	0.42; 0.69	Absorption in the green part of light
${\mathit{\beta }}_{\mathrm{3}}^I$	–	0.4; 0.7	Absorption in the red part of light
$K_{{\mathrm{PO}}_{\mathrm{4}}}^{I,\text{min}}$	$\mathrm{nmol}\mathrm{P}{\mathrm{L}}^{-\mathrm{1}}$	0.8; 2.4	Minimum half-saturation constant for phosphate
$K_{{\mathrm{NH}}_{\mathrm{4}}}^{I,\text{min}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}$	0.013; 0.039	Minimum half-saturation constant for ammonium
$K_{{\mathrm{NO}}_{\mathrm{3}}}^{I,\text{min}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}$	0.13; 0.39	Minimum half-saturation constant for nitrate
$K_{\mathrm{Si}}^{D,\text{min}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{Si}{\mathrm{L}}^{-\mathrm{1}}$	1	Minimum half-saturation constant for silicate
$K_{\mathrm{Si}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{Si}{\mathrm{L}}^{-\mathrm{1}}$	16.6	Parameter for the half-saturation constant
$K_{\mathrm{Si}}^I$	$\mathrm{\mu }\mathrm{mol}\mathrm{Si}{\mathrm{L}}^{-\mathrm{1}}$	2; 20	Parameters for Si $/$ C
$K_{\mathrm{Fe}}^{I,\text{min}}$	$\mathrm{nmol}\mathrm{Fe}{\mathrm{L}}^{-\mathrm{1}}$	1; 3	Minimum half-saturation constant for iron uptake
$S_{\text{rat}}^I$	–	3; 3	Size ratio of Phytoplankton
${\mathit{\theta }}_{\mathrm{m}}^{\mathrm{Si},D}$	$\mathrm{mol}\mathrm{Si}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	0.159	Optimal Si $/$ C uptake ratio of diatoms
${\mathit{\theta }}_{\text{opt}}^{\mathrm{Fe},I}$	$\mathrm{\mu }\mathrm{mol}\mathrm{Fe}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	7; 7	Optimal iron quota
${\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},I}$	$\mathrm{\mu }\mathrm{mol}\mathrm{Fe}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	40; 40	Maximum iron quota
$m^I$	${\mathrm{d}}^{-\mathrm{1}}$	0.01; 0.01	phytoplankton mortality rate
$w^P$	${\mathrm{d}}^{-\mathrm{1}}\mathrm{mol}{\mathrm{C}}^{-\mathrm{1}}$	0.01	Minimum quadratic mortality of phytoplankton
$w_{\text{max}}^D$	${\mathrm{d}}^{-\mathrm{1}}\mathrm{mol}{\mathrm{C}}^{-\mathrm{1}}$	0.03	Maximum quadratic mortality of diatoms
${\mathit{\theta }}_{\text{max}}^{\text{Chl},I}$	$\mathrm{mg}\mathrm{Chl}(\mathrm{mg}\mathrm{C}{)}^{-\mathrm{1}}$	0.033; 0.05	Maximum Chl $/$ C ratios of phytoplankton
${\mathit{\theta }}_{\text{min}}^{\text{Chl}}$	$\mathrm{mg}\mathrm{Chl}(\mathrm{mg}\mathrm{C}{)}^{-\mathrm{1}}$	0.0033	Minimum Chl $/$ C ratios of phytoplankton
$I_{\text{max}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	1; 1	Threshold concentration for size dependency

Continued.

(b)
Parameter	Units	Value	Description
$b_Z$	–	1.079; 1.079	Temperature sensitivity term
$e_{\text{max}}^I$	–	0.3; 0.35	Maximum growth efficiency of zooplankton
${\mathit{\sigma }}^I$	–	0.3; 0.3	Non-assimilated fraction
${\mathit{\gamma }}^I$	–	0.6; 0.6	Excretion as DOM
$g_{\mathrm{m}}^I$	${\mathrm{d}}^{-\mathrm{1}}$	3; 0.75	Maximum grazing rate
$g_{\text{FF}}^M$	$(\mathrm{m}\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$	$\mathrm{2}\times {10}^{\mathrm{3}}$	Flux feeding rate
$K_{\mathrm{G}}^I$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	20; 20	Half-saturation constant for grazing
$p_P^I$	–	1; 0.3	Preference for nanophytoplankton
$p_D^I$	–	0.5; 1	Preference for diatoms
$p_{\text{POC}}^I$	–	0.1,0.3	Preference for POC
$p_Z^M$	–	1.0	Preference for microzooplankton
$F_{\text{thresh}}^I$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	0.3; 0.3	Food threshold for zooplankton
$J_{\text{thres}}^Z$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	0.001	Specific food thresholds for microzooplankton
$J_{\text{thres}}^M$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	0.001	Specific food thresholds for mesozooplankton
$m^I$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	0.004; 0.03	Zooplankton quadratic mortality
$r^I$	${\mathrm{d}}^{-\mathrm{1}}$	0.03,0.005	Zooplankton linear mortality
$K_{\mathrm{m}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	0.2	Half-saturation constant for mortality
${\mathit{\nu }}^I$	–	0.5; 0.75	Fraction of calcite that does not dissolve in guts
${\mathit{\theta }}^{\mathrm{Fe},\mathrm{Zoo}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{Fe}\mathrm{mol}{\mathrm{C}}^{-\mathrm{1}}$	10	Fe $/$ C ratio of zooplankton

Continued.

(c)
Parameter	Units	Value	Description
${\mathit{\lambda }}_{\text{DOC}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.3	Remineralization rate of DOC
$K_{\text{DOC}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}$	417	Half-saturation constant for DOC remin.
$K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{Bact}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}$	0.03	${\mathrm{NO}}_{\mathrm{3}}$ half-saturation constant for DOC remin.
$K_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{Bact}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}$	0.003	${\mathrm{NH}}_{\mathrm{4}}$ half-saturation constant for DOC remin.
$K_{{\mathrm{PO}}_{\mathrm{4}}}^{\text{Bact}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{P}{\mathrm{L}}^{-\mathrm{1}}$	0.003	${\mathrm{PO}}_{\mathrm{4}}$ half-saturation constant for DOC remin.
$K_{\mathrm{Fe}}^{\text{Bact}}$	$\mathrm{nmol}\mathrm{Fe}{\mathrm{L}}^{-\mathrm{1}}$	0.01	Fe half-saturation constant for DOC remin.
$a_{\mathrm{1}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	0.37	Aggregation rate (turbulence) of DOC$\to$POC
$a_{\mathrm{2}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	102	Aggregation rate (turbulence) of DOC$\to$POC
$a_{\mathrm{3}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	3530	Aggregation rate (turbulence) of DOC$\to$GOC
$a_{\mathrm{4}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	5095	Aggregation rate (Brownian) of DOC$\to$POC
$a_{\mathrm{5}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	114	Aggregation rate (Brownian) of DOC$\to$POC

Continued.

(d)
Parameter	Units	Value	Description
${\mathit{\lambda }}_{\text{POC}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.025	Degradation rate of POC
$w_{\text{POC}}$	$\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$	2	Sinking speed of POC
$w_{\text{GOC}}^{\text{min}}$	$\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$	30	Minimum sinking speed of GOC${}_b$
$w_{\text{dust}}$	$\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$	2	Sinking speed of dust
$a_{\mathrm{6}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	25.9	Aggregation rate (turbulence) of POC$\to$GOC
$a_{\mathrm{7}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	4452	Aggregation rate (turbulence) of POC$\to$GOC
$a_{\mathrm{8}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	3.3	Aggregation rate (settling) of POC$\to$GOC
$a_{\mathrm{9}}$	$(\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}}{)}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	47.1	Aggregation rate (settling) of POC$\to$GOC
${\mathit{\lambda }}_{\mathrm{Fe}}^{\text{min}}$	${\mathrm{d}}^{-\mathrm{1}}$	$\mathrm{3}\times {10}^{-\mathrm{5}}$	Minimum scavenging rate of iron
${\mathit{\lambda }}_{\mathrm{Fe}}$	${\mathrm{d}}^{-\mathrm{1}}\mathrm{\mu }{\mathrm{mol}}^{-\mathrm{1}}\mathrm{L}$	0.005	Slope of the scavenging rate of iron xlam1
${\mathit{\lambda }}_{\mathrm{Fe}}^{\text{dust}}$	${\mathrm{d}}^{-\mathrm{1}}{\mathrm{mg}}^{-\mathrm{1}}\mathrm{L}$	150	Scavenging rate of iron by dust
${\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.197	Dissolution rate of calcite
nca	–	1	Exponent in the dissolution rate of calcite
${\mathit{\chi }}_{\text{lab}}^{\mathrm{0}}$	–	0.5	Proportion of the most labile phase in $\mathrm{PSi}$
${\mathit{\lambda }}_{\mathrm{PSi}}^{\text{slow}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.003	Slow dissolution rate of $\mathrm{BSi}$
${\mathit{\lambda }}_{\mathrm{PSi}}^{\text{fast}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.025	Fast dissolution rate of $\mathrm{BSi}$

Continued.

(e)
Parameter	Units	Value	Description
${\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}$	${\mathrm{d}}^{-\mathrm{1}}$	0.05	Maximum nitrification rate
${{\mathrm{O}}_{\mathrm{2}}}^{\text{min,1}}$	$\mathrm{\mu }\mathrm{mol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}$	1	Half-saturation constant for denitrification
${{\mathrm{O}}_{\mathrm{2}}}^{\text{min,2}}$	$\mathrm{\mu }\mathrm{mol}{\mathrm{O}}_{\mathrm{2}}{\mathrm{L}}^{-\mathrm{1}}$	6	Half-saturation constant for denitrification
$L_{\mathrm{T}}$	$\mathrm{nmol}{\mathrm{L}}^{-\mathrm{1}}$	0.6	Total concentration of iron ligands
$N_{\text{fix}}^{\mathrm{m}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{d}}^{-\mathrm{1}}$	0.013	Maximum rate of nitrogen fixation
$K_{\mathrm{Fe}}^{\text{Dz}}$	$\mathrm{nmol}\mathrm{Fe}{\mathrm{L}}^{-\mathrm{1}}$	0.1	Fe half-saturation constant of nitrogen fixation
$E_{\text{fix}}$	$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$	50	Photosynthetic parameter of nitrogen fixation
${\mathrm{Fe}}_{\text{ice}}$	$\mathrm{nmol}\mathrm{Fe}{\mathrm{L}}^{-\mathrm{1}}$	15	iron concentration in sea ice
$F_{\mathrm{Fe},\min }^{\text{sed}}$	$\mathrm{\mu }\mathrm{mol}\mathrm{Fe}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$	2	Maximum sediment flux of Fe
${\mathrm{Sol}}_{\mathrm{Fe}}$	–	0.02	Solubility of iron in dust
${\mathrm{O}}_{\mathrm{2}}^{\text{ut}}$	$\mathrm{mol}{\mathrm{O}}_{\mathrm{2}}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	133$/$122	O $/$ C for ammonium-based processes
${\mathrm{O}}_{\mathrm{2}}^{\text{nit}}$	$\mathrm{mol}{\mathrm{O}}_{\mathrm{2}}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	32$/$122	O $/$ C ratio of nitrification
$r_{{\mathrm{NH}}_{\mathrm{4}}}^{\star }$	$\mathrm{mol}\mathrm{N}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	3$/$5	C$/$N ratio of ammonification
$r_{{\mathrm{NO}}_{\mathrm{3}}}^{\star }$	$\mathrm{mol}\mathrm{N}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	105$/$16	C$/$N ratio of denitrification
${\mathit{\theta }}^{\mathrm{N},\mathrm{C}}$	$\mathrm{mol}\mathrm{N}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$	16/122	N $/$ C Redfield ratio
$r_{{\mathrm{CaCO}}_{\mathrm{3}}}$	–	0.3	Rain-ratio parameter

Phytoplankton

Nanophytoplankton

$$\begin{array}{cc}& \frac{\partial P}{\partial t}=(\mathrm{1}-{\mathit{\delta }}^P){\mathit{\mu }}^{P}P-m^P\frac{P}{K_{\mathrm{m}}+P}P-\text{sh}\times w^P{P}^{\mathrm{2}}\\ & & -g^Z\left(P\right)Z-g^M\left(P\right)M\end{array}$$

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 (${\mathit{\mu }}^P$) can be computed according to two different parameterizations: $$\begin{array}{cc}& {\mathit{\mu }}^P={\mathit{\mu }}_P{f}_{\mathrm{1}}\left(L_{\text{day}}\right)f_{\mathrm{2}}\left(z_{\text{mxl}}\right)\\ & & \left(\mathrm{1},-,\exp ,\left(\frac{-{\mathit{\alpha }}^P{\mathit{\theta }}^{\text{Chl}{,}^P}{\text{PAR}}^P}{L_{\text{day}}({\mathit{\mu }}_{\text{ref}}+b_{\text{resp}})}\right)\right)L_{\text{lim}}^P,& {\mathit{\mu }}^P={\mathit{\mu }}_P{f}_{\mathrm{1}}\left(L_{\text{day}}\right)f_{\mathrm{2}}\left(z_{\text{mxl}}\right)\\ & & \left(\mathrm{1},-,\exp ,\left(\frac{-{\mathit{\alpha }}^P{\mathit{\theta }}^{\text{Chl}{,}^P}{\text{PAR}}^P}{L_{\text{day}}{\mathit{\mu }}_P{L}_{\text{lim}}^P}\right)\right)L_{\text{lim}}^P,\end{array}$$ where $b_{\text{resp}}$ is a small respiration rate and ${\mathit{\mu }}_{\text{ref}}$ a reference growth rate, independent of temperature. All other 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. () is used. In the previous equations, $L_{\text{day}}$ is day length ($\in [\mathrm{0},\mathrm{1}]$). $f_{\mathrm{1}}\left(L_{\text{day}}\right)$ expresses the dependency of growth rate to the length of the day . $z_{\text{mxl}}$ is the depth of the mixed layer and $f_{\mathrm{2}}\left(z_{\text{mxl}}\right)$ imposes an additional reduction of the growth rate when the mixed layer depth exceeds the euphotic depth: $$\begin{array}{}& & f_{\mathrm{1}}\left(L_{\text{day}}\right)=1.5\frac{L_{\text{day}}}{0.5+L_{\text{day}}},& & \mathrm{\Delta }z=max\left(\mathrm{0},,,z_{\text{mxl}},-,z_{\text{eu}}\right),& & t_{\text{dark}}=(\mathrm{\Delta }z{)}^{\mathrm{2}}/\mathrm{86\hspace{0.17em}400},& & f_{\mathrm{2}}\left(z_{\text{mxl}}\right)=\mathrm{1}-\frac{t_{\text{dark}}}{t_{\text{dark}}^P+t_{\text{dark}}},\end{array}$$ where $z_{\text{eu}}$ is the depth of the euphotic zone defined as the depth at which there is 1 $\mathrm{‰}$ of surface PAR. $t_{\text{dark}}^P$ is set to 3 $\mathrm{days}$ for nanophytoplankton and 4 $\mathrm{days}$ for diatoms, as diatoms generally better cope with prolonged dark periods. $t_{\text{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 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$; 86 400 converts $t_{\text{dark}}$ from ${\mathrm{s}}^{-\mathrm{1}}$ to ${\mathrm{day}}^{-\mathrm{1}}$. Figure  displays $f_{\mathrm{2}}\left(z_{\text{mxl}}\right)$ as a function of $\mathrm{\Delta }z$.

Reduction of growth rate when the mixed layer depth exceeds the euphotic depth for nanophytoplankton (continuous line) and diatoms (dashed line). Depth corresponds to $\mathrm{\Delta }Z$.

[Figure omitted. See PDF]

${\mathit{\mu }}_P$ is defined as follows : $$\begin{array}{}& & f_P\left(T\right)=b_P^T,& & {\mathit{\mu }}_P={\mathit{\mu }}_{\text{max}}^{\mathrm{0}}{f}_P\left(T\right).\end{array}$$

In PISCES, vertical penetration of the photosynthetic available radiation (PAR) is based on a simplified version of the model by , which is described in . Visible light is split into three wavebands: blue (400–500 $\mathrm{nm}$), green (500–600 $\mathrm{nm}$) and red (600–700 $\mathrm{nm}$). For each waveband, the chlorophyll-dependent attenuation coefficients are fitted to the coefficients computed from the full spectral model of (as modified in ) 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). $$\begin{array}{}& & {\text{PAR}}_{\mathrm{1}}\left(\mathrm{0}\right)={\text{PAR}}_{\mathrm{2}}\left(\mathrm{0}\right)={\text{PAR}}_{\mathrm{3}}\left(\mathrm{0}\right)=\frac{{\mathit{\rho }}_{\text{par}}}{\mathrm{3}}\text{SW}& & {\text{PAR}}^P\left(z\right)={\mathit{\beta }}_{\mathrm{1}}^P{\text{PAR}}_{\mathrm{1}}\left(z\right)+{\mathit{\beta }}_{\mathrm{2}}^P{\text{PAR}}_{\mathrm{2}}\left(z\right)+{\mathit{\beta }}_{\mathrm{3}}^P{\text{PAR}}_{\mathrm{3}}\left(z\right)\end{array}$$

Light absorption by phytoplankton depends on the waveband and on the species. The normalized coefficients ${\mathit{\beta }}_i$ have been computed for each phytoplankton group by averaging and normalizing, for each waveband, the absorption coefficients published in .

In PISCES, the nutrient limitation terms are defined as follows: $$\begin{array}{}& & L_{\text{lim}}^P=min\left(L_{{\mathrm{PO}}_{\mathrm{4}}}^P,,,L_{\mathrm{N}}^P,,,L_{\mathrm{Fe}}^P\right),& & L_{{\mathrm{PO}}_{\mathrm{4}}}^P=\frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K_{{\mathrm{PO}}_{\mathrm{4}}}^P},& & L_{\mathrm{N}}^P=L_{{\mathrm{NO}}_{\mathrm{3}}}^P+L_{{\mathrm{NH}}_{\mathrm{4}}}^P,& & L_{{\mathrm{NH}}_{\mathrm{4}}}^P=\frac{K_{{\mathrm{NO}}_{\mathrm{3}}}^P{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NO}}_{\mathrm{3}}}^P{K}_{{\mathrm{NH}}_{\mathrm{4}}}^P+K_{{\mathrm{NH}}_{\mathrm{4}}}^P{\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}}^P{\mathrm{NH}}_{\mathrm{4}}},& & L_{{\mathrm{NO}}_{\mathrm{3}}}^P=\frac{K_{{\mathrm{NH}}_{\mathrm{4}}}^P{\mathrm{NO}}_{\mathrm{3}}}{K_{{\mathrm{NO}}_{\mathrm{3}}}^P{K}_{{\mathrm{NH}}_{\mathrm{4}}}^P+K_{{\mathrm{NH}}_{\mathrm{4}}}^P{\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}}^P{\mathrm{NH}}_{\mathrm{4}}},& & L_{\mathrm{Fe}}^P=min\left(\mathrm{1},,,max,\left(\mathrm{0},,,\frac{{\mathit{\theta }}^{\mathrm{Fe},P}-{\mathit{\theta }}_{\text{min}}^{\mathrm{Fe},P}}{{\mathit{\theta }}_{\text{opt}}^{\mathrm{Fe},P}}\right)\right).\end{array}$$

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 (${\mathit{\theta }}_{\text{opt}}^{\mathrm{Fe},P}$) is used in the denominator which allows luxury uptake as in the model proposed by .

The choice of the half-saturation constants is rather difficult as observations show that they can vary by several orders of magnitude e.g.,. 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) . 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.,: $$\begin{array}{}& & P_{\mathrm{1}}=min\left(P,,,P_{\text{max}}\right),& & P_{\mathrm{2}}=max\left(\mathrm{0},,,P,-,P_{\text{max}}\right),& & K_i^P=K_i^{P,\text{min}}\frac{P_{\mathrm{1}}+S_{\text{rat}}^P{P}_{\mathrm{2}}}{P_{\mathrm{1}}+P_{\mathrm{2}}},\end{array}$$ where $S_{\text{rat}}^P$ is the size ratio of the larger size class over the smaller size class. $K_i^{P,\text{min}}$ is the half-saturation constant of the smaller size class. This parameterization assumes that half-saturation constants increase linearly with size . 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 . However, in a recent review, 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_{\text{max}}$, $K_i^{P,\text{min}}$, and $S_{\text{rat}}^P$) 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 . This acclimation mechanism is not included in PISCES, except for the case of silicate (see Sect. ).

The distinction between new production based on nitrate and regenerated production based on ammonium is computed as follows : $${\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^P={\mathit{\mu }}^P\frac{L_{{\mathrm{NO}}_{\mathrm{3}}}^P}{L_{{\mathrm{NO}}_{\mathrm{3}}}^P+L_{{\mathrm{NH}}_{\mathrm{4}}}^P}{\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^P={\mathit{\mu }}^P\frac{L_{{\mathrm{NH}}_{\mathrm{4}}}^P}{L_{{\mathrm{NO}}_{\mathrm{3}}}^P+L_{{\mathrm{NH}}_{\mathrm{4}}}^P},$$ where ${\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^P$ and ${\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^P$ 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 ${\mathrm{s}}^{-\mathrm{1}}$ in the mixed layer and to 0.01 ${\mathrm{s}}^{-\mathrm{1}}$ below. This means that the aggregation is reduced by a factor of 100 below the mixed layer.

Diatoms

$$\begin{array}{cc}\frac{\partial D}{\partial t}=& (\mathrm{1}-{\mathit{\delta }}^D){\mathit{\mu }}^{D}D-m^D\frac{D}{K_{\mathrm{m}}+D}D-\text{sh}\times w^D{D}^{\mathrm{2}}\\ & & -g^Z\left(D\right)Z-g^M\left(D\right)M\end{array}$$

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: $${\text{PAR}}^D={\mathit{\beta }}_{\mathrm{1}}^D{\text{PAR}}_{\mathrm{1}}+{\mathit{\beta }}_{\mathrm{2}}^D{\text{PAR}}_{\mathrm{2}}+{\mathit{\beta }}_{\mathrm{3}}^D{\text{PAR}}_{\mathrm{3}}.$$

The production terms for diatoms are defined as for nanophytoplankton, except that the limitation terms also include Si: $$\begin{array}{}& & L_{\text{lim}}^D=min\left(L_{{\mathrm{PO}}_{\mathrm{4}}}^D,,,L_{\mathrm{N}}^D,,,L_{\mathrm{Fe}}^D,,,L_{\mathrm{Si}}^D\right),& & L_{\mathrm{Si}}^D=\frac{\mathrm{Si}}{\mathrm{Si}+K_{\mathrm{Si}}^D}.\end{array}$$

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 $\mathrm{\mu }\mathrm{M}$, whereas values as high as 88.7 $\mathrm{\mu }\mathrm{M}$ have been measured for Antarctic species . 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 : $$K_{\mathrm{Si}}^D=K_{\mathrm{Si}}^{D,\text{min}}+\frac{\mathrm{7}{\stackrel{\mathrm{˘}}{\mathrm{Si}}}^{\mathrm{2}}}{(K_{\mathrm{Si}}{)}^{\mathrm{2}}+{\stackrel{\mathrm{˘}}{\mathrm{Si}}}^{\mathrm{2}}},$$ where $\stackrel{\mathrm{˘}}{\mathrm{Si}}$ here is the maximum Si concentration over a year (note that during the first year of a pluriannual simulation, $\stackrel{\mathrm{˘}}{\mathrm{Si}}$ is set to a constant). For the other nutrients, we use the same parameterization as for nanophytoplankton (see Eq. ).

The diatoms aggregation term $w_p^D$ 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 : $$w^D=w^P+w_{\text{max}}^D(\mathrm{1}-L_{\text{lim}}^D).$$

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^{\text{Chl}}$ (where $I$ denotes $P$ or $D$, typical units are $\mathrm{\mu }\mathrm{g}$ Chl ${\mathrm{L}}^{-\mathrm{1}}$ or $\mathrm{mg}\mathrm{Chl}{\mathrm{m}}^{-\mathrm{3}}$) for both phytoplankton groups is parameterized using the photo-adaptive model of : $$\begin{array}{cc}\frac{\partial I^{\text{Chl}}}{\partial t}=& (\mathrm{1}-{\mathit{\delta }}^I)(12{\mathit{\theta }}_{\text{min}}^{\text{Chl}}+({\mathit{\theta }}_{\text{max}}^{\text{Chl},I}-{\mathit{\theta }}_{\text{min}}^{\text{Chl}}\left){\mathit{\rho }}^{I^{\text{Chl}}}\right){\mathit{\mu }}^{I}I-m^I\frac{I}{K_{\mathrm{m}}+I}{I}^{\text{Chl}}\\ & & -\text{sh}\times w^{I}I{I}^{\text{Chl}}-{\mathit{\theta }}^{\text{Chl},I}{g}^Z\left(I\right)Z-{\mathit{\theta }}^{\text{Chl},I}{g}^M\left(I\right)M,\end{array}$$ where $I$ is the phytoplankton group and ${\mathit{\theta }}^{\text{Chl},I}$ is the chlorophyll-to-carbon ratio of the considered phytoplankton class; 12 represents the molar mass of carbon; ${\mathit{\rho }}^{I^{\text{Chl}}}$ the ratio of energy assimilated to energy absorbed as defined by : $$\begin{array}{}& & {\mathit{\rho }}^{I^{\text{Chl}}}=\frac{144{\stackrel{\mathrm{˘}}{\mathit{\mu }}}^{I}I}{{\mathit{\alpha }}^I{I}^{\text{Chl}}\frac{{\text{PAR}}^I}{L_{\text{day}}}},& & {\stackrel{\mathrm{˘}}{\mathit{\mu }}}^I={\mathit{\mu }}_P{f}_{\mathrm{2}}\left(z_{\text{mxl}}\right)\left(\mathrm{1},-,\exp ,\left(\frac{-{\mathit{\alpha }}^I{\mathit{\theta }}^{\text{Chl},I}{\text{PAR}}^I}{L_{\text{day}}{\mathit{\mu }}_P{L}_{\text{lim}}^I}\right)\right)L_{\text{lim}}^I.\end{array}$$

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 $\mathrm{mg}\mathrm{Chl}{\mathrm{m}}^{-\mathrm{3}}$. It should be noted that for chlorophyll synthesis, the second parameterization of phytoplankton growth is used to compute ${\stackrel{\mathrm{˘}}{\mathit{\mu }}}^I$ (see Eq. ). This is necessary because of the expression for ${\mathit{\rho }}_{\text{Chl}}^I$.

Iron in nanophytoplankton and diatoms

The temporal evolution of the iron biomass of phytoplankton $I^{\mathrm{Fe}}$ (model units are mol Fe ${\mathrm{L}}^{-\mathrm{1}}$), where $I$ denotes $P$ or $D$, is driven by the following equation $$\begin{array}{cc}\frac{\partial I^{\mathrm{Fe}}}{\partial t}=& (\mathrm{1}-{\mathit{\delta }}^I){\mathit{\mu }}^{I^{\mathrm{Fe}}}I-m^I\frac{I}{K_{\mathrm{m}}+I}{I}^{\mathrm{Fe}}-\text{sh}\times w^{I}I{I}^{\mathrm{Fe}}\\ & & -{\mathit{\theta }}^{\mathrm{Fe},I}{g}^Z\left(I\right)Z-{\mathit{\theta }}^{\mathrm{Fe},I}{g}^M\left(I\right)M.\end{array}$$

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^{\mathrm{Fe}}$) rather than the iron quota (${\mathit{\theta }}^{\mathrm{Fe},I}$) directly. Growth rate of the iron biomass of phytoplankton is parameterized according to $${\mathit{\mu }}^{I^{\mathrm{Fe}}}={\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},I}{L}_{\text{lim,1}}^{I^{\mathrm{Fe}}}{L}_{\text{lim,2}}^{I^{\mathrm{Fe}}}\frac{\mathrm{1}-\frac{{\mathit{\theta }}^{\mathrm{Fe},I}}{{\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},I}}}{1.05-\frac{{\mathit{\theta }}^{\mathrm{Fe},I}}{{\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},I}}}{\mathit{\mu }}_P.$$

As in , iron uptake is also downregulated via a feedback from ${\mathit{\theta }}^{\mathrm{Fe},I}$ using a normalized inverse hyperbolic function with a small shape factor set to 0.05.

In the former equation, $L_{\text{lim,1}}^{I^{\mathrm{Fe}}}$ is the iron limitation term and is modeled as follows: $$\begin{array}{}& & L_{\text{lim,1}}^{I^{\mathrm{Fe}}}=\frac{\mathrm{bFe}}{\mathrm{bFe}+K_{\mathrm{Fe}}^{I^{\mathrm{Fe}}}},& & K_{\mathrm{Fe}}^{I^{\mathrm{Fe}}}=K_{\mathrm{Fe}}^{{\mathrm{I}}^{\mathrm{Fe},\min }}\frac{I_{\mathrm{1}}+S_{\text{rat}}^I{I}_{\mathrm{2}}}{I_{\mathrm{1}}+I_{\mathrm{2}}},& & I_{\mathrm{2}}=max\left(\mathrm{0},,,I,-,I_{\text{max}}\right)\text{, }{I}_{\mathrm{1}}=I-I_{\mathrm{2}},\end{array}$$ where $\mathrm{bFe}$ is the concentration of bioavailable iron (see Sect. ). The half-saturation constant for iron uptake is also increasing with phytoplankton biomass as for the other half-saturation constants (see Eq. ).

At low iron concentrations, observations suggest that iron uptake might be enhanced, at least for some species , giving surge uptake. 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 . In PISCES, a different parameterization has been chosen since downregulation is already included in Eq. (): $$L_{\text{lim,2}}^{I^{\mathrm{Fe}}}=\frac{\mathrm{4}-4.5L_{\mathrm{Fe}}^I}{L_{\mathrm{Fe}}^I+0.5}.$$

$L_{\text{lim,2}}$ equals 4 at very low iron concentrations and 1 at high iron concentration. Overall, the downregulation in Eq. () 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 .

The demands for iron in phytoplankton are for photosynthesis, respiration and nitrate/nitrite reduction. Following , 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: $$\begin{array}{cc}& {\mathit{\theta }}_{\text{min}}^{\mathrm{Fe},I}=\frac{0.0016}{55.85}{\mathit{\theta }}^{\text{Chl},I}+\frac{1.21\times {10}^{-\mathrm{5}}\times 14}{55.85\times 7.625}{L}_{\mathrm{N}}^P\\ & & \times 1.5+\frac{1.15\times {10}^{-\mathrm{4}}\times 14}{55.85\times 7.625}{L}_{{\mathrm{NO}}_{\mathrm{3}}}^P.\end{array}$$

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 . The modeled iron quota in PISCES varies thus between this minimum quota ${\mathit{\theta }}_{\text{min}}^{\mathrm{Fe},I}$ and the maximum quota ${\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},I}$ , i.e., between about 1 and 40 $\mathrm{\mu }\mathrm{mol}\mathrm{Fe}(\mathrm{mol}\mathrm{C}{)}^{-\mathrm{1}}$ when using the standard set of parameters (see Table ).

Silicon in diatoms

$$\begin{array}{cc}\frac{\partial D^{\mathrm{Si}}}{\partial t}=& {\mathit{\theta }}_{\text{opt}}^{\mathrm{Si},D}(\mathrm{1}-{\mathit{\delta }}^D){\mathit{\mu }}^{D}D-{\mathit{\theta }}^{\mathrm{Si},D}{g}^M\left(D\right)M\\ & -{\mathit{\theta }}^{\mathrm{Si},D}{g}^Z\left(D\right)Z-m^D\frac{D}{K_{\mathrm{m}}+D}{D}^{\mathrm{Si}}\\ & & -\text{sh}\times w^{D}D{D}^{\mathrm{Si}}\end{array}$$

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\pm 0.13$ $\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}$ . Light, N, P, or Fe stress has been demonstrated to lead to heavier silicification e.g.,. 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 . 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): $$\begin{array}{cc}& {\mathit{\theta }}_{\text{opt}}^{\mathrm{Si},D}={\mathit{\theta }}_{\mathrm{m}}^{\mathrm{Si},D}{L}_{\text{lim,1}}^{D^{\mathrm{Si}}}\\ & & min\left(5.4,,,\left(4.4,\exp ,\left(-,4.23,F_{\text{lim,1}}^{D^{\mathrm{Si}}}\right),F_{\text{lim,2}}^{D^{\mathrm{Si}}},+,\mathrm{1}\right),\left(\mathrm{1},+,\mathrm{2},L_{\text{lim,2}}^{D^{\mathrm{Si}}}\right)\right).\end{array}$$

Relative to the original parameterization, an additional limitation term by Si has been added ($F_{\text{lim,2}}^{D^{\mathrm{Si}}}$) to produce a lighter silicification in case of Si exhaustion.

The different terms in Eq. () are defined as follows: $$\begin{array}{}& & F_{\text{lim,1}}^{D^{\mathrm{Si}}}=min\left(\frac{{\mathit{\mu }}^D}{{\mathit{\mu }}_P{L}_{\text{lim}}^D},,,L_{{\mathrm{PO}}_{\mathrm{4}}}^D,,,L_{\mathrm{N}}^D,,,L_{\mathrm{Fe}}^D\right),& & F_{\text{lim,2}}^{D^{\mathrm{Si}}}=min\left(\mathrm{1},,,2.2,max,\left(\mathrm{0},,,L_{\text{lim,1}}^{D^{\mathrm{Si}}},-,0.5\right)\right),& & L_{\text{lim,1}}^{D^{\mathrm{Si}}}=\frac{\mathrm{Si}}{\mathrm{Si}+K_{\mathrm{Si}}^{\mathrm{1}}},& & L_{\text{lim,2}}^{D^{\mathrm{Si}}}=\left\{\begin{array}{ll}\frac{{\mathrm{Si}}^{\mathrm{3}}}{{\mathrm{Si}}^{\mathrm{3}}+(K_{\mathrm{Si}}^{\mathrm{2}}{)}^{\mathrm{3}}}& \text{if}\mathit{\phi }<\mathrm{0}\\ \mathrm{0}& \text{if}\mathit{\phi }>\mathrm{0},\end{array}\right.\end{array}$$ where $\mathit{\phi }$ 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 . To reproduce those high Si $/$ C ratios, we have introduced the term $L_{\text{lim,2}}^{D^{\mathrm{Si}}}$ 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_{\mathrm{Si}}^{\mathrm{2}}$. 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.

Zooplankton

Microzooplankton

$$\begin{array}{cc}\frac{\partial Z}{\partial t}=& e^Z\left(g^Z,\left(,P,\right),+,g^Z,\left(,D,\right),+,g^Z,\left(,\text{POC},\right)\right)Z\\ & -g^M\left(Z\right)M-m^Z{f}_Z\left(T\right)Z^{\mathrm{2}}\\ & & -r^Z{f}_Z\left(T\right)\left(\frac{Z}{K_{\mathrm{m}}+Z},+,\mathrm{3},\mathrm{\Delta },\left(,{\mathrm{O}}_{\mathrm{2}},\right)\right)Z\end{array}$$

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: $$\begin{array}{}& & g_m^Z=g_{\text{max}}^{\mathrm{0},Z}{f}_Z\left(T\right),& & f_Z\left(T\right)=b_Z^T,\end{array}$$ where $g_{\text{max}}^{\mathrm{0},Z}$ is the maximum grazing rate at 0 ${}^{\circ }\mathrm{C}$, $b_Z$ is the temperature dependence and $T$ is the temperature. In their review, have found a $Q_{10}$ ($Q_{10}=b_Z^{10}$) 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 $$\begin{array}{cc}& F=\sum _J{p}_J^{Z}max\left(\mathrm{0},,,J,-,J_{\text{thresh}}^Z\right)\\ & F_{\text{lim}}=max\left(\mathrm{0},,,F,-,min,\left(0.5,F,,,F_{\text{thresh}}^Z\right)\right)\\ & & g^Z\left(I\right)=g_m^Z\frac{F_{\text{lim}}}{F}\frac{p_I^{Z}max\left(\mathrm{0},,,I,-,I_{\text{thresh}}^Z\right)}{K_{\mathrm{G}}^Z+{\sum }_J{p}_J^{Z}J},\end{array}$$ where $J$ denotes all the species microzooplankton can graze upon ($P$, $D$ and POC) and $p_J^Z$ is the preference microzooplankton has for each $J$. In PISCES, we have chosen a Michaelis–Menten parameterization with no switching and a threshold ($F_{\text{thresh}}^Z$) . This choice is rather arbitrary. Another very popular formulation in models is the Michaelis–Menten parameterization with active switching introduced by . However, this parameterization exhibits anomalous dynamics such as sub-optimal feeding . In our parameterization, a threshold for each individual resource ($J_{\text{thresh}}^Z$) can be specified in addition to the global threshold ($F_{\text{thresh}}^Z$). 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 . A complete review of the different expected responses has been presented by . For instance, when confronted with poor food quality, zooplankton can increase their ingestion rate , or decrease it as the food can become deleterious . 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 ${\mathit{\theta }}^{\mathrm{Fe},I}$ or the N $/$ C ratio ${\mathit{\theta }}^{\mathrm{N},I}$ of the preys decreases), $e^Z$ decreases: $$\begin{array}{}& & e_{\mathrm{N}}^Z=min\left(\mathrm{1},,,\frac{{\sum }_I{\mathit{\theta }}^{\mathrm{N},I}{g}^Z\left(I\right)}{{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\sum }_I{g}^Z\left(I\right)},,,\frac{{\sum }_I{\mathit{\theta }}^{\mathrm{Fe},I}{g}^Z\left(I\right)}{{\mathit{\theta }}^{\mathrm{Fe},Z}{\sum }_I{g}^Z\left(I\right)}\right),& & e^Z=e_{\mathrm{N}}^{Z}min\left(e_{\text{max}}^Z,,,(,\mathrm{1},-,{\mathit{\sigma }}^Z,),\frac{{\sum }_I{\mathit{\theta }}^{\mathrm{Fe},I}{g}^Z\left(I\right)}{{\mathit{\theta }}^{\mathrm{Fe},Z}{\sum }_I{g}^Z\left(I\right)}\right).\end{array}$$

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. b). 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 ${\mathit{\theta }}^{\mathrm{N},I}$).

Mesozooplankton

$$\begin{array}{cc}\frac{\partial M}{\partial t}=& e^M\left(g^M,\left(,P,\right),+,g^M,\left(,D,\right),+,g^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{GOC},\right)\right.\\ & \left..+,g_{\text{FF}}^M,\left(,\text{POC},\right),+,g^M,\left(,Z,\right)\right)M\\ & -m^M{f}_M\left(T\right)M^{\mathrm{2}}-r^M{f}_M\left(T\right)\\ & & \left(\frac{M}{K_{\mathrm{m}}+M},+,\mathrm{3},\mathrm{\Delta },\left(,{\mathrm{O}}_{\mathrm{2}},\right)\right)M\end{array}$$

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 .

Parameterization of mesozooplankton grazing is similar to microzooplankton. In addition to the “conventional” concentration-dependent grazing described by Eq. (), 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 . 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: $$\begin{array}{}& & g_{\text{FF}}^M\left(\text{POC}\right)=g_{\text{FF}}{f}_M\left(T\right)w_{\text{POC}}\text{POC},& & g_{\text{FF}}^M\left(\text{GOC}\right)=g_{\text{FF}}{f}_M\left(T\right)w_{\text{GOC}}\text{GOC}.\end{array}$$

This importance of flux feeding has been analyzed in PISCES by . 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. (), the term with a quadratic dependency to mesozooplankton does not depict aggregation but grazing by the higher, non-resolved trophic levels. Following , 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: $$\begin{array}{}& & P_{\text{up}}^M={\mathit{\sigma }}^M{f}_{\text{up}}\left(e_{\text{max}}^M\right)m^M{f}_M\left(T\right)M^{\mathrm{2}},& & R_{\text{up}}^M=(\mathrm{1}-{\mathit{\sigma }}^M-e_{\text{max}}^M)f_{\text{up}}\left(e_{\text{max}}^M\right)m^M{f}_M\left(T\right)M^{\mathrm{2}},\end{array}$$ where function $f_{\text{up}}\left(x\right)$ is $$f_{\text{up}}\left(x\right)=\sum _{i=\mathrm{0}}^{\mathrm{\infty }}{x}^i=\frac{\mathrm{1}}{\mathrm{1}-x}\text{for}\mathrm{0}<x<1.$$

It should be noted here that a similar quadratic term is also included in the equation for microzooplankton (see Eq. ) 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 $$\begin{array}{cc}\frac{\partial \text{DOC}}{\partial t}=& (\mathrm{1}-{\mathit{\gamma }}^Z)(\mathrm{1}-e^Z-{\mathit{\sigma }}^Z)\sum _I{g}^Z\left(I\right)Z\\ & +(\mathrm{1}-{\mathit{\gamma }}^M)(\mathrm{1}-e^M-{\mathit{\sigma }}^M)\\ & \left(\sum _I,g^Z,\left(,I,\right),+,g_{\text{FF}}^M,\left(,\text{GOC},\right)\right)M+{\mathit{\delta }}^D{\mathit{\mu }}^{D}D\\ & +{\mathit{\delta }}^P{\mathit{\mu }}^{P}P+{\mathit{\lambda }}_{\text{POC}}^{\star }\text{POC}\\ & +(\mathrm{1}-{\mathit{\gamma }}^M)R_{\text{up}}^M-\text{Remin}-\text{Denit}\\ & & -{\mathrm{\Phi }}_{\mathrm{1}}^{\text{DOC}}-{\mathrm{\Phi }}_{\mathrm{2}}^{\text{DOC}}-{\mathrm{\Phi }}_{\mathrm{3}}^{\text{DOC}},\end{array}$$ where $I$ includes $P$, $D$ and POC for microzooplankton and $P$, $D$, $Z$ and POC for mesozooplankton (see Eqs.  and , 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 . Transport of this pool of dissolved organic matter can make a significant part of the carbon pump . 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: $$\begin{array}{}& & \mathrm{\text{Remin}}=min\left(\frac{{\mathrm{O}}_{\mathrm{2}}}{{\mathrm{O}}_{\mathrm{2}}^{\text{ut}}},,,{\mathit{\lambda }}_{\text{DOC}},f_P,\left(,T,\right),\left(\mathrm{1},-,\mathrm{\Delta },\left(,{\mathrm{O}}_{\mathrm{2}},\right)\right),L_{\text{lim}}^{\text{bact}},\frac{\text{Bact}}{{\text{Bact}}_{\text{ref}}},\text{DOC}\right),& & \mathrm{\text{Denit}}=min\left(\frac{{\mathrm{NO}}_{\mathrm{3}}}{r_{{\mathrm{NO}}_{\mathrm{3}}}^{\star }},,,{\mathit{\lambda }}_{\text{DOC}},f_P,\left(,T,\right),\mathrm{\Delta },\left(,{\mathrm{O}}_{\mathrm{2}},\right),L_{\text{lim}}^{\text{bact}},\frac{\text{Bact}}{{\text{Bact}}_{\text{ref}}},\text{DOC}\right).\end{array}$$

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 $\mathrm{\Delta }\left({\mathrm{O}}_{\mathrm{2}}\right)$ that varies between 0 and 1 (see Sect.  for the formulation of this factor). It is assumed that the specific rates of degradation (${\mathit{\lambda }}_{\text{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.,, 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^{\text{bact}}$) according to $$\begin{array}{}& & L^{\text{bact}}=L_{\text{lim}}^{\text{bact}}{L}_{\text{DOC}}^{\text{bact}},& & L_{\text{DOC}}^{\text{bact}}=\frac{\text{DOC}}{\text{DOC}+K_{\text{DOC}}},& & L_{\text{lim}}^{\text{bact}}=min\left(L_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}},,,L_{{\mathrm{PO}}_{\mathrm{4}}}^{\text{bact}},,,L_{\mathrm{Fe}}^{\text{bact}}\right),& & L_{\mathrm{Fe}}^{\text{bact}}=\frac{\mathrm{bFe}}{\mathrm{bFe}+K_{\mathrm{Fe}}^{\text{bact}}},& & L_{{\mathrm{PO}}_{\mathrm{4}}}^{\text{bact}}=\frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K_{{\mathrm{PO}}_{\mathrm{4}}}^{\text{bact}}},& & L_{\mathrm{N}}^{\text{bact}}=L_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}+L_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}},& & L_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}}=\frac{K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}{\mathrm{NH}}_{\mathrm{4}}}{K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}{K}_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}}+K_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}}{\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}{\mathrm{NH}}_{\mathrm{4}}},& & L_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}=\frac{K_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}}{\mathrm{NO}}_{\mathrm{3}}}{K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}{K}_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}}+K_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{bact}}{\mathrm{NO}}_{\mathrm{3}}+K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{bact}}{\mathrm{NH}}_{\mathrm{4}}}.\end{array}$$

The half-saturation constants of the P and N limitation terms ($K_i^{\text{bact}}$) are set in the namelist.

In PISCES, bacterial biomass is not explicitly modeled; Instead, we use the following formulation: $$\begin{array}{}& & z_{\text{max}}=max\left(z_{\text{mxl}},,,z_{\text{eu}}\right),& & \text{Bact}=\left\{\begin{array}{ll}min\left(0.7,(,Z,+,\mathrm{2},M\right),\mathrm{4}\mathrm{\mu }\mathrm{mol}\mathrm{C}{\mathrm{L}}^{-\mathrm{1}})& \text{if}z\le z_{\text{max}}\\ \text{Bact}\left(z_{\text{max}}\right){\left(\frac{z_{\text{max}}}{z}\right)}^{0.683}& \text{Otherwise}.\end{array}\right.\end{array}$$

In the previous equation, $0.7(Z+\mathrm{2}M)$ is a proxy for the bacterial concentration. This relationship has been constructed from an unpublished version of PISCES (already mentioned in ) that includes an explicit description of the bacterial biomass. Below a certain depth ($z_{\text{max}}$), this biomass decreases with depth via a power-law function .

In Eq. (), the terms ${\mathrm{\Phi }}^{\text{DOC}}$ denote aggregation processes and are described hereafter (see Sect. ). For DOM, we consider turbulence-induced as well as Brownian aggregation processes. $$\begin{array}{}& & {\mathrm{\Phi }}_{\mathrm{1}}^{\text{DOC}}=\text{sh}\times (a_{\mathrm{1}}\text{DOC}+a_{\mathrm{2}}\text{POC})\text{DOC}& & {\mathrm{\Phi }}_{\mathrm{2}}^{\text{DOC}}=\text{sh}\times a_{\mathrm{3}}\text{GOC}\times \text{DOC}& & {\mathrm{\Phi }}_{\mathrm{3}}^{\text{DOC}}=(a_{\mathrm{4}}\text{POC}+a_{\mathrm{5}}\text{DOC})\text{DOC}\end{array}$$

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 $\mathrm{\mu }\mathrm{m}$) and GOC for the larger class (100–5000 $\mathrm{\mu }\mathrm{m}$).

• A more complex model proposed by 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: $$\begin{array}{cc}\frac{\partial \text{POC}}{\partial t}=& {\mathit{\sigma }}^Z\sum g^Z\left(X\right)Z+0.5m^D\frac{D}{D+K_{\mathrm{m}}}D\\ & +r^Z{f}_Z\left(T\right)\frac{Z}{Z+K_{\mathrm{m}}}Z+m^Z{f}_Z\left(T\right)Z^{\mathrm{2}}\\ & +(\mathrm{1}-0.5R_{{\mathrm{CaCO}}_{\mathrm{3}}})(m^P\frac{P}{P+K_{\mathrm{m}}}P+w^P{P}^{\mathrm{2}})\\ & +{\mathit{\lambda }}_{\text{POC}}^{\star }\text{GOC}+{\mathrm{\Phi }}_{\mathrm{1}}^{\text{DOC}}+{\mathrm{\Phi }}_{\mathrm{3}}^{\text{DOC}}\\ & -\left(g^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right)\right)M-g^Z\left(\text{POC}\right)Z\\ & & -{\mathit{\lambda }}_{\text{POC}}^{\star }\text{POC}-\mathrm{\Phi }-w_{\text{POC}}\frac{\partial \text{POC}}{\partial z},\end{array}$$ where $w_{\text{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_{{\mathrm{CaCO}}_{\mathrm{3}}}$). We assume that 50 $\mathrm{\%}$ of the organic matter of the calcifiers is associated with the shell. Since calcite is significantly denser than organic matter, 50 $\mathrm{\%}$ 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 ${\mathit{\lambda }}_{\text{POC}}^{\star }$ 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 . In , the attenuation coefficient ($b$) for the flux was found to be about 0.4 instead of the standard value 0.86 . This corresponds to a 45 $\mathrm{\%}$ decrease of the degradation rate in anoxic waters relative to oxic waters, which is implemented as $${\mathit{\lambda }}_{\text{POC}}^{\star }={\mathit{\lambda }}_{\text{POC}}{f}_P\left(T\right)\left(\mathrm{1},-,0.45,\mathrm{\Delta },\left(,{\mathrm{O}}_{\mathrm{2}},\right)\right).$$

POC experiences aggregation due to turbulence and differential settling: $$\begin{array}{cc}& \mathrm{\Phi }=\text{sh}\times a_{\mathrm{6}}{\text{POC}}^{\mathrm{2}}+\text{sh}\times a_{\mathrm{7}}\text{POC}\times \text{GOC}+a_{\mathrm{8}}\text{POC}\\ & & \times \text{GOC}+a_{\mathrm{9}}{\text{POC}}^{\mathrm{2}}.\end{array}$$

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 steady-state power-law size spectrum for particles with an exponent of 3.6. Subsequently, the different coagulation kernels e.g., 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 $$\begin{array}{cc}\frac{\partial \text{GOC}}{\partial t}=& {\mathit{\sigma }}^M\left(\sum _I,g^M,\left(,I,\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{GOC},\right)\right)\\ & M+r^M{f}_M\left(T\right)\frac{M}{M+K_{\mathrm{m}}}M+P_{\text{up}}^M\\ & +0.5R_{{\mathrm{CaCO}}_{\mathrm{3}}}\left(m^P,\frac{P}{P+K_{\mathrm{m}}},P,+,w^P,P^{\mathrm{2}}\right)\\ & +0.5m^D\frac{D}{D+K_{\mathrm{m}}}D{w}^D{D}^{\mathrm{2}}\\ & +\mathrm{\Phi }+{\mathrm{\Phi }}_{\mathrm{2}}^{\text{DOC}}-g_{\text{FF}}^M\left(\text{GOC}\right)M-{\mathit{\lambda }}_{\text{POC}}^{\star }\text{GOC}\\ & & -w_{\text{GOC}}\frac{\partial \text{GOC}}{\partial z}.\end{array}$$

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.,. Such an increase is consistent with the power-law formulation proposed by . Such an increase in the settling speed is parameterized in PISCES for GOC as follows: $$\begin{array}{}& & z_{\text{max}}=max\left(z_{\text{eu}},,,z_{\text{mxl}}\right),& & w_{\text{GOC}}=w_{\text{GOC}}^{\text{min}}+(200-w_{\text{GOC}}^{\text{min}})\frac{max\left(\mathrm{0},,,z,-,z_{\text{max}}\right)}{5000}.\end{array}$$

The parameters in this equation have been adjusted using a model of aggregation/disaggregation with multiple size classes . The maximum sinking speed is set to 200 $\mathrm{m}{\mathrm{d}}^{-\mathrm{1}}$ and is reached at about 5000 m depth over most of the ocean since $z_{\text{max}}$ is generally less than 100 m. We have not included any ballasting effect due to the higher density of biogenic silica or calcite . In fact, observations are rather contradictory on this ballast effect . 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 .

Kriest model of particulate organic matter

Here we present a brief overview of the model of . The reader is referred to the literature, where the method has been presented e.g.,, for more detail. The model postulates that the carbon content ($m\left(d_i\right)$), the sinking speed ($w\left(d_i\right)$) and the abundance of the aggregates ($n\left(d_i\right)$) can be described by power-law functions of their diameters ($d_i$): $$\begin{array}{}& & m\left(d_i\right)=C{d}_i^{\mathit{\zeta }},& & w\left(d_i\right)=B{d}_i^{\mathit{\nu }},& & n\left(d_i\right)=A{d}_i^{\mathit{ϵ}}.\end{array}$$

It is also assumed, as in , 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: $$\mathit{ϵ}=\frac{(\mathit{\zeta }+\mathrm{1})\text{POC}-m_l\text{NUM}}{\text{POC}-m_l\text{NUM}},$$ where $m_l$ is the mass of the smallest aggregate (of size $l$).

Having $\mathit{ϵ}$, the average sinking speed of numbers ($w_{\text{NUM}}$) and mass ($w_{\text{POC}}$) can be computed following : $$\begin{array}{}& & w_{\text{POC}}=w_l\frac{\mathit{\zeta }+\mathrm{1}-\mathit{ϵ}+(\frac{L}{l}{)}^{\mathrm{1}+\mathit{\nu }+\mathit{\zeta }-\mathit{ϵ}}\mathit{\nu }}{\mathrm{1}+\mathit{\nu }+\mathit{\zeta }-\mathit{ϵ}},& & w_{\text{NUM}}=w_l\frac{\mathrm{1}-\mathit{ϵ}+(\frac{L}{l}{)}^{\mathrm{1}+\mathit{\nu }-\mathit{ϵ}}\mathit{\nu }}{\mathrm{1}+\mathit{\nu }-\mathit{ϵ}}.\end{array}$$

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 $\mathit{ϵ}$. As a result, the sinking speeds for both mass and number vary with space and time.

Aggregation ($\mathit{\xi }$) 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 ; see for the term $\mathit{\xi }$ and its computation. Currently it is assumed that turbulent shear rate is high in the mixed layer (1 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$), and low below (0.01 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). Summing up the number of collisions due to turbulent shear and differential settlement, $C_{\text{sh}}$ and $C_{\text{ds}}$, respectively, the decrease of the number of particles due to aggregation is then: $$\mathit{\xi }=\text{Stick}\times (C_{\text{sh}}+C_{\text{ds}}).$$

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 $$\begin{array}{cc}\frac{\partial \text{POC}}{\partial t}=& {\mathit{\sigma }}^Z\sum g^Z\left(X\right)Z+{\mathit{\sigma }}^M\left(\sum _I,g^M,\left(,I,\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right)\right)M\\ & +m^P\frac{P}{P+K_{\mathrm{m}}}P\\ & +w^P{P}^{\mathrm{2}}+m^D\frac{D}{D+K_{\mathrm{m}}}D+w^D{D}^{\mathrm{2}}\\ & +r^Z{f}_Z\left(T\right)\frac{Z}{Z+K_{\mathrm{m}}}Z+m^Z{f}_Z\left(T\right)Z^{\mathrm{2}}\\ & +r^M{f}_M\left(T\right)\frac{M}{M+K_{\mathrm{m}}}M+P_{\text{up}}^M-g^M\left(\text{POC}\right)M\\ & -g^Z\left(\text{POC}\right)Z-{\mathit{\lambda }}_{\text{POC}}^{\star }\text{POC}\\ & +{\mathrm{\Phi }}_{\mathrm{1}}^{\text{DOC}}+{\mathrm{\Phi }}_{\mathrm{3}}^{\text{DOC}}-\left(g^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right)\right)M\\ & -g^Z\left(\text{POC}\right)Z\\ & & -{\mathit{\lambda }}_{\text{POC}}^{\star }\text{POC}-w_{\text{POC}}\frac{\partial \text{POC}}{\partial z}.\end{array}$$

This is exactly equal to the sum of the two equations used for the temporal evolution of POC and GOC in the two-compartment model of PISCES (see Eqs.  and ). $$\begin{array}{cc}\frac{\partial \text{NUM}}{\partial t}=& \frac{{\mathit{\sigma }}^Z\sum g^Z\left(X\right)Z}{{\overline{m}}_Z}+\frac{{\mathit{\sigma }}^M\left({\sum }_I,g^M,\left(,I,\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right)\right)M}{{\overline{m}}_M}\\ & +\frac{m^P\frac{P}{P+K_{\mathrm{m}}}P+w^P{P}^{\mathrm{2}}}{{\overline{m}}_P}\\ & +\frac{m^D\frac{D}{D+K_{\mathrm{m}}}D+w^D{D}^{\mathrm{2}}}{{\overline{m}}_D}\\ & +\frac{r^Z{f}_Z\left(T\right)\frac{Z}{Z+K_{\mathrm{m}}}Z+m^Z{f}_Z\left(T\right)Z^{\mathrm{2}}}{{\overline{m}}_Z}\\ & +\frac{r^M{f}_M\left(T\right)\frac{M}{M+K_{\mathrm{m}}}M+P_{\text{up}}^M}{{\overline{m}}_M}\\ & -\frac{\left(g^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right)\right)M}{{\overline{m}}_M}\\ & -\frac{g^Z\left(\text{POC}\right)Z}{{\overline{m}}_Z}-{\mathit{\lambda }}_{\text{POC}}^{\star }\text{POC}+\frac{{\mathrm{\Phi }}_{\mathrm{1}}^{\text{DOC}}+{\mathrm{\Phi }}_{\mathrm{3}}^{\text{DOC}}}{m_l}\\ & & -\mathit{\xi }-w_{\text{NUM}}\frac{\partial \text{NUM}}{\partial z}\end{array}$$

In this equation, each process affecting the mass of the particles is divided by the mean mass ($\overline{m}$) of the compartment exerting this process to convert to numbers.

Iron in particles

In this subsection, the description corresponds to the two-compartment version of the model. To obtain the Kriest version, the equations for both $\mathrm{SFe}$ and $\mathrm{BFe}$, the iron content of the small and big particles, respectively, should be simply summed. $$\begin{array}{cc}\frac{\partial \mathrm{SFe}}{\partial t}=& {\mathit{\sigma }}^Z\sum _I{\mathit{\theta }}^{\mathrm{Fe},I}{g}^Z\left(I\right)Z+{\mathit{\theta }}^{\mathrm{Fe},Z}\left(r^Z{f}_Z\right(T)\frac{Z}{Z+K_{\mathrm{m}}}Z\\ & +m^Z{f}_Z\left(T\right)Z^{\mathrm{2}})\\ & +{\mathit{\lambda }}_{\text{GOC}}^{\star }\mathrm{BFe}+{\mathit{\theta }}^{\mathrm{Fe},P}(\mathrm{1}-0.5R_{{\mathrm{CaCO}}_{\mathrm{3}}})\\ & (m^P\frac{P}{P+K_{\mathrm{m}}}P+\text{sh}\times w^P{P}^{\mathrm{2}})\\ & +{\mathit{\theta }}^{\mathrm{Fe},\mathrm{D}}0.5m^D\frac{D}{D+K_{\mathrm{m}}}D+{\mathit{\lambda }}_{\mathrm{Fe}}\text{POC}{\mathrm{Fe}}^{\prime }\\ & +\text{Cgfe}\mathrm{1}-{\mathit{\lambda }}_{\text{POC}}^{\star }\mathrm{SFe}-{\mathit{\theta }}^{\mathrm{Fe},\text{POC}}\mathrm{\Phi }\\ & -{\mathit{\theta }}^{\mathrm{Fe},\text{POC}}\left(g^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right)\right)M+\\ & {\mathit{\kappa }}_{\text{Bact}}^{\mathrm{SFe}}\text{Bactfe}-{\mathit{\theta }}^{\mathrm{Fe},\text{POC}}{g}^Z\left(\text{POC}\right)\\ & & -w_{\text{POC}}\frac{\partial \mathrm{SFe}}{\partial z},\frac{\partial \mathrm{BFe}}{\partial t}=& {\mathit{\sigma }}^M\left(\sum _I,{\mathit{\theta }}^{\mathrm{Fe},I},g^M,\left(,I,\right)\right.\\ & \left..+,{\mathit{\theta }}^{\mathrm{Fe},\text{POC}},g_{\text{FF}}^M,\left(,\text{POC},\right),+,{\mathit{\theta }}^{\mathrm{Fe},\text{GOC}},g_{\text{FF}}^M,\left(,\text{GOC},\right)\right)M\\ & +{\mathit{\theta }}^{\mathrm{Fe},\mathrm{M}}\left(r^M{f}_M\right(T)\frac{M}{M+K_{\mathrm{m}}}M+P_{\text{up}}^M)+\\ & {\mathit{\theta }}^{\mathrm{Fe},P}0.5R_{{\mathrm{CaCO}}_{\mathrm{3}}}(m^P\frac{P}{P+K_{\mathrm{m}}}P\\ & +\text{sh}\times w^P{P}^{\mathrm{2}})\\ & +{\mathit{\theta }}^{\mathrm{Fe},\mathrm{D}}(0.5m^D\frac{D}{D+K_{\mathrm{m}}}D\\ & +\text{sh}\times w^D{D}^{\mathrm{2}})+{\mathit{\kappa }}_{\text{Bact}}^{\mathrm{BFe}}\text{Bactfe}\\ & +{\mathit{\lambda }}_{\mathrm{Fe}}\text{GOC}{\mathrm{Fe}}^{\prime }+{\mathit{\theta }}^{\mathrm{Fe},\text{POC}}\mathrm{\Phi }+\text{Cgfe}\mathrm{2}\\ & -{\mathit{\theta }}^{\mathrm{Fe},\text{GOC}}{g}_{\text{FF}}^M\left(\text{GOC}\right)M-{\mathit{\lambda }}_{\text{POC}}^{\star }\mathrm{BFe}\\ & & -w_{\text{GOC}}\frac{\partial \mathrm{BFe}}{\partial z},\end{array}$$ where Fe${}^{\prime }$ is the free form of dissolved iron. Its determination is detailed in Sect. . Bactfe is the amount of iron taken up by bacteria which is lost as particulate organic iron. Its computation is detailed in Sect. .

The free form of dissolved iron Fe${}^{\prime }$ 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.,: $$\begin{array}{cc}& {\mathit{\lambda }}_{\mathrm{Fe}}^{\star }={\mathit{\lambda }}_{\mathrm{Fe}}^{\text{min}}+{\mathit{\lambda }}_{\mathrm{Fe}}(\text{POC}+\text{GOC}+{\mathrm{CaCO}}_{\mathrm{3}}+\mathrm{BSi})\\ & & +{\mathit{\lambda }}_{\mathrm{Fe}}^{\text{dust}}\text{Dust},& & \text{Scav}={\mathit{\lambda }}_{\mathrm{Fe}}^{\star }{\mathrm{Fe}}^{\prime }.\end{array}$$

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 . The concentration of lithogenic particles is estimated as described in Eq. (). Model estimates suggest a different affinity for these particles compared to biogenic particles, which justifies the split between biogenic and lithogenic materials in Eq. (). The amount of iron that is scavenged by POC (${\mathit{\lambda }}_{\mathrm{Fe}}\text{POC}{\mathrm{Fe}}^{\prime }$) and GOC (${\mathit{\lambda }}_{\mathrm{Fe}}\text{GOC}{\mathrm{Fe}}^{\prime }$) is then allocated to $\mathrm{SFe}$ and $\mathrm{BFe}$, respectively.

PSi

$$\begin{array}{cc}\frac{\partial \mathrm{PSi}}{\partial t}=& {\mathit{\theta }}^{\mathrm{Si},D}{g}^M\left(D\right)M\\ & +{\mathit{\theta }}^{\mathrm{Si},D}{g}^Z\left(D\right)Z+{\mathit{\theta }}^{\mathrm{Si},D}{m}^D\frac{D}{K_{\mathrm{m}}+D}{D}^{\mathrm{Si}}\\ & +\text{sh}\times w^{D}D{D}^{\mathrm{Si}}\\ & & -{\mathit{\lambda }}_{\mathrm{PSi}}^{\star }{\text{Diss}}_{\mathrm{Si}}\mathrm{PSi}-w_{\text{GOC}}\frac{\partial {\mathrm{CaCO}}_{\mathrm{3}}}{\partial z}\end{array}$$

The dissolution rate of $\mathrm{PSi}$ depends on in situ temperature and on silicic acid saturation following the parameterization proposed by : $$\begin{array}{cc}& {\mathrm{Si}}_{\text{eq}}={10}^{6.44-\frac{968}{T+273.15}}\\ & {\mathrm{Si}}_{\text{sat}}=\frac{{\mathrm{Si}}_{\text{eq}}-\mathrm{Si}}{{\mathrm{Si}}_{\text{eq}}}\\ & {\mathit{\lambda }}_{\mathrm{PSi}}^{\star }={\mathit{\lambda }}_{\mathrm{PSi}}\left[0.225,\left(\mathrm{1},+,\frac{T}{15}\right),{\mathrm{Si}}_{\text{sat}}\right.\\ & & +0.775{\left({\left(\mathrm{1},+,\frac{T}{400}\right)}^{\mathrm{4}},{\mathrm{Si}}_{\text{sat}}\right)}^{\mathrm{9}}.\end{array}$$

The evolution of ${\mathit{\lambda }}_{\mathrm{PSi}}^{\star }$ as a function of Si and of temperature is shown on Fig. .

Laboratory experiments show that the diatom frustule is made of two biogenic silica phases which dissolve simultaneously, but at different rates e.g.,. 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 . 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 . In PISCES, despite this uncertainty, we model silica dissolution using two phases. The proportion of the most “labile” phase is set to a constant (${\mathit{\chi }}_{\text{lab}}$) in the upper ocean and is computed in the rest of the ocean assuming steady state: $$\begin{array}{cc}& z_{\text{max}}=max\left(z_{\text{eu}},,,z_{\text{mxl}}\right)\\ & & {\mathit{\chi }}_{\text{lab}}=\left\{\begin{array}{ll}{\mathit{\chi }}_{\text{lab}}^{\mathrm{0}}& \text{if}z\le z_{\text{max}}\\ {\mathit{\chi }}_{\text{lab}}^{\mathrm{0}}\exp \left(-,\left({\mathit{\lambda }}_{\mathrm{PSi}}^{\text{lab}},-,{\mathit{\lambda }}_{\mathrm{PSi}}^{\text{ref}}\right),\left(\frac{z-z_{\text{max}}}{w_{\text{GOC}}}\right)\right)& \text{Otherwise},\end{array}\right.& & {\mathit{\lambda }}_{\mathrm{PSi}}={\mathit{\chi }}_{\text{lab}}{\mathit{\lambda }}_{\mathrm{PSi}}^{\text{lab}}+(\mathrm{1}-{\mathit{\chi }}_{\text{lab}}){\mathit{\lambda }}_{\mathrm{PSi}}^{\text{ref}}.\end{array}$$

Nutrients

Nitrate and ammonium

$$\begin{array}{cc}& \frac{\partial {\mathrm{NO}}_{\mathrm{3}}}{\partial t}=\text{Nitrif}-{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^{P}P-{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^{D}D\\ & & -R_{{\mathrm{NH}}_{\mathrm{4}}}{\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}\mathrm{\Delta }\left({\mathrm{O}}_{\mathrm{2}}\right){\mathrm{NH}}_{\mathrm{4}}-R_{{\mathrm{NO}}_{\mathrm{3}}}\text{Denit}& \frac{\partial {\mathrm{NH}}_{\mathrm{4}}}{\partial t}={\mathit{\gamma }}^Z(\mathrm{1}-e^Z-{\mathit{\sigma }}^Z)\sum _I{g}^Z\left(I\right)Z+{\mathit{\gamma }}^M(\mathrm{1}-e^M-{\mathit{\sigma }}^M)\\ & \times \left(\sum _I,g^M,\left(,I,\right),+,g_{\text{FF}}^M,\left(,\text{POC},\right),+,g_{\text{FF}}^M,\left(,\text{GOC},\right)\right)M\\ & +{\mathit{\gamma }}^M{R}_{\text{up}}^M+\text{Remin}+\text{Denit}+N_{\text{fix}}\\ & -\text{Nitrif}-{\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}\mathrm{\Delta }\left({\mathrm{O}}_{\mathrm{2}}\right){\mathrm{NH}}_{\mathrm{4}}\\ & & -{\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^{P}P-{\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^{D}D\end{array}$$

Nitrification (Nitrif) corresponds to the conversion of ammonium to nitrate due to bacterial activity. It is assumed to be photo-inhibited e.g., and reduced in suboxic waters: $$\begin{array}{}& & \text{Nitrif}={\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}\frac{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{1}+\stackrel{\mathrm{‾}}{\text{PAR}}}\left(\mathrm{1},-,\mathrm{\Delta },\left(,{\mathrm{O}}_{\mathrm{2}},\right)\right),& & \stackrel{\mathrm{‾}}{\text{PAR}}=\stackrel{\mathrm{‾}}{{\text{PAR}}_{\mathrm{1}}+{\text{PAR}}_{\mathrm{2}}+{\text{PAR}}_{\mathrm{3}}},\end{array}$$ where $\stackrel{\mathrm{‾}}{\text{PAR}}$ is the PAR averaged over the mixed layer and $\mathrm{\Delta }\left({\mathrm{O}}_{\mathrm{2}}\right)$ varies between 0 (oxic conditions, ${\mathrm{O}}_{\mathrm{2}}>{{\mathrm{O}}_{\mathrm{2}}}^{\text{min,1}}$) and 1 (anoxia) according to $$\mathrm{\Delta }\left({\mathrm{O}}_{\mathrm{2}}\right)=min\left(\mathrm{1},,,max,\left(\mathrm{0},,,0.4,\frac{{{\mathrm{O}}_{\mathrm{2}}}^{\text{min,1}}-{\mathrm{O}}_{\mathrm{2}}}{{{\mathrm{O}}_{\mathrm{2}}}^{\text{min,2}}+{\mathrm{O}}_{\mathrm{2}}}\right)\right).$$

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_{{\mathrm{NO}}_{\mathrm{3}}}$ can be computed from $R_{-{\mathrm{O}}_{\mathrm{2}}/{\mathrm{NO}}_{\mathrm{3}}}$ and is found to be 0.86 . Equation (), implies that denitrification stops at oxygen concentration above 6 $\mathrm{\mu }\mathrm{M}$ . 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 (${\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}$).

Finally, nitrogen fixation is parameterized in PISCES as follows: $$\begin{array}{cc}& & L_N^{\text{Dz}}=\left\{\begin{array}{ll}0.01& \text{if}{L}_N^P\ge 0.8\\ \mathrm{1}-L_N^P& \text{Otherwise},\end{array}\right.& N_{\text{fix}}=N_{\text{fix}}^{m}max\left(\mathrm{0},,,{\mathit{\mu }}_P,-,2.15\right)L_N^{\text{Dz}}\\ & & min\left(\frac{\mathrm{bFe}}{K_{\mathrm{Fe}}^{\text{Dz}}+\mathrm{bFe}},,,\frac{{\mathrm{PO}}_{\mathrm{4}}}{K_{{\mathrm{PO}}_{\mathrm{4}}}^{P,\text{min}}+{\mathrm{PO}}_{\mathrm{4}}}\right)\left(\mathrm{1},-,e^{\frac{-\text{PAR}}{E_{\text{fix}}}}\right).\end{array}$$

This very crude parameterization is based on the following assumptions that have been inferred from studies of Trichodesmium e.g.,:

• Nitrogen fixation is restricted to warm waters above 20 ${}^{\circ }\mathrm{C}$ (${\mathit{\mu }}_P>{\mathit{\mu }}_P\left(20\right)=2.15$).

• Nitrogen fixation is restricted to areas with insufficient nitrogen ($L_N^P<0.8$).

• Nitrogen fixation requires iron and phosphorus.

• Nitrogen fixation needs high light levels, i.e., $E_{\text{fix}}$ is high.

The scaling factor $N_{\text{fix}}^m$ is set from the namelist and thus, may be chosen by the user.

Phosphate

$$\begin{array}{cc}\frac{\partial {\text{PO}}_{\mathrm{4}}}{\partial t}=& {\mathit{\gamma }}^Z(\mathrm{1}-e^Z-{\mathit{\sigma }}^Z)\sum _I{g}^Z\left(I\right)Z\\ & +{\mathit{\gamma }}^M(\mathrm{1}-e^M-{\mathit{\sigma }}^M)\left(\sum _I,g^M,\left(,I,\right),+,\sum _I,g_{\text{FF}}^M,\left(,I,\right)\right)M\\ & & +{\mathit{\gamma }}^M{R}_{\text{up}}^M+\text{Remin}+\text{Denit}-{\mathit{\mu }}^{P}P-{\mathit{\mu }}^{D}D\end{array}$$

All terms in this equation have been described previously.

Iron

$$\begin{array}{cc}\frac{\partial \mathrm{Fe}}{\partial t}=& max\left(\mathrm{0},,,(,\mathrm{1},-,{\mathit{\sigma }}^Z,),\frac{{\sum }_I{\mathit{\theta }}^{\mathrm{Fe},I}{g}^Z\left(I\right)}{{\sum }_I{g}^Z\left(I\right)},-,e_{\mathrm{N}}^Z,{\mathit{\theta }}^{\mathrm{Fe},Z}\right)\sum _I{g}^Z\left(I\right)Z\\ & +max\left(\mathrm{0},,,(,\mathrm{1},-,{\mathit{\sigma }}^M,),\frac{{\sum }_I{\mathit{\theta }}^{\mathrm{Fe},I}{g}^M\left(I\right)+{\sum }_I{\mathit{\theta }}^{\mathrm{Fe},I}{g}_{\text{FF}}^M\left(I\right)}{{\sum }_I{g}^M\left(I\right)+{\sum }_I{g}_{\text{FF}}^M\left(I\right)}\right.\\ & \left..-,e_{\mathrm{N}}^M,{\mathit{\theta }}^{\mathrm{Fe},Z}\right)\\ & \left(\sum _I,g^M,\left(,I,\right),+,\sum _I,g_{\text{FF}}^M,\left(,I,\right)\right)M+{\mathit{\gamma }}^M{\mathit{\theta }}^{\mathrm{Fe},Z}{R}_{\text{up}}^M\\ & +{\mathit{\lambda }}_{\text{POC}}^{\star }\mathrm{SFe}\\ & -(\mathrm{1}-{\mathit{\delta }}^P){\mathit{\mu }}^{P^{\mathrm{Fe}}}P-(\mathrm{1}-{\mathit{\delta }}^D){\mathit{\mu }}^{{\mathrm{D}}^{\mathrm{Fe}}}D-\text{Scav}-\\ & & \text{Cgfe}\mathrm{1}-\text{Cgfe}\mathrm{2}-\text{Aggfe}-\text{Bactfe}\end{array}$$

Iron scavenging (Scav) has been described previously in Sect. . Iron is present in seawater largely as colloids e.g.,. 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 . In PISCES, we model this process following the approach chosen for DOM (see Sect. ): $$\begin{array}{cc}& \text{Cgfe}\mathrm{1}=\left((,a_{\mathrm{1}},\text{DOC},+,a_{\mathrm{2}},\text{POC},),\times ,\text{sh},+,a_{\mathrm{4}},\text{POC}\right.\\ & & \left..+,a_{\mathrm{5}},\text{DOC}\right)\times {\mathrm{Fe}}_{\text{coll}},& & \text{Cgfe}\mathrm{2}=a_{\mathrm{3}}\text{GOC}\times \text{sh}\times {\mathrm{Fe}}_{\text{coll}},\end{array}$$ where ${\mathrm{Fe}}_{\text{coll}}$ is computed from the iron chemistry model (see below).

When dissolved iron concentration exceeds the total ligand concentration $L_{\mathrm{T}}$, scavenging is enhanced as it is done in many other biogeochemical models e.g.,: $$\text{Aggfe}=1000{\mathit{\lambda }}^{\mathrm{Fe}}max\left(\mathrm{0},,,\mathrm{Fe},-,L_{\mathrm{T}}\right){\mathrm{Fe}}^{\prime }.$$

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 . Observations show that they have quite elevated Fe $/$ C ratios and account for a significant fraction of the total biological uptake of iron . The bacterial uptake of iron is parameterized according to $$\text{Bactfe}={\mathit{\mu }}_P{L}_{\text{lim}}^{\text{Bact}}{\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},\mathrm{Bact}}\frac{\mathrm{Fe}}{K_{\mathrm{Fe}}^{\text{B,1}}+\mathrm{Fe}}\text{Bact},$$ where ${\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},\mathrm{Bact}}$ 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 (${\mathrm{Fe}}^{\prime }$) and dissolved complexed iron ($\mathrm{FeL}$).

• The complex chemistry model of as modified by . This model is based on two ligands (${\mathrm{L}}_{\mathrm{W}}$ and ${\mathrm{L}}_{\mathrm{S}}$) and five iron forms: free $\mathrm{Fe}\left(\mathrm{II}\right)$ (${\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime }$) and $\mathrm{Fe}\left(\mathrm{III}\right)$ (${\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }$), $\mathrm{Fe}\left(\mathrm{III}\right)$ bound to the weak ligand (${\mathrm{FeL}}_{\mathrm{W}}$), $\mathrm{Fe}\left(\mathrm{III}\right)$ bound to the strong ligand (${\mathrm{FeL}}_{\mathrm{S}}$) and solid iron (${\mathrm{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 , whereas the complex model is detailed in . For the complex model, all chemical constants have identical values to what was chosen in and are thus not listed in Table a–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 ${\mathrm{Fe}}^{\prime }$ and of “complexed” iron $\mathrm{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 . 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 $$\begin{array}{cc}& L_{\mathrm{T}}=\mathrm{FeL}+L^{\prime }\\ & \mathrm{Fe}=\mathrm{FeL}+{\mathrm{Fe}}^{\prime }\\ & & K_{\text{eq}}^{\mathrm{Fe}}=\frac{\mathrm{FeL}}{L^{\prime }{\mathrm{Fe}}^{\prime }}.\end{array}$$

The chemical equilibrium constant $K_{\text{eq}}^{\mathrm{Fe}}$ is computed from the formulation proposed by . Solving this set of equations is equivalent to solve a second-order polynomial equation in ${\mathrm{Fe}}^{\prime }$, whose solution is $$\begin{array}{cc}& \mathrm{\Delta }=\mathrm{1}+K_{\text{eq}}^{\mathrm{Fe}}{L}_{\mathrm{T}}-K_{\text{eq}}^{\mathrm{Fe}}\mathrm{Fe}\\ & & {\mathrm{Fe}}^{\prime }=\frac{-\mathrm{\Delta }+\sqrt{{\mathrm{\Delta }}^{\mathrm{2}}+\mathrm{4}{K}_{\text{eq}}^{\mathrm{Fe}}\mathrm{Fe}}}{\mathrm{2}{K}_{\text{eq}}^{\mathrm{Fe}}}.\end{array}$$

Colloidal iron is assumed to represent 50 % of FeL: $${\mathrm{Fe}}_{\text{coll}}=0.5\mathrm{FeL}.$$

The total ligand concentration $L_{\mathrm{T}}$ can be either constant over the ocean, using a value defined in the namelist or can be variable using the relationship proposed by : $$L_{\mathrm{T}}=max\left(0.09,(,\text{DOC},+,40\right)-\mathrm{3},0.6),$$ where $L_{\mathrm{T}}$ is in nmol ${\mathrm{L}}^{-\mathrm{1}}$ and DOC in $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{L}}^{-\mathrm{1}}$.

Complex chemistry model

The iron chemical system is governed by the following set of four equations $$\begin{array}{cc}& \mathrm{0}=k_{\text{lW}}{\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }{L}_{\mathrm{W}}-k_{\text{bW}}{\mathrm{FeL}}_{\mathrm{W}}-k_{\text{phW}}{\mathrm{FeL}}_{\mathrm{W}}\\ & & -k_{\text{th}}{\mathrm{FeL}}_{\mathrm{W}},& & \mathrm{0}=k_{\text{lS}}{\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }{L}_{\mathrm{S}}-k_{\text{bS}}{\mathrm{FeL}}_{\mathrm{S}}-k_{\text{phS}}{\mathrm{FeL}}_{\mathrm{S}},& & \mathrm{0}=k_{\text{phW}}{\mathrm{FeL}}_{\mathrm{W}}+k_{\text{phS}}{\mathrm{FeL}}_{\mathrm{S}}+k_{\text{th}}{\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }-k_{\text{ox}}{\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime },& & \mathrm{0}=k_{\text{pcp}}{\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }-k_r{\mathrm{Fe}}_P.\end{array}$$ A supplementary reaction has been added relative to the original set of equations. In the Pacific Ocean, thermal (dark) reduction of $\mathrm{Fe}\left(\mathrm{III}\right)$ organic complexes has been shown to produce the accumulation of a sizeable amount of ${\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime }$ in the mesopelagic zone .

Additional constraints are given by the conservation of total dissolved iron ($\mathrm{Fe}$), $L_{\text{WT}}$ and $L_{\text{ST}}$ over the fast timescale: $$\begin{array}{}& & \mathrm{Fe}={\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }+{\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime }+{\mathrm{FeL}}_{\mathrm{W}}+{\mathrm{FeL}}_{\mathrm{S}}+{\mathrm{Fe}}_P,& & L_{\text{WT}}={\mathrm{FeL}}_{\mathrm{W}}+L_{\mathrm{W}},& & L_{\text{ST}}={\mathrm{FeL}}_{\mathrm{S}}+L_{\mathrm{S}}.\end{array}$$

Solving this system of equations is equivalent to solving a third-order polynomial equation in ${\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }$ (Eq. 16 in ). Because thermal aphotic reduction of ${\mathrm{FeL}}_{\mathrm{W}}$ has been added here, the definition of some coefficients in the original study has changed: $$\begin{array}{}& & b=\mathrm{1}+\frac{k_{\text{phW}}+k_{\text{th}}}{k_{\text{ox}}},& & K_{\mathrm{W}}=\frac{k_{\text{phW}}+k_{\text{th}}+k_{\text{bW}}}{k_{\text{lW}}},\end{array}$$ where $k_{\text{th}}$ has been set to 0.0048 ${\mathrm{h}}^{-\mathrm{1}}$. Then, knowing ${\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }$, the other four iron species can be computed.

Observations suggest that the weak ligand ($L_{\mathrm{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.,. In PISCES, we assume that two-thirds of the total ligand concentration above 0.6 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{1}}$ is going to $L_{\mathrm{S}}$, the rest is attributed to $L_{\mathrm{W}}$: $$\begin{array}{}& & L_{\text{WT}}=0.6+\frac{\mathrm{1}}{\mathrm{3}}(L_{\mathrm{T}}-0.6),& & L_{\text{ST}}=\frac{\mathrm{2}}{\mathrm{3}}(L_{\mathrm{T}}-0.6).\end{array}$$

As in the simple chemistry model, the ligand concentration $L_{\mathrm{T}}$ can be either constant over the ocean, using a value defined in the namelist or can be variable using the relationship proposed by (see Eq. ).

The rate constants required by the model are identical to those described by as modified by . Furthermore, we have slightly changed the formulation of the oxidation rate constant used in the original model: $$k_{\text{ox}}=k_{\text{ox}}^{\prime }\frac{max\left({\mathrm{O}}_{\mathrm{2}},,,\mathrm{1},,\mathrm{\mu },\mathrm{mol},,{\mathrm{L}}^{-\mathrm{1}}\right)}{{\mathrm{O}}_{\mathrm{2}\text{sat}}}.$$

This avoids numerical problems in strongly anoxic areas where oxygen concentration is close to 0. Bioavailable iron can be defined either as ${\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime }+{\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }+{\mathrm{FeL}}_{\mathrm{S}}$ or as ${\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime }+{\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }+{\mathrm{FeL}}_{\mathrm{S}}+{\mathrm{FeL}}_{\mathrm{W}}$. $k_{\text{th}}$ has assigned the value computed from the observations by , consistent with the data of . Colloidal iron and dissolved inorganic iron are defined as $$\begin{array}{}& & {\mathrm{Fe}}_{\text{coll}}=0.5({\mathrm{Fe}}_p+{\mathrm{FeL}}_{\mathrm{W}}+{\mathrm{FeL}}_{\mathrm{S}}),& & {\mathrm{Fe}}^{\prime }={\mathrm{Fe}\left(\mathrm{III}\right)}^{\prime }+{\mathrm{Fe}\left(\mathrm{II}\right)}^{\prime }.\end{array}$$

We assumed that 50 % of the iron bound to ligands and of the particulate inorganic iron is colloidal iron.

Si

$$\frac{\partial \mathrm{Si}}{\partial t}={\mathit{\lambda }}_{\mathrm{PSi}}^{\star }{\text{Diss}}_{\mathrm{Si}}\mathrm{PSi}-{\mathit{\theta }}_{\text{opt}}^{D,\text{Si}}(\mathrm{1}-{\mathit{\delta }}^D){\mathit{\mu }}^{D}D$$

All terms in this equation have been already defined previously.

Calcite

$$\frac{\partial {\mathrm{CaCO}}_{\mathrm{3}}}{\partial t}=P_{{\mathrm{CaCO}}_{\mathrm{3}}}-{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}^{\star }{\mathrm{CaCO}}_{\mathrm{3}}-w_{\text{GOC}}\frac{\partial {\mathrm{CaCO}}_{\mathrm{3}}}{\partial z}$$

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 $$\begin{array}{cc}{P}_{{\mathrm{CaCO}}_{\mathrm{3}}}=& R_{{\mathrm{CaCO}}_{\mathrm{3}}}\left({\mathit{\eta }}^Z,g^Z,\left(,P,\right),Z,+,{\mathit{\eta }}^M,g^M,\left(,P,\right),M\right.\\ & & \left..+,0.5,(,m^P,\frac{P}{K_{\mathrm{m}}+P},P,+,\text{sh},\times ,w^P,P^{\mathrm{2}},)\right).\end{array}$$

The rain ratio $R_{{\mathrm{CaCO}}_{\mathrm{3}}}$ is variable. We propose the following parameterization for this ratio: $$\begin{array}{cc}{R}_{{\mathrm{CaCO}}_{\mathrm{3}}}=& r_{{\mathrm{CaCO}}_{\mathrm{3}}}{L}_{\text{lim}}^{{\mathrm{CaCO}}_{\mathrm{3}}}\frac{T}{0.1+T}max\left(\mathrm{1},,,\frac{P}{\mathrm{2}}\right)\\ & \times \frac{max\left(\mathrm{0},,,\text{PAR},-,\mathrm{1}\right)}{\mathrm{4}+\text{PAR}}\frac{30}{30+\text{PAR}}\\ & \times \left(\mathrm{1},+,\exp ,\left(\frac{-(T-10{)}^{\mathrm{2}}}{25}\right)\right)\\ & & \times min\left(\mathrm{1},,,\frac{50}{z_{\text{mxl}}}\right).\end{array}$$

This parameterization is based on a set of very simple assumptions, mainly inferred from the review by :

• 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 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, respectively.

• Coccolithophores are not found when the temperature of sea water is below 0 ${}^{\circ }\mathrm{C}$.

• Coccolithophores are found in stratified waters. Their abundance decreases when the mixed layer depth ($z_{\text{mxl}}$) exceeds 50 m.

• Maximum levels of coccolithophores are found in the mid-latitudes, where temperature is around 10 ${}^{\circ }\mathrm{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 (${\mathit{\eta }}^I$) of the grazed shells are routed to sinking calcite. The rest is taken to dissolve in the acidic guts of zooplankton . 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 .

The dissolution of calcite is modeled as in : $$\begin{array}{}& & \mathrm{\Delta }{\text{CO}}_{\mathrm{3}}^{\mathrm{2}-}=max\left(\mathrm{0},,,\mathrm{1},-,\frac{{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}}{{\mathrm{CO}}_{\mathrm{3},\mathrm{sat}}^{\mathrm{2}-}}\right),& & {\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}^{\star }={\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}(\mathrm{\Delta }{\text{CO}}_{\mathrm{3}}^{\mathrm{2}-}{)}^{\text{nca}}.\end{array}$$

The carbonate system

$$\begin{array}{cc}& \frac{\partial \text{DIC}}{\partial t}={\mathit{\gamma }}^Z(\mathrm{1}-e^Z-{\mathit{\sigma }}^Z)\sum _I{g}^Z\left(I\right)Z+{\mathit{\gamma }}^M(\mathrm{1}-e^M-{\mathit{\sigma }}^M)\\ & \left(\sum _I,g^M,\left(,I,\right),+,\sum _I,g_{\text{FF}}^M,\left(,I,\right)\right)M\\ & +{\mathit{\gamma }}^M{R}_{\text{up}}^M+\text{Remin}+\text{Denit}+{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}^{\star }{\mathrm{CaCO}}_{\mathrm{3}}-P_{{\mathrm{CaCO}}_{\mathrm{3}}}\\ & & -{\mathit{\mu }}^{D}D-{\mathit{\mu }}^{P}P& \frac{\partial \text{Alk}}{\partial t}={\mathit{\theta }}^{\mathrm{N},\mathrm{C}}\text{Remin}+{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}(r_{{\mathrm{NO}}_{\mathrm{3}}}^{\star }+\mathrm{1})\text{Denit}\\ & +{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\mathit{\gamma }}^Z(\mathrm{1}-e^Z-{\mathit{\sigma }}^Z)\sum _I{g}^Z\left(I\right)Z\\ & +{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\mathit{\gamma }}^M(\mathrm{1}-e^M-{\mathit{\sigma }}^M)\\ & \left(\sum _I,g^M,\left(,I,\right),+,\sum _I,g_{\text{FF}}^M,\left(,I,\right),+,{\mathit{\theta }}^{\mathrm{N},\mathrm{C}},{\mathit{\gamma }}^M,R_{\text{up}}^M\right)M\\ & +{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^{P}P+{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^{D}D+{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{N}_{\text{fix}}+\mathrm{2}{\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}^{\star }{\mathrm{CaCO}}_{\mathrm{3}}\\ & +{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}\mathrm{\Delta }\left({\mathrm{O}}_{\mathrm{2}}\right)(r_{{\mathrm{NH}}_{\mathrm{4}}}^{\star }-\mathrm{1}){\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}{\mathrm{NH}}_{\mathrm{4}}\\ & -{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^{P}P-{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}{\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^{D}D\\ & & -\mathrm{2}{\mathit{\theta }}^{\mathrm{N},\mathrm{C}}\text{Nitrif}-\mathrm{2}{P}_{{\mathrm{CaCO}}_{\mathrm{3}}}\end{array}$$

All terms in the above equations have been described previously in this document. In addition to these biogeochemical fluxes, the ocean exchanges CO${}_{\mathrm{2}}$ with the atmosphere at the sea surface. The gas exchange coefficient is computed from the relationship proposed by . No exchange is allowed with the atmosphere across sea ice: $$k_{{\mathrm{gCO}}_{\mathrm{2}}}=k_{{\mathrm{gCO}}_{\mathrm{2}}}^{\prime }\times (\mathrm{1}-{\mathit{\%}}_{\text{ice}}),$$ where ${\mathit{\%}}_{\text{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 $p$${\mathrm{CO}}_{\mathrm{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 $p$${\mathrm{CO}}_{\mathrm{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

$$\begin{array}{cc}& \frac{\partial {\mathrm{O}}_{\mathrm{2}}}{\partial t}={\mathrm{O}}_{\mathrm{2}}^{\text{ut}}({\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^{P}P+{\mathit{\mu }}_{{\mathrm{NH}}_{\mathrm{4}}}^{D}D)+({\mathrm{O}}_{\mathrm{2}}^{\text{ut}}+{\mathrm{O}}_{\mathrm{2}}^{\text{nit}})\\ & ({\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^{P}P+{\mathit{\mu }}_{{\mathrm{NO}}_{\mathrm{3}}}^{D}D)+{\mathrm{O}}_{\mathrm{2}}^{\text{nit}}{N}_{\text{fix}}\\ & -{\mathrm{O}}_{\mathrm{2}}^{\text{ut}}{\mathit{\gamma }}^Z(\mathrm{1}-e^Z-{\mathit{\sigma }}^Z)\sum _I{g}^Z\left(I\right)Z-{\mathrm{O}}_{\mathrm{2}}^{\text{ut}}{\mathit{\gamma }}^M\\ & (\mathrm{1}-e^M-{\mathit{\sigma }}^M)\left(\sum _I,g^M,\left(,I,\right),+,\sum _I,g_{\text{FF}}^M,\left(,I,\right)\right)M-{\mathrm{O}}_{\mathrm{2}}^{\text{ut}}{\mathit{\gamma }}^M{R}_{\text{up}}^M\\ & & -{\mathrm{O}}_{\mathrm{2}}^{\text{ut}}\text{Remin}-{\mathrm{O}}_{\mathrm{2}}^{\text{nit}}\text{Nitrif}\end{array}$$

In this equation, the stoichiometric ratio ${\mathrm{O}}_{\mathrm{2}}^{\text{ut}}$ represents the change in oxygen relative to carbon when ammonium is converted to organic matter, whereas ${\mathrm{O}}_{\mathrm{2}}^{\text{nit}}$ 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 .

Oxygen is exchanged with the atmosphere using the parameterization of 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${}_{\mathrm{2}}$, no air–sea fluxes are allowed when the ocean is covered by sea ice (see Eq. ).

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. ) 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 $\mathrm{\%}$ which is the average content of crustal material e.g.,. 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 $\mathrm{\%}$ of the particulate iron dissolves in a day . This is equivalent to a remineralization length scale of 20 000 $\mathrm{m}$ if the sinking speed is set to a typical value of 2 $\mathrm{m}{\mathrm{day}}^{-\mathrm{1}}$, on the same order as the length scale prescribed for the same process by . Atmospheric deposition of Si is also considered following 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 $\mathrm{ppm}$ and that the solubility in surface sea water is 10 $\mathrm{\%}$ . 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 $$\text{Dust}=\frac{D_{\text{dust}}}{w_{\text{dust}}},$$ where $D_{\text{dust}}$ is dust deposition at the surface and $w_{\text{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 . For DIC and alkalinity, we use results from the Global Erosion Model (GEM) of , neglecting the POC delivery as most of it is lost in the estuaries and in the coastal zone . 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 $\mathrm{Tg}\mathrm{Fe}{\mathrm{yr}}^{-\mathrm{1}}$ as estimated by .

Reductive mobilization of iron from marine sediments

Reductive mobilization of iron from marine sediments have been recognized as a significant source to the ocean . 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 . 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 and perhaps the magnitude of the organic carbon export to the sediments . Such potential relationships are not yet embedded in PISCES.

In a way similar to , 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 : $$\begin{array}{cc}& & z_{\text{Fsed}}=min\left(\mathrm{8},,,{\left(\frac{z}{500\mathrm{m}}\right)}^{-1.5}\right),& {\mathit{\zeta }}_{\text{Fsed}}=-0.9543+0.7662\ln \left(z_{\text{Fsed}}\right)\\ & & -0.235{\left(\ln ,\left(,z_{\text{Fsed}},\right)\right)}^{\mathrm{2}},& & \text{Fsed}=min\left(\mathrm{1},,,\exp ,\left(,{\mathit{\zeta }}_{\text{Fsed}},\right),/,0.5\right).\end{array}$$

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 . Our modulation factor is simply set equal to this relative contribution. The maximum iron flux from the sediments has been set by default to 2 $\mathrm{\mu }\mathrm{mol}\mathrm{Fe}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{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 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 $\mathrm{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 high-resolution 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 ${\mathit{\%}}_{\text{sed}}$ (which varies between 0 and 1): $$F_{\mathrm{Fe}}^{\text{sed}}=F_{\mathrm{Fe},\max }^{\text{sed}}\times \text{Fsed}\times {\mathit{\%}}_{\text{sed}}.$$

This corresponds to a global flux of 34 Gmol Fe yr${}^{-\mathrm{1}}$.

Iron from hydrothermalism

Recent studies have shown that hydrothermalism may deliver to the deep ocean a significant amount of dissolved iron e.g.,. Despite very large uncertainties, this source has been estimated, based on discrete data and a model, to 3 to $\mathrm{9}\times {10}^{\mathrm{8}}$ $\mathrm{mol}$ Fe ${\mathrm{yr}}^{-\mathrm{1}}$ globally . In PISCES, this source is included following the modeling study by 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 ${}^{\mathrm{3}}$He and dFe and using a data compilation of dFe $/$ ${}^{\mathrm{3}}$He (see the Supplement of ). Then, the spatial distribution of this flux has been derived from previous modeling works on ${}^{\mathrm{3}}$He, which relate the ${}^{\mathrm{3}}$He flux to the ridge-spreading rates ; 0.2 $\mathrm{\%}$ of the delivered iron is assumed to be soluble.

${\mathit{\theta }}_{\text{opt}}^{\mathrm{Si},D}$ as a function of Si concentration and $F_{\text{lim,1}}^{D^{\mathrm{Si}}}$. The vertical axis corresponds to $\log \left(\mathrm{Si}\right)$.

[Figure omitted. See PDF]

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 as well as the releases of iron accumulated in sea ice during winter . Measurements in sea ice have found iron concentrations of more than 1 order of magnitude higher than in adjacent sea water . About 90 $\mathrm{\%}$ of this iron has been shown to be of oceanic origin . Thus, iron is taken up from sea water when ice forms and is released back to the ocean when it melts. 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: $$\begin{array}{}& & F_{\mathrm{Fe}}^{\text{ice},-}=min\left(\mathrm{0},,,-,{\text{EP}}_{\text{oi}}\right)\times \mathrm{Fe},& & F_{\mathrm{Fe}}^{\text{ice},+}=max\left(\mathrm{0},,,-,{\text{EP}}_{\text{oi}}\right)\times {\mathrm{Fe}}_{\text{ice}},& & F_{\mathrm{Fe}}^{\text{ice}}=F_{\mathrm{Fe}}^{\text{ice},-}+F_{\text{ice}}^{\text{ice},+},\end{array}$$ where ${\text{EP}}_{\text{oi}}$ is the water flux (in kg ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) from the ice to the ocean and ${\mathrm{Fe}}_{\text{ice}}$ is the iron concentration in sea ice which has been found to be on the order of 10 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{1}}$. In this equation, $F_{\mathrm{Fe}}^{\text{ice},-}$ is thus the loss of iron from the ocean when sea ice forms and $F_{\mathrm{Fe}}^{\text{ice},+}$ is the release of iron to ocean 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${}^{-\mathrm{1}}$.

Dissolution rate of PSI (${\mathit{\lambda }}_{\mathrm{PSi}}^{\star }$) normalized to its value at 0 ${}^{\circ }\mathrm{C}$ with no silicate. Temperature is in ${}^{\circ }\mathrm{C}$.

[Figure omitted. See PDF]

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 with some modifications as described by . The main modification is the addition of denitrification to the set of early diagenetic reactions. Parameter values are identical to those in .

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 : $$E_{\text{burial}}=0.013+\frac{0.53F_{\text{OC}}^{\mathrm{2}}}{(7.0+F_{\text{OC}}{)}^{\mathrm{2}}},$$ where $E_{\text{burial}}$ is the burial efficiency and $F_{\text{OC}}$ is the flux of organic carbon at the bottom (in mmol C ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). We then use the metamodel by to determine the proportion of degradation of the remaining organic matter that is due to denitrification: $$\begin{array}{cc}\log \left(P_{\text{denit}}\right)=& -2.2567-1.185\log \left(F_{\text{OC}}\right)-0.221{\left(\log ,\left(,F_{\text{OC}},\right)\right)}^{\mathrm{2}}\\ & -0.3995\log \left({\mathrm{NO}}_{\mathrm{3}}\right)\log \left({\mathrm{O}}_{\mathrm{2}}\right)\\ & +1.25\log \left({\mathrm{NO}}_{\mathrm{3}}\right)+0.4721\log \left({\mathrm{O}}_{\mathrm{2}}\right)\\ & -0.0996\log \left(z\right)\\ & & +0.4256\log \left(F_{\text{OC}}\right)\log \left({\mathrm{O}}_{\mathrm{2}}\right),\end{array}$$ where the tracer concentrations are in $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{L}}^{-\mathrm{1}}$ and $F_{\text{OC}}$ is the flux of organic carbon at the bottom (in $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{cm}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). In this equation, oxygen and nitrate concentrations are not allowed to be below 10 $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{L}}^{-\mathrm{1}}$ and 1 $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{L}}^{-\mathrm{1}}$, respectively. Then, the fluxes of nitrate and oxygen to the sediment as a consequence of denitrification and oxic degradation, respectively, can be computed: $$\begin{array}{}& & F_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{denit}}=R_{{\mathrm{NO}}_{\mathrm{3}}}{P}_{\text{denit}}{F}_{\text{OC}},& & F_{{\mathrm{O}}_{\mathrm{2}}}^{\text{oxic}}={\mathrm{O}}_{\mathrm{2}}^{\text{ut}}(\mathrm{1}-P_{\text{denit}})F_{\text{OC}}.\end{array}$$

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: $${\mathit{\%}}_{{\mathrm{CaCO}}_{\mathrm{3}}}=min\left(\mathrm{1},,,1.3,\frac{0.2-\mathrm{\Omega }}{0.4-\mathrm{\Omega }}\right),$$ where $\mathrm{\Omega }$ is the calcite saturation level. This relationship has been deduced from the study by . The permanent burial of calcite is modulated by ${\mathit{\%}}_{{\mathrm{CaCO}}_{\mathrm{3}}}$. The amount of calcite that is not buried, instantaneously dissolves back to the ocean.

Boolean variables in the namelist. These variables activate functionalities of PISCES.

Boolean name	Description
ln_co2int	Read atmospheric pco2 from a file (T) or constant (F)
ln_presatm	Constant atmospheric pressure (F) or from a file (T)
ln_varpar	PAR made a variable fraction of shortwave (T) or not (F)
ln_newprod	Use Eq. () (T) or Eq. () for phytoplankton growth
ln_dust	Dust input from the atmosphere (T)
ln_solub	Variable solubility of iron in dust (T)
ln_river	River discharge of nutrient (T)
ln_ironsed	Sedimentary source of iron (T)
ln_ironice	iron input from sea ice (T)
ln_hydrofe	iron input from hydrothermalism (T)
ln_pisdmp	Relaxation of some tracers to a mean value (T) ${}^{*}$
ln_check_mass	Check mass conservation (T)

${}^{*}$ The frequency at which the restoring technique is applied is specified by the parameter nn_pisdmp.

Model parameters and their default values

Table a–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 control parameters are listed in Table . 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 .

Available CPP keys in PISCES.

CPP Key	Description
key_pisces	Activate the PISCES model
key_kriest	Activate the Kriest model (see Sect. )
key_sed	Activate the sediment model (see Sect. )

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 . This model is based on an ocean general circulation model OPA9, coupled with the sea ice model Louvain-la-Neuve Ice Model (LIM2) . The spatial resolution is about 2${}^{\circ }$ by 2${}^{\circ }$$\cos \mathrm{\Phi }$ (where $\mathrm{\Phi }$ is the latitude) with a focusing of the meridional resolution to 0.5${}^{\circ }$ in the equatorial domain. The model has 30 vertical layers, with an increased vertical thickness from 10 $\mathrm{m}$ at the surface to 500 $\mathrm{m}$ at 5000 $\mathrm{m}$. Representation of the topography is based on the partial step thicknesses . Lateral mixing along isopycnal surfaces is performed both on tracers and momentum as in . The parameterization of is applied poleward of 10${}^{\circ }$ to represent the effects of non-resolved mesoscale eddies. Vertical mixing is parameterized using the turbulent kinetic energy (TKE) scheme of , as modified by .

The fields used to drive the ocean are identical to those used by . 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 $\mathrm{m}$ atmospheric temperature averaged over 1948–2003 , monthly relative humidity , monthly ISCCP total cloudiness averaged over 1983–2001 , monthly precipitation averaged over 1979–2001 and weekly wind stress based on European Remote-Sensing Satellite (ERS) satellite product and TAO observations . Surface heat fluxes and evaporation are computed using empirical bulk formulas as described by . To avoid any strong model drift, modeled sea surface salinity is restored to the monthly WOA01 data set with a nudging timescale of 40 $\mathrm{days}$ applied through local freshwater forcing (thereby conserving salt). The ocean dynamical model has been spun-up for 200 $\mathrm{years}$, starting from rest and from the climatology of for temperature and salinity.

Phosphate, oxygen, nitrate and silicic acid distributions have been initialized at uniform concentrations inferred from observed climatologies . Initial values for dissolved inorganic carbon and alkalinity are taken from the OCMIP guidelines . The ecological tracers are initialized uniformly to arbitrary low values. Iron concentrations are set everywhere to 0.6 $\mathrm{nM}$. The model is then spun-up offline for 4000 $\mathrm{years}$ using the circulation state predicted by the dynamical model. Atmospheric $p$${\mathrm{CO}}_{\mathrm{2}}$ is set to a pre-industrial value of 278 $\mathrm{ppm}$. After this integration, primary productivity as well as ${\mathrm{CO}}_{\mathrm{2}}$ fluxes drift by less than 0.001 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{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–$\mathrm{5}\times {10}^{-\mathrm{5}}$ yr${}^{-\mathrm{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 ).

Global budget

Global annual budget of C in the top 150 $\mathrm{m}$ of the ocean.

Carbon budget ${}^{\mathrm{a}}$
Primary production in the top 150 m of the ocean
	7.5	Primary production by diatoms
	36.8	Primary production by nanophytoplankton
44.3		Global total primary production
Export from the top 150 m of the ocean
	3.9	Vertical flux due to sinking big POC
	2	Vertical flux due to sinking small POC
	1	Advective/diffusive vertical flux of organic matter
6.9 ${}^b$		Total vertical flux of organic matter
Various fluxes in the top 150 m of the ocean
	35.8	Grazing by microzooplankton on phytoplankton
	40.2	Total grazing by microzooplankton
	4	Grazing by mesozooplankton on phytoplankton
	11.2	Total grazing by mesozooplankton
51.2		Total grazing by zooplankton
22.3		Remineralization of DOC

${}^{\mathrm{a}}$ Carbon fluxes are all in $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$. ${}^{\mathrm{b}}$ The total vertical flux due to sinking POC is 7.3 Gt C yr${}^{-\mathrm{1}}$ at 100 m depth.

Table  presents the global carbon budget as simulated by PISCES, when embedded in ORCA2-LIM. The annual net predicted primary production is 44 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$. This value falls on the lower bound of the broad estimates given by satellite observations which give values between 37 and 67 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ . 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 , and coastal regions, characterized by a intense primary productivity, are not properly resolved by the coarse grid.

About 17 $\mathrm{\%}$ of the primary production is due to diatoms. Global estimates of the contribution of diatoms to total production are rather uncertain and broad. have suggested that diatoms may be responsible for up to 40 $\mathrm{\%}$ of the total primary production. However, as discussed by , this value is certainly overestimated. In recent years, algorithms, which attempt to retrieve the composition of phytoplankton from space, have been developed e.g.,. 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 . The estimated global contribution of diatoms from these methods ranges from as low as 7 $\mathrm{\%}$ to as high as 32 $\mathrm{\%}$ of the total phytoplankton (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 $\mathrm{\%}$ e.g.,.

Export production at 150 $\mathrm{m}$ is estimated to be 6.9 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$; 86 $\mathrm{\%}$ of this export is related to settling particles (one-third by the small sinking particles and two-third 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.,. Mesozooplankton grazes about 9 $\mathrm{\%}$ of total primary production. This value is close to other estimates either based on observations or models . Total gazing by mesozooplankton is predicted to be 11.2 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ by PISCES, quite similar to the value of $10.4\pm 3.7$ Gt C ${\mathrm{yr}}^{-\mathrm{1}}$ estimated by for the global respiration of mesozooplankton in the upper 200 $\mathrm{m}$ of the ocean. About 80 $\mathrm{\%}$ of total primary production, i.e., 35.8 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$, is consumed up by microzooplankton above the upper bound of the 25–33 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ given by when extrapolating observations. Despite estimates of grazing by microzooplankton are quite badly constrained, this might suggest that it is overestimated in the model.

Global annual budget of calcite and Si in the top 150 $\mathrm{m}$ of the ocean.

Calcite budget${}^{\mathrm{a}}$
	1.6	Production of calcite
	0.8	Dissolution of calcite
0.8		Vertical flux of sinking calcite particles
Biogenic silica budget${}^{\mathrm{b}}$
	145.8	Production of $\mathrm{BSi}$
	39.6	Dissolution of $\mathrm{BSi}$
106.2		Vertical flux of dissolved $\mathrm{BSi}$

${}^{\mathrm{a}}$ Calcite fluxes are all in $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$. ${}^{\mathrm{b}}$ Biogenic silica fluxes are all in Tmol Si ${\mathrm{yr}}^{-\mathrm{1}}$.

Table  shows the calcite and silicon budgets for the upper 150 $\mathrm{m}$ of the ocean. Production of calcite and export at 150 $\mathrm{m}$ are simulated to be, respectively, about 1.6 and 0.8 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ by PISCES. These numbers fall within the limits of the quite large range of 0.4–1.8 $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$ estimated either for global calcification or export of particulate inorganic carbon (PIC) . For silicate, the model predicts a vertical export of biogenic silicate of 106 Tmol Si ${\mathrm{yr}}^{-\mathrm{1}}$. This value is within the $105\pm 17$ Tmol Si ${\mathrm{yr}}^{-\mathrm{1}}$ estimated for the global ocean . Global production of biogenic silica by diatoms is 146 Tmol Si ${\mathrm{yr}}^{-\mathrm{1}}$ in our model. This value is quite low compared to the 239 Tmol Si ${\mathrm{yr}}^{-\mathrm{1}}$ given by . About 27 $\mathrm{\%}$ of biogenic silica dissolves in the top 150 $\mathrm{m}$ of the ocean, half the estimate of and . 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), estimated biogenic silica production to be about 103 Tmol Si ${\mathrm{yr}}^{-\mathrm{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 . 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 occurs within the Southern Ocean, a region where diatom cells are very heavily silicified .

Annual budget${}^{*}$ of N over the global ocean.

Sources of nitrogen to the ocean
	36	River discharge
	67	Atmospheric deposition
	111.8	Nitrogen fixation
214.8		Total input of nitrogen
Sinks of nitrogen from the ocean
	77.6	Denitrification in the water column
	92.8	Denitrification in the sediments
	23.2	Permanent burial in the sediments
193.6		Total loss of nitrogen
21.2		Net budget of nitrogen (Sources minus Sinks)

${}^{*}$ All nitrogen fluxes are in $\mathrm{Tg}\mathrm{N}{\mathrm{yr}}^{-\mathrm{1}}$.

Table  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.8 $\mathrm{Tg}\mathrm{N}{\mathrm{yr}}^{-\mathrm{1}}$. This value is close to the mean value of about 140 $\mathrm{Tg}\mathrm{N}{\mathrm{yr}}^{-\mathrm{1}}$ estimated from direct observations or nutrients analysis . Figure  shows a comparison between the spatial distribution of observed nitrogen fixation rates from the MARine 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 $\mathrm{Tg}\mathrm{N}{\mathrm{yr}}^{-\mathrm{1}}$, respectively. Sediment denitrification estimates are significantly higher, in the range of 130–300 $\mathrm{Tg}\mathrm{N}{\mathrm{yr}}^{-\mathrm{1}}$ . 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 $\mathrm{Tg}\mathrm{N}{\mathrm{yr}}^{-\mathrm{1}}$.

Modeled tracer distributions

Chlorophyll

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

[Figure omitted. See PDF]

The modeled chlorophyll distribution is compared to GLOBCOLOUR satellite observations for two seasons in Fig. . 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 either in phytoplankton or in organic matter . 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 . 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 . Chlorophyll concentrations are high in the eastern boundary upwelling systems. The sedimentary source of iron plays a critical role in these systems. When this iron source is not included in models, modeled chlorophyll concentrations are much lower .

Annual-mean depth averaged ${\mathrm{N}}_{\mathrm{2}}$ fixation rates in $\mathrm{\mu }\mathrm{mol}$ N ${\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$. (a) Database from the MARine Ecosystem Model Intercomparison Projec (MAREMIP) project (Luo et al., 2013); (b) model predictions.

[Figure omitted. See PDF]

Surface seasonal mean chlorophyll concentrations (mg chl $a$ ${\mathrm{m}}^{-\mathrm{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 omitted. See PDF]

In two of the 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.,. 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 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{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 ). 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 .

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.,. 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

Spatial distribution of annual-mean iron concentrations (in nmol ${\mathrm{L}}^{-\mathrm{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 $\mathrm{m}$ of the ocean. On panels (b) and (c), iron has been averaged over 200–1000 $\mathrm{m}$. The bottom two panels display the iron distributions average over the depth range 1000–5000 $\mathrm{m}$. Model values have been sampled at the same location and month as the data.

[Figure omitted. See PDF]

Figure  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 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 with enhanced scavenging by the dust particles deposited at the ocean surface . Very high iron concentrations, typically above 1 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{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 , in the eastern boundary upwelling systems but also downstream of the islands, especially in the Southern Ocean . In the rest of the open ocean, iron concentrations are typically low, generally below 0.2 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{1}}$, especially in the HNLC regions. PISCES tends to exaggerate these low concentrations.

Annual-mean ${\mathrm{NO}}_{\mathrm{3}}$ concentrations in $\mathrm{\mu }\mathrm{mol}$ N ${\mathrm{L}}^{-\mathrm{1}}$. Observations are from the World Ocean Atlas 2009 . (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 omitted. See PDF]

Iron concentrations increase with depth due to the remineralization of organic particles settling from the surface waters . However, except near the coasts, concentrations rarely exceed 1 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{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 . In the deepest waters of the Pacific and Indian oceans, iron concentrations tend to decrease to the bottom of the ocean and they often fall below 0.6 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{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.,. 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 . 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 $\mathrm{nmol}{\mathrm{L}}^{-\mathrm{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 , and the GLobal Ocean Data Analysis Project (GLODAP) database for DIC and alkalinity .

Annual-mean ${\mathrm{SiO}}_{\mathrm{3}}$ concentrations in $\mathrm{\mu }\mathrm{mol}$ Si ${\mathrm{L}}^{-\mathrm{1}}$. Observations are from the World Ocean Atlas 2009 . Panels are the same as on Fig. .

[Figure omitted. See PDF]

Figures  and  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 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 . 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.

Annual-mean ${\mathrm{O}}_{\mathrm{2}}$ concentrations in $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{L}}^{-\mathrm{1}}$. Observations are from the World Ocean Atlas 2009 . Panels are the same as on Fig. .

[Figure omitted. See PDF]

In Fig. , the modeled oxygen distribution is evaluated against observations. Not surprisingly, the surface distribution compares quite well to the observations as oxygen is 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 . The nearly homogeneous oxygen concentrations south of 60${}^{\circ }$ S are a consequence of this too intense winter mixing, which thus ventilates the deep ocean with too much oxygen.

Annual-mean natural DIC concentrations in $\mathrm{\mu }\mathrm{mol}$ ${\mathrm{L}}^{-\mathrm{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. .

[Figure omitted. See PDF]

Annual-mean alkalinity concentrations in $\mathrm{\mu }\mathrm{mol}$ eq ${\mathrm{L}}^{-\mathrm{1}}$. Observations are from GLODAP. Panels are the same as on Fig. .

[Figure omitted. See PDF]

Figures  and  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 ${\mathrm{CO}}_{\mathrm{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 . It should be also mentioned here that no observations were available north of 60${}^{\circ }$ 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 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 , Taylor diagrams have become quite popular as they synthesize, in a quite convenient way, several statistical diagnostics.

Figures  and  show Taylor diagrams for surface chlorophyll and mesozooplankton averaged over the top 150 $\mathrm{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 annual-mean 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. ). 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) produce a much greater spatial and temporal variability for mesozooplankton, especially in the high latitudes and along the continental margins.

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${}^{\circ }$ S).

[Figure omitted. See PDF]

Figure  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 re-inspection of Fig. , 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 . 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 . The agreement between the model and this database is much better (see Fig. ), thus confirming that interpolation in GLODAP potentially leads to a strong underestimate of the real spatial variability.

Sensitivity tests with some new parameterizations

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${}^{\circ }$ S).

[Figure omitted. See PDF]

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 $\mathrm{years}$ in which specific parameterizations have been either changed or removed. Table  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.

Sensitivity experiments performed with PISCES to evaluate the impact of specific parameterizations. Primary production (PP) and export production at 150 $\mathrm{m}$ (EP) are in $\mathrm{Gt}\mathrm{C}{\mathrm{yr}}^{-\mathrm{1}}$.

Experiment	Description	Parameterization choices	PP	EP
PAR	Impact of variable PAR fraction	ln_varpar $\mathrm{=}$ .false.	44.4	5.8
LIGHT	Impact of light limitation	Eq. ()	42.6	7.3
SIZE	Impact of variable cell sizes	xsizern, xsizerd $=\mathrm{1}$	44.8	6.2
FOOD	Impact of food quality	${\mathit{\theta }}^{\mathrm{N},I}=0.136$	43.4	6.1

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 $\mathrm{m}$ of the ocean. Black dot corresponds to ${\mathrm{NO}}_{\mathrm{3}}$, brown dot to ${\mathrm{O}}_{\mathrm{2}}$, red dot to ${\mathrm{PO}}_{\mathrm{4}}$, green dot to ${\mathrm{SiO}}_{\mathrm{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 to compare with the model.

[Figure omitted. See PDF]

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 $\mathrm{\%}$, 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 ). 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 (see Eq. ). 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 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 ${}^{\circ }\mathrm{C}$ for no nutrient limitation.

Surface seasonal mean chlorophyll anomaly (mg chl $a$ ${\mathrm{m}}^{-\mathrm{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 omitted. See PDF]

Figure a 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).

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 $\mathrm{mg}$ chl $a$ ${\mathrm{m}}^{-\mathrm{3}}$ are shown.

[Figure omitted. See PDF]

Figure  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 . 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 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. ). 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.,. 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 .

Figure c 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 phytoplankton 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 . 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 , 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 . 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. ). In the FOOD sensitivity experiment, the effect of food quality on the gross growth efficiency has been removed, i.e., $e_{\mathrm{N}}^Z$ is set to 1.

Surface chlorophyll concentrations are almost unaltered when the impact of food quality is removed (see Fig. e 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.

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.

[Figure omitted. See PDF]

Figure  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 $\mathrm{\%}$ because of more intense grazing by zooplankton. These changes are not perceptible in chlorophyll concentrations (at least with the color scale chosen on Fig. ) 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 $\mathrm{\%}$ in the subtropical gyres of all oceans and this increase even exceeds 200 $\mathrm{\%}$ 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 $\mathrm{mg}$ chl $a$ ${\mathrm{m}}^{-\mathrm{3}}$, whereas it is below 0.1 $\mathrm{mg}$ chl $a$ ${\mathrm{m}}^{-\mathrm{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 , 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 major 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. . 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 .

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) , and nitrous oxide (N${}_{\mathrm{2}}$O) . An explicit representation of paleo-proxies, such as $\mathit{\delta }{}^{13}\mathrm{C}$ , Pa $/$ Th , Nd , is also available.

PISCES is still in a phase of active development despite the fact that its development started more than 10 $\mathrm{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., . 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 . Second, almost all marine biogeochemical models 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.,. The conceptual schemes on which biogeochemical models, including PISCES, should then be revised, in particular the distinction between phytoplankton and microzooplankton.

Model structure

The model is being coded in FORTRAN 90. To activate PISCES, the cpp key key_pisces should be declared. Only the subroutines that compute the biological or chemical sources and sinks are considered to be part of PISCES. Thus, this excludes the computation of the advection–diffusion equation (the transport of the tracers), as it is not specific to PISCES. There are two types of subroutines: the initialization of the tracers and of the parameters and the computation of the various biogeochemical sources and sinks. The latter PISCES subroutines are called from within the ocean model time loop.

The objective here is not to precisely detail the PISCES code but rather to list the different modules and to briefly describe their role. All the subroutines that compute the biogeochemical sources/sinks are called from p4zsms which is then the main PISCES subroutine.

• p4zbio.F90: computation of the new tracer concentrations by summing up all the different sources and sinks;

• p4zche.F90: computation of the various chemical constants;

• p4zfechem.F90: computation of the iron chemistry. Scavenging of iron, aggregation of iron colloids;

• p4zflx.F90: air–sea fluxes of ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{2}}$;

• p4zint.F90: time interpolation of various terms (e.g., growth rate);

• p4zlim.F90: co-limitations of phytoplankton growth by the different nutrients;

• p4zlys.F90: calcite chemistry and dissolution;

• p4zmeso.F90: sources and sinks of mesozooplankton (mortality, grazing, etc.);

• p4zmicro.F90: sources and sinks of microzooplankton;

• p4zmort.F90: computation of the various mortality terms of nanophytoplankton and diatoms;

• p4zopt.F90: optical model and computation of the euphotic depth;

• p4zprod.F90: growth rate of the two phytoplankton groups;

• p4zrem.F90: remineralization of organic matter, dissolution of biogenic silica;

• p4zsed.F90: top and bottom boundary conditions of the biogeochemical tracers (deposition, sedimentary losses, etc.);

• p4zsink.F90: aggregation of organic matter, computation of the particles sinking speeds. Vertical sedimentation of particles using a MUSCL advection scheme;

• p4zsms.F90: main PISCES subroutine which calls the other subroutine.

Besides the subroutines listed above, several subroutines perform the model initialization. We will only discuss the initialization of the parameters necessary to PISCES. The tracers concentrations are excluded here as their initialization will of course vary with the ocean model.

• trcini.pisces.F90: initialization of various biogeochemical parameters. Allocation of the arrays used in PISCES. This subroutine also calls all the initialization subroutines included in the PISCES subroutines listed above.

• trcnam_pisces.F90: this subroutine reads the information necessary to write the netcdf files when IOM is not used.

• par_pisces.F90: it sets the PISCES parameters such as the number of tracers and the name of the indices, the number of additional diagnostics, etc.

• sms_pisces.F90: this subroutine defines some general PISCES variables and arrays and allocates them.

(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 organic and inorganic matter with their default values in PISCES. (d) Model parameters for various processes with their default values in PISCES.

(a)
Parameter	Coding name
$b_{\text{resp}}$	bresp
${\mathit{\alpha }}^I$	pislope; pislope2
${\mathit{\delta }}^I$	excret; excret2
$K_{{\mathrm{NH}}_{\mathrm{4}}}^{I,\text{min}}$	concnnh4; concdnh4
$K_{{\mathrm{NO}}_{\mathrm{3}}}^{I,\text{min}}$	concnno3; concdno3
$K_{\mathrm{Si}}^{\mathrm{1}}$	xksi1
$K_{\mathrm{Si}}^{\mathrm{2}}$	xksi2
$K_{\mathrm{Fe}}^{I,\text{min}}$	concnfer; concdfer
$S_{\text{rat}}^I$	xsizern; xsizerd
${\mathit{\theta }}_{\mathrm{m}}^{\mathrm{Si},D}$	grosip
${\mathit{\theta }}_{\text{opt}}^{\mathrm{Fe},I}$	qnfelim; qdfelim
${\mathit{\theta }}_{\text{max}}^{\mathrm{Fe},I}$	fecnm; fecdm
$m^I$	mprat; mprat2
$w^P$	wchl
$w_{\text{max}}^D$	wchld
$K_{\mathrm{m}}$	xkmort
${\mathit{\theta }}_{\text{max}}^{\text{Chl},I}$	chlcnm; chlcdm
${\mathit{\theta }}_{\text{min}}^{\text{Chl}}$	chlcmin
$I_{\text{max}}$	xsizephy; xsizedia

Continued.

(b)
Parameter	Coding name
$e_{\text{max}}^I$	epsher; epsher2
${\mathit{\sigma }}^I$	unass; unass2
${\mathit{\gamma }}^I$	sigma; sigma2
$g_{\mathrm{m}}^I$	graze; graze2
$g_{\text{FF}}^M$	grazflux
$K_{\mathrm{G}}^I$	xkgraz; xkgraz2
$p_P^I$	xpref2p; xprefp
$p_D^I$	xpref2d; xprefc
$p_{\text{POC}}^I$	xpref2c; xprefpoc
$p_Z^M$	xprefz
$F_{\text{thresh}}^I$	xthresh; xthresh2
$m^I$	mzrat; mzrat2
$r^I$	resrat; resrat2
${\mathit{\nu }}^I$	part; part2
${\mathit{\theta }}^{\mathrm{Fe},\mathrm{Zoo}}$	ferat3

Continued.

(c)
Parameter	Coding name
${\mathit{\lambda }}_{\text{DOC}}$	xremik
$K_{\text{DOC}}$	xkdoc
$K_{{\mathrm{NO}}_{\mathrm{3}}}^{\text{Bact}}$	concbno3
$K_{{\mathrm{NH}}_{\mathrm{4}}}^{\text{Bact}}$	concbnh4
$K_{\mathrm{Fe}}^{\text{Bact}}$	concbfe
${\mathit{\lambda }}_{\text{POC}}$	xremip
$w_{\text{POC}}$	wsbio
$w_{\text{GOC}}^{\text{min}}$	wsbio2
$w_{\text{dust}}$	wdust
${\mathit{\lambda }}_{\mathrm{Fe}}$	xlam1
${\mathit{\lambda }}_{\mathrm{Fe}}^{\text{dust}}$	xlamdust
${\mathit{\lambda }}_{{\mathrm{CaCO}}_{\mathrm{3}}}$	kdca
nca	nca
${\mathit{\chi }}_{\text{lab}}^{\mathrm{0}}$	xsilab
${\mathit{\lambda }}_{\mathrm{PSi}}^{\text{slow}}$	xsirem
${\mathit{\lambda }}_{\mathrm{PSi}}^{\text{fast}}$	xsiremlab

Continued.

(d)
Parameter	Coding name
${\mathit{\lambda }}_{{\mathrm{NH}}_{\mathrm{4}}}$	nitrif
${{\mathrm{O}}_{\mathrm{2}}}^{\text{min,1}}$	oxymin
$L_{\mathrm{T}}$	ligand
$N_{\text{fix}}^{\mathrm{m}}$	nitrfix
$K_{\mathrm{Fe}}^{\text{Dz}}$	concfediaz
$E_{\text{fix}}$	diazolight
${\mathrm{Fe}}_{\text{ice}}$	icefeinput
$F_{\mathrm{Fe},\min }^{\text{sed}}$	sedfeinput
${\mathrm{Sol}}_{\mathrm{Fe}}$	dustsolub
$r_{{\mathrm{CaCO}}_{\mathrm{3}}}$	caco3r

Many parameter values of the model can be specified from the namelist namelist_pisces. When such is the case, the corresponding parameter name in the namelist file is indicated in Table a–d.

Acknowledgements

The authors are grateful to the whole community of PISCES users. The model would be useless without them. In particular, we thank Thomas Gorguès, Keith B. Rodgers, Coralie Perruche, Christophe Menkès, Vincent Echevin, Gildas Cambon and the whole NEMO system team for their precious help, expertise and support which made the release of this version possible. O. Aumont, L. Bopp, M. Gehlen were supported by ANR-CEP09 MACROES. O. Aumont received additional support from the Labex Mer via grant ANR-10-LABX-19-01. Olivier Aumont will be eternally grateful to Ernst Maier-Reimer. Even if he did not directly participate in the design of PISCES, nothing would have been possible without him. PISCES was built upon HAMOCC3 and HAMOCC3.1 and still some parts of the current code of PISCES come from these models. This work is dedicated to his memory.Edited by: A. Ridgwell