Articles | Volume 16, issue 16
Model description paper
30 Aug 2023
Model description paper |  | 30 Aug 2023

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

Özgür Gürses, Laurent Oziel, Onur Karakuş, Dmitry Sidorenko, Christoph Völker, Ying Ye, Moritz Zeising, Martin Butzin, and Judith Hauck

The cycling of carbon in the oceans is affected by feedbacks driven by changes in climate and atmospheric CO2. 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 CO2 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 (CO2) 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 CO2 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−1 (26 % of the total CO2 emissions; Friedlingstein et al.2022b). A similar proportion was taken up by the terrestrial biosphere, amounting to 3.1 ± 0.6 PgC yr−1 (2012–2021), but the total air-to-land CO2 flux is substantially lower because of emissions from land use change, mainly deforestation, that amounted to 1.2 ± 0.7 PgC yr−1 (2012–2021). The ocean carbon sink has grown over the past few decades in response to the near-exponential rise in CO2 emissions (Friedlingstein et al.2022b). While the global ocean carbon sink estimate is assigned an uncertainty of 0.4 PgC yr−1 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 pCO2 and CO2 fluxes (Mongwe et al.2018). Both climate change and rising atmospheric CO2 will feed back on the fraction of CO2 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 CO2 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 CO2 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 CO2 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 CO2 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 Gruber2006). Without the biological carbon pump, atmospheric CO2 would be higher by 200 ppm (Maier-Reimer et al.1996), and perturbations thereof can have large effects on atmospheric CO2 (Kwon et al.2009; Lauderdale and Cael2021), as also known from paleo evidence (Galbraith and Skinner2020).

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 McKinley2021), 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; Bunsen2022). 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 (McWilliams2016), 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 Cael2021). 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 (Danilov2013). 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 (Hohn2009). 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 CO2 forcing and the associated model bias and drift from the experiment with a constant preindustrial (PI) CO2 level. We thus focus on evaluating model aspects with regard to the effects of climate change and CO2 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 1Scheme 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 (Redi1982) and the Gent-McWilliams (GM; Gent and McWilliams1990; Griffies1998) eddy stirring parameterization are applied. Both GM and the Redi (1982) scheme are scaled with horizontal resolution, with a maximum value of 2000 m2 s−1 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−4 m2 s−1 for momentum and 1 × 10−5 m2 s−1 for tracers. Furthermore, the Monin–Obukhov length-dependent vertical mixing parameterization is applied in the surface boundary layer south of 50 S (Timmermann and Beckmann2004).

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öhler2013), as well as the marine iron cycle (e.g. Völker and Tagliabue2015; Tagliabue et al.2016; Ye and Völker2017; 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 2Schematic 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 (C:N:Chl:CaCO3 ratios for small phytoplankton and C:N:Chl: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 Roche1994; Raven1988). 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. CaCO3 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−1 at the surface to 192 m d−1 at 6000 m depth; Kriest and Oschlies2008). The sinking speed of the second group (fast-sinking detritus) is constant throughout the water column (200 m d−1; 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.


  1. The routines for calculating carbonate chemistry and air–sea CO2 exchange used in FESOM1.4–REcoM2, which followed the guidelines provided by the Ocean Carbon-Cycle Model Intercomparison Project (Orr1999), were replaced by the mocsy 2.0 scheme of Orr and Epitalon (2015). While both use the same thermodynamic equilibrium to calculate surface pCO2 and CO2 flux, mocsy 2.0 uses the faster and more accurate algorithm SolveSAPHE (Munhoven2013). 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 O2 flux is calculated using the mocsy 2.0 routines (Orr and Epitalon2015). Photosynthesis, respiration, and remineralization change oxygen with a fixed O2:C ratio, and remineralization does not depend on the O2 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.


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.


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 Yeager2009) 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 1/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−3 (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 1List of simulations performed in this study.

Download Print Version | Download XLSX

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,, 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 CO2 levels (hereafter referred to as A). Atmospheric CO2 mixing ratio (xCO2) values are taken from the Global Carbon Budget (Joos and Spahni2008; Ballantyne et al.2012; Friedlingstein et al.2022a). A second simulation was forced by RYF61 atmospheric reanalysis fields and a preindustrial atmospheric CO2 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 CO2 mixing ratio of 278 ppm. This experiment is referred as D and is used to separate the effects of rising atmospheric CO2 and of climate change on the DIC inventory. Using the simulations A and B, the global ocean anthropogenic CO2 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 CO2 mixing ratio of 278 ppm until the air–sea CO2 reached a quasi-equilibrium state. The Aspinup and Bspinup simulations are a continuation of the pre‐spin-up simulation, with either increasing (Aspinup) or constant (Bspinup) atmospheric CO2, 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 Falkowski1997)(Buitenhuis et al.2010; Moriarty and O'Brien2013)(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 Gruber2021)

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

Download Print Version | Download XLSX

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 CO2 flux primarily under elevating CO2. 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 Vallis2007; 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 3Maps 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 4Vertical representation of the Atlantic meridional overturning circulation (Sv) in simulations A, B, and their difference (Sv).


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


Figure 6Maps 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. (2021ab 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 -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−3). 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−3 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−1 in Schourup-Kristensen et al. (2014) to 168 Tmol Si yr−1 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−1 (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−1 in Schourup-Kristensen et al. (2014) to 85.5 Tmol Si yr−1, which is higher than the range of 21–54 Tmol Si yr−1 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−3 for DIN and 0.47 and 0.54 mmol m−3 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 7Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (0–100 m) concentration of dissolved inorganic nitrogen (mmol m−3; a), dissolved inorganic silicon (mmol m−3; 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 8Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (0–50 m) concentration of dissolved iron (µmol m−3; a), and of the AI-based global reconstruction by Huang et al. (2022b). 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−1 (e.g. Hatta et al.2015), while modelled maxima are >3 nmol L−1. 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ölker2017; 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 Falkowski1997; 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 Marra2022). 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−3. 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−3 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−2 d−1), 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−1, the modelled global total NPP is slightly above the range of the earlier modelling studies (23.7–30.7 PgC yr−1; Schneider et al.2008) and within the range of recent Earth system models (24.5–57.3 PgC yr−1; Séférian et al.2020). It is lower than other satellite-based estimates of 47.3 PgC yr−1 (Behrenfeld and Falkowski1997), 52 PgC yr−1 (Westberry et al.2008), and 48.7–52.5 PgC yr−1 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; Mitchell1992; 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 Arrigo2020). 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 9Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (log10 transformed) chlorophyll a concentration (mg Chl m−3) 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 10Maps of simulated FESOM2.1–REcoM3 (simulation A) vertically integrated net primary production (mgC m−2 d−1) of small phytoplankton (a), diatoms (b), and the sum of both phytoplankton groups (c). The satellite-based Vertically Generalized Production Model (VGPM; Behrenfeld and Falkowski1997) 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 11Latitudinal distribution of vertically integrated and zonally averaged (a) chlorophyll a (mg Chl m−3) and (b) net primary production (mg C m−2 d−1) 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 Falkowski1997; 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-N2-fixing cyanobacteria such as Synechococcus and Prochlorococcus) and nitrogen fixers, and this might contribute to an underestimation of NPP.

Figure 12Maps 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'Brien2013). 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 13Maps 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−2) 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'Brien2013) 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 Falkowski1997)(Westberry et al.2008)(Schneider et al.2008)(Kulk et al.2020)Séférian et al.2020(Schlitzer2004)(Dunne et al.2007)(Henson et al.2011)(Siegel et al.2014)Dunne et al.2007(Tréguer et al.2021)Lee2001; 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 Falkowski1997)(Schlitzer2002; Nevison et al.2012)(Dunne et al.2007)(Dunne et al.2007)

Table 3Global 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.

Download Print Version | Download XLSX

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−3; 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−3). 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 CaCO3 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 CO2 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 CO2 in equilibrium than a low-alkalinity surface ocean. The range of biases is similar (±100 mmol m−3) to other ocean biogeochemical models (e.g. Tjiputra et al.2020; Long et al.2021a).

Figure 14Maps of simulated FESOM2.1–REcoM3 (simulation A) surface (0–100 m) concentration of dissolved inorganic carbon (mmol m−3; a) and alkalinity (mmol m−3; 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 pCO2 and air–sea CO2 flux

We compare the pattern of the temporal mean (2012–2021) surface ocean partial pressure of CO2 (pCO2; Fig. 15) and air–sea CO2 flux (Fig. 16) to the pCO2-based data product of Chau et al. (2022), with a seamless coverage from open ocean to the coasts (Fig. 16). Different pCO2 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 McKinley2021). 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 pCO2 in FESOM2.1–REcoM3 with a direct comparison to surface ocean pCO2 observations from the Surface Ocean CO2 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 15Maps of surface ocean pCO2 (µatm). The top row compares the (a) simulated FESOM2.1–REcoM3 (simulation A) surface partial pressure of CO2 to the (b) pCO2-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 16Maps of air–sea CO2 fluxes (mol C m−2 yr−1). The top row compares the (a) simulated FESOM2.1–REcoM3 (simulation A) CO2 flux to the (b) pCO2-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 pCO2 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 pCO2 product of Chau et al. (2022), the model pCO2 values are overestimated in the subtropical gyres (Fig. 15c). Furthermore, the North Atlantic pCO2 is on average lower than the pCO2 product, and the two data sets also differ in terms of the over- versus undersaturation of pCO2 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 pCO2 observations (e.g. Metzl et al.2006; Gregor et al.2019). FESOM2.1–REcoM3 also simulates very high pCO2 values on the Russian shelves in the Arctic, where hardly any observations exist. Similarly high pCO2 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 17Comparing the annual mean pCO2 (µ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 pCO2 reasonably well compared to SOCAT when accounting for where and when the pCO2 sampling took place (Fig. 17). The annual correlation coefficient and RMSE between the simulated and observed global mean pCO2 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 pCO2 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 pCO2 in the tropics and underestimates pCO2 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 CO2 flux spatial pattern was reasonably reproduced by FESOM2.1–REcoM3, with CO2 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−2 yr−1). Generally, the CO2 flux patterns mirror the pCO2 patterns (Fig. 15) but with the additional imprint of spatial variability in the wind speed. Hence, the CO2 uptake in the subpolar Southern Ocean may appear large compared to pCO2, 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 pCO2 product. While it is well established that the outgassing of CO2 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 CO2 flux and the pCO2-based data product generally mimics pCO2 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 CO2 fluxes with CO2 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 pCO2 products (Resplandy et al.2023). The large mismatch in pCO2 on the Siberian shelves does not show up in the CO2 flux, as sea ice prevents CO2 outgassing throughout most of the year.

Figure 18Time series of simulated annual mean global ocean–atmosphere CO2 flux (in PgC yr−1) 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 CO2 concentration of 278 ppm (dashed lines) and the spin-up under increasing CO2 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 CO2 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 CO2 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 CO2 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 CO2 flux in the FESOM2.1–REcoM3 control simulation amounts to −0.12 PgC yr−1. 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−1 at the end of the simulation. In addition to the bias, the drift is reduced from 0.00264 PgC yr−2 in FESOM1.4–REcoM2 to −0.00011 PgC yr−2 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 CO2 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 19Globally integrated annual air–sea CO2 flux from 10 global ocean biogeochemistry models (GOBMs) and seven pCO2-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−1 (Regnier et al.2022) to the pCO2 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 Gruber2021). The thick black line indicates the model ensemble mean, and the thick blue line shows the mean of the pCO2 product ensemble. Thin dashed lines are from individual GOBMs and pCO2 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 pCO2-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 CO2 concentration and without climate change forcing (simulation B). Although being consistent with the interannual variability, air–sea CO2 fluxes of FESOM1.4 are at the lower end of the range compared to other global ocean biogeochemistry models and pCO2-based estimates. In contrast, starting from the mid-1960s, FESOM2.1 shows a higher CO2 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 CO2 flux is mostly attributed to the level of CO2 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−1 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−1; 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 CO2 flux. Yet, the FESOM2.1–REcoM3 CO2 fluxes have been lower than the mean of the pCO2-based data products since about 2008 and thus affirm the growing discrepancy between global ocean biogeochemistry models and pCO2 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 pCO2 products are statistical models that interpolate and extrapolate sparse pCO2 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 CO2 flux trend (2000–2018) by 20 %–35 %, based on the current pCO2 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 (Sundquist1985; Keppler et al.2020). The DIC inventory is thought to change primarily in response to the rise in atmospheric CO2, 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 CO2. 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 4FESOM2.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 CO2 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).

