FOCI-MOPS v1 – integration of marine biogeochemistry within the Flexible Ocean and Climate Infrastructure version 1 (FOCI 1) Earth system model

The consideration of marine biogeochemistry is essential for simulating the carbon cycle in an Earth system model. Here we present the implementation and evaluation of a marine biogeochemical model, the Model of Oceanic Pelagic Stoichiometry (MOPS) in the Flexible Ocean and Climate Infrastructure (FOCI) climate model. FOCI-MOPS enables the simulation of marine biological processes, i.e. the marine carbon, nitrogen, and oxygen cycles with prescribed or prognostic atmospheric CO₂ concentration. A series of experiments covering the historical period (1850–2014) were performed following the DECK (Diagnostic, Evaluation and Characterization of Klima) and CMIP6 (Coupled Model Intercomparison Project 6) protocols. Overall, modelled biogeochemical tracer distributions and fluxes, transient evolution in surface air temperature, air–sea CO₂ fluxes, and changes in ocean carbon and heat contents are in good agreement with observations. Modelled inorganic and organic tracer distributions are quantitatively evaluated by statistically derived metrics. Results of the FOCI-MOPS model, including sea surface temperature, surface pH, oxygen (100–600 m), nitrate (0–100 m), and primary production, are within the range of other CMIP6 model results. Overall, the evaluation of FOCI-MOPS indicates its suitability for Earth climate system simulations.

1 Introduction

The strongest anthropogenic forcing on the Earth system during the last century has been a rise in atmospheric CO₂ concentrations due to anthropogenic CO₂ emissions (Shukla et al., 2019). About half of those emissions are currently taken up by the terrestrial biosphere and the ocean (Friedlingstein et al., 2020; Sabine et al., 2004; Gruber et al., 2019) in about equal proportions. The anthropogenic carbon is taken up by the ocean mostly due to the physical–chemical processes of the solubility pump (Sarmiento and Gruber, 2002) and is taken up on land by increased net primary productivity (Arneth et al., 2010). In addition to the uptake of anthropogenic carbon, the natural carbon fluxes are perturbed by climate change. In the ocean, the increase in seawater temperature directly decreases the solubility of CO₂. Global warming also leads to changes in ocean circulation; for instance, shifting wind patterns might change the Southern Ocean upwelling and increase natural CO₂ outgassing (Le Quéré et al., 2007). The increased natural outgassing in the Southern Ocean is also shown in modelling studies (Zickfeld et al., 2007; Tjiputra et al., 2010). The chemical capacity of the ocean to take up CO₂ decreases with increasing CO₂ concentrations in seawater (Friedlingstein et al., 2006; Riebesell et al., 2009; Fassbender et al., 2017). On land, increasing temperatures limit plant growth in low latitudes and enhance the decomposition of organic matter (Pugnaire et al., 2019; Lin et al., 2010; Sarmiento and Gruber, 2002). Those mechanisms are expected to lead to a weakening of the terrestrial and marine sinks for the extra carbon arising from human activity, but a detailed quantitative understanding is still lacking.

For a comprehensive investigation of climate–carbon cycle interactions and possible feedbacks, the implementation of ocean biogeochemistry in climate models is crucial. While existing global ocean biogeochemical models simulate surface ocean pCO₂ reasonably well, there are still discrepancies between the model results and data products concerning oceanic CO₂ sink estimates (Hauck et al., 2020). In order to improve our understanding of the Earth system, continuous development of the ocean in Earth system models is required. This includes an adequate representation of the marine carbon uptake variability on the finite atmospheric CO₂ pool and hence on climate.

A new climate model, the Flexible Ocean and Climate Infrastructure (FOCI), has been successfully developed (Matthes et al., 2020). The model consists of a fully coupled atmosphere–ocean–sea ice general circulation model and includes a land model, the Jena Scheme for Biosphere–Atmosphere Coupling in Hamburg (JSBACH; Brovkin et al., 2009; Reick et al., 2013); plus options for interactive stratospheric chemistry in the atmosphere, the ECHAM6.3-HAM2.3-MOZ1.0 (ECHAM6-HAMMOZ; Schultz et al., 2018); and the option for regional grid refinement in the ocean, the Adaptive Grid Refinement In Fortran package (AGRIF; Debreu et al., 2008) for the Nucleus for European Modelling of the Ocean (NEMO; Madec and the NEMO team, 2016). Here we present the implementation of the marine biogeochemical model component, the Model of Oceanic Pelagic Stoichiometry (MOPS; Kriest and Oschlies, 2015), into FOCI. MOPS enables the simulation of marine biological processes, i.e. the marine carbon, nitrate (NO₃), phosphate (PO₄), and oxygen (O₂) cycles. MOPS features a smaller number of prognostic variables than other Coupled Model Intercomparison Project Phase 6 (CMIP6) models (Séférian et al., 2020) (e.g. it does not simulate iron and silicate), which makes it computationally a comparatively more efficient model. Biogeochemical parameters in MOPS have been calibrated so that it accurately reproduces observed nutrient distributions and fluxes such as N₂ fixation and denitrification (Kriest and Oschlies, 2015). Because circulation patterns used in earlier calibration exercises of MOPS (Kriest et al., 2020) are similar to those of the physics-only FOCI model (Matthes et al., 2020), a similar performance of MOPS is expected here.

In this paper, we present the technical description of the marine biogeochemistry component in FOCI and its validation for the model mean state following a 500-year spin-up simulation, for historical simulations covering the period from 1850 to 2014, and control simulations with pre-industrial conditions. We also discuss the variability among ensemble members of each set-up and the differences between CO₂-concentration-driven and CO₂-emission-driven experiments.

2 Model description

2.1 Ocean circulation and the coupling to the atmosphere

The physical ocean model component in FOCI is detailed in Matthes et al. (2020). In brief, the ocean model is built on NEMO version 3.6 (Madec and the NEMO team, 2016) with a nominal global ocean resolution of $\mathrm{1}/\mathrm{2}{}^{\circ }$ on a tri-polar grid (ORCA05), with Louvain-la-Neuve sea Ice Model version 2 (LIM2) as the dynamic–thermodynamic sea ice model (Madec and the NEMO team, 2016). There are 46 vertical levels, with thicknesses varying from 6 m at the surface to 250 m in the deep ocean. A two-step flux-corrected transport, total variance dissipation scheme (TVD; Zalesak, 1979) is used for tracer advection to ensure positive definite values. Tracer diffusion is aligned along isopycnals, and viscosity is applied via a bi-Laplacian operator. The exchange of momentum, heat, freshwater fluxes, and sea ice properties between the ocean and the atmosphere is realized by the OASIS3-MCT coupler (Valcke, 2013). Note that all air–sea flux calculations are performed in the atmospheric module and have to be mapped from the coarser spatial grid of the atmospheric model (approximately 1.8^∘) to the finer one of the ocean ($\mathrm{1}/\mathrm{2}{}^{\circ }$).

2.2 Ocean biogeochemistry

MOPS (Model of Oceanic Pelagic Stoichiometry) simulates the elemental cycles of oceanic phosphorus, nitrogen, and oxygen and consists of seven compartments, namely phosphate, nitrate, oxygen, phytoplankton, zooplankton, detritus, and dissolved organic matter (DOM; Fig. 1). We here only provide a general overview of the model structure and the changes made for its implementation in FOCI. For further details of MOPS, we refer the reader to the original description of MOPS by Kriest and Oschlies (2015) and to the detailed model description in Appendix A.

Figure 1 Schematic of the ocean biogeochemistry model FOCI-MOPS. Arrows indicate processes and fluxes between model compartments (solid squares). The dashed lines between ALK and NO₃ and between ALK and PO₄ show that ALK is affected by changes in NO₃ and PO₄. CaCO₃ is represented as a dotted square as it is not a prognostic tracer in the model. See Appendix A for a detailed description. ALK stands for alkalinity, DIC stands for dissolved inorganic carbon, and DOM stands for dissolved organic matter.

For the implementation in FOCI, MOPS has been complemented with a carbon cycle that includes biological uptake and remineralization effects on dissolved inorganic carbon (DIC) and alkalinity (ALK; assuming fixed elemental ratios according to the stoichiometry by Paulmier et al., 2009) and the effects of formation and dissolution of calcite on these two tracers. For biogenic calcite production and dissolution, we implement an implicit approach (Schmittner et al., 2008) where the production of calcite is calculated from organic detritus production in a fixed ratio. Integrated over the entire water column and assuming a fixed molar CaCO₃:P ratio from the production of detritus at each time step, the newly produced vertically integrated calcite is then immediately distributed and released as DIC and ALK over the water column with an e-folding length scale. Air–sea gas exchange of CO₂ at the sea surface is calculated according to Orr et al. (2017). Altogether, the ocean biogeochemistry is simulated via nine prognostic tracers and five chemical elements (PO₄ (P), NO₃ (N), O₂ (O), DIC (C), ALK (C, N, P), phytoplankton (C, N, P), zooplankton (C, N, P), detritus (C, N, P, Ca), and DOM (C, N, P)).

In MOPS, phytoplankton growth depends on ambient PO₄, NO₃, temperature, and light. Phytoplankton is grazed by zooplankton and parameterized by a Holling-III function that uses a sigmoidal functional response of the grazing rate to increasing food (Holling and Buckingham, 1976). Subsequent zooplankton egestion and plankton mortality produce sinking detritus and neutrally buoyant DOM. The sinking speed of detritus increases linearly with depth, and the remineralization rate is constant and independent of temperature. In the absence of lateral or vertical exchange, these would result in a flux profile given by a power law of depth (the so-called “Martin” curve with exponent b; Martin et al., 1987). However, the model also simulates oxygen-dependent remineralization of organic matter. If oxygen falls below a threshold, denitrification sets in during anaerobic remineralization and reduces NO₃, albeit at a slower rate than aerobic remineralization. The decrease of remineralization caused by oxygen deficiency, especially in oxygen minimum zones, distorts the nominal value of b. There is no benthic denitrification in the model. The loss of fixed nitrogen due to denitrification affects the supply of NO₃ to the surface and influences nitrogen fixation at the sea surface, which in the model is diagnosed depending on temperature and the local ratio between simulated NO₃ and PO₄.

There is no sediment module in MOPS. Organic detritus arriving at the seafloor is partially buried and the burial fraction depends on the rain rate (Appendix A3). The resulting loss of phosphorus and nitrogen to the sediment is compensated for by adding an amount of PO₄ and NO₃ equivalent to the globally integrated burial, distributed homogeneously in the topmost model layer, rather than through river runoff. Likewise, we account for the burial loss of organic carbon and the associated virtual flux of ALK by a compensating supply of DIC and a decrease of ALK homogeneously distributed over the global sea surface. Supply and removal of these tracers to the surface layer ensures mass conservation with respect to the fluxes across the sea floor. Calcite arriving at the seafloor is not buried but immediately dissolved in the deepest model box, accommodating for the fact that a climate model like FOCI does not allow for the spin-up times needed for CaCO₃ sediments to reach a steady state.

