Calibration of a simple and a complex model of global marine biogeochemistry

The assessment of the ocean biota’s role in climate climate change is often carried out with global biogeochemical ocean models that contain many components, and involve a high level of parametric uncertainty. Examination the models’ fit to climatologies of inorganic tracers, after the models have been spun up to steady state, is a common, but computationally expensive procedure to assess model performance and reliability. Using new tools that have become available for global model 5 assessment and calibration in steady state, this paper examines two different model types a complex seven-component model (MOPS), and a very simple two-component model (RetroMOPS) for their fit to dissolved quantities. Before comparing the models, a subset of their biogeochemical parameters has been optimised against annual mean nutrients and oxygen. Both model types fit the observations almost equally well. The simple model, which contains only nutrients and dissolved organic phosphorus (DOP), is sensitive to the parameterisation of DOP production and decay. The spatio-temporal decoupling of 10 nitrogen and oxygen, and processes involved in their uptake and release, renders oxygen and nitrate valuable tracers for model calibration. In addition, the non-conservative nature of these tracers (with respect to their upper boundary condition) introduces the global bias as a useful additional constraint on model parameters. Dissolved organic phosphorous at the surface behaves antagonistically to phosphate, and suggests that observations of this tracer although difficult to measure may be an important asset for model calibration. 15

