• L4: I would suggest adding "different" to the realistic initial conditions to make it more clear. (as mentioned in L63)
• L22: Is it possible to have riverine boundary impact for the intrinsic effects? Obviously atmospheric boundary is the most important, but if the riverine boundary could have an impact on the intrinsic effects. If there is large riverine input, it might lead to a different intrinsic effect.
• I believe L24-25: The sentence regarding the mesoscale eddies and flow instabilities is breaking the flow of the paragraph regarding Hasselmann's study. I would suggest moving this sentence to/towards the end of the paragraph. Maybe connect with energy cascade of meso to submesoscale eddies/flow.
• L49-50: What is the basis of Tang's study in the scale/grid resolution? It would be helpful to mention that the capability of the high resolution ocean models to resolve the subgrid scale processes compared to the course resolution models.
• L80 mentions 0.1 degrees resolution for ECMWF. Coppini et al. (2023) mentions the same resolutions 0.125 degrees before 2020, and 0.1 degrees after 2020. Given the simulation period of Jan 2016 to Oct 2020, is there a possibility of a mismatch in the resolution?
• L85-90: It is not very clear that the initial condition is used from the previous simulation (Assuming to my understanding). If this is the case, wouldn't the model carry some of the intrinsic variability through initial conditions to the next simulation? I believe this should be clarified.
• L164: I think one of the most important (albeit expected) results of this study is this line. I would discuss or emphasize this result more.
• L168-170: should be rewritten to make it more clear. There seems to be a missing word or two.

• I would argue in some part of the manuscript the number of ensemble simulations. Tang uses 4 simulation ensemble and Penduff uses 50 as large ensemble. It would be helpful to argue how the number 20 came up for the ensemble simulations?
• Depending in this how many ensemble simulations would make a difference to be able to identify the intrinsic variability?

• In general it is an important first step analysis towards understanding the intrinsic variability in the Mediterranean Sea. In addition, the tides would definitely add an interesting approach to the study and the results as mentioned in L215-. Overall all it is a good manuscript and I would recommend it for publication after the minor revisions mentioned.

In the paper by Benincasa et al. 20 simulations of the operational forecasting system of the Mediterranean Sea are used, through an ensemble approach, to assess the internal/intrinsic ocean variability. It is shown that such a variability is associated with the mesoscale activity and that, with the exception of the Adriatic Sea and the Gulf of Gabes, is larger than the response to surface forcing in all the Mediterranean Sea. Internal variability has a clear season cycle for temperature in the surface layers while, for marine current velocities, it is always dominant and largest at the surface.

The paper is clearly written and well organized and I have only a major concern related to the model resolution and the ability of resolving mesoscale instabilities. The model horizontal resolution is about 4km (1/24 deg) which may be not enough for a full development of mesoscale instabilities: in many Mediterranean areas the Rossby radius (Rd) is less than 8km (i.e. less than 2 model grid cells, see Fig.1 of Beuvier et al. 2012). This seems to be confirmed also by the Rd estimates provided in Grilli and Pinardi (1998) that in some cases are even smaller and closer to the model grid resolution. One may argue that the significant presence of internal variability found in this paper in some areas (e.g. the southern parts of the basin, see Fig.5) is just imputable to the model ability to fully resolve there mesoscales features. Indeed the southern parts of the basins are characterized by larger Rd values (see always Fig.1 of Beuvier et al. 2012). I suggest that the authors discuss such an important limitation of their methodological setup and the sensitivity of their results to horizontal resolution. As the authors themselves report at L48-49, the higher is the horizontal resolution the larger is the intrinsic variability.

A process of revisions is suggested to address also the following minor concerns:

1. L9 “probably”: it should not be so hard to assess whether or not the peak at 30m is really connected to the thermocline formation. Monin-Okubov depth?
2. L42-43, just a curiosity: is there any quantification of the role of submesoscales in setting up intrinsic variability?
3. L77, no tides: this is also reported in the conclusions at L215-216 and may represent an important future extension. But tides are important also in other areas apart Gibraltar (e.g. in the North Adriatic) where the intrinsic variability of this study may be underestimated
4. L96-97: is there a specific reason to stop at 1000m and not to perform the analysis for deeper layers?
5. L108 \sigma_A formula: if \tau is the chosen period I am not sure I fully understand why t starts from 1 and there is a minus 1 in the denominator. My guess is that the total number of discrete timesteps making up the whole period should have been used instead.
6. L132: similar considerations on the \sqrt(\sigma^2) and RMSE formulae as the point above
7. L231: typo Eurosea project

Papers cited in this review and not present in the paper:

• Beuvier J., Béranger K., Lebeaupin Brossier C., Somot S., Sevault F., Drillet Y., Bourdallé-Badie R., Ferry N., Lyard F. (2012). Spreading of the Western Mediterranean Deep Water after winter 2005: Time scales and deep cyclone transport. Journal of Geophysical Research: Oceans, 117 (C7)
• Grilli F., Pinardi N. (1998). The computation of Rossby radii dynamical processes of deformation for the Mediterranean Sea, MTP News 6 (4)