With the implementation of MOPS, 1 year of FOCI-MOPS simulation takes 0.6 h and costs about 730 CPU hours with 1260 CPUs (Intel^® Xeon^® Platinum 9242 Processor) on 14 nodes on the North German Supercomputing Alliance (HLRN) complex LISE at the Zuse Institute Berlin (ZIB). For details of the hardware, please refer to the HLRN-IV documentation on HLRN web page (https://www.hlrn.de/, last access: 10 March 2022). With the same CPU configuration, the computing time for FOCI-MOPS increases by 26 % compared to a physics-only FOCI version. For comparison, an increase of 42 % computing cost was found for the biogeochemical component PISCESv2‐gas implemented in CNRM-ESM2-1 with an online grid-coarsening algorithm (Berthet et al., 2019).

A direct calibration by means of optimization of FOCI-MOPS, a computationally expensive Earth system model (ESM), is presently not feasible. Therefore, we selected MOPS parameters that resulted from a calibration using transport matrices that were derived from a circulation of a reanalysis dataset, the Estimating the Circulation and Climate of the Ocean (ECCO), which is computationally more efficient. In particular, six biogeochemical model parameters of MOPS were adjusted via an automatic calibration procedure against observed nutrient and oxygen distributions as described in Kriest et al. (2020, optimisation ECCO^*). While the circulation of the ECCO and FOCI are largely similar, to account for existing differences we manually adjusted three of the six parameters after initial tests as described in Appendix A5. The Appendix also describes the choice of parameters for regulating calcite formation and dissolution. The final full set of parameter values can be found in Table A1.

2.3 Model simulations and data used for model evaluation

Following the CMIP6 protocol (Eyring et al., 2016), we performed a series of experiments to evaluate FOCI-MOPS (Table 1). A 500-year spin-up with marine biogeochemistry (spinup) was restarted from the end of a 1500-year “physics-only” FOCI spin-up under year 1850 climate conditions (e.g. solar radiation, greenhouse gases, atmospheric nitrogen deposition, sulfate aerosol from volcanic eruptions, land usage, and population density; see Matthes et al., 2020, for details of the boundary conditions.) Details of physical characteristics of the 1500-year spin-up (FOCI1.3-SW038) are described in Matthes et al. (2020), including a cold bias in sea surface temperature (SST) and surface air temperature (SAT) in the North Atlantic and a warm bias in the Southern Ocean. The depth of the maximum transport of the Atlantic meridional overturning circulation (AMOC) is also shallower compared to the RAPID array observations (McCarthy et al., 2015). For the 500-year FOCI-MOPS spin-up, phosphate (PO₄), nitrate (NO₃), and oxygen (O₂) are initialized using the WOA2013 data set (Garcia et al., 2013a, b), and pre-industrial DIC and ALK are taken from GLODAPv2.2016b (Lauvset et al., 2016). The 480th, 490th, and 500th year of the spinup served as different initial conditions for an ensemble of three pre-industrial control (piControl) and transient historical (Hist) simulations (1850–2014) with prescribed atmospheric CO₂ concentrations. We also carried out a set of experiments where the model was forced with CO₂ emissions and the atmospheric CO₂ was calculated prognostically (instead of prescribing atmospheric CO₂ concentrations). For those experiments, we use the 480th year of the FOCI-MOPS spin-up as a starting point for a 250-year CO₂ zero-emission-driven spin-up (ESM-spinup) to allow for some equilibration between the atmosphere, land, and ocean carbon compartments. An ensemble of three ESM model pre-industrial control (ESM-piControl) and transient historical (ESM-Hist) simulations was then started from the 230th, 240th and 250th year of ESM-spinup.

Table 1 Overview of FOCI-MOPS simulations.

To evaluate the performance of FOCI-MOPS, we compared the distribution of inorganic tracers in Hist to interpolated and non-interpolated data of GLODAPv2.2016b (Lauvset et al., 2016; Olsen et al., 2016). For O₂, NO₃, PO₄, and ALK, model outputs are averaged over 1972 to 2013, and modelled year 2002 is used for DIC. For organic tracers, 10-year-mean (from 2005 to 2014) chlorophyll estimates from remote sensing (MODIS-Aqua, https://jeodpp.jrc.ec.europa.eu/ftp/public/JRC-OpenData/GMIS/satellite/9km/, last access: 20 January 2021, Melin, 2013), together with in situ observations of mesozooplankton, particulate organic nitrogen, and dissolved organic phosphorus, are used for comparison. Details of the biogeochemical data sets and model evaluation metrics are described in Appendix B.

3 Evaluation of model results

3.1 Temporal evolution of model simulations

3.1.1 Spin-up drift

Over the first 100 years of the FOCI-MOPS 500-year spin-up (spinup), a small (<1.5 %) decrease in inorganic nutrients (PO₄ and NO₃) is used to build up organic material (plankton biomass, DOM, and detritus), as they are initialized to zero at the start of spinup. After these initial adjustments, most of the global mean tracer concentrations and globally integrated fluxes showed drifts that were small relative to their mean concentrations (up to 0.1 % in the carbon flux at 2000 m) and reached or asymptotically approached a steady state at the end of spinup (Fig. 2). An exception is NO₃ because the marine nitrogen cycle did not reach a steady state, with N₂ fixation and denitrification not yet being in equilibrium after 500 years. This is commonly observed in models due to the spatial separation of these counteracting N-cycle processes. The equilibration of marine N₂ fixation and denitrification can take thousands of years (Falkowski, 1997; Oschlies et al., 2019). Owing to a net NO₃ loss due to a higher global denitrification than N₂ fixation rate, together with the build-up of organic matter in the beginning, global-average NO₃ is about 1.1 mmol N m⁻³ lower at the end of the spinup compared to the beginning, amounting to a decrease of 0.007 % per year. After a small positive spike in the beginning, O₂ concentrations continuously decrease throughout the spinup and in the end are close to a steady state. The loss of O₂ indicates that in the model the supply of O₂ via mixing is too slow and/or the consumption due to remineralization is overestimated. The changes in O₂ concentration are associated with the temporal drift of the carbon export across 2000 m. While the export flux at 100 m already reached a steady state after 100 years, the flux at 2000 m was still increasing at the end of the spin-up. This reflects that the remineralization of organic matter was slowing down in the upper 2000 m due to the decreasing O₂ concentration, allowing an increasing fraction of organic material to remineralize below 2000 m. After 500 years, global-average O₂ is about 8.4 mmol O₂ m⁻³ (5 %) lower than in the beginning. In the model, DIC varies mainly due to a build-up of organic matter and changes in air–sea CO₂ fluxes. Drift in the DIC inventory during the last 100 years in the spinup is −0.086 Pg C yr⁻¹, which meets the “acceptably small drift” (±0.1 Pg C yr⁻¹) suggested in Jones et al. (2016). Global CaCO₃ production equals dissolution as there is no dynamic CaCO₃ pool in the model. Therefore, the alkalinity inventory is only affected by the changes in PO₄ and NO₃. Since PO₄ and NO₃ decrease during the formation of organic matter and the nitrate reservoir additionally decreases resulting from the imbalance of denitrification and N₂ fixation (plus nitrification), the ALK inventory increases during the spinup in the stoichiometric ratio of $\mathrm{0.9914}:-\mathrm{1}$ ALK : NO₃ (Paulmier et al., 2009). The remaining small drifts in the spinup might exert an effect on the historical simulations, where they were accounted for by subtracting the respective control simulations.

Figure 2 Time series of tracer concentrations and fluxes in the 500-year spinup experiment of the inorganic tracers (PO₄, NO₃, O₂, DIC, and ALK), organic tracers (phytoplankton, zooplankton, detritus, and DOP), and fluxes (primary production, export production at 100 m, carbon flux at 2000 m (F2000), CaCo₃ production, denitrification (Denitr), and N₂ fixation (Nfix)).

3.1.2 Historical simulations

Even after the 500-year spin-up (see Sect. 3.1.1), some drifts in tracers and fluxes can still exist and can been seen in the piControl and ESM-piControl (Figs. S1–S2 in the Supplement). The drifts in Hist and ESM-Hist (Figs. S3–S4) runs are removed by subtracting the piControl and ESM-piControl simulation trends from the corresponding historical runs (Fig. 3). Drifts continue over the additional years of the CO₂-emission-driven spin-up (ESM-spinup); therefore, the absolute tracer concentrations and fluxes in ESM-Hist, which is initialized from the end of ESM-spinup, are partly different from those in the Hist runs, as is apparent in the time series of NO₃, O₂, DIC, and ALK discussed above. Nevertheless, the historical evolution and variability of the ESM-Hist simulation is qualitatively and quantitatively very similar to the one in Hist. In the Hist simulations, PO₄ increases slightly, this is mainly due to a decrease in dissolved organic phosphorus (DOP; the model also implicitly represents DOM in phosphorus units in a C : N : P molar ratio of $\mathrm{117}:\mathrm{16}:\mathrm{1}$), which declines by 0.007 % per year. Decreases in phytoplankton, zooplankton, and detritus are also small, ranging from 0.005 % to 0.01 % per year. Global-average NO₃ in 2014 is about 0.017 mmol N m⁻³ (0.06 %) higher than in 1850. One reason for the increase in NO₃ is a decrease of denitrification. The reduction of DOM and organic particle pools also contributes to the NO₃ increase (0.006 mmol N m⁻³). The changes in NO₃ are also reflected in a decreasing ALK inventory. The simulated relative decrease of the O₂ inventory is higher than that of PO₄ and NO₃. In the model, the ocean loses about 1.5 mmol O₂ m⁻³ (0.9 %) of oxygen between 1850 and 2014. Changes in circulation and mixing, together with the decrease in solubility due to a warming ocean, dominate the decline in marine O₂ content. The small decrease in export production would lead to reduced O₂ consumption via respiration and therefore elevate rather than reduce the oxygen inventory. Globally averaged DIC increases by 8 mmol C m⁻³ (0.35 %), mainly due to increased air–sea CO₂ fluxes into the ocean under rising atmospheric CO₂, i.e. the uptake of anthropogenic CO₂.

Figure 3 Time series of tracer concentrations and fluxes for the 165 years (1850–2014) Hist and ESM-Hist experiments, each for a three-member ensemble. Each ensemble member is corrected by the drifts in the respective piControl and ESM-piControl runs. Dark blue and dark red lines stand for the mean in Hist and ESM-Hist, respectively, and light blue and light red shaded areas represent 1 standard deviation in Hist and ESM-Hist, respectively. The arrangement of the panels is the same as in Fig. 2

3.2 Biogeochemical model performance

3.2.1 Spatial distribution of inorganic tracers

The surface concentration of O₂ is primarily determined by the temperature-dependent solubility. In the ocean interior, it is a sensitive balance between ocean circulation and oxygen consumption during remineralization of detritus and dissolved organic matter (DOM). In FOCI-MOPS, water column denitrification occurs in low-oxygen waters when O₂ concentration is below 36 mmol O₂ m⁻³. The low-oxygen waters develop where the physical supply of O₂ is sluggish and O₂ consumption via respiration is high. The model reproduces the low-oxygen waters at the eastern margins of tropical and subtropical ocean basins (Figs. 4a and S5). In the zonally averaged ocean basin means, the low-oxygen waters are situated between 100 and 1500 m in the Atlantic, Indian, and Pacific oceans (Fig. 4a and d–f). In the same depth range, the modelled O₂ is biased low, particularly in the Southern Hemisphere (Figs. 4i–k and S5). Another region with clear negative O₂ biases can be seen at depths below 1500 m and all the way to the bottom in the Pacific Ocean (Fig. 4k). The low biases in the zonal mean are due to an underestimated O₂ in the eastern equatorial Pacific (Fig. S5), a feature that is likely due to overestimated production in the euphotic zone above and imperfect physics and remineralization settings, which are commonly found in numerical models (Dietze and Loeptien, 2013; Ilyina et al., 2013; Cabré et al., 2015; Paulsen et al., 2018). In general, O₂ in the model is biased high between 1000 and 3000 m and is biased low below 3000 m. In addition to the biological processes, the difference can be explained by the ventilation of the water masses. The modelled AMOC has a shallow bias and is indeed weaker in the deep water (3000–5000 m; Fig. 5), as also remarked upon in Matthes et al. (2020) and common across climate models (e.g. Weijer et al., 2020). A more sluggishly ventilated deep water in latitude–depth structure (Fig. 5a) is consistent with the higher concentration in inorganic tracers and the negative biases in O₂.

Figure 4 O₂ climatology over 1972–2013 of the Hist simulations at (a) 300 m depth and zonally averaged across the (c) Southern, (d) Atlantic, (e) Indian, (f) Pacific, and (g) Arctic oceans in FOCI-MOPS. The corresponding difference of model minus observation are shown in (b) and (h)–(l), respectively. Observational data are from GLODAPv2.2016b, which covers 1972–2013 (Lauvset et al., 2016; Olsen et al., 2016). Red contour lines in (a) depict O₂ of 36 mmol O₂ m⁻³, below this concentration anaerobic remineralization kicks in in the model. Partitioning of ocean basins follows the World Ocean Atlas 2013 definitions.

Figure 5 Atlantic meridional overturning stream function climatology over 2004–2014 of the Hist simulations. (a) Latitude–depth structure, (b) vertical profile at 26.5^∘ N, and (c) z derivative of AMOC at 26.5^∘ N. Blue lines and the shaded areas are modelled results, and dashed black lines are observational values from the RAPID array (McCarthy et al., 2015).

Surface phosphate (PO₄) concentrations simulated by FOCI-MOPS generally agree with observations, with positive biases in the Atlantic and a negative bias in the subarctic gyre in the northern Pacific Ocean (Fig. 6a and b). In the interior, the model–data misfits are generally smallest in the Southern Ocean (SO) and the Arctic (Fig. 6h and l). In the Atlantic, Indian, and Pacific oceans, modelled PO₄ shows positive biases at around 100–1000 m and below 3000 m and negative biases between 1000 and 3000 m (Fig. 6i–k).

Figure 6 PO₄ climatology over 1972–2013 of the Hist simulations (a, c–g) and the model–observation difference (b, h–i). Observational data are from GLODAPv2.2016b, which covers 1972–2013 (Lauvset et al., 2016; Olsen et al., 2016). The order of the panels is the same as in Fig. 4.

In order to better understand the distribution of PO₄, we adopt the methodology of Duteil et al. (2012) and calculate regenerated (PO${{}_{\mathrm{4}}}_{\mathrm{reg}}$) and preformed PO₄ (PO${{}_{\mathrm{4}}}_{\mathrm{pre}}$) assuming oxygen saturation (O${{}_{\mathrm{2}}}_{\mathrm{sat}}$) at isopycnal outcrops (e.g. PO${{}_{\mathrm{4}}}_{\mathrm{reg}}=$ AOU$/r-$O₂: PO₄, AOU = O${{}_{\mathrm{2}}}_{\mathrm{sat}}-$ O${{}_{\mathrm{2}}}_{\mathrm{obs}}$, PO${{}_{\mathrm{4}}}_{\mathrm{pre}}=$ PO₄−PO${{}_{\mathrm{4}}}_{\mathrm{reg}}$). For observations, we apply the O₂: PO₄ ratio of 170 mol O₂: mol P estimated in Anderson and Sarmiento (1994). For model output, we apply the ratio of 165.08044 mol O₂: mol P used in the model.

Except for the deep Atlantic (deeper than 3000 m), preformed PO₄ is generally too low in the ocean interior compared to observations (Fig. 7h–k). A lack of iron limitation in the model could cause a too strong uptake of PO₄ by phytoplankton in iron-limited regions, such as the Southern Ocean, and lead to too low preformed PO₄. The overestimated regenerated PO₄ in the deep Pacific >3000 m (Fig. 8k) contributes to the overall high bias of PO₄ (Fig. 6k) and is consistent with the low O₂ biases in the same region (Figs. 4k and S5). The high bias of PO₄ may be at least partly explained by excessive remineralization of detritus and DOM in the eastern Pacific Ocean as a result of overestimated primary production (PP) and export production (EP; see Sect. 3.3). A too sluggish ventilation of these water masses could also explain the excessive accumulation of regenerated PO₄ and apparent oxygen utilization (Fig. 5), as the PO₄ is also overestimated in the deep Atlantic Ocean.

Figure 7 Preformed PO₄ climatology over 1972–2013 of the Hist simulations (a, c–g) and the difference with observations using GLODAPv2.2016b, which covers 1972–2013 (Lauvset et al., 2016; Olsen et al., 2016) (b, h–i). Preformed PO₄ is estimated as PO₄ minus regenerated PO₄. The order of the panels is the same as in Fig. 4.

Figure 8 Regenerated PO₄ climatology over 1972–2013 of the Hist simulations (a, c–g) and the difference with observations using GLODAPv2.2016b, which covers 1972–2013 (Lauvset et al., 2016; Olsen et al., 2016) (b, h–i). Regenerated PO₄ is estimated as apparent oxygen utilization; see text. The order of the panels is the same as in Fig. 4.

Biases of climatological nitrate NO₃ overall show a similar distribution as for phosphate (Fig. 9a). Simulated NO₃ differs from the pattern of PO₄ mostly in low latitudes at intermediate depth. Some positive biases can be seen in the Arctic and near the Southern Hemispheric Subtropical Front (Fig. 9b). In the ocean interior, NO₃ is removed by denitrification. Substantial rates of denitrification can be realized (depending on the supply of organic matter) when O₂ concentrations fall below a few mmol m⁻³ (Eqs. A11 and A16). Modelled denitrification occurs in the Atlantic Ocean, the northern Indian Ocean, and the eastern tropical Pacific Ocean and is most prominent in the Pacific Ocean (black contour lines in Fig. 9i–k). Globally integrated water column denitrification in the model amounts to 0.13 Pg N yr⁻¹, higher than observation-based estimates (ranging from 0.02 to 0.12 Pg N yr⁻¹; see Sect. 3.3 below). It is likely that excessive denitrification contributes to some of the negative biases between 300 and 1000 m in the low-latitude oceans (Figs. 9i–k and S5) that are not accompanied by similar negative biases in PO₄ (Fig. 6i–k).

Figure 9 NO₃ climatology over 1972–2013 of the Hist simulations and the difference model-observation. Black contour lines represent zonally integrated denitrification rates (mmol N m⁻² d⁻¹). Observational data are from GLODAPv2.2016b, which covers 1972–2013 (Lauvset et al., 2016; Olsen et al., 2016). Panels are ordered as in Fig. 4.

The spatial pattern of dissolved inorganic carbon (DIC) and alkalinity (ALK) in general agree with observed patterns at the sea surface, but both show negative biases (Figs. 10b and 11b). In the interior, the biases turn positive below 3000 m for both DIC and ALK, which is most obvious in the Atlantic and the Pacific Oceans, similar to the patterns of PO₄ and NO₃. In the ocean, potential differences in CO₂ uptake between model and observations can contribute to DIC biases. Both DIC and ALK are affected by the remineralization of detritus and the production of CaCO₃, via zooplankton egestion and plankton mortality in the surface layer and dissolution of CaCO₃ throughout the water column, respectively. There is a clear transition from a negative bias of DIC and ALK in the surface ocean to a positive bias in deep waters across all ocean basins, except the Arctic where plankton activity is limited (Fig. 11h–l). This indicates that the e-folding export and dissolution of CaCO₃ (Sect. A2.2) might lead to an excessive transport of ALK (and DIC) to the deep ocean. The positive biases in DIC and ALK in the deep water might also indicate a too sluggish ventilation of deep waters, which is also consistent with the positive biases in PO₄ and NO₃, as well as the negative bias in O₂.

Figure 10 Dissolved inorganic carbon (DIC) for the year 2002 in the Hist simulations and the model–observation difference. Observational data are sum of the preindustrial TCO₂ and the accumulated anthropogenic carbon in the year 2002 from GLODAPv2.2016b (Lauvset et al., 2016; Olsen et al., 2016). The order of panels is the same as in Fig. 4.

Figure 11 Alkalinity (ALK) climatology over 1972–2013 in the Hist simulations and the model–observations difference. Observational data are from GLODAPv2.2016b, which covers 1972–2013 (Lauvset et al., 2016; Olsen et al., 2016). The order of panels is the same as in Fig. 4.

Typically, 1 standard deviation of distribution of the inorganic tracers is smaller than 10 % of the mean value of the Hist ensemble simulations (Figs. S6–S12). Areas with higher variation usually occur where vertical mixing is stronger, such as at high latitudes.

3.2.2 Spatial distribution of organic tracers in the surface

In FOCI-MOPS, the growth of phytoplankton is determined by temperature, light, and ambient NO₃ and PO₄ concentrations (Eqs. A1–A5). The distribution of phytoplankton (Phy) in the surface layer (∼6 m) overall agrees relatively well with observational estimates (Fig. 12a, e, and i), the most distinct discrepancies occur in coastal regions, where the observed phytoplankton concentrations often exceed 0.06 mmol P m⁻³ but modelled phytoplankton do not (Fig. 12e). The difference might be due to the lack of terrestrial nutrient runoff in the current model version. In addition, shelves are not well resolved with a $\mathrm{1}/\mathrm{2}{}^{\circ }$ spatial resolution, a typical bias of global models. In the open-ocean pelagic regions, however, the modelled phytoplankton is usually higher than the observations, which might be due to the missing iron limitation in the model. The higher phytoplankton biomass might also explain some of the low biases in PO₄ in the equatorial regions.

Figure 12 Organic tracer concentrations for the first layer (∼6 m) of phytoplankton (Phy), 0–100 m of zooplankton (Zoo), particulate organic phosphorus (POP), and dissolved organic phosphorus (DOP) climatology over 2005–2014 in the Hist simulations (a–d), observations (e–h), and the zonally averaged values (i–l). In the zonally averaged panels, model results are subsampled according to the availability of observations. Observational data of Phy, Zoo, and POP are derived from chlorophyll (MODIS-Aqua; Melin, 2013), mesozooplankton (MAREDAT; Moriarty and O'Brien, 2013), and POM (Martiny et al., 2014), respectively. Observations of DOP are from multiple published (Torres-Valdés et al., 2009; Moutin et al., 2008; Landolfi et al., 2016; Yoshimura et al., 2007) and unpublished (Landolfi, unpublished) data sets. See Sects. B2.2, B2.3, B2.4, and B2.5 for details of the observations.

Zooplankton and particulate organic phosphorus (POP, sum of simulated phosphorus in detritus, phytoplankton, plus half of the zooplankton; see Sect. B2.4 for further details) are generally closely associated with the presence of phytoplankton; their distributions largely correlate with those of phytoplankton in the model (Fig. 12b and c). While the large-scale pattern is similar, Zooplankton and POP are less homogeneously distributed than phytoplankton. For example, modelled phytoplankton biomass is around 3-fold greater in the nutrient-replete eastern tropical Pacific than in the oligotrophic gyres, while the same ratio can be over 10-fold for zooplankton. Modelled zooplankton are overestimated around the Equator and the Southern Ocean around 40^∘ S in the surface (0–100 m; Fig. 12j), where phytoplankton is also overestimated (Fig. 12i). Observations show much higher POP at high latitudes than low latitudes at the surface (0–100 m; Fig. 12k), which is not captured by the model. Note that observational data of phytoplankton, zooplankton, and POP are from different types of estimates; therefore, the observations do not necessarily correlate to each other as in the model. As dissolved organic phosphorus (DOP) is neutrally buoyant, it is controlled by ocean circulation in addition to the production and remineralization. DOP is much lower in the deeper water than at the surface; therefore, in upwelling regions such as the eastern equatorial Pacific, DOP is lower. Meanwhile, nutrient supply from the upwelling promotes the growth of phytoplankton and the production of POP. As a result, the spatial distribution of surface DOP in the model in general is negatively correlated with POP and is more homogeneously distributed (Fig. 12d). Observations of DOP are very sparse but suggest that DOP tends to be overestimated in the model. This may indicate a too long DOP remineralization timescale; however, a shorter timescale may lead to an even stronger bias in surface phosphate concentrations (Fig. 12). Alternatively, changing the partitioning between DOP and POP production (e.g., via parameter σ_DOP in Eqs. A9 and A10), together with modified POP sinking speed and remineralization could potentially improve DOP without negative effects on other (surface) components. Further investigations regarding the simultaneous fit of MOPS to all inorganic and organic tracers are currently underway.

Similar to the inorganic tracers, 1 standard deviation of distribution of the organic tracers is smaller than 10 % of the mean value of the Hist ensemble simulations (Fig. S13).

3.2.3 Statistical performance of the simulated inorganic and organic tracers

The correlation coefficient r, standard deviation σ_M, and centred root-mean-squared error (RMSE^′, RMSE minus global bias) of modelled nutrients and oxygen show a very good match compared to the non-interpolated observations (Table 2, Fig. 13), although nutrients are biased somewhat low, especially at the surface (0–100 m), which is the only depth range where oxygen shows a high bias (note the shading of symbols in Fig. 13). The modelled ALK and DIC are biased slightly low above 2000 m. The modelled surface ALK averaged over 1972 to 2013 is 2286±1 mmol m⁻³, and surface DIC averaged over 1870 to 1899 (pre-industrial DIC) is 1985±0 mmol m⁻³. They are lower than the values among CMIP5 models (surface ALK and pre-industrial DIC are 2365–2500 and 2050–2170 mmol m⁻³, respectively; Oka, 2020) but are only slightly lower than the observations (surface ALK and pre-industrial DIC are 2362 and 2019 mmol m⁻³, respectively; Lauvset et al., 2016; Olsen et al., 2016). Note that the low bias in the full-domain ALK here is contrary to a higher inventory when compared to the initial state (Figs. 2 and 3). The low bias against individual (non-interpolated) data results from the fact that more observations are located in the upper 2000 m, where the model underestimates ALK concentrations. Compared to global model assessments of CMIP5 models with similar metrics evaluated by Ilyina et al. (2013), Séférian et al. (2013), and Kwiatkowski et al. (2014), the model performs well (Table 2) with respect to dissolved inorganic tracers.

Table 2 Correlation coefficient r, normalized standard deviation ${\mathit{\sigma }}_{\mathrm{M}}/{\mathit{\sigma }}_{\mathrm{O}}$, global bias normalized by observed mean (in squared brackets the non-normalized bias is shown, given in µg Chl L⁻¹ for phytoplankton and in mmol m⁻³ for all other tracers), and Bhattacharyya distance (BD) of experiment Hist presented in this study and from studies by Ilyina et al. (2013, historical run in “MR” model configuration, Table 5; “I2013”), Séférian et al. (2013, one biogeochemical model with three different circulations, Table 3; “S2013”) and Kwiatkowski et al. (2014, six biogeochemical models in one circulation, Table 3; “K2014”). Except for chlorophyll, the values in FOCI-MOPS are from comparison with non-interpolated data from GLODAPv2.2016b (Lauvset et al., 2016; Olsen et al., 2016). Round brackets are introduced for distinguishing hyphen and negative signs. For all model studies we report metrics for the surface (here refers to 0–100 m); additionally, for the metrics by Ilyina et al. (2013) we show the range over five different depth levels (surface, 100, 500, 1000, and 3000 m) in the “all domains” column. For the all domains column in this study we show the range over seven different vertical domains (0–100, 100–200, 200–500, 500–1000, 1000–2000, 2000–5000, 0–8000 m). For metrics by Séférian et al. (2013), we report the range over three different circulations and for metrics by Kwiatkowski et al. (2014) the range over five different models is reported. Note that, except for the Bhattacharyya distance, all metrics of Hist have been calculated from volume-weighted properties.

^* Computed from phosphorus units via a C : Chl according to Sathyendranath et al. (2009).

Figure 13 Taylor diagrams for phosphate (PO₄), nitrate (NO₃), oxygen (O₂), dissolved inorganic carbon (DIC), and alkalinity (ALK) (left to right). Colours indicate model bias (relative to observations; in percent). The model is subsampled according to the availability of observations. Observational data are from GLODAPv2.2016b (Lauvset et al., 2016; Olsen et al., 2016). Symbols indicate different vertical domains: star indicates 0–100 m, square indicates 100–200 m, diamond indicates 200–500 m, triangle indicates 500–1000 m, delta indicates 1000–2000 m, small delta indicates 2000–5000 m, and circle indicates full domain.

The model fit worsens for the organic tracers, with a correlation coefficient between 0.11 (DOP) and 0.42 (POP), and generally a too low spatial variance (between 41 and 61 % of the observations; Table B2 and Fig. 14). Except for POP, organic tracers are biased high, and RMSE and RMSE^′ are large. Thus, the statistical metrics indicate a lower performance of the model for the organic tracers chlorophyll, zooplankton, DOP, and POP. Nevertheless, the results for phytoplankton and chlorophyll are in line with earlier global model studies (Table 2). Note also that the uncertainty of the observational estimates of organic tracers is larger compared to inorganic tracers, the in situ observations are sparse in space and time, with sampling biases, e.g., towards summer and with less high-latitude observations (see also Appendix B2.).

Figure 14 Taylor diagrams for phytoplankton (PHY), zooplankton (ZOO), particulate organic phosphorus (POP), and dissolved organic phosphorus (DOP) as shown in Fig. 12 (left to right). Colours indicate model bias (relative to observations; in percent). Observational data of Phy, Zoo, and POP are derived from chlorophyll (MODIS-Aqua; Melin, 2013), mesozooplankton (MAREDAT; Moriarty and O'Brien, 2013), and POM (Martiny et al., 2014), respectively. Observations of DOP are from multiple published (Torres-Valdés et al., 2009; Moutin et al., 2008; Landolfi et al., 2016; Yoshimura et al., 2007) and unpublished (Landolfi, unpublished) data sets. See Sects. B2.2, B2.3, B2.4, and B2.5 for details of the observations. The model is subsampled according to the availability of observations. Symbols indicate different vertical domains: triangle indicates surface (∼6 m), star indicates 0–100 m, and circle indicates full domain (POP only).

Even if a biogeochemical model component is dynamically correct, a slight distortion in the physical model (e.g. a slight spatial shift in ocean current) can cause a large RMSE and thereby induce a large model error. To account for such effects we have calculated three metrics that do not rely on the exact spatial pattern of tracers, but on the concentration distribution, namely the Bhattacharyya distance (BD, Bhattacharyya, 1946), which evaluates the similarity between observed and simulated frequency distributions of tracers, the Hellinger distance (HD, Hellinger, 1909), which is related to BD via BD $=-\ln \left(\mathrm{1}-{\mathrm{HD}}^{\mathrm{2}}\right)$, and the L1 norm, which evaluates the absolute difference between observed and simulated distributions. For all three metrics, smaller model–data misfit is associated with larger area of overlap between the two distributions, which yields smaller values (see Appendix B3). Figures 15 and 16 show examples for the distributions of inorganic and organic tracers, and Table B2 lists the metrics for different tracers. Applying these metrics, it is evident that most organic tracers (phytoplankton, zooplankton, and POP), despite suffering from a large error with respect to r and RMSE, in general reflect the observed distribution (Fig. 16), and they thus exhibit values for, e.g., BD, that are in the same range as those of inorganic tracers. Only the strong overestimate of DOP by the model is also reflected in BD. Here the model shows a much higher misfit than for any other tracers for all metrics.

Figure 15 Frequency distribution of phosphate (PO₄), nitrate (NO₃), oxygen (O₂), DIC, and alkalinity (left to right) from observations (filled black bars) and model (red bars) for the surface (0100 m, top) and the full model domain (bottom). Numbers denote three different metrics for the similarity of the distributions, namely L1 (Eq. B3), HD (Eq. B2), and BD (Eq. B1). Observational data are from GLODAPv2.2016b. The data cover 1972–2013. except for DIC, which represents for the year 2002 (Lauvset et al., 2016; Olsen et al., 2016).

Figure 16 Frequency distribution of surface phytoplankton, zooplankton (0–100 m), POP (full domain), and DOP (0–100 m). Observational data of Phy, Zoo, and POP are derived from chlorophyll (MODIS-Aqua; Melin, 2013), mesozooplankton (MAREDAT; Moriarty and O'Brien, 2013), and POM (Martiny et al., 2014), respectively. Observations of DOP are from multiple published (Torres-Valdés et al., 2009; Moutin et al., 2008; Landolfi et al., 2016; Yoshimura et al., 2007) and unpublished (Landolfi, unpublished) data sets. See Sects. B2.2, B2.3, B2.4, and B2.5 for details of the observations. The model is subsampled according to the availability of observations. Colours and numbers are as described in Fig. 15.

To summarize, in terms of tracer distributions the performance of MOPS in FOCI is comparable to other global models and even generally better with respect to inorganic tracers. The fit deteriorates with regard to the organic components in the cases of plankton and POP and has a high bias for DOP.

3.3 Globally integrated biogeochemical fluxes and their spatial patterns

3.3.1 Oceanic biogeochemical fluxes

We average the globally integrated biogeochemical fluxes simulated by FOCI-MOPS over the time period 2005 to 2014 and over the ensemble of three Hist simulations (Table 3). Overall, primary production (37.8 Pg C yr⁻¹), export production, the sinking of detritus at 100 m (6.8 Pg C yr⁻¹), the flux of detritus at 2000 m (0.35 Pg C yr⁻¹), burial of detritus that sank to the bottom of the ocean (0.34 Pg C yr⁻¹), and N₂ fixation (0.12 Pg N yr⁻¹) are within the range of observation-based estimates. CaCO₃ production (0.77 Pg C yr⁻¹) is lower than previous observation-based estimates (0.8–4.7 Pg C yr⁻¹). Water column denitrification (0.13 Pg N yr⁻¹) is slightly higher than observational estimates.

Carr et al. (2006)Honjo et al. (2008)Buitenhuis et al. (2013)Dunne et al. (2007)Lutz et al. (2007)Honjo et al. (2008)Siegel et al. (2014)Dunne et al. (2007)Honjo et al. (2008)Guidi et al. (2015)Muller-Karger et al. (2005)Wallmann (2010)Iglesias-Rodriguez et al. (2002)Balch et al. (2007)Lee (2001)Buitenhuis et al. (2019)Eugster and Gruber (2012)Luo et al. (2012)Somes et al. (2013)Wang et al. (2019)Bianchi et al. (2012)Eugster and Gruber (2012)DeVries et al. (2012)DeVries et al. (2013)Somes et al. (2013)Wang et al. (2019)

Table 3 Primary production, export production, flux of detritus at 2000 m, global burial, CaCO₃ production (Pg C yr⁻¹), N₂ fixation, and water column denitrification (Pg N, yr⁻¹) in Hist simulations (averaged over 2005–2014) and observation-based estimates.

The highest annual primary production occurs in the eastern equatorial Pacific Ocean and some coastal regions (>1000 mg C m⁻² d⁻¹, Fig. 17a). In the subtropical oligotrophic gyres, it is typically lower than 200 mg C m⁻² d⁻¹. The spatial distributions of export production at 100 m and CaCO₃ production (Fig. 17b and c) are similar to that of primary production, suggesting that the latter is largely driving the former. The spatial pattern of the flux of sinking detritus at 2000 m is less similar to that of primary production due to the effects of advection and mixing (Fig. 17d). The distinct high flux at 2000 m depth in the eastern Pacific is partly due to the low O₂ levels between the surface and 2000 m depth mentioned earlier. Regions of elevated flux of detritus largely overlap with O₂ concentrations below 36 mmol O₂ m⁻³, which is the threshold for anaerobic remineralization in the model. Denitrification is shaped by low O₂ levels associated with high primary production, consequent detritus flux, and sluggish ventilation, particularly in the so-called shadow zones along the eastern margins of the subtropical oceans (Fig. 17e). These biological and physical effects together result in the most pronounced water column denitrification in the eastern tropical Pacific, which accounts for 71 % of the global denitrification rate in our model, which is within the range that was derived from nitrogen gas measurements (70 %–88 %; DeVries et al., 2012). Due to temperature limitation, N₂ fixation occurs mostly between 40^∘ N and 40^∘ S in the model. In the model, 60 % of N₂ fixation occurs in the Pacific Ocean, in agreement with (though at the lower end of) observational estimates, which range between 60 % and 80 % (Luo et al., 2012). Variation (1 standard deviation) among the three Hist simulations of the biogeochemical fluxes are smaller than 10 % of the mean fluxes (Fig. S14).

Figure 17 Patterns of primary production, export production at 100 m, CaCO₃ production, flux of detritus at 2000 m, denitrification, and N₂ fixation in the Hist simulations averaged over 2005–2014. The red contour line in (d) indicates an O₂ level of 36 mmol O₂ m⁻³. Values are column-integrated except for the export production and the flux of detritus.

3.3.2 Air–sea exchange of carbon

Figure 18 shows the ocean pCO₂ and the air–sea CO₂ fluxes (positive downwards) in the Hist simulations for the 2005–2014 period, as well as the observation-based estimate for the same period (Landschützer et al. , 2017). In general, sea surface pCO₂ and air–sea CO₂ flux are affected by various biological and physical processes (Dong et al., 2016; Qu et al., 2022). For example, a too high phytoplankton production could result in a too low DIC and lead to an underestimate of the magnitude of negative CO₂ flux (outgassing). Our model results are consistent with such a scenario in the eastern equatorial Pacific, northern Pacific, and parts of the Southern Ocean. The too high phytoplankton production and resultant underestimated DIC could be due to the lack of Fe limitation, as we can see surface PO₄ and NO₃ are underestimated in the same regions. On the other hand, physical factors such as sea surface temperature (SST; Fig. 19; see Fig. S16 for the standard deviation among the Hist simulations) can also affect pCO₂ and CO₂ flux. Overestimated SST can be associated with positive pCO₂ and negative air–sea CO₂ flux biases. Those biological and physical factors result in complex patterns in the pCO₂ and CO₂ flux, and it is often difficult to quantify the contribution of individual factors. On top of these factors, our model–data comparison is for the period of present-day data (2005–2014). Therefore, the mismatches in pCO₂ and CO₂ uptake can be from two components: the mismatch inherited from the spin-up at a steady state and the mismatch in the uptake of anthropogenic carbon accumulated during the historical period. Since there are no pre-industrial observations, we cannot differentiate the contributions from each part. Nevertheless, the modelled large-scale pattern of CO₂ fluxes and uptake agrees with observations, and the biases in air–sea CO₂ fluxes are in the range of other Earth system models participating in the CMIP5 intercomparison (Dong et al., 2016). When zonally averaged, much of the biases are averaged out, and both ocean pCO₂ and air–sea CO₂ fluxes in FOCI-MOPS match the observations rather well (Fig. 18d and h).

Figure 18 Ocean pCO₂ and air–sea CO₂ fluxes (positive downwards). (a, e) Hist ensemble mean averaged over 2005–2014. Panels (b) and (f) are observations (Landschützer et al. , 2017) over the same time period. Panels (c) and (g) are model biases with respect to the observations, and (d) and (h) are zonally averaged values of observations and model results where observations are available.

Figure 19 Sea surface temperature (SST). (a) Hist ensemble mean averaged over 2005–2014. (b) The difference between model results and observations (HadISST, Rayner et al., 2003).

At high latitudes, the variation in the ocean pCO₂ and the air–sea CO₂ fluxes are usually higher (Fig. S15) due to the physical mixing being stronger and having higher uncertainty.

3.4 Response to increasing atmospheric CO₂

3.4.1 Atmospheric CO₂, air–sea CO₂ flux, surface air temperature anomaly, and ocean heat content anomaly

Atmospheric CO₂ in ESM-piControl drifts slowly from 286.0±0.4 ppm between 1850 and 1879 to 287.1±0.7 ppm between 2005 and 2014, with an average of 286.4±0.14 ppm over the whole period, which is 2 ppm above the prescribed concentration in piControl. During the historical period, atmospheric CO₂ in ESM-Hist increases from 286.1±0.46 ppm in 1850 to 410.1±2.33 ppm in 2014, slightly higher than the historical value of 397.5 ppm (Fig. 20a). According to these results, the pre-industrial ocean in FOCI-MOPS acts as a small net CO₂ source to the atmosphere. The mean global oceanic CO₂ loss in piControl and ESM-piControl during the simulation period are 0.07±0.01 and 0.05±0.01 Pg C yr⁻¹, respectively. In the historical simulations, the globally integrated air–sea CO₂ flux gradually increases between 1850 and 1960 along with increasing anthropogenic CO₂ perturbation. After 1960, the rate of increase accelerates, and the flux reaches 2.11±0.10 and 2.21±0.14 Pg C yr⁻¹ over 2005–2014 period in Hist and ESM-Hist, respectively, both are slightly lower than the 2.34 Pg C yr⁻¹ estimated by the Global Carbon Project (Friedlingstein et al., 2020; see the summary tab; Fig. 20b).

Figure 20 Time series of global atmospheric CO₂, globally integrated air–sea CO₂ fluxes (positive downwards), globally averaged surface temperature, and heat content at 0–2000 m in the ocean of Hist and ESM-Hist. The piControl runs are subtracted to account for drifts. For surface temperature anomalies, the mean values over 1850–1879 are subtracted from the historical period, and for ocean heat content anomalies, the mean values over 1955–1964 are subtracted. (a) Atmospheric CO₂ and (b) air–sea CO₂ fluxes, (c) surface temperature anomaly, and (d) ocean heat anomaly (0–2000 m) are all given. Grey (blue) and brown (red) lines stand for mean values of concentration-driven and emission-driven pre-industrial control (historical) simulations, respectively. The purple line is air–sea CO₂ flux data from the Global Carbon Project (Friedlingstein et al., 2020, see summary tab in the data sheet). The light green line is surface temperature data from HadCRUT5 (Morice et al., 2021). Dark green and brown lines are ocean heat content data from Ishii et al. (2017) and Cheng et al. (2017).

The cumulative air–sea CO₂ flux (positive downwards) from 1850 to 2014 in ESM-Hist amounts to 139±2 Pg C, which is 10 Pg C higher than that of Hist, with most of the difference building up over the time period 1850 to 1994 (Table 4). In the land model JSBACH, the air–land carbon flux is fulfilled by the NPP (net primary productivity), which is determined by the difference between photosynthesis and autotrophic respiration of vegetation. The vegetation carbon is lost to grazing by herbivores, litter production, and crop harvests and is consequently transported back to the atmosphere. We have followed Liddicoat et al. (2021) and calculated the compatible emission in the historical simulations. The cumulative compatible emission from 1850 to 2014 in the Hist simulations is 367±3 Pg C, lower than the lower limit of observation (380 Pg C), but agrees with the mean of CMIP6 models when the two highest model values are excluded (367 Pg C, Liddicoat et al., 2021). The underestimated cumulative compatible emission implies missing processes including land carbon uptake during the Second World War, during that period land use might not be correctly accounted for in the model forcing (Bastos et al., 2016). Despite the different cumulative fluxes until 1994, the cumulative CO₂ flux is about the same (25±1 Pg C) over 1994 to 2007 in the Hist and ESM-Hist runs (Table 4). Both Hist and ESM-Hist simulate lower cumulative CO₂ fluxes when compared with observations and are comparable to or slightly lower than CMIP6 models (Table 4). The cumulative fluxes over 1994–2007 are within the range of observations. The latitudinal distributions of cumulative CO₂ fluxes and carbon storage in FOCI-MOPS are consistent with those of CMIP6 models (Fig. 21). The highest CO₂ uptake between 40 and 65^∘ S is associated with the wind-driven upwelling, which brings deep water with low anthropogenic CO₂ to the surface and is able to uptake a higher amount of CO₂ (Frölicher et al., 2015). The stronger upwelling seems to be related to the higher variation in CO₂ uptake around the similar latitudes within the ensemble and among the CMIP6 models. The anthropogenic CO₂ is then transported to the mid-latitudes around 25–40^∘ S. A similar pattern occurs in the Northern Hemisphere. The considerable difference between Hist and ESM-Hist simulations around 45–65^∘ S is because of the underestimated cumulative compatible emission in the Hist. Nevertheless, the cumulative air–sea CO₂ flux is about 10 Pg C higher in the ESM-Hist and is less than 0.03 % of the total DIC inventory, it has an insignificant contribution to the interior DIC distribution between ESM-Hist and Hist (Fig. S17).

(1800–1994; Gruber et al., 2019)(Gruber et al., 2019)(CMIP6 model results; Terhaar et al., 2021)

Table 4 Cumulative air–sea CO₂ fluxes (Pg C) in FOCI-MOPS and observation-based estimates and CIMP6 model results. FOCI-MOPS results are from mean of Hist and ESM-Hist runs. Changes in piControl runs are subtracted from the historical runs.

^a Including only half of the air–sea CO₂ flux in 1994 and in 2007 to be consistent with the data period (mid-1994 to mid-2007)
provided in Gruber et al. (2019).
^b Mean of some CMIP6 models, including ACCESS-ESM1-5, CanESM5-CanOE, CanESM5, CESM2, CESM2-WACCM, CNRM-ESM2-1, GFDL-CM4, GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, NorESM2-LM, UKESM1-0-LL, and MIROC-ES2L (Terhaar et al., 2021, unconstrained results).

Figure 21 Cumulative air–sea CO₂ flux (positive downwards) and carbon storage between 1850 and 2014 simulated by FOCI-MOPS and CMIP6 models. The cumulative changes in air–sea CO₂ flux over 1850 to 2014 and changes in the carbon storage between 1850 and 2014 in piControl runs are subtracted from the historical runs. Zonally integrated (a) cumulative air–sea CO₂ flux and (b) carbon storage. Blue and red lines stand for mean values of Hist and ESM-Hist simulation ensembles of FOCI-MOPS, respectively, and black lines and shaded areas represent the mean and ± standard deviation of some CMIP6 models, including ACCESS-ESM1-5, CanESM5-CanOE, CanESM5, CESM2, CESM2-WACCM, CNRM-ESM2-1, GFDL-CM4, GFDL-ESM4, IPSL-CM6A-LR, MPI-ESM1-2-HR, NorESM2-LM, UKESM1-0-LL, and MIROC-ES2L. Results from MPI-ESM1-2-LR and MRI-ESM2-0 are also included for the carbon storage (Frölicher et al., 2015; Terhaar et al., 2021).

In Hist and ESM-Hist, surface air temperature anomalies remain indistinguishable from piControl and ESM-piControl and only start to increase from the 1920s onwards. The temperature anomalies during 2005–2014 are 0.80±0.05, 1.04±0.12, and 0.93 ^∘C in Hist, ESM-Hist, and observations, respectively (Fig. 20c).

The heat content integrated from 0 to 2000 m during 2005–2014 in Hist increased about 257±6 zeta joules (ZJ) from 1955–1964 period, very close to the observations (260–267 ZJ). The increase in ESM-Hist is 317±3 ZJ, 60 ZJ higher than in the Hist, and reflects the higher surface temperature anomalies, which presumably are a consequence of the higher atmospheric CO₂ concentrations in experiment ESM-Hist.

In general, Hist and ESM-Hist have a similar temporal evolution and response to anthropogenic forcing, and both largely agree with observations in terms of air–sea CO₂ fluxes, temperature, and ocean heat content.

3.4.2 Environmental drivers of marine biogeochemical changes

In addition to the carbon cycle, of particular interest in this study is the impact of climate change on marine biogeochemistry. Several drivers are considered important in this respect, including sea surface temperature (SST), seawater pH, oxygen (O₂) and nitrate (NO₃) levels, and primary production (PP) (Kwiatkowski et al., 2020). Here we show the temporal evolution of the anomalies of these variables relative to the 1870–1899 reference period, as described in Kwiatkowski et al. (2020). Overall, the anomalies in both Hist and ESM-Hist simulations agree with the range of the CMIP6 mean (Fig. 22). SST anomalies in Hist and ESM-Hist are comparable during the historical period until the year 2000. Between 2005 and 2014, the SST anomaly in ESM-Hist amounts to 0.70±0.13 ^∘C, slightly higher than 0.54±0.05 ^∘C in Hist (Fig. 22a). In 2014, the 0.23 ^∘C higher SST in ESM-Hist compared to Hist is in line with the higher surface air temperature (Fig. 20c). Between 2005 and 2014, the sea surface pH values in Hist and ESM-Hist are 0.095±0.000 and 0.099±0.002 units lower than 1870–1899 (Fig. 22b), respectively. The pH is slightly lower in ESM-Hist than in the Hist simulations due to a higher air–sea CO₂ flux (Fig. 20b). The O₂ anomalies averaged between 100 and 600 m are $-\mathrm{2.92}\pm \mathrm{0.58}$ mmol O₂ m⁻³ in Hist and $-\mathrm{3.03}\pm \mathrm{1.05}$ mmol O₂ m⁻³ in ESM-Hist simulations between 2005 and 2014 (Fig. 22c), which is equivalent to a decrease of 2 % in both simulations in this depth range. When it comes to the changes in the global O₂ inventory, FOCI-MOPS simulated a decrease of 284±32 Tmol O₂ per decade over 1960–2010. While the rate is substantially lower than what observations suggest (961±429 Tmol O₂ per decade; Schmidtko et al., 2017), such an underestimation of observationally estimated changes in the marine O₂ inventory is a common feature among Earth system models (Oschlies et al., 2018). The NO₃ concentrations in the euphotic zone (0–100 m) are negatively correlated with SST (Fig. 22d), suggesting a reduced surface supply of NO₃ with increasing ocean stratification, a relationship that also exists in other global models (Fu et al., 2016). Similar to the temporal evolution of SST, the NO₃ remains relatively unchanged until the 2000s, and the anomalies between 2005 and 2014 are $-\mathrm{0.08}\pm \mathrm{0.07}$ and $-\mathrm{0.25}\pm \mathrm{0.11}$ mmol N m⁻³ in the Hist and ESM-Hist simulations, respectively (Fig. 22d). Temporal evolutions of PP in Hist and ESM-Hist are similar, and the averaged anomalies over 2005–2014 are $-\mathrm{1.3}\pm \mathrm{0.4}$ % and $-\mathrm{1.8}\pm \mathrm{0.8}$ % in the Hist and ESM-Hist simulations, respectively. These are at the upper end of the PP decrease projected with CMIP6 models (Fig. 22e).

Figure 22 Temporal evolution of global mean (a) sea surface temperature (SST), (b) sea surface pH, (c) oxygen (O₂) between 100 and 600 m, (d) nitrate (NO₃) in the upper 100 m, and (e) vertically integrated primary production (PP) in percentage. Values are relative to the mean of the time period 1870–1899 after subtracting the drifts of the piControl runs. Blue and red lines stand for mean values of Hist and ESM-Hist simulation ensembles of FOCI-MOPS, respectively, and black lines and shaded areas represent the mean and ± standard deviation of CMIP6 models, respectively, as described in Kwiatkowski et al. (2020).

4 Conclusions

In this study we present the implementation and evaluation of the marine biogeochemistry component coupled to FOCI. The resulting FOCI-MOPS model is based on MOPS (“Model of Oceanic Pelagic Stoichiometry”; Kriest and Oschlies, 2015), which simulates the elemental cycles of oceanic phosphorus, nitrogen, and oxygen between their dissolved (PO₄, NO₃, O₂, and DOM) and particulate (phytoplankton, zooplankton, and detritus) pools. DIC and ALK are included in the implementation, providing a fully coupled carbon cycle in FOCI-MOPS. Spin-up (spinup) and an ensemble of three pre-industrial control (piControl) and historical (Hist) experiments were performed with prescribed atmospheric CO₂ concentrations, following the CMIP6 protocols (Eyring et al., 2016). All tracers and fluxes approached steady state or only showed small drifts relative to their mean concentrations in the end of the 500-year spin-up. The marine carbon inventory decreased 0.086 Pg C per year during the last 100 years in the spin-up, which is smaller than the drift suggested as “acceptable” (<0.1 Pg C yr⁻¹) in the Coupled Climate–Carbon Cycle Model Intercomparison Project (C4MIP; Jones et al., 2016). Based on the applied metrics in Table 2, we conclude that FOCI-MOPS reproduces observed biogeochemical patterns of inorganic tracers and phytoplankton/chlorophyll well and that the model's large-scale performance is comparable with other CMIP models. Overall, globally integrated biogeochemical fluxes are in line with observational estimates (Table 3). The spatial patterns of surface ocean pCO₂ and air–sea CO₂ fluxes agree well with observations (Fig. 18). Simulated changes due to the increasing atmospheric CO₂ in globally integrated air–sea CO₂ flux and ocean heat content (0–2000 m), as well as globally averaged surface air temperature, also agree well with observations.

We also performed CO₂ emission-driven experiments, including ESM-spinup, ESM-piControl, and ESM-Hist. The air–sea CO₂ flux in these simulations is slightly higher in the historical emission-driven simulation (ESM-Hist) than in the atmospheric-CO₂-driven simulation (Hist), resulting in a higher cumulative air–sea CO₂ flux of 139 Pg C relative to 129 Pg C in the Hist between 1850 and 2014. The difference in CO₂ perturbation forcing contributes to a higher atmospheric CO₂ until the end of the simulation period in ESM-Hist than in the Hist simulations. This higher atmospheric CO₂ likely contributes to the higher cumulative air–sea CO₂ flux, as well as the higher surface temperatures and ocean heat content in ESM-Hist. Concerning the environmental drivers of marine biogeochemical changes (sea surface temperature, seawater pH, oxygen and nitrate levels, and primary production), their anomalies in both Hist and ESM-Hist simulations agree with the range of the CMIP6 model results (Fig. 22).

We plan to investigate possible shortcomings of the simulated ventilation by direct comparison of simulated and observed abiotic transient tracers (CFCs and SF6). This will better constrain shortcomings in the biogeochemical model component's parameterizations of export and remineralization. We will implement an iron model (Somes et al., 2021) in FOCI-MOPS to improve the potential model deficiency caused by lacking iron limitation in phytoplankton growth. Sensitivity runs with altered remineralization schemes are currently under way. In addition, evaluations of surface seasonal cycle of the model will also be carried out. We plan to use the model to investigate ocean-based CO₂ removal approaches for climate change mitigation, such as ocean alkalinity enhancement. While this can be done to some extent with the current model, new parameterizations and improvements in the simulation of ocean carbonate chemistry will likely be required. However, some of this work can only be done when more experimental data is available to constrain the model.

We note that FOCI-MOPS is less complex than most of the biogeochemistry components employed in other Earth system models (Séférian et al., 2020), for instance it does not explicitly resolve silicate or iron cycles and does not consider the CaCO₃ saturation state in the dissolution scheme. Nevertheless, as shown by our evaluation in this study, FOCI-MOPS overall shows an adequate performance that makes it an appropriate tool for studying marine biogeochemistry and biogeochemistry–climate interactions in different climate change scenarios, in order to inform the development of emission pathways that are consistent with the agreed climate targets.

Appendix A: Biogeochemical model equations and parameters

The biogeochemical model describes the cycling of phosphorus, nitrogen and oxygen in a stoichiometrically consistent manner (“Model of Oceanic Pelagic Stoichiometry”; Kriest and Oschlies, 2015). The model contains seven components, of which five are calculated in phosphorus units, namely phytoplankton (PHY), zooplankton (ZOO), detritus, dissolved organic matter (DOM), and phosphate (PO₄). Additionally, nitrate (NO₃) and oxygen (O₂) are simulated, with biogeochemical interactions among the different elements coupled via fixed stoichiometric ratios. As the model also simulates NO₃ loss and gain through denitrification and nitrogen fixation, respectively, the nitrate-to-phosphate ratio can vary. Details of the model can be found in Kriest and Oschlies (2015), and in Sect. A1 we briefly describe the equations. We have coupled a carbon cycle, which involves the effects of biogeochemical interactions, calcite formation and dissolution on the two additional components of dissolved inorganic carbon (DIC) and alkalinity (ALK), and air–sea gas exchange of CO₂ across the sea surface. The implementation of the carbon cycle is described in Sect. A2. Benthic–pelagic exchanges are described in Sect. A3, and a summary of all equations is given in Sect. A4. The biogeochemical model parameters and details on their calibration can be found in Sect. A5.

A1 The biogeochemical model MOPS

A1.1 Plankton dynamics

Phytoplankton (PHY) growth depends on temperature T (f₁(T)), daylight intensity (I, (W m⁻² d⁻¹)), day length τ (d) (f₂(I,τ)), and nutrients (f₃(PO₄,NO₃)). Temperature-dependent growth is formulated following Eppley (1972), in the notation by Schmittner et al. (2008):

$$\begin{array}{}\text{(A1)}& f_{\mathrm{1}}\left(T\right)={\mathit{\mu }}_{\mathrm{PHY}}^{\prime }{e}^{\frac{T}{\mathrm{15.65}}},\end{array}$$

where ${\mathit{\mu }}_{\mathrm{PHY}}^{\prime }=\mathrm{0.6}$ (d⁻¹) is the maximum growth rate of phytoplankton at T=0 ^∘C. Light limitation f(I) is parameterized following Smith (1936), integrated over vertical box thickness and day according to Evans and Parslow (1985), in the notation by Evans and Garçon (1997). We note that while the formulation by Evans and Garçon (1997) is based on maximum growth rate and the initial slope of the P–I curve, α ((W m⁻²)⁻¹ d⁻¹), we here calculate the light limitation based on the half-saturation constant for light, I_c=9.653 (W m⁻²), as expressed through $I_{\mathrm{c}}={\mathit{\mu }}_{\mathrm{PHY}}^{\prime }/\mathit{\alpha }$ (thereby resulting in α=0.0622 ((W m⁻²)⁻¹d⁻¹)):

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}{f}_{\mathrm{2}}\left(I\right)& ={\displaystyle \frac{\mathit{\tau }}{\mathrm{\Delta }z(k_{\mathrm{w}}+k_c\mathrm{PHY})}}\left(\mathit{\varphi }\left({\displaystyle \frac{\mathrm{2}I}{I_{\mathrm{c}}\mathit{\tau }}}\right)\right.\\ & \left.-\mathit{\varphi }\left({\displaystyle \frac{\mathrm{2}I}{I_{\mathrm{c}}\mathit{\tau }}}{e}^{-\mathrm{\Delta }z(k_w+k_c\mathrm{PHY})}\right)\right)\end{array},\end{array}$$

where I (W m⁻² d⁻¹) is the total light during a day at the top of each vertical layer, τ is the length of a day in (fraction of) days, Δz is the vertical box thickness (m), and k_w=0.04 (m⁻¹) and k_c=0.48 ((mmol P m⁻³)⁻¹ m⁻¹) are the attenuation coefficients of water and phytoplankton, respectively. The function ϕ is given by (Evans and Garçon, 1997)

$$\begin{array}{}\text{(A3)}& \mathit{\varphi }\left(u\right)=\ln \left(u+\sqrt{\mathrm{1}+u^{\mathrm{2}}}\right)-{\displaystyle \frac{\sqrt{\mathrm{1}+u^{\mathrm{2}}}-\mathrm{1}}{u}}.\end{array}$$

The dependence of growth on nutrients f₃(PO₄,NO₃) is based on a Monod function of the least available nutrient X, assuming a constant stoichiometry of phytoplankton given by d=16 (mol N : mol P):

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}{f}_{\mathrm{3}}({\mathrm{PO}}_{\mathrm{4}},{\mathrm{NO}}_{\mathrm{3}})& ={\displaystyle \frac{X}{k_{\mathrm{PHY}}+X}}\\ \mathrm{with}X& =\min \left({\mathrm{PO}}_{\mathrm{4}},{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{d}}\right),\end{array}\end{array}$$

where k_PHY=0.031 (mmol P m⁻³) is the half-saturation constant for phosphate (Kriest et al., 2017). Total phytoplankton growth PP (mmol P m⁻³ d⁻¹) is then given by the minimum of light and nutrient limitation, but this is only done if the most limiting nutrient is above a certain threshold ($X>{\mathrm{P}}^{*}={\mathrm{10}}^{-\mathrm{6}}$ (mmol P m⁻³)):

$$\begin{array}{}\text{(A5)}& \mathrm{PP}=f_{\mathrm{1}}\left(T\right)\mathrm{PHYmin}\left(f_{\mathrm{2}}\right(I),f_{\mathrm{3}}({\mathrm{PO}}_{\mathrm{4}},{\mathrm{NO}}_{\mathrm{3}}\left)\right).\end{array}$$

Phytoplankton experiences a linear loss given by λ_PHY=0.03 (d⁻¹). It is grazed by zooplankton (ZOO), where grazing G described by a Holling-III function, with a maximum grazing rate μ_ZOO=1.893 (d⁻¹) and half-saturation constant K_ZOO=0.086 (mmol P m⁻³):

$$\begin{array}{}\text{(A6)}& G={\mathit{\mu }}_{\mathrm{ZOO}}\mathrm{ZOO}{\displaystyle \frac{{\mathrm{PHY}}^{\mathrm{2}}}{K_{\mathrm{ZOO}}^{\mathrm{2}}+{\mathrm{PHY}}^{\mathrm{2}}}}.\end{array}$$

Only a fraction ϵ_ZOO=0.75 of grazing G (in (mmol P m⁻³ d⁻¹)) is effectively ingested, the rest is released again via egestion. Zooplankton further experiences a quadratic mortality κ_ZOO=4.548 ((mmol P m⁻³)⁻¹ d⁻¹) and a linear excretion rate given by λ_ZOO=0.03 (d⁻¹). Phytoplankton and zooplankton further die with a constant mortality rate of ${\mathit{\lambda }}_{\mathrm{PHY}}^{\prime }={\mathit{\lambda }}_{\mathrm{ZOO}}^{\prime }=\mathrm{0.01}$ (d⁻¹) but only when present above the lower concentration threshold P${}^{*}={\mathrm{10}}^{-\mathrm{6}}$ (mmol P m⁻³). Therefore, the source-minus-sink terms due to biogeochemical interactions for phytoplankton ($S_{\mathrm{PHY}}^{\mathrm{Bio}}$) and zooplankton ($S_{\mathrm{ZOO}}^{\mathrm{Bio}}$) are

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{PHY}}^{\mathrm{Bio}}& =\mathrm{PP}-G-{\mathit{\lambda }}_{\mathrm{PHY}}\mathrm{PHY}\\ & -{\mathit{\lambda }}_{\mathrm{PHY}}^{\prime }\max (\mathrm{0},\mathrm{PHY}-{\mathrm{P}}^{*})\end{array},\text{(A8)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{ZOO}}^{\mathrm{Bio}}& ={\mathit{ϵ}}_{\mathrm{ZOO}}G-{\mathit{\lambda }}_{\mathrm{ZOO}}\mathrm{ZOO}-{\mathit{\kappa }}_{\mathrm{ZOO}}{\mathrm{ZOO}}^{\mathrm{2}}\\ & -{\mathit{\lambda }}_{\mathrm{ZOO}}^{\prime }\max (\mathrm{0},\mathrm{ZOO}-{\mathrm{P}}^{*}).\end{array}\end{array}$$

Note that whereas Kriest and Oschlies (2015) assumed that plankton cycling only occurs in the upper, well-lit waters, we here skip this restriction and compute plankton interactions over the full water column. (Of course, production will cease in the aphotic zone because of the light limitation.) Further, we avoid possible negative concentrations because of the Eulerian time-stepping by computing biogeochemical fluxes only when plankton concentrations are positive.

A1.2 DOP and detritus

Dissolved organic matter is implicitly represented by the DOP in phosphorus units in a C : N : P molar ratio of $\mathrm{117}:\mathrm{16}:\mathrm{1}$ in the model. We assume that a fraction σ_DOP=0.15 of egestion, quadratic zooplankton mortality, and phytoplankton loss given by λ_PHY is released as DOP; the rest becomes detritus. Further, phytoplankton and zooplankton mortality leads immediately to the production of DOP. The DOP is remineralized in all layers with a constant rate ${\mathit{\lambda }}_{\mathrm{DOP}}^{\prime }=\mathrm{0.17}$ (yr⁻¹) but only when present above the lower limit of ${\mathrm{P}}^{*}={\mathrm{10}}^{-\mathrm{6}}$ (mmol P m⁻³). Its oxic and suboxic remineralization further depends on the availability of oxidants oxygen and nitrate, as described by the terms $S_{\mathrm{DOP}}^{\mathrm{Rox}}$ and $S_{\mathrm{DOP}}^{\mathrm{Rsubox}}$, which are described in detail in Sect. A1.3 below. Thus, the source-minus-sink term for DOP, $S_{\mathrm{DOP}}^{\mathrm{Bio}}$, is

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{S}_{\mathrm{DOP}}^{\mathrm{Bio}}& ={\mathit{\sigma }}_{\mathrm{DOP}}\left[\right(\mathrm{1}-{\mathit{ϵ}}_{\mathrm{ZOO}})G+{\mathit{\kappa }}_{\mathrm{ZOO}}{\mathrm{ZOO}}^{\mathrm{2}}\\ & +{\mathit{\lambda }}_{\mathrm{PHY}}\mathrm{PHY}]+{\mathit{\lambda }}_{\mathrm{PHY}}^{\prime }\max (\mathrm{0},\mathrm{PHY}-{\mathrm{P}}^{*})\\ & +{\mathit{\lambda }}_{\mathrm{ZOO}}^{\prime }\max (\mathrm{0},\mathrm{ZOO}-{\mathrm{P}}^{*})-S_{\mathrm{DOP}}^{\mathrm{Rox}}-S_{\mathrm{DOP}}^{\mathrm{Rsubox}}\end{array}.\end{array}$$

Like DOP, detritus (DET) is remineralized with a fixed rate ${\mathit{\lambda }}_{\mathrm{DET}}^{\prime }=\mathrm{0.05}$ (d⁻¹) to nutrients when present above the lower limit of ${\mathrm{P}}^{*}={\mathrm{10}}^{-\mathrm{6}}$ (mmol P m⁻³), and its decomposition depends on oxygen and nitrate, as described for the terms $S_{\mathrm{DET}}^{\mathrm{Rox}}$ and $S_{\mathrm{DET}}^{\mathrm{Rsubox}}$ in Sect. A1.3 below. In addition, detritus sinks through the water column. We assume that the sinking speed of detritus increases linearly with depth, according to $w=w\left(z\right)=w_{\mathrm{lin}}z$, where z is the centre of a layer. In steady state and in the absence of any other processes, this parameterization can be regarded as equivalent to the so-called “Martin” (power law) curve of particle flux, with the exponent b given by $b={\mathit{\lambda }}_{\mathrm{DET}}^{\prime }/w_{\mathrm{lin}}$ (see Kriest and Oschlies, 2008, for a detailed discussion). For easier comparison with other model studies, which explicitly define b, and for comparison with empirically observed values for this parameter, in our model experiments we prescribe b=1.41309 and evaluate w_lin from it via $w_{\mathrm{lin}}={\mathit{\lambda }}_{\mathrm{DET}}^{\prime }/b$. Note that in MOPS, due to reduction of remineralization by lack of oxidants (Sect. A1.3), the local effective “Martin” exponent b may be smaller than initially prescribed. The source-minus-sink term of detritus, $S_{\mathrm{DET}}^{\mathrm{Bio}}$, is therefore

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{S}_{\mathrm{DET}}^{\mathrm{Bio}}& =(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{DOP}})\left[\right(\mathrm{1}-{\mathit{ϵ}}_{\mathrm{ZOO}})G+{\mathit{\kappa }}_{\mathrm{ZOO}}{\mathrm{ZOO}}^{\mathrm{2}}\\ & +{\mathit{\lambda }}_{\mathrm{PHY}}\mathrm{PHY}]-S_{\mathrm{DET}}^{\mathrm{Rox}}-S_{\mathrm{DET}}^{\mathrm{Rsubox}}+{\displaystyle \frac{\partial w\mathrm{DET}}{\partial z}}\end{array}.\end{array}$$

A1.3 Oxic and suboxic remineralization

If oxygen is above a threshold defined by ${{\mathrm{O}}_{\mathrm{2}}}^{*}=max(\mathrm{0},{\mathrm{O}}_{\mathrm{2}}-{{\mathrm{O}}_{\mathrm{2}}}^{\min })$, with ${{\mathrm{O}}_{\mathrm{2}}}^{\min }=\mathrm{1}$ (mmol O₂ m−3), organic matter is remineralized aerobically according to a sigmoidal function:

$$\begin{array}{}\text{(A11)}& l_{{\mathrm{O}}_{\mathrm{2}}}={\displaystyle \frac{({{\mathrm{O}}_{\mathrm{2}}}^{*}{)}^{\mathrm{2}}}{({{\mathrm{O}}_{\mathrm{2}}}^{*}{)}^{\mathrm{2}}+K_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{2}}}},\end{array}$$

