Ocean biogeochemistry in the coupled ocean–sea ice–biogeochemistry model FESOM2.1–REcoM3

The cycling of carbon in the oceans is affected by feedbacks driven by changes in climate and atmospheric CO₂. Understanding these feedbacks is therefore an important prerequisite for projecting future climate. Marine biogeochemistry models are a useful tool but, as with any model, are a simplification and need to be continually improved. In this study, we coupled the Finite-volumE Sea ice–Ocean Model (FESOM2.1) to the Regulated Ecosystem Model version 3 (REcoM3). FESOM2.1 is an update of the Finite-Element Sea ice–Ocean Model (FESOM1.4) and operates on unstructured meshes. Unlike standard structured-mesh ocean models, the mesh flexibility allows for a realistic representation of small-scale dynamics in key regions at an affordable computational cost. Compared to the previous coupled model version of FESOM1.4–REcoM2, the model FESOM2.1–REcoM3 utilizes a new dynamical core, based on a finite-volume discretization instead of finite elements, and retains central parts of the biogeochemistry model. As a new feature, carbonate chemistry, including water vapour correction, is computed by mocsy 2.0. Moreover, REcoM3 has an extended food web that includes macrozooplankton and fast-sinking detritus. Dissolved oxygen is also added as a new tracer. In this study, we assess the ocean and biogeochemical state simulated with FESOM2.1–REcoM3 in a global set-up at relatively low spatial resolution forced with JRA55-do (Tsujino et al., 2018) atmospheric reanalysis. The focus is on the recent period (1958–2021) to assess how well the model can be used for present-day and future climate change scenarios on decadal to centennial timescales. A bias in the global ocean–atmosphere preindustrial CO₂ flux present in the previous model version (FESOM1.4–REcoM2) could be significantly reduced. In addition, the computational efficiency is 2–3 times higher than that of FESOM1.4–REcoM2. Overall, it is found that FESOM2.1–REcoM3 is a skilful tool for ocean biogeochemical modelling applications.

1 Introduction

There is an unequivocal consensus and concern about the effects of increasing greenhouse gases in the atmosphere due to human activities. Since the beginning of the industrial era (year 1750), the concentration of carbon dioxide (CO₂) in the air has substantially risen from 277 to 417.2 ppm (year 2022; Friedlingstein et al., 2022b). The ocean has taken up a remarkably constant fraction of 25 %–30 % of human CO₂ emissions from fossil fuel burning and land use change throughout time (Crisp et al., 2022). For the recent decade, 2012–2021, the rate of ocean anthropogenic carbon uptake (including the effects of climate change) amounted to 2.9 ± 0.4 PgC yr⁻¹ (26 % of the total CO₂ emissions; Friedlingstein et al., 2022b). A similar proportion was taken up by the terrestrial biosphere, amounting to 3.1 ± 0.6 PgC yr⁻¹ (2012–2021), but the total air-to-land CO₂ flux is substantially lower because of emissions from land use change, mainly deforestation, that amounted to 1.2 ± 0.7 PgC yr⁻¹ (2012–2021). The ocean carbon sink has grown over the past few decades in response to the near-exponential rise in CO₂ emissions (Friedlingstein et al., 2022b). While the global ocean carbon sink estimate is assigned an uncertainty of 0.4 PgC yr⁻¹ and medium confidence, regional patterns of the sink differ more strongly. This points to the balance between physical and biological processes, which are more difficult to model, as also illustrated in model deficiencies in accurately representing the seasonal cycle of pCO₂ and CO₂ fluxes (Mongwe et al., 2018). Both climate change and rising atmospheric CO₂ will feed back on the fraction of CO₂ emissions that will end up in the ocean over the next century (Friedlingstein et al., 2003; Canadell et al., 2021). Models are important tools for estimating how large these feedbacks are.

The flux of CO₂ between the atmosphere and ocean is controlled by two main mechanisms, namely the solubility pump and the biological pump. The solubility pump describes the air–sea CO₂ exchange that occurs to satisfy a thermodynamic equilibrium and the subsequent transport of carbon from the surface to the deep ocean with the overturning circulation. This leads to CO₂ uptake at mid- to high latitudes through high solubility in cold waters and large vertical motion in deep-water formation regions. In contrast, warm ocean regions in the tropics and subtropics and upwelling regions lose carbon to the atmosphere (Takahashi et al., 2009; Wanninkhof et al., 2013). The solubility pump is responsible for anthropogenic carbon uptake. The biological carbon pump comprises the fixation of CO₂ into biomass by phytoplankton and the subsequent downward transfer of dead organic material (Boyd et al., 2019). The biological carbon pump is responsible for 75 % of the natural vertical carbon gradient and for the large interbasin gradient between the deep Pacific and Atlantic (Sarmiento and Gruber, 2006). Without the biological carbon pump, atmospheric CO₂ would be higher by 200 ppm (Maier-Reimer et al., 1996), and perturbations thereof can have large effects on atmospheric CO₂ (Kwon et al., 2009; Lauderdale and Cael, 2021), as also known from paleo evidence (Galbraith and Skinner, 2020).

Global ocean biogeochemistry models (GOBMs; Fennel et al., 2022) are used to assess the global ocean carbon sink (Hauck et al., 2020), its regional patterns (Fay and McKinley, 2021), and effects of climate change and variability on the ocean carbon sink (Le Quéré et al., 2010; Hauck et al., 2013; DeVries et al., 2019; Bunsen, 2022). Through their representation of pH, the marine oxygen cycle, and phytoplankton primary production as the base of the marine food web, they also offer information about the environmental conditions for marine life and how these will develop under climate change (Bopp et al., 2013; Laufkötter et al., 2015; Kwiatkowski et al., 2020). However, modelling the marine biogeochemistry is subject to several sources of uncertainties. First, GOBMs are expensive with respect to the computational cost, due to the advection of a large number of tracers, and therefore often demand low spatial resolution. This leads to deficiencies in the representation of significant physical processes such as (sub)mesoscale currents (McWilliams, 2016), which can have large impacts on transport and mixing processes that strongly affect biological productivity (Lévy et al., 2018; Keerthi et al., 2022). Second, the descriptions of ecological interactions and of the physiology of primary and secondary producers in GOBMs are still mostly based on empirical or semi-empirical mathematical descriptions, such as the dependency of zooplankton grazing rates on prey abundance (Doney et al., 2001; Rohr et al., 2022). These contain a large number of parameters that are only partly constrained from observations, making it necessary to tune these parameters in GOBMs to some extent. Choices in these parameters can have strong effects on the biological carbon pump (e.g.  Lauderdale and Cael, 2021). It has been demonstrated that the largest source of uncertainty for projections of net primary production (NPP; Tagliabue et al., 2021) comes from model uncertainty and not scenario uncertainty (Frölicher et al., 2016).

Ocean circulation models formulated on unstructured meshes have become an alternative to the existing structured global ocean models (Danilov, 2013). The Finite-Element Sea ice–Ocean Model (hereafter FESOM1.4; Wang et al., 2014) is one of the first global models with multiple resolutions designed to simulate the large-scale ocean circulation. While it has already been used in numerous applications (Sidorenko et al., 2015; Wekerle et al., 2017), another dynamical core, the Finite-volumE Sea ice–Ocean Model version 2.1 (FESOM2.1), has been developed (Danilov et al., 2017). The advantages of a finite-volume formulation are (a) better throughput and scalability as a result of a more efficient data structure (Koldunov et al., 2019), (b) the availability of clearly defined fluxes, and (c) the possibility to choose from a selection of transport algorithms, which was very limited before (Danilov et al., 2017). Furthermore, the arbitrary Lagrangian–Eulerian (ALE) vertical coordinate is introduced, which provides different types of vertical coordinates (Scholz et al., 2019). The Regulated Ecosystem Model (REcoM) is an ocean biogeochemistry model that describes the lower trophic levels of the marine ecosystem, using the plankton functional type approach. It bases its description of primary production on a physiological model for phytoplankton growth that takes into account the nutrient availability effects on photoacclimation (Geider et al., 1998) and, for diatoms, on the relative frustule weight (Hohn, 2009). One specificity of REcoM is the representation of flexible stoichiometry, which leads to a description of elemental fluxes that can deviate from the fixed Redfield ratios often used in models (Redfield et al., 1963).

Here, we document the ocean biogeochemistry in the Regulated Ecosystem Model version 3 (REcoM3), coupled to the ocean and sea ice model FESOM2.1, and assess its performance in reproducing carbon and nutrient biogeochemical fluxes and the distribution of phytoplankton and zooplankton. Our aim is to analyse the new set-up regarding the coupled model state under historical atmospheric CO₂ forcing and the associated model bias and drift from the experiment with a constant preindustrial (PI) CO₂ level. We thus focus on evaluating model aspects with regard to the effects of climate change and CO₂ increase on carbon fluxes on century-scale timescales. We exclude in our analysis the deep-sea distribution of carbon and nutrients. This would require model runs over at least 500 to 2000 years (Séférian et al., 2020), which will be done in follow-up work.

2 Methods

2.1 Model description

We present the coupled ocean–sea ice–biogeochemistry model FESOM2.1–REcoM3. The previous model version (FESOM1.4–REcoM2) has been described by Schourup-Kristensen et al. (2014). Unlike its predecessor FESOM1.4, which uses a finite-element formulation, the ocean model is now based on a finite-volume discretization, which makes tracer conservation much easier to achieve. FESOM2.1 was described by Danilov et al. (2017) and evaluated in Scholz et al. (2019, 2022). The ocean biogeochemistry is simulated by the Regulated Ecosystem Model version 3 (REcoM3), which builds upon the previous version of REcoM2 (Hauck et al., 2013; Schourup-Kristensen et al., 2014). The advection and diffusion of 28 passive biogeochemical tracers is handled by FESOM2.1, whereas REcoM3 calculates the sources and sinks driven by biological interactions or biogeochemical exchange processes.

2.1.1 Ocean model FESOM2.1

FESOM2.1 solves the hydrostatic primitive equations under the Boussinesq approximation (Danilov et al., 2017). When set in a differential form, this equation is discretized on a finite set of points (nodes). As a first step of the mesh generation, a 2-dimensional grid is created by combining these nodes into triangular shapes (elements). At this stage, the mesh resolution (i.e. the size of triangles) can be adjusted in areas of interest without requiring a nesting approach. A 3-dimensional mesh is produced by projecting the triangles in a vertical direction and forming prisms. The scalar quantities (tracers and pressure) are located at the nodes, while the horizontal velocities are defined at centroids of the elements (see Figs. 1 and 2 in Danilov et al., 2017). A pair of control volumes are defined, where the vector control volumes are the prisms based on elements, and the scalar control volumes are formed by connecting cell centroids with edge midpoints (Fig. 1). Integration is carried out on a staggered Arakawa B-type mesh (Scholz et al., 2019).

Figure 1 Scheme of the cell–vertex discretization in 3-dimensional space. Blue dots correspond to scalar quantities, including REcoM3 state variables, located at the mid-layer vertices of triangles. Red dots represent horizontal velocities located at the mid-layer cell centres of the triangles. Yellow dots depict the vertical transfer velocities, which are placed at the layer boundaries aligned with the scalar quantities in the vertical.

We use FESOM2.1, an updated version of FESOM2.0. The updated model features include several developments, such as parallel and asynchronous output writing. An important new feature that we applied is the kinematic backscatter parameterization. This method takes into account the scales at which energy is scattered back to the resolved flow by introducing a negative viscosity term (Juricke et al., 2020). This greatly improves the simulation of eddy effects in coarse-resolution mesh set-ups (Juricke et al., 2020). The model code also includes the representation of ice shelf cavities (Timmermann et al., 2012), which has been used in regional studies with FESOM1.4–REcoM2 (Nissen et al., 2022). Ice shelf cavities are, however, not used in this study. Isoneutral tracer diffusion (Redi, 1982) and the Gent-McWilliams (GM; Gent and McWilliams, 1990; Griffies, 1998) eddy stirring parameterization are applied. Both GM and the Redi (1982) scheme are scaled with horizontal resolution, with a maximum value of 2000 m² s⁻¹ at 100 km horizontal resolution. The scaling decreases linearly below a resolution of 40 km to reach 0 at 30 km resolution, thus effectively switching the parameterization off. As a vertical mixing parameterization, the K-profile scheme is used (KPP; Large et al., 1994), with a background vertical diffusivity of 1 × 10⁻⁴ m² s⁻¹ for momentum and 1 × 10⁻⁵ m² s⁻¹ for tracers. Furthermore, the Monin–Obukhov length-dependent vertical mixing parameterization is applied in the surface boundary layer south of 50^∘ S (Timmermann and Beckmann, 2004).

Regarding the vertical discretization, FESOM2.1 is formulated with an arbitrary Lagrangian–Eulerian (ALE) scheme, which is a synthesis of different types of vertical coordinates. In the model configuration used here, we apply a full free-surface formulation and thus permit the vertical movement of the surface and of all other layers (referred to as zstar; Scholz et al., 2019). This drastically improves the tracer conservation properties (Campin et al., 2004). Partially filled cells are used at the ocean floor, resulting in a smoother representation of the bathymetry.

The sea ice component (Finite-Element Sea Ice Model, FESIM version 2) solves for sea ice concentration, ice and snow thickness, and ice drift velocity (Danilov et al., 2015). It is discretized on the same unstructured horizontal mesh as the ocean model. The elastic viscous plastic solver and flux-corrected transport scheme are used for sea ice advection (Danilov et al., 2015). The formulation of sea ice thermodynamics follows the work of Timmermann et al. (2009).

2.1.2 Biogeochemistry model REcoM3

REcoM3 is a water column biogeochemistry and ecosystem model, which incorporates cycles of carbon and nutrients (nitrogen, iron, and silicon) with varying intracellular stoichiometry in phytoplankton, zooplankton, and detritus (see the Appendix for detailed description and equations). Starting from the work by Schartau et al. (2007), REcoM was first used to describe carbon overconsumption in mesocosm experiments. After being coupled to the ocean and sea ice model MITgcm (Marshall et al., 1997), the previous version (REcoM2) with two phytoplankton classes, one zooplankton and one detritus class, was applied to study the cycling of marine carbon on present-day (Hauck et al., 2013, 2018) and glacial timescales (Du et al., 2022; Völker and Köhler, 2013), as well as the marine iron cycle (e.g. Völker and Tagliabue, 2015; Tagliabue et al., 2016; Ye and Völker, 2017; Pagnone et al., 2019). Moreover, REcoM2 was employed in assessments on the efficiency of ocean alkalinity enhancement (Köhler et al., 2013; Hauck et al., 2016), in data assimilation studies (Pradhan et al., 2019), and as a test bed for model development; e.g. for the development of a parameterization of iron-ligand binding based on pH (Ye et al., 2020), among others. Simultaneously, REcoM2 was coupled to FESOM1.4 (Schourup-Kristensen et al., 2014). These coupled model set-ups were used either in a global configuration (e.g. Schourup-Kristensen et al., 2014; Hauck et al., 2020), with a regional focus on the Arctic or the Antarctic (Hauck et al., 2015; Schourup-Kristensen et al., 2018; Oziel et al., 2022; Nissen et al., 2022) or in regional configurations (Taylor et al., 2013; Losch et al., 2014). Recently, the model has matured to include two groups of each classes of phytoplankton, zooplankton, and detritus (REcoM3; Fig. 2).

Figure 2 Schematic diagram of the components and interactions in the REcoM3 model. Small phytoplankton (phy) and diatoms (dia) take up inorganic nutrients (nut). Small zooplankton (small zoo.) and macrozooplankton (macrozoo.) consume phytoplankton and particles. Macrozooplankton feed on small zooplankton. Phytoplankton aggregation, zooplankton sloppy feeding, mortality, and fecal pellets generate sinking detritus (slow and fast). Sinking detritus degrades to dissolved organic carbon and nitrogen. Dissolved organic material (DOM) then remineralizes to dissolved inorganic carbon and nitrogen. The number of tracers and detailed processes are shown in the Appendix (Fig. A2).

Marine primary production is computed through representation of two phytoplankton functional types (PFTs), namely diatoms and small phytoplankton. The diverse group of small phytoplankton comprises a wide range of taxa, including, for instance, non-silicifying, calcifying, and non-calcifying haptophytes and green algae. The model allows PFTs to adapt their internal stoichiometry ($\mathrm{C}:\mathrm{N}:\mathrm{Chl}:{\mathrm{CaCO}}_{\mathrm{3}}$ ratios for small phytoplankton and $\mathrm{C}:\mathrm{N}:\mathrm{Chl}:\mathrm{Si}$ for diatoms) to nutrient levels, ambient light, and temperature, based on the photoacclimation model by Geider et al. (1998). Si uptake by diatoms is regulated as well, based on the internal Si:N quota, following Hohn (2009). This parameterization takes into account the strong decoupling between Si and N metabolism (e.g. Claquin et al., 2002) and prescribes the observed change in Si:N ratios under Fe and N limitation. The intracellular iron pool is derived from intracellular nitrogen with a fixed Fe:N ratio, based on the fact that intracellular iron is mostly associated with the photosynthetic electron transport chain and nitrogen metabolism (Geider and La Roche, 1994; Raven, 1988). REcoM3 also includes the photodamage parameterization by Álvarez et al. (2018). Calcium carbonate production is assumed to be linearly dependent on the gross small phytoplankton production. CaCO₃ dissolution is described by a depth-dependent dissolution rate.

Zooplankton is represented by two groups, namely small zooplankton and polar macrozooplankton (Karakuş et al., 2021), and each group has a carbon and nitrogen tracer. The small zooplankton group in the model is associated with relatively higher grazing rates compared to macrozooplankton and is widely spread in the global ocean. The polar macrozooplankton is mainly present in the Southern Ocean and northern high latitudes. The respiration rate is described mechanistically for macrozooplankton, taking into account the reduced metabolism in winter and increased metabolism at high grazing rates (Karakuş et al., 2021). For small zooplankton, respiration is calculated with a fixed respiration rate constant and biomass, which is in contrast to the previous version REcoM2, where respiration was used to drive zooplankton C:N back towards the Redfield ratio (Hauck et al., 2013; Schourup-Kristensen et al., 2014). Grazing is computed by applying a sigmoidal function, with variable preferences on both phytoplankton and detritus (Fasham et al., 1990).

Particulate organic matter (detritus) is split into two groups. The sinking speed of the first detritus group increases linearly with depth (from 20 m d⁻¹ at the surface to 192 m d⁻¹ at 6000 m depth; Kriest and Oschlies, 2008). The sinking speed of the second group (fast-sinking detritus) is constant throughout the water column (200 m d⁻¹; Karakuş et al., 2021). Remineralization of carbon and nitrogen occurs in two steps. Detrital material is first degraded to dissolved organic matter and then remineralized to the inorganic forms (dissolved inorganic carbon and nitrogen). For iron, it is assumed that the organic form is directly bioavailable, so it enters the dissolved iron pool in one step.

REcoM3 comprises a single-layer sediment pool for nitrogen, silicon, dissolved inorganic carbon, and calcium carbonate. The sinking detritus and associated minerals are accumulated in this layer when they reach the ocean floor. This material is subsequently returned back to the water column to the pools of dissolved inorganic nitrogen, carbon, and silicon, as well as alkalinity, with a fixed remineralization rate. The release of iron to the bottom layer of the ocean is assumed to be proportional to the release of inorganic nitrogen (Elrod et al., 2004).

2.1.3 Updates to previous REcoM version coupled to FESOM1.4

There are numerous improvements relative to the previously documented version of FESOM1.4–REcoM2 (Schourup-Kristensen et al., 2014), and the main changes are listed below.

REcoM

1. The routines for calculating carbonate chemistry and air–sea CO₂ exchange used in FESOM1.4–REcoM2, which followed the guidelines provided by the Ocean Carbon-Cycle Model Intercomparison Project (Orr, 1999), were replaced by the mocsy 2.0 scheme of Orr and Epitalon (2015). While both use the same thermodynamic equilibrium to calculate surface pCO₂ and CO₂ flux, mocsy 2.0 uses the faster and more accurate algorithm SolveSAPHE (Munhoven, 2013). Among other differences, it follows best-practice guides and uses recommended equilibrium constants. The gas exchange formulation is updated to that of Wanninkhof (2014), which is largely equivalent to that of Ho et al. (2006). The computed fluxes are scaled with the ice-free area.

