This paper presents a spatial Bayesian hierarchical model of sea-level extremes and uses it to analyse tide gauge observations along the Finish coastline. Estimates of extreme sea-level event probabilities, which typically are expressed in terms of return levels, are crucial to flood risk quantification. However, such estimates are often subject to large uncertainty owing to issues related to the small sample sizes and large data dispersion typical of tide gauge observations. Furthermore, when using traditional single-site approaches, estimates of event probabilities are only possible at gauged locations. These issues can be partly overcome by exploiting spatial dependencies in extreme sea levels, or simply by pooling information across data sites, which leads to estimates of return levels with reduced uncertainty and allows for estimation at unobserved locations. Despite the advantages of spatial modelling, most studies of sea-level extremes to date have analyses extremes on a site-by-site basis. In this regard, this paper represents a valuable contribution to the literature on sea-level extremes. The paper shows that pooling information across space leads to more robust estimates of event probabilities, though in this study all tide gauge records are relatively long and as a result the single-site model (‘Separate’) is still able to estimate the GEV parameters with high confidence. The benefits of spatial modelling are much larger in regions with short tide gauge records, and this should be more strongly emphasized in the paper. The paper is well written, the methods are valid, and overall the results are interesting. I do not have any major objections to the paper, but I do have some comments and suggestions, as outlined below, that would like to see addressed before the paper is published in NHESS.

General comments:

1. One of the motivations for using spatial modelling is the ability to make estimates at ungauged locations. However, other than in Fig. 2, the paper focuses on estimates at gauged sites and does not sufficiently assess the skill of the Bayesian models at ungauged sites. I would suggest the authors perform an experiment in which they leave one tide gauge out at a time, estimate the GEV parameters at the omitted site, and then compare the result with estimates based on all the data. I would also suggest the authors include a map of gridded estimates of 50-year return levels along the Finish coastline.
2. Another motivation for using a spatial model is the reduction in estimation uncertainty. I would suggest the authors quantify and discuss this reduction in more detail. By which factor is the uncertainty reduced? Figures 6 and 7 already provide a visual indication, but I think the discussion should be more quantitative. Perhaps a figure or a table showing posterior standard deviations for the 50-year return levels is all that is needed.
3. I think that authors should perform an analysis of sensitivity to prior choices, especially for the parameters defining the spline and GP models. It is well known that the GP parameters (standard deviation and length scale) are challenging to estimate. Also, please explain how and why these priors were chosen.
4. please show the posterior estimates (with uncertainty estimates) for all the scalar parameters (and hyperparameters) of the model, either as a plot or a table.

Specific comments:

Extraction of annual maxima. Was the tidal component removed prior to extracting the annual maxima from the tide gauge records?

Equation 7. The Greek letters used to denote the GP standard deviation and length scale are different between the article and the Supplementary Information.

It is unclear to me what the authors mean by ‘empirical estimates’. The estimates from the Bayesian hierarchical models are conditional on the observations, so they are ‘empirical’ too, aren’t they?

Please add either posterior SDs or credible intervals to Table 2.

Discussion: Line ~335. While I agree that it should be emphasized that to quantify flood risk one should include mean sea level changes, I do not think that excluding mean sea level influences is a limitation of your study, rather it is a choice to focus on the storm surge component of sea level. The actual limitation is to assume stationarity, but this is discussed in the next paragraph.

Discussion. Another limitation that is not mentioned is that the Bayesian hierarchical models used in this study assume conditional independence in the likelihood. In other words, they assume that, after accounting for dependence in the marginal GEV parameter, the annual maxima are independent across stations. However, this assumption is unlikely to hold because the stations are geographically close and thus they are going to be affected by the same extreme events, which means that the time series of annual maxima are going to be correlated between stations (what is called ‘residual dependence’). Ignoring residual dependence means that your uncertainty estimates are narrower than they should be (probably only slightly), but other than that it should not significantly affect the estimates presented in the paper. This limitation should be discussed. Calafat and Marcos (2020) provide a way for addressing residual dependence, but I recognize that this is beyond the scope of this paper.

The future “climate” of water levels is one of the core problems for low-lying areas. The manuscript addresses this problem by means of advanced statistical modeling of parameters of extreme value distributions for future water levels and a sort of ensemble projection of extreme water levels and their return periods.

The analysis is theoretically sound, relies on high-quality data sets, has been performed professionally and leads to an interesting set of results. The presentation is clear and well structured, uses correct English and brings enough details for understanding the material.

General comments:

I am thus generally happy to recommend the manuscript for publication.

Before sending to print, however, I recommend to expand the presentation a little bit to cover some aspects that may mislead inexperienced readers and to make a few adjustments that would make the interpretation more exact and the message clearer. The recommended changes and additions only address single wording features and interpretation aspects (most of which are technically acceptable as presented in the manuscript) and do not involve any large changes to the presentation.

A potential trap for some readers may be the interpretation of the limited set of arguments of the Weibull distribution. Even though the authors mention on lines 307–308 that the shown values [of the upper threshold for the argument of the reverse 3D Weibull distribution] should not be interpreted as actual limits for the sea level, I would recommend commenting the related aspects in more detail to make the situation clear. There are two aspects worth of mentioning.

