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 2 concentration. A series of experiments covering the historical period (1850– 2014) were performed following the DECK (Diagnostic, Evaluation and Characterization of Klima) and CMIP6 (Cou-pled Model Intercomparison Project 6) protocols. Overall, modelled biogeochemical tracer distributions and ﬂuxes, transient evolution in surface air temperature, air–sea CO 2 ﬂuxes, 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.


Introduction
The strongest anthropogenic forcing on the Earth system during the last century has been a rise in atmospheric CO 2 concentrations due to anthropogenic CO 2 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 2 . Global warming also leads to changes in ocean circulation; for instance, shifting wind patterns might change the Southern Ocean upwelling and increase natural CO 2 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 2 decreases with increasing CO 2 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 ter-Published by Copernicus Publications on behalf of the European Geosciences Union.

5988
C.-T. Chien: Presentation and evaluation of FOCI-MOPS v1 restrial 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 2 reasonably well, there are still discrepancies between the model results and data products concerning oceanic CO 2 sink estimates . 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 2 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 3 ), phosphate (PO 4 ), and oxygen (O 2 ) cycles. MOPS features a smaller number of prognostic variables than other Coupled Model Intercomparison Project Phase 6 (CMIP6) models ) (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 2 fixation and denitrification (Kriest and Oschlies, 2015). Because circulation patterns used in earlier calibration exercises of MOPS  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 spinup simulation, for historical simulations covering the period from 1850 to 2014, and control simulations with preindustrial conditions. We also discuss the variability among ensemble members of each set-up and the differences between CO 2 -concentration-driven and CO 2 -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 1/2 • 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 (1/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. 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 3 : 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 2 at the sea surface is calculated according to Orr et al. (2017). Altogether, the ocean biogeochemistry is simulated 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 3 and between ALK and PO 4 show that ALK is affected by changes in NO 3 and PO 4 . CaCO 3 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. via nine prognostic tracers and five chemical elements (PO 4 (P), NO 3 (N), O 2 (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 4 , NO 3 , 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 3 , 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 3 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 3 and PO 4 .
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 4 and NO 3 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 3 sediments to reach a steady state.
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.

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 "physicsonly" FOCI spin-up under year 1850 climate conditions (e.g. solar radiation, greenhouse gases, atmospheric nitrogen de-

5990
C.-T. Chien: Presentation and evaluation of FOCI-MOPS v1 position, 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 . For the 500-year FOCI-MOPS spin-up, phosphate (PO 4 ), nitrate (NO 3 ), and oxygen (O 2 ) are initialized using the WOA2013 data set (Garcia et al., 2013a, b), and preindustrial DIC and ALK are taken from GLODAPv2.2016b . 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 2 concentrations. We also carried out a set of experiments where the model was forced with CO 2 emissions and the atmospheric CO 2 was calculated prognostically (instead of prescribing atmospheric CO 2 concentrations). For those experiments, we use the 480th year of the FOCI-MOPS spin-up as a starting point for a 250-year CO 2 zero-emissiondriven 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.
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 Olsen et al., 2016). For O 2 , NO 3 , PO 4 , and ALK, model outputs are averaged over 1972 to 2013, and modelled year 2002 is used for DIC. For organic tracers, 10year-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. Over the first 100 years of the FOCI-MOPS 500-year spinup (spinup), a small (< 1.5 %) decrease in inorganic nutrients (PO 4 and NO 3 ) 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 3 because the marine nitrogen cycle did not reach a steady state, with N 2 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 2 fixation and denitrification can take thousands of years (Falkowski, 1997;Oschlies et al., 2019). Owing to a net NO 3 loss due to a higher global denitrification than N 2 fixation rate, together with the build-up of organic matter in the beginning, global-average NO 3 is about 1.1 mmol N m −3 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 2 concentrations continuously decrease throughout the spinup and in the end are close to a steady state. The loss of O 2 indicates that in the model the supply of O 2 via mixing is too slow and/or the consumption due to remineralization is overestimated. The changes in O 2 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 2 concentration, allowing an increasing fraction of organic material to remineralize below 2000 m. After 500 years, global-average O 2 is about 8.4 mmol O 2 m −3 (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 2 fluxes. Drift in the DIC inventory during the last 100 years in the spinup is −0.086 Pg C yr −1 , which meets the "acceptably small drift" (±0.1 Pg C yr −1 ) suggested in Jones et al. (2016). Global CaCO 3 production equals dissolution as there is no dynamic CaCO 3 pool in the model. Therefore, the alkalinity inventory is only affected by the changes in PO 4 and NO 3 . Since PO 4 and NO 3 decrease during the formation of organic matter and the nitrate reservoir additionally decreases resulting from the imbalance of denitrification and N 2 fixation (plus nitrification), the ALK inventory increases during the spinup in the stoichiometric ratio of 0.9914 : −1 ALK : NO 3 (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.

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 Sup-  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 2 -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 3 , O 2 , 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 4 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 117 : 16 : 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. Globalaverage NO 3 in 2014 is about 0.017 mmol N m −3 (0.06 %) higher than in 1850. One reason for the increase in NO 3 is a decrease of denitrification. The reduction of DOM and organic particle pools also contributes to the NO 3 increase (0.006 mmol N m −3 ). The changes in NO 3 are also reflected in a decreasing ALK inventory. The simulated relative decrease of the O 2 inventory is higher than that of PO 4 and NO 3 . In the model, the ocean loses about 1.5 mmol O 2 m −3 (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 2 content. The small decrease in export production would lead to reduced O 2 consumption via respiration and therefore elevate rather than reduce the oxygen inventory. Globally averaged DIC increases by 8 mmol C m −3 (0.35 %), mainly due to increased air-sea CO 2 fluxes into the ocean under rising atmospheric CO 2 , i.e. the uptake of anthropogenic CO 2 .
3.2 Biogeochemical model performance

Spatial distribution of inorganic tracers
The surface concentration of O 2 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 2 concentration is below 36 mmol O 2 m −3 . The low-oxygen waters develop where the physical supply of O 2 is sluggish and O 2 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 2 is biased low, particularly in the Southern Hemisphere (Figs. 4i-k and S5). Another region with clear negative O 2 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 2 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 2 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 2 . Surface phosphate (PO 4 ) 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 4 shows positive biases at around 100-1000 m and below 3000 m and negative biases between 1000 and 3000 m ( Fig. 6i-k).
In order to better understand the distribution of PO 4 , we adopt the methodology of Duteil et al. (2012) and calculate regenerated (PO 4reg ) and preformed PO 4 (PO 4pre ) assuming oxygen saturation (O 2sat ) at isopycnal outcrops (e.g. PO 4reg = AOU/r−O 2 : PO 4 , AOU = O 2sat − O 2obs , PO 4pre = PO 4 −PO 4reg ). For observations, we apply the O 2 : PO 4 ratio of 170 mol O 2 : mol P estimated in Anderson and Sarmiento (1994). For model output, we apply the ratio of 165.08044 mol O 2 : mol P used in the model.
Except for the deep Atlantic (deeper than 3000 m), preformed PO 4 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 4 by phytoplankton in iron-limited regions, such as the Southern Ocean, and lead to too low preformed PO 4 . The overestimated regenerated PO 4 in the deep Pacific > 3000 m (Fig. 8k) contributes to the overall high bias of PO 4 (Fig. 6k) and is consistent with the low O 2 biases in the same region (Figs. 4k and S5). The high bias of PO 4 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 4 and apparent oxygen utilization (Fig. 5), as the PO 4 is also overestimated in the deep Atlantic Ocean. Biases of climatological nitrate NO 3 overall show a similar distribution as for phosphate (Fig. 9a). Simulated NO 3 differs from the pattern of PO 4 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 3 is removed by denitrification. Substantial rates of denitrification can be realized (depending on the supply of organic matter) when O 2 concentrations fall below a few mmol m −3 (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 −1 , higher than observation-based estimates (ranging from 0.02 to 0.12 Pg N yr −1 ; 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 4 ( Fig. 6i-k).
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 4 and NO 3 . In the ocean, potential differences in CO 2 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 3 , via zooplankton egestion and plankton mortality in the surface layer and dissolution of CaCO 3 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 efolding export and dissolution of CaCO 3 (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 4 and NO 3 , as well as the negative bias in O 2 .
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.

Spatial distribution of organic tracers in the surface
In FOCI-MOPS, the growth of phytoplankton is determined by temperature, light, and ambient NO 3 and PO 4 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 −3 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 1/2 • 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 4 in the equatorial regions. 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 largescale 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).

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 −3 , and surface DIC averaged over 1870 to 1899 (pre-industrial DIC) is 1985 ± 0 mmol m −3 . They are lower than the values among CMIP5 models (surface ALK and pre-industrial DIC are 2365-2500 and 2050-2170 mmol m −3 , respectively; Oka, 2020) but are only slightly lower than the observations (surface ALK and pre-industrial DIC are 2362 and 2019 mmol m −3 , respectively; Lauvset et al., 2016; Olsen  (Table 2) with respect to dissolved inorganic tracers. 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.).
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 1 − HD 2 , 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 5998 C.-T. Chien: Presentation and evaluation of FOCI-MOPS v1 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.

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 −1 ), export production, the sinking of detritus at 100 m (6.8 Pg C yr −1 ), the flux of detritus at 2000 m (0.35 Pg C yr −1 ), burial of detritus that sank to the bottom of the ocean (0.34 Pg C yr −1 ), and N 2 fixation (0.12 Pg N yr −1 ) are within the range of observation-based estimates. CaCO 3 production (0.77 Pg C yr −1 ) is lower than previous observation-based estimates (0.8-4.7 Pg C yr −1 ). Water column denitrification (0.13 Pg N yr −1 ) is slightly higher than observational estimates. The highest annual primary production occurs in the eastern equatorial Pacific Ocean and some coastal regions (> 1000 mg C m −2 d −1 , Fig. 17a). In the subtropical oligotrophic gyres, it is typically lower than 200 mg C m −2 d −1 . The spatial distributions of export production at 100 m and CaCO 3 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 2 levels between the surface and 2000 m depth mentioned earlier. Regions of elevated flux of detritus largely overlap with O 2 concentrations below 36 mmol O 2 m −3 , which is the threshold for anaerobic remineralization in the model. Denitrification is shaped by low O 2 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  Olsen et al., 2016). The order of panels is the same as in Fig. 4. 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 2 fixation occurs mostly between 40 • N and 40 • S in the model. In the model, 60 % of N 2 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 18 shows the ocean pCO 2 and the air-sea CO 2 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 2 and air-sea CO 2 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 2 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 4 and NO 3 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 2 and CO 2 flux. Overestimated SST can be associated with positive pCO 2 and negative airsea CO 2 flux biases. Those biological and physical factors result in complex patterns in the pCO 2 and CO 2 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)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). Therefore, the mismatches in pCO 2 and CO 2 uptake can be from two components: the mismatch inherited from the spinup 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 2 fluxes and uptake agrees with observations, and the biases in air-sea CO 2 fluxes are in the range of other Earth system models participating in the CMIP5 intercomparison (Dong et al., 2016). When zon- ally averaged, much of the biases are averaged out, and both ocean pCO 2 and air-sea CO 2 fluxes in FOCI-MOPS match the observations rather well ( Fig. 18d and h).

Air-sea exchange of carbon
At high latitudes, the variation in the ocean pCO 2 and the air-sea CO 2 fluxes are usually higher (Fig. S15) due to the physical mixing being stronger and having higher uncertainty.
3.4 Response to increasing atmospheric CO 2 3.4.1 Atmospheric CO 2 , air-sea CO 2 flux, surface air temperature anomaly, and ocean heat content anomaly Atmospheric CO 2 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 2 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 2 source to the atmosphere. The mean global oceanic CO 2 loss in piControl and ESM-piControl during the simulation period are 0.07 ± 0.01 and 0.05 ± 0.01 Pg C yr −1 , respectively. In the historical simulations, the globally integrated air-sea CO 2 flux gradually increases between 1850 and 1960 along with increasing anthropogenic CO 2 perturbation. After 1960, the rate of increase accelerates, and the flux reaches 2.11 ± 0.10 and 2.21 ± 0.14 Pg C yr −1 over 2005-2014 period in Hist and ESM-Hist, respectively, both are slightly lower than the 2.34 Pg C yr −1 estimated by the Global Carbon Project (Friedlingstein et al., 2020; see the summary tab; Fig. 20b). The cumulative air-sea CO 2 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  Moriarty and O'Brien, 2013), and POM , 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. 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 2 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 2 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 2 fluxes and carbon storage in FOCI-MOPS are consistent with those of CMIP6 models (Fig. 21). The highest CO 2 uptake between 40 and 65 • S is associated with the wind-driven upwelling, which brings deep water with low anthropogenic CO 2 to the surface and is able to uptake a higher amount of CO 2 (Frölicher et al., 2015). The stronger upwelling seems to be related to the higher variation in CO 2 uptake around the similar latitudes within the ensemble and among the CMIP6 models. The anthropogenic CO 2 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 Table 2. Correlation coefficient r, normalized standard deviation σ M /σ O , global bias normalized by observed mean (in squared brackets the non-normalized bias is shown, given in µg Chl L −1 for phytoplankton and in mmol m −3 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 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. 45-65 • S is because of the underestimated cumulative compatible emission in the Hist. Nevertheless, the cumulative air-sea CO 2 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).

Metric
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 2 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 2 fluxes, temperature, and ocean heat content.   Landolfi et al., 2016;Yoshimura et al., 2007) and unpublished (Landolfi, unpublished)

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 2 ) and nitrate (NO 3 ) levels, and primary production (PP) . 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 sur-  face 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 2 flux (Fig. 20b) Schmidtko et al., 2017), such an underestimation of observationally estimated changes in the marine O 2 inventory is a common feature among Earth system models (Oschlies et al., 2018). The NO 3 concentrations in the euphotic zone (0-100 m) are negatively correlated with SST (Fig. 22d), suggesting a reduced surface supply of NO 3 with increasing ocean stratification, a relationship that also exists in other global models Table 3. Primary production, export production, flux of detritus at 2000 m, global burial, CaCO 3 production (Pg C yr −1 ), N 2 fixation, and water column denitrification (Pg N, yr −1 ) in Hist simulations ( (Fig. 22e).

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 4 , NO 3 , O 2 , 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 his-torical (Hist) experiments were performed with prescribed atmospheric CO 2 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 −1 ) in the Cou-  pled 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 2 and air-sea CO 2 fluxes agree well with observations (Fig. 18). Simulated changes due to the increasing  (Frölicher et al., 2015;Terhaar et al., 2021). atmospheric CO 2 in globally integrated air-sea CO 2 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 2 emission-driven experiments, including ESM-spinup, ESM-piControl, and ESM-Hist. The air-sea CO 2 flux in these simulations is slightly higher in the historical emission-driven simulation (ESM-Hist) than in the atmospheric-CO 2 -driven simulation (Hist), resulting in a higher cumulative air-sea CO 2 flux of 139 Pg C relative to 129 Pg C in the Hist between 1850 and 2014. The difference in CO 2 perturbation forcing contributes to a higher atmospheric CO 2 until the end of the simulation period in ESM-Hist than in the Hist simulations. This higher atmospheric CO 2 likely contributes to the higher cumulative airsea CO 2 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 ± standard deviation of CMIP6 models, respectively, as described in Kwiatkowski et al. (2020). 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 2 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 , for instance it does not explicitly resolve silicate or iron cycles and does not consider the CaCO 3 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 4 ). Additionally, nitrate (NO 3 ) and oxygen (O 2 ) are simulated, with biogeochemical interactions among the different elements coupled via fixed stoichiometric ratios. As the model also simulates NO 3 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 2 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.1 Plankton dynamics
Phytoplankton (PHY) growth depends on temperature T (f 1 (T )), daylight intensity (I , (W m −2 d −1 )), day length τ (d) (f 2 (I, τ )), and nutrients (f 3 (PO 4 , NO 3 )). Temperaturedependent growth is formulated following Eppley (1972), in the notation by Schmittner et al. (2008): where µ PHY = 0.6 (d −1 ) 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 −2 ) −1 d −1 ), we here calculate the light limitation based on the half-saturation constant for light, I c = 9.653 (W m −2 ), as expressed through I c = µ PHY /α (thereby resulting in α = 0.0622 ((W m −2 ) −1 d −1 )): where I (W m −2 d −1 ) 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 −1 ) and k c = 0.48 ((mmol P m −3 ) −1 m −1 ) are the attenuation coefficients of water and phytoplankton, respectively. The function φ is given by (Evans and Garçon, 1997) The dependence of growth on nutrients f 3 (PO 4 , NO 3 ) 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): where k PHY = 0.031 (mmol P m −3 ) is the half-saturation constant for phosphate . Total phytoplankton growth PP (mmol P m −3 d −1 ) 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 > P * = 10 −6 (mmol P m −3 )): Phytoplankton experiences a linear loss given by λ PHY = 0.03 (d −1 ). It is grazed by zooplankton (ZOO), where grazing G described by a Holling-III function, with a maximum grazing rate µ ZOO = 1.893 (d −1 ) and half-saturation constant K ZOO = 0.086 (mmol P m −3 ): Only a fraction ZOO = 0.75 of grazing G (in (mmol P m −3 d −1 )) is effectively ingested, the rest is released again via egestion. Zooplankton further experiences a quadratic mortality κ ZOO = 4.548 ((mmol P m −3 ) −1 d −1 ) 6010 C.-T. Chien: Presentation and evaluation of FOCI-MOPS v1 and a linear excretion rate given by λ ZOO = 0.03 (d −1 ). Phytoplankton and zooplankton further die with a constant mortality rate of λ PHY = λ ZOO = 0.01 (d −1 ) but only when present above the lower concentration threshold P * = 10 −6 (mmol P m −3 ). Therefore, the source-minus-sink terms due to biogeochemical interactions for phytoplankton (S Bio PHY ) and zooplankton (S Bio ZOO ) are 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 117 : 16 : 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 λ DOP = 0.17 (yr −1 ) but only when present above the lower limit of P * = 10 −6 (mmol P m −3 ). Its oxic and suboxic remineralization further depends on the availability of oxidants oxygen and nitrate, as described by the terms S Rox DOP and S Rsubox DOP , which are described in detail in Sect. A1.3 below. Thus, the sourceminus-sink term for DOP, S Bio DOP , is Like DOP, detritus (DET) is remineralized with a fixed rate λ DET = 0.05 (d −1 ) to nutrients when present above the lower limit of P * = 10 −6 (mmol P m −3 ), and its decomposition depends on oxygen and nitrate, as described for the terms S Rox DET and S Rsubox DET 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(z) = w 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 = λ DET /w 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 lin = λ DET /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-minussink term of detritus, S Bio DET , is therefore