2. Dissolved oxygen was added as a new tracer in REcoM3. The air–sea O₂ flux is calculated using the mocsy 2.0 routines (Orr and Epitalon, 2015). Photosynthesis, respiration, and remineralization change oxygen with a fixed O₂:C ratio, and remineralization does not depend on the O₂ levels in the current model version.

3. A second zooplankton group and a fast-sinking detritus class were added. The second zooplankton group represents a slow-growing polar macrozooplankton with a feeding preference on diatoms, which produces fast-sinking and carbon-rich fecal pellets (Karakuş et al., 2021).

4. The intracellular iron concentration is connected to intracellular nitrogen via a constant ratio Fe:N, leading to some variation in the Fe:C ratio, as briefly presented in Tagliabue et al. (2016) and Pagnone et al. (2019).

5. A sedimentary release of iron was added to the model (Tagliabue et al., 2016); this was done in addition to the previously considered Fe input with dust deposition.

FESOM

Biogeochemical fluxes returned back to the ocean from the benthos are treated with a specific bottom boundary condition. Variable bottom topography leads to a smaller scalar control volume located at the lowermost level. This is because the scalar control volumes are obtained by connecting the areas from the elements to which they are attached (see Fig. 1 in Danilov et al., 2017). Therefore, the number of elements around a single surface node may vary with depths when it meets non-flat topography. We thus computed the control volume and associated fluxes for each node by considering all surrounding elements at different depth levels.

Forcing

Our simulation was forced by the atmospheric reanalysis JRA55-do data set (Tsujino et al., 2018) instead of the CORE-II data set (Large and Yeager, 2009) that was used in previous assessments (Schourup-Kristensen et al., 2014). JRA55-do is a blend of reanalysis data and satellite observations and has the advantage to provide regularly updated near-real-time data up to the present day, with a higher temporal (3 h) resolution. The freshwater supplied by rivers is a climatology and provided by Large and Yeager (2004) as part of the CORE-II forcing. Nutrient, carbon and alkalinity supply via river discharge is not included in the experiments described here.

2.2 Experimental set-up and data

In this study, we used a mesh with a nominal resolution of 1^∘ as a background. The horizontal resolution is enhanced on the equatorial belt and in the region north of 50^∘ N to match $\mathrm{1}/\mathrm{3}$^∘ and 25 km, respectively. The mesh has 48 unevenly spaced vertical layers, where the layer thickness ranges from 5 m at the surface to 250 m in the deep ocean (Scholz et al., 2019).

Initial fields for temperature and salinity were taken from the winter statistical fields of the Polar science center Hydrographic Climatology (PHC3.0; updated from Steele et al., 2001) that ingests observations from the period 1900–1994. Total alkalinity (Alk) and preindustrial dissolved inorganic carbon (DIC) were initialized from version 2 of the Global Ocean Data Analysis Project (GLODAPv2) climatology centred on the year 2002 (Lauvset et al., 2016) and based on data collected between 1972 and 2013. Dissolved inorganic nitrogen (DIN) and dissolved silicic acid (DSi) were started with values from the World Ocean Atlas climatology of 2013 (Garcia et al., 2014) that occupied the period between 1955 and 2012. We used the World Ocean Atlas climatology of 2018 for dissolved oxygen (Garcia et al., 2019a, see Table 2), based on data for the time span 1955–2017.

Due to scarcity of observations, the iron field (DFe) was initialized with output from the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) model (Aumont et al., 2003), which was corrected using observed profiles for the Southern Ocean (de Baar et al., 1999; Boye et al., 2001). Sensitivity tests indicated that high values stemming from a hydrothermal vent in the eastern equatorial Pacific led to unreasonably large values in the interior Pacific Ocean due to advective fluxes. Therefore, the region spanning the latitudes of 9.5^∘ N–12.5^∘ S and longitudes 72–106^∘ W was masked to a maximum value of 0.3 µmol m⁻³ (below 2000 m). All other tracers were initialized with small values.

Iron was supplied to the ocean by dust deposition and from sediments. The sedimentary flux was assumed to scale with the organic nitrogen flux into the sediment, as found in Elrod et al. (2004). REcoM3 used monthly averages of dust deposition (Albani et al., 2014). We assumed that 3.5 % of the dust field consists of iron, of which 1.5 % dissolves into a bio-available form when deposited on the ocean surface. We did not include aeolian nitrogen deposition in our simulations.

Table 1 List of simulations performed in this study.

The atmospheric reanalysis data sets of JRA55‐do v.1.5.0 (Tsujino et al., 2018) were used to force the model for the period 1958–2021 (hereafter JRA55-do). A single repeating annual cycle of all forcing fields (year 1961) was used to perform the spin-up simulations and a control experiment. This is referred to as repeat-year forcing (hereafter called RYF61). We have deliberately chosen the year 1961, as it had rather neutral El Niño–Southern Oscillation conditions and also contained a low amount of anthropogenic perturbation, compared to the years of 1990 and 1991 recommended by Stewart et al. (2020).