Download Print Version | Download XLSX

Figure 20Maps of surface (0–10 m; ac) and intermediate depth (300–500 m; df) concentration of simulated FESOM2.1–REcoM3 (simulation A) dissolved O2 (mmol m−3; a, d), World Ocean Atlas 2018 climatology of dissolved O2 (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 O2/N2 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 CO2 flux in the tropics with too little ocean uptake (or too much release) of CO2.

3.4 Oxygen

The simulated distributions of global O2 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 m3 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 O2 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−3. Within the 100–600 m layer, FESOM2.1–REcoM3 performed remarkably well, with simulated values of about 160 ± 105 mmol m−3, which is very close to the observations from WOA 2018 (158 ± 103 mmol m3). 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 CO2 flux of the model simulation and the pCO2-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 CO2 (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 pCO2 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:

(A1) S t = - U S + κ S + SMS ( S ) ,

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 Oschlies2008). 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 A1Maps of simulated FESOM2.1–REcoM3 (simulation A) vertically integrated net primary production (mgC m−2 d−1) 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 A2Conceptual 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:

(A2) SMS ( DIN ) = ρ DON f T DON DON remineralization - V small N PhyC small N assimilation, small phytoplankton - V dia N PhyC dia N assimilation, diatoms .

The state variables DON, PhyCsmall, and PhyCdia are listed in Table A1. The value of the remineralization rate constant (ρDON) is given in Table A8. The temperature dependency of remineralization (fT) is calculated in Eq. (A33). See Sect. A3.4 for details on the carbon-specific nitrogen assimilation rates VsmallN and VdiaN (Table A5).

Dissolved silicic acid (DSi)

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

(A3) SMS ( DSi ) = ρ Si T DetSi Remineralization, slow-sinking detritus + ρ Si T DetZ 2 Si Remineralization, fast-sinking detritus - V Si PhyC dia Si assimilation, diatoms

The state variables PhyCdia, DetSi, and DetZ2Si are listed in Table A1. The temperature-dependent remineralization rate of silicon (ρSiT) and the carbon-specific Si assimilation rate (VSi) 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 (qFe: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.

(A4) SMS ( DFe ) = q Fe : N ( ϵ phy N f lim, small N : Cmax PhyN small Excretion, small phytoplankton + ϵ phy N f lim, dia N : Cmax PhyN dia Excretion, diatoms + ρ DetN f T DetN Remineralization, slow-sinking detritus + ρ DetN f T DetZ 2 N Remineralization, fast-sinking detritus + ϵ zoo N ZooN Excretion, small zooplankton + ϵ zoo 2 N Zoo 2 N Excretion, macrozooplankton - V small N PhyC small N assimilation, small phytoplankton - V dia N PhyC dia N assimilation, diatom ) - κ Fe DetC Fe Scavenging, slow-sinking detritus - κ Fe DetZ 2 C Fe Scavenging, fast-sinking detritus

The state variables PhyCsmall, PhyCdia, PhyNsmall, PhyNdia, DetC, DetN, DetZ2C, DetZ2N, ZooN, and Zoo2N are listed in Table A1. The intracellular Fe:N ratio (qFe:N) and scavenging rate of iron (κFe) are given in Table A4. Excretion rates (ϵphyN, ϵzooN, and ϵzoo2N) and the degradation rate for detritus N (ρDetN) are listed in Table A8. The temperature dependency (fT) is calculated in Eq. (A33). The limitation of intracellular nitrogen (flim, smallN:Cmax; flim, diaN:Cmax) is described in Eq. (A45). Scavenging is calculated following Parekh et al. (2004). The total concentration of dissolved iron (FeT) is separated into free iron (Fe) and iron complexed with organic ligands (FeL), which is not scavenged. Complexation reactions are fast (Tagliabue and Völker2011), so we assume instantaneous equilibrium between free iron and free ligand (L), which is computed using a constant KFeL=[Fe][L][FeL], by solving the following:

(A5) Fe T = Fe + Fe L L T = Fe L + L .

For simplicity, we assume here a constant total ligand concentration LT, 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 KFeL and LT 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 CaCO3 in zooplankton guts. Loss terms are carbon fixation by primary producers and the formation of calcium carbonate. In addition, the sea–air flux of CO2 leads to an exchange of carbon with the atmosphere, depending on the partial pressure difference in the CO2 between ocean and atmosphere. This exchange is treated separately as a boundary condition. The partial pressure of surface ocean CO2 is computed using the mocsy 2.0 routines (Orr and Epitalon2015).

(A6) SMS ( DIC ) = r small - P small PhyC small Net respiration, small phytoplankton + r dia - P dia PhyC dia Net respiration, diatom + ρ DOC f T DOC Remineralization of DOC + r zoo ZooC Respiration, small zoo. + r zoo 2 Zoo 2 C Respiration, macrozoo. + Diss calc DetCalc Calcite dissolution, slow-sinking detritus + G small zoo q small CaCO 3 : N Diss calc_guts CaCO 3 dissolution in guts, small zoo. - ψ P small PhyC small Calcification + Diss calc 2 DetZ 2 Calc Calcite dissolution, fast-sinking detritus + G small zoo 2 q small CaCO 3 : N Diss calc_guts CaCO 3 dissolution in guts, macrozoo.

The state variables PhyCsmall, PhyCdia, DOC, ZooC, Zoo2C, DetCalc, and DetZ2Calc are listed in Table A1. Respiration rate constants of small phytoplankton (rsmall), diatoms (rdia), and zooplankton groups (rzoo and rzoo2) are computed in Sects. A3.2 and A4.1, respectively. Photosynthesis terms (Psmall and Pdia) are calculated in Eq. (A36). The remineralization rate constant (ρDOC) is listed in Table A8, and the temperature dependency (fT) is given in Eq. (A33). Calcite dissolution by detritus (Disscalc; Disscalc2) is calculated in Eq. (A28). The constant for the dissolution of calcium carbonate in zooplankton guts (Disscalc_guts) is listed in Table A5. Gsmallzoo and Gsmallzoo2 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 CaCO3 in zooplankton guts. Alkalinity increases when nitrogen is assimilated and when CaCO3 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).

(A7) SMS ( Alk ) = ( 1 + 1 / 16 ) V small N PhyC small N assimilation, small phytoplankton + ( 1 + 1 / 16 ) V dia N PhyC dia N assimilation, diatom - ( 1 + 1 / 16 ) ρ DON f T DON Remineralization of DON - 2 ψ P small PhyC small Calcification + 2 Diss calc DetCalc Calcite dissolution, slow-sinking detritus + 2 G small zoo q small CaCO 3 : N Diss calc_guts CaCO 3 dissolution in guts, small zoo. + 2 Diss calc 2 DetZ 2 Calc Calcite dissolution, fast-sinking detritus + 2 G small zoo 2 q small CaCO 3 : N Diss calc_guts CaCO 3 dissolution in guts, macrozoo.

The state variables PhyCsmall, PhyCdia, DON, DetCalc, and DetZ2Calc are listed in Table A1. The N assimilation (VsmallN and VdiaN) is calculated in Sect. A3.4. The remineralization rate constant (ρDON) is given in Table A8. The temperature dependency (fT) is calculated in Eq. (A33). The value of the calcite production ratio (ψ) is given in Table A3. The photosynthesis term (Psmall) is calculated in Eq. (A36). The calcite dissolution by detritus (Disscalc, Disscalc2) is calculated in Eq. (A28). Dissolution of calcium carbonate in guts (Disscalc_guts) is listed in Table A5. Gsmallzoo and Gsmallzoo2 are grazing terms and explained in Sect. A4.2.

A1.3 Phytoplankton


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.

(A8)SMS(PhyNsmall)=VsmallNPhyCsmallN assimilation-ϵphyNflim, smallN:CmaxPhyNsmallDON excretion-AggPhyNsmallAggregation loss-GsmallzooGrazing loss by small zoo.-Gsmallzoo2Grazing loss by macrozoo.(A9)SMS(PhyNdia)=VdiaNPhyCdiaN assimilation-ϵphyNflim, diaN:CmaxPhyNdiaDON excretion-AggPhyNdiaAggregation loss-GdiazooGrazing loss by small zoo.-Gdiazoo2Grazing loss by macrozoo.

The state variables PhyCsmall, PhyNsmall, PhyCdia, and PhyNdia are listed in Table A1. The N assimilation (VsmallN and VdiaN) is explained in Sect. A3.4. The constant excretion rate constant (ϵphyN) 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 (flim, smallN:Cmax; flim, diaN: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 (Gsmallzoo, Gsmallzoo2, Gdiazoo, and Gdiazoo2) are explained in Sect. A4.2.


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.

(A10)SMS(PhyCsmall)=(Psmall-rsmall)PhyCsmallNet photosynthesis-AggPhyCsmallAggregation loss-ϵphyCflim, smallN:CmaxPhyCsmallExcretion of DOC-qsmallC:NGsmallzooGrazing loss by small zoo.-qsmallC:NGsmallzoo2Grazing loss by macrozoo.(A11)SMS(PhyCdia)=(Pdia-rdia)PhyCdiaNet photosynthesis-AggPhyCdiaAggregation loss-ϵphyCflim, diaN:CmaxPhyCdiaExcretion of DOC-qdiaC:NGdiazooGrazing loss by small zoo.-qdiaC:NGdiazoo2Grazing loss by macrozoo.

The state variables PhyCsmall and PhyCdia are listed in Table A1. The photosynthesis terms (Psmall and Pdia) are calculated in Eq. (A36). The rates of respiration by small phytoplankton (rsmall) and diatoms (rdia) are explained in Sect. A3.2. The constant for DOC excretion rate of phytoplankton (ϵphyC; Table A8) is downregulated by the limiter factor (flim, smallN:Cmax; flim, diaN:Cmax) when the N:C ratio becomes too high (Eq. A45). Phytoplankton aggregation (Agg) is calculated in Eq. (A42). Grazing terms (Gsmallzoo, Gsmallzoo2, Gdiazoo, and Gdiazoo2) are explained in Sect. A4.2. qC:N=PhyC/PhyN, is used to convert the grazing units from millimoles of N to millimoles of C.


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 CaCO3 by heterotrophs, such as foraminifera or pteropods is neglected. Biogenic CaCO3 produced by coccolithophores is transformed into detritus CaCO3 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.

(A12) SMS ( PhyCalc ) = ψ P small PhyC small Calcification - r small PhyCalc Respiration - G small zoo q small CaCO 3 : N Grazing loss, small zoo. - G small zoo 2 q small CaCO 3 : N Grazing loss, macrozoo. - ϵ phy C f lim, small N : Cmax PhyCalc Excretion loss - Agg PhyCalc Aggregation loss

The state variables PhyCsmall and PhyCalc are listed in Table A1. The value of the calcite production ratio (ψ) is given in Table A3. The constant excretion rate (ϵphyC; Table A8) is downregulated by the limiter factor flim, smallN:Cmax (Eq. A45) when the N:C ratio becomes too high. Photosynthesis (Psmall), respiration (rsmall), and the aggregation of phytoplankton (Agg) rates are calculated in Eqs. (A36), (A38), and (A42), respectively. Grazing terms (Gsmallzoo and Gsmallzoo2) are explained in Sect. A4.2. qsmallCaCO3:N=PhyCalc/PhyNsmall is used to convert the grazing units from millimoles of N to millimoles of CaCO3.

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.

(A13) SMS ( PhySi ) = V Si PhyC dia Diatom Si assimilation - ϵ phy N f lim, dia N : Cmax PhySi dia Excretion to detritus - Agg PhySi dia Aggregation loss - G dia zoo q Si : N Grazing loss, small zoo. - G dia zoo 2 q Si : N Grazing loss, macrozoo.

The state variables PhyCdia and PhySidia are described in Table A1. Si assimilation (VSi) and aggregation rates (Agg) are calculated in Eqs. (A41) and (A42), respectively. The constant excretion rate (ϵphyN; Table A8) is downregulated by the limiter factor flim, diaN:Cmax (Eq. A45) when the N:C ratio becomes too high. Grazing terms (Gdiazoo and Gdiazoo2) are explained in Sect. A4.2. The intracellular ratio between diatom silicon and nitrate is defined as qSi:N=PhySidia/PhyNdia.

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.

(A14)SMS(PhyChlsmall)=SsmallchlPhyCsmallChlorophyllasynthesis-GsmallzooqsmallChl:NGrazing loss, small zoo.-Gsmallzoo2qsmallChl:NGrazing loss, macrozoo.-degsmallchlPhyChlsmallDegradation loss-AggPhyChlsmallAggregation loss(A15)SMS(PhyChldia)=SdiachlPhyCdiaChlorophyllasynthesis-GdiazooqdiaChl:NGrazing loss, small zoo.-Gdiazoo2qdiaChl:NGrazing loss, macrozoo.-degdiachlPhyChldiaDegradation loss-AggPhyChldiaAggregation loss

The state variables PhyCsmall, PhyCdia, PhyChlsmall, and PhyChldia are listed in Table A1. The chlorophyll a synthesis (Ssmallchl; Sdiachl) and the aggregation (Agg) terms are calculated in Eqs. (A39) and (A42), respectively. The degradation parameters (degsmallchl; degdiachl) are given in Table A8. Grazing terms (Gsmallzoo, Gsmallzoo2, Gdiazoo, and Gdiazoo2) are explained in Sect. A4.2. The conversion factor from millimoles of N to milligrams of Chl a is defined as qChl:N=PhyChl/PhyN.

A1.4 Zooplankton


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.

(A16)SMS(ZooN)=γzooGtotzooGrazing-GzooGrazing loss, macrozoo.-mzooZooN2Mortality-ϵzooNZooNExcretion of DON(A17)SMS(Zoo2N)=γzoo2Gtotzoo2Grazing-mzoo2Zoo2N2Mortality-ϵzoo2NZoo2NExcretion of DON-fnGtotzoo2Fecal pellet

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 (Gtotzoo and Gtotzoo2) are calculated in Sect. A4.2. The mortality parameter (mzoo and mzoo2) and fecal pellet production rate constant (fn) are listed in Table A3. The DON excretion terms (ϵzooN and ϵzoo2N) are given in Table A8.


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.

(A18)SMS(ZooC)=γzoo(GsmallzooqsmallC:N+GdiazooqdiaC:N)Grazing on phytoplankton+γzoo(GdetzooqdetC:N+GdetZ2zooqdetZ2C:N)Grazing on detritus-GzooqzooC:NGrazing loss by macrozoo.-mzooZooN2qzooC:NZooplankton mortality-rzooZooCRespiration loss-ϵzooCZooCExcretion of DOC(A19)SMS(Zoo2C)=γzoo2(Gsmallzoo2qsmallC:N+Gdiazoo2qdiaC:N)Grazing on phytoplankton+γzoo2(Gdetzoo2qdetC:N+GdetZ2zoo2qdetZ2C:N)Grazing on detritus+γzoo2(GzooqzooC:N)Grazing on small zoo.-mzoo2Zoo2N2qzoo2C:NZooplankton mortality-rzoo2Zoo2CRespiration loss-ϵzoo2CZoo2CExcretion of DOC-fcGcfluxFecal pellet

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 (Gsmallzoo, Gdiazoo, Gdetzoo, GdetZ2zoo, Gsmallzoo2, Gdiazoo2, Gdetzoo2, GdetZ2zoo2, and Gzoo) are calculated in Sect. A4.2. The respiration terms of zooplankton (rzoo and rzoo2) are calculated in Eqs. (A50) and (A51). Mortality parameters (mzoo and mzoo2) are listed in Table A3. The DOC excretion terms (ϵzooC,ϵzoo2C) are in Table A8. The grazing flux in terms of nitrogen biomass is converted to carbon biomass using the respective intracellular C:N ratios (qsmallC:N, qdiaC:N, qdetC:N, qdetZ2C:N, qzooC:N, and qzoo2C:N), where qsmallC:N=PhyCsmall/PhyNsmall, qdiaC:N=PhyCdia/PhyNdia, qdetC:N=DetC/DetN, qdetZ2C:N=DetZ2C/DetZ2N, qzooC:N=ZooC/ZooN, and qzoo2C:N=Zoo2C/Zoo2N. Total grazed carbon biomass (Gcflux) and the fecal pellet production rate constant (fc; Table A3) together determine the fraction of carbon being lost to the large detritus carbon pool via fecal pellets.

A1.5 Detritus


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.

(A20)SMS(DetN)=(Gsmallzoo+Gdiazoo)(1-γzoo)Sloppy feeding+mzooZooN2Zooplankton mortality-γzoo(Gdetzoo+GdetZ2zoo)Grazing loss, small zoo.+Agg(PhyNsmall+PhyNdia)Phytoplankton aggregation-ρDetNfTDetNDegradation to DON(A21)SMS(DetZ2N)=(Gsmallzoo2+Gdiazoo2+Gzoo)(1-γzoo2)Sloppy feeding-γzoo2(Gdetzoo2+GdetZ2zoo2)Grazing loss, macrozoo.+mzoo2Zoo2N2Mortality+fnGtotFecal pellet-ρDetNfTDetZ2NDegradation to DON

The state variables PhyNsmall, PhyNdia, ZooN, DetN, Zoo2N, and DetZ2N are listed in Table A1. The grazing efficiency (γzoo and γzoo2), mortality (mzoo and mzoo2), and fecal pellet production rate constant (fn) are listed in Table A3. Grazing terms (Gsmallzoo, Gdiazoo, Gdetzoo, GdetZ2zoo, Gsmallzoo2, Gdiazoo2, Gdetzoo2, GdetZ2zoo2, and Gzoo) are calculated in Sect. A4.2. The remineralization rate constant of DON (ρDetN) is listed in Table A8. The temperature dependency fT is calculated in Eq. (A33). The aggregation (Agg) term is calculated in Eq. (A42).


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.

(A22)SMS(DetC)=(GsmallzooqsmallC:N+GdiazooqdiaC:N)(1-γzoo)Sloppy feeding+mzooZooN2qzooC:Nsmall zoo. mortality-γzoo(GdetzooqdetC:N+GdetZ2zooqdetZ2C:N)Grazing loss by macrozoo.+Agg(PhyCsmall+PhyCdia)Phytoplankton aggregation-ρDetCfTDetCDegradation to DOC(A23)SMS(DetZ2C)=(Gsmallzoo2qsmallC:N+Gdiazoo2qdiaC:N+GzooqzooC:N)(1-γzoo2)Sloppy feeding-γzoo2(Gdetzoo2qdetC:N+GdetZ2zoo2qdetZ2C:N)Grazing loss by macrozoo.+mzoo2Zoo2N2qzoo2C:NMortality+fcGcfluxFecal pellet-ρDetCfTDetZ2CDegradation to DOC

The state variables PhyCsmall, PhyCdia, ZooN, DetC, Zoo2N, and DetZ2C are listed in Table A1. The grazing efficiency (γzoo and γzoo2) and mortality (mzoo and mzoo2) parameters are listed in Table A3. Grazing terms (Gsmallzoo, Gdiazoo, Gdetzoo, GdetZ2zoo, Gsmallzoo2, Gdiazoo2, Gdetzoo2, GdetZ2zoo2, and Gzoo) are calculated in Sect. A4.2. The remineralization rate of DOC (ρDetC) is listed in Table A8. Temperature dependency fT is calculated in Eq. (A33). The aggregation (Agg) term is calculated in Eq. (A42). The total grazed carbon biomass (Gcflux) and the fecal pellet production rate constant (fc; Table A3) together determine the fraction of carbon being lost to the large detritus carbon pool via fecal pellets. The quotas qsmallC:N=PhyCsmall/PhyNsmall, qdiaC:N=PhyCdia/PhyNdia, qzooC:N=ZooC/ZooN, qzoo2C:N=Zoo2C/Zoo2N, qdetC:N=DetC/DetN, and qdetZ2C:N=DetZ2C/DetZ2N are used to convert the units from mmol N to mmol C.


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

(A24)SMS(DetSi)=(ϵphyNflim, diaN:CmaxDiatom excretion+Agg)AggregationDiaSi+GdiazooqSi:NSloppy feeding-ρSiTDetSiRemineralization to DSi(A25)SMS(DetZ2Si)=Gdiazoo2qSi:NSloppy feeding-ρSiTDetZ2SiRemineralization to DSi

The state variables DiaSi, DetSi, and DetZ2Si are listed in Table A1. The constant excretion rate (ϵphyN; Table A8) is downregulated by the limiter factor flim, diaN:Cmax (Eq. A45) when the N:C ratio becomes too high. The remineralization rates (ρSiT), the aggregation (Agg), and the grazing on diatoms (Gdiazoo and Gdiazoo2) are calculated in Eqs. (A35), (A42), and (A55), respectively. The intracellular ratio between diatom silicon and carbon is defined as qSi:N=PhySidia/PhyNdia.


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

(A26)SMS(DetCalc)=ϵCphyflim, smallN:CmaxPhyCalcSmall phytoplankton, excretion+(AggAggregation+rsmallRespiration)PhyCalc+GsmallzooqsmallCaCO3:NGrazing loss-GsmallzooqsmallCaCO3:NDisscalc-gutsCaCO3dissolution in guts-DisscalcDetCalcCaCO3dissolution, slow-sinking detritus(A27)SMS(DetZ2Calc)=Gsmallzoo2qsmallCaCO3:NGrazing loss-Gsmallzoo2qsmallCaCO3:NDisscalc-gutsCaCO3dissolution in guts-Disscalc2DetZ2CalcCaCO3dissolution, fast-sinking detritus

The state variables PhyCalc, DetCalc, and DetZ2Calc are listed in Table A1. The constant excretion rate (ϵphyC; Table A8) is downregulated by the limiter factor flim, smallN:Cmax (Eq. A45) when the N:C ratio becomes too high. The respiration (rsmall), the aggregation (Agg), and the grazing on small phytoplankton (Gsmallzoo and Gsmallzoo2) are calculated in Eqs. (A38), (A42), and (A54), respectively. The ratio qsmallCaCO3:N=PhyCalc/PhyNsmall.

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

(A28) Diss calc = Diss calc_rate w det Diss calc 2 = Diss calc_rate

Disscalc and Disscalc2 are the dissolution rate constants for slow- and fast-sinking detritus classes (Table A5). The reference dissolution rate (Disscalc_rate; Table A8) is based on a length scale of 3500 m and velocity of 20 m d−1. The sinking speed at depth z (wdet; Table A5) is calculated as follows:

(A29) w det = 0.0288 z + w 0 .

Here, z denotes the depth, and w0 is the sinking speed at the ocean surface (Table A3). The dissolution rate for a fast-sinking detritus class (Disscalc2) is assumed to be constant throughout the water column and is set to the value of Disscalc_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 O2 leads to an exchange of oxygen with the atmosphere, depending on the partial pressure difference in the O2 between ocean and atmosphere. This exchange is treated separately as a boundary condition. The partial pressure of surface ocean O2 is computed using the mocsy 2.0 routines (Orr and Epitalon2015).

(A30) SMS ( Oxy ) = P small - r small PhyC small Net production, small phytoplankton + P dia - r dia PhyC dia Net production, diatom - ρ DOC f T DOC Remineralization of DOC - r zoo ZooC Respiration, small zoo. - r zoo 2 Zoo 2 C Respiration, macrozoo.

The state variables PhyCsmall, PhyCdia, DOC, ZooC, and Zoo2C are listed in Table A1. Respiration rate constants of small phytoplankton (rsmall), diatoms (rdia), and zooplankton groups (rzoo and rzoo2) are computed in Sects. A3.2 and A4.1, respectively. Photosynthesis terms (Psmall and Pdia) are calculated in Eq. (A36). The remineralization rate constant ρDOC is listed in Table A8 and the temperature dependency (fT) 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.

(A31) SMS ( DON ) = ϵ phy N f lim, small N : Cmax PhyN small Excretion, small phytoplankton + ϵ dia N f lim, dia N : Cmax PhyN dia Excretion, diatom + ϵ zoo N ZooN Excretion, small zoo. + ϵ zoo 2 N Zoo 2 N Excretion, macrozoo. + ρ DetN f T DetN Detritus degradation, slow-sinking + ρ Det 2 ZN f T DetZ 2 N Detritus degradation, fast sinking - ρ DON f T DON Remineralization

The state variables PhyNsmall, PhyNdia, ZooN, DetN, Zoo2N, DetZ2N, and DON are listed in Table A1. The constant excretion rate of nitrogen from phytoplankton and zooplankton classes (ϵphyN, ϵdiaN, ϵzooN, and ϵzoo2N), 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 (flim, smallN:Cmax and flim, diaN:Cmax; Eq. A45) when the N:C ratio becomes too high. The temperature dependency fT 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.

(A32) SMS ( DOC ) = ϵ phy C f lim, small N : Cmax PhyC small Excretion, small phytoplankton + ϵ dia C f lim, dia N : Cmax PhyC dia Excretion, diatom + ϵ zoo C ZooC Excretion, small zoo. + ϵ zoo 2 C Zoo 2 C Excretion, macrozoo. + ρ DetC f T DetC Detritus degradation, slow sinking + ρ Det 2 C f T Det 2 C Detritus degradation, fast sinking - ρ DOC f T DOC Remineralization

The state variables PhyCsmall, PhyCdia, ZooC, DetC, Zoo2C, Det2C, and DOC are listed in Table A1. The constant excretion rate of nitrogen from phytoplankton and zooplankton classes (ϵphyC, ϵdiaC, ϵzooC, and ϵzoo2C), 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 (flim, smallN:Cmax and flim, diaN:Cmax; Eq. A45) when the N:C ratio becomes too high. Temperature dependency fT 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.

(A33) f T = exp - 4500 1 T - 1 T ref

T and Tref 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örtner2016) is used for the parameterization of the temperature dependency (fTzoo2; Table A5). Specifically, the following parameterization provides an optimum curve with a maximum at 0.5 C, as described in Karakuş et al. (2021).

(A34) f Tzoo 2 = exp Q a T r - Q a T 1 + exp Q h T h - Q h T

Tr is the intrinsic optimum temperature for development, and Th is the temperature above which inhibitive processes dominate. Qa and Qh 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 (ρSiT; Table A5) is calculated following Maerz et al. (2020) but with a minimum dissolution rate.

(A35) ρ Si T = max 0.023 2.6 T - 10 10 , ρ Si

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 (Psmall and Pdia) is parameterized as follows:

(A36) P small = P max small 1.0 - exp - α small q Chl : C PAR P max small , P dia = P max dia 1.0 - exp - α dia q Chl : C PAR P max dia .

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 (qChl:C) is defined as PhyChl/PhyC and varies as a result of photoacclimation. The apparent maximum photosynthetic rate (Pmaxsmall and Pmaxdia) is defined below.

(A37) P max small = μ C small max min f lim, small Fe , f lim, small N : Cmin f T , P max dia = μ C, dia max min f lim, dia Fe , f lim, dia N : Cmin , f lim, dia Si : Cmin f T

The value of μC smallmax and μC, diamax is listed in Table A7. The limitation terms (flim, smallN:Cmin, flim, diaN:Cmin, flim, diaSi:Cmin, flim, smallFe, and flim, diaFe) are presented in Sect. A3.6, and the temperature dependency (fT) is calculated in Eq. (A33).

A3.2 Respiration

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

(A38) r small = res small f lim, small N : Cmax Maintenance + ζ V small N N assim , r dia = res dia f lim, dia N : Cmax Maintenance + ζ V dia N N assim .

The values for the maintenance respiration rate (ressmall and resdia) 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 (Hohn2009). The limiter function (flim, smallN:Cmax and flim, diaN:Cmax) is described in Eq. (A45), and the N assimilation rate (VsmallN and VdiaN) is calculated in Eq. (A40).

A3.3 Chlorophyll a synthesis

The chlorophyll synthesis rate (Ssmallchl and Sdiachl; 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.

(A39) S small chl = V small N q max, small Chl : N min 1 , P small α small q Chl : C PAR , S dia chl = V dia N q max, dia Chl : N min 1 , P dia α dia q Chl : C PAR

The N assimilation (VsmallN and VdiaN) is computed in Eq. (A40). The conversion factor of the maximum Chl:N ratio (qmax, smallChl:N and qmax, diaChl:N) and the light-harvesting efficiency values (αsmall and αdia) are listed in Table A7. The C-specific photosynthesis (Psmall and Pdia) is given in Eq. (A36). PAR is the photosynthetically available radiation (Table A5), and the intracellular Chl to C ratio (qChl:C) is defined as PhyChl/PhyC.

A3.4 N and Si assimilation


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 (Pmax) 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 flimN:Cmax (see Eq. A45). Nitrogen assimilation is downregulated at high intracellular N:C ratio.

(A40) V small N = V cm small P max small σ N : C small f lim, small N : Cmax DIN K small N + DIN , V dia N = V cm dia P max dia σ N : C dia f lim, dia N : Cmax DIN K dia N + DIN

Vcmsmall, Vcmdia, σN:Csmall, σN:Cdia, KsmallN, and KdiaN are listed in Table A7. The maximum rate of photosynthesis (Pmaxsmall and Pmaxdia) is given in Eqs. (A37). flim, smallN:Cmax and flim, diaN:Cmax are described in Eq. (A45). DIN corresponds to the in situ concentration.


The building of a silica frustule of diatoms requires silicate uptake. The C-specific Si assimilation rate (VSi) is a function of a factor for C-specific N uptake (Vcmdia), a rate constant of C-specific photosynthesis (μC,diamax), 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 flim) further regulate the Si assimilation.

(A41) V Si = V cm dia μ C, dia max f T σ Si : C f lim Si : Cmax f lim, dia N : Cmax DSi K Si + DSi

The scaling factor for the N uptake (Vcmdia), the maximum rate constant of C-specific photosynthesis (μC, diamax), the uptake ratio of the maximum Si:C (σSi:C), and half-saturation constant for silicate uptake (KSi) are listed in Table A7. The temperature dependency (fT) is computed in Eq. (A33). The limitation by the intracellular ratios N:C and Si:C (flim, diaN:Cmax and flimSi: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 (1-qlimdia). When the nutrient limitation is high (i.e. low qlimdia), the aggregation rate increases in the model.

(A42)Agg=ϕphyPhyNsmall+(1-qlimdia)PhyNdia+ϕdetDetN+DetZ2N(A43)qlimdia=minflim, diaFe,flim, diaN:Cmin,flim, diaSi:Cmin

The state variables PhyNsmall, PhyNdia, 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 (flim, diaN:Cmin, flim, diaSi:Cmin, and flim, diaFe) 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).

(A44) f lim ( θ , q 1 , q 2 ) = 1 - exp ( - θ ( | Δ q | - Δ q ) 2 )

Here, Δq=q1-q2 is the difference between the current intracellular nutrient:C quota and a prescribed maximum or minimum quota. The dimensionless constant θ controls the limitation.


The limiter flimN:Cmax downregulates the metabolic processes such as nitrogen and Si assimilation, excretion, and maintenance respiration of phytoplankton when the intracellular nitrogen quota (qN:C) becomes too high. flimN:Cmax is one when the current qN:C<0.151 (i.e. Redfield ratio; 16N:106C) and zero for qN:C>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 CaCO3 processes of phytoplankton (see the discussion of CaCO3 in Sect. A1.5).

(A45) f lim, small N : Cmax = f lim ( θ max N , q small N : C , q small N : Cmax ) , f lim, dia N : Cmax = f lim ( θ max N , q dia N : C , q dia N : Cmax )

The limitation function for the quota regulation is calculated with Eq. (A44). qsmallN:C and qdiaN:C form the current intracellular nitrogen quota for small phytoplankton and diatoms, respectively. Dimensionless constants θmaxN, qsmallN:Cmax, and qdiaN:Cmax are listed in Table A6.


The limiter flimSi:Cmax downregulates the Si assimilation of diatoms when the intracellular silicon quota (Si:C) becomes too high. flimSi:Cmax is one when the current qN:C<0.76 and zero for qN:C>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:

(A46) f lim Si : Cmax = f lim ( θ max Si , q Si : C , q Si : Cmax )

Dimensionless constants θmaxSi and qSi:Cmax are listed in Table A6.


Carbon fixation and aggregation loss in diatoms are further downregulated by a factor (flim, diaSi:Cmin; see Eq. A44) when the intracellular silicon quota (qSi:C) approaches a minimum value (qSi:Cmin), mimicking the arrest of cellular division at low cellular Si (Claquin et al.2002). flim, diaSi:Cmin is zero when the current qSi:C<0.04 and one for qSi:C>0.08.

(A47) f lim, dia Si : Cmin = f lim ( θ min Si , q Si : Cmin , q Si : C )

The dimensionless constants θminSi and qSi:Cmin are listed in Table A6.


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

(A48) f lim, small Fe = DFe K small Fe + DFe , f lim, dia Fe = DFe K dia Fe + DFe

The variable DFe is listed in Table A1. The half-saturation constants (KsmallFe and KdiaFe) are given in Table A6.


In addition to iron limitation, photosynthesis is limited by nitrogen in small phytoplankton and diatoms using the Eq. (A44). Nitrogen limitation (flim, smallN:Cmin and flim, diaN:Cmin) is described as a function of the intracellular nitrogen quota (qsmallN:C and qdiaN:C) with growth ending at a minimum quota (qsmallN:Cmin and qdiaN:Cmin).

(A49) f lim, small N : Cmin = f lim ( θ min N , q small N : Cmin , q small N : C ) , f lim, dia N : Cmin = f lim ( θ min N , q dia N : Cmin , q dia N : C )

Dimensionless constants θminN, qsmallN:Cmin, and qdiaN: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 (rzoo; Table A5) is assumed to drive it back with a timescale τ.

(A50) r zoo = q zoo C : N - q standard C : N τ f T

The timescale for respiration (τ) is listed in Table A7. The temperature dependence (fT) is calculated in Eq. (A33). The ratios are defined as qzooC:N=ZooC/ZooN and qStandardC:N=106/16.


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

(A51) r zoo 2 = R s ( 1 + R f + R a )

The standard respiration rate (Rs) is listed in Table A3. The feeding activity factor (Rf; 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 (Ra; 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:

(A52) G tot zoo = ξ zoo i p i N i 2 σ zoo + i p i N i 2 f T ZooN .

Gtotzoo (Gtotzoo2) 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 (fT and fTzoo2) are given in Eqs. (A33) and (A34). In the model, relative grazing preferences are implemented following Fasham et al. (1990). Variable relative grazing preferences (pi) are calculated using the nominal preferences for small phytoplankton, diatoms, slow- and fast-sinking detritus, and small zooplankton (Table A10) as follows:

(A53) p i = p i N i i p i N i .

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; PhyNsmall), diatoms (i=2; PhyNdia), 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 (Gsmall, Gdia, Gdet, GdetZ2, and Gzoo).


where Gzoo is associated with macrozooplankton grazing on small zooplankton. PhyNsmall, PhyNdia, 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 CaCO3 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 wdet (Eq. A29) and wdetZ2 = 200 m d−1, respectively. This results in a detrital flux (BenFDetN, BenFDetZ2N, BenFDetC, BenFDetZ2C, BenFDetSi, BenFDetZ2Si, BenFDetCalc, and BenFDetZ2Calc; Table A11) from the water column to the benthos.


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.


BenFDIN, BenFDSi, BenFDIC, and BenFAlk (Table A11) denote the fluxes of DIN, DSi, DIC, and Alk returned into the bottom layer of the ocean. Constant remineralization rates (ρbenN, ρbenSi, and ρbenC) are listed in Table A8. The calcite dissolution rates Disscalc and Disscalc2 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 A1List of oceanic state variables in REcoM3.

Download Print Version | Download XLSX

Table A2List of benthic state variables in REcoM3.

Download Print Version | Download XLSX

Table A3Parameters for equations of sources minus sinks.

Download Print Version | Download XLSX

Table A4Parameters for iron calculations.

Download Print Version | Download XLSX

Table A5Model variables.

Download Print Version | Download XLSX

Table A6Parameters for limitation functions.

Download Print Version | Download XLSX

Table A7Parameters for phytoplankton processes.

Download Print Version | Download XLSX

Table A8Degradation parameters for sources-minus-sinks equations.

Download Print Version | Download XLSX

Table A9Parameters for macrozooplankton grazing.

Download Print Version | Download XLSX

Table A10Parameters for zooplankton processes.

Download Print Version | Download XLSX

Table A11Benthos variables.

Download Print Version | Download XLSX

Code availability

The FESOM2.1–REcoM3 source code is available at (last access: 31 December 2022). The version of FESOM2.1–REcoM3 used for this paper can be found at (Gürses2023). A manual is available at (last access: 23 July 2023).

Data availability

The GLODAPv2 mapped product of dissolved inorganic carbon and alkalinity can be found at (last access: 26 July 2023) (Lauvset et al.2016). Dissolved inorganic nitrogen and dissolved inorganic silicon from the World Ocean Atlas 2013 can be found at (Garcia et al.2013). The oxygen data set from the World Ocean Atlas 2018 is available at (last access: 26 July 2023) (Garcia et al.2019a). The OC-CCI data set is available at (last access 26 July 2023) (Sathyendranath et al.2019). The Southern Ocean chlorophyll a concentration data set is available at (last access 26 July 2023) (Johnson et al.2013). CbPM (Westberry et al.2008) and VGPM (Behrenfeld and Falkowski1997) products of net primary production are available at (last access 26 July 2023). MAREDAT products can be found at (Buitenhuis et al.2012), (O'Brien and Moriarty2012), and (Moriarty2012). Global maps of trends and climatological fields are available at (Sallée et al.2021b). The Surface Ocean CO2 Atlas (SOCAT) was taken from (last access: 26 July 2023) (Bakker et al.2016). The Global Ocean Surface Carbon is available at (last access: 26 July 2023) (Chau et al.2022). Annual air–sea CO2 flux data is freely available at (last access 26 July 2023) (Friedlingstein et al.2022a). The Polar science center Hydrographic Climatology (Steele et al.2001) data used for model initialization and the CORE-II freshwater runoff data (Large and Yeager2004) are available online. The JRA-55-based surface data set for driving ocean–sea ice models (JRA55-do) is taken from (last access: 23 August 2023, Japan Meteorological Agency2023). The model results are available at (Gürses et al.2023).

Author contributions

The conceptualization was done by ÖG, LO, CV, and JH. The data were prepared by ÖG, JH, and OK. All authors analysed the simulations and visualizations and contributed to the writing of the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


The work reflects only the authors' views; the European Commission and their executive agency are not responsible for any use that may be made of the information that the work contains.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank Sergey Danilov from AWI Bremerhaven, Germany, for his helpful support. We acknowledge the Global Carbon Project, which is responsible for the Global Carbon Budget, and we thank the ocean modelling and fCO2 mapping groups for producing and making available their model and fCO2 product output.

Financial support

This research has been funded by the Initiative and Networking Fund of the Helmholtz Association (Helmholtz Young Investigator Group Marine Carbon and Ecosystem Feedbacks in the Earth System, MarESys; grant no. VH-NG-1301) and the ERC-2022-STG OceanPeak (grant no. 101077209). Laurent Oziel and Ying Ye received funding from the European Union's Horizon 2020 research and innovation programme (grant no. 820989; COMFORT project). Laurent Oziel received funding from the Germany Ministry for Education and Research (BMBF; grant no. 03F0918A; nuArctic project). Moritz Zeising has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation; grant no. 268020496 – TRR 172) within the Transregional Collaborative Research Center “ArctiC Amplification: Climate Relevant Atmospheric and SurfaCe Processes, and Feedback Mechanisms (AC)3” project. Ying Ye and Martin Butzin have received funding by the German Federal Ministry of Education and Research (BMBF; grant no. 01LP1919A; PalMod project). Martin Butzin has additionally been funded through the DFG-ANR project of MARCARA.

The article processing charges for this open-access publication were covered by the Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research (AWI).

Review statement

This paper was edited by Riccardo Farneti and reviewed by two anonymous referees.


Albani, S., Mahowald, N. M., Perry, A. T., Scanza, R. A., Zender, C. S., Heavens, N. G., Maggi, V., Kok, J. F., and Otto-Bliesner, B. L.: Improved dust representation in the Community Atmosphere Model, J. Adv. Model. Earth Sy., 6, 541–570,, 2014. a, b, c

Álvarez, E., Thoms, S., and Völker, C.: Chlorophyll to Carbon Ratio Derived From a Global Ecosystem Model With Photodamage, Global Biogeochem. Cy., 32, 799–816,, 2018. a, b

Anderson, L. G., Jutterström, S., Hjalmarsson, S., Wåhlström, I., and Semiletov, I. P.: Out-gassing of CO2 from Siberian Shelf seas by terrestrial organic matter decomposition, Geophys. Res. Lett., 36, L20601,, 2009. a

Anderson, T. R., Gentleman, W. C., and Sinha, B.: Influence of grazing formulations on the emergent properties of a complex ecosystem model in a global ocean general circulation model, Prog. Oceanogr., 87, 201–213,, 2010. a

Aumont, O., Orr, J. C., Monfray, P., Ludwig, W., Amiotte-Suchet, P., and Probst, J.-L.: Riverine-driven interhemispheric transport of carbon, Global Biogeochem. Cy., 15, 393–405,, 2001. a, b

Aumont, O., Maier-Reimer, E., Blain, S., and Monfray, P.: An ecosystem model of the global ocean including Fe, Si, P colimitations, Global Biogeochem. Cy., 17, 1060,, 2003. a

Aumont, O., Ethé, C., Tagliabue, A., Bopp, L., and Gehlen, M.: PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies, Geosci. Model Dev., 8, 2465–2513,, 2015. a, b, c

Bakker, D., De Baar, H., and Bathmann, U.: Changes of carbon dioxide in surface waters during spring in the Southern Ocean, Deep-Sea Res. Pt. II, 44, 91–127,, 1997. a

Bakker, D. C. E., Pfeil, B., Landa, C. S., Metzl, N., O'Brien, K. M., Olsen, A., Smith, K., Cosca, C., Harasawa, S., Jones, S. D., Nakaoka, S., Nojiri, Y., Schuster, U., Steinhoff, T., Sweeney, C., Takahashi, T., Tilbrook, B., Wada, C., Wanninkhof, R., Alin, S. R., Balestrini, C. F., Barbero, L., Bates, N. R., Bianchi, A. A., Bonou, F., Boutin, J., Bozec, Y., Burger, E. F., Cai, W.-J., Castle, R. D., Chen, L., Chierici, M., Currie, K., Evans, W., Featherstone, C., Feely, R. A., Fransson, A., Goyet, C., Greenwood, N., Gregor, L., Hankin, S., Hardman-Mountford, N. J., Harlay, J., Hauck, J., Hoppema, M., Humphreys, M. P., Hunt, C. W., Huss, B., Ibánhez, J. S. P., Johannessen, T., Keeling, R., Kitidis, V., Körtzinger, A., Kozyr, A., Krasakopoulou, E., Kuwata, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lo Monaco, C., Manke, A., Mathis, J. T., Merlivat, L., Millero, F. J., Monteiro, P. M. S., Munro, D. R., Murata, A., Newberger, T., Omar, A. M., Ono, T., Paterson, K., Pearce, D., Pierrot, D., Robbins, L. L., Saito, S., Salisbury, J., Schlitzer, R., Schneider, B., Schweitzer, R., Sieger, R., Skjelvan, I., Sullivan, K. F., Sutherland, S. C., Sutton, A. J., Tadokoro, K., Telszewski, M., Tuma, M., van Heuven, S. M. A. C., Vandemark, D., Ward, B., Watson, A. J., and Xu, S.: A multi-decade record of high-quality fCO2 data in version 3 of the Surface Ocean CO2 Atlas (SOCAT), Earth Syst. Sci. Data, 8, 383–413,, 2016 (data available at:, last access: 26 July 2023). a, b, c, d

Ballantyne, A. P., Alden, C. B., Miller, J. B., Tans, P. P., and White, J. W. C.: Increase in observed net carbon dioxide uptake by land and oceans during the past 50 years, Nature, 488, 70–72,, 2012. a

Battaglia, G., Steinacher, M., and Joos, F.: A probabilistic assessment of calcium carbonate export and dissolution in the modern ocean, Biogeosciences, 13, 2823–2848,, 2016. a

Behrenfeld, M. J. and Falkowski, P. G.: A consumer’s guide to phytoplankton primary productivity models, Limnol. Oceanogr., 42, 1479–1491,, 1997. a, b, c, d, e, f, g, h

Berelson, W. M., Balch, W. M., Najjar, R., Feely, R. A., Sabine, C., and Lee, K.: Relating estimates of CaCO3 production, export, and dissolution in the water column to measurements of CaCO3 rain into sediment traps and dissolution on the sea floor: A revised global carbonate budget, Global Biogeochem. Cy., 21, GB1024,, 2007. a

Berthet, S., Séférian, R., Bricaud, C., Chevallier, M., Voldoire, A., and Ethé, C.: Evaluation of an Online Grid-Coarsening Algorithm in a Global Eddy-Admitting Ocean Biogeochemical Model, J. Adv. Model. Earth Sy., 11, 1759–1783,, 2019. a, b

Blondeau-Patissier, D., Gower, J. F., Dekker, A. G., Phinn, S. R., and Brando, V. E.: A review of ocean color remote sensing methods and statistical techniques for the detection, mapping and analysis of phytoplankton blooms in coastal and open oceans, Prog. Oceanogr., 123, 123–144,, 2014. a

Bopp, L., Resplandy, L., Orr, J. C., Doney, S. C., Dunne, J. P., Gehlen, M., Halloran, P., Heinze, C., Ilyina, T., Séférian, R., Tjiputra, J., and Vichi, M.: Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models, Biogeosciences, 10, 6225–6245,, 2013. a

Bourgeois, T., Goris, N., Schwinger, J., and Tjiputra, J. F.: Stratification constrains future heat and carbon uptake in the Southern Ocean between 30 S and 55 S, Nat. Commun., 13, 340,, 2022. a

Boyd, P., Claustre, H., Levy, M., Siegel, D. A., and Weber, T.: Multi-faceted particle pumps drive carbon sequestration in the ocean, Nature, 568, 327–335,, 2019. a

Boye, M., van den Berg, C. M., de Jong, J. T., Leach, H., Croot, P., and de Baar, H. J.: Organic complexation of iron in the Southern Ocean, Deep-Sea Res. Pt. I, 48, 1477–1497,, 2001. a

Buitenhuis, E., Rivkin, R. B., Séailley, S., and Le Quéré, C.: Biogeochemical fluxes through microzooplankton, Global Biogeochem. Cy., 24, GB4015,, 2010. a, b, c

Buitenhuis, E. T., Rivkin, R. B., Sailley, S., and Le Quéré, C.: Global distributions of microzooplankton abundance and biomass – Gridded data product (NetCDF) – Contribution to the MAREDAT World Ocean Atlas of Plankton Functional Types, PANGAEA [data set],, 2012. a

Bunsen, F.: Impact of recent climate variability on oceanic CO2 uptake in a global ocean biogeochemistry model, Master's thesis, Kiel University Christian-Albrechts-Universität, Faculty of Mathematics and Natural Sciences, (last access: 14 August 2023​​​​​​​), 2022. a

Bushinsky, S. M., Landschützer, P., Rödenbeck, C., Gray, A. R., Baker, D., Mazloff, M. R., Resplandy, L., Johnson, K. S., and Sarmiento, J. L.: Reassessing Southern Ocean air‐sea CO2 flux estimates with the addition of biogeochemical float observations, Global Biogeochem. Cy., 33, 1370–1388,, 2019. a

Butzin, M. and Pörtner, H. O.: Thermal growth potential of Atlantic cod by the end of the 21st century, Glob. Change Biol., 22, 4162–4168,, 2016. a

Campin, J.-M., Adcroft, A. J., Hill, C., and Marshall, J.: Conservation of properties in a free-surface model, Ocean Model., 6, 221–244,, 2004. a

Canadell, J., Monteiro, P., Costa, M., Cotrim da Cunha, L., Cox, P., Eliseev, A., Henson, S., Ishii, M., Jaccard, S., Koven, C., Lohila, A., Patra, P., Piao, S., Rogelj, J., Syampungani, S., Zaehle, S., and Zickfeld, K.: Global Carbon and other Biogeochemical Cycles and Feedbacks, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 673–816,, 2021. a

Carr, M.-E., Friedrichs, M. A., Schmeltz, M., Noguchi Aita, M., Antoine, D., Arrigo, K. R., Asanuma, I., Aumont, O., Barber, R., Behrenfeld, M., Bidigare, R., Buitenhuis, E. T., Campbell, J., Ciotti, A., Dierssen, H., Dowell, M., Dunne, J., Esaias, W., Gentili, B., Gregg, W., Groom, S., Hoepffner, N., Ishizaka, J., Kameda, T., Le Quéré, C., Lohrenz, S., Marra, J., Mélin, F., Moore, K., Morel, A., Reddy, T. E., Ryan, J., Scardi, M., Smyth, T., Turpie, K., Tilstone, G., Waters, K., and Yamanaka, Y.: A comparison of global estimates of marine primary production from ocean color, Deep-Sea Res. Pt. II, 53, 741–770,, 2006. a

Chau, T. T. T., Gehlen, M., and Chevallier, F.: A seamless ensemble-based reconstruction of surface ocean pCO2 and air–sea CO2 fluxes over the global coastal and open oceans, Biogeosciences, 19, 1087–1109,, 2022 (data available at:, last access: 26 July 2023). a, b, c, d, e, f, g, h, i, j

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., Chhabra, A., DeFries, R., Galloway, J., Heimann, M., Le Quéré, C., Myneni, R., Piao, S., and Thornton, P.: Carbon and Other Biogeochemical Cycles, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 465–570, 2014. a

Claquin, P., Martin-Jézéquel, V., Kromkamp, J. C., Veldhuis, M. J. W., and Kraay, G. W.: Uncoupling of Silicon Compared With Carbon and Nitrogen Metabolisms and the Role of the Cell Cycle in Continuous Cultures of Thalassiosira Pseudonana (Bacillariophyceae) Under Light, Nitrogen, and Phosphorus Control, J. Phycol., 38, 922–930,, 2002. a, b

Cocco, V., Joos, F., Steinacher, M., Frölicher, T. L., Bopp, L., Dunne, J., Gehlen, M., Heinze, C., Orr, J., Oschlies, A., Schneider, B., Segschneider, J., and Tjiputra, J.: Oxygen and indicators of stress for marine life in multi-model global warming projections, Biogeosciences, 10, 1849–1868,, 2013. a

Cram, J. A., Weber, T., Leung, S. W., McDonnell, A. M. P., Liang, J.-H., and Deutsch, C.: The Role of Particle Size, Ballast, Temperature, and Oxygen in the Sinking Flux to the Deep Sea, Global Biogeochem. Cy., 32, 858–876,, 2018. a

Crisp, D., Dolman, H., Tanhua, T., McKinley, G. A., Hauck, J., Bastos, A., Sitch, S., Eggleston, S., and Aich, V.: How Well Do We Understand the Land-Ocean-Atmosphere Carbon Cycle?, Rev. Geophys., 60, e2021RG000736,, 2022. a, b

Danabasoglu, G., Yeager, S. G., Bailey, D., Behrens, E., Bentsen, M., Bi, D., Biastoch, A., Böning, C., Bozec, A., Canuto, V. M., Cassou, C., Chassignet, E., Coward, A. C., Danilov, S., Diansky, N., Drange, H., Farneti, R., Fernandez, E., Fogli, P. G., Forget, G., Fujii, Y., Griffies, S. M., Gusev, A., Heimbach, P., Howard, A., Jung, T., Kelley, M., Large, W. G., Leboissetier, A., Lu, J., Madec, G., Marsland, S. J., Masina, S., Navarra, A., George Nurser, A., Pirani, A., y Mélia, D. S., Samuels, B. L., Scheinert, M., Sidorenko, D., Treguier, A.-M., Tsujino, H., Uotila, P., Valcke, S., Voldoire, A., and Wang, Q.: North Atlantic simulations in Coordinated Ocean-ice Reference Experiments phase II (CORE-II). Part I: Mean states, Ocean Model., 73, 76–107,, 2014. a

Danilov, S.: Ocean modeling on unstructured meshes, Ocean Model., 69, 195–210,, 2013. a

Danilov, S., Wang, Q., Timmermann, R., Iakovlev, N., Sidorenko, D., Kimmritz, M., Jung, T., and Schröter, J.: Finite-Element Sea Ice Model (FESIM), version 2, Geosci. Model Dev., 8, 1747–1761,, 2015. a, b

Danilov, S., Sidorenko, D., Wang, Q., and Jung, T.: The Finite-volumE Sea ice–Ocean Model (FESOM2), Geosci. Model Dev., 10, 765–789,, 2017. a, b, c, d, e, f, g

de Baar, H. J., de Jong, J. T., Nolting, R. F., Timmermans, K. R., van Leeuwe, M. A., Bathmann, U., Rutgers van der Loeff, M., and Sildam, J.: Low dissolved Fe and the absence of diatom blooms in remote Pacific waters of the Southern Ocean, Mar. Chem., 66, 1–34,, 1999. a

Denman, K. L., Brasseur, G., Chidthaisong, A., Ciais, P., Cox, P. M., Dickinson, R. E., Hauglustaine, D., Heinze, C., Holland, E., Jacob, D., Lohmann, U., Ramachandran, S., Leite da Silva Dias, P., Wofsy, S. C., and Zhang, X.: Couplings Between Changes in the Climate System and Biogeochemistry, in: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Solomon, S., Qin, D., Manning, M., Marquis, M., Averyt, K., Tignor, M. M. B., Miller, H. L., and Chen, Z. L., Cambridge University Press, Cambridge, UK and New York, USA, 499–587, 2007. a

DeVries, T., Le Quéré, C., Andrews, O., Berthet, S., Hauck, J., Ilyina, T., Landschützer, P., Lenton, A., Lima, I. D., Nowicki, M., Schwinger, J., and Séférian, R.: Decadal trends in the ocean carbon sink, P. Natl. Acad. Sci. USA, 116, 11646–11651,, 2019. a

Doney, S. C., Lindsay, K., Moore, J. K., Dutkiewicz, S., Friedrichs, M. A. M., and Matear, R. J.: Marine Biogeochemical Modeling: Recent Advances and Future Challenges, Oceanography, 14, 93–107,, 2001. a

Doney, S. C., Lima, I., Feely, R. A., Glover, D. M., Lindsay, K., Mahowald, N., Moore, J. K., and Wanninkhof, R.: Mechanisms governing interannual variability in upper-ocean inorganic carbon system and air–sea CO2 fluxes: Physical climate and atmospheric dust, Deep-Sea Res. Pt. II, 56, 640–655,, 2009. a, b

Du, J., Ye, Y., Zhang, X., Völker, C., and Tian, J.: Southern Control of Interhemispheric Synergy on Glacial Marine Carbon Sequestration, Geophys. Res. Lett., 49, e2022GL099048,, 2022. a

Dunne, J. P., Sarmiento, J. L., and Gnanadesikan, A.: A synthesis of global particle export from the surface ocean and cycling through the ocean interior and on the seafloor, Global Biogeochem. Cy., 21, GB4006,, 2007. a, b, c, d, e, f, g, h

Elrod, V. A., Berelson, W. M., Coale, K. H., and Johnson, K. S.: The flux of iron from continental shelf sediments: a missing source for global budgets, Geophys. Res. Lett., 31, L12307,, 2004. a, b

Fasham, M. J., Ducklow, H. W., and McKelvie, S. M.: A nitrogen-based model of plankton dynamics in the oceanic mixed layer, J. Mar. Res., 48, 591–639, 1990. a, b

Fay, A. R. and McKinley, G. A.: Observed Regional Fluxes to Constrain Modeled Estimates of the Ocean Carbon Sink, Geophys. Res. Lett., 48, e2021GL095325,, 2021. a, b

Fay, A. R., Gregor, L., Landschützer, P., McKinley, G. A., Gruber, N., Gehlen, M., Iida, Y., Laruelle, G. G., Rödenbeck, C., Roobaert, A., and Zeng, J.: SeaFlux: harmonization of air–sea CO2 fluxes from surface pCO2 data products using a standardized approach, Earth Syst. Sci. Data, 13, 4693–4710,, 2021. a

Fennel, K., Mattern, J. P., Doney, S. C., Bopp, L., Moore, A. M., Wang, B., and Yu, L.: Ocean biogeochemical modelling, Nature Reviews Methods Primers, 2, 76,, 2022. a

Friedlingstein, P., Dufresne, J.-L., Cox, P. M., and Rayner, P. J.: How positive is the feedback between climate change and the carbon cycle?, Tellus B, 55, 692–700,, 2003. a

Friedlingstein, P., Jones, M. W., O'Sullivan, M., Andrew, R. M., Bakker, D. C. E., Hauck, J., Le Quéré, C., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S. R., Anthoni, P., Bates, N. R., Becker, M., Bellouin, N., Bopp, L., Chau, T. T. T., Chevallier, F., Chini, L. P., Cronin, M., Currie, K. I., Decharme, B., Djeutchouang, L. M., Dou, X., Evans, W., Feely, R. A., Feng, L., Gasser, T., Gilfillan, D., Gkritzalis, T., Grassi, G., Gregor, L., Gruber, N., Gürses, Ö., Harris, I., Houghton, R. A., Hurtt, G. C., Iida, Y., Ilyina, T., Luijkx, I. T., Jain, A., Jones, S. D., Kato, E., Kennedy, D., Klein Goldewijk, K., Knauer, J., Korsbakken, J. I., Körtzinger, A., Landschützer, P., Lauvset, S. K., Lefèvre, N., Lienert, S., Liu, J., Marland, G., McGuire, P. C., Melton, J. R., Munro, D. R., Nabel, J. E. M. S., Nakaoka, S.-I., Niwa, Y., Ono, T., Pierrot, D., Poulter, B., Rehder, G., Resplandy, L., Robertson, E., Rödenbeck, C., Rosan, T. M., Schwinger, J., Schwingshackl, C., Séférian, R., Sutton, A. J., Sweeney, C., Tanhua, T., Tans, P. P., Tian, H., Tilbrook, B., Tubiello, F., van der Werf, G. R., Vuichard, N., Wada, C., Wanninkhof, R., Watson, A. J., Willis, D., Wiltshire, A. J., Yuan, W., Yue, C., Yue, X., Zaehle, S., and Zeng, J.: Global Carbon Budget 2021, Earth Syst. Sci. Data, 14, 1917–2005,, 2022a (data available at:, last access: 26 July 2023). a, b, c, d, e, f, g, h, i

Friedlingstein, P., O'Sullivan, M., Jones, M. W., Andrew, R. M., Gregor, L., Hauck, J., Le Quéré, C., Luijkx, I. T., Olsen, A., Peters, G. P., Peters, W., Pongratz, J., Schwingshackl, C., Sitch, S., Canadell, J. G., Ciais, P., Jackson, R. B., Alin, S. R., Alkama, R., Arneth, A., Arora, V. K., Bates, N. R., Becker, M., Bellouin, N., Bittig, H. C., Bopp, L., Chevallier, F., Chini, L. P., Cronin, M., Evans, W., Falk, S., Feely, R. A., Gasser, T., Gehlen, M., Gkritzalis, T., Gloege, L., Grassi, G., Gruber, N., Gürses, Ö., Harris, I., Hefner, M., Houghton, R. A., Hurtt, G. C., Iida, Y., Ilyina, T., Jain, A. K., Jersild, A., Kadono, K., Kato, E., Kennedy, D., Klein Goldewijk, K., Knauer, J., Korsbakken, J. I., Landschützer, P., Lefèvre, N., Lindsay, K., Liu, J., Liu, Z., Marland, G., Mayot, N., McGrath, M. J., Metzl, N., Monacci, N. M., Munro, D. R., Nakaoka, S.-I., Niwa, Y., O'Brien, K., Ono, T., Palmer, P. I., Pan, N., Pierrot, D., Pocock, K., Poulter, B., Resplandy, L., Robertson, E., Rödenbeck, C., Rodriguez, C., Rosan, T. M., Schwinger, J., Séférian, R., Shutler, J. D., Skjelvan, I., Steinhoff, T., Sun, Q., Sutton, A. J., Sweeney, C., Takao, S., Tanhua, T., Tans, P. P., Tian, X., Tian, H., Tilbrook, B., Tsujino, H., Tubiello, F., van der Werf, G. R., Walker, A. P., Wanninkhof, R., Whitehead, C., Willstrand Wranne, A., Wright, R., Yuan, W., Yue, C., Yue, X., Zaehle, S., Zeng, J., and Zheng, B.: Global Carbon Budget 2022, Earth Syst. Sci. Data, 14, 4811–4900,, 2022b. a, b, c, d, e, f, g

Frölicher, T. L., Rodgers, K. B., Stock, C. A., and Cheung, W. W. L.: Sources of uncertainties in 21st century projections of potential ocean ecosystem stressors, Global Biogeochem. Cy., 30, 1224–1243,, 2016. a

Galbraith, E. D. and Skinner, L. C.: The Biological Pump During the Last Glacial Maximum, Annu. Rev. Mar. Sci., 12, 559–586,, 2020. a

Gangstø, R., Gehlen, M., Schneider, B., Bopp, L., Aumont, O., and Joos, F.: Modeling the marine aragonite cycle: changes under rising carbon dioxide and its role in shallow water CaCO3 dissolution, Biogeosciences, 5, 1057–1072,, 2008. a

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Baranova, O. K., Zweng, M. M., Reagan, J. R., and Johnson, D. R.: World ocean atlas 2013. Volume 4, Dissolved inorganic nutrients (phosphate, nitrate, silicate), U.S. Department of Commerce, National Oceanic and Atmospheric Administration, National Environmental Satellite, Data and Information Service [data set],, 2013. a

Garcia, H. E., Locarnini, R. A., Boyer, T. P., Antonov, J. I., Baranova, O. K., Zweng, M. M., Reagan, J. R., and Johnson, D. R.: World Ocean Atlas 2013, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate, silicate), edited by: Levitus, S., Technical Ed.: Mishonov, A., Tech. Rep., NOAA Atlas NESDIS 76,, 2014. a, b, c

Garcia, H. E., Weathers, K. W., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 3: Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation, Technical Editor: Mishonov, A., Tech. Rep., NOAA Atlas NESDIS 83, (last access 26 July 2023), 2019a. a, b

Garcia, H. E., Weathers, K. W., Paver, C. R., Smolyar, I., Boyer, T. P., Locarnini, R. A., Zweng, M. M., Mishonov, A. V., Baranova, O. K., Seidov, D., and Reagan, J. R.: World Ocean Atlas 2018, Volume 4: Dissolved Inorganic Nutrients (phosphate, nitrate, silicate), Technical Editor: Mishonov, A., Tech. Rep., NOAA Atlas NESDIS 84, 2019b. a, b, c, d

Gehlen, M., Bopp, L., Emprin, N., Aumont, O., Heinze, C., and Ragueneau, O.: Reconciling surface ocean productivity, export fluxes and sediment composition in a global biogeochemical ocean model, Biogeosciences, 3, 521–537,, 2006. a

Geider, R. J. and La Roche, J.: The role of iron in phytoplankton photosynthesis, and the potential for iron-limitation of primary productivity in the sea, Photosynth. Res., 39, 275–301,, 1994. a

Geider, R. J., MacIntyre, H. L., and Kana, T. M.: A dynamic regulatory model of phytoplanktonic acclimation to light, nu- trients, and temperature, Limnol. Oceanogr., 43, 679–694,, 1998. a, b, c, d, e, f

Gent, P. and McWilliams, J.: Isopycnal Mixing in Ocean Circulation Models, J. Phys. Oceanogr., 20, 150–155, 1990. a

Gloege, L., McKinley, G. A., Landschützer, P., Fay, A. R., Frölicher, T. L., Fyfe, J. C., Ilyina, T., Jones, S., Lovenduski, N. S., Rodgers, K. B., Schlunegger, S., and Takano, Y.: Quantifying Errors in Observationally Based Estimates of Ocean Carbon Sink Variability, Global Biogeochem. Cy., 35, e2020GB006788,, 2021. a, b, c

Goris, N., Tjiputra, J. F., Olsen, A., Schwinger, J., Lauvset, S. K., and Jeansson, E.: Constraining Projection-Based Estimates of the Future North Atlantic Carbon Uptake, J. Climate, 31, 3959–3978,, 2018. a

Gray, A., Johnson, K. S., Bushinsky, S. M., Riser, S. C., Russell, J. L., Talley, L. D., Wanninkhof, R., Williams, N. L., and Sarmiento, J. L.: Autonomous biogeochemical floats detect significant carbon dioxide outgassing in the high-latitude Southern Ocean., Geophys. Res. Lett., 45, 9049–9057,, 2018. a

Gregor, L. and Gruber, N.: OceanSODA-ETHZ: a global gridded data set of the surface ocean carbonate system for seasonal to decadal studies of ocean acidification, Earth Syst. Sci. Data, 13, 777–808,, 2021. a, b

Gregor, L., Lebehot, A. D., Kok, S., and Scheel Monteiro, P. M.: A comparative assessment of the uncertainties of global surface ocean CO2 estimates using a machine-learning ensemble (CSIR-ML6 version 2019a) – have we hit the wall?, Geosci. Model Dev., 12, 5113–5136,, 2019. a

Griffies, S.: The Gent-McWilliams Skew Flux, J. Phys. Oceanogr., 28, 831–841,<0831:TGMSF>2.0.CO;2, 1998. a

Griffies, S. M., Biastoch, A., Böning, C., Bryan, F., Danabasoglu, G., Chassignet, E. P., England, M. H., Gerdes, R., Haak, H., Hallberg, R. W., Hazeleger, W., Jungclaus, J., Large, W. G., Madec, G., Pirani, A., Samuels, B. L., Scheinert, M., Gupta, A. S., Severijns, C. A., Simmons, H. L., Treguier, A. M., Winton, M., Yeager, S., and Yin, J.: Coordinated Ocean-ice Reference Experiments (COREs), Ocean Model., 26, 1–46,, 2009. a, b

Gruber, N., Gloor, M., Mikaloff Fletcher, S. E., Doney, S. C., Dutkiewicz, S., Follows, M. J., Gerber, M., Jacobson, A. R., Joos, F., Lindsay, K., Menemenlis, D., Mouchet, A., Müller, S. A., Sarmiento, J. L., and Takahashi, T.: Oceanic sources, sinks, and transport of atmospheric CO2, Global Biogeochem. Cy., 23, GB1005,, 2009. a

Gruber, N., Clement, D., Carter, B. R., Feely, R. A., van Heuven, S., Hoppema, M., Ishii, M., Key, R. M., Kozyr, A., Lauvset, S. K., Monaco, C. L., Mathis, J. T., Murata, A., Olsen, A., Perez, F. F., Sabine, C. L., Tanhua, T., and Wanninkhof, R.: The oceanic sink for anthropogenic CO2 from 1994 to 2007, Science, 363, 1193–1199,, 2019. a, b, c, d, e, f, g, h, i, j, k

Gürses, Ö.: Model code used in FESOM2.1-REcoM3 description paper, Zenodo [code],, 2023. a

Gürses, Ö., Oziel, L., Karakus, O., Sidorenko, D., Völker, C., Ye, Y., Zeising, M., Butzin, M., and Hauck, J.: Ocean biogeochemistry in the coupled ocean–sea ice–biogeochemistry model FESOM2.1–REcoM3 (0.1.0), Zenodo [data set],, 2023. a

Hatta, M., Measures, C. I., Wu, J., Roshan, S., Fitzsimmons, J. N., Sedwick, P., and Morton, P.: An overview of dissolved Fe and Mn distributions during the 2010–2011 U.S. GEOTRACES north Atlantic cruises: GEOTRACES GA03, Deep-Sea Res. Pt. II, 116, 117–129,, 2015. a

Hauck, J., Völker, C., Wang, T., Hoppema, M., Losch, M., and Wolf Gladrow, D. A.: Seasonally different carbon flux changes in the Southern Ocean in response to the southern annular mode, Global Biogeochem. Cy., 27, 1236–1245,, 2013. a, b, c, d

Hauck, J., Völker, C., Wolf‐Gladrow, D. A., Laufkötter, C., Vogt, M., Aumont, O., Bopp, L., Buitenhuis, E. T., Doney, S. C., Dunne, J., Gruber, N., Hashioka, T., John, J., Quéré, C. L., Lima, I. D., Nakano, H., Séférian, R., and Totterdell, I.: On the Southern Ocean CO2 uptake and the role of the biological carbon pump in the 21st century, Global Biogeochem. Cy., 29, 1451–1470,, 2015. a

Hauck, J., Köhler, P., Wolf-Gladrow, D., and Völker, C.: Iron fertilisation and century-scale effects of open ocean dissolution of olivine in a simulated CO2 removal experiment, Environ. Res. Lett., 11, 024007,, 2016. a

Hauck, J., Lenton, A., Langlais, C., and Matear, R.: The Fate of Carbon and Nutrients Exported Out of the Southern Ocean, Global Biogeochem. Cy., 32, 1556–1573,, 2018. a

Hauck, J., Zeising, M., Le Quéré, C., Gruber, N., Bakker, D. C. E., Bopp, L., Chau, T. T. T., Gürses, Ö., Ilyina, T., Landschützer, P., Lenton, A., Resplandy, L., Rödenbeck, C., Schwinger, J., and Séférian, R.: Consistency and Challenges in the Ocean Carbon Sink Estimate for the Global Carbon Budget, Frontiers in Marine Science, 7,, 2020. a, b, c, d, e, f, g, h, i, j

Hauck, J., Nissen, C., Landschützer, P., Rödenbeck, C., Bushinsky, S., and Olsen, A.: Sparse observations induce large biases in estimates of the global ocean CO2 sink: an ocean model subsampling experiment, Philos. T. Roy. Soc. A, 381, 20220063,, 2023. a, b

Henson, S. A., Sanders, R., Madsen, E., Morris, P. J., Le Moigne, F., and Quartly, G. D.: A reduced estimate of the strength of the ocean’s biological carbon pump, Geophys. Res. Lett., 38, L04606,, 2011. a

Henson, S. A., Laufkötter, C., Leung, S., Giering, S. L. C., Palevsky, H. I., and Cavan, E. L.: Uncertain response of ocean biological carbon export in a changing world, Nat. Geosci., 15, 248–254,, 2022. a

Ho, D. T., Law, C. S., Smith, M. J., Schlosser, P., Harvey, M., and Hill, P.: Measurements of air-sea gas exchange at high wind speeds in the Southern Ocean: Implications for global parameterizations, Geophys. Res. Lett., 33, L16611,, 2006. a

Hohn, S.: Coupling and decoupling of biogeochemical cycles in marine ecosystems, PhD thesis, Universität Bremen, Fachbereich Biologie, (last access: 14 August 2023​​​​​​​), 2009. a, b, c, d

Huang, Y., Tagliabue, A., and Cassar, N.: Data-Driven Modeling of Dissolved Iron in the Global Ocean, Frontiers in Marine Science, 9, 1–14,, 2022. a, b, c, d, e

Iida, Y., Kojima, A., Takatani, Y., Nakano, T., Sugimoto, H., Midorikawa, T., and Ishii, M.: Trends in pCO2 and sea–air CO2 flux over the global open oceans for the last two decades, J. Oceanogr., 71, 637–661,, 2015. a, b

Japan Meteorological Agency: JRA-55 based surface dataset for driving ocean-sea ice models (JRA55-do), Research Data Server [data set],, last access: 23 August 2023. a

Jin, X., Gruber, N., Dunne, J. P., Sarmiento, J. L., and Armstrong, R. A.: Diagnosing the contribution of phytoplankton functional groups to the production and export of particulate organic carbon, CaCO3, and opal from global nutrient and alkalinity distributions, Global Biogeochem. Cy., 20, GB2015,, 2006. a

Johnson, R., Strutton, P. G., Wright, S. W., McMinn, A., and Meiners, K. M.: Three improved Satellite Chlorophyll algorithms for the Southern Ocean, J. Geophys. Res.-Oceans, 118, 3694–3703,, 2013 (data available at:, last access: 26 July 2023). a, b, c, d, e

Joos, F. and Spahni, R.: Rates of change in natural and anthropogenic radiative forcing over the past 20,000 years, P. Natl. Acad. Sci. USA, 105, 1425–1430,, 2008. a

Jungclaus, J. H., Fischer, N., Haak, H., Lohmann, K., Marotzke, J., Matei, D., Mikolajewicz, U., Notz, D., and von Storch, J. S.: Characteristics of the ocean simulations in the Max Planck Institute Ocean Model (MPIOM) the ocean component of the MPI-Earth system model, J. Adv. Model. Earth Sy., 5, 422–446,, 2013. a

Juricke, S., Danilov, S., Koldunov, N., Oliver, M., and Sidorenko, D.: Ocean Kinetic Energy Backscatter Parametrization on Unstructured Grids: Impact on Global Eddy-Permitting Simulations, J. Adv. Model. Earth Sy., 12, e2019MS001855,, 2020. a, b

Karakuş, O., Völker, C., Iversen, M., Hagen, W., Wolf‐Gladrow, D., Fach, B., and Hauck, J.: Modeling the Impact of Macrozooplankton on Carbon Export Production in the Southern Ocean, J. Geophys. Res.-Oceans, 126, e2021JC017315,, 2021. a, b, c, d, e, f, g, h

Karakuş, O., Völker, C., Iversen, M., Hagen, W., and Hauck, J.: The Role of Zooplankton Grazing and Nutrient Recycling for Global Ocean Biogeochemistry and Phytoplankton Phenology, J. Geophys. Res.-Biogeo., 127, e2022JG006798,, 2022. a, b, c

Keerthi, M. G., Prend, C. J., Aumont, O., and Lévy, M.: Annual variations in phytoplankton biomass driven by small-scale physical processes, Nat. Geosci., 15, 1027–1033,, 2022. a

Keppler, L., Landschützer, P., Gruber, N., Lauvset, S. K., and Stemmler, I.: Seasonal Carbon Dynamics in the Near-Global Ocean, Global Biogeochem. Cy., 34, e2020GB006571,, 2020. a

Koeve, W., Duteil, O., Oschlies, A., Kähler, P., and Segschneider, J.: Methods to evaluate CaCO3 cycle modules in coupled global biogeochemical ocean models, Geosci. Model Dev., 7, 2393–2408,, 2014. a

Köhler, P., Abrams, J. F., Völker, C., Hauck, J., and Wolf-Gladrow, D. A.: Geoengineering impact of open ocean dissolution of olivine on atmospheric CO2, surface ocean pH and marine biology, Environ. Res. Lett., 8, 014009,, 2013. a

Koldunov, N. V., Aizinger, V., Rakowsky, N., Scholz, P., Sidorenko, D., Danilov, S., and Jung, T.: Scalability and some optimization of the Finite-volumE Sea ice–Ocean Model, Version 2.0 (FESOM2), Geosci. Model Dev., 12, 3991–4012,, 2019. a

Kriest, I. and Oschlies, A.: On the treatment of particulate organic matter sinking in large-scale models of marine biogeochemical cycles, Biogeosciences, 5, 55–72,, 2008. a, b, c

Kulk, G., Platt, T., Dingle, J., Jackson, T., Jönsson, B. F., Bouman, H. A., Babin, M., Brewin, R. J. W., Doblin, M., Estrada, M., Figueiras, F. G., Furuya, K., González-Benítez, N., Gudfinnsson, H. G., Gudmundsson, K., Huang, B., Isada, T., Kovac, Z., Lutz, V. A., Marañón, E., Raman, M., Richardson, K., Rozema, P. D., Poll, W. H. v. d., Segura, V., Tilstone, G. H., Uitz, J., Dongen-Vogels, V. v., Yoshikawa, T., and Sathyendranath, S.: Primary Production, an Index of Climate Change in the Ocean: Satellite-Based Estimates over Two Decades, Remote Sensing, 12, 826,, 2020. a, b

Kwiatkowski, L., Torres, O., Bopp, L., Aumont, O., Chamberlain, M., Christian, J. R., Dunne, J. P., Gehlen, M., Ilyina, T., John, J. G., Lenton, A., Li, H., Lovenduski, N. S., Orr, J. C., Palmieri, J., Santana-Falcón, Y., Schwinger, J., Séférian, R., Stock, C. A., Tagliabue, A., Takano, Y., Tjiputra, J., Toyama, K., Tsujino, H., Watanabe, M., Yamamoto, A., Yool, A., and Ziehn, T.: Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections, Biogeosciences, 17, 3439–3470,, 2020. a

Kwon, E. Y., Primeau, F., and Sarmiento, J. L.: The impact of remineralization depth on the air–sea carbon balance, Nat. Geosci., 2, 630–635,, 2009. a

Lacroix, F., Ilyina, T., and Hartmann, J.: Oceanic CO2 outgassing and biological production hotspots induced by pre-industrial river loads of nutrients and carbon in a global modeling approach, Biogeosciences, 17, 55–88,, 2020. a

Lacroix, F., Ilyina, T., Mathis, M., Laruelle, G. G., and Regnier, P.: Historical increases in land-derived nutrient inputs may alleviate effects of a changing physical climate on the oceanic carbon cycle, Glob. Change Biol., 27, 5491–5513,, 2021. a, b

Landschützer, P., Gruber, N., and Bakker, D. C. E.: Decadal variations and trends of the global ocean carbon sink, Global Biogeochem. Cy., 30, 1396–1417,, 2016. a, b

Large, W. G. and Yeager, S. G.: Diurnal to decadal global forcing for ocean and sea-ice models: the datasets and flux climatologies, Tech. Rep., CGD division of the National Center for Atmospheric Research, NCAR technical note, NCAR/TN-460+STR,, 2004. a, b

Large, W. G. and Yeager, S. G.: The global climatology of an interannually varying air-sea flux data set., Clim. Dynam., 33, 341–364,, 2009. a

Large, W. G., McWilliams, J. C., and Doney, S. C.: Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization, Rev. Geophys., 32, 363–403,, 1994. a

Lauderdale, J. M. and Cael, B. B.: Impact of remineralization profile shape on the air-sea carbon balance, Geophys. Res. Lett., 48, e2020GL091746,, 2021. a, b

Laufkötter, C., Vogt, M., Gruber, N., Aita-Noguchi, M., Aumont, O., Bopp, L., Buitenhuis, E., Doney, S. C., Dunne, J., Hashioka, T., Hauck, J., Hirata, T., John, J., Le Quéré, C., Lima, I. D., Nakano, H., Seferian, R., Totterdell, I., Vichi, M., and Völker, C.: Drivers and uncertainties of future global marine primary production in marine ecosystem models, Biogeosciences, 12, 6955–6984,, 2015. a, b

Lauvset, S. K., Key, R. M., Olsen, A., van Heuven, S., Velo, A., Lin, X., Schirnick, C., Kozyr, A., Tanhua, T., Hoppema, M., Jutterström, S., Steinfeldt, R., Jeansson, E., Ishii, M., Perez, F. F., Suzuki, T., and Watelet, S.: A new global interior ocean mapped climatology: the 1× 1 GLODAP version 2, Earth Syst. Sci. Data, 8, 325–340,, 2016 (data available at:, last access: 26 July 2023). a, b, c, d, e

Lee, K.: Global net community production estimated from the annual cycle of surface water total dissolved inorganic carbon, Limnol. Oceanogr., 46, 1287–1297, Wiley,, 2001. a

Lee, Z. and Marra, J. F.: The Use of VGPM to Estimate Oceanic Primary Production: A “Tango” Difficult to Dance, Journal of Remote Sensing, 2022, 9851013,, 2022. a

Lenton, A., Tilbrook, B., Law, R. M., Bakker, D., Doney, S. C., Gruber, N., Ishii, M., Hoppema, M., Lovenduski, N. S., Matear, R. J., McNeil, B. I., Metzl, N., Mikaloff Fletcher, S. E., Monteiro, P. M. S., Rödenbeck, C., Sweeney, C., and Takahashi, T.: Sea–air CO2 fluxes in the Southern Ocean for the period 1990–2009, Biogeosciences, 10, 4037–4054,, 2013. a

Le Quéré, C., Takahashi, T., Buitenhuis, E. T., Rödenbeck, C., and Sutherland, S. C.: Impact of climate change and variability on the global oceanic sink of CO2, Global Biogeochem. Cy., 24, GB4007,, 2010. a

Lévy, M., Franks, P. J. S., and Smith, K. S.: The role of submesoscale currents in structuring marine ecosystems, Nat. Commun., 9, 47–58,, 2018. a

Lewis, K. M. and Arrigo, K. R.: Ocean Color Algorithms for Estimating Chlorophyll a, CDOM Absorption, and Particle Backscattering in the Arctic Ocean, J. Geophys. Res.-Oceans, 125, e2019JC015706,, 2020. a

Lewis, K. M., van Dijken, G. L., and Arrigo, K. R.: Changes in phytoplankton concentration now drive increased Arctic Ocean primary production, Science, 369, 198–202,, 2020. a

Liao, E., Resplandy, L., Liu, J., and Bowman, K. W.: Amplification of the Ocean Carbon Sink During El Niños: Role of Poleward Ekman Transport and Influence on Atmospheric CO2, Global Biogeochem. Cy., 34, e2020GB006574,, 2020.