The article provides an estimate for one of the components of ocean eddy energy pathway (through relative wind stress), partly for improving our understanding of ocean energetics, and partly for informing the use of eddy energy constrained parameterisations. The topic is relevant to the journal and is of interest, although I am biased because I have publications in the related area... The methodology is reasonable, with an toy model that provides a analytical expression for the dissipation rate that can be applied to ocean datasets, with appropriate caveats stated; although see some technical comments I have below. Comments I have are mostly presentation related, but there are a few technical ones that I think the authors could provide a clarification for in a revised version (and yes those are probably going to be innocuous sounding but potentially annoying questions).

Technical comments

• There is the paper of Rai et al. (2021) [Rai et al., Sci. Adv. 2021; 7 : eabf4920] that looks into eddy killing, providing an estimate of energy loss from the ocean arising from relative wind stress. The paper is close enough in spirit that appropriate citations should be made, probably in multiple places throughout the text, and in particular highlighting how the present methodology is different. Would be good to see a comparison on how the present results are similar/different to the results there.
  
  
  
  (There is also a Rai et al (submitted) but that's not published yet, and may be not as relevant as the 2021 paper.)
• A two-layer shallow water model is employed as a way to have a representation of the first baroclinic mode, however the choice of "mode" (or the basis) is the "standard" choice assuming flat bottom essentially. That choice of basis is not unique and increasingly there are studies using surface modes (e.g. LaCasce 2017; Groeskamp et al. 2020; Ni et al., accepted at GRL or 2023?) I am not suggesting you redo your analysis using surface modes, but given there is a choice of basis, minimally it would be good to speculate how the results are going to be the same and/or different (although of course if you redo the analysis, instead of speculating you could then quantify the similarity and differences, and strengthening the conclusions).
• Regarding Eq. (19), I'll be honest and didn't think about this as much as I probably should have. Two points here:
  
  
  
  1) The notation is a little un-rigourous, because it's not clear which quantities are scalars, vectors and/or tensors, missing some contraction operators and the like. For example, the final term on the right hand is written $(\nabla^2 u)^2$, but you really mean $\nabla^2 u \cdot \nabla^2 u = |\nabla^2 u|^2$. Is the second term a tensor ($|\nabla u|^2$) the hit by the Laplacian operator? This is mostly notation but probably could be cleaned up a bit.
  
  
  
  2) Normally in shallow water it is known that simple choices of the Laplacian as a diffusion leads to sign indefinite energy dissipation (e.g. Peter Gent in 1993, "The Energetically Consistent Shallow-Water Equations"; Gilbert et al., 2014, "On the form of the viscous term for two dimensional Navier–Stokes flows"; in my PhD thesis). One offending reason is that the primitive variables in shallow water are (h, u, v), but the conservative variables are (h, U = hu, V = hv). If for example in the prognostic equation we have $-\nabla^2 u$, then multiplying by $hu$ then integrating by parts gives (abusing notation a bit)
  
  
  
  \int hu \cdot \nabla^2 u =  boundary terms - \int \nabla (hu) \cdot \nabla u,
  
  
  
  which gives the expected sign-definite term $-h|\nabla u|^2$, but there is a cross term involving $\nabla h$ floating around. The problem however doesn't exist in the conservative variables. It seems somehow you don't have this issue here? I am not seeing how your terms could be written as a flux given there is that annoying $h$ term floating around, so a clarification would be useful.
  
  
  
  Ultimately I assume it's quantitatively not going to be important, because your "hyperdiffusion" $\nabla^4 u$ (I put quotation marks because I am arguing $u$ is not the variable you should be hitting $\nabla^4$ with) is presumably going to be a small effect. Clarification and maybe appropriate references here would be helpful (e.g. Gent 1993; note also he makes the point about energetic consistency, which I am not convinced you have here, but it's probably not too important).