A series of experiments were carried out in a global set-up to investigate the performance of the coupled FESOM2.1–REcoM3 model. The experiments follow the definitions used in the Global Carbon Budget (Friedlingstein et al., 2022a) and in the RECCAP (REgional Carbon Cycle Assessment and Processes, https://reccap2-ocean.github.io/, last access: 14 August 2023) projects and are summarized in Table 1. Our first experiment was forced with varying climate from the JRA55-do data set and varying atmospheric CO₂ levels (hereafter referred to as A). Atmospheric CO₂ mixing ratio (xCO₂) values are taken from the Global Carbon Budget (Joos and Spahni, 2008; Ballantyne et al., 2012; Friedlingstein et al., 2022a). A second simulation was forced by RYF61 atmospheric reanalysis fields and a preindustrial atmospheric CO₂ mixing ratio of 278 ppm. This configuration, here termed B, is considered to be the control run. Our last simulation was forced with varying climate from the JRA55-do data set and a preindustrial atmospheric CO₂ mixing ratio of 278 ppm. This experiment is referred as D and is used to separate the effects of rising atmospheric CO₂ and of climate change on the DIC inventory. Using the simulations A and B, the global ocean anthropogenic CO₂ sink was estimated by taking the model biases and drift from the control run into account. We used a coupled system spin-up (i.e. a direct strategy; Séférian et al., 2016). Before starting simulations A, B, and D, we performed spin-up experiments in two stages. In the first stage, a 189-year long (equivalent to three cycles of JRA55-do forcing) preindustrial spin-up simulation (named the pre-spin-up) was conducted using RYF61 atmospheric forcing and a preindustrial atmospheric CO₂ mixing ratio of 278 ppm until the air–sea CO₂ reached a quasi-equilibrium state. The A_spinup and B_spinup simulations are a continuation of the pre‐spin-up simulation, with either increasing (A_spinup) or constant (B_spinup) atmospheric CO₂, and run from 1800–1957. From the spin-up simulations, A, B, and D were branched off in 1958 and run until the end of 2021. FESOM1.4–REcoM2 and FESOM2.1–REcoM3 reach a throughput of 6 simulated years per day (SYPD) and 16 SYPD using the same mesh configuration and the same experimental set-up (see Table 1) on 288 cores with time steps of 15 and 45 min, respectively. All modelled mean fields shown in this work are averaged over the period 2012–2021, unless stated otherwise.

(Lauvset et al., 2016)(Lauvset et al., 2016)(Garcia et al., 2014)(Garcia et al., 2014)(Garcia et al., 2019b)(Sathyendranath et al., 2019)(Johnson et al., 2013)(Westberry et al., 2008)(Behrenfeld and Falkowski, 1997)(Buitenhuis et al., 2010; Moriarty and O'Brien, 2013)(Moriarty et al., 2013)Sallée et al. (2021a)(Bakker et al., 2016)(Chau et al., 2022)Chau et al. (2022)(Friedlingstein et al., 2022b; Wright et al., 2021; Schwinger et al., 2016; Lacroix et al., 2021; Berthet et al., 2019; Hauck et al., 2020; Liao et al., 2020; Doney et al., 2009; Aumont et al., 2015; Nakano et al., 2011; Urakawa et al., 2020; Long et al., 2021a; Landschützer et al., 2016; Rödenbeck et al., 2022; Chau et al., 2022; Gloege et al., 2021; Zeng et al., 2014; Iida et al., 2015; Gregor and Gruber, 2021)

Table 2 List of the observational data sets used to initialize the biogeochemistry model and assess its performance.

3 Results and discussion

In this section, we assess the performance of FESOM2.1–REcoM3 in simulating the observed mean state of nutrients, chlorophyll a, net primary production, and export production in the near-surface ocean, as well as air–sea CO₂ flux primarily under elevating CO₂. Before assessing the biogeochemical variables, we analyse the key features of the ocean model.

3.1 Modelled hydrography, mixed layer, and Atlantic meridional overturning circulation

An extended analysis of analogous FESOM2.1 runs is presented in Scholz et al. (2019, 2022). Here we analyse only a few relevant diagnostics to prove the validity of the presented research. We start the analysis by inspecting the spatial distribution of the model bias in the surface hydrography, which is presented in Fig. 3 as the difference between modelled mean 2012–2021 and the PHC3.0 Climatology (Steele et al., 2001). For temperature and salinity, respectively, we found a global spatial correlation coefficient (r) of 0.99 and 0.99, with a root mean squared error (RMSE) of 0.82 ^∘C and 0.43 psu. In the northern North Atlantic, the bias is expressed by cold (∼ 4 ^∘C colder) and fresh (∼ 1 psu fresher) anomalies around Newfoundland, which is a typical bias for standalone and climate models at coarse resolutions (see, e.g., Scaife et al., 2011). Further south, the bias depicts a dipole anomaly associated with the Gulf Stream going too far north, which is a commonly addressed shortcoming for non-eddy-permitting models (see, e.g., Zhang and Vallis, 2007; Storkey et al., 2018). Similar issues are found in comparable current systems, such as the Kuroshio and Malvina systems. It is, however, surprising that in general FESOM is far too saline at the surface and is on average 0.3 psu saltier than the climatology. The reason for this bias could be imperfections in the river discharge from CORE-II forcing and the relatively low surface salinity restoring, which uses a piston velocity of 50 m/300 d in the simulations. In most of the ocean, the sea surface temperature (SST) and sea surface salinity (SSS) differences act in an opposite manner on buoyancy. Hence, the increase or decrease in SST is accompanied by an increase or decrease in SSS. The only exception is the Indian Ocean, where east and west in simulation A become less and more buoyant, respectively (Fig. 3).

Figure 3 Maps of simulated FESOM2.1–REcoM3 (simulation A) surface temperature (^∘C) (a), practical salinity (unitless) (d), with observations from the Polar science center Hydrographic Climatology (PHC3.0; updated from Steele et al., 2001), and corresponding differences (c, f), averaged over the time period 2012–2021.

Figure 4 Vertical representation of the Atlantic meridional overturning circulation (Sv) in simulations A, B, and their difference (Sv).

Figure 5 Time series of the annual mean Atlantic meridional overturning circulation (Sv) maxima in simulations A and B.

Figure 6 Maps of simulated FESOM2.1–REcoM3 (simulation A) maximum mixed layer depth (m) in March (a) and September (d), averaged over the time period 2012–2021, using the observation-based maximum mixed layer depth from Sallée et al. (2021a, b and e), which occupied the period between 1970 and 2018, and the corresponding differences (c, f).

In Fig. 4, we augment the diagnostic by inspecting the Atlantic meridional overturning circulation (AMOC), which provides the most general characteristic of water mass transformation and production. The mean AMOC in both runs is expressed with the basin-wide mid-depth cell, showing a maximum of ∼ 15 Sv at ca 40^∘ N. The bottom cell, induced by the circulation of the Antarctic bottom water, is also well reproduced, with a minimum of $\sim -\mathrm{5}$ Sv. Even though the runs depict large differences in temperature and salinity from the observed climatology, the simulated AMOC shows the canonical picture known from other standalone ocean and coupled climate models (Griffies et al., 2009; Jungclaus et al., 2013; Danabasoglu et al., 2014). This indicates that although biases in the representation of water mass properties and ventilation mechanisms are present, they still result in a reasonable density distribution which maintains realistic transport.

The difference between simulations A and B shows that the mid-depth and bottom cells are stronger in simulation B. Consequently, the difference A − B is expressed by a basin-wide positive anomaly, with a maximum of ∼ 3 Sv. We also show the time series, for both runs, of AMOC maxima for the years 1958–2021 (Fig. 5). In simulation A, the time series depicts a multidecadal variability, with a minimum of ∼ 9.5 Sv and a maximum of ∼ 13.5 Sv. Concurrently, the reference simulation B depicts a nearly constant value, with a small increase between 9.5 and 10 Sv, which is a result of the repeat-year forcing.

Finally, in Fig. 6, we present the simulated and observed (Sallée et al., 2021a; referred to as Atlas) annual maximum mixed layer depth (MLD) pattern for March and September, following the same methods (the depth at which the potential density referenced to the surface exceeds the density of the water by a threshold of 0.03 kg m⁻³). Overall, the modelled MLD fits well with the observations, although some common discrepancies remained in the deep mixing areas. In the Northern Hemisphere, the deepest MLD (>1000 m) is found in the Labrador and Greenland–Iceland–Norwegian seas. The magnitude is larger than in Sallée et al. (2021a) but is in the same range as in other modelling studies (Griffies et al., 2009; Sidorenko et al., 2011). In the Southern Hemisphere, winter deep mixing in high latitudes is also overestimated compared to the observations, especially in the Pacific sector of the Southern Ocean. From inspecting the model runs and their differences, we conclude that FESOM2.1 simulated a reasonable ocean state which can be used for further analysis.

3.2 Nutrients, ocean productivity, and ecosystem

3.2.1 Modelled versus in situ nutrients

We first compared the spatial distribution of the surface (averaged over the top 100 m depth layer) ocean dissolved inorganic nitrogen (DIN) and dissolved silicate (DSi) from REcoM3 with the World Ocean Atlas 2018 (Garcia et al., 2019b) climatologies (Fig. 7). While simulated surface DIN concentrations were lower than observations in the subpolar regions, a large positive DIN bias of up to 20 mmol m⁻³ was found in the subtropical South Pacific Ocean. The simulated DSi was overestimated in the Southern Ocean and underestimated in the northern Pacific. Exceptions are the Pacific and Atlantic sectors of the coastal Southern Ocean, where the modelled DSi concentrations are lower than the observations. These patterns were already present in FESOM1.4–REcoM2 (Schourup-Kristensen et al., 2014); however, two recent improvements should be noted. First, the large and positive DIN bias in the northern subtropical Pacific (Schourup-Kristensen et al., 2014) disappeared. This is caused by replacing the dust deposition input forcing field from Mahowald et al. (2003) with Albani et al. (2014), which results in a more realistic (i.e. less strong) iron limitation. Second, the silicate bias in the Southern Ocean is reduced in magnitude and extent compared to Schourup-Kristensen et al. (2014). This is related to tuning experiments (not shown), which resulted in a larger share of diatoms in the Southern Ocean (Fig. 11) compared to Schourup-Kristensen et al. (2014), thus drawing down more silicic acid. Along with the increased share of diatoms, the Southern Ocean and global opal export also increased. For the global ocean, the opal export increased from 74.5 Tmol Si yr⁻¹ in Schourup-Kristensen et al. (2014) to 168 Tmol Si yr⁻¹ in the present study and is thus at the upper end instead of the lower end of the range of 69–185 Tmol Si yr⁻¹ (Dunne et al., 2007) and higher than the best estimate of Tréguer et al. (2021, Table 3). In the Southern Ocean, opal export increased from 21.5 Tmol Si yr⁻¹ in Schourup-Kristensen et al. (2014) to 85.5 Tmol Si yr⁻¹, which is higher than the range of 21–54 Tmol Si yr⁻¹ reported by Dunne et al. (2007). The silicic acid bias is rather insensitive to the formulation and parameter choice of opal dissolution but very sensitive to the share of diatoms in the Southern Ocean. The correlation coefficient (r) and root mean squared error (RMSE) between the simulated and observed annual mean were 0.88 and 0.86 mmol m⁻³ for DIN and 0.47 and 0.54 mmol m⁻³ for DSi. The correlation with observed DIN is higher than in Schourup-Kristensen et al. (2014, 0.75), which we relate to the disappearance of the DIN bias in the northern subtropical Pacific. The correlation with observed DSi is lower than in FESOM-1.4–REcoM2, despite the reduction in magnitude and extent in the Southern Ocean DSi bias. Moderately high silicic acid values in the northern high latitudes are not reproduced. This may be related to mixing that is too sluggish or to overly strong silicic acid drawdown by diatoms (Figs. 9 and A1), which is possibly linked to an iron limitation that may be too weak (Fig. 12).

Figure 7 Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (0–100 m) concentration of dissolved inorganic nitrogen (mmol m⁻³; a), dissolved inorganic silicon (mmol m⁻³; d), with observations from the World Ocean Atlas 2018 climatology (b and e; Garcia et al., 2019b), and the corresponding differences (c, f) averaged over the time period 2012–2021.

Figure 8 Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (0–50 m) concentration of dissolved iron (µmol m⁻³; a), and of the AI-based global reconstruction by Huang et al. (2022, b). Note the different colour scales for the two plots.

Despite the enormous increase in the number of observations of dissolved iron with the GEOTRACES project, observations have not reached a global coverage that makes it possible to construct a global climatology. Therefore, the modelled dissolved iron is compared here to the global surface pattern of dissolved iron by Huang et al. (2022), which uses an artificial intelligence (AI) method (random forest) to construct a near-global iron field, based on the observations in the second intermediate GEOTRACES data product (Schlitzer et al., 2018), plus some older in situ iron observations compiled in Tagliabue et al. (2012), and on co-located hydrographic observations. The pattern of modelled dissolved iron (Fig. 8; averaged over the top 50 m) shows the expected pattern of high concentrations in regions with high dust deposition, mainly in the tropical Atlantic Ocean and the eastern part of the Arabian Sea and also to some extent in the southern subtropical Atlantic and Indian oceans. Concentrations are extremely low in the subpolar Southern Ocean and in almost the whole equatorial and South Pacific. Iron concentrations are also low in the subpolar North Pacific, and – less so, but still noticeable – in the subpolar North Atlantic. Oceanic regions adjacent to extended shelves, especially in the Arctic, show somewhat elevated iron concentrations. If we compare this to the AI-generated global pattern of dissolved iron from Huang et al. (2022), then we find qualitatively similar patterns, like the elevated iron concentration in the equatorial and subtropical Atlantic and the Arabian Sea or the low concentrations in the subpolar Southern Ocean, the equatorial Pacific, and the subpolar North Pacific, but the amplitude of the patterns is smaller overall. The largest discrepancy in amplitude is found under the Saharan dust plume in the tropical Atlantic, where the model produces maximum dissolved iron values that are almost 3 times as high as the reconstruction from Huang et al. (2022). Direct observations in the tropical Atlantic also show dissolved iron concentrations that reach 1.2 nmol L⁻¹ (e.g. Hatta et al., 2015), while modelled maxima are >3 nmol L⁻¹. A further important difference is that the distribution by Huang et al. (2022) shows slightly elevated iron concentrations in the centre of the subtropical South Pacific, where the model in contrast has extremely low values. This discrepancy causes an iron limitation that is too strong in this region in the model, probably explaining the overly high DIN concentrations in the model in the South Pacific. The fact that the amplitude of the patterns in modelled dissolved iron is too high, which is also found in other models, likely has a number of causes. The most important one is the assumption of a constant solubility in dust-deposited iron. Dust deposition close to the main source regions is on average coarser and has experienced less chemical processing during its transport, both of which would lead to a lower solubility. The opposite is true for particles deposited far from the source regions, such as in the South Pacific. A second contribution might be the missing source from pyrogenic aerosols, which are far more soluble. Also, the effect of dust particles as iron scavengers, which has not been included in this simulation, has been shown to reduce the overly high dissolved iron concentrations often found in models in the main dust deposition regions (Ye and Völker, 2017; Pagnone et al., 2019). Furthermore, the intensity and extension of dust plumes vary between modelled dust deposition fluxes (e.g. Myriokefalitakis et al., 2018). The field of dust deposition by Albani et al. (2014), used in our model to calculate aeolian iron input, is within the range of modern estimates but surely contains some uncertainties. Despite the overall amplitude of the patterns in the dissolved iron being too strong, especially in the regions of high dust deposition, the model is able to reproduce the main regions in which iron availability limits phytoplankton productivity (Moore et al., 2013), namely the subpolar Southern Ocean, the equatorial and North Pacific, and also to some extent the seasonal iron limitation in the subpolar North Atlantic (Nielsdóttir et al., 2009), but overestimates iron limitation in the subtropical South Pacific.

3.2.2 Modelled versus satellite-based phytoplankton biomass and productivity

We first evaluated the spatial distributions of the modelled chlorophyll a concentration (Fig. 9), which is an indicator for phytoplankton biomass, and vertically integrated the net primary production (NPP, Fig. 10). We compared chlorophyll a concentrations obtained from FESOM2.1–REcoM3 simulations, averaged from 2012 to 2021, with an ocean colour remote sensing merged data set (from the Ocean Colour Climate Change Initiative, OC-CCI; Sathyendranath et al., 2019), averaged from 1998 to 2019. We also compared the modelled NPP with satellite estimations, such as the Vertically Generalized Production Model (VGPM; Behrenfeld and Falkowski, 1997; Fig. 10) and the updated carbon-based productivity model (CbPM; Westberry et al., 2008; see Fig. A1 in the Appendix). VGPM is a chlorophyll-based algorithm that can be considered to be a standard NPP estimation from ocean colour for the last 20 years (Lee and Marra, 2022). CbPM uses spectrally resolved light attenuation and is based on a semi-analytical algorithm (Garver–Siegel–Maritorena, GSM; Maritorena et al., 2002). All data sets were also compared along a latitudinal distribution with an improved chlorophyll a algorithm for the Southern Ocean (Johnson et al., 2013; Fig. 11).

The results for chlorophyll a and NPP obtained here are comparable to those presented by Schourup-Kristensen et al. (2014). Over large parts of the global ocean, the mean surface chlorophyll a concentrations are in agreement with the observations (Fig. 9c and d). Yet, there are regional differences. In temperate latitudes, the modelled chlorophyll a concentrations are somewhat higher than observed, while the subtropical gyres show concentrations slightly lower than the observations. The comparison of modelled and observation-based satellite estimates of chlorophyll a yielded a correlation of 0.66 and an RMSE of 0.38 mg m⁻³. FESOM2.1–REcoM3 also shows a reasonably well-simulated latitudinal variation in chlorophyll a compared to satellite estimations (Fig. 11a). The model underestimates chlorophyll a concentrations in most of the coastal regions, especially in the high-latitude regions. In the southern high latitudes, FESOM2.1–REcoM3 follows the Southern-Ocean-adjusted chlorophyll a data set (Johnson et al., 2013) quite well, except for the coastal regions close to Antarctica (approximately south of 70^∘ S). In the Arctic Ocean, the model strongly underestimates the chlorophyll a concentrations, which is driven by negative biases reaching up to 3 mg chlorophyll a m⁻³ on the continental shelves (Fig. 9). Although FESOM2.1–REcoM3 did reproduce the NPP distribution at low latitudes well (Figs. 10 and 11), it also strongly underestimated the NPP at higher latitudes when compared to VGPM (r=0.43; RMSE = 0.34 mgC m⁻² d⁻¹), in particular in productive areas north of 50^∘ N and coastal areas (Fig. 10). For regional applications, further analysis and possibly tuning may be needed. When compared with VGPM, the model simulation generally underestimated the remotely sensed NPP estimations (Table 3), especially in the subtropical Pacific. Yet, with a value of 35.9 PgC yr⁻¹, the modelled global total NPP is slightly above the range of the earlier modelling studies (23.7–30.7 PgC yr⁻¹; Schneider et al., 2008) and within the range of recent Earth system models (24.5–57.3 PgC yr⁻¹; Séférian et al., 2020). It is lower than other satellite-based estimates of 47.3 PgC yr⁻¹ (Behrenfeld and Falkowski, 1997), 52 PgC yr⁻¹ (Westberry et al., 2008), and 48.7–52.5 PgC yr⁻¹ reported by Kulk et al. (2020).

Both simulated chlorophyll a concentrations and NPP from FESOM2.1–REcoM3 seemed to be underestimated in coastal regions. Primary production and chlorophyll a levels that are too low were particularly evident in coastal regions, which could be linked to deficiencies in either the chlorophyll a data set and/or in the model. For the former, the chlorophyll a OC-CCI data set and the CbPM primary production data set uses the GSM algorithm. GSM tries to distinguish the optical signatures from phytoplankton, particles, and dissolved organic matter but still requires regional tuning in coastal regions, where the presence of non-biotic optically active material (i.e. yellow substances and sediments) makes chlorophyll a retrieval challenging (Blondeau-Patissier et al., 2014). The overestimation of chlorophyll a in coastal waters is even more pronounced with the use of standard global chlorophyll algorithms in the VGPM primary production data set, such as OC4 that are only adapted to CASE-I waters (low influence of dissolved organic matter and non-algal particles). Therefore, both remotely sensed NPP estimations carry uncertainties related to the global algorithms. For example, turbid waters over the Arctic shelves are known to artificially increase both chlorophyll a and NPP estimates from remote sensing (Matsuoka et al., 2012; Mitchell, 1992; Mustapha et al., 2012). Some recent advances used local parameterizations with in situ data, which resulted in much lower productivity levels in those coastal areas (Lewis et al., 2020; Lewis and Arrigo, 2020). Generally, the NPP and chlorophyll differences compared to satellite-based estimates could also be linked to model deficiencies, such as a coarse model resolution and associated weak upwelling, missing complexity in the simulated phytoplankton classes, and also the so-far unconsidered nutrient input from terrigenous sources.

Figure 9 Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (log₁₀ transformed) chlorophyll a concentration (mg Chl m⁻³) of small phytoplankton (a), diatoms (b), and the sum of both phytoplankton groups (c). The satellite-based merged data set OC-CCI (Sathyendranath et al., 2019) is shown in panel (d), with the corresponding differences between FESOM2.1–REcoM3 and OC-CCI shown in panel (e). Note the different time periods for the simulation (2012–2021) and OC-CCI (1998–2019).

Figure 10 Maps of simulated FESOM2.1–REcoM3 (simulation A) vertically integrated net primary production (mgC m⁻² d⁻¹) of small phytoplankton (a), diatoms (b), and the sum of both phytoplankton groups (c). The satellite-based Vertically Generalized Production Model (VGPM; Behrenfeld and Falkowski, 1997) is shown in panel (d), with the corresponding differences between FESOM2.1–REcoM3 and VGPM shown in panel (e). All fields are averaged over the time period from 2012 to 2021.

The low values of primary production could be caused by several top-down and/or bottom-up effects. The nutrient dynamics that partly control NPP are the result of a delicate balance between physical (mixing, stratification, and upwelling systems) and biogeochemical processes. To investigate bottom-up controls on regional NPP dynamics, we derived the most limiting factor (either light or nutrients) of the growth of diatom and small phytoplankton. This factor ranges between zero (most limiting) and one (no limitation) and is based on the nutrient uptake Michaelis–Menten kinetics of REcoM. The Michaelis–Menten coefficient (MM) is computed as MM = [Nut] $/$ ([Nut] + KNut), with [Nut] being the nutrient concentration, and KNut a nutrient and phytoplankton-dependent half-saturation constant. The light limitation is defined as the carbon-specific photosynthesis rate divided by the maximum photosynthetic rate. We derived maps showing the most limiting factor (the factor closest to zero, with either the nutrients of DIN, DSi, or DFe or light) in the annual mean (Fig. 12).

Figure 11 Latitudinal distribution of vertically integrated and zonally averaged (a) chlorophyll a (mg Chl m⁻³) and (b) net primary production (mg C m⁻² d⁻¹) simulated by FESOM2.1–REcoM3 (blue line). The satellite-based merged chlorophyll a data sets of OC-CCI (Sathyendranath et al., 2019; orange line) and the improved chlorophyll a algorithm for the Southern Ocean (Johnson et al., 2013; green line) are shown in panel (a). The satellite-based data set of the Vertically Generalized Production Model (VGPM; Behrenfeld and Falkowski, 1997; purple line) and the carbon-based productivity model (CbPM; Westberry et al., 2008; red line) are shown in panel (b).

The spatial distribution of the dominant growth-limiting factor for diatoms and small phytoplankton over the time period 2012–2021 is shown in Fig. 12. Over large areas of the South Pacific and almost the entire Southern Ocean, diatoms were limited by iron availability. Elsewhere, except for the Arctic Ocean, where light was the most limiting factor, diatom growth was controlled by the abundance of dissolved silicic acid. Nutrient uptake of small phytoplankton was limited by iron in the South Pacific, DIN within the band of 45^∘ S–45^∘ N in the Atlantic and Indian oceans, and insufficient light at high latitudes (south of 45^∘ S and north of 45^∘ N).

The large-scale patterns of limitation were in general agreement with observations (Moore et al., 2013) and other modelling studies (Long et al., 2021a), although the degree of silicic acid limitation for diatoms (outside the iron-limited Southern Ocean) varied across models (Laufkötter et al., 2015). The severe iron limitation in most of the Pacific might contribute to the lower productivity levels than observed in the same regions (Fig. 10).

In addition to bottom-up explanations, the formulation and parameter choices for zooplankton grazing may be a reason for low primary production (Anderson et al., 2010; Prowe et al., 2012; Karakuş et al., 2021). In fact, Karakuş et al. (2022) demonstrated that a separation of the small zooplankton group in REcoM into micro- and mesozooplankton leads not only to a 25 % increase in NPP but also to a reduction in the overly strong iron limitation in the South Pacific, due to nutrient recycling by zooplankton. Furthermore, REcoM does not explicitly represent picophytoplankton (e.g. non-N₂-fixing cyanobacteria such as Synechococcus and Prochlorococcus) and nitrogen fixers, and this might contribute to an underestimation of NPP.

Figure 12 Maps showing the simulated spatial distribution of the most limiting factor in the surface water for diatoms (a) and small phytoplankton (b). Hatching denotes areas of weak limitation (all limiting factors >0.5). Note that Fe is for iron, DIN is for dissolved inorganic nitrogen, and DSi is for dissolved silicic acid.

3.2.3 Modelled versus MAREDAT zooplankton biomass

In REcoM3, the small zooplankton group is widely spread in the global ocean and the highest biomass occurs in high-productivity regions (Fig. 13a). The macrozooplankton is present in the high latitudes (Fig. 13b), since it is parameterized as a polar macrozooplankton group (Karakuş et al., 2021). We compare the latitudinal distribution of the integrated modelled zooplankton biomass with gridded global zooplankton biomass data from MAREDAT (Buitenhuis et al., 2010; Moriarty et al., 2013; Moriarty and O'Brien, 2013). The simulated biomass of small and total zooplankton reproduces MAREDAT-derived biomass reasonably well in low to midlatitudes but underestimates biomass in the polar regions (Fig. 13c). The underestimation of zooplankton biomass in the northern high latitudes may be related to an underestimation of primary production in the same region. In agreement with the MAREDAT data set (Moriarty et al., 2013), macrozooplankton is not present in low latitudes.

Figure 13 Maps of annual mean surface (a) small zooplankton and (b) macrozooplankton concentrations in FESOM2.1–REcoM3 (simulation A). Latitudinal distribution of vertically integrated zooplankton biomass (mg C m⁻²) of (c) modelled small zooplankton (solid blue line) and the sum of microzooplankton and mesozooplankton from MAREDAT (orange dots for individual observations and solid brown line for the zonal mean of the observations; Buitenhuis et al., 2010; Moriarty and O'Brien, 2013) and (d) modelled macrozooplankton (solid blue line) and macrozooplankton from MAREDAT (orange dots; Moriarty et al., 2013). The modelled zooplankton biomass is averaged over the time period from 2012 to 2021. The zonal mean of macrozooplankton is not calculated due to the low number of observations.

3.2.4 Synthesis

The modelled biogeochemical fluxes were compared to the previous version FESOM1.4–REcoM2 and observational studies (Table 3). Modelled global NPP is higher in FESOM2.1–REcoM3 than in FESOM1.4–REcoM2 but still lower than in satellite-based estimates. The estimate is comparable to other global modelling studies (Schneider et al., 2008; Séférian et al., 2020). Export production (EP) is slightly higher in FESOM2.1–REcoM3 than in the previous version and falls within the observational range previously documented in the literature for both the global ocean and the Southern Ocean. For the global ocean, FESOM2.1–REcoM3 NPP and EP estimations remained at the lower end of the range, despite a slight increase in NPP. A more detailed description of zooplankton results in more efficient nutrient recycling and can thus increase NPP by 25 % (see also explanation in Sect. 3.2.2; Karakuş et al., 2022). In the Southern Ocean, estimations of NPP and EP remained very close to observation-based estimates. Maybe the most noticeable change between the two model versions is the substantial increase in opal export, which increased by a factor of 4 in the Southern Ocean, passing from the lower to the higher end of the observational range of an earlier review (Dunne et al., 2007), and lies above an updated estimate (Tréguer et al., 2021). This is due to an increase in the relative contribution of diatoms to the total NPP in high latitudes (Fig. 11).

(Behrenfeld and Falkowski, 1997)(Westberry et al., 2008)(Schneider et al., 2008)(Kulk et al., 2020)Séférian et al., 2020(Schlitzer, 2004)(Dunne et al., 2007)(Henson et al., 2011)(Siegel et al., 2014)Dunne et al., 2007(Tréguer et al., 2021)Lee, 2001; Jin et al., 2006Gangstø et al., 2008; Berelson et al., 2007; Dunne et al., 2007; Battaglia et al., 2016; Gehlen et al., 2006(Carr et al., 2006)(Behrenfeld and Falkowski, 1997)(Schlitzer, 2002; Nevison et al., 2012)(Dunne et al., 2007)(Dunne et al., 2007)

Table 3 Global and Southern Ocean net primary production (NPP) and export production (EP) in FESOM2.1–REcoM3 and estimates from the literature. The Southern Ocean is here considered to be the region south of 50^∘ S. The numbers for VGPM and CbPM are recalculated after interpolation to the model mesh over the years 2012–2019.

3.3 Carbon cycle

3.3.1 Dissolved inorganic carbon and alkalinity

Insight into the carbonate system can be obtained by inspecting surface maps of modelled dissolved inorganic carbon and alkalinity and the corresponding observational Global Ocean Data Analysis Project (GLODAPv2) climatologies (Fig. 14). Global patterns of simulated concentrations resemble the observed fields reasonably well (r=0.99; RMSE = 36.5 mmol m⁻³; calculated from annual means), with the highest DIC values in the subtropical gyres of the Atlantic and south Pacific and the subpolar North Atlantic and the Southern Ocean. Similar to GLODAP, the highest alkalinity values are found in the subtropical gyres of the Atlantic and south Pacific, with a good agreement with global spatial features (r=0.99; RMSE = 33.9 mmol m⁻³). Yet, simulated surface DIC and alkalinity concentrations were slightly overestimated in most of the global ocean. Two major exceptions are the Arctic Ocean and the North Atlantic, where the concentrations were underestimated, and the tropical upwelling regions and the northern Indian Ocean for DIC. The bias patterns differ relative to FESOM-1.4–REcoM (which was too low for DIC and alkalinity in the tropics and subtropics and too high in the high latitudes; not shown), which indicates that different realizations in circulation or mixing may drive these bias patterns. This is in line with an overestimation of surface salinity in most of the global ocean, with the exception of the North Atlantic and the Arctic Ocean (see Fig. 3). Also, surface alkalinity biases are generally attributed to a dominant physical (preformed) signal, with a smaller contribution from the calcium carbonate cycle and a negligible contribution from organic matter remineralization (Koeve et al., 2014). However, tuning the model to result in a higher CaCO₃ production could possibly also counteract the positive alkalinity bias. Similarly, a higher NPP in the South Pacific could regionally ameliorate the high DIC bias. A positive bias in alkalinity at constant atmospheric CO₂ in the spin-up (not shown) and simulation A (Fig. 14) leads to a positive bias in DIC, as surface water with a higher alkalinity can hold more CO₂ in equilibrium than a low-alkalinity surface ocean. The range of biases is similar (±100 mmol m⁻³) to other ocean biogeochemical models (e.g. Tjiputra et al., 2020; Long et al., 2021a).

Figure 14 Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (0–100 m) concentration of dissolved inorganic carbon (mmol m⁻³; a) and alkalinity (mmol m⁻³; d). Corresponding data from GLODAPv2 (b, e; Lauvset et al., 2016) and model–data differences (c, f) are shown. The comparison is for the period 1998–2006, which is centred on the year 2002 to be comparable to GLODAP.

3.3.2 Surface ocean pCO₂ and air–sea CO₂ flux

We compare the pattern of the temporal mean (2012–2021) surface ocean partial pressure of CO₂ (pCO₂; Fig. 15) and air–sea CO₂ flux (Fig. 16) to the pCO₂-based data product of Chau et al. (2022), with a seamless coverage from open ocean to the coasts (Fig. 16). Different pCO₂ products largely agree with each other in terms of spatial patterns, although they differ with respect to amplitude and timing of variability in the regionally or globally integrated fluxes (Fay et al., 2021; Fay and McKinley, 2021). Therefore, we chose one product (Chau et al., 2022) and focus on the comparison of the spatial pattern with our model. We further evaluate the temporal evolution of pCO₂ in FESOM2.1–REcoM3 with a direct comparison to surface ocean pCO₂ observations from the Surface Ocean CO₂ Atlas (SOCAT; Bakker et al., 2016), where we subsampled the model output for spatiotemporal locations where observations exist, following Hauck et al. (2020) and Friedlingstein et al. (2022b) in Fig. 17.

Figure 15 Maps of surface ocean pCO₂ (µatm). The top row compares the (a) simulated FESOM2.1–REcoM3 (simulation A) surface partial pressure of CO₂ to the (b) pCO₂-based data product (Chau et al., 2022) and both are averaged over 2012–2020. Panel (c) shows model–data differences. Note that simulation and observations are masked every month with the sea ice concentration >15 %.

Figure 16 Maps of air–sea CO₂ fluxes (mol C m⁻² yr⁻¹). The top row compares the (a) simulated FESOM2.1–REcoM3 (simulation A) CO₂ flux to the (b) pCO₂-based data product (Chau et al., 2022) and both are averaged over 2012–2020. Panel (c) shows model–data differences. Negative numbers indicate a flux into the ocean. Note that simulation and observations are masked every month with the sea–ice concentration >15 %.

The large-scale spatial patterns of pCO₂ are well reproduced (Fig. 15) with high values in the tropics that are typically higher than atmospheric values (red colours) and lower values in the subpolar Southern and Pacific oceans and the high-latitude North Atlantic (r=0.75 and RMSE = 29.2 µatm; ice-free areas). However, compared to the pCO₂ product of Chau et al. (2022), the model pCO₂ values are overestimated in the subtropical gyres (Fig. 15c). Furthermore, the North Atlantic pCO₂ is on average lower than the pCO₂ product, and the two data sets also differ in terms of the over- versus undersaturation of pCO₂ relative to the atmosphere in the polar Southern Ocean (higher values in FESOM2.1–REcoM3). The latter may well be explained by a known summer bias in the Southern Ocean pCO₂ observations (e.g. Metzl et al., 2006; Gregor et al., 2019). FESOM2.1–REcoM3 also simulates very high pCO₂ values on the Russian shelves in the Arctic, where hardly any observations exist. Similarly high pCO₂ values were reported for this region by Anderson et al. (2009), but missing repeat observations prevent a conclusion on whether this is a robust signal and what its extent in time and space is.

Figure 17 Comparing the annual mean pCO₂ (µatm) from FESOM2.1–REcoM3 (simulation A; subsampled for spatiotemporal locations of observations in SOCAT; red) with observations from SOCATv2022 (light blue; updated from Bakker et al., 2016). Results are shown as spatially averaged for (a) the global ocean, (b) the north (>30^∘ N), (c) the tropics (30^∘ S–30^∘ N), and (d) the south (<30^∘ S). The time series are shown for all observations in SOCAT (since 1970), but the correlation coefficient r (unitless) and root mean squared error (RMSE; µatm) are indicated in the panels for the time period 1990–2021.

FESOM2.1–REcoM3 reproduced the temporal evolution of surface ocean pCO₂ reasonably well compared to SOCAT when accounting for where and when the pCO₂ sampling took place (Fig. 17). The annual correlation coefficient and RMSE between the simulated and observed global mean pCO₂ are 0.93 and 4.6 µatm, respectively. The subsampled model follows the SOCAT time series closely, including its variability due to a sampling distribution in space and time. The global mismatch with SOCAT pCO₂ as measured by the RMSE is comparable to or slightly below the value for FESOM-1.4–REcoM2 (see Fig. S9 in the Supplement of Hauck et al., 2020, for 1985–2018) and comparable to but at the high end of the range of other models in GCB 2022 (1990–2021; Friedlingstein et al., 2022b). On a monthly scale, the RMSE is higher (38 µatm), as the models capture the large-scale patterns better than smaller-scale variability, according to a previous assessment (Hauck et al., 2020). An analysis of large-scale regional patterns (north, tropics, and south; Fig. 17) reveals that the model reproduced the trend well but overestimates the mean pCO₂ in the tropics and underestimates pCO₂ in the northern extratropics and to a lesser extent in the southern extratropics in recent decades, as also indicated in the maps (Fig. 15).

The air–sea CO₂ flux spatial pattern was reasonably reproduced by FESOM2.1–REcoM3, with CO₂ uptake in the subpolar regions of both hemispheres and outgassing in the tropics and north Pacific (Fig. 16; r=0.72; RMSE = 1.45 mol C m⁻² yr⁻¹). Generally, the CO₂ flux patterns mirror the pCO₂ patterns (Fig. 15) but with the additional imprint of spatial variability in the wind speed. Hence, the CO₂ uptake in the subpolar Southern Ocean may appear large compared to pCO₂, which is not as strongly undersaturated in the South Atlantic as in the North Atlantic. Regions of mean outgassing in the Southern Ocean are seen to a smaller extent in the model than in the pCO₂ product. While it is well established that the outgassing of CO₂ in the polar Southern Ocean occurs in winter (e.g. Bakker et al., 1997), its magnitude and timing varies between estimates and is being debated (Gruber et al., 2009; Lenton et al., 2013; Gray et al., 2018; Bushinsky et al., 2019; Sutton et al., 2021; Long et al., 2021b). The misfit between the annual mean modelled CO₂ flux and the pCO₂-based data product generally mimics pCO₂ biases and thus shows small positive differences (less uptake or more outgassing) in the subtropical gyres and small negative differences (stronger uptake or less outgassing) in the equatorial Pacific and the Southern Ocean (Fig. 16; bottom panel). The strongest biases were found in the northern high latitudes (negative bias) and the upwelling zone of the eastern tropical Pacific (positive bias). FESOM2.1–REcoM3 also generally captures the large-scale patterns of coastal CO₂ fluxes with CO₂ uptake in the mid- and high latitudes (poleward of 25^∘ N/S) and outgassing in the tropical coastal ocean, as described in a recent synthesis based on low- and high-resolution models and pCO₂ products (Resplandy et al., 2023). The large mismatch in pCO₂ on the Siberian shelves does not show up in the CO₂ flux, as sea ice prevents CO₂ outgassing throughout most of the year.

Figure 18 Time series of simulated annual mean global ocean–atmosphere CO₂ flux (in PgC yr⁻¹) in the experiments conducted in this study. FESOM2.1–REcoM3 spin-up was conducted for 347 years (including 189 years of pre-spin-up; not shown in the plot) under repeat-year forcing taken from the year 1961 (RYF61). Here we show the spin-up since 1800 that has been continued as the control simulation B after 1958 for FESOM-1.4–REcoM2 (yellow) and FESOM2.1–REcoM3 (magenta), with a constant CO₂ concentration of 278 ppm (dashed lines) and the spin-up under increasing CO₂ that has been continued as simulation A after 1958 (solid lines). The control simulation B started in the year 1958 and was conducted for 64 years with RYF61 (dashed lines). Simulation A also started in 1958 and was forced with interannual varying forcing JRA55-do-1.5.0 (solid lines). Please note that the spin-up periods for FESOM1.4–REcoM2 and FESOM2.1–REcoM3 differ from each other, with the latter being longer than the former.

We continue our investigation with the analysis of the global ocean–atmosphere CO₂ flux time series (Fig. 18). In 1800, the first year of spin-up after the first 189 years of pre-spin-up of simulation B, the global ocean–atmosphere CO₂ flux was already relatively stable and converged towards a value close to zero. Under the assumption that the ocean and atmosphere were in equilibrium at constant preindustrial CO₂ values and without riverine carbon being transported into the ocean (Aumont et al., 2001; Resplandy et al., 2018; Regnier et al., 2022), an equilibrium flux of zero is expected for simulation B. Any deviation from this can be considered a bias (Hauck et al., 2020). The global bias of the annual air–sea CO₂ flux in the FESOM2.1–REcoM3 control simulation amounts to −0.12 PgC yr⁻¹. The control simulation conducted with the older model version FESOM1.4 had a larger bias, with a positive flux of around 0.4 PgC yr⁻¹ at the end of the simulation. In addition to the bias, the drift is reduced from 0.00264 PgC yr⁻² in FESOM1.4–REcoM2 to −0.00011 PgC yr⁻² in FESOM2.1–REcoM3 with a longer spin-up. Despite different spin-up procedures (FESOM1.4 has a shorter spin-up period), simulation A with both FESOM2.1 and FESOM1.4 reveals similar CO₂ fluxes under interannually varying forcing after 1980, which indicates a dominance of the forcing over the initial conditions. This also questions the common assumption that the same bias occurs in the control and historical simulations.

Figure 19 Globally integrated annual air–sea CO₂ flux from 10 global ocean biogeochemistry models (GOBMs) and seven pCO₂-based data products used in the Global Carbon Budget 2022 (Table 4 in Friedlingstein et al., 2022a), namely after applying bias correction to the models and river flux adjustment of 0.65 PgC yr⁻¹ (Regnier et al., 2022) to the pCO₂ products (Wright et al., 2021; Schwinger et al., 2016; Lacroix et al., 2021; Berthet et al., 2019; Hauck et al., 2020; Liao et al., 2020; Doney et al., 2009; Aumont et al., 2015; Nakano et al., 2011; Urakawa et al., 2020; Long et al., 2021a; Landschützer et al., 2016; Rödenbeck et al., 2022; Chau et al., 2022; Gloege et al., 2021; Zeng et al., 2014; Iida et al., 2015; Gregor and Gruber, 2021). The thick black line indicates the model ensemble mean, and the thick blue line shows the mean of the pCO₂ product ensemble. Thin dashed lines are from individual GOBMs and pCO₂ products. FESOM2.1–REcoM3 (magenta) shows the ocean carbon flux for the period of 1959–2021, whereas FESOM1.4–REcoM2 (yellow) covers the period from 1959–2019. Positive numbers indicate a flux into the ocean.

We next assess the model performance of the interannually varying simulation A by comparing it with the Global Carbon Budget's ensemble of pCO₂-based data products and other ocean biogeochemistry models (Fig. 19). Note that all model time series shown in Fig. 19 are referenced relative to their control simulations, with a constant atmospheric CO₂ concentration and without climate change forcing (simulation B). Although being consistent with the interannual variability, air–sea CO₂ fluxes of FESOM1.4 are at the lower end of the range compared to other global ocean biogeochemistry models and pCO₂-based estimates. In contrast, starting from the mid-1960s, FESOM2.1 shows a higher CO₂ flux in comparison to FESOM1.4. Considering the fact that both model versions do not differ much from each other in simulation A, the increase in the net CO₂ flux is mostly attributed to the level of CO₂ fluxes in their control simulations (Fig. 18).

After accounting for the bias in simulation B, the simulated ocean carbon sink (1990–1999) is 1.74 ± 0.11 and 2.17 ± 0.13 PgC yr⁻¹ for FESOM1.4–REcoM2 and FESOM2.1–REcoM3, respectively. Hence, FESOM2.1–REcoM3 is closer to the best estimate for the 1990s (2.2 ± 0.4 PgC yr⁻¹; IPCC; based on seven different methodologies; Denman et al., 2007; Ciais et al., 2014) than FESOM1.4–REcoM2. The cumulative uptake over the period of 1959–2019 amounts to 93.4 PgC (FESOM1.4–REcoM2) and 116.6 PgC (FESOM2.1–REcoM3), which is a 25 % increase in CO₂ flux. Yet, the FESOM2.1–REcoM3 CO₂ fluxes have been lower than the mean of the pCO₂-based data products since about 2008 and thus affirm the growing discrepancy between global ocean biogeochemistry models and pCO₂ products (Friedlingstein et al., 2022a). It is likely that the models underestimate the mean ocean carbon uptake (Friedlingstein et al., 2022a), which is linked to biases in ventilation (Goris et al., 2018; Terhaar et al., 2021, 2022; Bourgeois et al., 2022) and surface ocean buffer capacity (Vaittinada Ayar et al., 2022; Terhaar et al., 2022), and it is thus encouraging that FESOM2.1–REcoM3 has a comparatively high mean flux (Fig. 19). The pCO₂ products are statistical models that interpolate and extrapolate sparse pCO₂ observations and have substantial uncertainties themselves. In particular, they are sensitive to sparse and unevenly distributed observations (e.g. Gloege et al., 2021; Hauck et al., 2023), and it was shown that two of these methods overestimate the decadal CO₂ flux trend (2000–2018) by 20 %–35 %, based on the current pCO₂ observation distribution using a synthetic data set (Hauck et al., 2023).

3.3.3 DIC inventory changes

The interior ocean DIC inventory in FESOM2.1–REcoM3 amounts to about 38 200 PgC, which is in the reported range of 37 200 ± 200 to 39 000 PgC (Sundquist, 1985; Keppler et al., 2020). The DIC inventory is thought to change primarily in response to the rise in atmospheric CO₂, and the resulting DIC inventory change is often referred to as anthropogenic carbon. Effects of climate change (warming and circulation changes) are thought to be 1 order of magnitude smaller, have the opposite sign, and affect both the natural carbon cycle and the anthropogenic carbon uptake and inventory (see Hauck et al., 2020; Friedlingstein et al., 2022a; Crisp et al., 2022, for more details on the different simulations and carbon components). The observation-based estimates of the anthropogenic DIC inventory change use back-calculation techniques to separate anthropogenic carbon changes from the vast natural carbon reservoir. Thus, we here analyse the DIC inventory change in FESOM–REcoM for simulation A minus simulation B, which quantifies the total DIC inventory change while accounting for model drift. In addition, to derive the comparable DIC inventory change component as in observation-based studies, we make use of a third simulation (called simulation D), which is forced by interannually varying climate and preindustrial atmospheric CO₂. Quantifying the DIC inventory change over a specific period from simulation A minus D is then coherent with the anthropogenic carbon definition used in Gruber et al. (2019).

Sabine et al. (2004)Gruber et al. (2019)Gruber et al. (2019)

Table 4 FESOM2.1–REcoM3 DIC inventory for simulation A (in PgC) in 1994 and change in DIC inventory between 1800–1994 and 1994–2007 calculated from simulation A minus simulation B to account for model drift. Thus, the FESOM2.1–REcoM3 numbers encompass anthropogenic carbon cycle processes and the effect of climate change on the natural carbon cycle. Gruber et al. (2019) estimate the anthropogenic carbon inventory change, which is equivalent to simulation A minus simulation D (constant atmospheric CO₂ and variable climate). We have given the Gruber et al. (2019) anthropogenic plus the back-of-the-envelope natural carbon inventory changes in parenthesis, which is roughly comparable to simulation A minus simulation B (only available for the global ocean).

Figure 20 Maps of surface (0–10 m; a–c) and intermediate depth (300–500 m; d–f) concentration of simulated FESOM2.1–REcoM3 (simulation A) dissolved O₂ (mmol m⁻³; a, d), World Ocean Atlas 2018 climatology of dissolved O₂ (b, e; Garcia et al., 2019b) and corresponding differences (c, f) over the time period from 2012–2021.

The DIC inventory grew over time, in accordance with observation-based estimates (Table 4; Sabine et al., 2004; Gruber et al., 2019). The simulated anthropogenic DIC inventory change 1800–1994 (119 PgC) is in good agreement with the observation-based anthropogenic DIC inventory change (118 ± 19 PgC). For the total DIC inventory change, FESOM2.1–REcoM3 estimates a somewhat higher number (121 PgC; compared to 111 ± 21 PgC) but is well within the reported uncertainty in the observation-based estimate. For the period 1994–2007, the total DIC increase is 29.9 PgC (simulation A minus simulation B; i.e. drift corrected) and thus slightly lower than the estimate by Gruber et al. (2019). However, Gruber et al. (2019) only quantify the ocean anthropogenic DIC inventory increase (34 ± 4 PgC) and neglect the counter-effect of climate change. When considering the poorly constrained response of the natural carbon inventory to climate change, the model estimate falls within the uncertainty range of Gruber et al. (2019, 29 ± 5 PgC). Alternatively, estimating the anthropogenic DIC inventory change from simulations A minus D leads to 30.9 PgC, which is within the uncertainty range of Gruber et al. (2019)'s 34 ± 4 PgC. In GCB 2022, only four models simulated an anthropogenic DIC inventory change (1994–2007; simulation A minus D) ≥30 PgC (i.e. within the Gruber et al., 2019, uncertainty range). The other six models ranged between 25.5 and 28.3 PgC, and the model ensemble mean was 28.3 ± 2.6 PgC (Friedlingstein et al., 2022a). FESOM2.1–REcoM3 is thus one of the few ocean biogeochemistry models that falls within the range of interior ocean anthropogenic carbon accumulation that is also supported by ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratios (Tohjima et al., 2019) and atmospheric inversions (also see the discussion in Friedlingstein et al., 2022b). Notably, FESOM2.1–REcoM3 reproduces the latitudinal distribution of anthropogenic carbon accumulation in 1994–2007, with the maximum in the tropics (30^∘ S–30^∘ N), followed by the Southern Ocean (south of 30^∘ S), and the north (north of 30^∘ N). However, it also underestimates the accumulation in the tropics, as most other models do (Friedlingstein et al., 2022a). If the observation-based assessment of the DIC inventory changes in the north, tropics, and south is correct, then this may indicate transport of anthropogenic carbon from the Southern Ocean into the tropics that is too weak or an air–sea CO₂ flux in the tropics with too little ocean uptake (or too much release) of CO₂.

3.4 Oxygen

The simulated distributions of global O₂ concentration at the surface ocean and intermediate depths was consistent with the observed patterns in the World Ocean Atlas (WOA) 2018 (Fig. 20; with r=0.98/0.91 and RMSE = 19.6/38.4 mmol m³ for surface (0–10 m) and intermediate depths (300–500 m), respectively). The model successfully reproduced the typical spatial patterns (Schmidtko et al., 2017), including (1) oxygen minimum zones in the western boundary upwelling systems, where old deoxygenated waters are brought to the surface, (2) high concentrations in the high-latitude regions, where cold temperature increases oxygen solubility (Arctic and Southern oceans), and (3) moderate oxygen concentrations in the more stratified tropical gyres. Nevertheless, there were regional discrepancies. At the surface, the model slightly underestimated O₂ concentrations in the high-latitude surface ocean. At intermediate depth, the model generally overestimated the oxygen levels, especially in the Pacific Ocean and the subpolar Southern Ocean, with biases exceeding 100 mmol m⁻³. Within the 100–600 m layer, FESOM2.1–REcoM3 performed remarkably well, with simulated values of about 160 ± 105 mmol m⁻³, which is very close to the observations from WOA 2018 (158 ± 103 mmol m³). A previous model intercomparison study of oxygen concentrations within the 100–600 m layer (Cocco et al., 2013) showed that such performances are common, as all evaluated models fell within the error range of observations.

4 Conclusions and outlook

We have presented the new coupled ocean biogeochemistry model FESOM2.1–REcoM3. Building upon finite volumes for the ocean component improves the numerical efficiency and leads to higher numerical throughput of the coupled model (Danilov et al., 2017). Furthermore, the biogeochemistry component was extended to incorporate state-of-the-art carbonate chemistry routines, a second zooplankton and detritus group, and simulate the cycling of oxygen in the ocean. In its present configuration, the overall realism of FESOM2.1–REcoM3 in simulating the observed mean biogeochemical state is comparable to that of most GOBMs, while being among the more realistic models for estimating the global ocean anthropogenic carbon uptake. There are still a number of model shortcomings, such as a lower simulated NPP and regional misfit between the annual mean CO₂ flux of the model simulation and the pCO₂-based data product that will be addressed in the future.

This model set-up provides the basis for further model development, for example, the inclusion of coccolithophores as an additional phytoplankton functional type and the sensitivity of phytoplankton growth to rising CO₂ (Seifert et al., 2022) and the separation of the generic small zooplankton group into micro- and mesozooplankton that reduces model biases in nutrient fields, increases net primary production, and better captures the top-down control on phytoplankton bloom phenology (Karakuş et al., 2022). We further plan to incorporate more detailed iron biogeochemistry, as developed in REcoM coupled to MITgcm (e.g. Ye et al., 2020), and the explicit representation of the effects of viscosity and ballasting on the particle sinking speed, as well as oxygen-dependent remineralization, following Cram et al. (2018) to address knowledge gaps in carbon export and transfer to depth (Henson et al., 2022). Other on-going work addresses the role of rivers for carbon and nutrient transport into the ocean and the remineralization timescale of this river-derived organic material (Aumont et al., 2001; Lacroix et al., 2020; Regnier et al., 2022) and thus tackles a major uncertainty in the ocean carbon cycle and that complicates the comparison of ocean carbon sink estimates based on pCO₂ products and ocean biogeochemistry models (e.g. Hauck et al., 2020).

Appendix A: Equations

This appendix provides an overview of the underlying model equations and lists all biogeochemical variables of FESOM2.1–REcoM3. Changes in state variables in REcoM3 are controlled by biological and chemical processes, in addition to the changes induced by ocean circulation, mixing, diffusion, and advection computed by FESOM2.1. While some variables exchange across the ocean surface and/or the sea floor, others, like dead organic matter (detritus), sink through the water column. The concentration change for a state variable S is formulated as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{\partial S}{\partial t}}=-\mathit{U}\cdot \mathbf{\nabla }S+\mathbf{\nabla }\cdot \left(\mathit{\kappa }\cdot \mathbf{\nabla }S\right)+\mathrm{SMS}\left(S\right),\end{array}$$

where S is the volumetric concentration of a state variable, U is the 3-dimensional advection velocity, and κ is the diffusivity. The term SMS(S) represents the biogeochemical sources minus sinks. The slow-sinking detritus class is assumed to sink with a velocity, which increases linearly with depth as a first-order description of the shift to larger and faster-sinking particles with increasing depth (Kriest and Oschlies, 2008). A constant sinking rate is applied to the fast-sinking detritus class. REcoM3 has 28 oceanic and four explicit benthic state variables (Tables A1 and A2).

Figure A1 Maps of simulated FESOM2.1–REcoM3 (simulation A) vertically integrated net primary production (mgC m⁻² d⁻¹) of small phytoplankton (a), diatoms (b) and the sum of both phytoplankton groups (c). The satellite-based carbon-based productivity model (CbPM) is shown (d; Westberry et al., 2008), with corresponding differences between FESOM2.1–REcoM3 and VGPM (e). All fields are averaged over the time period from 2012 to 2021.

Figure A2 Conceptual diagram of the ocean biogeochemical model REcoM3. The 28 tracers can be grouped (indicated by boxes) into dissolved nutrients, carbonate system parameters and oxygen (upper left), phytoplankton functional types (centre), zooplankton functional types (upper right), two detritus classes (lower right), and dissolved organic material (lower left). Source and sink terms are depicted by arrows. For reasons of diagrammatic clarity, connections of dissolved oxygen (Oxy) to other state variables are omitted here. Similarly, the release of alkalinity, dissolved inorganic nutrients, and organic matter from the sediment are not shown.

A1 Sources minus sinks

A1.1 Nutrients

Dissolved inorganic nitrate (DIN)

The simulated DIN conceptually represents the concentrations of nitrate, nitrite, and ammonia, while in practice only nitrate is considered. The concentration of DIN in the water column rises when dissolved organic nitrogen (DON) is remineralized and diminishes as a consequence of assimilation by small phytoplankton and diatoms as follows:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{DIN}\right)& =\underset{\text{DON remineralization}}{\underbrace{{\mathit{\rho }}_{\mathrm{DON}}\cdot f_{\mathrm{T}}\cdot \mathrm{DON}}}\\ & -\underset{\text{N assimilation, small phytoplankton}}{\underbrace{V_{\mathrm{small}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}-\underset{\text{N assimilation, diatoms}}{\underbrace{V_{\mathrm{dia}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}.\end{array}\end{array}$$

The state variables DON, PhyC_small, and PhyC_dia are listed in Table A1. The value of the remineralization rate constant (ρ_DON) is given in Table A8. The temperature dependency of remineralization (f_T) is calculated in Eq. (A33). See Sect. A3.4 for details on the carbon-specific nitrogen assimilation rates $V_{\mathrm{small}}^{\mathrm{N}}$ and $V_{\mathrm{dia}}^{\mathrm{N}}$ (Table A5).

Dissolved silicic acid (DSi)

Silicon assimilation (Si assimilation) increases when biogenic silica from one of the two detritus classes dissolves.

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{DSi}\right)& =\underset{\text{Remineralization, slow-sinking detritus}}{\underbrace{{\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}\cdot \mathrm{DetSi}}}\\ & +\underset{\text{Remineralization, fast-sinking detritus}}{\underbrace{{\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Si}}}-\underset{\text{Si assimilation, diatoms}}{\underbrace{V^{\mathrm{Si}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\end{array}\end{array}$$

The state variables PhyC_dia, DetSi, and DetZ2Si are listed in Table A1. The temperature-dependent remineralization rate of silicon (${\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}$) and the carbon-specific Si assimilation rate (V^Si) are calculated in Eqs. (A35) and (A41), respectively (Table A5).

Dissolved iron (DFe)

Excretion of phyto- and zooplankton and remineralization of detritus release iron with a fixed iron:nitrate ratio (q^Fe:N). Unlike for nitrogen, which is released as dissolved organic nitrogen and needs to be remineralized further to become available as nutrient again, the released iron is put directly into the dissolved pool iron, basically assuming that all dissolved iron is ultimately bio-available. Iron assimilation (again assumed to be proportional to nitrogen assimilation; hereafter N assimilation) by both phytoplankton classes lower the level of dissolved iron. In addition, free inorganic iron Fe^′ is scavenged onto sinking particles, with a rate that is proportional to particle concentration. We take detrital carbon as a proxy for the mass of sinking particles.

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{DFe}\right)& =q^{\mathrm{Fe}:\mathrm{N}}\cdot (\underset{\text{Excretion, small phytoplankton}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyN}}_{\mathrm{small}}}}+\underset{\text{Excretion, diatoms}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyN}}_{\mathrm{dia}}}}\\ & +\underset{\text{Remineralization, slow-sinking detritus}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetN}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetN}}}+\underset{\text{Remineralization, fast-sinking detritus}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetN}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{N}}}\\ & +\underset{\text{Excretion, small zooplankton}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{N}}\cdot \mathrm{ZooN}}}+\underset{\text{Excretion, macrozooplankton}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{N}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{N}}}\\ & -\underset{\text{N assimilation, small phytoplankton}}{\underbrace{V_{\mathrm{small}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}-\underset{\text{N assimilation, diatom}}{\underbrace{V_{\mathrm{dia}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}})\\ & -\underset{\text{Scavenging, slow-sinking detritus}}{\underbrace{{\mathit{\kappa }}_{\mathrm{Fe}}\cdot \mathrm{DetC}\cdot {\mathrm{Fe}}^{\prime }}}-\underset{\text{Scavenging, fast-sinking detritus}}{\underbrace{{\mathit{\kappa }}_{\mathrm{Fe}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{C}\cdot {\mathrm{Fe}}^{\prime }}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, PhyN_small, PhyN_dia, DetC, DetN, DetZ2C, DetZ2N, ZooN, and Zoo2N are listed in Table A1. The intracellular Fe:N ratio (q^Fe:N) and scavenging rate of iron (κ_Fe) are given in Table A4. Excretion rates (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}$, ${\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{N}}$, and ${\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{N}}$) and the degradation rate for detritus N (ρ_DetN) are listed in Table A8. The temperature dependency (f_T) is calculated in Eq. (A33). The limitation of intracellular nitrogen ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$; $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$) is described in Eq. (A45). Scavenging is calculated following Parekh et al. (2004). The total concentration of dissolved iron (Fe_T) is separated into free iron (Fe^′) and iron complexed with organic ligands (Fe_L), which is not scavenged. Complexation reactions are fast (Tagliabue and Völker, 2011), so we assume instantaneous equilibrium between free iron and free ligand (L^′), which is computed using a constant $K_{{\mathrm{Fe}}_{\mathrm{L}}}=\frac{\left[{\mathrm{Fe}}^{\prime }\right]\cdot \left[{\mathrm{L}}^{\prime }\right]}{\left[{\mathrm{Fe}}_{\mathrm{L}}\right]}$, by solving the following:

$$\begin{array}{}\text{(A5)}& {\mathrm{Fe}}_{\mathrm{T}}={\mathrm{Fe}}^{\prime }+{\mathrm{Fe}}_{\mathrm{L}}{\mathrm{L}}_{\mathrm{T}}={\mathrm{Fe}}_{\mathrm{L}}+{\mathrm{L}}^{\prime }.\end{array}$$

For simplicity, we assume here a constant total ligand concentration L_T, unlike in Völker and Tagliabue (2015). Variable ligand concentration, like in Misumi et al. (2011) or Völker and Tagliabue (2015), or variable ligand binding strength, like in Ye et al. (2020), will be explored in the future. The values for $K_{{\mathrm{Fe}}_{\mathrm{L}}}$ and L_T are listed in Table A4.

A1.2 Carbon cycle

Dissolved inorganic carbon (DIC)

DIC concentration increases with respiration of phyto- and zooplankton, remineralization of semi-labile dissolved organic carbon, dissolution of calcitic detritus, and dissolution of CaCO₃ in zooplankton guts. Loss terms are carbon fixation by primary producers and the formation of calcium carbonate. In addition, the sea–air flux of CO₂ leads to an exchange of carbon with the atmosphere, depending on the partial pressure difference in the CO₂ between ocean and atmosphere. This exchange is treated separately as a boundary condition. The partial pressure of surface ocean CO₂ is computed using the mocsy 2.0 routines (Orr and Epitalon, 2015).

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{DIC}\right)& =\underset{\text{Net respiration, small phytoplankton}}{\underbrace{\left(r_{\mathrm{small}}-P_{\mathrm{small}}\right)\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}+\underset{\text{Net respiration, diatom}}{\underbrace{\left(r_{\mathrm{dia}}-P_{\mathrm{dia}}\right)\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & +\underset{\text{Remineralization of DOC}}{\underbrace{{\mathit{\rho }}_{\mathrm{DOC}}\cdot f_{\mathrm{T}}\cdot \mathrm{DOC}}}+\underset{\text{Respiration, small zoo.}}{\underbrace{r_{\mathrm{zoo}}\cdot \mathrm{ZooC}}}\\ & +\underset{\text{Respiration, macrozoo.}}{\underbrace{r_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{C}}}+\underset{\text{Calcite dissolution, slow-sinking detritus}}{\underbrace{{\mathrm{Diss}}_{\mathrm{calc}}\cdot \mathrm{DetCalc}}}\\ & +\underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution in guts, small zoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}\cdot {\mathrm{Diss}}_{\text{calc\_guts}}}}-\underset{\mathrm{Calcification}}{\underbrace{\mathit{\psi }\cdot P_{\mathrm{small}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & +\underset{\text{Calcite dissolution, fast-sinking detritus}}{\underbrace{{\mathrm{Diss}}_{\mathrm{calc}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Calc}}}\\ & +\underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution in guts, macrozoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}\cdot {\mathrm{Diss}}_{\text{calc\_guts}}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, DOC, ZooC, Zoo2C, DetCalc, and DetZ2Calc are listed in Table A1. Respiration rate constants of small phytoplankton (r_small), diatoms (r_dia), and zooplankton groups (r_zoo and r_zoo2) are computed in Sects. A3.2 and A4.1, respectively. Photosynthesis terms (P_small and P_dia) are calculated in Eq. (A36). The remineralization rate constant (ρ_DOC) is listed in Table A8, and the temperature dependency (f_T) is given in Eq. (A33). Calcite dissolution by detritus (Diss_calc; Diss_calc2) is calculated in Eq. (A28). The constant for the dissolution of calcium carbonate in zooplankton guts (Diss_{calc_guts}) is listed in Table A5. $G_{\mathrm{small}}^{\mathrm{zoo}}$ and $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$ are grazing terms and explained in Sect. A4.2. The value of the calcite production ratio (ψ) is given in Table A3.

Total alkalinity (Alk)

The balance of alkalinity is affected by primary production, remineralization of dissolved organic matter, dissolution of calcitic detritus, and dissolution of CaCO₃ in zooplankton guts. Alkalinity increases when nitrogen is assimilated and when CaCO₃ is dissolved (Wolf-Gladrow et al., 2007). Simultaneously, it is reduced by the calcification and remineralization of dissolved organic nitrogen. The effect of phosphate assimilation and remineralization onto alkalinity is taken into account, assuming a constant N:P Redfield ratio (16:1).

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{Alk}\right)& =(\mathrm{1}+\mathrm{1}/\mathrm{16})\cdot \underset{\text{N assimilation, small phytoplankton}}{\underbrace{V_{\mathrm{small}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & +(\mathrm{1}+\mathrm{1}/\mathrm{16})\cdot \underset{\text{N assimilation, diatom}}{\underbrace{V_{\mathrm{dia}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & -(\mathrm{1}+\mathrm{1}/\mathrm{16})\cdot \underset{\text{Remineralization of DON}}{\underbrace{{\mathit{\rho }}_{\mathrm{DON}}\cdot f_{\mathrm{T}}\cdot \mathrm{DON}}}\\ & -\mathrm{2}\cdot \underset{\mathrm{Calcification}}{\underbrace{\mathit{\psi }\cdot P_{\mathrm{small}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & +\mathrm{2}\cdot \underset{\text{Calcite dissolution, slow-sinking detritus}}{\underbrace{{\mathrm{Diss}}_{\mathrm{calc}}\cdot \mathrm{DetCalc}}}\\ & +\mathrm{2}\cdot \underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution in guts, small zoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}\cdot {\mathrm{Diss}}_{\text{calc\_guts}}}}\\ & +\mathrm{2}\cdot \underset{\text{Calcite dissolution, fast-sinking detritus}}{\underbrace{{\mathrm{Diss}}_{\mathrm{calc}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Calc}}}\\ & +\mathrm{2}\cdot \underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution in guts, macrozoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}\cdot {\mathrm{Diss}}_{\text{calc\_guts}}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, DON, DetCalc, and DetZ2Calc are listed in Table A1. The N assimilation ($V_{\mathrm{small}}^{\mathrm{N}}$ and $V_{\mathrm{dia}}^{\mathrm{N}}$) is calculated in Sect. A3.4. The remineralization rate constant (ρ_DON) is given in Table A8. The temperature dependency (f_T) is calculated in Eq. (A33). The value of the calcite production ratio (ψ) is given in Table A3. The photosynthesis term (P_small) is calculated in Eq. (A36). The calcite dissolution by detritus (Diss_calc, Diss_calc2) is calculated in Eq. (A28). Dissolution of calcium carbonate in guts (Diss_{calc_guts}) is listed in Table A5. $G_{\mathrm{small}}^{\mathrm{zoo}}$ and $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$ are grazing terms and explained in Sect. A4.2.

A1.3 Phytoplankton

Nitrogen

The phytoplankton nitrogen pools increase through N assimilation. The assimilation process is assumed to be proportional to carbon biomass, with a carbon-specific uptake rate that depends on the C:N ratio of phytoplankton and the external DIN concentration (Geider et al., 1998). Excretion of biogenic nitrogen to semi-labile DON drains the pool. At a high intracellular C:N ratio, the excretion is downregulated. Aggregation and grazing by the two zooplankton groups transfer nitrogen to the zooplankton and detritus pools.

$$\begin{array}{}\text{(A8)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left({\mathrm{PhyN}}_{\mathrm{small}}\right)& =\underset{\text{N assimilation}}{\underbrace{V_{\mathrm{small}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & -\underset{\text{DON excretion}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyN}}_{\mathrm{small}}}}-\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhyN}}_{\mathrm{small}}}}\\ & -\underset{\text{Grazing loss by small zoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}}}-\underset{\text{Grazing loss by macrozoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}}}\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left({\mathrm{PhyN}}_{\mathrm{dia}}\right)& =\underset{\text{N assimilation}}{\underbrace{V_{\mathrm{dia}}^{\mathrm{N}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & -\underset{\text{DON excretion}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyN}}_{\mathrm{dia}}}}-\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhyN}}_{\mathrm{dia}}}}\\ & -\underset{\text{Grazing loss by small zoo.}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}}}}-\underset{\text{Grazing loss by macrozoo.}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyN_small, PhyC_dia, and PhyN_dia are listed in Table A1. The N assimilation ($V_{\mathrm{small}}^{\mathrm{N}}$ and $V_{\mathrm{dia}}^{\mathrm{N}}$) is explained in Sect. A3.4. The constant excretion rate constant (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}$) is given in Table A8. When the C:N ratio of the cells becomes too high, excretion of DON is downregulated by the limiter function ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$; $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$) that is described in Eq. (A45). Phytoplankton aggregation (Agg) defines the transfer of nitrogen into the detritus pools, which depends quadratically on detritus and phytoplankton concentrations (Eq. A42). Grazing loss terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$, $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}}$, and $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$) are explained in Sect. A4.2.

