The ¹⁷Δ_O2 preserved in ice cores is a key proxy for reconstructing past global biosphere productivity, yet interpreting this signal requires a robust understanding of both the oceanic and terrestrial end-members. Previous studies have incorporated oxygen isotopes into the ocean biogeochemical component of CESM, but a fully coupled Earth system modeling framework has remained unavailable. This study addresses this gap for the oceanic component by implementing the triple oxygen isotope composition of dissolved oxygen (¹⁷Δ_ocean) into the intermediate-complexity Earth System Model iLOVECLIM.

The authors incorporate the three principal processes controlling the isotopic composition of dissolved O₂—photosynthesis, respiration, and air–sea gas exchange—and evaluate model performance against a compilation of 2482 observational measurements under pre-industrial boundary conditions. The manuscript is generally well structured, the methods are described in sufficient detail, and the model–data comparison is comprehensive. I believe this work represents a valuable contribution and is suitable for publication in Geoscientific Model Development after the authors address the comments below.

Major comments

1. Equations (2), (5), and (8): If I understand the derivation correctly, a factor of 2 may be missing before the isotope ratio term (R) in Eqs. (2), (5), and (8). The conversion between isotopologue abundances and bulk isotope ratios for O₂ requires a factor of 2 arising from the symmetry of the O₂ molecule (Yeung et al., 2012). Although this omission may have only a minor impact on the numerical results if implemented consistently throughout the model, the equations should be verified and clarified to ensure that the formulation is physically and chemically correct.
2. Sensitivity tests of key parameters: The model evaluation relies on a number of prescribed parameters, including the coefficient of 2 in Eq. (10), the respiration fractionation factors (¹⁸α_respiration = 0.980 and ^17/18θ_respiration = 0.518), the atmospheric ¹⁷Δ_O2 boundary condition, and the initial dissolved O₂ concentrations and isotope compositions. Several of these parameters remain incompletely constrained and may exert a first-order control on the simulated distributions of δ¹⁸O_ocean and ¹⁷Δ_ocean. Moreover, the manuscript repeatedly suggests that some of the model–data discrepancies, particularly the overestimation of ¹⁷Δ_ocean in deep waters and the Southern Ocean, may result from the choice of respiration fractionation parameters. However, this hypothesis is not evaluated quantitatively. Given that the implementation of oxygen isotopes is the primary development presented in this study, the manuscript would be substantially strengthened by including sensitivity experiments that explore plausible ranges of the key isotope parameters. Such analyses would help determine whether the remaining model–data mismatches primarily reflect parameter uncertainty or structural limitations of the model itself, while also providing a useful assessment of the robustness of the simulated isotope fields. 
3. Air–sea gas-exchange isotope fractionation: The manuscript does not explicitly state the oxygen isotope fractionation factors used for air–sea gas exchange. Because air–sea exchange is one of the three fundamental processes controlling ¹⁷Δ_ocean, these values should be clearly documented. In addition, experimental studies have demonstrated that both kinetic and equilibrium isotope fractionation during O₂ gas transfer can substantially influence dissolved O₂ isotopologue compositions and, consequently, gross primary productivity estimates derived from ¹⁷Δ_ocean (Li et al., 2019). This effect may be particularly important in low-productivity regions where the biological signal is relatively weak. I therefore recommend that the authors explicitly report the fractionation parameters used in the model and discuss their potential influence on the simulated isotope fields.

Minor comments

1, In Eq. (7), should the asterisks on O_surf and O_eq be removed? Please verify the notation.

2, Line 184: “the the model” contains a duplicated word.

3, Table 1 uses inconsistent formats for reporting parameter uncertainty. Some entries are presented as ranges (e.g., [0.982–0.990]), whereas others are reported as mean ± uncertainty (e.g., 0.980 ± 0.003). Consider adopting a consistent format throughout the table.

4, In Table 2, the simulated global GPPO₂ value (9.68 × 10¹⁵ mol O₂ yr⁻¹) is compared with estimates from Huang et al. (2021). However, the Huang et al. estimates represent the modern ocean, whereas the model is run under pre-industrial boundary conditions. Although the difference in global productivity between these two states may be modest, this mismatch should be acknowledged. The same caveat applies to the comparison of isotope tracers with modern observational datasets.

References

1, Yeung, L. Y., Young, E. D., & Schauble, E. A. (2012). Measurements of 18O18O and 17O18O in the atmosphere and the role of isotope-exchange reactions. Journal of Geophysical Research, 117(D18).

2, Li, B., Yeung, L. Y., Hu, H., & Ash, J. L. (2019). Kinetic and equilibrium fractionation of O2 isotopologues during air–water gas transfer and implications for tracing oxygen cycling in the ocean. Marine Chemistry, 210, 61–71.

Overall, I think this is a useful advance to be published, and brings together understanding implemented in models of different parts of the Earth System into a coupled ocean-atmosphere framework.

General points to address:

More justification is needed for the simplification 'GPPO2 = 2 × PP'. You state from Juranek and Quay that the range in the ocean is 18-4.5 and that the ratio is close to 2.7, then chose a value of 2. If going for a single number, why not 2.7? Then whet is the implication of not considering the spatial variability? For example, if you were to consider spatial correlation between the GPPO2/PP ratio and e.g. spatial variability in O2 fluxes, does that lead to a bias in from the mean impact of this assumption that might suggest that you should be using a higher or lower value? As I say, simplification is quite understandable, it just needs more justification and explanation.

I'm not 100% clear as to whether oxygen is prognostic or not in the model. You state that you are inferring O2 fluxes from CO2 fluxes, which in my view may be acceptable in steady state system, but O2 and CO2 have vastly different equilibration timescales due to the slow hydration reaction of aqueous CO2. But you also describe the O2 flux calculation in the paragraph starting on line 176. I would suggest having a very clearly written summary at the start of the implementation section so the reader is clear when they go into this section what assumptions you are making, then you can fill them in on the detail subsequently. Please then explore with clear justification, how if these assumptions were not made the results would change, and in the context of this, explain why the choices made here are (or are not) acceptable.

The observational comparison is hard to take a lot from in this stud, as the model physics and background O2 concertation are not very good (for understandable reasons given the nature of the model). I wonder if it would be more useful to go the observational comparison in water-mass space rather than geographical coordinates, so the reader can assess the d17O more meaningfully?

Specific points to address:



'…provides a global integrator of past biospheric oxygen fluxes, and by extension carbon fluxes'. Not all processes change C and O together, particularly on longer timescale - moderate this statement.

Wording flow: 'However, deconvolving the signal of 17Δ requires to isolate the oceanic biosphere productivity (17Δocean).' -> ‘However, deconvolving the signal of 17Δ requires one to isolate the oceanic biosphere productivity (17Δocean) from …..[state what it is deconvolved from]’

Line 39, "have employ isotope-enabled." -> "have employed isotope-enabled."

Line 49, missing a space in ").T"