where $K_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{1.066}$ (mmol O₂ m−3) is the half-saturation constant for the heterotroph's uptake of oxygen. To prevent total oxygen consumption per time step from exceeding available oxygen, we first calculate the theoretical oxygen demand for aerobic remineralization of detritus and DOP $u_{{\mathrm{O}}_{\mathrm{2}}}$:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{u}_{{\mathrm{O}}_{\mathrm{2}}}& =l_{{\mathrm{O}}_{\mathrm{2}}}\left[{\mathit{\lambda }}_{\mathrm{DET}}^{\prime }max(\mathrm{0},\mathrm{DET}-{\mathrm{P}}^{*})\right.\\ & +\left.{\mathit{\lambda }}_{\mathrm{DOP}}^{\prime }max(\mathrm{0},\mathrm{DOP}-{\mathrm{P}}^{*})\right]R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}\mathrm{\Delta }t\end{array},\end{array}$$

where ${\mathit{\lambda }}_{\mathrm{DOP}}^{\prime }$ and ${\mathit{\lambda }}_{\mathrm{DET}}^{\prime }$ are the remineralization rates of DOP and detritus, respectively, and Δt is the time step length of the biogeochemical model. $R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}=\mathrm{165.08044}$ denotes the stoichiometric oxygen demand of aerobic remineralization. The aerobic decay rate limitation is then

