Contribution of high-frequency (T < 2 h) sea level oscillations to the Adriatic sea level maxima

Sea level time series of up to 17.5 years length, recorded with a 1 min sampling interval at 18 tide gauges, evenly distributed along the eastern and western coast of the Adriatic Sea (Mediterranean), were analysed in order to quantify contribution of high-frequency sea level oscillations to the positive sea level extremes of the Adriatic Sea. Two types of sea level extremes were defined and identified: (1) residual extremes which are mostly related to storm surges and (2) high-frequency (T < 2 h) extremes, strongest of which are meteotsunamis. The detailed analysis of extremes led to the following conclusions: (1) high-frequency sea level oscillations can dominate positive sea level extremes; (2) even when not dominating them, high-frequency oscillations can considerably contribute to extreme sea levels; (3) contribution of high-frequency oscillations to total signal is governed by a combination of bathymetry and atmospheric forcing, resulting in the strongest high-frequency oscillations over the middle Adriatic; (4) residual extremes mostly happen from October to January when they are also the strongest, while high-frequency extremes spread more evenly throughout the year, with the strongest events peaking during May to September; (5) tide gauge stations can be divided into three distinct groups depending on the characteristics of high-frequency oscillations which they record. Conclusively, both low-frequency and high-frequency sea level components must be considered when assessing hazards related to sea level extremes, implying that availability and analysis of 1 min sea level data are a must.