Carbon

The carbon biomass of small phytoplankton and diatoms increases as a result of carbon assimilation during photosynthesis. Loss terms include excretion of DOC, which is limited by the availability of proteins as in the nitrogen pool, respiration, aggregation, and grazing.

$$\begin{array}{}\text{(A10)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left({\mathrm{PhyC}}_{\mathrm{small}}\right)& =\underset{\text{Net photosynthesis}}{\underbrace{(P_{\mathrm{small}}-r_{\mathrm{small}})\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & -\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}-\underset{\text{Excretion of DOC}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & -\underset{\text{Grazing loss by small zoo.}}{\underbrace{q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}\cdot G_{\mathrm{small}}^{\mathrm{zoo}}}}-\underset{\text{Grazing loss by macrozoo.}}{\underbrace{q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}\cdot G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}}}\end{array}\text{(A11)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left({\mathrm{PhyC}}_{\mathrm{dia}}\right)& =\underset{\text{Net photosynthesis}}{\underbrace{(P_{\mathrm{dia}}-r_{\mathrm{dia}})\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & -\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}-\underset{\text{Excretion of DOC}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & -\underset{\text{Grazing loss by small zoo.}}{\underbrace{q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}}\cdot G_{\mathrm{dia}}^{\mathrm{zoo}}}}-\underset{\text{Grazing loss by macrozoo.}}{\underbrace{q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}}\cdot G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}}}\end{array}\end{array}$$