$$\begin{array}{}\text{(A13)}& s_{{\mathrm{O}}_{\mathrm{2}}}=l_{{\mathrm{O}}_{\mathrm{2}}}{\displaystyle \frac{min({\mathrm{O}}_{\mathrm{2}}^{*},u_{{\mathrm{O}}_{\mathrm{2}}})}{u_{{\mathrm{O}}_{\mathrm{2}}}}}.\end{array}$$

Therefore, the sinks of oxygen due to remineralization of DOP and detritus are defined by

$$\begin{array}{}\text{(A14)}& {\displaystyle }{S}_{\mathrm{DOP}}^{\mathrm{Rox}}=-{\mathit{\lambda }}_{\mathrm{DOP}}^{\prime }max(\mathrm{0},\mathrm{DOP}-{\mathrm{P}}^{*})s_{{\mathrm{O}}_{\mathrm{2}}},\text{(A15)}& {\displaystyle }{S}_{\mathrm{DET}}^{\mathrm{Rox}}=-{\mathit{\lambda }}_{\mathrm{DET}}^{\prime }max(\mathrm{0},\mathrm{DET}-{\mathrm{P}}^{*})s_{{\mathrm{O}}_{\mathrm{2}}}.\end{array}$$

If ${{\mathrm{O}}_{\mathrm{2}}}^{*}$ is lower than 36 mmol O₂ m⁻³, denitrification additionally sets in. As for oxygen, we first define a quadratic rate limitation of this process based on a minimum concentration of nitrate, ${{\mathrm{NO}}_{\mathrm{3}}}^{\min }=\mathrm{15.978}$ (mmol N m⁻³), with ${{\mathrm{NO}}_{\mathrm{3}}}^{*}=max(\mathrm{0},{\mathrm{NO}}_{\mathrm{3}}-{{\mathrm{NO}}_{\mathrm{3}}}^{\min })$. To account for inhibition of denitrification by oxygen we further reduce this rate by the inverse oxygen consumption rate:

