See attachment

*General comments

This study examines the abyssal circulation of the Ionian Sea, specifically how vertical coupling between surface, intermediate, and abyssal layers might produce some observed decadal-scale changes in the circulation. The approach is based on sets of four-layer quasigeostrophic models with different stratifications (layer thicknesses) initialized with a Gaussian vortex. Some of the model results are qualitatively consistent with features found in current meter data. The results are very sensitive to the prescribed stratification, and this sensitivity is examined.

The text reads generally well, and most figures are clear. I think there are two major points (detailed below) and several other specific ones that need to be addressed.

I think a central issue in this work is that there is important sensitivity to the choices of layer thicknesses (and densities), as thoroughly examined in the manuscript. This is expected in QG systems with a few layers. Because the primary motivation of the manuscript is to understand decadal changes in the abyssal circulation with the simplest possible model, it would make sense to have experiments with more layers representing the observed density profile (Figure 2) more accurately. This would better approximate the continuously-stratified limit (see, e.g., Gulliver & Radko's 2022 Figure 4 with a ten-layer model based on an observed profile), and help determine at which point the discrete representation of the stratification is sufficiently realistic to reproduce the observed processes.

Another major point is the choice of a single vortex as initial condition. If the motivation is to explain some of the changes in the observed abyssal currents (Figure 1), decaying turbulence experiments with random broadband initial conditions would be more relevant. Even considering that the Intermediate Water eddies found in the Mediterranean have a Gaussian velocity profile, It seems unlikely that changes in individual eddy structure and propagation could be responsible for the observed changes in a real, broad-banded flow with a developed turbulent cascade. Diagnosing the changes in the baroclinic/barotropic energy fluxes in a turbulent four-layer system (with different stratifications like the authors do) could therefore be more helpful to find quantitative links with the observations.

*Specific comments

Figure 1c,d and paragraph starting at line 84: As seen in the hodographs and noted in the text, the amplitudes of the mean flow are similar, but there seems to be much less energy in the subinertial band in the SMO-1 rotary spectra than in GNDT-1. Is this low frequency/mesoscale kinetic energy drop reported in other observations in the literature, and if so, is the reason understood?

Line 47: I think it is important to note here that adding bottom topography with realistic roughness does produce more stable vortices with much longer lifetimes (across different density stratifications), as often observed in real vortices (Gulliver & Radko, 2022).

Lines 113-114: It is said here that compressibility effects were corrected for in the stratification frequency, so I think it would make more sense to have potential density in Figure 2b and in the layer density values reported. Using in situ density rather than potential density is very unusual, and I do not see a case for this choice here. Because only lateral density gradients matter in the QGPV equations, the extra density term due to compressibility effects should make no difference dynamically (apart from any potential numerical error propagation). But I still see no reason to use in situ density instead, especially since it obscures vertical structure in the observed profile (Figure 2b).

Lines 226-234: The relationship between the interfacial deformation radii in a layered QG model (as the one in this work and in Carton et al., 2014) and the modal deformation radii in a continuously stratified system (as in, e.g., Nittis et al. 1993) needs to be clarified here. This is important because the layer-wise deformation radii increase towards the bottom, and that seems to agree with the scale of the circulation features. But the modal deformation radii decrease with increasing mode number, and are associated with vertical modes with more zero crossings (i.e., shorter vertical wavelength). The modal analogues of the layered deformation radii involve combinations of the layer thicknesses (e.g., sqrt(g' * H1*H2/(H1 + H2))/f for the baroclinic mode in a two-layer model). It seems that the abyssal gyres in question are better described in the layered sense as locally equivalent-barotropic features, so I suggest trying to clarify that point.

Figure 5: The axis labels are very difficult to read without zooming in, the font sizes need to be increased.

*Minor corrections and editorial suggestions

Line 21: "in the deep layers process" -> "in abyssal processes", "in deep circulation", or something similar.

Line 31: catching -> observing.

Line 60: wants to -> aims to.

Line 85: to which this study refers -> which this study refers to.

Lines 124, 125, 167: adimensional -> nondimensional.

Line 168: was -> were.

Line 169: be of course worked around adding -> of course be worked around by adding.

Line 200: oscillation -> oscillations.

Line 253: among -> between.

Line 276: can be seen -> it can be seen.

Lines 331, 348: vortexes -> vortices.

Line 334: below ones -> ones below.

Lines 349-350: This last sentence is a bit confusing, perhaps a rewrite?

Line 358: the perturbation changes however the modal structure -> however, the perturbation changes the modal structure.

Line 359: Repeated closing parenthesis.

Line 361: emerges a symmetry in the modal structure around the critical value -> a symmetry in the modal structure around the critical value emerges.

Line 363: appear -> appears.

Line 372: variation -> variations.

Line 401: may be is too simplified for a completed representation of all the processes of the ocean water column -> is too simplified to represent the deep potential vorticity propagation considered here.

Line 405: critical thicknesses ratios values -> critical values of thickness ratios

Line 410: layers -> layer

Line 413: created the preconditions for -> favored

Line 414: The acronym IW is not defined, better to spell out internal waves since it is used only one other time (line 99)

Line 432: help -> helped

*Gulliver & Radko (2022) reference: https://doi.org/10.1029/2021GL097686