The state variables PhyC_small and PhyC_dia are listed in Table A1. The photosynthesis terms (P_small and P_dia) are calculated in Eq. (A36). The rates of respiration by small phytoplankton (r_small) and diatoms (r_dia) are explained in Sect. A3.2. The constant for DOC excretion rate of phytoplankton (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}$; Table A8) is downregulated by the limiter factor ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$; $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$) when the N:C ratio becomes too high (Eq. A45). Phytoplankton aggregation (Agg) is calculated in Eq. (A42). Grazing terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$, $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}}$, and $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$) are explained in Sect. A4.2. $q^{\mathrm{C}:\mathrm{N}}=\mathrm{PhyC}/\mathrm{PhyN}$, is used to convert the grazing units from millimoles of N to millimoles of C.

CaCO₃

The formation of biogenic calcium carbonate in our model is limited to coccolithophores only, which are assumed to form a constant fraction of the non-diatom phytoplankton. Formation of CaCO₃ by heterotrophs, such as foraminifera or pteropods is neglected. Biogenic CaCO₃ produced by coccolithophores is transformed into detritus CaCO₃ with all forms of organic carbon loss, including organic matter excretion, respiration, aggregation, and grazing. Calcifiers are assumed to comprise a certain fraction of the total small phytoplankton concentration, which is specified by the parameter ψ (Table A3), thus tying the calcite production of calcifiers to the growth of small phytoplankton.

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{PhyCalc}\right)& =\underset{\mathrm{Calcification}}{\underbrace{\mathit{\psi }\cdot P_{\mathrm{small}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}-\underset{\mathrm{Respiration}}{\underbrace{r_{\mathrm{small}}\cdot \mathrm{PhyCalc}}}\\ & -\underset{\text{Grazing loss, small zoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}}}-\underset{\text{Grazing loss, macrozoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}}}\\ & -\underset{\text{Excretion loss}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot \mathrm{PhyCalc}}}-\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot \mathrm{PhyCalc}}}\end{array}\end{array}$$

The state variables PhyC_small and PhyCalc are listed in Table A1. The value of the calcite production ratio (ψ) is given in Table A3. The constant excretion rate (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}$; Table A8) is downregulated by the limiter factor $f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$ (Eq. A45) when the N:C ratio becomes too high. Photosynthesis (P_small), respiration (r_small), and the aggregation of phytoplankton (Agg) rates are calculated in Eqs. (A36), (A38), and (A42), respectively. Grazing terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$ and $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$) are explained in Sect. A4.2. $q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}=\mathrm{PhyCalc}/{\mathrm{PhyN}}_{\mathrm{small}}$ is used to convert the grazing units from millimoles of N to millimoles of CaCO₃.

Diatom silicon

The silica frustule of diatoms is built through Si assimilation, which we assume to be carbon specific and regulated by cellular quotas (see below). Any decrease in N biomass through excretion, grazing, or aggregation leads to a corresponding transfer of silica to the detritus silica pool.

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{PhySi}\right)& =\underset{\text{Diatom Si assimilation}}{\underbrace{V^{\mathrm{Si}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & -\underset{\text{Excretion to detritus}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhySi}}_{\mathrm{dia}}}}-\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhySi}}_{\mathrm{dia}}}}\\ & -\underset{\text{Grazing loss, small zoo.}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}}\cdot q^{\mathrm{Si}:\mathrm{N}}}}-\underset{\text{Grazing loss, macrozoo.}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}\cdot q^{\mathrm{Si}:\mathrm{N}}}}\end{array}\end{array}$$

The state variables PhyC_dia and PhySi_dia are described in Table A1. Si assimilation (V^Si) and aggregation rates (Agg) are calculated in Eqs. (A41) and (A42), respectively. The constant excretion rate (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}$; Table A8) is downregulated by the limiter factor $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$ (Eq. A45) when the N:C ratio becomes too high. Grazing terms ($G_{\mathrm{dia}}^{\mathrm{zoo}}$ and $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$) are explained in Sect. A4.2. The intracellular ratio between diatom silicon and nitrate is defined as $q^{\mathrm{Si}:\mathrm{N}}={\mathrm{PhySi}}_{\mathrm{dia}}/{\mathrm{PhyN}}_{\mathrm{dia}}$.

Chlorophyll a

Chlorophyll a synthesis is structured as a function of irradiance and of N assimilation, following Geider et al. (1998). Chlorophyll a is degraded at a light-dependent rate (see Álvarez et al., 2018) and lost via aggregation and grazing. The grazing losses in terms of nitrogen biomass are converted to chlorophyll loss using the intracellular Chl:N ratio.