be required for a fair assessment of the virtues of models of different complexity. Such a scan usually requires many model evaluations, which, given long equilibration time scales of coupled global models (Khatiwala, 2008;Wunsch and Heimbach, 2008;Primeau and Deleersnijder, 2009;Siberlin and Wunsch, 2011), is difficult to carry out. For assessment of only surface properties and processes a short model spinup may be sufficient; however, on a global scale, many centuries to millennia of coupled model simulations are necessary, in order to remove the drift in biogeochemical tracer fields and fit to observed 5 properties (Kriest and Oschlies, 2015;Seferian et al., 2016).
Only recently tools have become available that allow for a reduction in simulation times, such as the Transport Matrix Method (TMM, Khatiwala et al., 2005), Newton's method, which requires the inverse of the Jacobian (Kwon and Primeau, 2006Primeau, , 2008, matrix-free Newton-Krylov (MFNK, Khatiwala, 2008;Li and Primeau, 2008), or surrogate-based optimisation (Priess et al., 2013). The gain in computational efficiency resulting from these methods can then be used for a systematic 10 calibration of global biogeochemical models. For example, Kwon andPrimeau (2006, 2008) used global climatological data sets of phosphate, inorganic carbon, and alkalinity to calibrate a simple global biogeochemical model. The misfit between observed and simulated phosphate was used by DeVries et al. (2014) to calibrate parameters related to particle properties in a simple two-component, nutrient-restoring model. In a similar approach Holzer et al. (2014) optimised parameters for opal production and dissolution against observed silicate. Letscher et al. (2015) switched between a complex and a simple model of 15 ocean biogeochemistry to estimate production and decay rates of dissolved organic phosphorus on a global scale.
All these biogeochemical models employed in global parameter estimates were of a low level of biogeochemical complexity.
One reason for this restriction might be the large number of tracers in more complex models, which increases simulation time.
Another problem is associated with the variety of time scales associated with more complex models. Piwonski and Slawig (2016) used MFNK to evaluate the steady state of simple and complex biogeochemical models. They noted that " [...] for more 20 complex models the Newton method requires more attention to solver parameter settings [...]" (Piwonski and Slawig, 2016), which may be related to the highly nonlinear structure of these models. The nonlinearity, and the large number of parameters, also complicates their simultaneous optimisation (Ward et al., 2010). On a global scale, these problems are amplified by the sparsity of observations of organism groups, particularly of higher trophic levels. Observations of dissolved inorganic constituents, on the other hand, are much more frequent, and therefore provide more information on the spatio-temporal variability 25 of these tracers.
Recently, Kriest et al. (2017) combined the TMM with an estimation of distribution algorithm (Covariance Matrix Adaption Evolution Strategy, CMA-ES), to optimize six biogeochemical parameters of a seven component model against global climatologies of annual mean phosphate, nitrate and oxygen. They showed that annual mean tracer concentrations do not provide much information on parameters related to the dynamic biological processes taking place in the euphotic zone, but that 30 parameters related to larger time-and space scales could be estimated from these observations. To follow up on this, I here investigate if parameters related to oxidant-dependent decay in the mesopelagial are better constrained by this type of misfit function. This is done by replacing four parameters of the optimisation carried out by Kriest et al. (2017) by parameters related oxidant-affinity of remineralisation, and -to account for the possible alterations in fixed nitrogen turnover -by the maximum nitrogen fixation rate. Given the successful parameter optimisation of simpler models noted above, and also to acknowledge the fact, that these models have been popular and quite successful in global simulations of ocean biogeochemistry, this paper further presents an optimised model, which has been derived from downscaling the seven-component model MOPS (Kriest and Oschlies, 2015;Kriest et al., 2017) to a model that retains only three abiotic dissolved tracers (phosphate, nitrate, and oxygen) and one biotic tracer (dissolved organic phosphorus; DOP). This new model, which I refer to as "RetroMOPS", includes the oxidant-5 dependency of MOPS, but is otherwise very similar to models applied earlier in global models. In contrast to some of these models (Marchal et al., 1998;Najjar et al., 2007) it assumes no relaxation to observed tracer fields, but simulates changes in surface production fully prognostic, as in Bacastow and Maier-Reimer (1991);Maier-Reimer (1993); Matear and Hirst (2003); Parekh et al. (2005).
After a brief presentation of model MOPS (Kriest and Oschlies, 2015), the downscaled model RetroMOPS is introduced, 10 followed by an outline on circulation, optimisation and experimental design (section 2). In section 3 results from optimisation of both MOPS and RetroMOPS are presented and discussed. The paper closes with conclusions drawn from these experiments. plankton, zooplankton, dissolved organic phosphorus (DOP) and detritus are calculated in units of mmol P m 3 . Oxygen is coupled to the P-cycle with a constant stoichiometry given by R O2:P . Aerobic remineralisation of organic matter follows a saturation curve, with half-saturation constant K O2 . Aerobic remineralisation ceases when oxygen declines; at the same time, denitrification takes over, as long as nitrate is available above a defined threshold, DIN min . Like the oxic process, suboxic remineralisation follows a saturation curve for oxidant nitrate, with half-saturation constant K DIN . MOPS does not explicitly 20 resolve the different oxidation states of inorganic nitrogen (nitrite, N 2 O, ammonium), but assumes immediate coupling of the the different processes involved in nitrate reduction, the end-product being dinitrogen (see also Paulmier et al., 2009;Kriest and Oschlies, 2015). All organic components are characterised by a constant N:P stoichiometry of d = 16. Loss of fixed nitrogen is balanced by a simple parameterisation of nitrogen fixation by cyanobacteria, which relaxes the nitrate-to-phosphate ratio to d with a time constant, µ ⇤ NFix . Detritus sinks with a vertically increasing sinking speed: w = a z. Assuming a constant degrada-25 tion rate r, in equilibrium this would result in a particle flux curve given by F (z) / z b , with b = r/a. For better comparison with values of b derived from observations (e.g. Martin et al., 1987;Van Mooy et al., 2002;Buesseler et al., 2007), and with the simpler model RetroMOPS (see below), a is expressed in terms of b (assuming constant, nominal r = 0.05 d 1 ). A fraction of detritus deposited at the sea floor (at the bottom of the deepest vertical box) is buried instantaneously in some hypothetical sediment. Non-buried detritus is resuspended into the water column, where it is treated as regular detritus. The phosphorus 30 budget is closed on an annual time scale through resupply via river runoff. More details about the biogeochemical model and parameters can be found in Kriest and Oschlies (2013) and Kriest and Oschlies (2015).

Model RetroMOPS
MOPS' structure has been simplified by skipping the explicit simulation of phytoplankton, zooplankton, and detritus. The remaining equations of, and functional relationships between, phosphate, nitrate, oxygen and DOP have been parameterised similar to MOPS. Because the downscaled model resembles so many features of earlier biogeochemical models simulated in a global context (e.g., Bacastow and Maier-Reimer, 1991;Maier-Reimer, 1993;Matear and Hirst, 2003), but keeps the oxidant 5 dependency of MOPS, the model is named "RetroMOPS".

