the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Ocean biogeochemistry in the coupled ocean–sea ice–biogeochemistry model FESOM2.1–REcoM3
Laurent Oziel
Onur Karakuş
Dmitry Sidorenko
Christoph Völker
Moritz Zeising
Martin Butzin
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.
- Article
(18361 KB) - Full-text XML
- BibTeX
- EndNote
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 Gruber, 2006). 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 Cael, 2021), as also known from paleo evidence (Galbraith and Skinner, 2020).
Global ocean biogeochemistry models (GOBMs; Fennel et al., 2022) are used to assess the global ocean carbon sink (Hauck et al., 2020), its regional patterns (Fay and McKinley, 2021), and effects of climate change and variability on the ocean carbon sink (Le Quéré et al., 2010; Hauck et al., 2013; DeVries et al., 2019; Bunsen, 2022). Through their representation of pH, the marine oxygen cycle, and phytoplankton primary production as the base of the marine food web, they also offer information about the environmental conditions for marine life and how these will develop under climate change (Bopp et al., 2013; Laufkötter et al., 2015; Kwiatkowski et al., 2020). However, modelling the marine biogeochemistry is subject to several sources of uncertainties. First, GOBMs are expensive with respect to the computational cost, due to the advection of a large number of tracers, and therefore often demand low spatial resolution. This leads to deficiencies in the representation of significant physical processes such as (sub)mesoscale currents (McWilliams, 2016), which can have large impacts on transport and mixing processes that strongly affect biological productivity (Lévy et al., 2018; Keerthi et al., 2022). Second, the descriptions of ecological interactions and of the physiology of primary and secondary producers in GOBMs are still mostly based on empirical or semi-empirical mathematical descriptions, such as the dependency of zooplankton grazing rates on prey abundance (Doney et al., 2001; Rohr et al., 2022). These contain a large number of parameters that are only partly constrained from observations, making it necessary to tune these parameters in GOBMs to some extent. Choices in these parameters can have strong effects on the biological carbon pump (e.g. Lauderdale and Cael, 2021). It has been demonstrated that the largest source of uncertainty for projections of net primary production (NPP; Tagliabue et al., 2021) comes from model uncertainty and not scenario uncertainty (Frölicher et al., 2016).
Ocean circulation models formulated on unstructured meshes have become an alternative to the existing structured global ocean models (Danilov, 2013). The Finite-Element Sea ice–Ocean Model (hereafter FESOM1.4; Wang et al., 2014) is one of the first global models with multiple resolutions designed to simulate the large-scale ocean circulation. While it has already been used in numerous applications (Sidorenko et al., 2015; Wekerle et al., 2017), another dynamical core, the Finite-volumE Sea ice–Ocean Model version 2.1 (FESOM2.1), has been developed (Danilov et al., 2017). The advantages of a finite-volume formulation are (a) better throughput and scalability as a result of a more efficient data structure (Koldunov et al., 2019), (b) the availability of clearly defined fluxes, and (c) the possibility to choose from a selection of transport algorithms, which was very limited before (Danilov et al., 2017). Furthermore, the arbitrary Lagrangian–Eulerian (ALE) vertical coordinate is introduced, which provides different types of vertical coordinates (Scholz et al., 2019). The Regulated Ecosystem Model (REcoM) is an ocean biogeochemistry model that describes the lower trophic levels of the marine ecosystem, using the plankton functional type approach. It bases its description of primary production on a physiological model for phytoplankton growth that takes into account the nutrient availability effects on photoacclimation (Geider et al., 1998) and, for diatoms, on the relative frustule weight (Hohn, 2009). One specificity of REcoM is the representation of flexible stoichiometry, which leads to a description of elemental fluxes that can deviate from the fixed Redfield ratios often used in models (Redfield et al., 1963).
Here, we document the ocean biogeochemistry in the Regulated Ecosystem Model version 3 (REcoM3), coupled to the ocean and sea ice model FESOM2.1, and assess its performance in reproducing carbon and nutrient biogeochemical fluxes and the distribution of phytoplankton and zooplankton. Our aim is to analyse the new set-up regarding the coupled model state under historical atmospheric 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.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).
We use FESOM2.1, an updated version of FESOM2.0. The updated model features include several developments, such as parallel and asynchronous output writing. An important new feature that we applied is the kinematic backscatter parameterization. This method takes into account the scales at which energy is scattered back to the resolved flow by introducing a negative viscosity term (Juricke et al., 2020). This greatly improves the simulation of eddy effects in coarse-resolution mesh set-ups (Juricke et al., 2020). The model code also includes the representation of ice shelf cavities (Timmermann et al., 2012), which has been used in regional studies with FESOM1.4–REcoM2 (Nissen et al., 2022). Ice shelf cavities are, however, not used in this study. Isoneutral tracer diffusion (Redi, 1982) and the Gent-McWilliams (GM; Gent and McWilliams, 1990; Griffies, 1998) eddy stirring parameterization are applied. Both GM and the Redi (1982) scheme are scaled with horizontal resolution, with a maximum value of 2000 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 Beckmann, 2004).
Regarding the vertical discretization, FESOM2.1 is formulated with an arbitrary Lagrangian–Eulerian (ALE) scheme, which is a synthesis of different types of vertical coordinates. In the model configuration used here, we apply a full free-surface formulation and thus permit the vertical movement of the surface and of all other layers (referred to as zstar; Scholz et al., 2019). This drastically improves the tracer conservation properties (Campin et al., 2004). Partially filled cells are used at the ocean floor, resulting in a smoother representation of the bathymetry.
The sea ice component (Finite-Element Sea Ice Model, FESIM version 2) solves for sea ice concentration, ice and snow thickness, and ice drift velocity (Danilov et al., 2015). It is discretized on the same unstructured horizontal mesh as the ocean model. The elastic viscous plastic solver and flux-corrected transport scheme are used for sea ice advection (Danilov et al., 2015). The formulation of sea ice thermodynamics follows the work of Timmermann et al. (2009).
2.1.2 Biogeochemistry model REcoM3
REcoM3 is a water column biogeochemistry and ecosystem model, which incorporates cycles of carbon and nutrients (nitrogen, iron, and silicon) with varying intracellular stoichiometry in phytoplankton, zooplankton, and detritus (see the Appendix for detailed description and equations). Starting from the work by Schartau et al. (2007), REcoM was first used to describe carbon overconsumption in mesocosm experiments. After being coupled to the ocean and sea ice model MITgcm (Marshall et al., 1997), the previous version (REcoM2) with two phytoplankton classes, one zooplankton and one detritus class, was applied to study the cycling of marine carbon on present-day (Hauck et al., 2013, 2018) and glacial timescales (Du et al., 2022; Völker and Köhler, 2013), as well as the marine iron cycle (e.g. Völker and Tagliabue, 2015; Tagliabue et al., 2016; Ye and Völker, 2017; Pagnone et al., 2019). Moreover, REcoM2 was employed in assessments on the efficiency of ocean alkalinity enhancement (Köhler et al., 2013; Hauck et al., 2016), in data assimilation studies (Pradhan et al., 2019), and as a test bed for model development; e.g. for the development of a parameterization of iron-ligand binding based on pH (Ye et al., 2020), among others. Simultaneously, REcoM2 was coupled to FESOM1.4 (Schourup-Kristensen et al., 2014). These coupled model set-ups were used either in a global configuration (e.g. Schourup-Kristensen et al., 2014; Hauck et al., 2020), with a regional focus on the Arctic or the Antarctic (Hauck et al., 2015; Schourup-Kristensen et al., 2018; Oziel et al., 2022; Nissen et al., 2022) or in regional configurations (Taylor et al., 2013; Losch et al., 2014). Recently, the model has matured to include two groups of each classes of phytoplankton, zooplankton, and detritus (REcoM3; Fig. 2).
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 ( ratios for small phytoplankton and for diatoms) to nutrient levels, ambient light, and temperature, based on the photoacclimation model by Geider et al. (1998). Si uptake by diatoms is regulated as well, based on the internal Si:N quota, following Hohn (2009). This parameterization takes into account the strong decoupling between Si and N metabolism (e.g. Claquin et al., 2002) and prescribes the observed change in Si:N ratios under Fe and N limitation. The intracellular iron pool is derived from intracellular nitrogen with a fixed Fe:N ratio, based on the fact that intracellular iron is mostly associated with the photosynthetic electron transport chain and nitrogen metabolism (Geider and La Roche, 1994; Raven, 1988). REcoM3 also includes the photodamage parameterization by Álvarez et al. (2018). Calcium carbonate production is assumed to be linearly dependent on the gross small phytoplankton production. 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 Oschlies, 2008). 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.
REcoM
-
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 (Orr, 1999), were replaced by the mocsy 2.0 scheme of Orr and Epitalon (2015). While both use the same thermodynamic equilibrium to calculate surface pCO2 and CO2 flux, mocsy 2.0 uses the faster and more accurate algorithm SolveSAPHE (Munhoven, 2013). Among other differences, it follows best-practice guides and uses recommended equilibrium constants. The gas exchange formulation is updated to that of Wanninkhof (2014), which is largely equivalent to that of Ho et al. (2006). The computed fluxes are scaled with the ice-free area.
-
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 Epitalon, 2015). 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.
-
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).
-
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).
-
A sedimentary release of iron was added to the model (Tagliabue et al., 2016); this was done in addition to the previously considered Fe input with dust deposition.
FESOM
Biogeochemical fluxes returned back to the ocean from the benthos are treated with a specific bottom boundary condition. Variable bottom topography leads to a smaller scalar control volume located at the lowermost level. This is because the scalar control volumes are obtained by connecting the areas from the elements to which they are attached (see Fig. 1 in Danilov et al., 2017). Therefore, the number of elements around a single surface node may vary with depths when it meets non-flat topography. We thus computed the control volume and associated fluxes for each node by considering all surrounding elements at different depth levels.
Forcing
Our simulation was forced by the atmospheric reanalysis JRA55-do data set (Tsujino et al., 2018) instead of the CORE-II data set (Large and Yeager, 2009) that was used in previous assessments (Schourup-Kristensen et al., 2014). JRA55-do is a blend of reanalysis data and satellite observations and has the advantage to provide regularly updated near-real-time data up to the present day, with a higher temporal (3 h) resolution. The freshwater supplied by rivers is a climatology and provided by Large and Yeager (2004) as part of the CORE-II forcing. Nutrient, carbon and alkalinity supply via river discharge is not included in the experiments described here.
2.2 Experimental set-up and data
In this study, we used a mesh with a nominal resolution of 1∘ as a background. The horizontal resolution is enhanced on the equatorial belt and in the region north of 50∘ N to match ∘ 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.
The atmospheric reanalysis data sets of JRA55‐do v.1.5.0 (Tsujino et al., 2018) were used to force the model for the period 1958–2021 (hereafter JRA55-do). A single repeating annual cycle of all forcing fields (year 1961) was used to perform the spin-up simulations and a control experiment. This is referred to as repeat-year forcing (hereafter called RYF61). We have deliberately chosen the year 1961, as it had rather neutral El Niño–Southern Oscillation conditions and also contained a low amount of anthropogenic perturbation, compared to the years of 1990 and 1991 recommended by Stewart et al. (2020).
A series of experiments were carried out in a global set-up to investigate the performance of the coupled FESOM2.1–REcoM3 model. The experiments follow the definitions used in the Global Carbon Budget (Friedlingstein et al., 2022a) and in the RECCAP (REgional Carbon Cycle Assessment and Processes, https://reccap2-ocean.github.io/, last access: 14 August 2023) projects and are summarized in Table 1. Our first experiment was forced with varying climate from the JRA55-do data set and varying atmospheric CO2 levels (hereafter referred to as A). Atmospheric CO2 mixing ratio (xCO2) values are taken from the Global Carbon Budget (Joos and Spahni, 2008; Ballantyne et al., 2012; Friedlingstein et al., 2022a). A second simulation was forced by RYF61 atmospheric reanalysis fields and a preindustrial atmospheric 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 Falkowski, 1997)(Buitenhuis et al., 2010; Moriarty and O'Brien, 2013)(Moriarty et al., 2013)Sallée et al. (2021a)(Bakker et al., 2016)(Chau et al., 2022)Chau et al. (2022)(Friedlingstein et al., 2022b; Wright et al., 2021; Schwinger et al., 2016; Lacroix et al., 2021; Berthet et al., 2019; Hauck et al., 2020; Liao et al., 2020; Doney et al., 2009; Aumont et al., 2015; Nakano et al., 2011; Urakawa et al., 2020; Long et al., 2021a; Landschützer et al., 2016; Rödenbeck et al., 2022; Chau et al., 2022; Gloege et al., 2021; Zeng et al., 2014; Iida et al., 2015; Gregor and Gruber, 2021)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 Vallis, 2007; Storkey et al., 2018). Similar issues are found in comparable current systems, such as the Kuroshio and Malvina systems. It is, however, surprising that in general FESOM is far too saline at the surface and is on average 0.3 psu saltier than the climatology. The reason for this bias could be imperfections in the river discharge from CORE-II forcing and the relatively low surface salinity restoring, which uses a piston velocity of 50 m/300 d in the simulations. In most of the ocean, the sea surface temperature (SST) and sea surface salinity (SSS) differences act in an opposite manner on buoyancy. Hence, the increase or decrease in SST is accompanied by an increase or decrease in SSS. The only exception is the Indian Ocean, where east and west in simulation A become less and more buoyant, respectively (Fig. 3).
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 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).
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ölker, 2017; Pagnone et al., 2019). Furthermore, the intensity and extension of dust plumes vary between modelled dust deposition fluxes (e.g. Myriokefalitakis et al., 2018). The field of dust deposition by Albani et al. (2014), used in our model to calculate aeolian iron input, is within the range of modern estimates but surely contains some uncertainties. Despite the overall amplitude of the patterns in the dissolved iron being too strong, especially in the regions of high dust deposition, the model is able to reproduce the main regions in which iron availability limits phytoplankton productivity (Moore et al., 2013), namely the subpolar Southern Ocean, the equatorial and North Pacific, and also to some extent the seasonal iron limitation in the subpolar North Atlantic (Nielsdóttir et al., 2009), but overestimates iron limitation in the subtropical South Pacific.
3.2.2 Modelled versus satellite-based phytoplankton biomass and productivity
We first evaluated the spatial distributions of the modelled chlorophyll a concentration (Fig. 9), which is an indicator for phytoplankton biomass, and vertically integrated the net primary production (NPP, Fig. 10). We compared chlorophyll a concentrations obtained from FESOM2.1–REcoM3 simulations, averaged from 2012 to 2021, with an ocean colour remote sensing merged data set (from the Ocean Colour Climate Change Initiative, OC-CCI; Sathyendranath et al., 2019), averaged from 1998 to 2019. We also compared the modelled NPP with satellite estimations, such as the Vertically Generalized Production Model (VGPM; Behrenfeld and Falkowski, 1997; Fig. 10) and the updated carbon-based productivity model (CbPM; Westberry et al., 2008; see Fig. A1 in the Appendix). VGPM is a chlorophyll-based algorithm that can be considered to be a standard NPP estimation from ocean colour for the last 20 years (Lee and Marra, 2022). CbPM uses spectrally resolved light attenuation and is based on a semi-analytical algorithm (Garver–Siegel–Maritorena, GSM; Maritorena et al., 2002). All data sets were also compared along a latitudinal distribution with an improved chlorophyll a algorithm for the Southern Ocean (Johnson et al., 2013; Fig. 11).
The results for chlorophyll a and NPP obtained here are comparable to those presented by Schourup-Kristensen et al. (2014). Over large parts of the global ocean, the mean surface chlorophyll a concentrations are in agreement with the observations (Fig. 9c and d). Yet, there are regional differences. In temperate latitudes, the modelled chlorophyll a concentrations are somewhat higher than observed, while the subtropical gyres show concentrations slightly lower than the observations. The comparison of modelled and observation-based satellite estimates of chlorophyll a yielded a correlation of 0.66 and an RMSE of 0.38 mg m−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 Falkowski, 1997), 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; Mitchell, 1992; Mustapha et al., 2012). Some recent advances used local parameterizations with in situ data, which resulted in much lower productivity levels in those coastal areas (Lewis et al., 2020; Lewis and Arrigo, 2020). Generally, the NPP and chlorophyll differences compared to satellite-based estimates could also be linked to model deficiencies, such as a coarse model resolution and associated weak upwelling, missing complexity in the simulated phytoplankton classes, and also the so-far unconsidered nutrient input from terrigenous sources.
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).
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.
3.2.3 Modelled versus MAREDAT zooplankton biomass
In REcoM3, the small zooplankton group is widely spread in the global ocean and the highest biomass occurs in high-productivity regions (Fig. 13a). The macrozooplankton is present in the high latitudes (Fig. 13b), since it is parameterized as a polar macrozooplankton group (Karakuş et al., 2021). We compare the latitudinal distribution of the integrated modelled zooplankton biomass with gridded global zooplankton biomass data from MAREDAT (Buitenhuis et al., 2010; Moriarty et al., 2013; Moriarty and O'Brien, 2013). The simulated biomass of small and total zooplankton reproduces MAREDAT-derived biomass reasonably well in low to midlatitudes but underestimates biomass in the polar regions (Fig. 13c). The underestimation of zooplankton biomass in the northern high latitudes may be related to an underestimation of primary production in the same region. In agreement with the MAREDAT data set (Moriarty et al., 2013), macrozooplankton is not present in low latitudes.
3.2.4 Synthesis
The modelled biogeochemical fluxes were compared to the previous version FESOM1.4–REcoM2 and observational studies (Table 3). Modelled global NPP is higher in FESOM2.1–REcoM3 than in FESOM1.4–REcoM2 but still lower than in satellite-based estimates. The estimate is comparable to other global modelling studies (Schneider et al., 2008; Séférian et al., 2020). Export production (EP) is slightly higher in FESOM2.1–REcoM3 than in the previous version and falls within the observational range previously documented in the literature for both the global ocean and the Southern Ocean. For the global ocean, FESOM2.1–REcoM3 NPP and EP estimations remained at the lower end of the range, despite a slight increase in NPP. A more detailed description of zooplankton results in more efficient nutrient recycling and can thus increase NPP by 25 % (see also explanation in Sect. 3.2.2; Karakuş et al., 2022). In the Southern Ocean, estimations of NPP and EP remained very close to observation-based estimates. Maybe the most noticeable change between the two model versions is the substantial increase in opal export, which increased by a factor of 4 in the Southern Ocean, passing from the lower to the higher end of the observational range of an earlier review (Dunne et al., 2007), and lies above an updated estimate (Tréguer et al., 2021). This is due to an increase in the relative contribution of diatoms to the total NPP in high latitudes (Fig. 11).
(Behrenfeld and Falkowski, 1997)(Westberry et al., 2008)(Schneider et al., 2008)(Kulk et al., 2020)Séférian et al., 2020(Schlitzer, 2004)(Dunne et al., 2007)(Henson et al., 2011)(Siegel et al., 2014)Dunne et al., 2007(Tréguer et al., 2021)Lee, 2001; Jin et al., 2006Gangstø et al., 2008; Berelson et al., 2007; Dunne et al., 2007; Battaglia et al., 2016; Gehlen et al., 2006(Carr et al., 2006)(Behrenfeld and Falkowski, 1997)(Schlitzer, 2002; Nevison et al., 2012)(Dunne et al., 2007)(Dunne et al., 2007)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).
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 McKinley, 2021). Therefore, we chose one product (Chau et al., 2022) and focus on the comparison of the spatial pattern with our model. We further evaluate the temporal evolution of 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.
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.
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.
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.
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 (Sundquist, 1985; 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)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 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.
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).
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:
where S is the volumetric concentration of a state variable, U is the 3-dimensional advection velocity, and κ is the diffusivity. The term SMS(S) represents the biogeochemical sources minus sinks. The slow-sinking detritus class is assumed to sink with a velocity, which increases linearly with depth as a first-order description of the shift to larger and faster-sinking particles with increasing depth (Kriest and Oschlies, 2008). A constant sinking rate is applied to the fast-sinking detritus class. REcoM3 has 28 oceanic and four explicit benthic state variables (Tables A1 and A2).
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:
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 and (Table A5).
Dissolved silicic acid (DSi)
Silicon assimilation (Si assimilation) increases when biogenic silica from one of the two detritus classes dissolves.
The state variables PhyCdia, DetSi, and DetZ2Si are listed in Table A1. The temperature-dependent remineralization rate of silicon () 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.
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 (, , and ) 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 (; ) 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ölker, 2011), so we assume instantaneous equilibrium between free iron and free ligand (L′), which is computed using a constant , by solving the following:
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 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 Epitalon, 2015).
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. and 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).
The state variables PhyCsmall, PhyCdia, DON, DetCalc, and DetZ2Calc are listed in Table A1. The N assimilation ( and ) 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. and are grazing terms and explained in Sect. A4.2.
A1.3 Phytoplankton
Nitrogen
The phytoplankton nitrogen pools increase through N assimilation. The assimilation process is assumed to be proportional to carbon biomass, with a carbon-specific uptake rate that depends on the C:N ratio of phytoplankton and the external DIN concentration (Geider et al., 1998). Excretion of biogenic nitrogen to semi-labile DON drains the pool. At a high intracellular C:N ratio, the excretion is downregulated. Aggregation and grazing by the two zooplankton groups transfer nitrogen to the zooplankton and detritus pools.
The state variables PhyCsmall, PhyNsmall, PhyCdia, and PhyNdia are listed in Table A1. The N assimilation ( and ) is explained in Sect. A3.4. The constant excretion rate constant () 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 (; ) 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 (, , , and ) are explained in Sect. A4.2.
Carbon
The carbon biomass of small phytoplankton and diatoms increases as a result of carbon assimilation during photosynthesis. Loss terms include excretion of DOC, which is limited by the availability of proteins as in the nitrogen pool, respiration, aggregation, and grazing.
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 (; Table A8) is downregulated by the limiter factor (; ) when the N:C ratio becomes too high (Eq. A45). Phytoplankton aggregation (Agg) is calculated in Eq. (A42). Grazing terms (, , , and ) are explained in Sect. A4.2. , is used to convert the grazing units from millimoles of N to millimoles of C.
CaCO3
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.
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 (; Table A8) is downregulated by the limiter factor (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 ( and ) are explained in Sect. A4.2. 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.
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 (; Table A8) is downregulated by the limiter factor (Eq. A45) when the N:C ratio becomes too high. Grazing terms ( and ) are explained in Sect. A4.2. The intracellular ratio between diatom silicon and nitrate is defined as .
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.
The state variables PhyCsmall, PhyCdia, PhyChlsmall, and PhyChldia are listed in Table A1. The chlorophyll a synthesis (; ) and the aggregation (Agg) terms are calculated in Eqs. (A39) and (A42), respectively. The degradation parameters (; ) are given in Table A8. Grazing terms (, , , and ) are explained in Sect. A4.2. The conversion factor from millimoles of N to milligrams of Chl a is defined as .
A1.4 Zooplankton
Nitrogen
Both zooplankton classes increase their nitrogen biomass via grazing on phytoplankton and detritus, while mortality and excretion of DON reduce it. Macrozooplankton further feeds on small zooplankton and releases nitrogen via fecal pellet production.
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 ( and ) 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 ( and ) are given in Table A8.
Carbon
The zooplankton carbon biomass increases with carbon uptake via grazing and decreases through carbon losses through mortality, respiration, and carbon excretion to the semi-labile DOC pool. Macrozooplankton further gains carbon by grazing on small zooplankton and loses it via fecal pellet production.
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 (, , , , , , , , 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 () are in Table A8. The grazing flux in terms of nitrogen biomass is converted to carbon biomass using the respective intracellular C:N ratios (, , , , , and ), where , , , , , and . 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
Nitrogen
Detrital nitrogen pool increases as a result of sloppy feeding and mortality. Sloppy feeding is outlined as a function of grazing fluxes and grazing efficiency of macrozooplankton. In other words, the grazed phytoplankton partly goes to the macrozooplankton biomass, depending on the grazing efficiency. The phytoplankton aggregation contributes only to slow-sinking detritus. Fecal pellet production is defined only for macrozooplankton group. Detritus is degraded to DON, based on temperature and a remineralization rate.
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 (, , , , , , , , 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).
Carbon
Detrital carbon sources are associated with sloppy feeding, aggregation of phytoplankton, mortality of small zooplankton, and fecal pellet production by macrozooplankton. Degradation of DetC and DetZ2C to DOC is the only loss term.
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 (, , , , , , , , 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 , , , , , and are used to convert the units from mmol N to mmol C.
Silica
Biogenic detrital silica increases with excretion fluxes from diatoms to detritus, aggregation, and grazing and decreases with silica dissolution from DetSi and DetZ2Si.
The state variables DiaSi, DetSi, and DetZ2Si are listed in Table A1. The constant excretion rate (; Table A8) is downregulated by the limiter factor (Eq. A45) when the N:C ratio becomes too high. The remineralization rates (), the aggregation (Agg), and the grazing on diatoms ( and ) are calculated in Eqs. (A35), (A42), and (A55), respectively. The intracellular ratio between diatom silicon and carbon is defined as .
CaCO3
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.
The state variables PhyCalc, DetCalc, and DetZ2Calc are listed in Table A1. The constant excretion rate (; Table A8) is downregulated by the limiter factor (Eq. A45) when the N:C ratio becomes too high. The respiration (rsmall), the aggregation (Agg), and the grazing on small phytoplankton ( and ) are calculated in Eqs. (A38), (A42), and (A54), respectively. The ratio .
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.
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:
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 Epitalon, 2015).
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.
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 (, , , and ), 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 ( and ; 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.
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 (, , , and ), 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 ( and ; 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.
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örtner, 2016) 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).
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 (; Table A5) is calculated following Maerz et al. (2020) but with a minimum dissolution rate.
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:
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 and varies as a result of photoacclimation. The apparent maximum photosynthetic rate ( and ) is defined below.
The value of and is listed in Table A7. The limitation terms (, , , , and ) 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:
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 (Hohn, 2009). The limiter function ( and ) is described in Eq. (A45), and the N assimilation rate ( and ) is calculated in Eq. (A40).
A3.3 Chlorophyll a synthesis
The chlorophyll synthesis rate ( and ; 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.
The N assimilation ( and ) is computed in Eq. (A40). The conversion factor of the maximum Chl:N ratio ( and ) 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.
A3.4 N and Si assimilation
Nitrogen
The C-specific N assimilation rate is a function of the maximum rate of C-specific photosynthesis and DIN concentration. N assimilation depends on the DIN concentration in seawater via Michaelis–Menten kinetics. The maximum photosynthetic rate (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 (see Eq. A45). Nitrogen assimilation is downregulated at high intracellular N:C ratio.
, , , , , and are listed in Table A7. The maximum rate of photosynthesis ( and ) is given in Eqs. (A37). and are described in Eq. (A45). DIN corresponds to the in situ concentration.
Silicon
The building of a silica frustule of diatoms requires silicate uptake. The C-specific Si assimilation rate (VSi) is a function of a factor for C-specific N uptake (), a rate constant of C-specific photosynthesis (), 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.
The scaling factor for the N uptake (), the maximum rate constant of C-specific photosynthesis (), 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 ( and ) 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 (). When the nutrient limitation is high (i.e. low ), the aggregation rate increases in the model.
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 (, , and ) 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).
Here, 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 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. is one when the current (i.e. Redfield ratio; 16N:106C) and zero for (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).
The limitation function for the quota regulation is calculated with Eq. (A44). and form the current intracellular nitrogen quota for small phytoplankton and diatoms, respectively. Dimensionless constants , , and are listed in Table A6.
The limiter downregulates the Si assimilation of diatoms when the intracellular silicon quota (Si:C) becomes too high. is one when the current and zero for . 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:
Dimensionless constants and qSi:Cmax are listed in Table A6.
Carbon fixation and aggregation loss in diatoms are further downregulated by a factor (; 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). is zero when the current and one for .
The dimensionless constants 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.
The variable DFe is listed in Table A1. The half-saturation constants ( and ) 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 ( and ) is described as a function of the intracellular nitrogen quota ( and ) with growth ending at a minimum quota ( and ).
Dimensionless constants , , and 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 τ.
The timescale for respiration (τ) is listed in Table A7. The temperature dependence (fT) is calculated in Eq. (A33). The ratios are defined as and .
Macrozooplankton
The daily respiration rate constant of macrozooplankton (rzoo2; Table A5) is modelled following Karakuş et al. (2021).
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:
() 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:
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 (, , and ) 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.
The FESOM2.1–REcoM3 source code is available at https://github.com/FESOM/fesom2/tree/fesom2.1_recom (last access: 31 December 2022). The version of FESOM2.1–REcoM3 used for this paper can be found at https://doi.org/10.5281/zenodo.7502419 (Gürses, 2023). A manual is available at https://recom.readthedocs.io/en/latest/ (last access: 23 July 2023).
The GLODAPv2 mapped product of dissolved inorganic carbon and alkalinity can be found at https://glodap.info/index.php/mapped-data-product/ (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 https://doi.org/10.7289/V5J67DWD (Garcia et al., 2013). The oxygen data set from the World Ocean Atlas 2018 is available at https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/bin/woa18oxnu.pl?parameter=o (last access: 26 July 2023) (Garcia et al., 2019a). The OC-CCI data set is available at https://www.oceancolour.org/thredds/catalog-cci.html (last access 26 July 2023) (Sathyendranath et al., 2019). The Southern Ocean chlorophyll a concentration data set is available at https://imos.org.au/facilities/srs/oceancolour (last access 26 July 2023) (Johnson et al., 2013). CbPM (Westberry et al., 2008) and VGPM (Behrenfeld and Falkowski, 1997) products of net primary production are available at http://sites.science.oregonstate.edu/ocean.productivity/inputData.php (last access 26 July 2023). MAREDAT products can be found at https://doi.org/10.1594/PANGAEA.779970 (Buitenhuis et al., 2012), https://doi.org/10.1594/PANGAEA.785501 (O'Brien and Moriarty, 2012), and https://doi.org/10.1594/PANGAEA.777398 (Moriarty, 2012). Global maps of trends and climatological fields are available at https://doi.org/10.5281/zenodo.5776180 (Sallée et al., 2021b). The Surface Ocean CO2 Atlas (SOCAT) was taken from https://socat.info/index.php/previous-versions/ (last access: 26 July 2023) (Bakker et al., 2016). The Global Ocean Surface Carbon is available at https://data.marine.copernicus.eu/product/MULTIOBS_GLO_BIO_CARBON_SURFACE_REP_015_008/description (last access: 26 July 2023) (Chau et al., 2022). Annual air–sea CO2 flux data is freely available at https://globalcarbonbudgetdata.org (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 Yeager, 2004) are available online. The JRA-55-based surface data set for driving ocean–sea ice models (JRA55-do) is taken from https://climate.mri-jma.go.jp/pub/ocean/JRA55-do/ (last access: 23 August 2023, Japan Meteorological Agency, 2023). The model results are available at https://doi.org/10.5281/zenodo.8276875 (Gürses et al., 2023).
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.
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.
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).
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, https://doi.org/10.1002/2013MS000279, 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, https://doi.org/10.1029/2017GB005850, 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, https://doi.org/10.1029/2009GL040046, 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, https://doi.org/10.1016/j.pocean.2010.06.003, 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, https://doi.org/10.1029/1999GB001238, 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, https://doi.org/10.1029/2001GB001745, 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, https://doi.org/10.5194/gmd-8-2465-2015, 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, https://doi.org/10.1016/S0967-0645(96)00075-6, 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, https://doi.org/10.5194/essd-8-383-2016, 2016 (data available at: https://socat.info/index.php/previous-versions/, 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, https://doi.org/10.1038/nature11299, 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, https://doi.org/10.5194/bg-13-2823-2016, 2016. a
Behrenfeld, M. J. and Falkowski, P. G.: A consumer’s guide to phytoplankton primary productivity models, Limnol. Oceanogr., 42, 1479–1491, https://doi.org/10.4319/lo.1997.42.1.0001, 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, https://doi.org/10.1029/2006gb002803, 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, https://doi.org/10.1029/2019MS001644, 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, https://doi.org/10.1016/j.pocean.2013.12.008, 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, https://doi.org/10.5194/bg-10-6225-2013, 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, https://doi.org/10.1038/s41467-022-27979-5, 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, https://doi.org/10.1038/s41586-019-1098-2, 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, https://doi.org/10.1016/S0967-0637(00)00099-6, 2001. a
Buitenhuis, E., Rivkin, R. B., Séailley, S., and Le Quéré, C.: Biogeochemical fluxes through microzooplankton, Global Biogeochem. Cy., 24, GB4015, https://doi.org/10.1029/2009GB003601, 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], https://doi.org/10.1594/PANGAEA.779970, 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, https://epic.awi.de/id/eprint/54788/ (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, https://doi.org/10.1029/2019GB006176, 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, https://doi.org/10.1111/gcb.13375, 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, https://doi.org/10.1016/S1463-5003(03)00009-X, 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, https://doi.org/10.1017/9781009157896.007, 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, https://doi.org/10.1016/j.dsr2.2006.01.028, 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, https://doi.org/10.5194/bg-19-1087-2022, 2022 (data available at: https://data.marine.copernicus.eu/product/MULTIOBS_GLO_BIO_CARBON_SURFACE_REP_015_008/description, 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, https://doi.org/10.1046/j.1529-8817.2002.t01-1-01220.x, 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, https://doi.org/10.5194/bg-10-1849-2013, 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, https://doi.org/10.1029/2017GB005710, 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, https://doi.org/10.1029/2021RG000736, 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, https://doi.org/10.1016/j.ocemod.2013.10.005, 2014. a
Danilov, S.: Ocean modeling on unstructured meshes, Ocean Model., 69, 195–210, https://doi.org/10.1016/j.ocemod.2013.05.005, 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, https://doi.org/10.5194/gmd-8-1747-2015, 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, https://doi.org/10.5194/gmd-10-765-2017, 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, https://doi.org/10.1016/S0304-4203(99)00022-5, 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, https://doi.org/10.1073/pnas.1900371116, 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, https://doi.org/10.5670/oceanog.2001.10, 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, https://doi.org/10.1016/j.dsr2.2008.12.006, 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, https://doi.org/10.1029/2022GL099048, 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, https://doi.org/10.1029/2006GB002907, 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, https://doi.org/10.1029/2004GL020216, 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, https://doi.org/10.1029/2021GL095325, 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, https://doi.org/10.5194/essd-13-4693-2021, 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, https://doi.org/10.1038/s43586-022-00154-2, 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, https://doi.org/10.1034/j.1600-0889.2003.01461.x, 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, https://doi.org/10.5194/essd-14-1917-2022, 2022a (data available at: https://globalcarbonbudgetdata.org, 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, https://doi.org/10.5194/essd-14-4811-2022, 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, https://doi.org/10.1002/2015GB005338, 2016. a
Galbraith, E. D. and Skinner, L. C.: The Biological Pump During the Last Glacial Maximum, Annu. Rev. Mar. Sci., 12, 559–586, https://doi.org/10.1146/annurev-marine-010419-010906, 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, https://doi.org/10.5194/bg-5-1057-2008, 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], https://doi.org/10.7289/V5J67DWD, 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, https://doi.org/10.7289/V5J67DWD, 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, https://www.ncei.noaa.gov/access/world-ocean-atlas-2018/bin/woa18oxnu.pl?parameter=o (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, https://doi.org/10.5194/bg-3-521-2006, 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, https://doi.org/10.1007/BF00014588, 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, https://doi.org/10.4319/lo.1998.43.4.0679, 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, https://doi.org/10.1029/2020GB006788, 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, https://doi.org/10.1175/JCLI-D-17-0564.1, 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, https://doi.org/10.1029/2018GL078013, 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, https://doi.org/10.5194/essd-13-777-2021, 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, https://doi.org/10.5194/gmd-12-5113-2019, 2019. a
Griffies, S.: The Gent-McWilliams Skew Flux, J. Phys. Oceanogr., 28, 831–841, https://doi.org/10.1175/1520-0485(1998)028<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, https://doi.org/10.1016/j.ocemod.2008.08.007, 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, https://doi.org/10.1029/2008GB003349, 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, https://doi.org/10.1126/science.aau5153, 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], https://doi.org/10.5281/zenodo.7502419, 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], https://doi.org/10.5281/zenodo.8276875, 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, https://doi.org/10.1016/j.dsr2.2014.07.005, 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, https://doi.org/10.1002/2013GB004600, 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, https://doi.org/10.1002/2015GB005140, 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, https://doi.org/10.1088/1748-9326/11/2/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, https://doi.org/10.1029/2018GB005977, 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, https://doi.org/10.3389/fmars.2020.571720, 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, https://doi.org/10.1098/rsta.2022.0063, 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, https://doi.org/10.1029/2011GL046735, 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, https://doi.org/10.1038/s41561-022-00927-0, 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 param