$$\begin{array}{}\text{(A14)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left({\mathrm{PhyChl}}_{\mathrm{small}}\right)& =\underset{\mathrm{Chlorophyll}a\mathrm{synthesis}}{\underbrace{S_{\mathrm{small}}^{\mathrm{chl}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & -\underset{\text{Grazing loss, small zoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{\mathrm{Chl}:\mathrm{N}}}}-\underset{\text{Grazing loss, macrozoo.}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{\mathrm{Chl}:\mathrm{N}}}}\\ & -\underset{\text{Degradation loss}}{\underbrace{{\mathrm{deg}}_{\mathrm{small}}^{\mathrm{chl}}\cdot {\mathrm{PhyChl}}_{\mathrm{small}}}}-\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhyChl}}_{\mathrm{small}}}}\end{array}\text{(A15)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left({\mathrm{PhyChl}}_{\mathrm{dia}}\right)& =\underset{\mathrm{Chlorophyll}a\mathrm{synthesis}}{\underbrace{S_{\mathrm{dia}}^{\mathrm{chl}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & -\underset{\text{Grazing loss, small zoo.}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}}\cdot q_{\mathrm{dia}}^{\mathrm{Chl}:\mathrm{N}}}}-\underset{\text{Grazing loss, macrozoo.}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{dia}}^{\mathrm{Chl}:\mathrm{N}}}}\\ & -\underset{\text{Degradation loss}}{\underbrace{{\mathrm{deg}}_{\mathrm{dia}}^{\mathrm{chl}}\cdot {\mathrm{PhyChl}}_{\mathrm{dia}}}}-\underset{\text{Aggregation loss}}{\underbrace{\mathrm{Agg}\cdot {\mathrm{PhyChl}}_{\mathrm{dia}}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, PhyChl_small, and PhyChl_dia are listed in Table A1. The chlorophyll a synthesis ($S_{\mathrm{small}}^{\mathrm{chl}}$; $S_{\mathrm{dia}}^{\mathrm{chl}}$) and the aggregation (Agg) terms are calculated in Eqs. (A39) and (A42), respectively. The degradation parameters (${\mathrm{deg}}_{\mathrm{small}}^{\mathrm{chl}}$; ${\mathrm{deg}}_{\mathrm{dia}}^{\mathrm{chl}}$) are given in Table A8. Grazing terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$, $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}}$, and $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$) are explained in Sect. A4.2. The conversion factor from millimoles of N to milligrams of Chl a is defined as $q^{\mathrm{Chl}:\mathrm{N}}=\mathrm{PhyChl}/\mathrm{PhyN}$.

A1.4 Zooplankton

Nitrogen

Both zooplankton classes increase their nitrogen biomass via grazing on phytoplankton and detritus, while mortality and excretion of DON reduce it. Macrozooplankton further feeds on small zooplankton and releases nitrogen via fecal pellet production.

$$\begin{array}{}\text{(A16)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{ZooN}\right)& =\underset{\mathrm{Grazing}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}}\cdot G_{\mathrm{tot}}^{\mathrm{zoo}}}}-\underset{\text{Grazing loss, macrozoo.}}{\underbrace{G_{\mathrm{zoo}}}}\\ & -\underset{\mathrm{Mortality}}{\underbrace{m_{\mathrm{zoo}}\cdot {\mathrm{ZooN}}^{\mathrm{2}}}}-\underset{\text{Excretion of DON}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{N}}\cdot \mathrm{ZooN}}}\end{array}\text{(A17)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{Zoo}\mathrm{2}\mathrm{N}\right)& =\underset{\mathrm{Grazing}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}}\cdot G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}}}-\underset{\mathrm{Mortality}}{\underbrace{m_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}{\mathrm{N}}^{\mathrm{2}}}}\\ & -\underset{\text{Excretion of DON}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{N}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{N}}}-\underset{\text{Fecal pellet}}{\underbrace{f_{\mathrm{n}}\cdot G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}}}\end{array}\end{array}$$

The state variables ZooN and Zoo2N are listed in Table A1. Only a fraction of the grazed phytoplankton (γ_zoo and γ_zoo2; Table A3) enters the zooplankton biomass. The rest is transferred to detritus due to sloppy feeding. The grazing terms ($G_{\mathrm{tot}}^{\mathrm{zoo}}$ and $G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}$) are calculated in Sect. A4.2. The mortality parameter (m_zoo and m_zoo2) and fecal pellet production rate constant (f_n) are listed in Table A3. The DON excretion terms (${\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{N}}$ and ${\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{N}}$) are given in Table A8.

Carbon

The zooplankton carbon biomass increases with carbon uptake via grazing and decreases through carbon losses through mortality, respiration, and carbon excretion to the semi-labile DOC pool. Macrozooplankton further gains carbon by grazing on small zooplankton and loses it via fecal pellet production.

$$\begin{array}{}\text{(A18)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{ZooC}\right)& =\underset{\text{Grazing on phytoplankton}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}}\cdot (G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{dia}}^{\mathrm{zoo}}\cdot q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}})}}\\ & +\underset{\text{Grazing on detritus}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}}\cdot (G_{\det }^{\mathrm{zoo}}\cdot q_{\det }^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}}\cdot q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}})}}\\ & -\underset{\text{Grazing loss by macrozoo.}}{\underbrace{G_{\mathrm{zoo}}\cdot q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}}}\\ & -\underset{\text{Zooplankton mortality}}{\underbrace{m_{\mathrm{zoo}}\cdot {\mathrm{ZooN}}^{\mathrm{2}}\cdot q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}}}-\underset{\text{Respiration loss}}{\underbrace{r_{\mathrm{zoo}}\cdot \mathrm{ZooC}}}\\ & -\underset{\text{Excretion of DOC}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{C}}\cdot \mathrm{ZooC}}}\end{array}\text{(A19)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{Zoo}\mathrm{2}\mathrm{C}\right)& =\underset{\text{Grazing on phytoplankton}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}}\cdot (G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}})}}\\ & +\underset{\text{Grazing on detritus}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}}\cdot (G_{\det }^{\mathrm{zoo}\mathrm{2}}\cdot q_{\det }^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}})}}\\ & +\underset{\text{Grazing on small zoo.}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}}\cdot (G_{\mathrm{zoo}}\cdot q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}})}}\\ & -\underset{\text{Zooplankton mortality}}{\underbrace{m_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}{\mathrm{N}}^{\mathrm{2}}\cdot q_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}}}-\underset{\text{Respiration loss}}{\underbrace{r_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{C}}}\\ & -\underset{\text{Excretion of DOC}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{C}}}-\underset{\text{Fecal pellet}}{\underbrace{f_{\mathrm{c}}\cdot G_{\mathrm{cflux}}}}\end{array}\end{array}$$

The state variables ZooN, ZooC, Zoo2N, and Zoo2C are listed in Table A1. A fraction of the grazed phytoplankton (γ_zoo and γ_zoo2; Table A3) is kept in the zooplankton biomass, while the remainder is returned back to detritus pool as a consequence of sloppy feeding. Grazing terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}}$, $G_{\det }^{\mathrm{zoo}}$, $G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}}$, $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\det }^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}$, and G_zoo) are calculated in Sect. A4.2. The respiration terms of zooplankton (r_zoo and r_zoo2) are calculated in Eqs. (A50) and (A51). Mortality parameters (m_zoo and m_zoo2) are listed in Table A3. The DOC excretion terms (${\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{C}},{\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}}$) are in Table A8. The grazing flux in terms of nitrogen biomass is converted to carbon biomass using the respective intracellular C:N ratios ($q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}$, $q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}}$, $q_{\det }^{\mathrm{C}:\mathrm{N}}$, $q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}$, $q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}$, and $q_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}$), where $q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}={\mathrm{PhyC}}_{\mathrm{small}}/{\mathrm{PhyN}}_{\mathrm{small}}$, $q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}}={\mathrm{PhyC}}_{\mathrm{dia}}/{\mathrm{PhyN}}_{\mathrm{dia}}$, $q_{\det }^{\mathrm{C}:\mathrm{N}}=\mathrm{DetC}/\mathrm{DetN}$, $q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}=\mathrm{DetZ}\mathrm{2}\mathrm{C}/\mathrm{DetZ}\mathrm{2}\mathrm{N}$, $q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}=\mathrm{ZooC}/\mathrm{ZooN}$, and $q_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}=\mathrm{Zoo}\mathrm{2}\mathrm{C}/\mathrm{Zoo}\mathrm{2}\mathrm{N}$. Total grazed carbon biomass (G_cflux) and the fecal pellet production rate constant (f_c; Table A3) together determine the fraction of carbon being lost to the large detritus carbon pool via fecal pellets.

A1.5 Detritus

Nitrogen

Detrital nitrogen pool increases as a result of sloppy feeding and mortality. Sloppy feeding is outlined as a function of grazing fluxes and grazing efficiency of macrozooplankton. In other words, the grazed phytoplankton partly goes to the macrozooplankton biomass, depending on the grazing efficiency. The phytoplankton aggregation contributes only to slow-sinking detritus. Fecal pellet production is defined only for macrozooplankton group. Detritus is degraded to DON, based on temperature and a remineralization rate.

$$\begin{array}{}\text{(A20)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetN}\right)& =\underset{\text{Sloppy feeding}}{\underbrace{(G_{\mathrm{small}}^{\mathrm{zoo}}+G_{\mathrm{dia}}^{\mathrm{zoo}})\cdot (\mathrm{1}-{\mathit{\gamma }}_{\mathrm{zoo}})}}+\underset{\text{Zooplankton mortality}}{\underbrace{m_{\mathrm{zoo}}\cdot {\mathrm{ZooN}}^{\mathrm{2}}}}\\ & -\underset{\text{Grazing loss, small zoo.}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}}\cdot (G_{\det }^{\mathrm{zoo}}+G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}})}}\\ & +\underset{\text{Phytoplankton aggregation}}{\underbrace{\mathrm{Agg}\cdot ({\mathrm{PhyN}}_{\mathrm{small}}+{\mathrm{PhyN}}_{\mathrm{dia}})}}-\underset{\text{Degradation to DON}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetN}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetN}}}\end{array}\text{(A21)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetZ}\mathrm{2}\mathrm{N}\right)& =\underset{\text{Sloppy feeding}}{\underbrace{(G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}+G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}+G_{\mathrm{zoo}})\cdot (\mathrm{1}-{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}})}}\\ & -\underset{\text{Grazing loss, macrozoo.}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}}\cdot (G_{\det }^{\mathrm{zoo}\mathrm{2}}+G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}\cdot )}}\\ & +\underset{\mathrm{Mortality}}{\underbrace{m_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}{\mathrm{N}}^{\mathrm{2}}}}+\underset{\text{Fecal pellet}}{\underbrace{f_{\mathrm{n}}\cdot G_{\mathrm{tot}}}}\\ & -\underset{\text{Degradation to DON}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetN}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{N}}}\end{array}\end{array}$$

The state variables PhyN_small, PhyN_dia, ZooN, DetN, Zoo2N, and DetZ2N are listed in Table A1. The grazing efficiency (γ_zoo and γ_zoo2), mortality (m_zoo and m_zoo2), and fecal pellet production rate constant (f_n) are listed in Table A3. Grazing terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}}$, $G_{\det }^{\mathrm{zoo}}$, $G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}}$, $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\det }^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}$, and G_zoo) are calculated in Sect. A4.2. The remineralization rate constant of DON (ρ_DetN) is listed in Table A8. The temperature dependency f_T is calculated in Eq. (A33). The aggregation (Agg) term is calculated in Eq. (A42).

Carbon

Detrital carbon sources are associated with sloppy feeding, aggregation of phytoplankton, mortality of small zooplankton, and fecal pellet production by macrozooplankton. Degradation of DetC and DetZ2C to DOC is the only loss term.

$$\begin{array}{}\text{(A22)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetC}\right)& =\underset{\text{Sloppy feeding}}{\underbrace{(G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{dia}}^{\mathrm{zoo}}\cdot q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}})\cdot (\mathrm{1}-{\mathit{\gamma }}_{\mathrm{zoo}})}}\\ & +\underset{\text{small zoo. mortality}}{\underbrace{m_{\mathrm{zoo}}\cdot {\mathrm{ZooN}}^{\mathrm{2}}\cdot q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}}}\\ & -\underset{\text{Grazing loss by macrozoo.}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}}\cdot (G_{\det }^{\mathrm{zoo}}\cdot q_{\det }^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}}\cdot q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}})}}\\ & +\underset{\text{Phytoplankton aggregation}}{\underbrace{\mathrm{Agg}\cdot ({\mathrm{PhyC}}_{\mathrm{small}}+{\mathrm{PhyC}}_{\mathrm{dia}})}}-\underset{\text{Degradation to DOC}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetC}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetC}}}\end{array}\text{(A23)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetZ}\mathrm{2}\mathrm{C}\right)& =\underset{\text{Sloppy feeding}}{\underbrace{\begin{array}{c}(G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{zoo}}\cdot q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}})\\ \cdot (\mathrm{1}-{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}})\end{array}}}\\ & -\underset{\text{Grazing loss by macrozoo.}}{\underbrace{{\mathit{\gamma }}_{\mathrm{zoo}\mathrm{2}}\cdot (G_{\det }^{\mathrm{zoo}\mathrm{2}}\cdot q_{\det }^{\mathrm{C}:\mathrm{N}}+G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}})}}\\ & +\underset{\mathrm{Mortality}}{\underbrace{m_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}{\mathrm{N}}^{\mathrm{2}}\cdot q_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}}}+\underset{\text{Fecal pellet}}{\underbrace{f_{\mathrm{c}}\cdot G_{\mathrm{cflux}}}}\\ & -\underset{\text{Degradation to DOC}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetC}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{C}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, ZooN, DetC, Zoo2N, and DetZ2C are listed in Table A1. The grazing efficiency (γ_zoo and γ_zoo2) and mortality (m_zoo and m_zoo2) parameters are listed in Table A3. Grazing terms ($G_{\mathrm{small}}^{\mathrm{zoo}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}}$, $G_{\det }^{\mathrm{zoo}}$, $G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}}$, $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$, $G_{\det }^{\mathrm{zoo}\mathrm{2}}$, $G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}$, and G_zoo) are calculated in Sect. A4.2. The remineralization rate of DOC (ρ_DetC) is listed in Table A8. Temperature dependency f_T is calculated in Eq. (A33). The aggregation (Agg) term is calculated in Eq. (A42). The total grazed carbon biomass (G_cflux) and the fecal pellet production rate constant (f_c; Table A3) together determine the fraction of carbon being lost to the large detritus carbon pool via fecal pellets. The quotas $q_{\mathrm{small}}^{\mathrm{C}:\mathrm{N}}={\mathrm{PhyC}}_{\mathrm{small}}/{\mathrm{PhyN}}_{\mathrm{small}}$, $q_{\mathrm{dia}}^{\mathrm{C}:\mathrm{N}}={\mathrm{PhyC}}_{\mathrm{dia}}/{\mathrm{PhyN}}_{\mathrm{dia}}$, $q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}=\mathrm{ZooC}/\mathrm{ZooN}$, $q_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}=\mathrm{Zoo}\mathrm{2}\mathrm{C}/\mathrm{Zoo}\mathrm{2}\mathrm{N}$, $q_{\det }^{\mathrm{C}:\mathrm{N}}=\mathrm{DetC}/\mathrm{DetN}$, and $q_{\mathrm{detZ}\mathrm{2}}^{\mathrm{C}:\mathrm{N}}=\mathrm{DetZ}\mathrm{2}\mathrm{C}/\mathrm{DetZ}\mathrm{2}\mathrm{N}$ are used to convert the units from mmol N to mmol C.

Silica

Biogenic detrital silica increases with excretion fluxes from diatoms to detritus, aggregation, and grazing and decreases with silica dissolution from DetSi and DetZ2Si.

$$\begin{array}{}\text{(A24)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetSi}\right)& =(\underset{\text{Diatom excretion}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}}}+\underset{\text{Aggregation}}{\underbrace{\mathrm{Agg})}}\cdot \mathrm{DiaSi}\\ & +\underset{\text{Sloppy feeding}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}}\cdot q^{\mathrm{Si}:\mathrm{N}}}}-\underset{\text{Remineralization to DSi}}{\underbrace{{\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}\cdot \mathrm{DetSi}}}\end{array}\text{(A25)}& {\displaystyle }\mathrm{SMS}\left(\mathrm{DetZ}\mathrm{2}\mathrm{Si}\right)=\underset{\text{Sloppy feeding}}{\underbrace{G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}\cdot q^{\mathrm{Si}:\mathrm{N}}}}-\underset{\text{Remineralization to DSi}}{\underbrace{{\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Si}}}\end{array}$$

The state variables DiaSi, DetSi, and DetZ2Si are listed in Table A1. The constant excretion rate (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}$; Table A8) is downregulated by the limiter factor $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$ (Eq. A45) when the N:C ratio becomes too high. The remineralization rates (${\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}$), the aggregation (Agg), and the grazing on diatoms ($G_{\mathrm{dia}}^{\mathrm{zoo}}$ and $G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}$) are calculated in Eqs. (A35), (A42), and (A55), respectively. The intracellular ratio between diatom silicon and carbon is defined as $q^{\mathrm{Si}:\mathrm{N}}={\mathrm{PhySi}}_{\mathrm{dia}}/{\mathrm{PhyN}}_{\mathrm{dia}}$.

CaCO₃

The coccolithophore fraction of small phytoplankton loses biogenic CaCO₃ to the detrital CaCO₃ pool along with excretion, aggregation, respiration, and grazing. Dissolution of CaCO₃ leads to an increase in DIC and alkalinity.

$$\begin{array}{}\text{(A26)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetCalc}\right)& =\underset{\text{Small phytoplankton, excretion}}{\underbrace{{\mathit{ϵ}}_{\mathrm{C}}^{\mathrm{phy}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot \mathrm{PhyCalc}}}\\ & +(\underset{\text{Aggregation}}{\underbrace{\mathrm{Agg}}}+\underset{\mathrm{Respiration}}{\underbrace{r_{\mathrm{small}}}})\cdot \mathrm{PhyCalc}\\ & +\underset{\text{Grazing loss}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}}}\\ & -\underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution in guts}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}\cdot {\mathrm{Diss}}_{\text{calc-guts}}}}\\ & -\underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution, slow-sinking detritus}}{\underbrace{{\mathrm{Diss}}_{\mathrm{calc}}\cdot \mathrm{DetCalc}}}\end{array}\text{(A27)}& {\displaystyle }\begin{array}{rl}\mathrm{SMS}\left(\mathrm{DetZ}\mathrm{2}\mathrm{Calc}\right)& =\underset{\text{Grazing loss}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}}}\\ & -\underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution in guts}}{\underbrace{G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}\cdot q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}\cdot {\mathrm{Diss}}_{\text{calc-guts}}}}\\ & -\underset{{\mathrm{CaCO}}_{\mathrm{3}}\text{dissolution, fast-sinking detritus}}{\underbrace{{\mathrm{Diss}}_{\mathrm{calc}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Calc}}}\end{array}\end{array}$$

The state variables PhyCalc, DetCalc, and DetZ2Calc are listed in Table A1. The constant excretion rate (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}$; Table A8) is downregulated by the limiter factor $f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$ (Eq. A45) when the N:C ratio becomes too high. The respiration (r_small), the aggregation (Agg), and the grazing on small phytoplankton ($G_{\mathrm{small}}^{\mathrm{zoo}}$ and $G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}$) are calculated in Eqs. (A38), (A42), and (A54), respectively. The ratio $q_{\mathrm{small}}^{{\mathrm{CaCO}}_{\mathrm{3}}:\mathrm{N}}=\mathrm{PhyCalc}/{\mathrm{PhyN}}_{\mathrm{small}}$.

Calcite dissolution

As the detritus calcite sinks through the water column, it is subject to dissolution. We follow Yamanaka and Tajika (1996), assuming an exponential decrease in the CaCO₃ flux with depth. As we also assume an increasing sinking speed of small detritus with depth, following Kriest and Oschlies (2008), the dissolution rate is scaled with the sinking velocity.

$$\begin{array}{}\text{(A28)}& {\mathrm{Diss}}_{\mathrm{calc}}={\mathrm{Diss}}_{\text{calc\_rate}}\cdot w_{\det }{\mathrm{Diss}}_{\mathrm{calc}\mathrm{2}}={\mathrm{Diss}}_{\text{calc\_rate}}\end{array}$$

Diss_calc and Diss_calc2 are the dissolution rate constants for slow- and fast-sinking detritus classes (Table A5). The reference dissolution rate (Diss_{calc_rate}; Table A8) is based on a length scale of 3500 m and velocity of 20 m d⁻¹. The sinking speed at depth z (w_det; Table A5) is calculated as follows:

$$\begin{array}{}\text{(A29)}& w_{\det }=\mathrm{0.0288}\cdot z+w_{\mathrm{0}}.\end{array}$$