$$\begin{array}{}\text{(A16)}& l_{{\mathrm{NO}}_{\mathrm{3}}}={\displaystyle \frac{({{\mathrm{NO}}_{\mathrm{3}}}^{*}{)}^{\mathrm{2}}}{({{\mathrm{NO}}_{\mathrm{3}}}^{*}{)}^{\mathrm{2}}+K_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{2}}}}\left(\mathrm{1}-l_{{\mathrm{O}}_{\mathrm{2}}}\right),\end{array}$$

where $K_{{\mathrm{NO}}_{\mathrm{3}}}=\mathrm{23.104}$ (mmol N m⁻³) is the half-saturation constant for the denitrifiers' uptake of nitrate. As for oxygen, we restrict the use of nitrate to the amount available:

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}{u}_{{\mathrm{NO}}_{\mathrm{3}}}& =l_{{\mathrm{NO}}_{\mathrm{3}}}\left[{\mathit{\lambda }}_{\mathrm{DET}}^{\prime }max(\mathrm{0},\mathrm{DET}-{\mathrm{P}}^{*})\right.\\ & \left.+{\mathit{\lambda }}_{\mathrm{DOP}}^{\prime }max(\mathrm{0},\mathrm{DOP}-{\mathrm{P}}^{*})\right]R_{-{\mathrm{NO}}_{\mathrm{3}}:\mathrm{P}}\mathrm{\Delta }t\end{array},\end{array}$$

