Probabilistic Tsunami Risk Assessment from Incomplete and Uncertain Historical Impact Records: Mediterranean and Connected Seas

Tsunami risk is considered as the probability of a particular coastline being struck by a tsunami that may cause a certain level of impact (destructiveness). The impact metric of a tsunami is expressed in terms of tsunami intensity values, K, assigned on a 12-degree scale. To calculate tsunami risk we are based on the tsunami history of the region codified in tsunami catalogues. The probabilistic model adopted was used successfully for hazard assessment of earthquakes (Kijko et al. in Bull Seismol Soc Am 79:645–654, 2016) and of tsunamis (Smit et al. in Environmetrics 30:e2566, 2019) by considering seismic magnitude and tsunami height as metrics of the respective hazards. In this model, instead of hazard metrics we inserted risk metric, i.e. wave impact in terms of intensity values. The procedure allows utilization of the entire data set consisting not only from the complete (recent) part of tsunami catalogue but also from the highly incomplete and uncertain historical part of the catalogue. Risk is assessed in terms of probabilities of exceedance and return periods of certain intensity values in specific time frames. We applied the model using catalogues for the Mediterranean and connected seas. Sensitivity analysis showed that using complete data sets generally provided more realistic results than using entire data sets. Results indicated that the risk level depends on the seismicity level and not on the size of individual ocean basin. The highest tsunami risk level was found in the eastern Mediterranean (EM), with a significantly lower risk in the western Mediterranean (WM). In the Marmara Sea (MS), the tsunami risk was low, and the lowest was in the Black Sea (BS). The risk in the small Corinth Gulf (CG, Central Greece) was comparable to that of WM. The return period of damaging tsunamis (i.e. K ≥ 7) was 22 years in the entire Mediterranean basin and 31, 118, 135, 424, and 1660 years in the EM, WM, CG, MS, and BS basins, respectively.