A1.3 Oxic and suboxic remineralization
If oxygen is above a threshold defined by , organic matter is remineralized aerobically according to a sigmoidal function: where K O 2 = 1.066 (mmol O 2 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 O 2 : where λ DOP and λ DET are the remineralization rates of DOP and detritus, respectively, and t is the time step length of the biogeochemical model. R −O 2 :P = 165.08044 denotes the stoichiometric oxygen demand of aerobic remineralization. The aerobic decay rate limitation is then Therefore, the sinks of oxygen due to remineralization of DOP and detritus are defined by If O 2 * is lower than 36 mmol O 2 m −3 , denitrification additionally sets in. As for oxygen, we first define a quadratic rate limitation of this process based on a minimum concentration of nitrate, NO 3 min = 15.978 (mmol N m −3 ), with NO 3 * = max(0, NO 3 −NO 3 min ). To account for inhibition of denitrification by oxygen we further reduce this rate by the inverse oxygen consumption rate: where K NO 3 = 23.104 (mmol N m −3 ) is the half-saturation constant for the denitrifiers' uptake of nitrate. As for oxygen, we restrict the use of nitrate to the amount available: with R −NO 3 :P = 0.8 R −O 2 :P − d = 116.064352 (mmol NO 3 : mmol P),following the stoichiometry of Paulmier et al. (2009). The rate limitation of anaerobic decay is then Therefore, the sinks of oxygen and nitrate due to denitrification of DOP and detritus are defined by

A1.4 Nitrogen fixation
As in Kriest and Oschlies (2015), nitrogen fixation depends on temperature and nutrient ratio: with µ NFix = 1.88924 (µmol N m −3 d −1 ) being the maximum nitrogen fixation of the (implicit) cyanobacteria and t 2 , t 1 , t 0 , 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 4 > 10 −6 mmol P m −3 .

A1.5 Nutrients and oxygen
Phosphate (PO 4 ) 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 Bio PO 4 , of The loss and gain of nitrate (NO 3 ) 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): Finally, oxygen (O 2 ) increases due to photosynthesis and decreases because of aerobic remineralization and respiration by zooplankton.
In addition, oxygen exchanges with the atmosphere at the sea surface (i.e. for layer 1) are given as S Air O 2 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: Changes in phosphate and nitrate further affect alkalinity via