with $R_{-{\mathrm{NO}}_{\mathrm{3}}:\mathrm{P}}$ =0.8 $R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}-d=\mathrm{116.064352}$ (mmol NO₃: mmol P),following the stoichiometry of Paulmier et al. (2009). The rate limitation of anaerobic decay is then

$$\begin{array}{}\text{(A18)}& s_{{\mathrm{NO}}_{\mathrm{3}}}=l_{{\mathrm{NO}}_{\mathrm{3}}}{\displaystyle \frac{min({{\mathrm{NO}}_{\mathrm{3}}}^{*},u_{{\mathrm{NO}}_{\mathrm{3}}})}{u_{{\mathrm{NO}}_{\mathrm{3}}}}}.\end{array}$$

Therefore, the sinks of oxygen and nitrate due to denitrification of DOP and detritus are defined by

$$\begin{array}{}\text{(A19)}& {\displaystyle }{S}_{\mathrm{DOP}}^{\mathrm{Rsubox}}=-{\mathit{\lambda }}_{\mathrm{DOP}}^{\prime }max(\mathrm{0},\mathrm{DOP}-{\mathrm{P}}^{*})s_{{\mathrm{NO}}_{\mathrm{3}}},\text{(A20)}& {\displaystyle }{S}_{\mathrm{DET}}^{\mathrm{Rsubox}}=-{\mathit{\lambda }}_{\mathrm{DET}}^{\prime }max(\mathrm{0},\mathrm{DET}-{\mathrm{P}}^{*})s_{{\mathrm{NO}}_{\mathrm{3}}}.\end{array}$$

A1.4 Nitrogen fixation

As in Kriest and Oschlies (2015), nitrogen fixation depends on temperature and nutrient ratio:

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}{S}_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{NFix}}& ={\mathit{\mu }}_{\mathrm{NFix}}max\left(\mathrm{0},{\displaystyle \frac{t_{\mathrm{2}}{T}^{\mathrm{2}}+t_{\mathrm{1}}T-t_{\mathrm{0}}}{t_f}}\right)\\ & max\left(\mathrm{0},\mathrm{1}-{\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{d{\mathrm{PO}}_{\mathrm{4}}}}\right),\end{array}\end{array}$$

with μ_NFix=1.88924 (μmol N m⁻³ d⁻¹) being the maximum nitrogen fixation of the (implicit) cyanobacteria and t₂, t₁, t₀, and t_f coefficients that describe the temperature dependency of nitrogen fixation of Trichodesmium spp. (Breitbarth et al., 2007), using the approximation by Kriest and Oschlies (2015). We note that in the model nitrogen fixation only occurs when PO${}_{\mathrm{4}}>{\mathrm{10}}^{-\mathrm{6}}$ mmol P m⁻³.

A1.5 Nutrients and oxygen

Phosphate (PO₄) is affected by primary production, excretion by zooplankton, and the decay of dissolved and particulate organic matter (POM), as explained above, leading to a source-minus-sink term due to biogeochemical interactions, $S_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Bio}}$, of

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}{S}_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Bio}}& =\left(-\mathrm{PP}+{\mathit{\lambda }}_{\mathrm{ZOO}}\mathrm{ZOO}\right)+S_{\mathrm{DOP}}^{\mathrm{Rox}}+S_{\mathrm{DOP}}^{\mathrm{Rsubox}}\\ & +S_{\mathrm{DET}}^{\mathrm{Rox}}+S_{\mathrm{DET}}^{\mathrm{Rsubox}}.\end{array}\end{array}$$

The loss and gain of nitrate (NO₃) follows that of phosphate for aerobic processes and production. In addition, this tracer is affected by denitrification (fixed nitrogen loss) and nitrogen fixation (fixed nitrogen gain):

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}{S}_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{Bio}}& =S^{\mathrm{NFix}}-d\mathrm{PP}+d\left({\mathit{\lambda }}_{\mathrm{ZOO}}\mathrm{ZOO}+S_{\mathrm{DOP}}^{\mathrm{Rox}}\right.\\ & \left.+S_{\mathrm{DET}}^{\mathrm{Rox}}\right)-R_{-{\mathrm{NO}}_{\mathrm{3}}:\mathrm{P}}\left(S_{\mathrm{DOP}}^{\mathrm{Rsubox}}+S_{\mathrm{DET}}^{\mathrm{Rsubox}}\right)\end{array}.\end{array}$$

Finally, oxygen (O₂) increases due to photosynthesis and decreases because of aerobic remineralization and respiration by zooplankton.

$$\begin{array}{}\text{(A24)}& S_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{Bio}}=R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}\left(\mathrm{PP}-{\mathit{\lambda }}_{\mathrm{ZOO}}\mathrm{ZOO}-S_{\mathrm{DOP}}^{\mathrm{Rox}}-S_{\mathrm{DET}}^{\mathrm{Rox}}\right)\end{array}$$

In addition, oxygen exchanges with the atmosphere at the sea surface (i.e. for layer 1) are given as $S_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{Air}}$ following Orr et al. (2017).

A2 The carbon cycle

A2.1 Coupling to the biogeochemical core

Photosynthesis decreases DIC, whereas remineralization of organic matter to phosphate increases it. We assume a constant stoichiometry between phosphorus and carbon, a=117 (mol C : mol P); DIC thus changes in proportion to changes in phosphate:

$$\begin{array}{}\text{(A25)}& S_{\mathrm{DIC}}^{\mathrm{Bio}}=a{S}_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Bio}}.\end{array}$$

Changes in phosphate and nitrate further affect alkalinity via