Introduction
Tsunami waves can have a significant effect on human life, communities, property, and the natural environment. The devastating tsunamis of 26 December 2004 in the Indian Ocean and 11 March 2011 in Japan, both generated by subduction-related earthquakes of moment magnitudes M w 9.3 and M w 9.0, respectively, clearly showed the lethal and destructive potential of mega or large tsunamis. In the Mediterranean region and its connected seas, at least 44 impactful tsunamis, with an assigned intensity equal to or larger than 6, according to the 12-grade intensity scale of Papadopoulos and Imamura (2001), have been documented from the fifth century BC up to the present (Papadopoulos et al., 2014). Most of these tsunamis were produced by earthquakes, but some have been associated with other causes, such as volcanic eruptions and coastal or submarine landslides (Maramai et al., 2014). Considerable potential for tsunami generation exists in the Mediterranean region, although the risk varies between the tsunamigenic zones (Papadopoulos & Fokaefs, 2005).
An early definition of the term ''risk,'' as related to earthquakes, derives from UNESCO initiatives (Algermissen et al., 1979;Fournier d'Albe, 1982). In this definition, risk is considered as the convolution of three main attributes, namely hazard, vulnerability, and value exposed to hazard. This definition was later generalized and applied to other risks including tsunamis (e.g., Papadopoulos & Dermentzopoulos, Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00024-023-03262-6) contains supplementary material, which is available to authorized users. 1998; Curtis & Pelinovsky, 1999). This definition of risk is consistent with definitions in more recent glossaries (EC-Working Paper, 2010;UNISDR, 2009). According to the tsunami glossary (IOC, 2019), tsunami risk is the probability of a particular coastline being struck by a tsunami multiplied by the likely destructive effects of the tsunami and by the number of potential victims. However, the mode of convolution is not standardized, as it depends on the type of risk being considered and the availability of relevant data. For example, the expected tsunami impact on the population requires treatment that is quite different from that for assets, such as buildings, infrastructure, or financial values. Consequently, the general concept of tsunami risk assessment has been used with various alternatives. For example, vulnerability could be replaced either by exposure, or exposure could be included as an additional component in the risk definition.
At the current research stage, both probabilistic and scenario-based methods for tsunami hazard assessment are subject to various aleatory and epistemic uncertainties owing to gaps in the available tsunami datasets and the applied hazard assessment techniques (e.g., Behrens et al., 2021). Tsunami hazard estimates rely heavily on historical records owing to the small amount of tsunami data contained in the brief instrumental period. Catalogues incorporating paleo-earthquakes and paleo-tsunamis, determined from geological signatures, can provide helpful information for estimating return periods of large events (e.g., Papadopoulos et al., 2014;Paris et al., 2020;Pinegina et al., 2003;Priest et al., 2017). However, the possibility of overweighting of large earlier tsunamis should be treated with caution (Geist & Parsons, 2006).
Tsunami risk assessment traditionally includes qualitative or quantitative scenario-based methods (e.g., Pararas-Carayannis, 1988;Qinghai & Adams, 1988;Sato et al., 2003;Strunz et al., 2011). Such methods have been tested in the Mediterranean region by, e.g., Papadopoulos & Dermentzopoulos, 1998;Charalambakis et al., 2007;Jelínek et al., 2012;Pagnoni et al., 2015;Pagnoni & Tinti, 2016;Triantafyllou et al., 2019). Methods for Probabilistic Tsunami Risk Assessment (PTRA) have been proposed based on various approaches (e.g., Attary et al., 2019;Dominey-Howes et al., 2010;Goda & Song, 2019;Grezio et al., 2012;Jaimes et al., 2016;Løvholt et al., 2015;Renou et al., 2011). Compared with hazard assessment data, for which only the physical parameters of the events are required, less progress has generally been noted in PTRA with respect to Probabilistic Tsunami Hazard Assessment (PTHA). Such lack of progress is ascribed to vulnerability, exposure, and value data being required for the risk assessment. However, relevant data sets are limited or even lacking for many coastal areas worldwide. Therefore, no standard methods for the PTRA assessment have been concluded yet. To support future research, various types of gaps in research and in the data required to conduct PTRA need to be considered (Behrens et al., 2021). A recent analysis of various tsunami data sources and repositories has shown the potential contribution of ontology engineering as an example of computer science methods that enable improvements in tsunami-related data management (Nacházel et al., 2021). However, currently, only a limited number of search tools can be implemented practically for ontology learning (Babič et al., 2022).
In many ways, PTHA, which has gained ground over the last few years, follows the experience obtained from the Probabilistic Seismic Hazard Assessment (PSHA) (see reviews by Geist & Parsons, 2006;Smit et al., 2017;Basili et al., 2021;Behrens et al., 2021). The PTHA methods are becoming a standard basis for the PTRA. However, PTRA methods are still less abundant and standardized compared with those employed for PTHA, as no broadly accepted approaches have been defined (Behrens et al., 2021).
In the Mediterranean region ( Fig. 1), PTRA has been attempted by some authors who considered tsunamis as random events, the size of which, expressed in terms of tsunami intensity, is distributed according to the exponential-law model, similar to that traditionally applied for the frequency-magnitude relationship of earthquakes (Gutenberg & Richter, 1944). This approach has been employed for the entire Mediterranean basin and smaller parts of the basin Papadopoulos, 2003Papadopoulos, , 2009Papadopoulos et al., 2007Papadopoulos et al., , 2010. Other approaches to PTRA in the Mediterranean region and connected seas have been published by Hancilar (2012), Grezio et al. (2012) and Yavuz et al. (2020).
The various models applied for PTHA provide highly uncertain hazard estimates owing to incomplete, heterogeneous, and doubtful tsunami datasets. From this standpoint, techniques developed for seismic hazard assessment by utilizing incomplete and uncertain earthquake catalogues are of particular interest (Kijko & Sellevoll, 1989, 1992Kijko et al., 2016). Such approaches are suitable for hazard assessment based on rare tsunami events and have been applied by, e.g., Smit et al., (2017Smit et al., ( , 2019. In these applications, the historical Oceanographic basins examined in this study. Rectangle indicates the Corinth Gulf basin Vol. 180, (2023) Probabilistic Tsunami Risk Assessment 1787 tsunami wave height, h, was adopted as a metric of tsunami hazard.
In the current study, we introduce a data-driven methodology on the PTRA. The tsunami risk metric adopted was tsunami intensity, which directly estimates the destructive effects of a tsunami. This option allows probabilisitic estimation of the expected risk in terms of destructiveness by avoiding the need to estimate hazard, vulnerability and exposure separately. The mathematical model used follows the approach introduced successfully for the PSHA from incomplete and uncertain earthquake catalogues (Kijko & Sellevoll, 1989, 1992Kijko et al., 2016) and later extended for the PTHA from incomplete and uncertain tsunami catalogues (Smit et al., , 2019. Intensity data sets have been compiled from tsunami catalogues covering both the historical and instrumental record of tsunamis in the Mediterranean and the connected seas ( Fig. 1). Our approach was initially introduced in the EC-FP7 ASTARTE tsunami research project for test sites in Greece (Papadopoulos, 2016;Triantafyllou, 2021). Here, the technique was implemented fully using tsunami data for the regional PTRA in the Mediterranean.
Although most tsunamis are produced by earthquakes, the adopted method is not restricted to seismic tsunamis, therefore, we considered PTRA regardless of the tsunami generation mechanisms, i.e., seismic or others. In the main basins of the Mediterranean region and in certain time intervals, PTRA has been calculated in terms of exceedance probabilities and tsunami return periods, with various intensity levels. The approach developed and the results obtained are of strategic interest to initiatives aimed at the mitigation of tsunami risk, such as the Tsunami Early Warning and Mitigation System in the North-eastern Atlantic, Mediterranean, and Connected Seas (NEAMTWS/IOC/UNESCO).

Tsunami Impact and Intensity
Tsunamis have the potential to wreak havoc wherever they strike, causing death and destroying various assets such as buildings, infrastructure, lifelines, vessels, materials and objects, and coastal cultivated areas. The destructive force of tsunamis can cause social disruption and financial and environmental losses. These environmental effects include the destruction of flora and fauna, ground erosion, sediment transport, and the deposition of medium-to fine-grained material, boulders, and megaclasts.
''Intensity'' is a concept often used to describe the effects of various natural processes, e.g., the 12-degree Beaufort Wind Force scale is a well-known example. For earthquakes, several macroseismic intensity scales have been in use since the end of the nineteenth century, including the 12-degree Mercalli-Sieberg scale (Sieberg, 1912), which was later transformed into the Modified Mercalli (MM) scale (Wood & Neumann, 1931), and the European Macroseismic Scale (EMS-98) (Grünthal, 1998). Sieberg (1923) introduced a 6-degree tsunami intensity scale that was later refined by Ambraseys (1962). A relevant effort in Japan (e.g., Iida et al., 1967;Imamura, 1949) concluded with the concept of the Imamura-Iida 6-degree magnitude scale based on the maximum wave height observed at the coast (see review in Papadopoulos et al., 2020). This scale was modified by Soloviev (1970), who proposed the concept of tsunami intensity, S, expressed by: where H (in m) is the mean tsunami height measured at the coast closest to the source. However, S is a metric of tsunami magnitude rather than tsunami intensity, as H is still a physical feature and does not describe impact at all. The Soloviev scale is, therefore, preferred for tsunami hazard assessment (e.g., Smit et al., 2017) but not for tsunami risk assessment. Following the longstanding seismological tradition, Papadopoulos and Imamura (2001) introduced a 12-degree tsunami intensity scale that describes the various levels of tsunami destructive effects in the built and natural environments (e.g., damage in buildings or in vessels, people washed away etc.) (Table S1). From the description of tsunami effects a numerical (integer) intensity value is assigned scaling from 1 to 12. Tsunami intensity is not a scaled height metric. The intensity of a tsunami event in a particular inundation place is assessed only on the basis of the wave impact, it is not determined from the wave height at all. Therefore, tsunami intensity, K, is a metric of damage, i.e. a metric of risk and not a metric of hazard. On the other hand, there is a rough correlation of K with wave height. This is exactly the long tradition gained from the use of macroseismic intensity scales in seismology since the late of nineteenth century. In seismology, intensity roughly scales with ground acceleration but is not determined from the ground acceleration.

PTRA from Incomplete and Uncertain Impact Records
In most instances, instrumental records of earthquake occurrences are more complete, homogeneous, and accurate than historical earthquake data. Unfortunately, the instrumental records are relatively short, covering, at best, a period of * 120 years. However, the return periods of strong earthquakes in specific fault segments or zones are significantly longer than 120 years. Accordingly, reliable seismic hazard assessment requires incorporating data from earthquake catalogues preceding the instrumental era. However, including historical events in the analysis increases the incompleteness and heterogeneity of the catalogue and the errors involved in earthquake parameters, such as location, depth, and magnitude. Moreover, the dates of pre-instrumentally recorded earthquakes are often not well known.
In view of these challenging issues, a maximum likelihood approach has been introduced that utilizes incomplete (extreme) and complete catalogue segments with different magnitude thresholds (Kijko & Sellevoll, 1989) (Fig. 2). This approach also incorporates the calculations of magnitude heterogeneity and an estimation of the maximum magnitude expected in a fault or fault zone (Kijko & Sellevoll, 1992), as well as the uncertainties involved in the earthquake occurrence model (Kijko et al., 2016). This method appears quite promising for assessing PTHA (Smit et al., , 2019, as tsunami catalogues comprise highly incomplete, heterogeneous, and uncertain data files, which, as a rule, contain a relatively small number of events. Regardless the quality of the data, prehistoric and historical information can be combined with the more complete data files of the instrumental period (Fig. 2). In PTHA, the approach focuses on the empirical assessment of three main tsunami recurrence parameters, namely (1) the mean tsunami activity or occurrence rate, k; (2) the b-value, i.e., the slope of the log-linear domain of the frequency-size relationship; and (3) the coastline-characteristic maximum expected tsunami wave height, h max . In this approach, the wave height is a metric of the tsunami hazard. Replacement of tsunami height with an appropriate metric of destructiveness, such as tsunami intensity, K, transforms the PTHA method into a PTRA method. The mean tsunami activity rate, k, describes the number of tsunamis expected to occur within a specified stretch of coastline with damaging effects corresponding to a threshold value equal to or larger than K min within a unit of time, e.g., one year.
Among the various intensity scales available, the 12-degree scale proposed by Papadopoulos and Imamura (2001) has been implemented in the current study, as this scale provides a superior resolution for a more detailed description of the intensity levels. Based on worldwide data (Shuto, 1993), this scale also provides an empirical relationship of the intensity with the wave height, h. Such empirical relationships provide a rough idea on how K scales with h. Figure 3 shows the scaling of K with h for various intensity measures worldwide, as well as in the Mediterranean region. Compared with the other intensity measures considered, a superior K/h correlation was found for the EM basin (formula 2). The relation and correlation coefficient, r, are given by:

Temporal Distribution
The temporal distribution of tsunamis is assumed to follow a random (Poisson) model (e.g., Geist & Parsons, 2006), which is the most simplified and popular model adopted for earthquake occurrence (e.g., Lomnitz, 1973). For the Poisson model, the probability that n tsunami events will be observed along a certain stretch of coastline within a specified time interval, t, is where k : k(J min ) is the mean (annual) activity rate of tsunamis, with an intensity greater than or equal to J min .

Size-Frequency Relationship
In the tsunami size distribution model, parameter b plays the same role as parameter b in the frequency-magnitude relationship for earthquakes (Gutenberg & Richter, 1944): N is either the discrete frequency of magnitudes, M, in each magnitude bin or the cumulative frequency of magnitudes C M; a and b are parameters determined by the earthquake catalogue data. Parameter b is an expression of the relative number of small magnitude earthquakes to the large magnitude ones, whereas a is a measure of the seismicity level. Considering tsunami intensity, K, instead of earthquake magnitude, M, and setting b = b ln(10) and A = ln a, relation (4) is transformed to From (5), the tsunami size expressed in terms of intensity, K, scales according to an exponential law. The Hard Bounds Model (Kijko & Sellevoll, 1992) was followed here to account for the uncertainty involved in tsunami intensity determination, implying that the real, unknown tsunami intensity is assumed to occur within fixed boundary limits. In the current analysis, these limits were determined independently by an empirical expert judgment for each tsunami event.

Figure 2
Generalized illustration of data used to assess model recurrence parameters based on prehistoric (paleo-), historical, and instrumental datasets. Horizontal dimension expresses time. Vertical dimension expresses event size (e.g., magnitude or intensity). Vertical length and horizontal width of black rectangles represent uncertainty in size and time of occurrence of the observed wave, respectively. Prehistoric data are subject to uncertainty relevant to the time of occurrence, exact size of the event, and incompleteness in terms of the probability of detecting an event. Due to uncertainties in the time of occurrence of prehistorical data these times are represented by gray bars. Historical data, consisting of the largest observed events, as well as instrumental datasets are subject to incompleteness and uncertainty relevant to the observed event size, varying levels of certainty regarding the exact location of an event, and varying probabilities of all events above a certain minimum size being observed. White rectangles in the instrumental group represent time intervals within which data are available, while Time gaps (Tgs)

Maximum Likelihood Estimation
The coastline-characteristic tsunami risk parameters h = (k,b,J max ) are estimated using the maximum likelihood estimation (MLE) technique. Likelihood functions are defined for both the historical and the instrumental parts of the catalogue, LH(h) and LI(h), respectively, with LI(h) consisting of the likelihood function of each of the s sub-catalogues. Utilizing the additive property of the log-likelihood function L(h) (Rao, 1973), the MLE function for the entire catalogue takes the form Aleatory and epistemic uncertainties are introduced in the tsunami models using mixing distributions. For example, a mixing Poisson-gamma distribution based purely on empirical data was chosen for the uncertainty in the occurrence model  as an alternative formalism to the classic Bayesian method. The uncertainties involved in the recurrence parameters are estimated through formalism called Weighted Maximum Likelihood (Schoenberg, 1997).
In PTRA, assuming that the tsunami intensities observed along a stretch of coastline are continuous and independent random values distributed according to an exponential law, the probability distribution function (PDF), f, (K) and the cumulative distribution function (CDF), F(K), of tsunami intensity take the same functional forms (e.g., Burroughs & Tebbens, 2005) as those for earthquake magnitude (e.g., Cosentino et al., 1977;Kagan, 2002). Consequently, both PDF and CDF are defined for a minimum intensity threshold value, K min , and are truncated from the top by the coastline-characteristic maximum possible tsunami intensity, K max , so that: In both (7) and (8), for K \ K min we have f(K) = 0 and F(K) = 0; for K [ K max , f(K) = 0 and F(K) = 1. A similar formalism is applied to the distribution of the tsunami wave height, h, if PTHA is considered.
The most difficult parameter to assess is the coastline-characteristic maximum possible tsunami intensity, K max . In PTHA, parameter h max was assessed  in a way similar to that of the regional characteristic maximum possible earthquake magnitude, M max , in PSHA (Kijko & Sellevoll, 1989, 1992. As K max cannot be calculated using the sample likelihood function L(h), the procedure suggested by Cooke (1979) was followed, where the possible maximum intensity is determined as the maximum observed intensity in the entire catalogue, plus various additional terms (for more details and mathematical formulation see Kijko, 2004;Kijko & Singh, 2011;Smit et al., 2017).
The three coastline-characteristic tsunami risk parameters (k, b, K max ) were estimated by maximizing the likelihood function expressed by Eq. (6). The likelihood functions were derived separately for the historical and instrumental parts of the entire tsunami catalogue. Details on the mathematical formulation can be found in Smit et al. (2017).

Basins Examined
In this study, we conducted PTRA in the Mediterranean Sea region and connected seas. We discriminated between three large and three small basins in this geographic region (Fig. 1). The entire Mediterranean Sea, with a surface area of approximately 2,500,000 km 2 , is the largest basin and is separated into two main basins, namely the eastern (EM) and the western (WM) basins. However, the complex physiography of the Mediterranean Sea complicates the allocation of tsunami sources between the eastern and western basins. We assumed that the area to the west of the Tunisia-Sicily passage is part of the WM basin and includes the Tyrrhenian Sea. We considered the Adriatic Sea as belonging geographically to the EM basin, which, in this Vol. 180, (2023) Probabilistic Tsunami Risk Assessment allocation, has a surface area equivalent to * 3/5 of the entire Mediterranean basin surface area. Tsunamis produced in the Strait of Messina, at the southern tip of the Italian peninsula, were included in the geographic area of the EM basin. This decision was based on the good propagation of these tsunamis toward the south Ionian Sea, such as the strong tsunami of 28 December 1908. Careful consideration was given to the distinction of tsunamis that occurred in the Tyrrhenian Calabria or Ionian Calabria, Italy, for allocation either to the WM or to the EM basins, respectively.
The Black Sea, including the Azov Sea, the Marmara Sea, and the Gulf of Corinth, are the three smaller basins, with surface areas of 436,402 km 2 , 11,350 km 2 , and 2500 km 2 , respectively (Fig. 1). The Gulf of Corinth, Central Greece, is a nearly closed basin of small dimensions and, therefore, constitutes a special case. This basin has an elongated shape, stretching from east to west, with a length of * 115 km and maximum width of * 25 km (Fig. 1). At its eastern end, the Gulf of Corinth is practically closed and its only connection with the Aegean Sea is the narrow Isthmus of Corinth, with a length of 6.3 km. At its western end, the Gulf of Corinth has a narrow opening to the Ionian Sea through the 2.7 km long Rio-Antirrion Strait. The rate of tsunami occurrence is significant in the Gulf of Corinth because of the high level of seismicity in the area and its susceptibility to coastal and submarine landslides that also produce tsunamis (Papadopoulos, 2003). Although powerful and destructive tsunamis have been reported from both the northern and southern coasts of the Gulf of Corinth, even the strongest tsunamis cannot propagate with significant amplitudes outside the basin (e.g., Stefatos et al., 2005).

Catalogue of Tsunami Intensities
Various tsunami catalogues contain tsunami impact data deriving from numerous studies, reports, and books. Non-parametric or semi-parametric tsunami catalogues for the Mediterranean Sea and its connected seas have been published systematically since the 1960s (e.g., Ambraseys, 1962;Antonopoulos, 1979;Galanopoulos, 1960). The cataloguing continued in the 1980s, with the introduction of parametric catalogues (e.g., Papadopoulos & Chalkis, 1984). The cataloguing of tsunamis followed the standard format proposed by EC-supported research projects, leading to the establishment of a standardized European Tsunami Catalogue from the beginning of the 1990s (e.g., Maramai et al., 2014;Papadopoulos, 2001Papadopoulos, , 2003Papadopoulos et al., 2007Papadopoulos et al., , 2010Papadopoulos et al., , 2011Tinti & Maramai, 1996;Tinti et al., 2004). These catalogues list information on tsunami parameters, parameters of the causal sources (e.g., earthquakes), and tsunami effects on the anthropogenic and natural environments.
We consulted numerous descriptive and/or parametric catalogues and other relevant information with the aim to compile a catalogue of tsunami impact and to assign tsunami intensities for the largest possible number of tsunami events. In addition to the catalogues and studies mentioned before, we also consulted those by Guidoboni et al. (1994), Yalçıner et al. (2002), Papadopoulos (2009)  , and .
In total, 256 instrumental and pre-instrumental tsunami events were included in our PTRA analysis. Table 1 shows the allocation of these events to the six basins examined. Table S2 shows the full chronological list of tsunamis and their respective intensities. The data set we compiled covers the period from prehistoric times up to 2021. The quasi-exponential  Table S2 shows the catalogue list of 256 tsunami events, each with an estimated maximum tsunami intensity, K max , and respective standard deviation, which is a measure of the uncertainty related to K max . As the Hard Bounds Model was adopted for the uncertainty related to tsunami intensity, we introduced values of the form K max ± j, where j = 0, 0.5, 1, 2. In a few instances, incomplete event dates are available for pre-instrumental tsunamis, e.g., only the year of occurrence of the tsunami is known. In such instances, 1 July of the given year is considered conventionally as the event date. The 15th of a given month was adopted when only the day is missing from a certain date.
The instrumental period of the tsunami record conventionally starts in 1900, following the start of the instrumental period of seismology, implying determination of more accurate earthquake focal parameters. The instrumental period of tsunami recording in the Mediterranean region practically started only after the establishment of the NEAMTWS/IOC/UNESCO, by the end of 2005. The creation of NEAMTWS accelerated the deployment of modern, digital tide-gauge networks in the Mediterranean region. However, analog instruments, not suitable for accurately recording tsunamis, were operated for several years before 2005. Nevertheless, the macroscopic observation of tsunamis gradually improved after 1900.

Remarks on Data
One of the most important tsunamis listed in the catalogues consulted is the tsunami associated with . Therefore, this tsunami event was assigned a reliability of only 2 in contrast with the higher reliability assigned by previous authors. Further, we carefully examined the catastrophic earthquake and tsunami that reportedly occurred in the northern Evoikos Gulf, Central Greece, in AD 552. The event was reported by the contemporary Byzantine historian, Procopius, which is the only literary source. However, thorough examination indicated that the report of a tsunami was false ; therefore, this event was not included in the list of tsunamis we considered. The date of the AD mid-eighth century earthquake and tsunami that occurred in the Levantine Sea has been an issue of dispute (e.g., Ambraseys, 2009;Guidoboni et al., 1994;Salamon et al., 2007). Following Salamon et al. (2007), we adopted the date of 18 January 746 ).

Sensitivity Analysis
Initial calculations of the probabilities of exceedance, P, and return periods, R, of certain tsunami intensity values for each of the six basins examined were performed using the full data set. However, we also conducted a sensitivity analysis procedure for each basin by removing the less complete segment of data from the full data set and recalculating P and R using only the more complete part of the data set. Recalculation was performed for the three large basins, i.e., the entire Mediterranean, EM, and WM, as well as for the Black Sea basin using only the instrumental data set after removing the pre-instrumental events from the respective tsunami catalogues (Tables 1, S2). As regards the Marmara Sea, only four events are included in the instrumental part of the catalogue; therefore, sensitivity analysis was conducted using the more complete catalogue segment covering the period from AD 1265 onward (Table S2). As only six events are listed in the instrumental period of the catalogue of the Gulf of Corinth, for sensitivity analysis we used the more complete segment of the catalogue, i.e., AD 1742 onward (Table S2).

Results
Our results are summarized in Tables 2 and 3 for the return periods and exceedance probabilities, respectively. An essential common feature characterizes the results we obtained for the various basins, namely the level of the estimated tsunami risk is driven mainly by the level of seismic activity within the basin rather than the size of the basin. This finding is independent of the data sets employed, i.e., the full data sets or the complete data sets. Although tsunamis are not only generated by earthquakes but also by other causes, e.g., landslides and volcanic eruptions, seismic activity was the main cause of tsunami generation in all the basins examined. Compared with the entire Mediterranean Sea basin, the EM basin provides clear evidence of that feature ( Table 2, Fig. 6). In the two basins, the return period, R, for the same level of intensity, K, is quite close, which is ascribed to the seismic activity in the EM being much higher than in the rest of the Mediterranean. Therefore, the results obtained for the entire Mediterranean basin are driven mainly by the data sets for the EM basin. This explanation is relevant also for the significantly longer return periods of intensities in the WM basin compared with those in the EM. The return periods calculated in all the basins from the entire data sets and for several levels of K are systematically longer than those calculated from the complete data sets (Table 2, Fig. 6). This discrepancy is caused by the nature of the applied data when historical and prehistoric data are present in addition to complete data. As these data sets are highly incomplete and inaccurate compared with the complete data sets, the return periods calculated from the complete data sets are more realistic. For instance, in the entire Mediterranean Sea basin, R of 89, 189, and 1 250 years were calculated for intensities 7, 8, and 10, respectively, using the full data set (Table 2). However, the respective R values obtained using only the complete (instrumental) data sets are much shorter, i.e., 22, 38, and 180 years for the respective K levels. These values are more realistic, as they are close to the values of R = n/T for the Table 2 Mean return period, R, of tsunami intensity, K, in the basins examined h is the respective wave height from the empirical formula (2). The activity rate, k (events/year), is k = 1/R. J max ± r max is the expected maximum intensity along with the respective estimated standard deviation r max . Tsunami b is the estimated value of parameter b inserted in formula (4). In each basin the upper and lower row of return period, R, has been calculated from the full and complete data sets, respectively Vol. 180, (2023) Probabilistic Tsunami Risk Assessment 1795 complete data sets, i.e., 30, 40, and 122 years, where T = 122 years is the time length of the instrumental catalogue and n is the number of tsunami events with intensity 7 (n = 4), 8 (n = 3), and 10 (n = 1). Details on our results are presented in the next lines.

Mediterranean Basins
In the entire Mediterranean Sea basin, as well as in the EM basin, the maximum observed intensity and maximum expected intensity are K mobs = 10 and K max = 11, respectively (Table 2). However, the respective values in the WM basin are significantly lower, namely K mobs = 7 and K max = 8. In the entire Mediterranean Sea, the calculated return periods, R, of the various tsunami intensity levels from the complete data set are 14 and 22 years for slightly damaging (K = 6) and damaging (K = 7) tsunamis, respectively (Table 2, Fig. 7). The R values for heavily damaging (K = 8), destructive (K = 9), and very destructive (K = 10) tsunamis increase to 38, 72, and 180 years, respectively. However, as regards the full data set that covers both the instrumental and the pre-instrumental periods, only five tsunamis of estimated intensity 10 ± 1 are known, all having occurred in the EM basin (Table 2). One event occurred during the 122 years of the complete data period , implying that one or more Table 3 Probabilities of exceedance of tsunami intensity level, K, calculated for the most complete time segments of the tsunami catalogues (see text for sensitivity analysis) in the basins examined in time intervals T (years) powerful tsunamis that likely occurred in the preinstrumental period had not been recorded. The calculated R for devastating tsunamis (K = 11) is much longer, namely 2 150 years, but the associated uncertainty is significant (Fig. 7). No tsunami of such extremely high intensity is included in our database.
In the EM basin, the intensity return periods are longer but comparable to those calculated for the entire Mediterranean Sea (Table 2, Fig. 7). This result is ascribed to most tsunami events occurring in the EM basin, including the larger events contained in the data set for the entire Mediterranean. In the WM basin, however, the return periods are significantly longer, namely 30, 55, 135, and 1 600 years for K = 5, 6, 7, and 8, respectively (Table 2, Fig. 8 left).

Black Sea Basin
The tsunami risk is significantly lower in this basin. The values of K mobs = 6 and K max = 7 are considerably lower than are those in the Mediterranean Sea. For tsunami intensities 5, 6, and 7, the respective return periods in this basin are long, namely 46, 128, and 1 660 years (Table 2, Fig. 8 right).

Marmara Sea Basin
In this basin, the level of tsunami risk is the lowest of all the considered areas for K B 6 (Table 2, Fig. 9 left). The calculated R values are 115, 164, and 249 years for intensities 4, 5, and 6, respectively, i.e., longer than the R values in the other basins. However, R = 424 years is indicated for K = 7, which is much shorter than that in the Black Sea, whereas R for K = 8 is much shorter (955 years) than that in the WM basin. These results are ascribed to the maximum observed intensity K mobs = 8 and the maximum estimated intensity K max = 9 in the Marmara Sea basin, which exceed the respective values in the WM and Black Sea basins by one and two intensity degrees, respectively.

Gulf of Corinth
A significant level of tsunami risk characterizes the Gulf of Corinth, although it is a small basin, with a surface area, S, of only 2500 km 2 . The maximum intensity assigned is as high as K mobs = 9, whereas the maximum estimated intensity is K max = 10. In addition, the calculated R for several K levels is much shorter than that found in the much larger basins of the Black Sea (S = 436,402 km 2 ) and the Marmara Sea (S = 11,350 km 2 ) ( Table 2) (Figs. 8,9). Tsunami intensity return periods in the Corinth Gulf are comparable to those in the WM basin for K B 7, but much shorter for K = 8. The high tsunami activity in the Corinth Gulf is ascribed to both a high level of seismicity and the landslide susceptibility along the shores of the area (Papadopoulos, 2003).

Probabilities of Exceedance
The probabilities of exceedance, P, of certain tsunami intensity levels in the basins examined are illustrated in Figs. 10, 11, 12 and summarized in Table 3. Figure 13 shows a comparison of P calculated for various basins. In the entire Mediterranean basin, high probabilities of exceedance have found for slightly damaging (K = 6, P = 0.96) and damaging (K = 7, P = 0.87) tsunamis in a time frame of 50 years (Fig. 10). The exceedance probabilities of heavily damaging (K = 8), destructive (K = 9), and very destructive (K = 10) tsunamis within the time frame of 100 years are 0.91, 0.73, and 0.42, respectively. The respective probabilities in the EM basin are lower but comparable (Fig. 10). However, in the rest of the basins, the probabilities of exceedance of certain tsunami intensity values are significantly lower than those in the EM basin (Figs. 11,12). As expected, the probabilities are relatively high in the Gulf of Corinth basin despite its small size, e.g., P = 0.56 for exceeding K = 6 and K = 7 in time frames of 50 and 100 years, respectively.

Discussion
The tsunami risk assessment performed with the described methodology considered several oceanographic basins in the Mediterranean and connected seas and is, therefore, of interest in the context of regional scale. However, it is important for developing strategic plans that require the prioritization of actions aimed at long-term tsunami risk mitigation, including further development of monitoring networks and relevant infrastructures that can support early warning systems, such as the NEAMTWS/IOC/ UNESCO. More local or site-specific tsunami risk assessments employing the proposed methodology are also possible, provided that tsunami impact data are available. It is our intention to present relevant results in a future study.
Our results are directly comparable to only a few relevant studies published for PTRA in the oceanographic basins examined. For example, the approach to PTRA for the Corinth Gulf and the entire Mediterranean basin (Papadopoulos, 2003(Papadopoulos, , 2009) was a model combining the random distribution of tsunami events in time with the tsunami frequencyintensity distribution following exponential-law. However, the 6-degree tsunami intensity scale was considered in those studies. Therefore, for a direct comparison of results, we converted intensities, k, on the 6-degree scale to intensities, K, on the 12-degree scale, according to the empirical conversion relationship proposed for the entire Mediterranean basin by Papadopoulos and Fokaefs (2005). Namely, k of 3, 4, 5, and 6 were converted to K of 5, 7, 9, and 11, respectively. After these conversions, we found that the tsunami risk in terms of the mean return period, R, and the probabilities of exceedance of certain tsunami intensity levels, K, in the entire Uncertainty intervals were computed as noted in Fig. 7 Vol. 180, (2023) Probabilistic Tsunami Risk Assessment Mediterranean and the Gulf of Corinth is generally overestimated by the simple random/exponential-law model adopted by Papadopoulos (2003Papadopoulos ( , 2009 in comparison with the results obtained by the methodology applied in the present study. For example, in the entire Mediterranean basin, the previous authors found R equal to 4, 13, 41, and 132 years for K C j, with j = 5, 7, 9, and 11, respectively, whereas we calculated R equal to 8.5, 22, 72, and 2 150 years. An indirect comparison of our results can be conducted with PTHA studies performed in the Mediterranean basin. For instance, Sørensen et al. (2012) found that for a wave height h = 10 m somewhere in the entire Mediterranean basin, the probability of exceedance, P, is approximately 0.3, 0.4, and 0.83 in time frames, T, of 1, 30, and 100 years, respectively. These results are comparable to ours (Table 3), which indicated P equal to 0.3, 0.72, and 0.91 for the exceedance of tsunami intensity 8, with the corresponding h = 9.2, in T of 1, 50, and 100 years, respectively. However, the results of Sørensen et al. (2012) are not consistent with our P for low wave heights in short time frames. For Figure 10 Probabilities of exceedance of certain tsunami intensity levels in the entire Mediterranean and the Eastern Mediterranean basins, calculated from the complete data set (C) in time frames of 1 year (left) and 10, 50, and 100 years (right). Uncertainty intervals were computed as noted in Fig. 7 example, for h = 1 m and T = 1 year, they found P = 0.73, whereas ours is only P = 0.17 for the exceedance of intensity 4, corresponding to h = 1.2 m. For h = 2 m and T = 1 year, their result is P = 0.37, and ours is only P = 0.11 for the exceedance of intensity 5, which corresponds to h = 2.3 m. This inconsistency implies that the Monte Carlo procedure followed by Sørensen et al. (2012) for producing a synthetic tsunami catalogue, combined with their tsunami simulations, tended to overestimate the probability of exceedance of relatively low or moderate tsunami heights in short time frames. However, our result, i.e., tsunami risk is highest in the EM basin among the rest basins examined, is consistent with their conclusion that the tsunami hazard in this basin is undoubtedly higher than that in the WM basin.
In the entire NEAM Region, Basili et al. (2021) found that maximum inundation heights (MIHs) of several meters are rare but not impossible. Considering a 2% probability of exceedance in 50 years, corresponding to R&2 475 years, the points of interest with MIH [ 5 m are fewer than 1% and are all in the Mediterranean and particularly on the coasts of Libya, Egypt, Cyprus, and Greece, i.e. within the Eastern Mediterranean basin.
Since our analysis covers a long span of time, from the prehistoric to instrumental eras, one point of concern related to our method could be how the tsunami intensity metric (K) varies with historical time for equivalent tsunamis given that the descriptors of the tsunami intensity scale include observations from shore, felt reports on vessels, and various damage indicators. Historical changes in population density at the coast, vessel and building construction, coastal development, etc. should affect the impact for a given tsunami event (i.e., fixed source location, parameters). For this reason, a sensitivity test was performed, i.e. calculations of the probabilities of exceedance, P, and return periods, R, of certain tsunami intensity values for each basin have been made not only from the respective full data set but also from only the respective complete data set.
The PTRA methodology proposed in this study has the advantage of being based solely on observational impact data in terms of tsunami intensity, and, at the same time, incorporating data incompleteness, heterogeneity, and uncertainty. Intensity is an appropriate risk metric, as it integrates the overall destructive effects of tsunami events. Moreover, the methodology accepts the aggregation of data for tsunamis generated by earthquakes as well as by other causes, such as landslides, volcanic eruptions, and more. The disadvantage is that no data are available in many specific coastal areas. However, an alternative could be the organization of synthetic tsunami impact catalogues.

Conclusions
We successfully developed methodology for Probabilistic Tsunami Risk Assessment (PTRA) in the oceanographic basins of the entire Mediterranean, as well as of the Eastern Mediterranean (EM), Western Mediterranean (WM), Black Sea (BS), Marmara (MS) seas, and of the Gulf of Corinth in Central Greece. Tsunami intensity, K, according to the 12-grade Papadopoulos and Imamura (2001) scale was used as a risk metric employing tsunami impact data collected from instrumental and pre-instrumental observation periods regardless of the causes of the tsunami events, i.e., seismic or other. The methodology incorporates the incompleteness and uncertainties involved in the data files. The PTRA was approached in terms of return periods, R, and probabilities of exceedance, P, of certain K values in certain time frames, T. After sensitivity analysis, for each basin examined more realistic results were Figure 13 Comparison of the probabilities of exceedance, P, of certain tsunami intensity, K, in the large (black symbols) and small basins in time frames of 1, 10, 50, and 100 years, calculated using the respective complete data sets. Basin notation: Entire Mediterranean (M), Eastern Mediterranean ( obtained by using the respective complete data file compared with results based on the full data sets incorporating both instrumental and pre-instrumental data. The tsunami risk level in the various basins examined was summarized based on the more realistic results obtained by employing the complete data files.
In the entire Mediterranean basin, the calculated P values were 0.04, 0.34, 0.87, and 0.98 for exceedance of intensity K C 7, i.e., a damaging tsunami, in T of 1, 10, 50, and 100 years, respectively. The P values found for the EM basin were quite close, namely 0.03, 0.23, 0.78, and 0.94, ascribed to most strong submarine earthquakes capable of tsunami generation occurring in the EM. For K C 7, the R values in the entire area and the EM were 22 and 31 years, respectively.
In the WM basin, the tsunami risk was significantly lower compared with that in the EM. The P values found were 0.007, 0.07, 0.31, and 0.52 for K C 7 in T of 1, 10, 50, and 100 years, respectively. The return period of K C 7 was as long as 135 years.
Among the three smaller basins, the Gulf of Corinth, which is by far the smallest one, had the highest risk, as probability, P, was 0.009, 0.08, 0.34, and 0.56 for K C 7 in T of 1, 10, 50, and 100 years, respectively. The return period of K C 7 was as long as 118 years. These findings indicated that the tsunami risk in this small area was comparable to that in the much larger WM basin. The result is ascribed to the Gulf of Corinth being characterized by a high rate of tsunamigenic earthquakes and high susceptibility to coastal and submarine landslides, which also often cause tsunamis (Papadopoulos, 2003).
The tsunami risk was low in the MS, given that for K C 7, the value of R = 424 years was large enough, whereas P was as low as 0.002, 0.02, 0.11, and 0.21 for K C 7 in T of 1, 10, 50, and 100 years, respectively. In the BS, the lowest tsunami risk indicator, P, was 0.0006, 0.006, 0.03, and 0.06 for K C 7 in T of 1, 10, 50, and 100 years, respectively, while the return period of K C 7 is extremely long, i.e., R = 1 660 years. Accordingly, the tsunami risk level does not depend on the size of the basin but on the level of seismicity, which is the main cause of tsunamis.

Acknowledgements
This research is a contribution to the European Cooperation in Science and Technology COST project ''AGITHAR-Accelerating Global Science in Tsunami Hazard and Risk Analysis''. We thank Dr. Ruan Vermeulen for his constructive comments on the initial version of the manuscript. Thanks are also due to three anonymous reviewers for their suggestions that helped to improve the paper. Figure 1 was created with the use of Arcmap 10.8. Figures from 7 to 12 were created with the use of Matlab graphics.

Funding
Open access funding provided by HEAL-Link Greece. No funding has been provided for this research.

Data availability
The data sets used in this paper and the relevant citations are explained in the main text and summarized in Supplementary file 1.

Declarations
Conflict of interest The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.