Here, z denotes the depth, and w₀ is the sinking speed at the ocean surface (Table A3). The dissolution rate for a fast-sinking detritus class (Diss_calc2) is assumed to be constant throughout the water column and is set to the value of Diss_{calc_rate} (Table A8).

A1.6 Dissolved oxygen (Oxy)

Oxy concentration increases with the carbon fixation by primary producers. It decreases with the respiration of phytoplankton and zooplankton and the remineralization of dissolved organic carbon. In addition, sea–air flux of O₂ leads to an exchange of oxygen with the atmosphere, depending on the partial pressure difference in the O₂ between ocean and atmosphere. This exchange is treated separately as a boundary condition. The partial pressure of surface ocean O₂ is computed using the mocsy 2.0 routines (Orr and Epitalon, 2015).

$$\begin{array}{}\text{(A30)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{Oxy}\right)& =\underset{\text{Net production, small phytoplankton}}{\underbrace{\left(P_{\mathrm{small}}-r_{\mathrm{small}}\right)\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}\\ & +\underset{\text{Net production, diatom}}{\underbrace{\left(P_{\mathrm{dia}}-r_{\mathrm{dia}}\right)\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}-\underset{\text{Remineralization of DOC}}{\underbrace{{\mathit{\rho }}_{\mathrm{DOC}}\cdot f_{\mathrm{T}}\cdot \mathrm{DOC}}}\\ & -\underset{\text{Respiration, small zoo.}}{\underbrace{r_{\mathrm{zoo}}\cdot \mathrm{ZooC}}}-\underset{\text{Respiration, macrozoo.}}{\underbrace{r_{\mathrm{zoo}\mathrm{2}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{C}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, DOC, ZooC, and Zoo2C are listed in Table A1. Respiration rate constants of small phytoplankton (r_small), diatoms (r_dia), and zooplankton groups (r_zoo and r_zoo2) are computed in Sects. A3.2 and A4.1, respectively. Photosynthesis terms (P_small and P_dia) are calculated in Eq. (A36). The remineralization rate constant ρ_DOC is listed in Table A8 and the temperature dependency (f_T) is given in Eq. (A33).

A1.7 Dissolved organic material

Dissolved organic matter in our model is a representation of the semi-labile fraction only; the refractory and labile fractions are not included.

Dissolved organic nitrogen (DON)

DON is produced via nitrogen excretion by phytoplankton, by zooplankton, and by the degradation of detrital nitrogen. DON is turned into DIN by remineralization, which is the only sink term.

$$\begin{array}{}\text{(A31)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{DON}\right)& =\underset{\text{Excretion, small phytoplankton}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyN}}_{\mathrm{small}}}}+\underset{\text{Excretion, diatom}}{\underbrace{{\mathit{ϵ}}_{\mathrm{dia}}^{\mathrm{N}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyN}}_{\mathrm{dia}}}}\\ & +\underset{\text{Excretion, small zoo.}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{N}}\cdot \mathrm{ZooN}}}+\underset{\text{Excretion, macrozoo.}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{N}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{N}}}\\ & +\underset{\text{Detritus degradation, slow-sinking}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetN}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetN}}}\\ & +\underset{\text{Detritus degradation, fast sinking}}{\underbrace{{\mathit{\rho }}_{\mathrm{Det}\mathrm{2}\mathrm{ZN}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{N}}}-\underset{\mathrm{Remineralization}}{\underbrace{{\mathit{\rho }}_{\mathrm{DON}}\cdot f_{\mathrm{T}}\cdot \mathrm{DON}}}\end{array}\end{array}$$

The state variables PhyN_small, PhyN_dia, ZooN, DetN, Zoo2N, DetZ2N, and DON are listed in Table A1. The constant excretion rate of nitrogen from phytoplankton and zooplankton classes (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{N}}$, ${\mathit{ϵ}}_{\mathrm{dia}}^{\mathrm{N}}$, ${\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{N}}$, and ${\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{N}}$), the degradation rate of detritus (ρ_DetN and ρ_DetZ2N) and the remineralization rate of DON (ρ_DON) are listed in Table A8. The constant excretion rate of phytoplankton is downregulated by the limiter function ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$ and $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$; Eq. A45) when the N:C ratio becomes too high. The temperature dependency f_T is calculated in Eq. (A33).

Dissolved organic carbon (DOC)

DOC is produced via carbon excretion by phytoplankton, by zooplankton, and by the degradation of detrital carbon. DOC is turned into DIC by remineralization, which is the only sink term.

$$\begin{array}{}\text{(A32)}& \begin{array}{rl}\mathrm{SMS}\left(\mathrm{DOC}\right)& =\underset{\text{Excretion, small phytoplankton}}{\underbrace{{\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyC}}_{\mathrm{small}}}}+\underset{\text{Excretion, diatom}}{\underbrace{{\mathit{ϵ}}_{\mathrm{dia}}^{\mathrm{C}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\mathrm{PhyC}}_{\mathrm{dia}}}}\\ & +\underset{\text{Excretion, small zoo.}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{C}}\cdot \mathrm{ZooC}}}+\underset{\text{Excretion, macrozoo.}}{\underbrace{{\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}}\cdot \mathrm{Zoo}\mathrm{2}\mathrm{C}}}\\ & +\underset{\text{Detritus degradation, slow sinking}}{\underbrace{{\mathit{\rho }}_{\mathrm{DetC}}\cdot f_{\mathrm{T}}\cdot \mathrm{DetC}}}+\underset{\text{Detritus degradation, fast sinking}}{\underbrace{{\mathit{\rho }}_{\mathrm{Det}\mathrm{2}\mathrm{C}}\cdot f_{\mathrm{T}}\cdot \mathrm{Det}\mathrm{2}\mathrm{C}}}\\ & -\underset{\mathrm{Remineralization}}{\underbrace{{\mathit{\rho }}_{\mathrm{DOC}}\cdot f_{\mathrm{T}}\cdot \mathrm{DOC}}}\end{array}\end{array}$$

The state variables PhyC_small, PhyC_dia, ZooC, DetC, Zoo2C, Det2C, and DOC are listed in Table A1. The constant excretion rate of nitrogen from phytoplankton and zooplankton classes (${\mathit{ϵ}}_{\mathrm{phy}}^{\mathrm{C}}$, ${\mathit{ϵ}}_{\mathrm{dia}}^{\mathrm{C}}$, ${\mathit{ϵ}}_{\mathrm{zoo}}^{\mathrm{C}}$, and ${\mathit{ϵ}}_{\mathrm{zoo}\mathrm{2}}^{\mathrm{C}}$), the degradation rate of detritus (ρ_DetC and ρ_Det2C), and the remineralization rate of DOC (ρ_DOC) are listed in Table A8. The constant excretion rate of phytoplankton is downregulated by the limiter factor ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$ and $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$; Eq. A45) when the N:C ratio becomes too high. Temperature dependency f_T is calculated in Eq. (A33).

A2 Temperature dependence of rates

Arrhenius function

Most metabolic processes are faster at higher temperatures. This temperature dependence is defined relative to a reference temperature.

$$\begin{array}{}\text{(A33)}& f_{\mathrm{T}}=\exp \left(-\mathrm{4500}\cdot \left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{ref}}}}\right)\right)\end{array}$$

T and T_ref are the local and reference temperature in Kelvin, respectively (Table A6).

Macrozooplankton grazing

Macrozooplankton grazing is temperature dependent. A dimensionless exponential temperature function (Butzin and Pörtner, 2016) is used for the parameterization of the temperature dependency (f_Tzoo2; Table A5). Specifically, the following parameterization provides an optimum curve with a maximum at 0.5 ^∘C, as described in Karakuş et al. (2021).

$$\begin{array}{}\text{(A34)}& f_{\mathrm{Tzoo}\mathrm{2}}={\displaystyle \frac{\exp \left(\frac{Q_a}{T_r}-\frac{Q_a}{T}\right)}{\mathrm{1}+\exp \left(\frac{Q_h}{T_h}-\frac{Q_h}{T}\right)}}\end{array}$$

T_r is the intrinsic optimum temperature for development, and T_h is the temperature above which inhibitive processes dominate. Q_a and Q_h are the temperatures for the uninhibited and inhibited reaction kinetics, respectively (Table A9). T is the local temperature in Kelvin.

Silicon dissolution

The temperature-dependent dissolution rate of silicon (${\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}$; Table A5) is calculated following Maerz et al. (2020) but with a minimum dissolution rate.

$$\begin{array}{}\text{(A35)}& {\mathit{\rho }}_{\mathrm{Si}}^{\mathrm{T}}=max\left(\mathrm{0.023}\cdot {\mathrm{2.6}}^{\frac{T-\mathrm{10}}{\mathrm{10}}},{\mathit{\rho }}_{\mathrm{Si}}\right)\end{array}$$

T is the local temperature in degrees Celsius. The minimum dissolution rate (ρ_Si) is listed in Table A8.

A3 Phytoplankton processes

Phytoplankton growth equations are based on Geider et al. (1998), with small modifications for diatom silicon uptake, following Hohn (2009).

A3.1 Photosynthesis

The rate of the carbon-specific (C-specific from now on) photosynthesis for phytoplankton (P_small and P_dia) is parameterized as follows:

$$\begin{array}{}\text{(A36)}& \begin{array}{rl}{P}_{\mathrm{small}}& =P_{\max }^{\mathrm{small}}\cdot \left(\mathrm{1.0}-\exp \left({\displaystyle \frac{-{\mathit{\alpha }}_{\mathrm{small}}\cdot q^{\mathrm{Chl}:\mathrm{C}}\cdot \mathrm{PAR}}{P_{\max }^{\mathrm{small}}}}\right)\right),\\ & P_{\mathrm{dia}}=P_{\max }^{\mathrm{dia}}\cdot \left(\mathrm{1.0}-\exp \left({\displaystyle \frac{-{\mathit{\alpha }}_{\mathrm{dia}}\cdot q^{\mathrm{Chl}:\mathrm{C}}\cdot \mathrm{PAR}}{P_{\max }^{\mathrm{dia}}}}\right)\right).\end{array}\end{array}$$

The light-harvesting efficiency (α_small and α_dia) per chlorophyll is listed in Table A7. PAR is the photosynthetically available radiation (Table A5). The intracellular Chl to C ratio (q^Chl:C) is defined as PhyChl$/\mathrm{PhyC}$ and varies as a result of photoacclimation. The apparent maximum photosynthetic rate ($P_{\max }^{\mathrm{small}}$ and $P_{\max }^{\mathrm{dia}}$) is defined below.

$$\begin{array}{}\text{(A37)}& \begin{array}{rl}{P}_{\max }^{\mathrm{small}}& ={\mathit{\mu }}_{\text{C small}}^{\max }\cdot min\left(f_{\text{lim, small}}^{\mathrm{Fe}},f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmin}}\right)\cdot f_{\mathrm{T}},\\ & P_{\max }^{\mathrm{dia}}={\mathit{\mu }}_{\text{C, dia}}^{\max }\cdot min\left(f_{\text{lim, dia}}^{\mathrm{Fe}},f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmin}},f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}\right)\cdot f_{\mathrm{T}}\end{array}\end{array}$$

The value of ${\mathit{\mu }}_{\text{C small}}^{\max }$ and ${\mathit{\mu }}_{\text{C, dia}}^{\max }$ is listed in Table A7. The limitation terms ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmin}}$, $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmin}}$, $f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}$, $f_{\text{lim, small}}^{\mathrm{Fe}}$, and $f_{\text{lim, dia}}^{\mathrm{Fe}}$) are presented in Sect. A3.6, and the temperature dependency (f_T) is calculated in Eq. (A33).

A3.2 Respiration

The phytoplankton respiration rate (r_small and r_dia; Table A5) is calculated as a base respiration plus a second term proportional to N assimilation and as a measure of biosynthesis:

$$\begin{array}{}\text{(A38)}& \begin{array}{rl}{r}_{\mathrm{small}}& =\underset{\mathrm{Maintenance}}{\underbrace{{\mathrm{res}}_{\mathrm{small}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}}}+\underset{\text{N assim}}{\underbrace{\mathit{\zeta }\cdot V_{\mathrm{small}}^{\mathrm{N}}}},\\ & r_{\mathrm{dia}}=\underset{\mathrm{Maintenance}}{\underbrace{{\mathrm{res}}_{\mathrm{dia}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}}}+\underset{\text{N assim}}{\underbrace{\mathit{\zeta }\cdot V_{\mathrm{dia}}^{\mathrm{N}}}}.\end{array}\end{array}$$

The values for the maintenance respiration rate (res_small and res_dia) and the cost of biosynthesis (ζ) are listed in Table A7. Si assimilation is assumed to be inexpensive, so it is not included as additional cost in the respiration (Hohn, 2009). The limiter function ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$ and $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$) is described in Eq. (A45), and the N assimilation rate ($V_{\mathrm{small}}^{\mathrm{N}}$ and $V_{\mathrm{dia}}^{\mathrm{N}}$) is calculated in Eq. (A40).

A3.3 Chlorophyll a synthesis

The chlorophyll synthesis rate ($S_{\mathrm{small}}^{\mathrm{chl}}$ and $S_{\mathrm{dia}}^{\mathrm{chl}}$; Table A5) is proportional to N assimilation, with the proportionality factor varying as a function of the C-specific photosynthesis rate, relative to the maximum possible photosynthetic rate at the current Chl:C ratio of the cell, which depends on photosynthetically available radiation and light-harvesting efficiency.

$$\begin{array}{}\text{(A39)}& \begin{array}{rl}{S}_{\mathrm{small}}^{\mathrm{chl}}& =V_{\mathrm{small}}^{\mathrm{N}}\cdot q_{\text{max, small}}^{\mathrm{Chl}:\mathrm{N}}\cdot min\left(\mathrm{1},{\displaystyle \frac{P_{\mathrm{small}}}{{\mathit{\alpha }}_{\mathrm{small}}\cdot q^{\mathrm{Chl}:\mathrm{C}}\cdot \mathrm{PAR}}}\right),\\ & S_{\mathrm{dia}}^{\mathrm{chl}}=V_{\mathrm{dia}}^{\mathrm{N}}\cdot q_{\text{max, dia}}^{\mathrm{Chl}:\mathrm{N}}\cdot min\left(\mathrm{1},{\displaystyle \frac{P_{\mathrm{dia}}}{{\mathit{\alpha }}_{\mathrm{dia}}\cdot q^{\mathrm{Chl}:\mathrm{C}}\cdot \mathrm{PAR}}}\right)\end{array}\end{array}$$

The N assimilation ($V_{\mathrm{small}}^{\mathrm{N}}$ and $V_{\mathrm{dia}}^{\mathrm{N}}$) is computed in Eq. (A40). The conversion factor of the maximum Chl:N ratio ($q_{\text{max, small}}^{\mathrm{Chl}:\mathrm{N}}$ and $q_{\text{max, dia}}^{\mathrm{Chl}:\mathrm{N}}$) and the light-harvesting efficiency values (α_small and α_dia) are listed in Table A7. The C-specific photosynthesis (P_small and P_dia) is given in Eq. (A36). PAR is the photosynthetically available radiation (Table A5), and the intracellular Chl to C ratio (q^Chl:C) is defined as PhyChl$/\mathrm{PhyC}$.

A3.4 N and Si assimilation

Nitrogen

The C-specific N assimilation rate is a function of the maximum rate of C-specific photosynthesis and DIN concentration. N assimilation depends on the DIN concentration in seawater via Michaelis–Menten kinetics. The maximum photosynthetic rate (P_max) is converted to nitrogen units by multiplication with an optimal N:C uptake ratio (σ). Nitrogen uptake rates are further affected by the intracellular nitrogen status q through the limiter functions $f_{\text{lim}}^{\mathrm{N}:\mathrm{Cmax}}$ (see Eq. A45). Nitrogen assimilation is downregulated at high intracellular N:C ratio.

$$\begin{array}{}\text{(A40)}& \begin{array}{rl}{V}_{\mathrm{small}}^{\mathrm{N}}& =V_{\mathrm{cm}}^{\mathrm{small}}\cdot P_{\max }^{\mathrm{small}}\cdot {\mathit{\sigma }}_{\mathrm{N}:\mathrm{C}}^{\mathrm{small}}\cdot f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\displaystyle \frac{\mathrm{DIN}}{K_{\mathrm{small}}^{\mathrm{N}}+\mathrm{DIN}}},\\ & V_{\mathrm{dia}}^{\mathrm{N}}=V_{\mathrm{cm}}^{\mathrm{dia}}\cdot P_{\max }^{\mathrm{dia}}\cdot {\mathit{\sigma }}_{\mathrm{N}:\mathrm{C}}^{\mathrm{dia}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\displaystyle \frac{\mathrm{DIN}}{K_{\mathrm{dia}}^{\mathrm{N}}+\mathrm{DIN}}}\end{array}\end{array}$$

$V_{\mathrm{cm}}^{\mathrm{small}}$, $V_{\mathrm{cm}}^{\mathrm{dia}}$, ${\mathit{\sigma }}_{\mathrm{N}:\mathrm{C}}^{\mathrm{small}}$, ${\mathit{\sigma }}_{\mathrm{N}:\mathrm{C}}^{\mathrm{dia}}$, $K_{\mathrm{small}}^{\mathrm{N}}$, and $K_{\mathrm{dia}}^{\mathrm{N}}$ are listed in Table A7. The maximum rate of photosynthesis ($P_{\max }^{\mathrm{small}}$ and $P_{\max }^{\mathrm{dia}}$) is given in Eqs. (A37). $f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}$ and $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$ are described in Eq. (A45). DIN corresponds to the in situ concentration.

Silicon

The building of a silica frustule of diatoms requires silicate uptake. The C-specific Si assimilation rate (V^Si) is a function of a factor for C-specific N uptake ($V_{\mathrm{cm}}^{\mathrm{dia}}$), a rate constant of C-specific photosynthesis (${\mathit{\mu }}_{\mathrm{C},\mathrm{dia}}^{\max }$), the maximum Si:C uptake ratio (σ_Si:C), and DSi concentration. Temperature and the intracellular Si:C and N:C ratios (via the limiter functions f_lim) further regulate the Si assimilation.

$$\begin{array}{}\text{(A41)}& V^{\mathrm{Si}}=V_{\mathrm{cm}}^{\mathrm{dia}}\cdot {\mathit{\mu }}_{\text{C, dia}}^{\max }\cdot f_{\mathrm{T}}\cdot {\mathit{\sigma }}_{\mathrm{Si}:\mathrm{C}}\cdot f_{\lim }^{\mathrm{Si}:\mathrm{Cmax}}\cdot f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}\cdot {\displaystyle \frac{\mathrm{DSi}}{K_{\mathrm{Si}}+\mathrm{DSi}}}\end{array}$$

The scaling factor for the N uptake ($V_{\mathrm{cm}}^{\mathrm{dia}}$), the maximum rate constant of C-specific photosynthesis (${\mathit{\mu }}_{\text{C, dia}}^{\max }$), the uptake ratio of the maximum Si:C (σ_Si:C), and half-saturation constant for silicate uptake (K_Si) are listed in Table A7. The temperature dependency (f_T) is computed in Eq. (A33). The limitation by the intracellular ratios N:C and Si:C ($f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}$ and $f_{\lim }^{\mathrm{Si}:\mathrm{Cmax}}$) is described in Eqs. (A45) and (A46), respectively. DSi corresponds to the in situ concentration.

A3.5 Aggregation loss

The aggregation rate (Agg; Table A5) is proportional to the concentration of small phytoplankton, diatoms, and detritus. The effect of increased stickiness of diatoms under nutrient limitation (Waite et al., 1992; Aumont et al., 2015) is taken into account by multiplying the diatom biomass with ($\mathrm{1}-q_{\lim }^{\mathrm{dia}}$). When the nutrient limitation is high (i.e. low $q_{\lim }^{\mathrm{dia}}$), the aggregation rate increases in the model.