$$\begin{array}{}\text{(A26)}& S_{\mathrm{ALK}}^{\mathrm{Bio}}=-S_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Bio}}-S_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{Bio}}.\end{array}$$

A2.2 Calcite production and dissolution

We assume that a constant fraction $p^{{\mathrm{CaCO}}_{\mathrm{3}}}$ of detritus production via zooplankton egestion and plankton mortality to detritus is in the form of calcite.

$$\begin{array}{}\text{(A27)}& \begin{array}{rl}{p}^{{\mathrm{CaCO}}_{\mathrm{3}}}& =a{\mathit{\sigma }}_{{\mathrm{CaCO}}_{\mathrm{3}}}(\mathrm{1}-{\mathit{\sigma }}_{\mathrm{DOP}})\left[(\mathrm{1}-{\mathit{ϵ}}_{\mathrm{ZOO}})G\right.\\ & \left.+{\mathit{\kappa }}_{\mathrm{ZOO}}{\mathrm{ZOO}}^{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{PHY}}\mathrm{PHY}\right],\end{array}\end{array}$$

where a=117 (mol C:mol P) is the molar ratio of C : P in organic matter, and ${\mathit{\sigma }}_{{\mathrm{CaCO}}_{\mathrm{3}}}=\mathrm{0.032}$ (mol CaCO₃:mol C) is the calcite-to-organic-carbon ratio. Calcite production reduces DIC by 1 and alkalinity by 2:

$$\begin{array}{}\text{(A28)}& {\displaystyle }{S}_{\mathrm{DIC}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{P}}=-p^{{\mathrm{CaCO}}_{\mathrm{3}}},\text{(A29)}& {\displaystyle }{S}_{\mathrm{ALK}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{P}}=-\mathrm{2}{p}^{{\mathrm{CaCO}}_{\mathrm{3}}}.\end{array}$$

Following Schmittner et al. (2008), we integrate the production of calcite over the entire water column:

$$\begin{array}{}\text{(A30)}& P^{{\mathrm{CaCO}}_{\mathrm{3}}}=\underset{\mathrm{0}}{\overset{\mathrm{Bottom}}{\int }}{p}^{{\mathrm{CaCO}}_{\mathrm{3}}}\mathrm{d}z.\end{array}$$

The total production of calcite is then distributed and dissolved immediately over the entire water column with an e-folding profile $D=\exp (-z/l_{{\mathrm{CaCO}}_{\mathrm{3}}})$, with $l_{{\mathrm{CaCO}}_{\mathrm{3}}}=\mathrm{4289.4}$ m, thereby affecting alkalinity and DIC:

$$\begin{array}{}\text{(A31)}& {\displaystyle }{S}_{\mathrm{DIC}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{D}}=P^{{\mathrm{CaCO}}_{\mathrm{3}}}{\displaystyle \frac{\partial D}{\partial z}},\text{(A32)}& {\displaystyle }{S}_{\mathrm{ALK}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{D}}=\mathrm{2}{P}^{{\mathrm{CaCO}}_{\mathrm{3}}}{\displaystyle \frac{\partial D}{\partial z}},\end{array}$$

and thus the total source-minus-sink terms for DIC and alkalinity due to calcite formation and dissolution are

$$\begin{array}{}\text{(A33)}& {\displaystyle }{S}_{\mathrm{DIC}}^{{\mathrm{CaCO}}_{\mathrm{3}}}=S_{\mathrm{DIC}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{P}}+S_{\mathrm{DIC}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{D}},\text{(A34)}& {\displaystyle }{S}_{\mathrm{ALK}}^{{\mathrm{CaCO}}_{\mathrm{3}}}=S_{\mathrm{ALK}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{P}}+S_{\mathrm{ALK}}^{{\mathrm{CaCO}}_{\mathrm{3}}\mathrm{D}}.\end{array}$$

We note that the redistribution with an e-folding profile over all layers can result in some upward transport of the alkalinity gain caused by calcite dissolution if calcite-bearing detritus produced at greater depths dissolves further up in the water column. It may, however, be just a small problem, as most detritus will likely be produced in shallow layers, where zooplankton grazing and mortality is high.

A2.3 Air–sea gas exchange of CO₂ and solution of the carbonate system

The simulation of the carbonate chemistry system and of the air–sea gas exchange of CO₂ in MOPS follows an OCMIP-type formulation (Orr et al., 1999) with updates from Orr et al. (2017) for the air–sea gas exchange of CO₂.

The air–sea CO₂ flux is computed as follows:

$$\begin{array}{}\text{(A35)}& F_{{\mathrm{CO}}_{\mathrm{2}}}=k_{\mathrm{w}}\left([{{\mathrm{CO}}_{\mathrm{2}}}^{*}{]}_{\mathrm{sat}}-[{{\mathrm{CO}}_{\mathrm{2}}}^{*}]\right),\end{array}$$

where k_w (in m s⁻¹) is the gas transfer velocity, (CO${}_{\mathrm{2}}^{*}$)_sat (in mol kg⁻¹) is the saturation concentration of CO₂, and (CO${}_{\mathrm{2}}^{*}$) (in mol kg⁻¹) is the surface ocean dissolved CO₂ concentration (see below for further details). The instantaneous gas transfer velocity k_w is parameterized based on Wanninkhof (2014) as a quadratic function of the 10 m wind speed u:

$$\begin{array}{}\text{(A36)}& k_{\mathrm{w}}=\mathit{\xi }{\left({\displaystyle \frac{Sc}{\mathrm{660}}}\right)}^{-\mathrm{0.5}}{u}^{\mathrm{2}}\left(\mathrm{1}-f_{\mathrm{ice}}\right),\end{array}$$

where ξ is a constant, Sc is the Schmidt number, and f_ice is the fraction of the grid cell covered by sea ice.

The saturation concentration of CO₂ in equilibrium with the water-vapour-saturated atmosphere at a total atmospheric pressure of 1 atm (i.e. (CO${}_{\mathrm{2}}^{*}$)_sat) is computed as follows:

$$\begin{array}{}\text{(A37)}& {\left[{{\mathrm{CO}}_{\mathrm{2}}}^{*}\right]}_{\mathrm{sat}}=Fx{\mathrm{CO}}_{\mathrm{2}},\end{array}$$

where xCO₂ is the mole fraction of CO₂ in dry air and F is the solubility function (Weiss and Price, 1980, Eq. (13), Table 6, column 3), which includes all non-ideal effects and fits the effects of water vapour pressure for a total atmospheric pressure of 1 atm.

Once dissolved, CO₂ reacts with seawater, forming carbonic acid (H₂CO₃), most of which dissociates into two other inorganic species, bicarbonate (HCO${}_{\mathrm{3}}^{-}$) and carbonate (CO${}_{\mathrm{3}}^{\mathrm{2}-}$) ions. CO${}_{\mathrm{2}}^{*}$ refers to the sum of dissolved CO₂ and the much less abundant H₂CO₃. The sum of the three species CO${}_{\mathrm{2}}^{*}$ + HCO${}_{\mathrm{3}}^{-}$ + CO${}_{\mathrm{3}}^{\mathrm{2}-}$ is referred to as total dissolved inorganic carbon (DIC). Their partitioning depends on seawater pH, temperature, salinity, and pressure. The pH may be calculated from DIC and seawater's ionic charge balance, formalized as total alkalinity (ALK). DIC and ALK are carried as prognostic tracers in the ocean model, and both are used, along with temperature, salinity, and nutrient concentrations, to compute CO${}_{\mathrm{2}}^{*}$ at the ocean surface.

The carbonate chemistry system is solved using the equilibrium constants recommended for best practices (Dickson et al., 2007). The total pH scale is used for all constants except Ks (which uses the “free” scale following Dickson, 1990). The model does not explicitly simulate boron, sulfate, and fluoride, which are instead estimated as function of chlorinity based on Lee et al. (2010) (for boron), Morris and Riley (1966) (for sulfate), and Riley (1965) (for fluoride). Silicate, which is also not explicitly simulated, is computed as a function of density following Orr et al. (1999).

For reasons of computational efficiency, to solve the carbonate chemistry system we use the approximate and non-iterative method proposed by Follows et al. (2006) (Eqs. 8, 11, and 12 therein). The algorithm has been shown to provide a sufficiently accurate solution in the context of a three-dimensional global ocean carbon cycle model. The algorithm uses DIC, ALK, dissolved inorganic phosphorus, silica, boron, and thermodynamic equilibrium coefficients as inputs. The species retained in the expression for alkalinity are the phosphoric, silicic, carbonic, boric, sulfuric, fluoridic, and water acid systems. Starting from an initial guess of pH deriving from the previous time step, the formula provides an updated value of pH as output, which is then used to compute CO${}_{\mathrm{2}}^{*}$ (Eq. 8 in Follows et al., 2006) and surface CO₂ fugacity as follows:

$$\begin{array}{}\text{(A38)}& f{\mathrm{CO}}_{\mathrm{2}}={\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}^{*}}{K_{\mathrm{0}}}},\end{array}$$

where K₀ is the solubility coefficient of CO₂ in seawater (Eq. 12 and column 3 of Table 1 in Weiss, 1974).

A3 Benthic burial and nutrient re-supply

As in Kriest and Oschlies (2013) and Kriest and Oschlies (2015), we assume that once sinking detritus arrives at the seafloor a fraction of it is buried in the sediments. The amount buried, $S_{\mathrm{DET}}^{\mathrm{BUR}}$ (mmol P m⁻² d⁻¹), depends on the rain rate of detritus to the sea floor (F_B, (mmol P m⁻² d⁻¹)) via

$$\begin{array}{}\text{(A39)}& S_{\mathrm{DET}}^{\mathrm{BUR}}=\left\{\begin{array}{lr}\mathrm{1.6828}{F}_{\mathrm{B}}^{\mathrm{1.799}}& k=k_{\mathrm{b}}\\ \mathrm{0}& k<k_{\mathrm{b}}\end{array},\right.\end{array}$$

where k is the index of each vertical box (counting downwards) and k_b is the last ocean box above the sea floor at each horizontal model grid point.

In contrast to Kriest and Oschlies (2013) and Kriest and Oschlies (2015), the amount buried is not integrated globally and over a year and then resupplied to the ocean via river runoff; instead, in every time step the global amount of organic phosphorus, nitrogen, and buried content is resupplied at the sea surface (in the surface layer) as phosphate, nitrate, and DIC. Thus, if $B={\int }_A{F}_{\mathrm{BUR}}\mathrm{d}a$ is the global burial in a time step integrated over the sea floor area A and V₁ is the global volume of the surface layer,

$$\begin{array}{}\text{(A40)}& S_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Supply}}=\left\{\begin{array}{lr}\frac{B}{V_{\mathrm{1}}}& k=\mathrm{1}\\ \mathrm{0}& k>\mathrm{1}\end{array},\right.\end{array}$$

We note that this approach has a slightly “fertilizing” effect even away from the river mouths. To account for the implicit, simultaneous burial of particulate organic carbon and nitrogen, the respective supply of nitrate and DIC is parameterized using constant stoichiometry:

$$\begin{array}{}\text{(A41)}& {\displaystyle }{S}_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{Supply}}=d{S}_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Supply}},\text{(A42)}& {\displaystyle }{S}_{\mathrm{DIC}}^{\mathrm{Supply}}=a{S}_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Supply}}.\end{array}$$

We also account for the equivalent of negative alkalinity (that would be consumed when organic matter was remineralized instead of being buried) by subtracting it from alkalinity at the sea surface:

$$\begin{array}{}\text{(A43)}& S_{\mathrm{ALK}}^{\mathrm{Supply}}=-S_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Supply}}-S_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{Supply}}.\end{array}$$

In contrast to organic material, we assume no burial of calcite in the sediment but dissolve all calcite arriving in the bottom box immediately.

A4 Total source-minus-sinks

Summing up the effects of all processes, we arrive at the following equations for the source-minus-sink terms:

$$\begin{array}{}\text{(A44)}& {\displaystyle }{S}_{\mathrm{PHY}}=S_{\mathrm{PHY}}^{\mathrm{Bio}}(\mathrm{Eq}.\mathrm{A}\mathrm{7}),\text{(A45)}& {\displaystyle }{S}_{\mathrm{ZOO}}=S_{\mathrm{ZOO}}^{\mathrm{Bio}}(\mathrm{Eq}.\mathrm{A}\mathrm{8}),\text{(A46)}& {\displaystyle }{S}_{\mathrm{DOP}}=S_{\mathrm{DOP}}^{\mathrm{Bio}}(\mathrm{Eq}.\mathrm{A}\mathrm{9}),\text{(A47)}& {\displaystyle }{S}_{\mathrm{DET}}=S_{\mathrm{DET}}^{\mathrm{Bio}}-S_{\mathrm{DET}}^{\mathrm{BUR}}(\mathrm{Eqs}.\mathrm{A}\mathrm{10}\mathrm{and}\mathrm{A}\mathrm{39}),\text{(A48)}& {\displaystyle }{S}_{{\mathrm{PO}}_{\mathrm{4}}}=S_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Bio}}+S_{{\mathrm{PO}}_{\mathrm{4}}}^{\mathrm{Supply}}(\mathrm{Eqs}.\mathrm{A}\mathrm{22}\mathrm{and}\mathrm{A}\mathrm{40}),\text{(A49)}& {\displaystyle }{S}_{{\mathrm{NO}}_{\mathrm{3}}}=S_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{Bio}}+S_{{\mathrm{NO}}_{\mathrm{3}}}^{\mathrm{Supply}}(\mathrm{Eqs}.\mathrm{A}\mathrm{23}\mathrm{and}\mathrm{A}\mathrm{41}),\text{(A50)}& {\displaystyle }{S}_{{\mathrm{O}}_{\mathrm{2}}}=S_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{Bio}}+S_{{\mathrm{O}}_{\mathrm{2}}}^{\mathrm{Air}}(\mathrm{Eq}.\mathrm{A}\mathrm{24}\mathrm{and}\mathrm{Orr}\mathrm{et}\mathrm{al}.,\mathrm{2017}),\end{array}$$

$$\begin{array}{}\text{(A51)}& \begin{array}{rl}{S}_{\mathrm{DIC}}& =S_{\mathrm{DIC}}^{\mathrm{Bio}}+S_{\mathrm{DIC}}^{{\mathrm{CaCO}}_{\mathrm{3}}}+S_{\mathrm{DIC}}^{\mathrm{Supply}}+S_{\mathrm{DIC}}^{\mathrm{Air}}\\ & (\mathrm{Eqs}.\mathrm{A}\mathrm{25},\mathrm{A}\mathrm{33},\mathrm{A}\mathrm{42},\mathrm{and}\mathrm{Orr}\mathrm{et}\mathrm{al}.,\mathrm{2017}),\end{array}\end{array}$$

$$\begin{array}{}\text{(A52)}& \begin{array}{rl}{S}_{\mathrm{ALK}}& =S_{\mathrm{ALK}}^{\mathrm{Bio}}+S_{\mathrm{ALK}}^{{\mathrm{CaCO}}_{\mathrm{3}}}+S_{\mathrm{ALK}}^{\mathrm{Supply}}\\ & (\mathrm{Eqs}.\mathrm{A}\mathrm{26},\mathrm{A}\mathrm{34},\mathrm{and}\mathrm{A}\mathrm{43}).\end{array}\end{array}$$