Firstly, the limited region of validity of the 3-parameter Weibull distribution could be interpreted differently. On the one hand, there is a temptation to think that this distribution provides the final truth about some properties of the described processes. On the other hand, the existence of this kind of threshold is not really physical and could be interpreted as showing that the entire GEV approach loses its validity near and behind this threshold.

Secondly, the set of block maxima may contain elements of different water level “populations” of the Baltic Sea. The reason is the well-known property of the Baltic Sea: its water volume may increase or decrease considerably for several weeks by water exchange through the Danish straits. The “population” of the background water level of the Baltic Sea roughly follows a Gaussian distribution whereas the local storm-driven surges roughly follow an exponential distribution [Soomere, T., Eelsalu, M., Kurkin, A., Rybin, A., 2015. Separation of the Baltic Sea water level into daily and multi-weekly components. Continental Shelf Research, 103, 23–32, doi: 10.1016/j.csr.2015.04.018]. It may thus easily happen that the block (annual) maxima do not necessarily come from the same distribution. In this case the GEV distribution is just a passable approximation of the distribution of the block maxima and nothing more. It may easily be that the large scatter of the threshold of question is a reflection of this feature.

In this sense it is better to remove the conjecture “From theoretical perspective, this suggests that there might be an upper limit that the sea level extremes can reach along the Finnish coast” on lines 369–370 from the manuscript and also to modify the sentence “This also suggests that the hierarchical models can be used to estimate theoretical upper limits of the extremes of short-term sea level variations along the Finnish coast” to make sure that the unexperienced readers are not mislead.

Specific comments:

The Abstract seems too long, e.g., the sentence on lines 3–5 could be removed without any loss to the message and the material on lines 11–13 could be made more compact and smooth.

Line 21: probably “associated WITH” or similar.

Lines 23–24: even though the increase in the mean sea level has exceeded the global average during the past 50 years in the Baltic Sea in many locations, there are opposite examples, e.g., the sea level on the Latvian shores [Männikus, R., Soomere, T., Viška, M. 2020. Variations in the mean, seasonal and extreme water level on the Latvian coast, the eastern Baltic Sea, during 1961–2018. Estuarine Coastal and Shelf Science, 245, Art. No. 106827, https://doi.org/10.1016/j.ecss.2020.106827]. This feature very shortly reflected in (Weisse et al., 2021) and may easily be overlooked. Also, it seems to have local character.

Line 33: it is recommended to insert a reference to the analysis of meteotsunamis in the study area even though such a reference appears later.

Lines 38–39: while piling up water in the ends of the Bay of Bothnia and Gulf of Finland for sure is one of the main reasons for very high water level in these locations, the role of piling and emptying the entire subbasin is probably minor there compared to harbor-type oscillations. See, for, example [Jonsson, B., Döös, K., Nycander, J., Lundberg, P. 2008. Standing waves in the Gulf of Finland and their relationship to the basin-wide Baltic seiches. Journal of Geophysical Research-Oceans, 113 (C3), C03004, doi: 10.1029/2006JC003862]. Still, this effect seems to be a decisive one in some other basins, such as the Gulf of Riga [Männikus, R., Soomere, T., Kudryavtseva, N. 2019. Identification of mechanisms that drive water level extremes from in situ measurements in the Gulf of Riga during 1961–2017. Continental Shelf Research, 182, 22–36, doi: 10.1016/j.csr.2019.05.014.].

Line 61: “However, they did not consider spatial dependencies explicitly in their analysis” is ambiguous and is better to be removed.

Line 80: consider replacing “extends” by “applies”.

Line 116: “which should reduce the correlation between the annual maximum values.” is of course correct but this operation most likely almost totally removes this correlation.

Line 139: What is the meaning of the plus sign at the end of square brackets?

Line 141: consider replacing “y has bounded upper tail” (that is mathematically nonsense for an argument) by perhaps a longer explanation that the GEV distribution function is only defined until a specific value of y which is often associated with the theoretical maximum or minimum value of the process under consideration.

Line 173 and in several locations below: the simple use of “Common” (or similar) makes reading fairly complicated. Consider using “The COMM version/approach/ model” etc., e.g., as on line 246.

Line 199: Do you have a specific reason for using norm when evaluating the expression in square brackets?

Line 234: consider treating “elpdLOO” as a variable, e.g., $elpd_{LOO}$ unless you have reasons for using the text mode. Anyway, unify the use of $P_{LOO}$ in Table 1 and as text on line 243.

Line 248: probably “location and scale PARAMETERS” are meant.

Line 253: it is not recommended to start the sentence from a symbol or expression.

Line 267: consider expanding the expression “Weibull-type distribution” towards explanation that the GEV approach uses so-called reversed 3-parameter Weibull distribution (and not, e.g., the 2-parameter Weibull distribution that is common in the description of wind speed, wave heights, etc.). Just to make clear the scene for the reader.

Line 269–270: “they used ocean model output instead of observations in their analysis” is only partially true. They also used measured data from five locations and noted strong variations in the shape parameter depending on both the particular location and the method for evaluation of the parameters of the GEV distribution.

Line 276: “separate fits”: see comment to line 173.

Section 6 Conclusions: it is recommended to remove the short names of scenarios from the text in order to make the section readable on its own.

References:

Coles 2001/2004 is missing from the list