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 re-exposure 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 ACCESS-M (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 re-exposure times) as compared to our best estimate from the data-assimilated model. This is a general feature of optimizing biological-pump 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 spin-ups for each set of parameters to be explored. By contrast, the transport-matrix approach allows us to solve directly for the steady-state, which is computationally orders or magnitude more efficient than running a forward model to equilibrium (steady-state).

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 (1993-present) ?

It is certainly possible to embed the PC or PCO2 biogeochemistry model into any other circulation, including CMEMS-based 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 matrix-free Newton-Krylov), (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.

We thank the referee for the positive feedback. Our point-by-point response follows below:

• (1) Eqn. 3 and lines 84-85: 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 Green-function 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 advection-diffusion, and that this does not impact the solutions significantly while benefitting computational efficiency.

• (2) Eqn. 3 and lines 88-89: "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 steady-state phosphate equation would be singular (rank-deficient). 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 98-100: "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₂ falls below O₂^lim at the same rate as would occur if O₂ was equal to O₂^lim, but without utilizing O₂. 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₂ 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 organic-matter respiration continues when O₂ falls below O₂^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/bg-14-4965-2017)?

For simplicity, DOP is remineralized with a fixed 2-yr timescale at rateD_DOP = [POP]/τ_POM that is independent of [O₂] (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 semi-labile 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 646-647)? Is this a gradient-based 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 gradient-based.  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 zonal-mean preformed and regenerated DIP and preformed O₂ 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: Cross-validation experiments are mentioned (simulating ACCESS with OCIM-optimal 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 cross-validation experiments. 

In a supplementary file, we will also provide basin zonal means of the ACCESS-M tracers (DIP, DIC, O₂) obtained with OCIM2 parameters, and their difference with the optimized ACCESS-M tracers.

• Lines 221ff: "Instead, the Southern Ocean POC respiration rate is weaker for ACCESS-M, 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 ACCESS-M as further discussed in subsection 3.3.2.)" - Why is the respiration rate weaker in ACCESS-M? 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 ACCESS-M are larger. Is this because of lower temperatures in the Southern Ocean in ACCESS-M?

We thank the referee for catching this mistake. We will revise the paragraph to state:

Instead, enhanced vertical mixing dramatically reduces O₂ residence times such that total oxygen is closer to preformed oxygen in the Southern Ocean for ACCESS-M compared to OCIM. This occurs despite the larger ACCESS-M respiration rate coefficients (γ_s and γ_f, Table 1) and lower q₁₀ presumably because of the relatively low organic-matter 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 self-shading 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 ACCESS-M matrix. Hence, the effects of self-shading 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 self-shading.

• Eqn D1: Is there a factor (1-sigma) 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.