The manuscript aims at assessing the performance of the Metastatistical Extreme Value Distribution in estimating high quantile of extreme sea level also including future sea-level rise. The topic discussed is relevant and important to improve the resilience of coastal systems facing the effect of a changing climate. However, few aspects discussed in the manuscript need to be revised and discussed more in-depth. The terminology and notation used need also improvement to ensure consistency throughout the manuscript and avoid confusion in the readers.

From the title of the manuscript, the reader expects to read a study about extreme storm surge. However, the study’s objectives (Lines 66-68) refer to extreme sea level. Later on, Line 150, the Authors say that they will investigate the variable h(t) being the sum of tide and storm surge, so sea-level without mean sea level. I would encourage the Authors to clearly state the variable of interest and the variable used when performing the analyses, see also other comments below.

Information regarding MEVD, which is the main method investigated in the manuscript, is limited. The Authors say that this method guarantees “the least amount of a-priori assumption” (line 56). However, the following assumption must be made: F(x,θ) in Eq. 2, the threshold for the ordinary values, the estimation window for parameter estimation, the time-lag to ensure independence between ordinary values. How then is this method the one with the least amount of a-priori assumptions? I suggest clarifying further the advantages of the MEVD compared to the other two methods investigated. Moreover, additional information should be discussed: how the threshold for the ordinary value was selected (line 121 says “as small as possible”); how the 5-year estimation window was selected; why the 30-day lag time for the independence of the ordinary value is so different compared to the values found in the literature (lines 173-179); and how F(x,θ), which turns out to be a GDP (Line 267), is different compared to the classical GDP

I do see the value in implementing the cross-validation procedure to assess the predictability power of the distribution selected as representative of the observations. At the same time, I see the cross-validation as an additional measure of goodness of fit rather than the main one. The NDE only tests if the one quantile associated with the return period Tr of interest is well captured. What about the other quantiles? Is the distribution representative of the entire sample? Also, how the observed quantile h(obs,p) is calculated? Which sample (M,S, or V) is used? The Q-Q plots are mentioned only in the results section and they are only performed for the 30 years in-sample test. In my opinion, the Q-Q plots put the NDE into perspective and should be included as goodness-of-fit method. Also, it would be useful to have them in the main manuscript. I do understand that the space is limited, maybe the Authors could consider including in the main manuscript only the ones related to the MEDV.

In the section Return Period, the definition of Equation 4 needs to be further discussed. Even if the Authors replace (h) with (z-msl), Equation 4 is still the return period of (h), and not the return period of the (z), as indicated by the Authors. Mean sea level (msl) shows a clear linear trend and such trend is recognizable in (z). Similarly, in Equation 5, the distribution G is the distribution of the variable (h) and not the variable (z) as reported in line 341. This has an implication in Figure 5. I assume that the y-axis in Figure 5 “water level” refers to the variable (z). This variable (z) is time-dependent, while in Figure 5 it seems like the statistical properties of (z) are constant. I would have expected something similar to the effective return level plots, to show the effect of sea-level rise. How (msl), which is time-dependent, is added to (h), which is not time-dependent, to derive Figure 5? I suggest clarifying the transition from the analysis on the variable (h), a random variable, to (z), which presents a linear trend due to (msl). I also suggest being more precise with the notation and the terms used throughout the manuscript. It is very difficult to understand the variables the Authors refer to because are often called with many different terms, e.g., total sea level, water level, extreme sea level...

The Authors say that “MEVD proves to be a good model for the extreme sea levels” (line 288) and that “MEVD-based estimates outperform the traditional approaches” (line 301). I do fail to see what the Authors describe. In the QQ-plots Figure S2-6, MEVD in the in-sample analysis has, in general, the highest variability, especially compared to the GEV. In the out-of-sample, MEVD looks better for lower quantiles, but it has quite a large variability for higher quantiles, compared to the other distributions. Overall, it is difficult to quantify which distribution performs best. This is also reflected in the NDE plots, Figure 3, where the differences between distributions are minimal.

Point by point comments:

