the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Optimal parameters for the ocean's nutrient, carbon, and oxygen cycles compensate for circulation biases but replumb the biological pump
Mark Holzer
Matthew A. Chamberlain
Richard J. Matear
Nathaniel L. Bindoff
François W. Primeau
Download
 Final revised paper (published on 26 Jul 2023)
 Supplement to the final revised paper
 Preprint (discussion started on 06 Mar 2023)
Interactive discussion
Status: closed

RC1: 'Comment on egusphere2023363', Anonymous Referee #1, 03 Apr 2023
The paper titled "Optimal parameters for the ocean's nutrient, carbon, and oxygen cycles compensate for circulation biases but replumb the biological pump" presents the tuning of parameters of the PCO2 (or PC) models in 2 cases where it is constrained by 2 different circulations, OCIM2 and ACCESSM.
The paper is very well written and clear The figures are also clear and relevant. Generally, the method is conviincing. The result correspond to what one may expect: different circulations lead to different "optimal" parameters. But the authors go further than that, and show that different pathways are actually relevant in the 2 simuations. Furthermore, the idea of running the model with parameters optimized for 1 forcing, with the other forcing, is relevant. Then the study of potential generalization of the optimized parameters in a different situation (in this case, nutrient depletion) is very interesting.
Thus I have no major comments or criticism, only a few clarifications would be appreciated.
* would it make sense to compare the circulations using simpler tracer models, in order to separate the hydrodynamic and biogeochemical problems ? For example, would it be relevant to circulate an age tracer ? The authors seem to have done something in this direction (around lines 165) ? Another comparison could simply concern the vertical velocities or transports accross particular depths between the 2 circulations ? Also, please specify the vertical discretisation of the 2 circulations (z or sigma etc) to ease the reader's experience.
* in the authors' opinion, would the situation be different be similar when using more complex biogeochemical models ? Could the authors explain how their approach would be different using forward models ?
* the authors use at different places the word "bias" (e.g. line 280) as opposed to the word "errors". I suppose the circulation is indeed biased; can the authors explain whether this is systematic ? Can the present study be used as valuable information for future version of the circulations ?
* As the authors focus on the 1990's, would it be possible to couple PCO2 (or PC) to the Copernicus Marine (CMEMS) GLObal model product (1993present) ?
Citation: https://doi.org/10.5194/egusphere2023363RC1 
AC1: 'Reply on RC1', Benoit Pasquier, 06 Apr 2023
We thank Referee #1 for the positive feedback. Our clarifications to the Referee’s points are as follows:
 would it make sense to compare the circulations using simpler tracer models, in order to separate the hydrodynamic and biogeochemical problems ? For example, would it be relevant to circulate an age tracer ? The authors seem to have done something in this direction (around lines 165) ? Another comparison could simply concern the vertical velocities or transports accross particular depths between the 2 circulations ? Also, please specify the vertical discretisation of the 2 circulations (z or sigma etc) to ease the reader's experience.
Yes, an age tracer can be used to disentangle the effects of circulation and biogeochemistry. Detailed analyses of the model circulations, which included ideal mean age (mean time since surface contact) and mean reexposure time (mean time until next surface contact), have been carried out for ACCESS by Chamberlain et al. (2019) and Holzer et al. (2020), and for OCIM2 by DeVries and Holzer (2019). These references were already cited at various points in the manuscript, but in the revisions, we will also cite them where the circulations are first introduced in the Methods section.
(We mention the ideal mean age in the manuscript around line 165 because we used it to quantify the effect of the horizontal coarse graining (1°×1° to 2°×2°), which guided our reduction of the vertical mixing to compensate for the increased numerical diffusion.)
Following your suggestion, we will add to the revised manuscript that both matrix models use a vertical depth (as opposed to density) coordinate with layer thicknesses that increase with depth and span the range 10–335 m for ACCESSM (50 levels) and 36–634 m for OCIM2 (24 levels).
 in the authors' opinion, would the situation be different be similar when using more complex biogeochemical models? Could the authors explain how their approach would be different using forward models ?
