Review of Dietze and Löptien’s “Argon Saturation in a Suite of Coupled General Ocean Circulation Biogeochemical Models off Mauretania” (manuscript #: 2024-918)

General comments:

The authors of this manuscript examine the model-simulated inaccuracies in diapycnal mixing that cannot be explained by the horizontal resolution near Mauretania and find that it’s comparable to advection numerics and choices of background diffusivity by making use of argon saturation as a proxy for effective mixing. I can appreciate how much time and effort that went into this work but I have some reservations about publishing this manuscript as is. One concern that I have is that there is no comparison of the effective diapycnal mixing with other methods, like that of Holmes et al. (2021) or potentially separating out the spurious contribution to compare with a method like that developed in Ilıcak (2016)'s "Quantifying spatial distribution of spurious mixing in ocean models”.  I understand the authors just want to have a relative measure to rank their model configurations with their effective mixing against one another but without knowing whether their method is properly quantifying the effective mixing, the ranking could be misleading. For example, the authors include the mixed layer in their calculations of the effective mixing, despite how there are known confounding effects on the effective mixing. The second issue is that the authors try to connect Ekman pumping/suction to effective diapycnal mixing without explaining the relationship between the diapycnal advection and diapycnal mixing in the Ekman layer. The third issue is that the model spin-up time may suffice for many applications of model experiments, but they are investigating a signal that is so small that the rate at which their model diagnostics are changing at the time they are evaluating them may be relevant to the significance of the small differences they’re finding. I suggest major revisions. Specific comments are listed below:

Specific comments:

Line 3: The authors say “unacquainted vertical” by which they are referring to “diapycnal” but that way of describing diapycnal sounds more literary than it sounds like it has specificity; did the authors use AI to help them write this abstract?

Lines 38-39: The canonical reference for uncertainties in representation of diapycnal mixing (not just vertical) is MacKinnon et al. (2017)’s “Climate Process Team on Internal-Wave Driven Ocean Mixing” even though there’s also shear-driven and other mixing

Line 45: This seems like more “coupled” reasoning than circular because the biogeochemical modules can potentially be developed to be more realistic without developing the diapycnal diffusivity modules, but of course the biogeochemical variables depend upon the diapycnal diffusivities. It’s difficult to point to new model innovations in ESMs that don’t have this coupling problem. Your argument seems analogous to saying that assimilating sea-ice concentrations with the goal of improving the sea surface properties is circular reasoning because our sea-ice modules are imperfect. In any case, without improvements in the biogeochemical modules, the information propagated by the assimilation of biogeochemical variables may not inform the diapycnal diffusivities as much as the biogeochemical parameters themselves, which is the problem you’re pointing out. This would certainly be the case with the Green’s function-based approach of ECCO-Darwin (see Carroll et al. (2020)’s “The ECCO-Darwin Data-assimilative Global Ocean Biogeochemistry Model: Estimates of Seasonal to Multi-decadal Surface Ocean pCO2 and Air-sea CO2 Flux”) because the physical variables aren’t changed in their formulation; only the biogeochemical parameters are altered. I believe the proposal by Trossman et al. (2022) is that the ECCO-Darwin and ECCO setups could be used in an iterative process to incrementally improve the biogeochemical parameters and the physical parameters, respectively, but this is even more computationally expensive than the current process and may require a different process for sequential data assimilation systems, which they present unique issues with in their appendix. Additional evidence that this could work comes from Skakala et al. (2022)’s “The impact of ocean biogeochemistry on physics and its consequences for modelling shelf seas,” where it was shown that assimilation of a suite of biogeochemical variables had an impact on the diapycnal diffusivities beyond that of physical variables.

Lines 63-64: Will this always be the case? The atmospheric community seems to have come up with a solution to advection numerics in McGraw et al. (2024)’s “Preserving Tracer Correlations in Moment‐Based Atmospheric Transport Models”. Also, is it possible that spurious diapycnal mixing is the result of spurious vertical advection sometimes? This certainly happens when you perform sequential assimilation (see: Pilo et al. (2018)’s “Impact of data assimilation on vertical velocities in an eddy resolving ocean model”).

Lines 81-82: Isn’t there also the effect of surface disequilibrium at the time of water mass formation (e.g., when subduction occurs) that needs to be accounted for? Both the nonlinearity of solubility in the presence of mixing and surface disequilibrium effect cause apparent oxygen utilization to be different from true oxygen utilization, for example. So the argon tracer needs to be injected below the mixed layer such that it doesn’t get obducted throughout the experiment. Also, argon can’t be injected into the ocean to derive diapycnal diffusivities from observations (unless the background levels are well-known) because there are significant background levels of argon in the ocean, which will eventually render it another exotic tracer on top of the ones Ledwell used. I mention this because it was noted that argon could be used as a proxy for effective diapycnal mixing in both models and observations in the text.

Line 150: I understand not wanting to compete with Jason Momoa in Google searches by dubbing your model’s acronym MOMOA, but you’re going to compete with the Museum of Modern Art in New York City instead. In any case, the model and its tracer conservation properties won’t be as much of a problem as they would be if the authors chose another model like HYCOM. And the model domain is reasonable, considering problems with boundary conditions that can occur.

Line 169: You say that there are 55 geopotential levels in all model configurations here but in Table 1, it says there are 72 grid points in an undefined dimension. The caption of Table 1 only mentions zonal and meridional directions but there is a third dimension listed for the number of grid points. Is that dimension time or is one of your instances of vertical dimensions inaccurate (i.e., in Table 1 or on Line 169)?

Line 186: For future equilibrations of biogeochemical variables, you could consider using anther method like the one by Khatiwala (2024)’s “Efficient spin-up of Earth System Models using sequence acceleration”

Line 189: I would argue that use of COREv2 forcing is no longer considered standard, as JRA55-do has supplanted COREv2 in OMIP2 and ERA5 is also commonly used

Lines 201-203: If you’re interested in the evolution of each of the model configurations, then you should be assessing metrics that are related to the time rate of change, not averages

Lines 205-206: One year of spin-up is common amongst published high-resolution model simulations (based on the global kinetic energy metric you seem to be using) because the initial shock is essentially adjusted away but the model’s climatological state is still being approached asymptotically. That would need at least several more years to achieve. It may take much longer to equilibrate potential energy, which includes available potential energy relevant to the eddying circulation, but to first order using kinetic energy may be good enough.

Table 2: Continuing on from my previous comments, please include the time rate of change (linear slope of a regression against time, for example) of the Delta Ar[%] metric you include the average of in this table. That would help the reader see whether these values are close enough to what they would be in a more well-spun-up state. Also, why are you including the mixed layer in your calculations when Argon is not a conservative tracer near the sea surface and diapycnal mixing is not very meaningful within a bulk mixed layer? Are you not using a bulk mixed layer?

Lines 234-236: Again, because penetrating solar radiation, bubble entrainment, and possibly other processes affect Delta Ar[%], why are you including the mixed layer in your calculations?

Lines 266-268: This is qualitatively accurate. However, it can be noted that there are in-situ estimates of eddy kinetic energy compared with high-resolution model simulations in the literature, like those presented in Luecke et al. (2020)’s “Statistical Comparisons of Temperature Variance and Kinetic Energy in Global Ocean Models and Observations: Results From Mesoscale to Internal Wave Frequencies”. Also, along-track SWOT data resolves the higher wavenumbers and there are preliminary results that much more of the spectrum is resolved with SWOT. I’m not suggesting that you perform a comparison with in-situ (e.g., moored) or SWOT data. This is just something to be aware of.

Figure 3: It’s curious that the relatively warm SST feature near 18W, 21N in the High configuration is completely absent in Coarse, except for what seems like a separate merging feature near the upper part of the domain plotted here. The rest of the figure comparison looks like the High configuration’s SST was coarse-grained.

Table 3 and Figure 6: Okay, the higher the resolution, the greater the Ekman pumping/suction there is in its mean and variability, but how come you’re spending time on this when you’re not relating it to the effective diapycnal mixing? You’re relating how resolving the wind stress curl and upper-ocean features to diapycnal fluxes of tracer properties within the Ekman layer, which has more to do with advection. It’s possible, however, that, by conversation of volume in your model, diapycnal mixing will counter advective diapycnal fluxes but you are not showing this.

Lines 313-319: The differences in domain-averaged Delta Ar[%] amounting to <0.1% with the same back diffusivities (equivalent to increasing the background diffusivity from 20% to 40%) may not be significant/detectable here once you consider the rate at which Delta Ar[%] and/or global kinetic energy is changing as the simulations are spinning-up. It appears that in Figure 7d that the seasonal variability and trend in (Delta?) Ar[%] may be larger than the difference across High and Coarse. So it’s unclear whether the comparisons you make with the comparable effect of varying the background diffusivity and/or numerical advection scheme will hold up if your model was spun up more. A decrease in Delta Ar[%] from using a background diffusivity of 0.14 cm^2 s^-1 to 0.16 cm^2 s^-1 suggests varying rates of spin-up across different background diffusivities, for instance.

Minor corrections:

Line 4: “effected” should be “affected”

Lines 43 and 668: Trossmann should be Trossman

Line 209: “… effects associated to our…” should be “… effects associated with our…”

Line 233: “… that would allows for…” should be “… that would allow for…”

Line 261: “… starts to contributes to…” should be “… starts to contribute to…”

Line 295: “looses” should be “loses”

Line 419: Steven Griffith should be Stephen Griffies, I believe; Eric Galbreith should be Eric Galbraith

Review of “Argon Saturation in a Suite of Coupled General Ocean Circulation

Biogeochemical Models off Mauretania” by Dietze and Löptien

Summary:

The main finding of this paper is that diapycnal mixing, a key source of uncertainty in ocean models, is strongly influenced by factors such as model resolution, advection schemes, and background diffusivity. Using Argon saturation as a proxy, the study shows that the effect of differences in resolution is comparable to that of other factors. This study highlights the ongoing challenge of accurately representing ocean mixing in models used for climate projections. I found this manuscript concise and well-written. However, some clarifications in methods and a few minor corrections are needed before publication. I suggest minor revisions. Please see the general comments and specific comments below.

General comments:

1. Consider changing the title to highlight the key findings or purpose of this study. The current title is concise, but it could be refined for better focus and impact.
2. While the rationale for using a 3-year simulation period is clear in terms of model spin-up and avoiding equilibrium, I have concerns about the potential influence of spurious effects from the closed boundaries on the region of interest. Could the authors provide diagnostics or an analysis that evaluates how far spurious effects, such as reflected waves or boundary-induced disturbances, propagate into the domain over the 3-year period? This would help confirm whether the region of interest remains unaffected.

Detailed comments:

1. Line 134: What is the horizontal resolution of MOMSO? Is it 50 km as shown in Fig. 1?
2. Line 165: “resolutions” should be “resolution.”
3. Line 169: I have a question regarding the vertical grid description in line 169, where it is mentioned that the model uses 55 vertical levels but 72 vertical grid points in Table 1. Could you kindly clarify this aspect? Is the difference in grid points due to the use of staggered grids, boundary conditions, or another feature of the model's vertical discretization? If there is a specific reason for this configuration, I would appreciate further details as a reader.
4. Figure 1 caption: Panel (c) and panel (d) don’t show vertical resolution, right? Am I missing anything?
5. Line 188: I believe “affect” would be better than “effect” here.
6. Line 225: It seems to me that equation (1) is incorrect. Since delta_Ar = Ar - Ar_sat, why do we need to subtract 100% in this equation? This might be a typo, but if not, please check the calculations presented in this study. Or am I missing anything?
7. Line 229: Could the authors elaborate on why mixing always causes an increase in saturation? Wouldn't the sign of the change in saturation depend on the temperature before and after mixing, given that Argon saturation is a non-linear function of temperature?
8. Line 229: There is a typo. Change “his” to “is.”
9. Line 233: “Allows” should be “allow.”
10. Line 266: There is an extra parenthesis.
11. Line 295: “loses” should replace “looses.”