1. Line 92. Please revise the notation. Pr(Mn<= x) = F(x)^n where Mn is the maximum of a sequence of independent random variable X. See also Coles 2001 (line 415)
2. Line 154. Additional discussion is needed concerning the fact that h(t) can be considered a stochastic variable even though a determinist component is included. Also, a literature review on indirect and direct methods (Line 149) for extreme sea level is missing.
3. Lines 133. The Authors discuss the negligibility of tide-surge interaction. Does this condition hold in the case of Punta della Salute which is located within the Venice Lagoon?
4. How the GDP threshold is selected and tested?
5. It would be very interesting and useful to appreciate the difference between the performance of the distribution functions to see the sample of maxima used for fitting the distributions.
6. Lines 205-209. My suggestion is to revise this paragraph. The terminology is confusing. I believe the Authors here are discussing the variable (z), in which storm surge is a component.
7. Lines 220-221. The Authors say that the tidal and storm components do not change over time as mean sea level. How did the Author check that no trend is detected in the variable h?
8. Section 3: Was the trend test performed only on the annual maxima or also on the samples of maxima used to compute the GPD and the MEVD?
9. Line 281: Storm surge or storm surge and tide?
10. Line 285: what is L?

The paper “Extreme Storm Surge estimation and projection through the Metastatistical Extreme Value Distribution” by Caruso and Marani compare different approaches to perform extreme value analysis of storm surge data at four sites along the European coasts.

The goal of the paper is relevant for the journal and the analysis seem to be correct. Results are well presented and exhaustively explained, I particularly appreciated the quality of the figures.

However, in my opinion there are some aspects the Authors should clarify before the manuscript is deemed suitable for publication. A list of comments is reported below.

• Three approaches are compared in this research, i.e., GEV distribution on annual peak maxima (GEV-BM), the Metastatistical Extreme Value Distribution (MEVD), and GEV distribution on peaks over a higher threshold (GEV-POT). With respect to the latter approach, wouldn’t it be better to rely on a Generalized Pareto Distribution (GPD) when threshold exceedances are considered? As far as I remember, GPD is a derivation of GEV for POT data; as such, is it conceptually correct to test a GEV distribution rather than a GPD on POT data? Please comment on this in the Methods section and/or extend the explanation in the Introduction (e.g. lines 37-39).

• Lines 15-19 in the Introduction. As you speak of “active field” as for the modeling of extreme value probability of occurrence, you could reference more recent works.

• Line 29 in the Introduction. The list of reference is rather long; perhaps it would be enough to cite a few works and the “references therein”.

• Line 48 in the Introduction. You can also cite Solari et al. (2017).

Solari, S., Egüen, M., Polo, M. J., & Losada, M. A. (2017). Peaks Over Threshold (POT): A methodology for automatic threshold estimation using goodness of fit p‐value. Water Resources Research, 53(4), 2833-2849.

• Page 4, Fig. 1. Please reduce the y-axis range for Marseille plot.

• Page 5, line 107. If I understood correctly, “year” in the following line should be replaced with “block”.

• I would swap Section 2.2.1 and Section 2.2.2. First explain how you pre-processed the data, then the distribution used to model them.

• Section 2.2.1. I think you should explain what are the cumulative distributions F you tested for the ordinary values, and which one did you choose.

• Section 2.2.2, line 134. The fact that you neglect the interactions between tides and surges means that gauges are placed in deep waters. Is that true? Please add the respective water depths in Table 1 if such info are available.

• Section 2.2.2. Please number the equations.

• Section 2.2.3, lines 171-173. This paragraph is unclear. Indeed, looking at the correlograms (Fig. S1) it seems that independent events are achieved for no lags. This aspect is crucial so it should be better explained. Correlograms also reveal that tides are relevant (negligible) in Venice and Newly (Marseille and Hornbaek). Perhaps you could comment on this in the paper.
• Section 2.2.3, line 177. Please use consistent tenses throughout the paper when referencing other works. For instance, here you say “Bernardara et al. (2011) adopted”, while previously you use the present tense (e.g. page 6, line 117 or later in the paper at page 8, line 212).

• Page 9, line 239. I do not understand why the return period is expressed for annual maxima (AM). Apologies but I am not familiar with the MEVD, however it is clear that it allows to select multiple events per year. Then, why Eq. (3) is defined with respect to AM data?

• It is not clear the purpose of Section 3.1, given that no non-stationary distributions are subsequently employed. However, if you want to keep it, I suggest adding the confidence intervals of the slopes fitted to the data (and perhaps comment them with respect to the p-values of the Mann-Kendall test).

• Section 3.2, lines 277-279. What does it mean that thresholds are selected “based on local tidal ranges”? This is a pivotal step of the study, please extend the explanation (you could also add it to the Methods section).

• Figures S2-S6. Please use a 1:1 axis ratio. This would help to assess the quality of the fit.

• Figure 5. Levels in the return period plots refer to z or h? I find the terminology rather confusing throughout the whole manuscript, e.g. sometimes you talk about storm surge, some other time about extreme sea levels. Please be consistent.

• Conclusions, line 349. Please specify that MEVD outperforms the other distributions for long enough calibration periods.