Primary production
Like MOPS, RetroMOPS calculates primary production only in the euphotic zone, which, in the current configuration, is confined to the upper two numerical layers (k EZ = 2, z = 0 120 m). Phytoplankton is parameterised with a constant concentration of P HY = 0.02 mmol P m 3 , which is the mean phytoplankton concentration in the upper 120 m of two optimised model 10 setups MOPS oS and MOPS oD (see below). Using this constant phytoplankton concentration RetroMOPS calculates light-and nutrient dependent primary production P in each layer k as: where f (I(k)) defines light-limitation, µ PHY is the temperature-dependent maximum growth rate of phytoplankton, and L determines the limiting nutrient: L(k) = min(PO 4 (k), DIN(k)/d) (see Kriest and Oschlies, 2015, for more details). 15 2.2.2 The fate of primary production: Export, DOP production and remineralisation Instead of resolving heterotrophic processes (zooplankton grazing, excretion and egestion) at the sea surface explicitly, in RetroMOPS a fraction DOP of organic matter fixed photosynthetically is immediately released as dissolved organic phosphorous, DOP. DOP then decays to phosphate and nitrate with a constant rate DOP . To allow for a potential, fast recycling loop at the surface, RetroMOPS parameterises an additional decay rate, sDOP , that affects DOP only in the first two layers. 20 The remaining fraction of production, 1 DOP , of each layer in the euphotic zone is exported to the layers below, where it immediately remineralises to nutrients, following a power-law of depth. The discretised form for flux F into box j from all (surface) source layers k, with 1  k  k EZ is then given by: where z(k) denotes the thickness of a numerical (source) layer, and z(j) is the depth of the upper boundary of layer j.
The flux divergence, D = dF/dz, for any box j in discretised form is defined by Neglecting oxidant dependency of decay, the entire flux divergence D(j) would be released as phosphate and nitrate, with equivalent oxidant consumption. It is, however, possible that oxidants become depleted at some location. Earlier models in 5 this case continued the degradation of organic matter, thereby implicitly assuming unspecified oxidants (e.g., Marchal et al., 1998;Matear and Hirst, 2003;Najjar et al., 2007;Kriest et al., 2010Kriest et al., , 2012. In contrast, RetroMOPS, like MOPS, accounts for suppression of remineralisation (oxic/suboxic) in the absence of sufficient oxidants, by assuming saturation curves for the limitation by either oxygen or nitrate. The amount of organic matter available for oxidation is given by the decay of dissolved organic matter, DOP DOP , and by the flux divergence, D(j) (Eqn. 3). The discretised flux divergence, that can actually be 10 remineralised to phosphate and nitrate with available oxidants (oxygen and/or nitrate), D e↵ (j), is then determined by where s O2 (j) and s DIN (j) represent the oxidant limitation terms, as detailed in Kriest and Oschlies (2015, Equations 15-27).
The remaining flux divergence that cannot be remineralised under the given concentrations of oxidants is added as additional flux divergence to the layer below: where again D e↵ (j + 1) is evaluated. In the bottom layer the remaining flux that has not been remineralised in the water column eventually enters the sediment.

Benthic exchanges
Models that implicitly assume unspecified oxidants often prescribe a zero boundary flux, i.e. all organic matter in the last 20 bottom box is degraded instantaneously (e.g., Marchal et al., 1998;Matear and Hirst, 2003;Najjar et al., 2007;Yool et al., 2011). Both MOPS and RetroMOPS have to take "leftover" organic matter flux into account, that arrives undegraded at the sea floor because of incomplete remineralisation in the water column. The explicit detritus compartment in MOPS allows for only partial burial at the sea floor, which may result in detritus accumulation in the deepest model box (see Kriest and Oschlies, 2013). Because there is no such detritus compartment in RetroMOPS, all flux arriving at the sea floor is buried immediately. 25 Therefore, MOPS and RetroMOPS differ with respect to their lower boundary condition.

Nitrogen fixation
In RetroMOPS and MOPS nitrogen fixation balances the simulated loss of fixed nitrogen via denitrification. Both models do not explicitly simulate cyanobacteria, but assume zero net growth of these organisms, parameterised as an immediate release of fixed nitrogen as nitrate: f 1 parameterises the temperature dependence of nitrogen fixation with a second order polynomial approximation of the function by Breitbarth et al. (2007). f 2 regulates the relaxation of the nitrate:phosphate ratio towards the global observed stoichiometric ratio of d = 16. µ ⇤ NFix is the maximum nitrogen fixation of the parameterized cyanobacteria population (mmol N m 3 d 1 ; see Kriest and Oschlies, 2015, for more details).

Source-minus-sinks 10
Combining the above mentioned processes and interactions, the time rate of change for phosphate, nitrate, oxygen, and DOP due to biogeochemical processes are Summarising, RetroMOPS is similar to model "N-DOP" of Kriest et al. (2010Kriest et al. ( , 2012, to the phosphorus component of the model presented by Parekh et al. (2005), or to the models presented by Bacastow and Maier-Reimer (1991) and Maier-Reimer (1993), the exception being details of primary production at the sea surface, and the explicit parameterisation of oxidantdependent remineralisation. By assuming constant cyanobacteria biomass, and a relaxation of the nitrate:phosphate ratio via

Circulation and physical transport
All model simulations apply the Transport Matrix Method (TMM; Khatiwala, 2007, github.com/samarkhatiwala/tmm) for tracer transport, with monthly mean transport matrices (TMs) derived from a 2.8 global configuration of the MIT ocean 25 model, with 15 levels in the vertical (Marshall et al., 1997). Using this efficient offline approach, a time step length of 1/2 day for tracer transport and 1/16 day for biogeochemical interactions, simulation of 3000 years requires about 0.5-1.5 hrs on

Optimisation algorithm
Optimisation of parameters is carried out using an Estimation of Distribution Algorithm, namely the Covariance Matrix Adaption Evolution Strategy (CMA-ES; Hansen and Ostermeier, 2001;Hansen, 2006). The application of this algorithm to the 5 coupled biogeochemistry-TMM framework has shown good performance with respect to quality and efficiency (in terms of function evaluations), and is described only briefly below. More details about the algorithm, its setup and coupling to the global biogeochemical model can be found in Kriest et al. (2017).
Let n be the number of biogeochemical parameters to be estimated. In each iteration ("generation") the algorithm defines a population of individuals (biogeochemical parameter vectors of length n), with = 10 (derived from the default parameter 10 = 4 + 3 ln(n), Hansen and Ostermeier, 2001). The candidate vectors are sampled from a multi-variate normal-distribution, which generalizes the usual normal distribution, also known as Gaussian distribution, from R to the vector space R n .
Following the simulation of these individual model setups to steady state (3000 years), the misfit function is evaluated, and information of the current, as well as previous generations is used to update the probability distribution in R n such that the likelihood to sample good solutions increases. Usually, the realisation of the probability distribution update ensures that

Optimisation of MOPS
Building upon the optimisation of mostly surface-related parameters presented in Kriest et al. (2017, "OBS-NARROW", hereafter referred to as MOPS oS ), optimisation MOPS oD presented here aims at calibrating parameters related to processes that directly affect the oxidants nitrate and oxygen in subsurface layers. In MOPS oD the optimal parameters of MOPS oS for light and nutrient affinity of phytoplankton, zooplankton grazing and its mortality are retained, and parameters relevant for deep 5 aerobic and anaerobic remineralisation are subject to change during optimisation (Table 1).
Parameter K O2 determines the affinity of the aerobic remineralisation to oxygen, and the gradual transition from this process to denitrification (see Eqns. 15 and 20 of Kriest and Oschlies, 2015). K DIN determines the affinity of denitrification to nitrate.
Parameter DIN min defines the lower threshold for the onset of denitrification. MOPS oD also optimises the maximum rate of nitrogen fixation, µ ⇤ NF ix , which balances fixed nitrogen loss. The fifth and sixth parameter to be estimated are the oxygen 10 requirement per mole phosphorus remineralised, R O2:P , and the flux (or remineralisation) length scale, b.
To investigate the influence of observations entering the misfit function, MOPS oD is repeated with a reduced data set, that excludes the eastern equatorial Pacific (east of 140 W, between 10 S and 10 N) from the misfit function. This optimisation is named MOPS oD ⇤ . In the following, results from the optimised models MOPS oS , MOPS oD and MOPS oD ⇤ are compared to a reference experiment, MOPS r , which represents a "hand-tuned", a priori setup of this model.

Optimisation of RetroMOPS
In model RetroMOPS processes such as grazing of phytoplankton, and its subsequent release of organic or inorganic phosphorus are parameterised via a single component, DOP. Because DOP production and decay regulate the partitioning between sinking and dissolved organic matter, optimisation RetroMOPS o targets at these parameters, namely DOP , sDOP and DOP .
While DOP , as parameter that regulates the export ratio, may be more or less well constrained, sDOP and DOP both 20 include a variety of processes, which may act on time scales of days to years. In a set of nine a priori sensitivity experiments the effect of these parameters on the misfit function is explored by varying DOP between 0.18 y 1 and 0.72 y 1 , and sDOP between 0 y 1 and 0.36 y 1 (see table 2). The results of these sensitivity experiments provide a guidance for upper and lower boundaries of optimisation (table 1). The sensitivity experiment with the lowest misfit ( sDOP = 0, DOP = 0.36) is used for comparison with the optimised RetroMOPS, and referred to as RetroMOPS r . 25 The explicit representation of detritus in MOPS may result in considerable numerical diffusion (particularly on coarse vertical grids as used here; see also Kriest and Oschlies, 2011) and thus in a different estimate of optimal b then when applying a direct flux curve, such as in RetroMOPS. Therefore, b is included as fourth parameter to be optimised.  (Fig. S1), and the narrow, almost gaussian distribution of the best 10% to 1‰ of parameters (Fig. 1, Table 3). The good determination of b by dissolved inorganic tracers is agreement with earlier studies (Kwon and Primeau, 2006;Kriest 5 et al., 2017;Schartau et al., 2016). Parameters related to oxidant-dependent remineralisation approach the lower (K O2 ) or upper (K DIN , DIN min ) boundary, with a rather wide, skewed distribution. The rate for maximum nitrogen fixation shows a slightly skewed distribution, but suggests an overall good estimate of this parameter.
Fixed nitrogen loss and gain depend on parameters for oxidant-dependency of remineralisation: In MOPS oS , both fluxes are very high (Fig. 2), and outside the observed range (Table 4). Because optimisation MOPS oD results in a strongly reduced 10 affinity to, and higher threshold of, nitrate, its pelagic fixed nitrogen loss is almost halved, and now agrees with observed global estimates (Table 4). Further, as a result of reduced denitrification, the nitrate deficit in the eastern equatorial Pacific is smaller; however, at the cost of a small underestimate of observed oxygen in this region (Fig. 3). The latter is a consequence of the now very low half-saturation constant for oxygen uptake (Table 3). Overall, optimisation of parameters related to the oxidant affinity of oxic and suboxic remineralisation leads to a slightly improved fit to tracer concentrations, to J ⇤ = 98% of that of 15 MOPS oS (Table 3), and to a better agreement with observed estimates of global biogeochemical fluxes (Table 4).
Optimisation MOPS oD results in a high treshold for the limitation of denitrification, with K DIN and DIN min close to their upper boundaries. The increase protects nitrate from becoming depleted in the upwelling regions, particularly the eastern equatorial Pacific, and resembles results obtained by Moore and Doney (2007): To prevent their model from reproducing unrealistically low nitrate values in this region, they had to impose a threshold of 32 mmol NO 3 m 3 for the occurrence 20 of denitrification. An explanation for this requirement of a high nitrate threshold might be found in the representation of the equatorial intermediate current system in coarse resolution models, which can result in spurious zonal oxygen gradients (Dietze and Loeptien, 2013;Getzlaff and Dietze, 2013). It is possible that the optimisation of biogeochemical parameters attempts to ameliorate these effects, which are in fact caused by the parameterisation of physics.
To investigate the impact of this region on the parameter estimate, an additional optimisation was carried out, that applies the 25 same set of parameters to be optimised, but omits the eastern equatorial Pacific from the calculation of the misfit function. This optimisation MOPS oD ⇤ generates a lower threshold of nitrate for the onset of denitrification, and a higher maximum nitrogen fixation rate (Table 3), resulting in slightly enhanced fixed nitrogen turnover, particularly in the eastern equatorial Pacific (Fig. 2). Global fixed nitrogen loss increases by about 20%, towards the upper limit of observed estimates (Table 4). Compared to MOPS oD the estimates of K DIN and DIN min become more uncertain with respect to the best 10% to 1‰ individuals, and 30 show an even bimodal distribution (Fig. S2, Table 3). The uncertainty in parameter estimates can be related to the missing data in regions of simulated denitrification. Summarising, although the eastern equatorial Pacific, and potential unresolved processes in simulated circulation, evoke only relatively small effects on some parameter estimates, these nevertheless result in an increase in global fixed nitrogen loss of about 20%.

Sensitivity of RetroMOPS to DOP production and decay
In RetroMOPS fast DOP recycling leads to the higher primary production, export production, and deep organic particle flux, especially in the equatorial upwelling regions (Fig. 4). While this has only a small effect on vertically or globally averaged phosphate concentrations (Figures 5 and 6), it leads to a a large underestimate of nitrate in the ocean (Figures S4 and 6). The underestimate can be explained by the tight coupling between production, export and denitrification, which leads to higher 5 denitrification and global fixed N-loss (Fig. 4), and thus a larger nitrate deficit (Fig S4) in the eastern equatorial Pacific. This is in agreement with effects hypothesised and investigated by Landolfi et al. (2013).
In the model nitrogen fixation counteracts fixed nitrogen loss through denitrification. In contrast to nitrogen fixation, which is not much affected by DOP turnover rates, global fixed nitrogen loss increases with increasing DOP decay rate (Fig. 4).
The imbalance between nitrogen losses and gains suggests that the models even after 3000 years of simulation are not yet in The effect of DOP recycling on oxygen concentrations differs from its effect on nitrate. With fast recycling DOP is remineralised mostly at its place of production, and does not contribute much to oxygen consumption in deep waters (see also Fig S3).
As a consequence, deep oxygen concentrations are high, particularly in the northern North Pacific (Fig. 5), and global average 20 oxygen is overestimated by more than 10% (Fig. 6). Slow DOP recycling, in contrast, leads less organic matter remineralisation in preformed, well-ventilated waters, but more remineralisation in deep waters. This in turn results in an underestimate of global mean oxygen of almost 10% (for DOP = 0.18 y 1 and DOP = 0), which is somewhat surprising, given that production and export in this scenario are the lowest of all simulations (Fig. 4).
Overall, the best fit to observed inorganic tracer concentrations is achieved with moderate DOP recycling (Table 2, Fig. 5). 25 Most likely because of its fixed inventory, phosphate contributes to less than 1/3 of the misfit function, and is quite insensitive to changes in DOP recycling rate (Fig. 6). Nitrate and oxygen play a larger role for model fit, because their inventory can adapt to changing biogeochemistry. The misfit to nitrate and oxygen more or less increases in concert with their bias (Fig. 6).
Therefore, these tracers with their flexible inventory provide some very useful constraints on DOP recycling rates.
Slow DOP recycling increases DOP concentrations, particularly in the ACC and in the northern North Atlantic (Fig. 5), 30 and simulated concentrations largely exceed the observations (Yoshimura et al., 2007;Raimbault et al., 2008;Torres-Valdes et al., 2009;Letscher and Moore, 2015). Only the simulation with quite fast DOP recycling of DOP = 0.72 y 1 and sDOP = 0.36 y 1 results in reasonable concentrations of DOP -but at the cost of too high phosphate concentrations along these sections, and a too high global misfit (Table 2), a too low nitrate and too high oxygen inventory (Figures 5 and 6).

Optimisation of RetroMOPS
All four parameters of RetroMOPS o are quite well constrained by the observations, as indicated by the narrow, almost gaussian distribution around the optimal parameter (Figures 7, S5, and Table 3). Optimisation reduces the decay rate for surface DOP, sDOP , to almost zero, i.e., in RetroMOPS there seems to be no requirement for fast DOP turnover at the surface. The optimal DOP total remineralisation rate is about 0.5 y 1 , and the optimal fraction of primary production released as DOP is 73%, 5 resulting in a slightly higher turnover as compared to the reference scenario RetroMOPS r . Optimal DOP agrees very well with = 0.74 obtained by Kwon and Primeau (2006); however, their optimal DOP decay rate was twice as high (1 y 1 ).
When optimising a simple biogeochemical model similar to RetroMOPS against observed phosphate, Kwon and Primeau (2006) noted a correlation between DOP production fraction and decay rate, impeding the simultaneous estimation of these parameters. On the contrary, in optimisation RetroMOPS o both DOP and the DOP decay rates seem to be rather well con-10 strained. An analysis of the different components of the misfit function, similar to Fig. 4 of Kwon and Primeau (2006), helps to resolve this apparent contradiction. For this, in Fig. 8 misfit (J and J(j) of Eqn. 11) and bias of the best 5% of all individual is mapped against DOP and DOP decay timescale ⌧ = 1/( DOP + sDOP ).
Note that the analysis depicted in Fig. 8 differs from that of Kwon and Primeau (2006) in several aspects: Firstly, their global biogeochemical model was fully equilibrated (due to their direct evaluation of steady state), whereas simulations of RetroMOPS 15 may still exhibit some drift in nitrogen inventory (see subsection 3.2). Second, Kwon and Primeau (2006) evaluated model sensitivity at b = 1, while Fig. 8 displays a region ±5% around optimal b. Thirdly, Fig. 8 maps only the misfit of solutions realised by the optimisation routine, while Kwon and Primeau (2006) analysed the entire parameter space at b = 1. Finally, the misfit function applied here is based on three components, with very different properties and associated time scales (see above), which can be of advantage for parameter estimation. 20 The misfit to phosphate (Fig. 8, lower left panel) shows an elongated valley in the two-dimensional projection on DOP decay timescale ⌧ (years) and DOP production fraction DOP , and resembles Fig. 4 of Kwon and Primeau (2006). Indeed, one of the lowest misfits to phosphate is achieved with about the same set of parameters as in Kwon and Primeau (2006), namely ⌧ ⇡ 1, DOP ⇡ 0.73. However, nitrate and oxygen show a different, and, partly, antagonistic, pattern: the best fit to observed nitrate is achieved with rather high values of DOP ⇡ 0.8 and ⌧ between about 1-2 years, while the best fit to oxygen is obtained with 25 DOP ⇡ 0.7 and ⌧ ⇡ 1.5 years. The superposition of the different components of the misfit function leads to a unique optimum of ⌧ = 2 ( DOP = 0.47 and DOP = 0.02) and DOP = 0.73 (Table 3). Thus, oxygen and nitrate can provide some useful, independent information on these parameters.
As noted above, the advantage of including nitrate and oxygen in the misfit function is that, in contrast to phosphate, the inventory of these tracers may change freely according to model parameterisation. The resulting bias to observations thus adds 30 two important components to the misfit function, both of which are independent: while high DOP turnover (as simulated by low ⌧ ) biases nitrate low (Fig. 8, upper mid panel), the same value leads to an overestimate of oxygen (Fig. 8, upper right panel; see also Fig. 6). This behaviour can be explained with the different processes and boundary conditions for the two tracers already noted above: a high DOP turnover leads to higher fluxes and a tighter coupling of production and denitrification in upwelling waters, causing a nitrate deficit in the model (see above, and Fig S4). On the contrary, it reduces DOP in preformed waters e.g., in the Southern Ocean, thereby decreasing aerobic remineralisation and oxygen consumption in these waters on their passage towards, e.g., the northern North Pacific. The latter process increases oxygen particularly in deep waters (Fig S3).
The optimal b = 0.98 of RetroMOPS o is lower than that of MOPS oS and MOPS oD . This may be partially explained with the absence of numerical diffusion of detritus in RetroMOPS. As shown by Kriest and Oschlies (2011), in models that explicitly 5 simulate detritus sinking with an upstream scheme, the assumption of homogenous distribution of detritus in each vertical grid box causes an additional, usually downward transport of detritus. This results in an effective b which is about 10-20% smaller (corresponding to faster sinking) than the nominally prescribed b. Optimisation of MOPS accounts for this additional numerical transport by increasing b (= reducing sinking velocity) by some amount. Therefore, optimal b of MOPS without any influence numerical of diffusion would likely be around 1.1-1.2, i.e. closer to b = 0.98 of RetroMOPS o . Considering this effect, optimal 10 b of MOPS oD and, in particular, RetroMOPS o agree with the optimal value of b = 1 found by Kwon and Primeau (2006).
Despite its generally lower fluxes, fixed nitrogen loss in the eastern equatorial Pacific is higher in RetroMOPS o than in MOPS oD (Fig. 2), resulting in a nitrate deficit in this region. The deficit is comparable to that of MOPS oS , i.e. of a model simulation with default parameters for oxidant dependent processes (Fig. 3). Likely, the instantaneous remineralisation of sinking material inherent in the direct flux parameterisation of RetroMOPS, which causes a tighter spatial coupling between 15 production, sinking, remineralisation and upwelling. It has been suggested earlier that the production of slowly degradable organic matter above upwelling regions and/or oxygen minimum zones may help to decouple these processes, and avoid a runaway effect of nitrate loss (Landolfi et al., 2013;Dietze and Loeptien, 2013). The very low optimal value for surface DOP turnover sDOP found in this study supports this finding.
The total misfit to observed dissolved tracer concentrations of RetroMOPS o is only about 4% higher than that of MOPS oD , 20 i.e. RetroMOPS can perform almost as well as MOPS, with respect to annual mean phosphate, nitrate, and oxygen. Simulated biogeochemical fluxes of RetroMOPS o are generally lower than those of MOPS oD , and their horizontal pattern is less pronounced (Fig. 2). This likely arises from the prescribed, constant phytoplankton concentration of RetroMOPS o , which mutes biogeochemical dynamics in productive regions of the high latitudes and upwelling areas. Because RetroMOPS o applies the same parameters as MOPS oD for oxidant-dependent processes, its global fixed nitrogen loss and gain is comparable to that 25 of the more complex model. As for MOPS, optimisation of RetroMOPS against dissolved tracer concentrations improves the fit to global estimates of biogeochemical fluxes (Table 4), and indicates, that these tracers can provide means to calibrate biogeochemical model fluxes on a global scale, even -or especially -for a model as simple as RetroMOPS.

Conclusions
Optimisation of parameters for oxidant-dependent processes results in a slightly better fit to observed tracers, and in a much (e.g., using TMs extracted from the UVic model; Kvale et al., 2017) will help to examine, if a different parameterisation alters the current requirement for very high nitrate threshold of denitrification, that currently helps to prevent nitrate from depletion.
Oxygen and nitrate add important additional constraints on the estimation of biogeochemical parameters. Of particular importance is that, in addition to the spatial information they provide, their flexible inventory introduces the bias as additional information for model calibration. The different time and space scales of processes relavant for their inventory may help to 5 constrain parameters that govern dissolved organic matter production and decay. The effect of these tracers on parameter estimates is of particular importance for models such as RetroMOPS and MOPS, that aim at conserving all oxidants. It may be weaker for models that continue remineralisation even under suboxic and/or low nitrate conditions, thereby implicitly assuming some "hidden" oxidants. In these models it could be useful to track and examine potential oxidant deficits for model evaluation.
The DOP recycling rate affects surface DOP and phosphate concentrations conversely: either the model performs relatively      Table 3. Optimisation results: minimum misfit J ⇤ , optimum parameters and their uncertainties. To determine parameter uncertainty, we selected a group ⌦ of the 1‰ best individuals, i.e. individuals defined by a misfit Ji : Ji/J ⇤ 1  J , with J = 0.001. The number of these individuals N (⌦) is also denoted as fraction n(⌦) of all individuals of the optimisation ⇥ N , where N is the number of generations, and = 10 the population size. For each parameter ⇥ the first column gives the optimal parameter ⇥ ⇤ (i.e., the average parameter of the last generation). The second and third column present the parameter range of all individuals of ⌦, expressed as absolute value (R⇥(⌦)), and normalised by the a priori range of parameters (R A ⇥ ; see Table 1): r⇥(⌦) = R⇥(⌦)/R A ⇥ value.
Experiment:  Table 4. Global annual fluxes of primary production (P), grazing (GRAZ), fixed nitrogen loss through pelagic denitrification (NLOSS), export production (F120, flux through 120 m), flux through 2250 m (F2250), and benthic burial (BUR), in Pg N y 1 , for the reference experiment of MOPS r , MOPS oS , MOPS oD , MOPS oD ⇤ and RetroMOPS, for which we show the fluxes of the (best) reference experiment, RetroMOPS r , the range of all sensitivity experiments, and the optimised run, RetroMOPS o . Also shown are some globally derived, observed estimates. Conversion between different elements was carried out via N:P=16, and C:P=122.