A2.2 Calcite production and dissolution
We assume that a constant fraction p CaCO 3 of detritus production via zooplankton egestion and plankton mortality to detritus is in the form of calcite.
where a = 117 (mol C:mol P) is the molar ratio of C : P in organic matter, and σ CaCO 3 = 0.032 (mol CaCO 3 :mol C) is the calcite-to-organic-carbon ratio. Calcite production reduces DIC by 1 and alkalinity by 2: Following Schmittner et al. (2008), we integrate the production of calcite over the entire water column: 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 CaCO 3 ), with l CaCO 3 = 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 2 and solution of the carbonate system The simulation of the carbonate chemistry system and of the air-sea gas exchange of CO 2 in MOPS follows an OCMIPtype formulation (Orr et al., 1999) with updates from Orr et al. (2017) for the air-sea gas exchange of CO 2 . The air-sea CO 2 flux is computed as follows: where k w (in m s −1 ) is the gas transfer velocity, (CO * 2 ) sat (in mol kg −1 ) is the saturation concentration of CO 2 , and (CO * 2 ) (in mol kg −1 ) is the surface ocean dissolved CO 2 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: 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 2 in equilibrium with the water-vapour-saturated atmosphere at a total atmospheric pressure of 1 atm (i.e. (CO * 2 ) sat ) is computed as follows: where xCO 2 is the mole fraction of CO 2 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 2 reacts with seawater, forming carbonic acid (H 2 CO 3 ), most of which dissociates into two other inorganic species, bicarbonate (HCO − 3 ) and carbonate (CO 2− 3 ) ions. CO * 2 refers to the sum of dissolved CO 2 and the much less abundant H 2 CO 3 . The sum of the three species CO * 2 + HCO − 3 + CO 2− 3 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 * 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 noniterative 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 threedimensional 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 * 2 (Eq. 8 in Follows et al., 2006) and surface CO 2 fugacity as follows: where K 0 is the solubility coefficient of CO 2 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 BUR DET (mmol P m −2 d −1 ), depends on the rain rate of detritus to the sea floor (F B , (mmol P m −2 d −1 )) via 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 = A F BUR da is the global burial in a time step integrated over the sea floor area A and V 1 is the global volume of the surface layer, 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: 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: 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:

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 . In particular, Kriest et al. (2020) optimized the parameters for oxidant-dependent remineralization NO min 3 , K NO 3 , K O 2 , the maximum nitrogen fixation rate µ NFix , the oxygen requirement for aerobic remineralization, R −O 2 :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 −O 2 :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 −O 2 :P (165.08044 instead of 151.1 mol O 2 : 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 −3 d −1 ; see also Table A1).
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 σ CaCO 3 and l CaCo 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 σ CaCo 3 = 0.032 mol CaCO 3 : mol C org (instead of 0.035 mol CaCO 3 : mol C org ) and l CaCo 3 = 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 1 • × 1 • 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.  Kriest et al., 2020) 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 −1 to mmol m −3 using in situ density com-puted 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).

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 1 • × 1 • . 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 −3 , respectively, an unweighted mean of 0.0161 mmol P m −3 , and a standard deviation of 0.0122 mmol P m −3 (see Table B1).