Our main results should apply to more complex biogeochemical models, although our detailed quantitative findings are of course model specific. Optimizing biogeochemical parameters to match observations will compensate for the effects of circulation biases on the biogeochemistry, but the optimal parameters are almost certain to alter the inner workings of the biological pump (with production and export efficiency adjusting to the model’s reexposure times) as compared to our best estimate from the dataassimilated model. This is a general feature of optimizing biologicalpump parameters in the presence of circulation biases independent of the specifics of the model and its complexity.
Using a forward model would lead to qualitatively similar results, although a forward model would have the benefit of capturing seasonal and interannual variability. However, optimizing biogeochemical parameters with a forward model is generally prohibitively expensive requiring lengthy model spinups for each set of parameters to be explored. By contrast, the transportmatrix approach allows us to solve directly for the steadystate, which is computationally orders or magnitude more efficient than running a forward model to equilibrium (steadystate).
In response, we will add some discussion along these lines to the revised manuscript.
 the authors use at different places the word "bias" (e.g. line 280) as opposed to the word "errors". I suppose the circulation is indeed biased; can the authors explain whether this is systematic ? Can the present study be used as valuable information for future version of the circulations ?
We use the term “bias” deliberately because the ACCESS ocean circulation errors are systematic and characterized by coherent patterns as has been reported by Bi et al. (2013) (who also use the language of “biases”). We do hope that our work will be useful for future versions of the ACCESS ocean model, which can then be assessed not just in terms of the fidelity of the physical variables but also in terms of key biogeochemical metrics. In response, we will add a sentence to the revised manuscript inviting other models to be similarly benchmarked using diagnostics of the biological pump.
 As the authors focus on the 1990's, would it be possible to couple PCO2 (or PC) to the Copernicus Marine (CMEMS) GLObal model product (1993present) ?
It is certainly possible to embed the PC or PCO2 biogeochemistry model into any other circulation, including CMEMSbased circulations. However, this is beyond the scope of our study and would require significant work. Specifically, one would at a minimum need to (i) extract a transport matrix as was done for the ACCESS model from a relevant output, (ii) significantly coarsen in the horizontal to make the model numerically affordable or, alternatively, use different solvers suitable for much larger nonlinear systems (e.g., using matrixfree NewtonKrylov), (iii) determine appropriate eddy diffusion schemes, and (iv) carefully assess the resulting deep circulation (e.g., as quantified by ventilation tracers) to ensure fidelity to the parent model. We anticipate no change to the manuscript in response to this point will be required.
Citation: https://doi.org/10.5194/egusphere2023363AC1

AC1: 'Reply on RC1', Benoit Pasquier, 06 Apr 2023