A5 Parameter calibration

The biogeochemical model parameterization is based upon a previous objective calibration of model MOPS coupled to the Transport Matrix Method (Khatiwala, 2007; Khatiwala et al., 2018), using Transport Matrices derived from the ECCO project (Kriest et al., 2020). In particular, Kriest et al. (2020) optimized the parameters for oxidant-dependent remineralization NO${}_{\mathrm{3}}^{\min }$, $K_{{\mathrm{NO}}_{\mathrm{3}}}$, $K_{{\mathrm{O}}_{\mathrm{2}}}$, the maximum nitrogen fixation rate μ_NFix, the oxygen requirement for aerobic remineralization, $R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}$ and the parameter determining the particle flux profile b against observed nutrients and oxygen at a global scale, while all other parameters were kept constant.

However, when comparing parameters optimized for different circulations, Kriest et al. (2020) noted that three of the six parameters optimized were sensitive to characteristic features of the applied circulation, as expressed through the maximum mixed layer depth, age of North Atlantic Deep Water (NADW), and outcrop area of Sub-Antarctic Mode Water (SAMW) and Antarctic Intermediate Water (AAIW). We therefore adjusted $R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}$, b, and μ_NFix to the values for the respective physical diagnostic of NEMO by using the regression shown in Fig. 6 of Kriest et al. (2020). The adjustment led to a higher value for $R_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{P}}$ (165.08044 instead of 151.1 mol O₂: mol P), a slightly lower value for b (1.41309 instead of 1.46), and a lower value of μ_NFix (1.88924 instead of 2.29 μmol N m⁻³ d⁻¹; see also Table A1).

Table A1 Biogeochemical model parameters of MOPS (see also Kriest et al., 2020) and for the carbon cycle.

To adjust the parameters regarding the calcite cycle, we have extended MOPS in the setup ECCO described by Kriest et al. (2020) to include the carbon cycle described above but with a slightly different air–sea gas exchange and computation of the carbonate system, and we have optimized ${\mathit{\sigma }}_{{\mathrm{CaCO}}_{\mathrm{3}}}$ and $l_{{\mathrm{CaCo}}_{\mathrm{3}}}$ after a spin-up of 10 years against a data set of alkalinity and preindustrial DIC. The resulting parameters are only slightly different from those applied by Schmittner et al. (2008), namely ${\mathit{\sigma }}_{{\mathrm{CaCo}}_{\mathrm{3}}}=\mathrm{0.032}$ mol CaCO₃: mol C_org (instead of 0.035 mol CaCO₃: mol C_org) and $l_{{\mathrm{CaCo}}_{\mathrm{3}}}=\mathrm{4289.4}$ m (instead of 3500 m).

Appendix B: Biogeochemical model evaluation

B1 Model postprocessing

The model geometry is based on a curvilinear grid, complicating a direct comparison to observations, which are mostly available on regular (rectangular) grids. Thus, we have mapped the model output onto a horizontal grid defined by $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ using Ferret's functions curv_to_rect_map (with a radius of 2^∘ for map creation) and curv_to_rect. The vertical grid of NEMO was maintained. All further analysis of model fit was carried out on these remapped quantities unless stated otherwise.

B2 Data sets for model validation

For a complete biogeochemical model evaluation we used data sets of phosphate, nitrate, oxygen, DIC, total alkalinity, surface chlorophyll, mesozooplankton, and particulate and dissolved organic matter. Because many of the observed quantities are available in different spatial resolutions, they were gridded onto the rectangular model grid. Further details, including conversion from different units, are described below.

B2.1 Nutrients, oxygen, DIC, and alkalinity

For the spatial distribution (Sect. 3.2.1) comparisons of dissolved inorganic tracers we used the interpolated data, and for model evaluation with statistical metrics (Sect. 3.2.3) we used the non-interpolated data of the Global Ocean Data Analysis Project version 2 mapped climatology (GLODAPv2.2016b,  Lauvset et al., 2016; Olsen et al., 2016), as available at https://www.nodc.noaa.gov/ocads/oceans/GLODAPv2_2020/ (last access: 12 May 2016). Observed concentrations of all inorganic tracers have been converted from μmol kg⁻¹ to mmol m⁻³ using in situ density computed from GLODAP's temperature and salinity. Originally the non-interpolated data set contained between 158 401 and 252 808 data points for the different tracers because the NEMO grid, onto which the data are interpolated, has a higher vertical resolution; the final data set contains between 183 213 and 295 603 data points (see Table B1).

Lauvset et al. (2016)Olsen et al. (2016)Melin (2013)Moriarty and O'Brien (2013)Martiny et al. (2014)Torres-Valdés et al. (2009)Moutin et al. (2008)Yoshimura et al. (2007)

Table B1 Statistics for observations, regridded onto model grid: number of observations, minimum and maximum concentration, volume-weighted mean, and standard deviation. See Sect. B2 for further details.

B2.2 Phytoplankton

For the assessment of simulated phytoplankton we use chlorophyll data derived from remote sensing (MODIS-Aqua; Melin, 2013, https://jeodpp.jrc.ec.europa.eu/ftp/public/JRC-OpenData/GMIS/satellite/9km/climatology/, last access: 20 January 2021). The original surface data are available as a monthly climatology on a 9 km grid. After averaging to annual mean chlorophyll, the data are gridded (by averaging) onto a horizontal grid defined by $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$. Chlorophyll was converted to carbon using the algorithm derived by Sathyendranath et al. (2009) and then to phosphorus using a C : P ratio of 117 mol C : mol P. The resulting data set contains 36.669 data points, which are all located in the surface layer, with minimum and maximum values of 0 and 0.25 mmol P m⁻³, respectively, an unweighted mean of 0.0161 mmol P m⁻³, and a standard deviation of 0.0122 mmol P m⁻³ (see Table B1).

B2.3 Zooplankton

For the evaluation of simulated zooplankton we use the MAREDAT data set of mesozooplankton (Moriarty and O'Brien, 2013). This sparse, quasi-climatological data set contains 42.245 data points of monthly mean mesozooplankton (in mg C m⁻³) on a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ degree grid. After averaging over a year, and mapping onto the spatial grid of NEMO, we obtained a total of 37.838 data points for the upper 100 m. Conversion to phosphorus was carried out by assuming a C : P ratio of 117 mol C : mol P.

Many groups of mesozooplankton carry out diurnal vertical migration, meaning that they descend to depths between ≈200–500 m depth at dawn and ascend to the surface layers for feeding at dusk (e.g., Kiko et al., 2017, 2020). This process is so far not included in the model (see Aumont et al., 2017), and may lead to an overestimation of simulated zooplankton biomass when compared only against surface data. Therefore, to account for the maximum total potential biomass of grazers in the observations, we integrated all observed biomass within the upper 500 m and distributed it evenly over the upper 100 m for model comparison.

Further, the biogeochemical model does not distinguish between microzooplankton and mesozooplankton but aggregates both types into one single component. Unfortunately, samples for microzooplankton are much more sparse (only 2029 monthly data in the data set by Buitenhuis et al., 2012) than those of mesozooplankton and are often taken at other locations and during other times. Based on an analysis at stations where both small and large zooplankton observations are available, we estimated an approximate micro-to-mesozooplankton ratio of 1. For comparison with the model we therefore doubled the observed concentrations of mesozooplankton, resulting in minimum and maximum concentrations of 0 and 0.348 mmol P m⁻³, an average of 0.0062 mmol P m⁻³, and a standard deviation of 0.0101 mmol P m⁻³ (see Table B1).

B2.4 POP

There is no direct observational equivalent to simulated detritus; the closest type of observation are those of particulate organic phosphorus (POP), nitrogen (PON), or carbon (POC). For model evaluation we downloaded the data set by Martiny et al. (2014, https://datadryad.org/stash/downloads/file_stream/87679, last access: 16 April 2020), which contains observations of POP, PON, and POC. After omitting some entries where depth was not given, we obtained 6940 data points for POP and 46 705 data points for PON. Because of the much higher data frequency for PON, we used this variable as further diagnostic, and converted it to POP using a stoichiometric ratio of 16 mol N : mol P.

To map data onto a regular grid we averaged all data that fall within boxes defined by a horizontal resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ and 23 depth intervals centred at 5, 15, 27.5, 45, 65, 87.5, 117.5, 160, 222.5, 310, 435, 610, 847.5, 1160, 1542.5, 1975, 2450, 2950, 3450, 3950, 4450, 4950 and 5450 m. The resulting data set contained 6887 data points for PON, with minimum and maximum concentrations (after conversion to POP) of 0 and 1.69 mmol P m⁻³, respectively, an average of 0.0266 mmol P m⁻³, and a standard deviation of 0.0499 mmol P m⁻³. Following this mapping, we interpolated the data onto the NEMO grid. Because of the higher vertical resolution of NEMO, this data set contains more observations but similar statistics (see Table B1).

Particulate organic matter is usually collected with Niskin bottles and then filtered; thus it entails not only detritus (dead organic particles) but also phytoplankton and possibly a fraction of zooplankton. We therefore compare these observations to the sum of simulated detritus and phytoplankton plus half of the zooplankton, thereby assuming that microzooplankton is caught in the Niskin bottles and remains on filters for PON analysis.

B2.5 DOP

Most observations of dissolved organic phosphorus applied for model evaluation have been compiled by Angela Landolfi. They include data from cruises 36N, AMT10, AMT12, AMT14, AMT15, AMT16, and AMT17 (Torres-Valdés et al., 2009); the BIOSOPE cruise (Moutin et al., 2008); a cruise of the North Atlantic (cruise D279, April–May 2004; Landolfi et al., 2016); and the unpublished data from the Indian Ocean (cruise CD139, March–April 2002; Landolfi, unpublished). The coefficient of variation of DOP measurements is about 10 % in Landolfi et al. (2016). In the compilation we only included data with a positive (good) quality flag. We further included data from the northern Pacific from Fig. 2 of Yoshimura et al. (2007). Data were gridded onto a $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid with 23 vertical layers, as applied for mapping of PON. After mapping onto this grid we obtained 814 data points in the upper 100 m, with minimum and maximum values of 0 and 3.92 mmol P m⁻³, an average of 0.178 mmol P m⁻³, and a standard deviation of 0.21 mmol P m⁻³. Gridding onto the finer NEMO grid increases the number of observations but largely maintains the statistics (see Table B1). Note that the majority of the surface observational data (phytoplankton, zooplankton, POP, and DOP) were collected between spring and autumn, and potential biases could exist when compared with annual mean model results.

B3 Metrics

To assess the performance of the biogeochemical model we apply six statistical measures and metrics that account in different ways for potential errors in average concentrations (biases), the spatial variability of observations, and the match to spatial patterns. In particular we evaluated the model bias (both absolute and normalized by the observed mean), the model's standard deviation σ_M normalized by the standard deviation of observations σ_O, the root-mean-squared error (RMSE), the pattern error or centred rms difference RMSE^′ (RMSE minus bias), and the correlation coefficient r between model and observations. All calculations take into account the spatial dimensions of the model, thereby emphasizing deviations in the deep ocean, where box thicknesses are large. To investigate the model's representation of dissolved inorganic tracers in different vertical domains, beside the global model fit we also evaluated the metrics in different vertical domains, namely for 0–100, 100–200, 200–500, 500–1000, 1000–2000, and 2000–5000 m. Vertical domains of the organic tracers are surface for phytoplankton, the upper 100 m for zooplankton and DOP (see above for treatment of zooplankton observations), and the full domain for POP.

A slight distortion in the physical model (e.g. a current being located slightly off) may cause a large RMSE and thereby induce a large model error, even if the biogeochemical model is dynamically correct. To account for this, and to compare this model with the results by Ilyina et al. (2013), we further added a seventh metric, namely the Bhattacharyya distance (BD), which evaluates the similarity between observed and simulated frequency distributions of tracers. In particular, to evaluate BD we binned simulated and observed tracer concentrations into N=50 concentration classes over a typical range of concentrations (phosphate: 0–4 mmol P m⁻³; nitrate: 0–50 mmol N m⁻³; oxygen: 0–400 mmol O₂ m⁻³; DIC and alkalinity: 1700–2500 mmol m⁻³; phytoplankton and zooplankton: 0–0.1 mmol P m⁻³; POP: 0–0.5 mmol P m⁻³; DOP: 0–1 mmol P m⁻³). Examples for different distributions are given in Figs. B1 and 16. Following binning, BD was evaluated as

$$\begin{array}{}\text{(B1)}& \mathrm{BD}=-\ln \left(\sum _{i=\mathrm{1}}^N{m}_i{o}_i\right),\end{array}$$

where m_i and o_i are the frequencies of simulated and observed model boxes with concentration of class i. We note that BD relates to the Hellinger distance (HD) via BD $=-\ln (\mathrm{1}-{\mathrm{HD}}^{\mathrm{2}})$ as follows:

$$\begin{array}{}\text{(B2)}& \mathrm{HD}=\sqrt{\mathrm{0.5}\sum _{i=\mathrm{1}}^N(\sqrt{m_i}-\sqrt{o_i}{)}^{\mathrm{2}}}.\end{array}$$

In contrast to the BD, the HD is bounded by $\mathrm{0}\le \mathrm{HD}\le \mathrm{1}$. Finally, the L1 norm of distributions as given by the following equation is evaluated (note that this norm is bounded by $\mathrm{0}\le L\mathrm{1}\le \mathrm{2}$).

$$\begin{array}{}\text{(B3)}& L\mathrm{1}=\sum _{i=\mathrm{1}}^N|m_i-o_i|\end{array}$$

Figures B1 and 16 list the different metrics of inorganic and organic tracers in different vertical domains. In all cases, the smaller the area of overlap between the two distributions, the larger the metric will be.

Figure B1 Frequency distribution of phosphate, nitrate, oxygen, DIC, and alkalinity (left to right) from observations (black filled bars) and the model (red bars) for different vertical domains (top to bottom). Numbers denote three different metrics for the similarity of the distributions, namely L1 (Eq. B3), HD (Eq. B2), and BD (Eq. B1).

Table B2 Different metrics for surface (0–100 m), total domain (0–8000 m) and range of metrics over seven different vertical domains (0–100, 100–200, 200–500, 500–1000, 1000–2000, 2000–5000, 0–8000 m) of historical experiments ensemble mean from 1972 to 2013 for PO₄, NO₃, O₂, ALK; 2002 for DIC; and 2005–2014 for Phy, Zoo, POP, and DOP. Different metrics include the correlation coefficient (r), the root-mean-squared error (RMSE), the pattern error (RMSE^′, RMSE minus global bias), the global bias (in mmol m⁻³), bias normalized by observed global mean, normalized standard deviation, and the Bhattacharyya distance (BD) for each model component. Except for BD, all metrics have been calculated on a volume-weighted basis. See the text for further details.