Introduction
The rise in sea level due to climate change has led to an increase in number of floods and extreme sea level events that affect coastal ecosystems and communities (Vitousek et al. 2017;Taherkhani et al. 2020). These extremes have become more common and more violent during the end of the 20th and beginning of the 21st century (Woodworth et al. 2011;Sweet et al. 2014), thus sparking significant interest of the scientific community and general public. Numerous studies of extreme sea levels, as well as the processes that generate them, have been carried out (Hamlington et al. 2020;Ferrarin et al. 2022) but only a handful incorporated statistical analysis of sub 2 h oscillations Dusek et al. 2019;Williams et al. 2021;Zemunik et al. 2022). Lack of statistical analysis of shorter-period oscillations is primarily since 1 h sampling of sea level has been a norm in oceanography for a long time. Sea level time series sampled with a 1 min time step became widely available only recently (Zemunik et al. 2021), thus allowing for assessment of contribution of sub-hourly sea level oscillations (2 min < T < 2 h) to total flooding heights (Suursaar et al. 2006;Heidarzadeh and Rabinovich 2021;Pérez-Gómez et al. 2021;Medvedev et al. 2022). This sampling step is, however, still too long to study wind waves (T < 20 s), and shorter period infragravity waves (25 s < T < 300 s).
The Adriatic Sea, which is a semi-enclosed basin in the Mediterranean Sea ( Fig. 1), has a long history of sea level monitoring starting from the 18th century onwards . A lot of phenomena (sea level trends, decadal variability, tidal oscillations, storm surges, basin-wide and local seiches, wind waves, etc.) have been topics of research in the Adriatic Sea. The Adriatic Sea is prone to floodings due to storm surges, tsunamis and wind waves. The term tsunami incorporates different types of tsunamis: (i) seismic tsunamis (earthquake generated); (ii) landslide generated tsunamis; (iii) volcanic tsunamis; (iv) tsunamis generated by massive explosions; (v) tsunamis generated by meteorite fall; and (vi) meteorological tsunamis (meteotsunamis) (tsunamis of meteorological origin) (Tappin 2017). In the Introduction, we discuss two types of tsunamis which are known to cause significant flooding in the Adriatic Sea: seismic tsunamis and meteotsunamis, while in the rest of the paper, we focus only on meteotsunamis. Published analyses of the Adriatic storm surges, seismic tsunamis and meteotsunamis include case studies (e.g. Orlić 1980 for meteotsunamis, Maramai et al. 2007 for seismic tsunami, Međugorac et al. 2015 for storm surges), statistical analysis (e.g. Međugorac et al. 2018;Ferrarin et al. 2022;Šepić et al. 2022a), catalogues (Tinti et al. 2004, andPasarić et al. 2012 for seismic tsunamis; Orlić 2015, andŠepić et al. 2022b for meteotsunamis) and numerical modelling studies (e.g. Šepić et al. 2015a;Bajo et al. 2019).
Storm surges usually occur in autumn (with a peak occurrence in November) (Lionello et al. 2012) during episodes of strong sirocco wind which blows water from the south to the north of the Adriatic Sea. Most notable storm surges occur in Venice which is positioned at the very north of the shallow northern Adriatic and has the opening of the bay facing sirocco ). The highest measured hourly sea level in Venice (Punta della Salute tide gauge) occurred on 4 November 1966 reaching 167.5 cm above the 65-year mean sea level (Šepić et al. 2022a) and resulting in flooding of most of the historic town centre (Canestrelli et al. 2001). Storm surges can also happen during episodes of sirocco wind blowing over middle and southern areas of the Adriatic and combining with strong northeasterly bora wind over the northernmost part of the basin (Međugorac et al. 2018;Ferrarin et al. 2021). Apart from Venice, analyses of flood heights during storm surges were done for other coastal towns in the Adriatic (Šepić et al. 2022a). The highest hourly values include 182.8 cm above the 65-year mean value in Trieste, 114.9 cm in Rovinj, 111.5 cm in Bakar, 85.4 cm in Split and 74.0 cm in Dubrovnik (for location of stations see Fig. 1). In addition, distribution of the northern Adriatic storm surges has been estimated for present climate and for future climate (Conte and Lionello 2014). Storm surges have also been extensively numerically modelled (Međugorac et al. 2018;Bajo et al. 2019) and are also forecasted for the northern Adriatic .
Meteotsunamis are long ocean waves generated by the atmospheric processes such as atmospheric gravity waves, pressure jumps and frontal passages and further intensified by Fig. 1 Map of the area with indicated positions of tide gauges (white circles with black outline). Blue circles mark locations and amplitudes (mostly determined from eyewitness reports and photos) of historic meteotsunamis (Šepić et al. 2022b), and red circles maximum sea level heights (from mean sea level to maximum) related to storm surges (obtained from hourly measurements) (Pérez- Gómez et al. 2022;Šepić et al. 2022a, b). Position of the Adriatic Sea relative to the Mediterranean and data availability are given in the insets. Division of the Adriatic Sea to the northern, middle, and southern Adriatic is also depicted different mechanisms: Proudman or Greenspan and shelf and harbour resonances (Monserrat et al. 2006). Meteotsunamis are observed all around the world, especially at the so called meteotsunami "hot spots", such as the Balearic Islands (Spain), Nagasaki Bay (Japan), the Baltic, the Black Sea, and the Adriatic Sea (Monserrat et al. 2006;Rabinovich 2020;Vilibić et al. 2021b). The Adriatic Sea meteotsunamis started to attract scientific attention after the Great Vela Luka Flood of 1978 (Orlić 1980). With estimated throughto-crest height of short period (T < 2 h) sea level oscillations of ~ 6 m, the Great Vela Luka flood, remains, up to this date, the strongest known World meteotsunami and the event itself represents one of the turning points in global research of meteotsunamis (Vilibić et al. 2021a). Research of the Adriatic meteotsunamis includes numerous case studies (Vilibić et al. 2004;Šepić et al. 2009a, b, 2016, comprehensive overviews of events (Vilibić and Šepić 2009;Vilibić et al. 2021a) and numerical modelling efforts (Vilibić et al. 2004;Orlić et al. 2010;Denamiel et al. 2019;Bubalo et al. 2021). Further on, a catalogue of the Adriatic meteotsunamis has been created by Orlić (2015). The published version of the catalogue contains 21 events. However, regularly updated on-line version of the catalogue now contains 36 events up to the year 2021 (Šepić et al. 2022b).
Tsunamis, both of meteorological and seismic origin, have been observed in the Adriatic Sea for a long time. The most extensive Adriatic seismic tsunami catalogue is the Italian tsunami catalogue which has events from as early as the 1st century (Tinti et al. 2004). For older events, information comes from carefully studied bibliographic sources and for newer events from sea level data as well. By expanding the Italian seismic tsunami catalogue, which mostly lists events from the western side of the Adriatic, with events from the eastern side, Pasarić et al. (2012) created a newer catalogue of events after extensive bibliographic research that included events from the 16th century up to the present, resulting with a list of 16 seismic tsunamis which hit either eastern or western Adriatic coast from 1511 to 1979 (Pasarić et al. 2012). Three best-known Adriatic seismic tsunamis that also generated significant damage were the 1627 Gargano tsunami, the 1667 Dubrovnik tsunami and the 1979 Montenegro tsunami (Tinti et al. 2004;Pasarić et al. 2012). The historical meteotsunamis were occasionally wrongly interpreted as seismic tsunamis. The best-known example is the Great Vela Luka flood which was initially considered to be a seismic tsunami (Bedosti et al. 1980;Zore-Armanda 1988).
In addition, wind waves, with periods up to 10 s, also present a significant threat to the Adriatic coastline, especially to the northern part (Lionello et al. 2012). Wind waves are strongest during episodes of sirocco wind but also appear during periods of bora wind (Pomaro et al. 2017). Infragravity waves, with periods from 25 to 300 s, have been observed in the Adriatic and even measured propagating through river delta inland signalling that they could add to a potential flooding event (Melito et al. 2018).
The extreme sea level events caused by the discussed phenomena may pose a serious problem to various areas of human activity. Moreover, special sites of the UNESCO heritage (e.g. historical centres of Venice, Split, Dubrovnik; UNESCO list: https://whc.unesco. org/en/list/), that form a core of some marine cities, are also under the impact of these events (Fig. 1). So far, over the Adriatic but also globally, different types of positive sea level extremes have almost exclusively been studied separately. Nonetheless, during recent years case studies combining more types of extremes started to emerge: examples include studies of Venice flood of November 2019 ) during which hazardous heights were reached due to both storm surge (long period) and meteotsunami (short period), hurricane-generated floodings in Japan (Heidarzadeh and Rabinovich 2021;Medvedev et al. 2022), storm Gloria in the western Mediterranean (Pérez- Gómez et al. 2021), storm Gudrun in Gulf of Finland (Suursaar et al. 2006) and others. All of these analyses are case studies relating to single events and are not giving a true statistical sense of importance and contribution of high-frequency sea level oscillations to flooding events. The aim of our paper is to overcome this lack for the Adriatic Sea and to evaluate contribution of high-frequency sea level oscillations (2 min < T < 2 h) to the total Adriatic positive sea level extremes, resulting with an example-pioneering study which could be expanded for other global locations.
Paper is structured as follows: in Chapter 2, available data and analysis are presented; in Chapter 3, results of analysis are given; Chapter 4 brings the discussion and Chapter 5 the conclusion.

Materials and methods
We have analysed sea level data measured with a 1 min time step, originating from 18 tide gauges spread evenly along the eastern and western Adriatic coast (Fig. 1). Six tide gauges are in the northern Adriatic, nine in the middle, and three in the southern Adriatic. Division of the Adriatic Sea to the northern, middle and southern Adriatic is based on the bathymetric properties ( Fig. 1): the northern Adriatic encompasses shallow northernmost parts, from the Po river mouth to the north of the 270 m deep Jabuka Pit; the middle Adriatic spans area from the Jabuka Pit to the northwestern outskirts of the South Adriatic Pit; and the southern Adriatic spreads further on, to the Strait of Otranto. Length of available time series is from 3 up to 17.5 years with the longest records available for the Croatian stations of Rovinj, Bakar, Zadar, Split, Ploče and Dubrovnik ( Fig. 1; inset).
Data from the Italian stations were measured by the Istituto Superiore per la Protezione e la Ricerca Ambientale (ISPRA) and provided via the internet site: http:// www. ioc-seale velmo nitor ing. org/. ISPRA stations are Trieste, Venice, Ancona, San Benedetto del Tronto (also referred to as S.B. del Tronto in the paper), Ortona, Isole Tremiti (I. Tremiti), Vieste, Bari and Otranto. Data from the Croatian stations were measured and provided by three Institutions: for Bakar ) by the Department of Geophysics (Faculty of Science, University of Zagreb), for Rovinj, Zadar, Split, Ploče and Dubrovnik by the Hydrographic Institute of the Republic of Croatia (Split), and for Vela Luka (V. Luka), Sobra and Stari Grad (S. Grad) by the Institute of Oceanography and Fisheries (Split) (link to data repository: http:// faust. izor. hr/ autod atapub/ mjesu stdoh vatpod? jezik= eng). Data series from each tide gauge were quality controlled-with quality control including automatic spike detection procedures and additional visual checking of series. All nonphysical spikes and nonphysical outliers were removed. Means of data series were calculated and subtracted from the original series. The sea level time series were de-tided using MAT-LAB T-TIDE package (Pawlowicz et al. 2002) removing seven strongest tidal constituents (Međugorac et al. 2015) in the Adriatic (K1, O1, P1, K2, S2, M2 and N2). The de-tided signal (residual) was filtered using a Kaiser-Bessel window of 2 h length to separate highfrequency (HF) (T < 2 h) and low-frequency (LF) (T > 2 h) component of time series (e.g. Thomson and Emery 2014). Both LF and HF components encompass various types of sea level oscillations. LF processes are, starting from the shortest: storm surges and other synoptic scale oscillations, planetary waves, seasonal oscillations, annual and interannual oscillations and trends (e.g. Šepić et al. 2022a). HF processes are, starting from the shortest: wind waves and swells, infragravity waves, local seiches, and long ocean waves including meteotsunamis and seismic tsunamis (e.g. Medvedev et al. 2022). Due to their short periods, wind waves, swells and shorter period infragravity waves are not detectable 1 3 in 1 min time series; further on, no seismic tsunamis were recorded in the Adriatic Sea during the investigated period. Our HF series, thus, can contain longer period infragravity waves (T > 2 min), local seiches and other long ocean waves of which the strongest ones are meteotsunamis.
As the first step in the analysis, we defined and extracted two types of extreme positive sea levels episodes at each station: (i) extremes of the residual signal-"Residual extremes" from now on; (ii) extremes of the HF series-"HF extremes" from now on. We define Residual extremes as points in time during which residual signal surpasses its 99.85 percentile value (different value at each station) (Fig. 2, Table 1). We considered all points above the given threshold within a continuous time interval to belong to the same episode. Since maximum of residual signal does not have to coincide with the high tide, Residual extremes may or may not be flooding events (Pugh and Woodworth 2014). We define HF extremes as points in time during which HF series (T < 2 h) attain values higher than their respective 99.993 percentile (Fig. 3). Again, we consider all points above the given threshold within a continuous time interval to belong to the same episode. We have chosen the threshold percentiles (99.85 for Residual extremes, and 99.993 for HF extremes) so that, for most stations, 2.5-4.5 extreme episodes per year (in median, for all stations, 3 episodes per year) are extracted for both types of episodes. This is in line with number of yearly events (1-5) commonly chosen for analysis of return periods of sea level extremes (Tsimplis and Blackman 1997; Marcos et al. 2009;Šepić et al. 2022a). Consequently, at tide gauge stations which have the longest data sets and have fewer data gaps, there are more episodes (Tables 1 and 2).
If two, or more, episodes happened within a span of three days, we assumed that they were dependent and that they are the same episode. We have chosen this restriction because  1 3 of the Adriatic-wide seiche that has a fundamental period of 21.1 h . The active seiche is usually generated during a storm surge (Raicich et al. 1999;Vilibić 2000), decaying to 1/e of its initial amplitude after 3.2 days (Cerovečki et al. 1997), if no additional forcing is present. Thus, oscillations within a three-day period are likely to be dependant, especially when it comes to Residual extremes. HF extremes occurring within a 3 d period might on the other hand be independent-nonetheless, to be consistent, we have decided to apply the same 3 d limit for HF extremes as well. We calculated power spectra following a common methodology (Thomson and Emery 2014). For each tide gauge and each type of extreme episodes (Residual and HF extremes), we calculated power spectrum of each episode from the residual series. After the calculation for each episode, we averaged power spectra over all episodes (measured at one tide gauge) of one type (Residual or HF extremes). We calculated background spectra for each station using sea level data from the selected calm periods. By evaluating event and background spectra and comparing them, we were able to assess energy increase in sea level oscillations during two types of episodes. We set the number of data points used for estimation of background spectra to 2 13 (8192 points, i.e. ~ 136 h), and for both episode types to 1536 points (~ 26 h), centred around the episode. We used half-overlapping Kaiser-Bessel window, for both episodes and background spectra, with 2 9 (512 points, i.e. ~ 8.5 h) length, because we were primarily interested in short period (i.e. HF) component of oscillations. This choice led to 62 degrees of freedom for background spectra (same for each tide gauge) and a different number of degrees of freedom for each station and for each type of episode (from 120 for Vela Luka to 660 for Ploče for Residual extremes; and from 20 for Isole Tremiti to 640 for Zadar for HF extremes). The different number of degrees of freedom is related to different number of episodes per station. Finally, we estimated return values for  periods of 2, 5, 10, 50 and 100 years for both types of extreme episodes using generalized Pareto distribution (Marcos et al. 2009) applied to data from stations that have relatively long time series (~ 17.5 y): Dubrovnik, Ploče, Split, Zadar, Bakar and Rovinj. It should be noted that, throughout the manuscript (unless otherwise indicated) term "height" refers to sea level above the mean sea level.

Residual extremes
Selected parameters, related to Residual extremes, are given in Table 1 Average contribution (over all episodes) of the LF and the HF oscillations to Residual extremes is shown in Fig. 4a. Average LF contribution is highest over the northern Adriatic (Venice, Trieste, Rovinj and Bakar) from where it decays southward ranging, from 86 cm at Trieste to 39 cm at Otranto. Average HF contribution is highest over the middle Adriatic, both along its western and eastern coast, reaching values up to 23 cm (at Vela Luka station). The lowest average HF contribution is found at Dubrovnik reaching only 3 cm. Low average HF contributions are also observed over the northern Adriatic stations reaching values of 6 cm at Venice.
Next, we look at the strongest Residual extreme episodes (not the same episode at each station; Fig. 4b; Table 1 Finally, we take a look at the Residual extreme episodes (not the same episode at each station) during which amplitude of the HF component was the largest (Fig. 4c). Among the northern Adriatic stations, only Venice had a residual value (136 cm) higher than average residual value during Residual extremes, whereas at the other northern Adriatic stations (Trieste, Rovinj and Bakar), this value was lower than average residual value (65-75 cm) (comparison of Fig. 4a and c). Contribution of the HF component was smallest at Venice (10%; 14 cm), and larger than 50% at the remaining northern Adriatic stations, reaching heights from 37 to 45 cm. Shifting the focus towards the middle Adriatic, we notice that, at most stations, values of the residual signal were larger than values of the average residual signal, with the exception of Zadar, Split, Vieste and Stari Grad. Over all of the middle Adriatic stations, aside for Isole Tremiti (10 cm; 20%), contribution of the HF component was above 50%, reaching values of 91% at Vieste, 89% at Stari Grad, 79% at Ortona, etc. Over the southern Adriatic, contribution of the HF component dropped to smaller values, reaching 12 cm (27%) at Dubrovnik and 16 cm (20%) at Bari. At Otranto, the HF component again dominated the signal, with its value reaching 35 cm (76%).
Comparison of the three subplots of Fig Contribution of the HF oscillations to Residual extremes at each station is further detailed in Fig. 5 in which we show box plots of HF contribution to Residual extremes, and box plots of height of HF signal during Residual extremes, with outliers defined as plausible values higher than q3 + 1.5 × (q3-q1), where q1 and q3 are the 25th and 75th percentiles of the sample data. For a clearer overview, we grouped stations according to the median value of the HF contribution into: (1) the low HF contribution group (the median HF contribution lower than 10%) "Group 1", (2) the moderate HF contribution group (the median HF contribution from 10 to 20%) "Group 2" and (3) the high HF contribution group (the median HF contribution above 20%) "Group 3". Stations situated in the northern Adriatic (Venice, Trieste and Rovinj) and three stations from the middle and southern Adriatic (Isole Tremiti, Split and Dubrovnik) belong to Group 1. Stations falling into Group 2 are: Zadar, Bakar, Ploče, Ancona and Bari. The border station between Group 2 and Group 3 is Stari Grad. According to the median value of the HF contribution (18.9%) Stari Grad should belong to Group 2, but its total appearance (maximum, outliers) is more like stations of Group 3. That is why it is listed in Group 3 together with San Benedetto del Tronto, Sobra, Otranto, Vela Luka, Vieste and Ortona. It should be noted that at some stations extreme events (recognized as outlier values denoted with red crosses in Fig. 5) have maximum values that are double or triple of the median value, i.e. at Stari Grad and Vela Luka. At some other stations, outliers are also significantly higher than maximum values, e.g. at Zadar and Vieste. It can be noticed again that at most stations (aside for Dubrovnik, Isole Tremiti, Venice, Split and Bari) there are events during which the HF contribution to Residual extremes is close to or higher than 50%, dominating the residual signal and driving the extreme event.
We can get a further insight into characteristic of different groups from Fig. 6 where we show contribution of the LF and the HF components to Residual extremes, for two stations per each group. For Group 1 (first row) we show stations Dubrovnik and Venice; for Group 2 (second row) Zadar and Bakar, and for Group 3 (third row), Vela Luka and Vieste. As representative stations, we have chosen, for each group, two stations with the longest time series. The one exception is Vela Luka, which has a relatively short time series, but, because it is a location prone to strong meteotsunamis (Orlić 2015;Šepić et al. 2022b), we have chosen it to be one of the representative stations for Group 3. The corresponding plots for all the other stations are presented in Supplementary Information in Fig. S1.
Stations belonging to the same group have a similar distribution of the HF, the LF and the residual signal heights (Fig. 6). For Group 1, there is a clear division between histograms showing the HF height (red hatched) and histograms showing the LF (blue hatched) and the residual (grey) heights, pointing to the fact that Residual extremes of Group 1 are dominated by the LF signal. For Group 2 stations, histograms showing the HF level are "shifted to the right" (in comparison with Group 1 stations) and are closer to the LF and the residual histograms, with some bins overlapping, indicating that Residual extremes of Group 2 can occasionally be dominated by the HF oscillations. Finally, there is a clear overlapping of the HF, the LF and the residual height histograms at stations belonging to Group 3. This overlapping implies that Residual extremes of Group 3 are strongly influenced by the HF signal, as is also suggested in Fig. 4.
Residual extremes occur during the whole year, but the peak of their appearance and strength is during autumn months (November and December) (Fig. 7). Their seasonal distribution is, expectedly, the same as seasonal distribution of storm surges over the Adriatic Sea (Lionello et al. 2012). The contribution of the HF oscillations to Residual extremes is highest from May to September, when number of Residual extremes is lowest (Fig. 7). Nonetheless, it is evident that HF oscillations can contribute significantly to Residual extremes during the whole year. Seasonal distributions of number of Residual extremes and the HF contribution for each individual station are shown in Supplementary Information in Fig. S2. Individual distributions follow distribution like the one shown in Fig. 7. It should be noted, though, that there is at least one month during the warmer part of the year (April to September), at each station, during which no Residual extremes were recorded. Further on, at Ancona, Trieste, Bari, Isole Tremiti, Venice, Sobra and Split, no Residual extremes were recorded during 4-5 consecutive warmer period months.

HF extremes
From the previous section, the HF sea level oscillations can contribute significantly to Residual extremes. Further on, from previous research work (Fig. 1), it is evident that at many Adriatic Sea locations and, in particular, the middle Adriatic, it is the HF oscillations that represent the real danger when it comes to extreme positive sea levels, as this is a location where strongest Adriatic meteotsunamis occur (Orlić 2015;Šepić et al. 2022b). Thus, we now present results of analysis of HF extremes. As it was explained earlier, HF extremes are episodes during which the HF signal surpasses its 99.993 percentile value. Percentile thresholds reach values from 5 cm at station Dubrovnik, to 24.6 cm at station Vela Luka, showing no clear south-to-north gradient, but some intensification in the middle Adriatic ( Table 2). Number of HF extremes per year ranges from 0.7 at Isole Tremiti to 5.4 at Otranto, with median value of around 3.3 episodes per year, similar to what was obtained for Residual extremes. Difference in number of episodes at various stations can be related to definition of events: to satisfy criteria of independent events, we consider all HF episodes occurring within 3 days to represent one HF episode-at those stations where strong episodes are closer to each other, we are bound to have less episodes per year. For example, the station Isole Tremiti has all points exceeding a given threshold grouped around the two events meaning that all of them fall into two defined episodes.
Additional basic variables including values related to the maximum HF episode at each station are given in Table 2. At each station maximum HF event occurred at different date, except for the 9 July 2019 event which occurred simultaneously at Stari Grad and Ortona, and the 2 August 2019 event which occurred simultaneously at Ancona and Zadar. This opposes behaviour of Residual extremes which dominantly appeared  (Table 1) and further suggests that while Residual extremes are basin-wide events, HF extremes are much more localized, often peaking at one individual location (e.g. Šepić et al. 2022b). Maximum residual value, during the maximum HF extreme, was recorded in San Benedetto del Tronto reaching 87.8 cm. During that event, the LF component contributed to 44.2 cm and the HF component with 43.6 cm to the maximum residual level. Ortona, with the HF height of 60.7 cm, had the highest HF value during all HF extremes. During that event the residual signal was 75.9 cm meaning that almost 80% of the entire signal was due to the HF component. It should be noted that 10 of HF extremes given in Table 2 match the dates (but not necessarily locations) of the known Adriatic meteotsunamis listed in the Adriatic catalogue (Šepić et al. 2022b). Adriatic catalogue of meteotsunamis contains only events during which through-to-crest height of HF oscillations was above 100 cm; implying that their height, as we define it herein (from the mean sea level to maximum), was at least 50 cm.
Presenting similar statistics as in Fig. 4, but using HF extremes, average residual sea level height during HF extremes is much lower than during Residual extremes (Fig. 8a), reaching values from around 21 cm at station Dubrovnik to 54 cm at Vela Luka station, with a clear intensification of average HF oscillations for stations in the middle Adriatic. Also, average contribution of the HF oscillations to the residual signal is close to or larger than the contribution of the LF component-meaning that the residual signal is either evenly dominated by the HF and the LF signal or dominated by the HF signal. Exception is Isole Tremiti station, where the average HF contribution is around 9 cm (21%), and Dubrovnik where it is 6.5 cm (30%).
Among all the HF episodes, the ones with maximum residual signal are presented in Fig. 8b. The maximum residual height at most stations (Fig. 8b) is larger than average residual height of Residual extremes (Fig. 4a), i.e. the largest residual events among the HF extremes are also Residual extremes. Finally, we consider those HF extremes during which maximum value of the HF signals was measured (Fig. 8c). Northern stations (Venice, Trieste, Rovinj and Bakar) exhibited low residual values, reaching 47-73 cm, with the residual signal dominated by the HF component (> 56%). Larger residuals were recorded over the middle Adriatic, reaching values from 35 to around 88 cm, at Isole Tremiti and San Benedetto del Tronto, respectively. Over the middle Adriatic, the HF component contributed to residual signal with 37 to 89%. Moving towards the south (tide gauges at Dubrovnik, Bari and Otranto), the situation is reminiscent of that from the northern stations with a low residual value (reaching from around 21-53 cm, at Dubrovnik and Bari, respectively) and a great percentage of that value coming from the HF component.
The distribution of the HF, the LF and the residual heights during HF extremes is shown in Fig. 9 for six representative stations (the same stations as in Fig. 6). Although there are only six stations chosen, they are good representatives of their respective groups, and distributions for all the other stations are shown in Supplementary Information in Fig. S3. Looking at the distribution of the HF, the LF and the residual height histograms, we again notice distinction among three groups. The Group 1 stations (Dubrovnik and Venice) have overlapping HF and LF histograms suggesting similar values of the two, while residual histograms tend to have higher values. For the stations in Group 2, Zadar and Bakar, the division between histograms of LF and HF values starts to appear: the LF histograms move to the left (reaching lower values), and the HF histogram to the right (reaching higher values), indicating that HF oscillations are becoming more important. Furthermore, residual value histograms spread out more from the HF and LF ones. In the final group, Group 3, the LF histograms move further left, and the HF histograms further right, with the least 1 3 Fig. 8 a Average height of the residual (red circle; Z r ) and the LF (blue circle; Z LF ) time series during HF extremes; b the residual and the LF height during the HF extreme episode in which maximum residual sea level was measured; c the residual and the LF height during the episode in which HF component had its maximum value. Station names are coloured to distinguish different groups defined in Fig. 5: names of Group 1 stations are given in blue, of Group 2 stations in bold black, and of Group 3 stations in red overlapping between the two; the HF and residual histograms, though, overlap the mostimplying that the HF signal is strongest for Group 3 stations, controlling the residual levels.
Distributions of heights of HF component during HF extremes at each station are presented in Fig. 10. The stations generally keep the properties of previously defined groups (Fig. 5), but with some exceptions. Group 1 (stations with median contributions of the HF signal to Residual extremes lower than 10%; Fig. 5) generally contains stations with low median values of HF heights (< 20 cm)-however, at some stations (Trieste and Rovinj) HF heights are distributed more widely reaching outliers which are 2-3 times higher than median values. In Group 2 we find three stations with median HF values higher than 20 cm; and Zadar and Bari with lower median values (< 20 cm at both), with the former containing exceptional outliers. Finally, Group 3 contains stations with the highest median values of HF contribution (up to 35 cm at Vela Luka), but with Sobra and Otranto having median values lower than 20 cm.
We have also examined seasonal distribution of HF extremes. The frequency of occurrence is highest during months May to October, and these are also the months during which HF extremes reach highest values (Fig. 11): medians reach values from 15 cm in April to almost 23 cm in July. Analysis of the seasonal distribution and height of HF   Fig. S4) yielding, for most stations, a similar distribution as shown in Fig. 11. Two exceptions are station Isole Tremiti at which we extracted only two episodes and station Otranto at which there are less HF episodes during summertime than during winter.

Spectral properties of extreme events
We estimated spectral properties for all stations and plotted respective spectra for two stations from each of the three groups (Group 1-low HF contribution, Group 2-medium HF contribution, Group 3-high HF contribution) (Fig. 12). Chosen stations are the same as in Figs. 6 and 9; Dubrovnik and Venice for Group 1; Bakar and Zadar for Group 2; Fig. 12 Averaged power spectra density over: all Residual extremes (red line), all HF extremes (dashed black line) and background periods (grey line). The confidence intervals are given for each spectrum, 95% for Residual and HF averaged episodes and 99% for background. Vertical grey lines denote the times (periods labelled with "T") when power spectra densities intersect for the first and the second (last) time. Periods of significant peaks are given as well (in min) and Vela Luka and Vieste for Group 3. For each station, we show three power spectra, one estimated from Residual extremes, one estimated from HF extremes and one from the background signal (following methodology described in Materials and methods). We note here again that spectra are estimated using residual signal for both types of episodes and for background periods.
Throughout the stations and different groups, there are a few characteristics of the power spectra that are similar. For lower-frequency parts, up until the first vertical line in each plot, sea level oscillations related to Residual extremes contain more power than those related to HF extremes. This is in line with the physical processes responsible for oscillations in each of the extreme episode types. Residual extremes are mostly related to surge events, thus high power is contained in the part of the spectrum where storm surges are observed (from a few hours up to 3-4 days, (Rabinovich 2020)). The first vertical line in each plot marks an average period at which HF extremes become more energetic than Residual extremes. This vertical line, for most stations, is at 1-3 h.
The next part of the spectrum, from a couple of hours down to a couple of minutes, is normally dominated by seismic tsunamis, meteotsunamis, seiches and infragravity waves. Seiches are detected as peaks in the power spectrum with periods ranging from couple of minutes to tens of minutes. For most of our stations, there are spectral peaks which appear in all three types of power spectra (Residual extremes, HF extremes, background), indicating that seiches of certain period are quasi-constantly present in sea level oscillations, but have varying strength. Matching of peaks is seen at all periods at Vieste and Zadar. At Vela Luka and Bakar, the shortest HF peaks (periods shorter than 8.5 and 4 min respectively) are not seen in the background spectra-implying that these oscillations are generated only during more energetic events. At Venice and Dubrovnik, none of broad peaks detected during extreme events are found in the background spectra. Interestingly, at Dubrovnik, Venice and Vieste, there is a point (around the period of 2.5 and 3 min; marked by the second vertical grey line) in the power spectrum at which energy of HF extremes again becomes lower than energy of Residual extremes. These periods are normally occupied by infragravity waves-which are known to develop during storm surges (Hattori et al. 2020) and which are more likely to be detected at stations under greater influence of open sea wind waves. These waves were observed in the Adriatic in front of a river mouths, and deeper inland, during storms as reported by Melito et al. 2020.
On top of the noticed similarities, there are notable differences between selected stations, especially between the groups of stations. Starting with Group 1 we notice that power spectra of these stations (Dubrovnik and Venice) contain mostly broad and notwell defined peaks. Furthermore, these peaks are not present in the background spectra. Similar spectra are obtained for other stations in this group ( Supplementary Information  Fig. S5). Energy is generally increased during both types of episodes, assuming form of white noise over the shortest periods (T < 4 min) likely pointing to the increased infragravity wave activity (Rabinovich 2009). Following terminology from Medvedev et al. (2022) we can conclude that stations within Group 1 are dominantly "white-noise" stations, i.e. stations at which there are no pronounced spectral peaks, and whose spectra at the highest frequencies resemble white-noise spectra, especially during the most energetic Residual extremes, likely due to influence of storm waves, swells and infragravity waves. Group 2 is represented by stations Bakar and Zadar (Fig. 12). These spectra contain numerous welldefined short-period peaks of comparable strength and we consider these stations to be "multi-resonant" (polychromatic) stations ). Bakar has pronounced peaks at around 20, 8, 4 and 3 min, and Zadar at 11, 4 and 3 min. Finally, within Group 3, which is characterized by the strongest HF oscillations, spectra of stations have "resonant" (monochromatic) properties , meaning that one peak is significantly stronger than others. Nice example of such stations is Vela Luka where a peak centred at 18 min clearly surpasses all other peaks in strength, and Vieste at which similar peak is found centred around 12 min. At all other stations from Group 3 there is one such dominant peak, centred at 22 min at Stari Grad, at 18 min at San Benedetto del Tronto, at 21 min at Ortona, 9 min at Otranto and 85 min at Sobra (additional spectra are visible in Supplementary Information in Fig. S5). If we were to group our stations according to spectral properties (as opposed to the median grouping) minor shift between groups would occur, Split station has a multi-resonant (polychromatic) spectrum and would thus belong to Group 2 (instead of Group 1), and Bari and Ancona stations have resonant (monochromatic) spectra and would transfer to Group 3 (instead of Group 2). It should be noted, however, that these three stations are right next to borderlines according to the division into groups presented in Fig. 5.

Discussion
Numerous papers dealing with statistical analysis of extreme sea levels are focused on the LF sea level component (T > 2 h) (e.g. Menéndez and Woodworth 2010;Marcos et al. 2015). There are two reasons for this: the fundamental one is that the deadliest storm surges, associated with tropical cyclones, occur at periods from a couple of hours to half a day (von Storch and Woth 2008); the technical one is that sea level data have traditionally been sampled with hourly time step, which does not allow for analysis of extreme shortperiod oscillations such as seismic tsunamis and meteotsunamis. Nonetheless, the strongest of the HF events have been studied using either higher resolution data (Rabinovich 2020;Carvajal et al. 2021) which mostly started to become available after the 2004 Indian Ocean tsunami (Titov et al. 2005) or digitized versions of older analogue tide gauge charts (Šepić et al. 2015a;Pellikka et al. 2018). Statistical analyses of strength and frequency of extreme sea levels related to the HF oscillations are, however, sparse, because long enough 1 min sea level time series have only recently become available. Nonetheless, initial statistical studies of extreme HF oscillations were made during the past 10 years.  and Zemunik et al. (2022) provide global analysis of the HF oscillations done on time series of 1 to 12 years lengths. There are also some regional (covering stretches of coast longer than ~ 500 km) papers, like Šepić et al. (2015b), Dusek et al. (2019), Williams et al. (2021) or local (covering one individual location) like Ozsoy et al. (2016), providing similar analyses. However, all of the mentioned papers focus exclusively on the HF component of time series, i.e. on events comparable to our HF extremes. This paper represents a pioneering effort when it comes to analysis of contribution of the HF oscillations (T < 2 h) to total sea level extremes.
The first important finding of our study is that the HF oscillations represent an important, and occasionally even dominating contribution (> 50%) to extreme residual sea levels. This result is valid at most of our locations (excluding Venice, Isole Tremiti, Split, Dubrovnik and Bari-although time series at Bari and Isole Tremiti are relatively short, likely impacting the analysis) (Fig. 4). To better understand this contribution, we did an exercise of estimating return values (for periods of 2, 5, 10, 50 and 100 years) of Residual extremes and HF extremes but only for stations with data length of ~ 17.5 years (Rovinj, Bakar, Zadar, Split, Ploče and Dubrovnik; Fig. 13). Expectedly, Residual extremes have higher return values, reaching from 53 cm at Dubrovnik and 96 cm at Bakar for a 2-year return period to 69 cm at Dubrovnik and 132 cm at Bakar for a 100-year return period. HF extremes are similarly distributed, reaching from around 10 cm at Dubrovnik and 40 cm at Bakar for a 2-year return period to 14 cm at Dubrovnik and 54 cm at Bakar for a 100-year return period. Both distribution of Residual extremes return values and distribution of HF extremes return values seem to have south-to-north gradient; with the former related to south-to-north gradient of storm surge heights, and the latter likely occurring by chanceas our stations with the longest time series are such that stronger HF episodes are recorded at the northernmost stations Rovinj and Bakar (recognized as outliers at those stations), and weaker at the southernmost stations Split and Dubrovnik (Fig. 10). What should be noted here is that return values of the HF extremes amount to, depending on station, around 18 to 48 per cent of return values of Residual extremes, pointing again to importance and necessity of studying the HF extremes. If we had long enough time series measured at meteotsunami prone locations, such as Vela Luka and Stari Grad, where the HF oscillations can reach levels of more than 200 cm (Fig. 1), these percentages would be much higher.
There are two key elements governing appearance of sea level extremes: bathymetry and atmospheric forcing, with both elements important for Residual and HF extremes. Over the northern Adriatic, atmospheric forcing and bathymetry act jointly to generate extreme storm surge sea levels. Deep cyclones often pass over the northern Adriatic with related sirocco wind blowing towards the northern Adriatic, causing additional piling of water therein (Međugorac et al. 2015). Effect of shallow depths is superimposed to the atmospheric pressure effect as wind affects sea level more over shallow than over deep areas (sea level gradient is proportional to W/h, where W is wind stress, and h is depth) (Pugh and Woodworth 2014). The HF contribution is, on contrary, much weaker over the northern Adriatic-presumably because bathymetry of bays and harbours does not allow for dramatic amplification of the long ocean waves. This is confirmed by spectral analysis Fig. 13 Return values of (a) Residual and (b) HF extremes. The 95% confidence intervals are presented with whiskers which indicates that most of the northern Adriatic stations belong to Group 1 characterized by the spectra without significant spectral peaks at periods shorter than 2 h. Nonetheless, it should be noted that the 2019 flood of Venice was due to combined effect of storm surge and a smaller scale atmospheric disturbance , and that numerous meteotsunamis in Rovinj higher than 28 cm have been detected (Šepić et al. 2012). It should be also noted here that periods of Venice (2-3 h, Umgiesser et al. 2021) andTrieste (2.7-4.2 h, Šepić et al. 2022a) seiches are longer than what we consider to be HF oscillations in this paper. Also, a study has recently been published for Venice in which contribution of various sea level processes, acting on time scales from 20 min to decades, to total sea level extremes has been estimated (Ferrarin et al. 2022). In this study, authors used 10 min data for 2010-2019 period and considered all sea level oscillations at periods shorter than 10 h to be HF. While envisioning our research, we considered a higher cutoff period (e.g. 6 h) for separation of the LF and the HF component of sea level oscillations-in this way, seiches of larger bays/areas (like the abovementioned seiches recorded in Venice and Trieste) would be included in our HF signal. However, we decided to keep the 2 h limit, as it allows us to compare our Residual extremes with similar values from other studies which used 1 h sea level measurements (e.g. Pérez Gómez et al. 2021;Šepić et al. 2022a).
The middle Adriatic emerged as a "hot" location when it comes to the HF contribution to sea level extremes. Storm surges are generally weaker over this area, due to bathymetry characterized by greater depths than over the northern Adriatic and due to less significant "piling of water" effect. On the other hand, complex bathymetry of the eastern middle Adriatic provides an excellent background for the occurrence of strong HF oscillations with spectral analysis implying that most of the eastern middle Adriatic stations belong to Group 2 or Group 3 (Stari Grad, Vela Luka, Sobra, Ploče). Spectra of these stations are either polychromatic or monochromatic. Eastern middle Adriatic station Split is in Group 1, based on contribution of the HF oscillations to Residual extremes. However, as stated before, according to its spectral properties, it should be in Group 2 (Supplementary Information Fig. S5). Low contribution of the HF oscillations to Residual extremes in Split is likely because Split is well sheltered (by islands Brač, Šolta and Čiovo) from long ocean waves coming from the open sea. The Italian middle Adriatic stations (Ancona, Bari, San Benedetto del Tronto, Ortona and Vieste) also fall into either Group 2 or Group 3, and the strength of oscillation at the western (Italian) coast seems to be comparable to the strength of oscillations at the eastern coast, at least during the investigated period ( Figs. 4 and 8). This came as a certain surprise to us, as bathymetry in front of the Italian coast is not as complex as in front of the Croatian coast, and all Italian tide gauges are located within harbours which should not have, based on their shape and bathymetry, pronounced amplification factors (c.f., Rabinovich 2009 for amplification factors). Most of the middle Adriatic Italian tide gauges are located within rectangular harbours of dimensions of (500-1000 m) × (500-1000 m), which have depths of ~ 5-15 m. A rough estimate of eigenperiod of these harbours can be made, assuming a rectangular shape of harbour and flat (or linearly or parabolically sloping) bottom, following expressions from Rabinovich (2009). Resulting periods are from ~ 3 to 10 min-and such periods can indeed be found in spectra of the middle Adriatic Italian stations ( Fig. 12; Supplementary Information Fig. S5). However, at all these stations, another peak of ~ 13-27 min can be noticed. Since a peak of this period appears at all the middle Adriatic Italian stations, it is possible that this peak corresponds to edge waves (Ursell 1952) travelling along the gently sloping western Adriatic Sea. A more detailed analysis of onset times and similarity of HF oscillations at the western middle Adriatic coast should be done to assess this assumption. In addition to complex bathymetry, a reason for the increased strength of the middle Adriatic HF extremes might be that atmospheric forcing responsible for generating HF extremes over the middle Adriatic is stronger there than over other parts of the Adriatic Sea. The straightforward way to check this hypothesis would be to analyse long enough time series of 1 min air pressure measurements of spatial distribution similar to the one of our sea level measurements. It should be further noted that strong middle Adriatic HF oscillations often coincide with storm surges, i.e. with Residual extremes. Additionally, strong HF oscillation on the western side of the middle Adriatic can be a consequence of reflected long ocean waves coming from the eastern side during strong HF events, as was observed and modelled for the 1978 Vela Luka meteotsunami (Orlić et al. 2010).
Finally, the very deep southern Adriatic of not too complex coastal bathymetry is where both storm surge and HF component of sea level oscillations are rather weak. Nonetheless, the southern Adriatic station Otranto is classified in Group 3-although Otranto HF oscillations are weak, they are of comparable strength to the Otranto LF oscillations-which, at Otranto-being placed at the "mouth-of-bay" (i.e. at the entrance to the Adriatic Sea from much larger and deeper Mediterranean) are expectedly weak. Finally, we discuss Dubrovnik tide gauge-which is located at border of the southern and middle Adriatic at the entrance of a wide bay (Šepić et al. 2009b). Such tide gauge placement at the mouth of the bay results with the tide gauge recording only open ocean waves, and not seiches of the bay since the latter have a node line at the bay mouth. However, in vicinity of this tide gauge (Gruž, Komolac-both within 20 km) strong meteotsunamis have been observed (Fig. 1).
We summarize the main properties of the three discussed groups in Fig. 14 in which we present schemes of: (i) tide gauge locations; (ii) spectrums; and (iii)-(iv) histograms of residual, HF and LF series during (iii) Residual and (iv) HF extremes. Stations belonging to different groups are characterized as follows: • Group 1 stations are characterized by the lowest contributions of the HF oscillations to extreme events. These stations are commonly placed along the open coast. Their spectrum has no significant peaks and often assumes a white-noise form at the highest frequencies. Contribution of the HF signal to Residual extremes is extremely low, implying that Residual extremes are dominated by the LF signal. Even during HF extremes, the HF signal is usually not strong enough to dominate residual series. • Group 2 stations are characterized by a moderate contribution of the HF oscillations to extremes events. These stations are found either at coastlines placed along gently sloping bottoms, or within bays and harbours protected from incoming long ocean waves by nearby islands. Contribution of the HF signal to Residual extremes is moderate, with the HF oscillations dominating Residual extremes on rare occasions. During HF extremes, the HF signal can become strong enough to dominate the residual series. • Group 3 stations are characterized by the highest contribution of the HF oscillations to extremes events. These stations are placed within harbours/bays of high amplification factors, within channels prone to seiches, or in harbours/bays at coastlines placed along gently sloping bottoms. Contribution of the HF signal to Residual extremes is high, with the HF oscillations occasionally (more often than for Group 2 stations) dominating Residual extremes. During HF extremes, the HF signal is strong enough to dominate the residual series.
Finally, we did not study atmospheric background of extreme sea level events within this paper, but we can offer an insight into dominating atmospheric processes by examining 1 3 seasonal distribution of Residual and HF extremes. Residual extremes mostly occur during autumn months, with peak occurrence in November and December, i.e. during two months when cyclonic activity over the Adriatic is strongest (Pasarić and Orlić 2001). HF extremes have peak occurrence during summer months when meteotsunamigenic synoptic conditions are usually found over the Adriatic (Šepić et al. 2015c). These conditions include: (i) Fig. 14 A summary of the group properties is given. In the first row, a scheme of the shape of coastline/ harbour characteristic for tide gauge locations in three groups; spectral properties are sketched in the second row; histograms of heights that residual, HF and LF signals have during Residual extremes are sketched in the third row; sketched histograms of residual, HF and LF signals during HF extremes are in the final, forth row dry and warm air coming from Africa at altitude of around 1500 m (850 hPa); (ii) southwesterly winds having speeds of at least 20 m/s at around 5000 m (500 hPa); (iii) a dynamically unstable layer of air at heights between 600 and 400 hPa. In line with the above, HF oscillations contribute to Residual extremes most during the summer months. Further on, HF oscillations are often so strong during summer months that a HF extreme episode is often a Residual extreme as well. More precisely, we have examined a total of 528 Residual and 610 HF extremes. The number of Residual and HF extremes that coincide, i.e. happen in a span of 24 h is 122. This equals to ~ 23% of the total number of Residual and 20% of HF extremes. Focusing on summer months (i.e. May to August) when contribution of HF oscillations was the strongest, almost every Residual extreme is also a HF extreme, 27 out of 30 events in that time period to be exact. This abundance of simultaneous Residual and HF extremes during summertime, when meteotsunamis are more likely to occur because of favourable synoptic conditions (as stated before), signals that most summertime Residual extremes are due to the HF sea level oscillations. Additionally, summertime storms are more likely to produce infragravity waves that could penetrate deeper inland (Melito et al. 2018) meaning that their influence on the coast could be greater during that time. This presents an additional hazard for coastal communities, as during summer months there are more people involved in near-coast day-to-day activities. Furthermore, there are events (95 of them), when Residual and HF extremes happen simultaneously during other months. Superimposing a strong Residual extreme with any HF extreme can produce an extremely dangerous event for coastal communities. A recent example is the 2019 Venice flood ) during which HF oscillations significantly contributed to flooding.

Conclusion
Summarizing, main conclusions of our research are: 1. HF sea level oscillations can dominate positive sea level extremes. 2. Even when not dominating them, HF oscillations can considerably contribute to extreme sea levels. 3. Contribution of HF oscillations to total signal is governed by a combination of bathymetry and atmospheric forcing, resulting in the strongest HF oscillations over the middle Adriatic. 4. Residual extremes mostly happen from October to January when they are also the strongest, while HF extremes spread more evenly throughout the year, with the strongest events peaking during May to September. 5. Tide gauge stations can be divided into three distinct groups depending on the characteristics of HF oscillations which they record.
Additional studies can be performed to better investigate the HF oscillations and their contribution to either Residual or HF extreme events. These studies can incorporate research of prevailing synoptic conditions during which events occur and especially atmospheric conditions during the strongest events. Further on, acquiring time series of at least 20-year length at a significant number of stations will allow for a more precise estimation of return periods of hazardous HF oscillations. Index of likely strengths of HF oscillations and corresponding sea level height could be developed for locations of interest. This index