$$\begin{array}{}\text{(A42)}& {\displaystyle }\begin{array}{rl}\mathrm{Agg}& =\mathit{\varphi }{}_{\mathrm{phy}}\cdot \left({\mathrm{PhyN}}_{\mathrm{small}}+(\mathrm{1}-q_{\lim }^{\mathrm{dia}})\cdot {\mathrm{PhyN}}_{\mathrm{dia}}\right)\\ & +{\mathit{\varphi }}_{\det }\cdot \left(\mathrm{DetN}+\mathrm{DetZ}\mathrm{2}\mathrm{N}\right)\end{array}\text{(A43)}& {\displaystyle }{q}_{\lim }^{\mathrm{dia}}=min\left(f_{\text{lim, dia}}^{\mathrm{Fe}},f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmin}},f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}\right)\end{array}$$

The state variables PhyN_small, PhyN_dia, DetN, and DetZ2N are described in Table A1. The values of the maximum aggregation loss parameters (ϕ_phy and ϕ_det) are listed in Table A3. The limitation terms ($f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmin}}$, $f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}$, and $f_{\text{lim, dia}}^{\mathrm{Fe}}$) are presented below (Sect. A3.6).

A3.6 Nutrient limitation

The metabolic processes such as C-specific photosynthesis, respiration rate, and excretion losses are treated as functions of the intracellular nitrogen status (i.e. N:C ratios q), following Geider et al. (1998). Intracellular ratios between nutrients and carbon limit the uptake of nitrogen and silicon, which is modelled via a non-linear function, as in Schourup-Kristensen et al. (2014).

$$\begin{array}{}\text{(A44)}& f_{\lim }(\mathit{\theta },q_{\mathrm{1}},q_{\mathrm{2}})=\mathrm{1}-\exp (-\mathit{\theta }(\left|\mathrm{\Delta }q\right|-\mathrm{\Delta }q{)}^{\mathrm{2}})\end{array}$$

Here, $\mathrm{\Delta }q=q_{\mathrm{1}}-q_{\mathrm{2}}$ is the difference between the current intracellular nutrient:C quota and a prescribed maximum or minimum quota. The dimensionless constant θ controls the limitation.

$f_{\lim }^{\mathrm{N}:\mathrm{Cmax}}$

The limiter $f_{\lim }^{\mathrm{N}:\mathrm{Cmax}}$ downregulates the metabolic processes such as nitrogen and Si assimilation, excretion, and maintenance respiration of phytoplankton when the intracellular nitrogen quota (q^N:C) becomes too high. $f_{\lim }^{\mathrm{N}:\mathrm{Cmax}}$ is one when the current $q^{\mathrm{N}:\mathrm{C}}<\mathrm{0.151}$ (i.e. Redfield ratio; 16N:106C) and zero for $q^{\mathrm{N}:\mathrm{C}}>\mathrm{0.2}$ (i.e. 21.2N:106C). It determines the end of the uptake of nitrogen and silicon in assimilation processes and the cease of carbon and nitrogen release during the respiration and excretion of DON and/or DOC and the CaCO₃ processes of phytoplankton (see the discussion of CaCO₃ in Sect. A1.5).

$$\begin{array}{}\text{(A45)}& \begin{array}{rl}{f}_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmax}}& =f_{\lim }({\mathit{\theta }}_{\max }^{\mathrm{N}},q_{\mathrm{small}}^{\mathrm{N}:\mathrm{C}},q_{\mathrm{small}}^{\mathrm{N}:\mathrm{Cmax}}),\\ & f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmax}}=f_{\lim }({\mathit{\theta }}_{\max }^{\mathrm{N}},q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{C}},q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{Cmax}})\end{array}\end{array}$$

The limitation function for the quota regulation is calculated with Eq. (A44). $q_{\mathrm{small}}^{\mathrm{N}:\mathrm{C}}$ and $q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{C}}$ form the current intracellular nitrogen quota for small phytoplankton and diatoms, respectively. Dimensionless constants ${\mathit{\theta }}_{\max }^{\mathrm{N}}$, $q_{\mathrm{small}}^{\mathrm{N}:\mathrm{Cmax}}$, and $q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{Cmax}}$ are listed in Table A6.

$f_{\lim }^{\mathrm{Si}:\mathrm{Cmax}}$

The limiter $f_{\lim }^{\mathrm{Si}:\mathrm{Cmax}}$ downregulates the Si assimilation of diatoms when the intracellular silicon quota (Si:C) becomes too high. $f_{\lim }^{\mathrm{Si}:\mathrm{Cmax}}$ is one when the current $q^{\mathrm{N}:\mathrm{C}}<\mathrm{0.76}$ and zero for $q^{\mathrm{N}:\mathrm{C}}>\mathrm{0.8}$. It determines the end of the uptake of silicon in the assimilation processes. The limiter function is described in Eq. (A44) and is calculated as follows:

$$\begin{array}{}\text{(A46)}& f_{\lim }^{\mathrm{Si}:\mathrm{Cmax}}=f_{\lim }({\mathit{\theta }}_{\max }^{\mathrm{Si}},q^{\mathrm{Si}:\mathrm{C}},q^{\mathrm{Si}:\mathrm{Cmax}})\end{array}$$

Dimensionless constants ${\mathit{\theta }}_{\max }^{\mathrm{Si}}$ and q^Si:Cmax are listed in Table A6.

$f_{\lim }^{\mathrm{Si}:\mathrm{Cmin}}$

Carbon fixation and aggregation loss in diatoms are further downregulated by a factor ($f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}$; see Eq. A44) when the intracellular silicon quota (q^Si:C) approaches a minimum value (q^Si:Cmin), mimicking the arrest of cellular division at low cellular Si (Claquin et al., 2002). $f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}$ is zero when the current $q^{\mathrm{Si}:\mathrm{C}}<\mathrm{0.04}$ and one for $q^{\mathrm{Si}:\mathrm{C}}>\mathrm{0.08}$.

$$\begin{array}{}\text{(A47)}& f_{\text{lim, dia}}^{\mathrm{Si}:\mathrm{Cmin}}=f_{\lim }({\mathit{\theta }}_{\min }^{\mathrm{Si}},q^{\mathrm{Si}:\mathrm{Cmin}},q^{\mathrm{Si}:\mathrm{C}})\end{array}$$

The dimensionless constants ${\mathit{\theta }}_{\min }^{\mathrm{Si}}$ and q^Si:Cmin are listed in Table A6.

$f_{\lim }^{\mathrm{Fe}}$

Growth limitation by iron is modelled with Michaelis–Menten kinetics, implicitly assuming that all dissolved iron is ultimately bioavailable.

$$\begin{array}{}\text{(A48)}& f_{\text{lim, small}}^{\mathrm{Fe}}={\displaystyle \frac{\mathrm{DFe}}{K_{\mathrm{small}}^{\mathrm{Fe}}+\mathrm{DFe}}},f_{\text{lim, dia}}^{\mathrm{Fe}}={\displaystyle \frac{\mathrm{DFe}}{K_{\mathrm{dia}}^{\mathrm{Fe}}+\mathrm{DFe}}}\end{array}$$

The variable DFe is listed in Table A1. The half-saturation constants ($K_{\mathrm{small}}^{\mathrm{Fe}}$ and $K_{\mathrm{dia}}^{\mathrm{Fe}}$) are given in Table A6.

$f_{\lim }^{\mathrm{N}:\mathrm{Cmin}}$

In addition to iron limitation, photosynthesis is limited by nitrogen in small phytoplankton and diatoms using the Eq. (A44). Nitrogen limitation ($f_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmin}}$ and $f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmin}}$) is described as a function of the intracellular nitrogen quota ($q_{\mathrm{small}}^{\mathrm{N}:\mathrm{C}}$ and $q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{C}}$) with growth ending at a minimum quota ($q_{\mathrm{small}}^{\mathrm{N}:\mathrm{Cmin}}$ and $q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{Cmin}}$).

$$\begin{array}{}\text{(A49)}& \begin{array}{rl}{f}_{\text{lim, small}}^{\mathrm{N}:\mathrm{Cmin}}& =f_{\lim }({\mathit{\theta }}_{\min }^{\mathrm{N}},q_{\mathrm{small}}^{\mathrm{N}:\mathrm{Cmin}},q_{\mathrm{small}}^{\mathrm{N}:\mathrm{C}}),\\ & f_{\text{lim, dia}}^{\mathrm{N}:\mathrm{Cmin}}=f_{\lim }({\mathit{\theta }}_{\min }^{\mathrm{N}},q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{Cmin}},q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{C}})\end{array}\end{array}$$

Dimensionless constants ${\mathit{\theta }}_{\min }^{\mathrm{N}}$, $q_{\mathrm{small}}^{\mathrm{N}:\mathrm{Cmin}}$, and $q_{\mathrm{dia}}^{\mathrm{N}:\mathrm{Cmin}}$ are listed in Table A6.

A4 Zooplankton processes

A4.1 Zooplankton respiration

Small zooplankton

When the intracellular C:N ratio in zooplankton exceeds the Redfield ratio, a temperature-dependent respiration (r_zoo; Table A5) is assumed to drive it back with a timescale τ.

$$\begin{array}{}\text{(A50)}& r_{\mathrm{zoo}}={\displaystyle \frac{q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}-q_{\mathrm{standard}}^{\mathrm{C}:\mathrm{N}}}{\mathit{\tau }}}\cdot f_{\mathrm{T}}\end{array}$$

The timescale for respiration (τ) is listed in Table A7. The temperature dependence (f_T) is calculated in Eq. (A33). The ratios are defined as $q_{\mathrm{zoo}}^{\mathrm{C}:\mathrm{N}}=\mathrm{ZooC}/\mathrm{ZooN}$ and $q_{\mathrm{Standard}}^{\mathrm{C}:\mathrm{N}}=\mathrm{106}/\mathrm{16}$.

Macrozooplankton

The daily respiration rate constant of macrozooplankton (r_zoo2; Table A5) is modelled following Karakuş et al. (2021).

$$\begin{array}{}\text{(A51)}& r_{\mathrm{zoo}\mathrm{2}}=R_{\mathrm{s}}\cdot (\mathrm{1}+R_{\mathrm{f}}+R_{\mathrm{a}})\end{array}$$

The standard respiration rate (R_s) is listed in Table A3. The feeding activity factor (R_f; Table A5) is defined as the ratio of grazing flux to carbon biomass of macrozooplankton, which increases linearly from zero to one for ratio between 0 % and 10 % and is one otherwise. The respiration activity factor (R_a; Table A5) defines the reduced macrozooplankton respiration rate in austral or boreal winters with the value of −0.5.

A4.2 Grazing

In REcoM3, there are two zooplankton classes, namely small zooplankton (<2 cm) and macrozooplankton (2–20 cm). The small zooplankton group grazes on small phytoplankton, diatoms, and on fast- and slow-sinking detrital particles. Similarly, while macrozooplankton grazes on both phytoplankton classes and detritus groups, it further grazes on small zooplankton. The total grazing of both zooplankton groups is based on the Holling type III ingestion function as follows:

$$\begin{array}{}\text{(A52)}& G_{\mathrm{tot}}^{\mathrm{zoo}}={\mathit{\xi }}_{\mathrm{zoo}}\cdot {\displaystyle \frac{{\left({\sum }_i{p}_i\cdot N_i\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{\mathrm{zoo}}+{\left({\sum }_i{p}_i\cdot N_i\right)}^{\mathrm{2}}}}\cdot f_{\mathrm{T}}\cdot \mathrm{ZooN}.\end{array}$$

$G_{\mathrm{tot}}^{\mathrm{zoo}}$ ($G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}$) is the total grazing flux which is calculated for small (macro) zooplankton. ZooN (Zoo2N) is listed in Table A1. The maximum grazing rate (ξ_zoo and ξ_zoo2) and the half-saturation constants (σ_zoo and σ_zoo2) are listed in Table A10. The temperature dependency terms (f_T and f_Tzoo2) are given in Eqs. (A33) and (A34). In the model, relative grazing preferences are implemented following Fasham et al. (1990). Variable relative grazing preferences (p_i) are calculated using the nominal preferences for small phytoplankton, diatoms, slow- and fast-sinking detritus, and small zooplankton (Table A10) as follows:

$$\begin{array}{}\text{(A53)}& p_i={\displaystyle \frac{p_i^{\prime }\cdot N_i}{{\sum }_i{p}_i^{\prime }\cdot N_i}}.\end{array}$$

Here, summation i is done over each food source to calculate the relative proportion of the food. Total grazing is used to calculate the grazing of zooplankton groups on individual food source; i.e. small phytoplankton (i=1; PhyN_small), diatoms (i=2; PhyN_dia), and both detritus classes (i=3, DetN; i=4, DetZ2N) and (i=5, ZooN) in the case of macrozooplankton as the ratio of each food source to total food source (G_small, G_dia, G_det, G_detZ2, and G_zoo).

$$\begin{array}{}\text{(A54)}& {\displaystyle }\begin{array}{rl}{G}_{\mathrm{small}}^{\mathrm{zoo}}& =G_{\mathrm{tot}}^{\mathrm{zoo}}\cdot {\displaystyle \frac{p_{\mathrm{small}}\cdot {\mathrm{PhyN}}_{\mathrm{small}}}{{\sum }_i{p}_i\cdot N_i}},\\ & G_{\mathrm{small}}^{\mathrm{zoo}\mathrm{2}}=G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}\cdot {\displaystyle \frac{p_{\mathrm{small}}\cdot {\mathrm{PhyN}}_{\mathrm{small}}}{{\sum }_i{p}_i\cdot N_i}}\end{array}\text{(A55)}& {\displaystyle }\begin{array}{rl}{G}_{\mathrm{dia}}^{\mathrm{zoo}}& =G_{\mathrm{tot}}^{\mathrm{zoo}}\cdot {\displaystyle \frac{p_{\mathrm{dia}}\cdot {\mathrm{PhyN}}_{\mathrm{dia}}}{{\sum }_i{p}_i\cdot N_i}},\\ & G_{\mathrm{dia}}^{\mathrm{zoo}\mathrm{2}}=G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}\cdot {\displaystyle \frac{p_{\mathrm{dia}}\cdot {\mathrm{PhyN}}_{\mathrm{dia}}}{{\sum }_i{p}_i\cdot N_i}}\end{array}\text{(A56)}& {\displaystyle }\begin{array}{rl}{G}_{\det }^{\mathrm{zoo}}& =G_{\mathrm{tot}}^{\mathrm{zoo}}\cdot {\displaystyle \frac{p_{\det }\cdot \mathrm{DetN}}{{\sum }_i{p}_i\cdot N_i}},\\ & G_{\det }^{\mathrm{zoo}\mathrm{2}}=G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}\cdot {\displaystyle \frac{p_{\det }\cdot \mathrm{DetN}}{{\sum }_i{p}_i\cdot N_i}}\end{array}\text{(A57)}& {\displaystyle }\begin{array}{rl}{G}_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}}& =G_{\mathrm{tot}}^{\mathrm{zoo}}\cdot {\displaystyle \frac{p_{\mathrm{detZ}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{N}}{{\sum }_i{p}_i\cdot N_i}},\\ & G_{\mathrm{detZ}\mathrm{2}}^{\mathrm{zoo}\mathrm{2}}=G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}\cdot {\displaystyle \frac{p_{\mathrm{detZ}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{N}}{{\sum }_i{p}_i\cdot N_i}}\end{array}\text{(A58)}& {\displaystyle }{G}_{\mathrm{zoo}}=G_{\mathrm{tot}}^{\mathrm{zoo}\mathrm{2}}\cdot {\displaystyle \frac{p_{\mathrm{ZooN}}\cdot \mathrm{ZooN}}{{\sum }_i{p}_i\cdot N_i}},\end{array}$$

where G_zoo is associated with macrozooplankton grazing on small zooplankton. PhyN_small, PhyN_dia, ZooN, DetN, and DetZ2N are listed in Table A1.

A5 Bottom boundary fluxes

The model contains a benthic layer at the sea floor. Within this benthic layer, the total amounts of organic carbon, organic nitrogen, biogenic silica, and CaCO₃ are modelled.

Loss to benthos

When the slow- and fast-sinking detritus reach the ocean bottom, they continue to sink into the benthic layer with the speeds of w_det (Eq. A29) and w_detZ2 = 200 m d⁻¹, respectively. This results in a detrital flux (BenF_DetN, BenF_DetZ2N, BenF_DetC, BenF_DetZ2C, BenF_DetSi, BenF_DetZ2Si, BenF_DetCalc, and BenF_DetZ2Calc; Table A11) from the water column to the benthos.

$$\begin{array}{}\text{(A59)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetN}}=-w_{\det }\cdot \mathrm{DetN}\text{(A60)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetC}}=-w_{\det }\cdot \mathrm{DetC}\text{(A61)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetSi}}=-w_{\det }\cdot \mathrm{DetSi}\text{(A62)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetCalc}}=-w_{\det }\cdot \mathrm{DetCalc}\text{(A63)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetZ}\mathrm{2}\mathrm{N}}=-w_{\mathrm{detZ}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{N}\text{(A64)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetZ}\mathrm{2}\mathrm{C}}=-w_{\mathrm{detZ}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{C}\text{(A65)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetZ}\mathrm{2}\mathrm{Si}}=-w_{\mathrm{detZ}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Si}\text{(A66)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DetZ}\mathrm{2}\mathrm{Calc}}=-w_{\mathrm{detZ}\mathrm{2}}\cdot \mathrm{DetZ}\mathrm{2}\mathrm{Calc}\end{array}$$

These fluxes increase the total amount of the different benthic state variables. The state variables of DetN, DetC, DetSi, DetCalc, DetZ2N, DetZ2C, DetZ2Si, and DetZ2Calc are described in Table A1.

Input from benthos

The lowermost ocean layer located next to the benthic layer receives remineralized inorganic matter back from the benthos. At the same time, these fluxes reduce the number of the benthic variables. In addition, a sediment flux of Fe from the sediment is calculated from the nitrogen flux but assuming a Fe:N ratio that is higher than in the biomass. This parameterization models a scenario in which the release of iron from the sediment is driven by redox processes, which are ultimately tied to their remineralization of organic matter.

$$\begin{array}{}\text{(A67)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DIN}}={\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{N}}\cdot \mathrm{BenthosN}\text{(A68)}& {\displaystyle }{\mathrm{BenF}}_{\mathrm{DSi}}={\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{Si}}\cdot \mathrm{BenthosSi}\text{(A69)}& {\displaystyle }\begin{array}{rl}{\mathrm{BenF}}_{\mathrm{DIC}}& ={\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{C}}\cdot \mathrm{BenthosC}+{\mathrm{Diss}}_{\mathrm{calc}}\cdot \mathrm{BenthosCalc}\\ & +{\mathrm{Diss}}_{\mathrm{calc}\mathrm{2}}\cdot \mathrm{BenthosCalc}\mathrm{2}\end{array}\text{(A70)}& {\displaystyle }\begin{array}{rl}{\mathrm{BenF}}_{\mathrm{Alk}}& =(\mathrm{1}+\mathrm{1}/\mathrm{16})\cdot {\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{N}}\cdot \mathrm{BenthosN}+\mathrm{2}\cdot {\mathrm{Diss}}_{\mathrm{calc}}\\ & \cdot \mathrm{BenthosCalc}\end{array}\end{array}$$

BenF_DIN, BenF_DSi, BenF_DIC, and BenF_Alk (Table A11) denote the fluxes of DIN, DSi, DIC, and Alk returned into the bottom layer of the ocean. Constant remineralization rates (${\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{N}}$, ${\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{Si}}$, and ${\mathit{\rho }}_{\mathrm{ben}}^{\mathrm{C}}$) are listed in Table A8. The calcite dissolution rates Diss_calc and Diss_calc2 are calculated in Eq. (A28). BenthosN, BenthosSi, BenthosC, and BenthosCalc denote the vertically integrated benthos concentration of dissolved nitrogen, silicate, carbon, and calcium carbonate, respectively (Table A2). The alkalinity of the lowermost ocean layer located next to the benthic layer is changed by the remineralization of DIN, dissolved inorganic phosphate converted from DIN with Redfield ratio, and the dissolution of calcite from the benthos.

Table A1 List of oceanic state variables in REcoM3.

Table A2 List of benthic state variables in REcoM3.

Table A3 Parameters for equations of sources minus sinks.

Table A4 Parameters for iron calculations.

Table A5 Model variables.

Table A6 Parameters for limitation functions.

Table A7 Parameters for phytoplankton processes.

Table A8 Degradation parameters for sources-minus-sinks equations.

Table A9 Parameters for macrozooplankton grazing.

Table A10 Parameters for zooplankton processes.

Table A11 Benthos variables.