RC2: 'Comment on egusphere2023363', Anonymous Referee #2, 11 Apr 2023
The paper presents the outcome of a simple biogeochemical model that was optimised against observations of dissolved inorganic tracers in two different circulations. Optimal model parameters differ considerably, resulting in very different shallow POC fluxes and export production, but quite similar fractions of regenerated nutrients, and hence similar strengths of the biological pump. The authors show that depending on circulation, optimisation of the biogeochemical model parameters compensates for errors in circulation by inducing different particle transfer efficiencies and sequestration times.
The paper in generally wellwritten, and the results and conclusions provide a very interesting and thorough insight into aspects of largescale biogeochemical cycling and the biological carbon pump on decadal to millennial cycles in global models.
However, I have some open questions regarding the methods (see below), that I think could be addressed in more detail. I further think that the paper would benefit from some more details in the results section, either in the main part, or as a supplement  this relates mainly to preformed and/or regenerated tracers, as well as the results of the crossvalidation experiment (see below).(Methods:)
(1) Eqn. 3 and lines 8485: left hand term S_f[POP_f] and S_s[POP_s] and "is the divergence of the flux of the POP_k tracer with sinking speed w_k"  Does this mean that the POP tracers are not subject to advection and diffusion (as for DOP and DIP)? If so (if particulate and dissolved tracers are treated differently) would this in any way affect the solution/evaluation of the Green function (appendix D1 and D2)?
(2) Eqn. 3 and lines 8889: "The global phosphate inventory is prescribed by weakly restoring DIP to the global observed mean [DIP]"  Why is restoring necessary? Is there any loss or gain of phosphate in the model? In particular: what happens to sinking POP at the lower model boundary (given the " unrealistic POM accumulation in the bottom grid boxes under anoxic conditions" mentioned in line 100)? Is there any loss to the sediment, that necessitates the restoring? How big are the restoring terms for DIP (and TA)?
(3) Lines 98100: "Note, however, that Eq. (2) differs from previous parameterizations in that we implicitly include the effect of microbes switching to nitrate for organic matter oxidization (denitrification) by disallowing respiration rates to decline in anoxic waters below [Olim ] = 5 μM."  Does this mean respiration continues in the absence of oxygen? If so, is the oxygen dept accounted for? Would it play a large role, if it were accounted for in the cost function?
(4) In the model POP may sink very deeply, and dissolves to DOP. Is DOP remineralisation (to DIP) also oxygen dependent as for POP? Does DOP also accumulate at the bottom or elsewhere? Could it be that any improvement in PO4 (through optimisation of sigma) comes at the cost of DOP mismatches (e.g., Kriest, 2017, https://doi.org/10.5194/bg1449652017)?
(5) Optimisation: Could you briefly explain what is behind "MATLAB's unconstrained minimizer" (appendix B4, lines 646647)? Is this a gradientbased method? If so, is there the possibility that the optimisation especially of ACCESS became trapped in a local minimum?
(Results:)
(6) Preformed oxygen simulated by the models is mentioned several times (lines 220, 228, 307)  For me it would be very interesting and helpful to see plots (e.g., zonal means) of this and/or regenerated nutrients, especially to understand the following results on different pathways in the model better.
(7) Lines 239ff "When O2 is prescribed from observations (PC model) ..."  This is one of few references to the PC model, which later serves as an alternative model. Are there any differences in the spatial distribution (zonal means) of DIP? A few words on this (or a plot of zonal means in the supplement) could be helpful, as this model later (Table 1) serves as an alternative biogeochemistry.
(8) Lines 269ff: Crossvalidation experiments are mentioned (simulating ACCESS with OCIMoptimal parameters). As the two different circulations seem to require very different parameters (namely those that affect remineralisation), I think it would be very interesting to see zonal plots of DIP and oxygen for the crossvalidation experiments.
Minor comments:
 Lines 221ff: "Instead, the Southern Ocean POC respiration rate is weaker for ACCESSM, allowing more oxygen to be mixed throughout the water column. (The reduced respiration rate is largely driven by a lower optimal value of γs (Table 1), with POCs dominating respiration for ACCESSM as further discussed in subsection 3.3.2.)"  Why is the respiration rate weaker in ACCESSM? From what I see in Table 1, (basic) optimal respiration rates are 1.05 and 0.535 for large and small POP, and thus much larger than those of OCIM? Likewise, the Q10 values of ACCESSM are larger. Is this because of lower temperatures in the Southern Ocean in ACCESSM?
 Eqn 2: Is temperature given in K or °C? (because of "K" in the denominator of the exponent)
 Eqn A1 and line 578: The same PAR (light) was used for both models  does this relate only to the surface light, or also the light at depth?? If the latter, then the approach neglects selfshading by phytoplankton, and hence any feedback effects of phytoplankton on production, correct? This might be an important difference to more complex BGC models, together with the assumed constant phytoplankton concentration.
 Eqn D1: Is there a factor (1sigma) missing in this equation (as in Eqn 3)? Would this also affect the following equations (D4, D6, D7)?
 Table B1 (and Table 1): It would be very helpful to have the meaning of the parameters (e.g., "fraction of production routed to DOP" for sigma) in at least one of the tables.
Citation: https://doi.org/10.5194/egusphere2023363RC2 
AC2: 'Reply on RC2', Benoit Pasquier, 08 May 2023
We thank the referee for the positive feedback. Our pointbypoint response follows below:
 (1) Eqn. 3 and lines 8485: left hand term S_f[POP_f] and S_s[POP_s] and "is the divergence of the flux of the POP_k tracer with sinking speed w_k"  Does this mean that the POP tracers are not subject to advection and diffusion (as for DOP and DIP)? If so (if particulate and dissolved tracers are treated differently) would this in any way affect the solution/evaluation of the Green function (appendix D1 and D2)?
Yes, the particles are not advected and diffused along with the water. This approximation is justified because the transport of particles is generally dominated by their sinking. The approximation is desirable because adding the full advection–diffusion operator for the particles would reduce the sparsity of the system of equations, making it computationally more expensive to solve. This does not affect the Greenfunction approach of our diagnostics in any way. In the revisions, we will explicitly state that the particles are subjected only to gravitational sinking and not to water advectiondiffusion, and that this does not impact the solutions significantly while benefitting computational efficiency.
 (2) Eqn. 3 and lines 8889: "The global phosphate inventory is prescribed by weakly restoring DIP to the global observed mean [DIP]"  Why is restoring necessary? Is there any loss or gain of phosphate in the model? In particular: what happens to sinking POP at the lower model boundary (given the " unrealistic POM accumulation in the bottom grid boxes under anoxic conditions" mentioned in line 100)? Is there any loss to the sediment, that necessitates the restoring? How big are the restoring terms for DIP (and TA)?
The total amount of phosphate is conserved in our model. Because we are directly solving for the steady state, we cannot specify an initial phosphate inventory. Without the restoring term, the steadystate phosphate equation would be singular (rankdeficient). The weak (“geological”) restoring removes this deficiency (null space) from the system. (This has been discussed in several previous publications by the authors, e.g., Kwon and Primeau, 2005; Primeau, Holzer, and DeVries, 2013; Holzer, Kwon, and Pasquier, 2021.) POP that reaches the bottom grid box superjacent to the seafloor accumulates there (no loss to the sediments) until it is eventually remineralized to DIP or solubilized to DOP, both of which are then transported with water and mixed back into surrounding grid boxes. The DIP restoring term is multiple orders of magnitude smaller than all other DIP tendencies. In response, we will add a sentence stating that the total amount of phosphate is conserved in our model and explaining what is accomplished by the “geological” DIP restoring.
 (3) Lines 98100: "Note, however, that Eq. (2) differs from previous parameterizations in that we implicitly include the effect of microbes switching to nitrate for organic matter oxidization (denitrification) by disallowing respiration rates to decline in anoxic waters below [Olim ] = 5 μM."  Does this mean respiration continues in the absence of oxygen? If so, is the oxygen dept accounted for? Would it play a large role, if it were accounted for in the cost function?
Yes, this means that respiration continues when O_{2} falls below O_{2}^{lim} at the same rate as would occur if O_{2} was equal to O_{2}^{lim}, but without utilizing O_{2}. The effect of this on oxygen deficit (as measured by TOU for example) is consistently accounted for by the oxygen equation without the need for extra modelling because we optimize the mismatch with O_{2} not oxygen deficit. We experimented with various other parameterizations of denitrification and found the solutions to be very similar. In response, we will add a sentence making it explicit that organicmatter respiration continues when O_{2 }falls below O_{2}^{lim}.
 (4) In the model POP may sink very deeply, and dissolves to DOP. Is DOP remineralisation (to DIP) also oxygen dependent as for POP? Does DOP also accumulate at the bottom or elsewhere? Could it be that any improvement in PO4 (through optimisation of sigma) comes at the cost of DOP mismatches (e.g., Kriest, 2017, https://doi.org/10.5194/bg1449652017)?
For simplicity, DOP is remineralized with a fixed 2yr timescale at rate_{}D_{DOP} = [POP]/τ_{POM} that is independent of [O_{2}] (see also lines 101–103). DOP does not accumulate on the seafloor, but is transported with seawater. Only POP is allowed to accumulate at the seafloor, in which case the accumulated POP is eventually remineralized to DIP or dissolved into DOP, both of which are transported with water into the neighbouring boxes. It is entirely that improvements in the DIP mismatch come at the cost of DOP realism. However, we deliberately do not use DOP observations as a constraint on our model because our representation of DOP with a single semilabile pool is too simplistic to be properly constrained by observations. No changes to the manuscript.
 (5) Optimisation: Could you briefly explain what is behind "MATLAB's unconstrained minimizer" (appendix B4, lines 646647)? Is this a gradientbased method? If so, is there the possibility that the optimisation especially of ACCESS became trapped in a local minimum?
We used the MATLAB routine “fminunc”, which is indeed gradientbased. While it is theoretically possible that the optimization is trapped in a local minimum, during this research we performed many optimizations for slightly different versions of the model and they generally approached a similar state. In response, we will provide the name of the MATLAB routine in the revised manuscript and mention that this optimization method does not guarantee convergence to the global minimum.
 (6) Preformed oxygen simulated by the models is mentioned several times (lines 220, 228, 307)  For me it would be very interesting and helpful to see plots (e.g., zonal means) of this and/or regenerated nutrients, especially to understand the following results on different pathways in the model better.
We will provide plots of basin zonalmean preformed and regenerated DIP and preformed O_{2} and TOU in a supplementary file (not part of the Appendices).
 (7) Lines 239ff "When O2 is prescribed from observations (PC model) ..."  This is one of few references to the PC model, which later serves as an alternative model. Are there any differences in the spatial distribution (zonal means) of DIP? A few words on this (or a plot of zonal means in the supplement) could be helpful, as this model later (Table 1) serves as an alternative biogeochemistry.
There are minor differences in the zonal means of DIP (±0.1μM) between the PCO2 and PC models. In response, we will provide plots of the basin zonal means of DIP in a supplementary file.
 (8) Lines 269ff: Crossvalidation experiments are mentioned (simulating ACCESS with OCIMoptimal parameters). As the two different circulations seem to require very different parameters (namely those that affect remineralisation), I think it would be very interesting to see zonal plots of DIP and oxygen for the crossvalidation experiments.
In a supplementary file, we will also provide basin zonal means of the ACCESSM tracers (DIP, DIC, O_{2}) obtained with OCIM2 parameters, and their difference with the optimized ACCESSM tracers.
 Lines 221ff: "Instead, the Southern Ocean POC respiration rate is weaker for ACCESSM, allowing more oxygen to be mixed throughout the water column. (The reduced respiration rate is largely driven by a lower optimal value of γs (Table 1), with POCs dominating respiration for ACCESSM as further discussed in subsection 3.3.2.)"  Why is the respiration rate weaker in ACCESSM? From what I see in Table 1, (basic) optimal respiration rates are 1.05 and 0.535 for large and small POP, and thus much larger than those of OCIM? Likewise, the Q10 values of ACCESSM are larger. Is this because of lower temperatures in the Southern Ocean in ACCESSM?
We thank the referee for catching this mistake. We will revise the paragraph to state:
Instead, enhanced vertical mixing dramatically reduces O_{2} residence times such that total oxygen is closer to preformed oxygen in the Southern Ocean for ACCESSM compared to OCIM. This occurs despite the larger ACCESSM respiration rate coefficients (γ_{s} and γ_{f}, Table 1) and lower q_{10} presumably because of the relatively low organicmatter production in the polar Southern Ocean (Figure 4 below).
 Eqn 2: Is temperature given in K or °C? (because of "K" in the denominator of the exponent)
This equation involves the temperature difference T − T_{ref}, for which it does not matter if the unit is °C or K. Having the denominator as 10K has the advantage of being unambiguous — if 10°C was used, one might misinterpret the denominator as (10+273.15)K, which would be incorrect. We think the equation is clear as written.
 Eqn A1 and line 578: The same PAR (light) was used for both models  does this relate only to the surface light, or also the light at depth?? If the latter, then the approach neglects selfshading by phytoplankton, and hence any feedback effects of phytoplankton on production, correct? This might be an important difference to more complex BGC models, together with the assumed constant phytoplankton concentration.
Yes, we prescribe the same approximately exponentially attenuating PAR throughout for both models, taken from the same ACCESS 1.3 run used to build the ACCESSM matrix. Hence, the effects of selfshading are not modelled for simplicity. We agree that a more complex forward model could account for this. In the revisions, we will explicitly mention that the same PAR is used throughout the water column for both models and that PAR is not coupled to the plankton concentration, precluding effects from selfshading.
 Eqn D1: Is there a factor (1sigma) missing in this equation (as in Eqn 3)? Would this also affect the following equations (D4, D6, D7)?
We thank the Referee for catching this typo – a factor (1σ) was indeed missing, which will be corrected in the revised version.
 Table B1 (and Table 1): It would be very helpful to have the meaning of the parameters (e.g., "fraction of production routed to DOP" for sigma) in at least one of the tables.
We agree and will add a descriptive column to Table 1 in the revised manuscript.
Citation: https://doi.org/10.5194/egusphere2023363AC2

AC2: 'Reply on RC2', Benoit Pasquier, 08 May 2023