B2.3 Zooplankton
For the evaluation of simulated zooplankton we use the MAREDAT data set of mesozooplankton . This sparse, quasi-climatological data set contains 42.245 data points of monthly mean mesozooplankton (in mg C m −3 ) on a 1 • × 1 • 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., 2017Kiko et al., , 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 −3 , an average of 0.0062 mmol P m −3 , and a standard deviation of 0.0101 mmol P m −3 (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 1 • × 1 • 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 −3 , respectively, an average of 0.0266 mmol P m −3 , and a standard deviation of 0.0499 mmol P m −3 . 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.   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 1 • × 1 • 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 −3 , an average of 0.178 mmol P m −3 , and a standard deviation of 0.21 mmol P m −3 . 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 differ-ent 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 −3 ; nitrate: 0-50 mmol N m −3 ; oxygen: 0-400 mmol O 2 m −3 ; DIC and alkalinity: 1700-2500 mmol m −3 ; phytoplankton and zooplankton: 0-0.1 mmol P m −3 ; POP: 0-0.5 mmol P m −3 ; DOP: 0-1 mmol P m −3 ). Examples for different distributions are given in Figs. B1 and 16. Following binning, BD was evaluated as 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(1 − HD 2 ) as follows: In contrast to the BD, the HD is bounded by 0 ≤ HD ≤ 1. Finally, the L1 norm of distributions as given by the following equation is evaluated (note that this norm is bounded by 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). Code and data availability. The FOCI-MOPS is build based on a published FOCI version (Matthes et al., 2020). All modifications to the published version code and full runtime environment are provided at https://doi.org/10.5281/zenodo.6772175 (Chien et al., 2022). Model data and codes necessary to reproduce the figures present here are available at the same location. All data used in this paper are publicly accessible, including World Ocean Atlas 2013 (https://www.ncei.noaa. gov/products/ocean-climate-laboratory, last access: 15 August 2020, Garcia et al., 2013a, b), Global Ocean Data Analysis Project version 2 mapped climatology (GLODAPv2.2016b, https://www.glodap.info/index.php/mapped-data-product/, last access: 12 May 2016, Lauvset et al., 2016;Olsen et al., 2016), surface pCO 2 and air-sea CO 2 fluxes (https://www.ncei.noaa.gov/access/ ocean-carbon-acidification-