• So intuitively the action of the wind is at the surface, but that is being argued to be felt over the structure of the baroclinic mode, so there is implicitly some assumption about the time-scale of communication. Do we actually know what that is, and is that small enough for the conclusions here to internally consistent? I think I wouldn't raise this question so much if surface modes were used for example. Clarifications on this point would be welcome, maybe commenting on the evolution of the MITgcm model results.
• The toy model uses an eddy as a noun (quasi-circular object), but results are then applied generically to a velocity with no specific mention to the eddy, so I assume there is no eddy detection that is being done here, and the application is then considering eddies as the verb (the fluctuation, a special case being the quasi-circular object). Could the authors comment on this distinction, and which one would be the more appropriate thing to do? One could for example argue that we might want to identify the coherent eddies (so the noun, via Eulerian or Lagrangian approaches), then consider the impact of relative wind-stress on those identified eddies for consistency of the theory and application.
  
  
  
  (I don't personally believe you should do what I just suggested. I am just raising the point that the word "eddy" is sometimes used to mean different things by different people, and sometimes the intention is not as precise as it should be.)
• Maybe a hypocrite for asking this (because I didn't do it either in the 2022a paper), but a dissipation rate is estimated here but no estimation of a power? Could you estimate a power, and how might that compare with the results of Rai et al (2021) say? Or if you refrain from doing so, give a reason on why you don't?

Presentation comments

• (line 6): unless you can definitively say and show the "nonlinear baroclinic processes" causality, I would lessen the strength of the wording and say "there is divergence from the analytical model at around day 150, likely due to the presence of nonlinear baroclinic processes" (because it's really more an observation at the moment).
• (line 12 and later): the 10^-7 is mentioned but its significance (or lack of) is never given explicitly. Easiest to say up front why that that reference value is chosen.
• (two sentences spanning line 10-12): reads a bit clunky, could do with a re-write.
• (line 34-35): I would argue that's not a good comparison, because the low explicit eddy energy is to do with the coarse resolution model and much less on the GM parameterisation itself. I would personally just remove that sentence.
• (line 45-47): Jumpy sentence, consider rewrite (eddy saturation <-> GEOMETRIC while "other" <-> turbulent energy cascade, and as written the is ambiguity in how the sub-clauses are related to each other).
• (line 59): formatting of reference, brackets.
• (paragraph of line 60): Rai et al. (2021) should be cited and results compared accordingly in this paragraph.
• (line 95): $g$ should be the gravitational acceleration constant
• (line 95): "$\nabla_h$ IS THE horizontal gradient operator" or similar
• (sentence of line 114-115): Could be read as there is current feedback onto the wind profile, which I assume is not what was intended.
• (equation 9): consider using \begin{align} \end{align} with some & according to break the lines, probably
  
  \begin{align}
  
      W_{rel} &= \tau_rel \cdot u_g \\\
  
                    &= etc. \\\
  
                    &=
  
  \end{align}
• (line 143-145): citation to Rai et al. (2021) here also probably.
• (section 2.2 opening paragraph): probably comment about surface modes here, or say the discussion will be given in the conclusion section
• (line 163): "...where $\cdot_{1,2}$ denotes the upper and lower layer variables,..."
• (line 202, 203): might consider swapping $u$ and $\eta$ ordering so you can have equation reference ordering as (5) and (6)
• (line 223): not sure why you wouldn't just RK4 the whole thing, comment to clarify? (Because piggy-backing on the MITgcm AB3 time-stepper?)
• equation (23): \left( and \right) for brackets
• (liner 262): so you have sponge layers or diffusion to soak it up? clarification would be useful (if discussed in previous work, citation here would also be appropriate).
• (line 273): I might have opted to define $n_0(z)$ as the NEGATIVE vertical gradient to soak up that negative sign in $PE$ partly because it is never used again anyway (I'm not remotely attached to this suggestion).
• (line 318): remove comma maybe, "...but with vertical diffusion did result..."
• (line 343): "...zonal wind velocity, because wind patterns..."
• (line 345): "...is on a 2 degree horizontal grid."
• (line 355): comment on differences / similarities with Rai et al. (2021)
• (text below equation 31): comment on how surface modes might change results here, or say this will be talked about in conclusions section (I would minimally speculate what is expected to change with surface modes here though).
• (line 370-372): "...while there is a slow down...between seasons, with the largest absolute values..."
• (line 388): against, may want to be explicit about significance of 10^-7
• (line 390): extra space after $0^\circ$
• (line 395): for completeness it is also in Mak et al (2022a), Fig 7 I suppose.
• (line 405): is negative dissipation a problem?
• (line 407): would recommend $\widehat{lon}$, or even better $\widehat{\mbox{lon}}$
• (line 428): how significant are the uncertainties? a quantification would be useful
• (line 438-439): is the issue of "time-scale of response" an issue?
• (line 451): Mak et al (2022a) the one actually intended? (2022b is the prognostic calculation, while the 2022a is the inverse calculation). How are the results here compared to that say?
• (line 467-468): Rewrite? I assume you want "A further advantage of this work is having a simple analytical expression for this dissipation rate that can be applied to ocean datasets" or something similar. It doesn't flow very well at the moment.

The authors make a valiant attempt to estimate the eddy energy dissipation rate, and its spatial distribution, due to relative wind stress. The attempt is both useful and insightful. But I have a number of problems, especially with Section 5, and especially concerning  how to relate continuous stratification to the 2 layer model. Since Section 5 is central to the manuscript, I am recommending a major revision. I set out my comments below, starting with Section 5.

Section 5:

1. Line 367: There is a formula given here for the reduced gravity, g’. The first problem I have with this formula is the term involving the eigenmode, phi_1. The eigenmodes can always be scaled by a constant multiplying factor and since we are not told how the eigenmodes are normalized, the formula for g’, as written, is not sensible. It follows that we need to be told how the eigenmodes are normalized. Nevertheless, I cannot make sense of where this expression for g’ comes from. The fact the value of the eigenmode at z=0 appears is mysterious; the fact it appears squared even more so. Some explanation is required.
2. In the 2-layer model, the wave speed, c, for the first baroclinic mode is related to the reduced gravity by c^2 = g’ H_1 H_2/H where H = H_1 + H_2 and H_1, H_2 are the undistiurbed depths of the upper and lower layers respectively (see Gill(1982), equation (6.3.7)). Presumably the term involving the eigenmode in the expression for g’ above corresponds to the factor H/H_1?
3. Line 366: Here, it is written that H_1 is taken to be the depth of the zero crossing of the eigenmode. Why write this in terms of the minimum of the absolute value of the eigenmode when this is obviously zero at the zero crossing? I also wonder if this is the appropriate choice for H_1? It seems reasonable but might also turn out to be on the deep side? On the other hand, it is clear that the choice of H_1 is of great importance since once H_1 is known, so is H_2 and hence g’, given that c and H are known (see above).
4. Line 369: How is mu computed? We are not told anything about this. On line 186, mu is given by mu = -g’H_2/gH. How does this connect to the mu being used here?
5. I found myself wondering if a simpler model for an eddy than the 2 layer model used here (i.e. allowing the lower layer to be in motion and selecting only the baroclinic model) might be to assume that the lower layer is at rest, corresponding to a projection onto both the baroclinic and barotropic modes. Such a set-up would be equivalent to the 1 ½ layer model. The 1 ½ layer model is recovered from the 2 layer model in the limit that the depth of the lower layer, H_2, goes to infinity. The mathematics for the 1 ½ layer model is simpler than for the 2-layer model and I feel sure the expression in equation (29) would be the same as given in the manuscript with mu = g’/g (since H_2/H tends to one as H_2 tends to infinity). This shows very clearly the importance g’ for determining the dissipation rate. The authors might want to consider briefly discussing the 1 ½ layer model as an alternative to their single baroclinic mode setup.
6. An obvious question that is not addressed is what the authors do about the fact the ocean actually has variable bottom topography. I assume that the eigenvalue problem stated at the beginning of Section 5.1.1 uses H = 4000m? This should be made clear. How do the authors deal with the temperature and salinity data if the local depth is less than 4000m? Again, this needs to be explained.
7. Another issue throughout section 5 and in the figures is how the authors deal with the Rossby radius of deformation near the equator where the Coriolis parameter goes to zero? I assume a band around the equator must be left out. But we do not see this in Figure 7 and neither is anything ever said about it in the text.
8. Line 334: We need to be told what expression is being used for the energy in (29). In particular, we need to be told which term in (17) comes from the thermocline displacement? I also feel confident that it is easy to show that this term dominates the expression in (17) by substituting appropriate values for the parameters. Showing this would justify the use of the “key assumption” made here and make it clear it is not an ad hoc assumption.
9. It seems that the wind speed used in (29) is the seasonal mean wind speed. This also needs to be stated clearly. What about the fact that the winds are not steady? How could this affect the result? At least in (29), wind speed appears linearly (although the speed itself is not linear). In reality, however, the drag coefficient also depends on the wind speed. Indeed, I am reminded of the paper by Thompson, Marsden and Wright (1983, JPO) from 40 years ago…
10. Line 385: Is L_e being used or what is written on line 355?
11. Lines 389-390: To talk about the influence of bottom topography is a bit glib since what is plotted comes from equation (29) and hence must arise from the terms in (29). Maybe R, maybe the wind speed? Of course, it is true that this is where the Antarctic Circumpolar Current takes a turn to the north, consistent with the form drag effect across the Drake Passage sill, so bottom topography may play an indirect role.
12. Line 394: Likewise, here, the authors could (should) be able to say which terms in (29) are making the important contribution.
13. Line 409: The dissipation rates plotted in Figure 9 have to be consistent with Figure 8 by construction!

Other Comments:

1. Lines 20-21: There are lots of examples of how eddies modulate volume transport, going back to the early quasi-geostrophic models of Holland et al., e.g. Holland, Rhines and Keffer (1984), or even Holland (1978), to more recent papers such as Wang et al. (2017, GRL.
2. Lines 33-35: Tandon and Garrett (1996, JPO) were perhaps the first authors to ask what happens to the energy released by the GM scheme. The discussion here reminds me of the backscatter scheme…
3. Line 37: Although Jansen et al. gets referenced, there is no explicit mention of the backscatter scheme in the text?
4. Line 63: The drag coefficient is known to be a function of wind speed – there is no “could be” about it! Think of Large and Pond (1981, JPO) and all the subsequent updates.
5. Equation (6) should not include “A” explicitly since the amplitude is already contained in eta as given by (5).
6. Equation (8): It should be noted that this approximation assumes that that the absolute value of u_a is much bigger than the absolute value of u_g.
7. Equation (10): I assume that the drag coefficient is assumed to be a uniform constant here, independent of u_a?
8. Equation (10): Furthermore, if I understand correctly, it is not W_rel that is being integrated but only the last 4 terms in the expression for W_rel in (9c)?
9. Line 142: It is already clear from the last four terms in (9c) that P_rel is negative. This is because the only term that is not negative-definite is the (u_g)^3 term and this integrates out for a circular eddy. I am also noting here that the absolute value of u_a, plus u_a itself, is always positive or zero.
10. Lines 151-152: Better to write “zero net vertically-integrated flow”.
11. Line 163: Should say that eta_2 is measured positive upward. This is because interface displacement in the 2 layer model is often measured positive downward.
12. As I mentioned earlier, the authors could consider also discussing the simpler 1 ½ layer model in addition to their 2 layer model.
13. Line 186: It should be noted that these solutions for lambda and mu are only valid in the limit g’/g goes to zero – see Gill (1982), Section 6.2. From (13), for the barotropic mode this is obvious because lamda = 1 implies g’ = 0. But it is also true for the baroclinic mode.
14. Lines 198-199: It is not correct to say that the terms in the middle represent the redistribution of energy by nonlinear advection. The pressure work term is also contained in these terms! Think of the energy equation for the equations linearized about a state of rest, e.g. Gill (1982), Section 5.7.
15. Line 215: Surely the higher order boundary conditions on biharmonic viscosity must be used for the first two terms on the right hand side of (19) to drop out on integration? The no normal flow condition is not enough on its own.
16. Section 3.1: I am assuming that the imposed atmospheric wind is uniform, that the drag coefficient is independent of wind speed and that the model domain has closed boundaries? Please be clear about these things! Also, what depth is used for the model ocean, 4000m?
17. Line 260 and thereabouts: It might be helpful to show a plot of the vertical profile of the eddy used for initialization?
18. Section 4: I do not like the title. How about “Verifying the analytic model”?
19. Section 4, first paragraph: It should be stated in the text that this paragraph refers only to the first 150 days of integration, otherwise the impression one gets from Figure 3 is quite different.
20. Lines 318-319: What is the vertical mixing scheme used in the model? Vertical mixing by itself will always increase potential energy – vertical mixing does work against gravity by mixing dense water upwards. Perhaps the greater increase in potential energy due to vertical mixing in the cyclonic case is because the stratification is weaker than in that case?
21. Figure 4: I assume what is plotted is geostrophic relative vorticity, not total relative vorticity? Please make clear.
22. Line 454: The 2 layer model used by the authors is not “linearized”, as written in the text, although the nonlinear terms do not play a role in the analysis.