## 2 Materials and Methods
### Psha  
Concerning PSHA, the approach proposed by [PERSON] and [PERSON] [59; 60] was followed, utilizing the instrumental earthquake catalogue from [PERSON] et al. [25]. Given that this catalogue ends in 2009, it was deemed necessary to extend the period up to 2019 in a homogenous manner, meaning that only reviewed events from the International Seismological Centre (ISC) earthquake catalogue were taken into consideration. Therefore, we refrained from expanding the database for post-2019 events in order to retain homogeneity. It should be noted that only few \(\mathrm{M}_{\mathrm{w}}\geq 4\) magnitude earthquakes have been recorded in the study area since 2019, the largest being the \(\mathrm{M}_{\mathrm{w}}=5.3\) earthquake of 17 February 2021 [7; 47]. Thus, we considered that the annual rate of exceedance of the magnitude of completeness would not be altered significantly and would not impact the results of PSHA. According to the Cornell-McGuire approach to PSHA, the results follow a normal distribution, being essentially time-independent. Given that there are both aftershocks and foreshocks included in the utilized earthquake catalogue, a declustering procedure was conducted to retain only the mainshocks. However, it must be noted that ground motion exceedances can be caused randomly by any earthquake occurrence and a significant amount of data could be omitted by applying a declustering procedure. This would have an impact on the obtained results because annual rates of exceedance for the magnitude of completeness would be lower, resulting to the underestimation of the maximum expected ground motion parameters [58; 61]. The latter is strengthening the implementation of a non-declustered earth-quake catalogue in a Poissonian process, as ground motion at a site caused either by foreshocks or aftershocks may exceed a certain level [61]. Three seismotectonic models were considered herein; namely, the Euro-Mediterranean Seismic Hazard Model 2013 (ESHM13) [62; 63], the updated Euro-Mediterranean Seismic Hazard Model 2020 (ESHM20) [64], and the model from [PERSON] et al. [65] proposed for the Greek territory. The boundaries of the Area Source Zones are selected to include areas with similar seismological and tectonic characteristics, avoiding dividing fault systems, whenever possible. Each model divides the study area into polygons with common seismotectonic characteristics (Area Source Zones, ASZ). Individual faults were not implemented in the PSHA, as a significant part of the total seismicity of the WGoC is related to offshore faults for which certain characteristics are not known; i.e., mean slip rates. The decision to employ more than one seismotectonic model was made in order to acquire differently parameterized earthquake occurrence models so that the epistemic uncertainties would be reduced in a qualitative manner. The Modified Gutenberg-Richter (MG-R) earthquake occurrence model was applied to characterize the seismic potential of every area source zone for each seismotectonic model [66]. The MG-R was parameterized by importing the following seismicity data for each zone: (a) the b-value, (b) the threshold magnitude (M\({}_{0}\); also considered as the magnitude of completeness (M\({}_{\rm c}\))), (c) the average annual exceedance rate of M\({}_{\rm c}\) (\(\lambda\)(M\({}_{\rm c}\))), and (d) the maximum expected magnitude (M\({}_{\rm u}\)). The M\({}_{\rm c}\) and b-values were estimated for each zone using the maximum curvature (MAXC) method from [PERSON] and [PERSON] [67] and the maximum likelihood function from [PERSON] [68], respectively, with ZMAP software [69]. This method has been proved to be stable, even in cases when a small number of earthquakes within each zone are used [70; 71]. The \(\lambda\)(M\({}_{\rm c}\)) was calculated through the analysis of the earthquake catalogue and M\({}_{\rm u}\) was estimated via the Robson-Whitlock-Cooke (R-W-C) procedure [72; 73].