*Article* **How Useful Are Strain Rates for Estimating the Long-Term Spatial Distribution of Earthquakes?**

**Sepideh J. Rastin \*, David A. Rhoades, Christopher Rollins and Matthew C. Gerstenberger**

GNS Science, 1 Fairway Drive, Avalon, P.O. Box 30-368, Lower Hutt 5040, New Zealand; d.rhoades@gns.cri.nz (D.A.R.); c.rollins@gns.cri.nz (C.R.); m.gerstenberger@gns.cri.nz (M.C.G.)

**\*** Correspondence: s.rastin@gns.cri.nz; Tel.: +64-21-2810336

**Abstract:** Strain rates have been included in multiplicative hybrid modelling of the long-term spatial distribution of earthquakes in New Zealand (NZ) since 2017. Previous modelling has shown a strain rate model to be the most informative input to explain earthquake locations over a fitting period from 1987 to 2006 and a testing period from 2012 to 2015. In the present study, three different shear strain rate models have been included separately as covariates in NZ multiplicative hybrid models, along with other covariates based on known fault locations, their associated slip rates, and proximity to the plate interface. Although the strain rate models differ in their details, there are similarities in their contributions to the performance of hybrid models in terms of information gain per earthquake (IGPE). The inclusion of each strain rate model improves the performance of hybrid models during the previously adopted fitting and testing periods. However, the hybrid models, including strain rates, perform poorly in a reverse testing period from 1951 to 1986. Molchan error diagrams show that the correlations of the strain rate models with earthquake locations are lower over the reverse testing period than from 1987 onwards. Smoothed scatter plots of the strain rate covariates associated with target earthquakes versus time confirm the relatively low correlations before 1987. Moreover, these analyses show that other covariates of the multiplicative models, such as proximity to the plate interface and proximity to mapped faults, were better correlated with earthquake locations prior to 1987. These results suggest that strain rate models based on only a few decades of available geodetic data from a limited network of GNSS stations may not be good indicators of where earthquakes occur over a long time frame.

**Keywords:** earthquake forecasting; strain rates; multiplicative hybrid models; reverse testing; spatial distribution

#### **1. Introduction**

Estimation of the long-term spatial distribution of earthquakes is an essential component of Probabilistic Seismic Hazard Analysis (PSHA). In the traditional approach to PSHA, the estimation is based mainly on the combined information from studies of the long-term slip rates and rupture histories of mapped faults and spatial smoothing of the locations of past earthquakes. However, other emerging data streams, including geodetic observations from the global navigation satellite system (GNSS) and derived strain-rate models, seem to offer promise for improving the estimation [1–5]. In a revision of the New Zealand (NZ) National Seismic Hazard Model (NZNSHM) currently being undertaken, we are interested in assessing the relative value of different available data inputs, including information from fault studies, tectonics, the earthquake catalogue, and strain rate models, for estimating the long-term spatial distribution of earthquakes. In particular, we are interested in assessing how strongly strain rate estimates based on a few decades of geodetic data should be weighted when estimating the long-term spatial distribution of earthquakes.

A multiplicative hybrid modelling framework [6] was developed to optimally combine a set of spatially gridded forecasts submitted to the five-year regional earthquake likelihood

**Citation:** Rastin, S.J.; Rhoades, D.A.; Rollins, C.; Gerstenberger, M.C. How Useful Are Strain Rates for Estimating the Long-Term Spatial Distribution of Earthquakes?*Appl. Sci.* **2022**, *12*, 6804. https://doi.org/10.3390/ app12136804

Academic Editors: Eleftheria E. Papadimitriou and Alexey Zavyalov

Received: 25 May 2022 Accepted: 2 July 2022 Published: 5 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

models (RELM) experiment in California [7,8]. It was found that, although one best model strongly outperformed all others individually [9], it was possible to find several hybrid models that could outperform the best individual model over the 5-year period of the experiment. Such hybrid models were formed from pairs of models, including the best model as the baseline.

Multiplicative hybrid modelling has since been used to combine gridded information from fault, earthquake, and strain rate covariates in NZ [10,11] with a simple Stationary Uniform Poisson (SUP) baseline earthquake likelihood model. When hybrids only of fault- and earthquake-based variables were considered, any hybrid model that included one particular covariate called "Proximity to Mapped Faults" (PMF) outperformed all hybrids without PMF in both a 20-year fitting period and an independent testing period [10]. The PMF covariate was derived from the locations of NZ mapped faults and their estimated slip rates. Subsequently, Rhoades et al. [11] used models for maximum shear, rotational and dilatational strain rates by Beavan [12], based on GNSS data from 1996–2011, as additional covariates. They found that any hybrid including shear strain rate (BSS) outperformed all hybrids not including BSS in both the 20-year fitting period and an independent testing period.

In the present study, we restrict our attention to (maximum) shear strain rate models. In addition to BSS, we consider the global shear strain rate model of Kreemer et al. [13] based on GNSS data over the period 1996–2013 (GSS) and an updated NZ strain rate model by Haines and Wallace [14] based on GNSS data from 1995–2013 (HWS). We analyse the performance of multiplicative hybrids during the fitting and testing periods used by Rhoades et al. [11] with GSS and HWS as alternative strain rate covariates. We then examine the stability of the multiplicative modelling results by testing them against earthquakes from 1951 to 1986. Furthermore, we evaluate the strain rate and other covariates as earthquake predictors using Molchan error diagrams. We also produce smoothed scatter plots of covariate values corresponding to targeted earthquakes against time from 1951 to 2020. These results are relevant to the question of how useful strain rate estimates are for the assessment of long-term seismic hazards.

#### **2. Method and Data**

#### *2.1. Multiplicative Hybrid Model*

We follow the method used by Rhoades et al. in [6,10,11] to produce hybrid earthquake likelihood models from a baseline model and a set of covariates. As in [10,11], the stationary uniform Poisson (SUP) model is used as the baseline model, and subsets of the available spatially varying covariates are assimilated into it to form a range of hybrid models.

The hybrid and baseline models are regional earthquake likelihood models [8] defined by a set of expected numbers of earthquakes {*λ*(*j,k*)} in cells of location, indexed by *j*, and magnitude, indexed by *k*, where *j* = 1, ... , *nj*, and *k* = 1, ... , *nk*. The baseline model is denoted by {*λ*0(*j,k*)}, and the hybrid model by {*λH*(*j,k*)}. The *ni* covariates take values in the location cells only and are denoted by {*ρi*(*j*), *j* = 1, ... , *nj*; *i* = 1, ... , *ni*}. The covariates may be real-valued or categorical (e.g., binary) variables. The multiplicative hybrid model is of the form

$$
\lambda\_H(j,k) = \lambda\_0(j,k) \exp(a) \Pi\_{i=1}^{n\_i} f\_i[\rho\_i(j)]. \tag{1}
$$

Here, *a* is a normalising parameter that ensures the total expected numbers of earthquakes in the hybrid model match with the actual number in the catalogue to which it is fitted. The function *fi* converts a cell value of the *i*th covariate into a multiplier to be applied to the corresponding cells of the baseline model. For a real-valued covariate, the multiplier *fi* is a non-negative, monotone, non-decreasing and nonlinear function of the form

$$f\_i[\rho] = \exp(b\_i[\ln(1+\rho)]^{c\_i}), \ (b\_i \ge 0; c\_i \ge 0). \tag{2}$$

It thus preserves the ordering of values in cells but not the ratio of values between cells. For a binary covariate, where *ρ*(*j*) takes only the values 0 and 1, *fi* is of the form

$$f\_i[\rho] = \exp(d\_i \rho)\_\prime \tag{3}$$

where adjustable parameters *di* are allowed to take positive or negative values.

The adjustable parameters in the hybrid model are the normalising parameter *a*, the shape parameters *bi* and *ci* for each real-valued covariate, and *di* for each binary covariate. In fitting a hybrid, the parameters are chosen to optimize the log-likelihood, ln *L*, of the target earthquakes, given by Rhoades et al. [15] as

$$\ln L = \sum\_{n=1}^{N} \ln \lambda\_H(j\_{n\prime} k\_n) - \hat{N}\_{H\prime} \tag{4}$$

in which the target earthquakes occur in cells {(*jn*, *kn*), *<sup>n</sup>* <sup>=</sup> 1, ··· , *<sup>N</sup>*}, and *<sup>N</sup>*<sup>ˆ</sup> *<sup>H</sup>* denotes the total expected number of target earthquakes, i.e.,

$$
\hat{N}\_H = \Sigma\_{j=1}^{n\_j} \Sigma\_{k=1}^{n\_k} \lambda\_H(j,k). \tag{5}
$$

Following Rhoades et al. [6], the corrected information gain per earthquake (IGPEc) of a fitted hybrid model {*λH*(*j*,*k*)} over the baseline model {*λ*0(*j*,*k*)} is calculated using

$$\text{IGPEc} = \frac{-\Delta}{2N} \,\text{}^\prime\tag{6}$$

in which Δ is the change in the corrected Akaike Information Criterion (AICc) between the hybrid model and the baseline model. AICc is defined by

$$\text{AICc} = -2\ln L + 2p + \frac{p+1}{N-p-1},\tag{7}$$

in which *p* is the number of fitted parameters in the hybrid model, and *N* is the number of target earthquakes [16].

Confidence limits on IGPEc are estimated using an adaptation of the *T*-test [15] for comparing one earthquake likelihood model to another based on independent data, as described by Rhoades et al. [6].

For testing of the fitted models on independent data, in which there are no adjustable parameters, the information gain per earthquake (IGPE) is simply the change in loglikelihood between the hybrid and baseline models divided by the number of target earthquakes, and the *T*-test given in [15] applies.

The covariates are defined on the NZ Collaboratory for the Study of Earthquake Predictability (CSEP) testing grid of 0.1-degree rectangular cells [17]. The covariates used in this study are strain rates (BSS, GSS, and HWS), Proximity to the Plate Interface (PPI), Proximity to Mapped Faults (PMF), and Fault in cell (FLT). We have provided details of strain rate covariates in Appendix A. The non-strain-rate covariates are slightly modified compared to those used by Rhoades et al. [10]. We have provided information on these modifications in Appendix B. It is not meant to include the smoothed seismicity covariates from Rhoades et al. [10] because they were fitted to the earthquakes prior to 1986 that are now used for reverse testing. Figure 1 displays the spatial distributions of the non-strainrate covariates.

**Figure 1.** Spatial distributions of the modified (**a**) Proximity to plate interface (PPI) and (**b**) Proximity to mapped faults (PMF) covariates (log of grid cell values). The colour scale shows the logarithm of the relative rate. (**c**) Spatial distribution of updated Fault in cell (FLT) binary covariate using all faults from the New Zealand Community Fault Model [18].

The strain rate covariates (BSS, GSS, and HWS) are included separately in hybrid models as single covariates and with all possible combinations of the other (non-strain-rate) covariates. Hybrid models are fitted to the period 1987–2006, with forward testing over the period 2012–2015 and reverse testing over the period 1951–1986. The spatial distributions of the strain rate covariates are shown in Figure 2.

**Figure 2.** Spatial distribution of (**a**) GSS (1996–2013), (**b**) BSS (1996–2011), and (**c**) HWS (1995–2013) strain rate models in the New Zealand CSEP testing region. The colour scale shows the logarithm of the relative rate.

#### *2.2. Molchan Error Diagram*

We obtain Molchan error diagrams [19,20] for all covariates at different time intervals. Molchan error diagrams are equivalent to Receiver Operating Characteristic curves, which have been widely used in other contexts [21,22]. They can be used to investigate the correlation between a spatially distributed covariate and the spatial distribution of target earthquakes.

The error diagram is based on the concept of an alarm being declared wherever a covariate exceeds a certain threshold. As the threshold is decreased, the proportion of space occupied by alarms increases, and the proportion of target earthquakes outside of the alarm space (i.e., the proportion of "unpredicted earthquakes") decreases. In the error diagram, the proportion of unpredicted earthquakes is plotted against the proportion of space occupied by alarms, but the threshold for declaring an alarm is not shown. The error diagram is, therefore, a monotone decreasing curve.

For an ideal predictor, all target earthquakes are predicted by an alarm occupying almost zero area. The error diagram curve then consists of the line segments joining the point (0,1) to (0,0) and the point (0,0) to (1,0). For a totally uninformative predictor, the proportion of space occupied would be the same as the proportion of predicted earthquakes. The error diagram curve would then be the diagonal line joining (1,0) to (0,1).

The area above the error diagram curve within the unit box is called the Area Skill Score (ASS) [23]. The ASS is a rough measure of the overall potential the covariate has for predicting the locations of target earthquakes. An ideal predictor has an ASS of 1, and a totally uninformative predictor has an ASS of 0.5.

If a covariate is better correlated with the locations of target earthquakes during one time period than another, then the error diagram will vary with time accordingly. Therefore, we plot and compare error diagrams for different periods.

#### *2.3. Scatter Plot of Covariates for the Target Earthquakes*

We have investigated gradual changes over time in the correlation of a covariate with earthquake locations using a simple scatter plot. The scatter plot displays values of a covariate in cells corresponding to target earthquakes versus time. Trends in these values, whether decreasing or increasing, are revealed by a smoothed local regression fit to the data.

#### *2.4. Data*

The target earthquakes in this study are events from 1951 to 2020 with a hypocentral depth between 0 and 40 km and a magnitude greater than 4.95 (GeoNet preferred magnitude) within the CSEP testing region [17]. The events of interest were extracted from the GeoNet catalogue of NZ earthquakes (www.geonet.org.nz; last accessed on 1 July 2022). Magnitudes assigned by GeoNet are mostly local magnitudes, but moment magnitudes are preferred for larger events. The catalogue is consistent with that used by Rhoades et al. [10,11] within the fitting and testing periods of those studies.

#### **3. Results**

#### *3.1. Hybrid Models*

IGPE estimates and 95% confidence intervals of the hybrid models with BSS for the fitting period (1987–2006), forward testing period (2012–2015), and reverse testing period (1951–1986) are presented in Figure 3. Corresponding results with GSS and HWS are presented in Figures 4 and 5, respectively. The IGPEs are relative to the SUP baseline model. To compare the IGPEs of two hybrid models, one needs to obtain the differences between their log-likelihoods for the target earthquakes and perform a *T*-test on those differences. In this study, we have only presented *T*-tests for comparing each hybrid model to the baseline model.

The combination of covariates contributing to each hybrid is indicated by a model number composite. This composite combines the baseline model identifier with those of the selected covariates: "1" for SUP, "2" for PPI, "3" for PMF, "4" for FLT, and "5" for strain rate covariate (BSS, GSS, or HWS). Table 1 contains the name, acronym, and identifier for the baseline model and covariates. For example, the composite "123" represents a hybrid with SUP as baseline and PPI and PMF as covariates.


**Table 1.** Name, acronym, and identifier for the baseline model and covariates.

**Figure 3.** Information gain per earthquake (IGPE) of hybrid models relative to SUP baseline model with the Beavan shear strain rate (BSS) covariate. The combination of covariates contributing to each hybrid is indicated by the model number, which is a composite of the numbers corresponding to the baseline model and selected covariates, as follows: 1-SUP; 2-PPI; 3-PMF; 4-FLT; 5-BSS. (**a**) Fitting period (1987–2006); (**b**) Forward testing period (2012–2015); (**c**) Reverse testing period (1951–1986). Error bars are 95% confidence intervals.

**Figure 4.** Information gain per earthquake (IGPE) of hybrid models relative to SUP baseline model with the GSS strain rate covariate. The combination of covariates contributing to each hybrid is indicated by the model number, which is a composite of the numbers corresponding to the baseline model and selected covariates, as follows: 1-SUP; 2-PPI; 3-PMF; 4-FLT; 5-GSS. (**a**) Fitting period (1987–2006); (**b**) Forward testing period (2012–2015); (**c**) Reverse testing period (1951–1986). Error bars are 95% confidence intervals.

**Figure 5.** Information gain per earthquake of hybrid models relative to SUP baseline model with the HWS strain rate covariate. The combination of covariates contributing to each hybrid is indicated by the model number, which is a composite of the numbers corresponding to the baseline model and selected covariates, as follows: 1-SUP; 2-PPI; 3-PMF; 4-FLT; 5-HWS. (**a**) Fitting period (1987–2006); (**b**) Forward testing period (2012–2015); (**c**) Reverse testing period (1951–1986). Error bars are 95% confidence intervals.

Results for the fitting and forward testing period with the BSS covariate were previously reported by Rhoades et al. [11] and showed that the IGPE was higher for all hybrid models including BSS than for all models not including BSS. Here, with modified non-strainrate covariates for the fitting period (Figure 3a), all hybrids including BSS (composites containing "5") have IGPE > 1.1. The largest IGPE is 1.48 for "1345". On the other hand, all hybrids excluding BSS have IGPE < 1.1. For these, the largest IGPE is 1.00 for "134".

When the BSS covariate is replaced by either GSS (Figure 4a) or HWS (Figure 5a), the results for the fitting period are different. With BSS, "15" (SUP + BSS) has the largest IGPE amongst all the simple hybrids. However, with GSS and HWS, "13" (SUP + PMF) outperforms "15". However, similar to BSS, the largest IGPE values for GSS and HWS are obtained for "1345" (1.08 and 1.01, respectively).

For the forward testing period, the hybrids again perform differently with different strain rate covariates. With BSS (Figure 3b), the hybrids including strain rates are not as dominant as in Rhoades et al. [11], due to the revision of the non-strain-rate covariates (Appendix B). Hybrid "13" performs almost as well as "15", i.e., the IGPEs of both hybrids are similar. With GSS (Figure 4b), "13" performs much better than "15". With HWS, "15" outperforms "13" and all other simple hybrids. Moreover, all hybrid models including HWS outperform most hybrids without HWS except "123".

Hybrids with HWS perform not quite as well as those with BSS in the fitting period (Figures 4a and 5a) but better than those with BSS in the forward testing period (Figures 4b and 5b). Larger IGPEs in the forward testing period may be due to a partial overlap between this period (2012–2015) and the time range of GNSS data used to derive HWS (1995–2013). Hybrids with GSS perform not as well in either the fitting or forward testing period as their counterparts with BSS and HWS (Figure 4a, Figure 5a, and Figure 6a, and Figure 4b, Figure 5b, and Figure 6b, respectively).

**Figure 6.** Molchan error diagrams for covariates (**a**) HWS, (**b**) BSS, (**c**) GSS, (**d**) PPI, (**e**) PMF, and (**f**) FLT, and target earthquakes in the forward test period 2012–2015.

Hybrids with BSS perform poorly in the reverse testing period 1951–1986 (Figure 3c) and even worse than those without BSS, including the SUP baseline (since their IGPEs are negative). For all models without BSS, IGPE is mostly lower than in the fitting and forward testing periods. The largest IGPE for the reverse testing period was 0.27, obtained by hybrid "12", and is only half of the corresponding IGPE for the fitting period.

With GSS, the IGPE of the "15" hybrid is again negative in the reverse testing period, as are the IGPEs for "12345" and "1345". With HWS, the "15" hybrid has a small positive IGPE of 0.03 for the reverse testing period, and all the hybrids have small positive IGPE values (Figure 5c). Thus, hybrid models with HWS perform slightly better than the corresponding models with BSS or GSS in the reverse testing period (c.f. Figures 3c, 4c and 5c). However, overall, the performance of hybrids including strain rates is poor in the reverse testing period.

#### *3.2. Error Diagrams*

A hybrid model can perform poorly in an independent testing period either because the covariates that contribute to it are not well correlated with earthquake locations or because the model is overfitted to a training period in which the covariates are better correlated with earthquake locations than in the testing period. Molchan error diagrams give an indication of the correlation between covariates and earthquake locations in different time periods independent of any model fitting. For the error diagram analysis, we have split the 36-year reverse testing period into three separate 12-year periods.

In the forward testing period 2012–2015, the non-strain-rate covariates PPI and PMF and strain covariates BSS and HWS are well-correlated with the locations of target earthquakes, as shown by the error diagrams deviating far below the diagonal (Figure 6). The ASS for each of these covariates is rather high, between 0.85 and 0.88. HWS has the highest ASS of 0.88. Again, this may be due to the partial overlap between the forward testing

period and the time range of GNSS data used to derive HWS. In contrast, GSS and FLT have smaller ASS values of 0.78 and 0.72, respectively, which implies they are not as strongly correlated with the locations of target earthquakes. This is reflected in the error diagrams not reaching as far below the diagonal as those of the other covariates.

The correlation between the covariates and target earthquake locations varies appreciably over time (Figure 7). The strain rate covariates, especially BSS, show this variation more than the other covariates. Figure 7 compares error diagrams and ASS values for each covariate in five selected time periods—1951–1962; 1963–1974; 1975–1986; 1987–2006; and 2012–2015.

The range of variations in ASS values for a single covariate is an indicator of how widely the correlation varies between time periods. For the PMF covariate (Figure 7a), the correlation is equally high during fitting and forward testing periods (*ASS*1987–2006 = *ASS*2012–2015 = 0.85) and lowest during 1951–1962 (*ASS*1951–1962 = 0.72). For the PPI covariate (Figure 7b), the correlation with earthquake locations is highest during 2012–2015 (*ASS*2012–2015 = 0.84) and lowest during 1963–1974 (*ASS*1963–1974 = 0.70). For the FLT covariate (Figure 7c), the correlation is highest during the forward testing period (*ASS*2012–2015 = 0.72) and lowest during 1975–2006 (*ASS*1975–1986 = 0.56 , *ASS*1987–2006 = 0.57). For the HWS covariate (Figure 7d), the correlation is highest during 2012–2015 (*ASS*2012–2015 = 0.88) and lowest during 1963–1974 (*ASS*1963–1974 = 0.68). For the BSS covariate (Figure 7e), the correlation is highest during 2012–2015 (*ASS*1987–2006 = 0.87) and lowest during 1963–1974 (*ASS*1963–1974 = 0.65). For the GSS covariate (Figure 7f), the correlation is highest during 1987–2006 (*ASS*1987–2006 = 0.81) and lowest during 1951–1962 (*ASS*1951–1962 = 0.60). The low ASS can be explained by the fact that many of the target events in 1951–1962 were located where GSS had a value of zero (Figure 8).

**Figure 8.** Spatial distribution of GSS and target earthquakes during 1951–1962.

The variation of ASS values with time period is plotted for each covariate in Figure 9. No covariate has a consistently higher ASS than the other covariates across all time periods. As can be seen, the highest ASS among the covariates in 1951–1962 is that of PPI, in 1963–1974 that of PMF, in 1975–1986 that of GSS, in 1987–2006 that of BSS, and in 2012–2015 that of HWS. At the other end of the scale, the FLT covariate has the lowest ASS value in all time periods except 1963–1974, in which it shares the lowest value with BSS. Among the strain rate covariates, the highest ASS value is that of the HWS from 1951 to 1962 and 2012 to 2015, GSS from 1963 to 1974 and 1975 to 1986, and BSS from 1987 to 2006. Figure 9 also shows that the HWS covariate has a higher ASS value than BSS in all time periods except 1987–2006—the fitting period for the multiplicative models.

#### *3.3. Smoothed Scatter Plots*

To give a better insight into changes over time in the correlation of covariates with earthquake locations, we plotted the covariate values for the cells corresponding to target earthquakes versus time over the period 1951–2020. We then fitted a trend curve using locally weighted scatter plot smoothing (lowess) based on Cleveland's method in [24]. We observed clear trends in the strain rate covariate values corresponding to target earthquake locations over the period 1951–2020, as shown by Figure 10.

**Figure 9.** Area Skill Score (ASS) of Molchan error diagram for covariates over different times.

**Figure 10.** Smoothed scatter plots of (**a**) GSS, (**b**) BSS, and (**c**) HWS strain rate covariates for targeted earthquakes against time from 1951 to 2020.

The trend is stronger for BSS and HWS than for GSS. However, for all three covariates, the smoothed trend line starts from its lowest level in 1951, rises gradually to a peak in 1999, and then tails off somewhat through to 2020. However, the HWS curve also tails off less than the others. These results confirm and elaborate on the trends over time that was observed in Figures 8 and 10. The HWS model was extracted from the longest period of recordings (1995–2013) and a larger number of GNSS stations compared with the GSS and BSS models. In addition, unlike BSS and GSS that were both based on an older method of Haines and Holt [25–27], the HWS model was derived using the novel VDOHS technique [14]. Such differences could explain the differences observed in the smoothed scatter plots versus time.

We observed that FLT was the only covariate to show a clear trend with the magnitude of target earthquakes. Figure 11 plots the FLT values versus the magnitudes of the target earthquakes from 1951 to 2020. The locally weighted smoothed curve, shown by the solid black line, reveals an increasing trend with magnitude. This may not necessarily imply that future large earthquakes are more likely to occur on known faults than smaller earthquakes. It may instead reflect the fact that large earthquakes cause surface ruptures more frequently than small ones; consequently, new ruptures over the past 70 years have revealed the associated faults and allowed them to be mapped.

**Figure 11.** Fault covariate values for targeted earthquakes against magnitude, 1951–2020. The solid black line shows the smoothed scatter plot. The blue triangles depict the values of the FLT covariate for each target earthquake.

#### **4. Discussion**

The low information gain of hybrid models including strain rate covariates in reverse testing is consistent with the relatively low ASS values of the strain rate covariates and low strain rate values associated with target earthquakes from 1951 to 1986. The ASS values of the strain rate covariates in the reverse testing period are not only lower than in the fitting period but are also lower than the ASS values of the PPI and PMF covariates during the reverse testing period (Figures 7 and 8). The contrast of ASS values between the fitting and reverse testing periods is particularly strong for BSS. This explains the negative IGPEs for hybrids with BSS in the reverse testing period (Figure 3c).

Similarly, the large IGPE for hybrids with HWS in forward testing is consistent with the high ASS for HWS in 2012–2015 and with the slow tailing off between 2010 and 2020 of the trend in HWS associated with target earthquakes. Evidently, if the forward testing period were extended through to 2020, hybrids including HWS would have continued to perform well.

Differences between the three present strain rate models can be attributed to the differences in the geodetic datasets on which they were based and the evolving techniques for estimating strain rates from the geodetic data. The most recently developed of the three models is the HWS model, which was published in 2020. HWS was derived from the most complete data and up-to-date techniques and showed the best correlation with the locations of the most recent earthquakes (Figure 7). However, all three strain rate models exhibit a disappointingly low correlation with the locations of earthquakes from 1951 to 1986.

The implicit assumption behind the fitting of hybrid models to estimate long-term earthquake rates is that the correlation between the covariates and the locations of target earthquakes is stable over time. Alternatively, if not stable over time, the correlation during the fitting period should be indicative of the long-term correlation. The correlations for PPI and PMF appear to be more stable over time than those for the strain rate covariates, based on the ASS variations. PPI and PMF represent features that are expected to affect earthquake occurrence over a long time frame. The greater instability of the correlation for strain rate covariates may be attributed to the short-term non-stationarity of earthquake occurrence. The strain rate covariates are most highly correlated with earthquake occurrence within the time period of the data on which they are based or within a few years after (Figure 7).

Without a dense GNSS network and extensive historical and paleo-seismicity records, it is not possible to pinpoint temporal and spatial variations of strain rates. Recently, Iezzi et al. obtained high temporal histories of slip rates for three parallel normal faults in central Greece with progressive clusters using in situ 36Cl cosmogenic dating [28]. They also had access to associated GNSS recordings of a dense network. They compared regional decadal geodetic strain rates with strain rates derived from slip-rate pulses over a few thousand years on each fault. They found that strain builds up across a fault system over thousands of years, but not all individual faults actively slip at the same time. They discussed how the available geodetic data alone was not capable of pinpointing the switchover of slip from one fault to another over time. In comparison, New Zealand has neither a dense GNSS network nor a comparable record of fault slip histories. The density of GNSS stations has been progressively improving but is still sparse in some parts and limited to onshore regions [14]. Obviously, this affects the accuracy of the strain rate models in NZ offshore regions, where many of the earthquakes occur.

Given the low earthquake location–strain rate correlations in the reverse testing period here, it is expected that there will be future periods during which such correlations are similarly low.

In fitting hybrid models for estimating long-term seismicity rates, it seems best to use the longest possible fitting period subject to earthquake catalogue quality considerations. The NZ catalogue covers about 180 years of historical and instrumental records. With the improvement of the NZ seismographic network, the magnitude of completeness has lowered to 5.0 for the recent 70 years of the GeoNet catalogue. Such a level of magnitude completeness fulfils the magnitude threshold of 4.95 applied to the target earthquakes for this study. However, we cannot be confident that a 70-year period, or any other fixed period, is long enough to represent the long-term spatial distribution of earthquakes. The main purpose of the hybrid fitting is to realistically weight the information offered by the available covariates, none of which depend directly on the catalogue. Here we showed that it would be unwise to base this weighting on a 20-year subset of the available catalogue. Performing tests with the 180-year catalogue could potentially help us to understand how adequate a 70-year catalogue is for fitting such hybrid models.

The robustness of our results against time-varying errors in earthquake locations could be questioned. However, since the strain rate models and all other covariates were defined on a rather coarse grid of 0.1 × 0.1-degree cells, they naturally can accommodate location errors without having much effect on the correlation with earthquake locations. We examined the robustness of the correlations between all covariates and target earthquake locations by

introducing synthetic location errors within the range reported by Bondár et al. [29] to be the most probable values for the ISC-GEM catalogue. The location errors were added to the epicentres of target events from 1987 onwards. We tested the hypothesis that events prior to 1987 had larger location errors due to the seismograph network's sparsity, and that could be the reason for low correlations between the covariates and target earthquake locations. However, this hypothesis was rejected; despite the added location errors for all the covariates, the main features of the smoothed curves were preserved.

#### **5. Conclusions**

We analysed the performance of multiplicative hybrid models in forward and reverse testing periods. To complement the performance analysis, we also obtained Molchan error diagrams and smoothed scatter plots of strain rate covariates versus time. All the analyses reveal strong variability over time in the correlation of strain rates with earthquake locations. These results suggest that strain rate models based on the available geodetic data for only a few decades may not be good indicators of where earthquakes occur over a long time frame. Therefore, the longest available high-quality earthquake catalogue should be used in constructing hybrid models to estimate long-term earthquake rates for seismic hazard analysis.

**Author Contributions:** Conceptualization, S.J.R., D.A.R. and M.C.G.; methodology, S.J.R. and D.A.R.; software, S.J.R., D.A.R. and C.R.; formal analysis, S.J.R. and D.A.R.; writing—original draft preparation, S.J.R. and D.A.R.; writing—review and editing, S.J.R., D.A.R., C.R. and M.C.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Strategic Science Investment Fund (SSIF) and the National Seismic Hazard Model Revision Project of the Ministry of Business, Innovation and Employment of New Zealand.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The New Zealand Earthquake Catalogue was obtained from GeoNet. Available online: http://www.geonet.org.nz (accessed on 30 May 2022). Data on earthquake faults and slip rates were obtained from the New Zealand Community Fault Model v1.0 (Data Set), GNS Science, https://doi.org/10.21420/NMSX-WP67 (accessed on 30 May 2022).

**Acknowledgments:** We acknowledge the New Zealand GeoNet project and its sponsors EQC, GNS Science, LINZ, NEMA, and MBIE for providing the data used in this study. Kiran Kumar Thingbaijam and David Burbidge provided helpful internal reviews of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A**

#### *Appendix A.1 Strain Rate Covariates*

Beavan and Haines [30] utilized GNSS data from 9 geodetic networks across NZ, comprising 362 GNSS stations. They calculated continuous horizontal velocity and associated strain rate fields at the Earth's surface in NZ between 1991 and 1998. They adopted the method proposed by Haines et al. [26], which enhanced NZ's known deformations and discovered new features in both the North and South Islands.

John Beavan (1950–2012) directed most of the GNSS data acquisition in NZ [14]. His last shear strain rate model of NZ was based on survey-mode GNSS data collected in NZ between January 1996 and February 2011. This model was based on the method of Beavan and Haines [30] and was submitted to Land Information New Zealand (LINZ) to be used in updating the geodetic and cadastral systems following the 4 September 2010 Darfield earthquake. This was the shear strain rate (BSS) used by Rhoades et al. [11].

The Global Strain Rate Model (GSRM v.2.1) of Kreemer et al. [13] is an improved version of the GSRM-1 Kreemer et al. model in [2] with a larger input dataset and improved spatial resolution. The maximum shear strain rate from GSRM v.2.1 is referred to here as the GSS model. Like the BSS model, the GSS model was obtained using the Haines and Holt method [25,27,30] to model the strain rate field. GSS is derived from worldwide GNSS data between 1 January 1996 and 31 December 2013, including velocities from 233 published and unpublished regionally focused studies. Table A1 lists the NZ-focused studies used in the GNSS compilation of the GSS model [13]. Apart from in the southeast South Island region, the GSS model uses less up-to-date data than the BSS model.

The NZ strain rate model of Haines and Wallace [14] was produced from GNSS interseismic velocities at 918 stations between 1995 and 2013 [31]. This model was based on the novel method of Haines et al. [32]. They used a physics-based approach to invert the GNSS displacement fields for Vertical Derivatives of Horizontal Stress (VDoHS) rates. This method uses constraints from linear elasticity in place of artificial smoothing. This is to ensure physical plausibility, higher resolution, and an improved signal-to-noise ratio compared to the commonly used methods of Haines and Holt [25] and Beavan and Haines [30]. The one modification made to the published Haines and Wallace model of [14] is one implemented by the NSHM Geodesy Working Group: a sill model of rapid magmatic contraction in the Taupo Volcanic Zone (TVZ) was removed from the GNSS velocities, and then the velocities were re-inverted using the Haines et al. [32] method [31]. The maximum shear strain rate from this modified model is the scalar quantity that we use in our hybrid modelling and is henceforth referred to as HWS [33].

**Table A1.** Data from NZ GNSS campaigns used in constructing the GSS rate model.


#### **Appendix B**

*Appendix B.1 Proximity to the Plate Interface (PPI)*

A PPI covariate was initially introduced by Rhoades et al. [10]. Here we use a modified version of PPI that differs in some minor details ([33]; Figure 1a). PPI is based on the notion of proximity to the plate interface. A map of the interface between the Australian and Pacific plates in the Hikurangi subduction zone [37] is used to define the North Island plate boundary. South of the Hikurangi subduction zone, the plate boundary is defined by the Hope Fault, the Alpine Fault, and the Fiordland subduction interface, as given in the 2010 NZNSHM fault model [38]. In estimating the proximity of a cell to the interface, no allowance is made for the possibility of different slip rates on different patches of the plate boundary. The plate boundary is represented by a set of *nI* points on the interface at approximately 1 km spacing, and a cell is represented by a set of *nc* points at its middle latitude and longitude coordinates with depths from 1 to 39 km at a spacing of 2 km. The PPI covariate *ρ*(*j*) in the *j*th grid cell is defined by

$$\rho(j) = \frac{1}{n\_{\mathcal{C}}} \sum\_{l=1}^{n\_{\mathcal{C}}} h(j, l), \tag{A1a}$$

where *h*(*i*,*l*) is an index of the proximity of the *l*th point in the *j*th cell to *i*th point on the interface:

$$h(j,l) = \frac{1}{\left[d\_I^2 + \min\left(\Delta(i,l)^2\right)\_{i=1,\dots,n\_l}\right]^{3/2}}\tag{A1b}$$

where Δ(*i*,*l*) is the distance in km from [*xI*(*i*), *yI*(*i*), *zI*(*i*)] to *xj*(*l*), *yj*(*l*), *zj*(*l*) and *dI* is a constant smoothing distance, taken as 10 km.

#### *Appendix B.2 Proximity to Mapped Faults (PMF)*

The concept of proximity to mapped faults was initially introduced by Rhoades and Stirling [39]. It takes account of the distance from the *j*th cell to all fault elements and their associated slip rates, assumed uniform on each fault segment. As in [39], each fault plane is divided into numerous point sources closely spaced at intervals of 1 km. Thus the set of fault segment planes and associated slip rates is expanded into a much larger set of point sources, at longitude, latitude, and depth coordinates [*xF*(*i*), *yF*(*i*), *zF*(*i*)], *i* = 1, ··· , *nF*, and with associated slip rates *ri*, *i* = 1, ··· , *nF*. Similarly, a cell with dimensions of 0.1 degrees in latitude and longitude and 40 km in depth is expanded into a representative set of points within it *xj*(*l*), *yj*(*l*), *zj*(*l*) , *l* = 1, ··· , *nc* on a cuboidal grid spaced at 0.033 degrees in latitude and longitude and 5 km in depth. Following [10], the PMF covariate *ρ*(*j*) in the *j*th grid cell is defined by

$$\rho(j) = \frac{1}{n\_{\text{ct}}n\_f} \sum\_{i=1}^{n\_f} \sum\_{l=1}^{n\_c} h(i, l), \tag{A2}$$

where *h*(*i*,*l*) is an index of the proximity of the *l*th point in the *j*th cell to slip on the *i*th fault point source:

$$h(i,l) = \frac{r\_i}{\left[d\_F^2 + \Delta(i,l)^2\right]^{\frac{3}{2}}} \; \prime \tag{A3}$$

where *ri* is the slip rate in mm/yr, Δ(*i*,*l*) is the distance in km from [*xF*(*i*), *yF*(*i*), *zF*(*i*)] to *xj*(*l*), *yj*(*l*), *zj*(*l*) , and *dF* is a constant smoothing distance. The value of *ρ*(*j*) is, therefore, high for cells near to, or intersected by, mapped faults and low far away from mapped faults. Further, *ρ*(*j*) is greater when the slip rates of nearby mapped faults are greater.

Here we used a modified version [33]. In the modified PMF fault, sources are planes rather than linear traces and distances are measured in three dimensions rather than two. Rhoades and Stirling [39] found an optimal smoothing distance of 1.0 km for a PMF model defined in continuous space. Rhoades et al. used a smoothing distance of 2.5 km, but their PMF has since been found to have small-scale variations between adjacent cells due to the smoothing distance being too small [10]. Here we use a larger smoothing distance of *dF* = 10 km, so that the denominators in Equation (5) are not overly sensitive to the exact positions of points on the fault surfaces and within the cells. The PMF covariate used here (Figure 1b) also uses updated fault information from the New Zealand Community Fault Model [18].

#### *Appendix B.3 Fault in Cell (FLT)*

The fault in cell (FLT) covariate [10] is defined to have the value 1 if a fault intersects the cell anywhere between the surface and a depth of 40 km, given the surface trace of the fault, its dip angle and estimated width. Otherwise, it takes the value 0. The FLT covariate is shown in Figure 1c [33] and has been updated using all faults from the New Zealand Community Fault Model [18].

#### **References**


**George Kaviris 1,\*, Angelos Zymvragakis 1, Pavlos Bonatis 2, Vasilis Kapetanidis <sup>1</sup> and Nicholas Voulgaris <sup>1</sup>**


**Abstract:** The Gulf of Corinth (Central Greece) is one of the most rapidly extending rifts worldwide, with its western part being the most seismically active, hosting numerous strong (M ≥ 6.0) earthquakes that have caused significant damage. The main objective of this study was the evaluation of seismic hazard through a probabilistic and stochastic methodology. The implementation of three seismotectonic models in the form of area source zones via a logic tree framework revealed the expected level of peak ground acceleration and velocity for return periods of 475 and 950 years. Moreover, PGA values were obtained through the stochastic simulation of strong ground motion by adopting worst-case seismic scenarios of potential earthquake occurrences for known active faults in the area. Site-specific analysis of the most populated urban areas (Patras, Aigion, Nafpaktos) was performed by constructing uniform hazard spectra in terms of spectral acceleration. The relative contribution of each selected fault segment to the seismic hazard characterizing each site was evaluated through response spectra obtained for the adopted scenarios. Almost all parts of the study area were found to exceed the reference value proposed by the current Greek National Building Code; however, the three urban areas are covered by the Eurocode 8 regulations.

**Keywords:** PGA; PGV; return period; PSHA; stochastic seismic hazard assessment

#### **1. Introduction**

The Gulf of Corinth, located in Central Greece, is a continental rift with an extension rate of 6–16 mm/yr [1], one of the highest values known worldwide for this type of tectonic structure. Its main axis is aligned in an approximately E–W direction, flanked by major normal faults at its southern and northern shoulders. Major earthquakes (Mw ≥ 6.0) have occurred in the past, both at the eastern [2] and western parts of the gulf. However, the Western Gulf of Corinth (WGoC) is the most active in terms of observed microseismicity [3,4], as it exhibits a higher extension rate—i.e., ~10.8 mm/yr near Aigion—than its eastern counterpart (~5.5 mm/yr) [5]. This has led to the WGoC being studied extensively during past decades, and it is closely monitored by the Corinth Rift Laboratory (CRL) local seismological network [6]. The CRL has also obtained the status of an international EPOS Near Fault Observatory (NFO), providing high-resolution multidisciplinary data and products [7–9].

The WGoC is bordered by major north-dipping normal faults in the south, with measured dipping angles on the surface ranging from 50◦ (e.g., Psathopyrgos fault [10]) to ~60◦ (e.g., Aigion fault [11]) and up to 65–70◦ (e.g., Lakka fault [12]), with average slip rates in the Late Quaternary of the order of 2.4–5.0 mm/yr for the Aigion and East and West Helike faults and 2.0–3.5mm/yr for the Psathopyrgos fault [10]. On the northern coast, steep south-dipping antithetic normal faults have been mapped, such as the Marathias fault, dipping at 55◦ [13], and the Trizonia fault, dipping at 64–72◦, the latter with an average slip rate of 0.36–0.44 mm/yr over the last ~130 ka [14]. Offshore E–W-trending normal faults, along with smaller structures, oblique to the main axis of the rift with significant

**Citation:** Kaviris, G.; Zymvragakis, A.; Bonatis, P.; Kapetanidis, V.; Voulgaris, N. Probabilistic and Scenario-Based Seismic Hazard Assessment on the Western Gulf of Corinth (Central Greece). *Appl. Sci.* **2022**, *12*, 11152. https://doi.org/ 10.3390/app122111152

Academic Editor: José A. Peláez

Received: 14 October 2022 Accepted: 1 November 2022 Published: 3 November 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

strike-slip components, have also been documented [10,14,15]. Average slip rates during the Quaternary for major offshore faults in the WGoC—i.e., the North and South Eratini and West and East Channel faults—are estimated to be between 0.4 mm/yr and over 1.4 mm/yr [10].

Seismicity in the WGoC is mainly concentrated at focal depths between 5 and 10 km [6], with the vast majority of fault plane solutions indicating dominant E–W normal faulting [4,16]. However, oblique-normal slip has also been reported, usually related to seismic swarms [7,17] or to the activation of older structures [18]; i.e., the External Hellenides, which crosscut the gulf in a roughly N–S direction. Although such faults are not optimally oriented to the regional crustal stress field [19], slip can be facilitated by the intrusion of high pore-pressure fluids in the fault network, which play a significant role in the evolution of swarms in the area [20,21]. Towards the western edge of the rift, the crustal stress field changes [19], favoring a dextral strike-slip regime along the SW–NE-trending Rion–Patras fault. The relocated seismicity presented in Figure 1 shows that the bulk of activity mainly occurs along the rift axis, between Trizonia Island and Aigion, with hypocenters deepening towards the north, developing a north-dipping low-angle detachment zone in which the major faults of the WGoC are rooting [22]. Several earthquakes have been located at intermediate depths of 50–60 km, with resolved focal mechanisms representing reverse faulting (Figure 1). These are indicative of seismic activity related to the oceanic slab at the northwestern end of the Hellenic subduction zone [23], with ascending fluids due to its dehydration possibly playing a role in the triggering of earthquake swarms at shallow depths [22,24].

**Figure 1.** Seismotectonic map of the western Gulf of Corinth (WGoC). Solid circles and stars (for Mw ≥ 5.0) denote earthquakes with Mw ≥ 4.0 from the 1900–2009 catalogue of Makropoulos et al. [25], extended in this study up to 2019, with symbol size proportional to the magnitude (Mw) and color according to the focal depth. Selected microseismicity, presented with small, hollow, red circles, is adopted from the 2000–2015 relocated catalogue of Duverger et al. [3]. Focal mechanisms of significant events, presented as beachballs, are adopted from [26–32], as well as the databases of the NKUA, NOA, CMT, and ISC. Blue squares mark the macroseismic epicenters of historical earthquakes (1000–1899) acquired from the SHEEC database [33]. Fault lines are after [10,14,34,35]. Faults marked with red labels were used for the stochastic seismic hazard modeling in this study.

Despite the constant presence of microseismicity, certain parts of major faults in the WGoC appear to be locked, accumulating stress that is released in major events. The last destructive earthquake in the area was the 15 June 1995 Ms = 6.2 event, which severely damaged the town of Aigion [29]. The source of this earthquake was a north-dipping, low-angle (33◦) normal fault, with more recent data excluding its association with the East Helike fault [22]. The rupture surface of the 1995 event appears to be locked, with the observed microseismicity stopping abruptly at the eastern end of the study area (Figure 1). Other major earthquakes have also occurred in the past at the eastern end of the WGoC, with epicenters towards the northern coasts. Significant earthquakes have struck the WGoC in the historical era (i.e., before 1900). One of the more notable was the 373 BCE earthquake, with an estimated minimum magnitude of M = 6.6 [36], which caused a tsunami that destroyed the ancient city of Helike [37]. An earthquake of similar macroseismic characteristics struck the same area on 26 December 1861 (Mw = 6.5 according to Boiselet [38]), triggering a tsunami and causing a strip of plain to submerge under the sea [37]. Both earthquakes are attributed to the north-dipping Helike normal fault, most likely its eastern branch [22]. Recently, the contemporary town of Helike was in the epicentral area of a swarm that started in May 2013, persisting for several months and involving a few earthquakes with magnitudes ranging from 3.3 to 3.7, but which was attributed to activity near the root of Pirgaki fault [17,39]. Other historical earthquakes that struck the town of Aigion occurred in 1748, 1817, and 1888, with estimated magnitudes of 6.3–6.6 [33]. The area of Patras was affected by historical earthquakes in 1785, 1804, and 1806, with magnitudes between 6.1 and 6.4, while the town of Nafpaktos was struck by earthquakes in 1462, 1703, 1714, 1756, and 1831, with magnitudes in the range from 5.9 to 6.5 [33]. In the instrumental era (1900 to today), before the 15 June 1995 earthquake, the previous major event with a magnitude ~6 occurred in the WGoC on 30 May 1909, while intermediate-depth events with magnitudes of 6.4 and 5.9 occurred in 1965 and 1993, respectively; the epicenters of these events were located at the eastern end of the study area.

Although no Mw ≥ 6.0 earthquake has struck the WGoC since 1995, moderate magnitude events (Mw ≈ 5) have occurred in recent years. On 14 July 1993, an Mw = 5.5 event with a strike-slip focal mechanism [30,40] struck near Patras, causing considerable damage. Although the locally mapped structures trend SW–NE, it has been suggested, based on the aftershock distribution, that the earthquake ruptured a NW–SE sinistral-slip fault [41]. A notable case is the earthquake doublet of January 2010 near Efpalio, with an Mw = 5.3 event being followed by an Mw = 5.2 four days later, both related to blind E–W-trending normal faults [32,42,43]. This activity likely played a role in the triggering of another Mw = 5.1 earthquake on 7 August 2011 to the west, near Nafpaktos, characterized by a similar normal faulting style. Later, in 2013–2014, the WGoC presented strong signs of microseismic activity, with an earthquake swarm initiating in September 2013 at the offshore region between Nafpaktos and Psathopyrgos that slowly migrated eastwards, triggering earthquake clusters and culminating in an Mw = 5.0 event, offshore of Aigion, on 7 November 2014 [21,44]. A second seismic swarm excitation, characterized by bilateral spatiotemporal migration in the offshore area between Nafpaktos and Psathopyrgos, involved an Mw = 4.9 earthquake on 21 September 2014, activating a patch related to the westernmost edge of the Psathopyrgos fault, near its junction with the Rion–Patras fault zone [3]. The most recent significant activity in the WGoC was the 2020–2021 seismic sequence [7,45,46], which started on 23 December 2020 with a moderate event near Marathias and evolved in three stages, first triggering earthquake clusters towards Trizonia Island involving some strikeslip events and later producing an Mw = 5.3 earthquake offshore of Psathopyrgos that exhibited peculiarities in terms of its complicated rupture characteristics [47].

In this study, we reassess the seismic hazard in the WGoC utilizing a well-established methodology and incorporating recent data. The term "seismic hazard" describes the expected ground motion level for potential earthquake occurrences within a study area. The methodologies used to evaluate seismic hazard are grouped into two major categories: probabilistic and scenario-based. The first exploits data from historical and instrumental

earthquake catalogues and integrates a seismotectonic model for its application. Subsequently, given an earthquake occurrence model, results for the maximum expected ground motion for selected return periods are generated. The second methodology is based on individual earthquake scenarios, with which synthetic waveforms are produced and analyzed. Strong ground motion may be a catalyst for other secondary hazards, the most important of which are rock falls, landslides, and liquefaction phenomena. Therefore, the visualization of the spatial distribution of peak ground motion parameters can be of great importance for scientists and decision makers. The advantage that the study area offers is the wealth of earthquake data. The latter is critical because seismic hazard can be assessed without adopting algorithms that cope with this difficulty, such as the one proposed by Meirova et al. [48]. A plethora of seismic hazard studies have been conducted for the Greek territory, either as a whole or for specific areas, such as [49–58]. The applied methods for the assessment of seismic hazard are described in detail in Section 2.

Patras, the third largest city in Greece and one of its largest harbors, is situated at the western edge of the study area, on the southeastern shores of the Patras Gulf. The nearby town of Rion hosts important infrastructure, such as the University of Patras and the University General Hospital of Patras, while the Rion–Antirrion bridge connects Peloponnese with Central Greece. Other important towns on the shores of the WGoC include Aigion and Nafpaktos. As the seismic risk strongly depends on the exposure to seismically prone areas, the aforementioned city, towns, and infrastructure may suffer increased economic losses from potentially destructive earthquakes, both in terms of human lives and structural damage. As a consequence, the proper assessment of seismic hazard in the WGoC area, both probabilistic and with specific earthquake scenarios, is required for the subsequent estimation of seismic risk and consideration of means for its mitigation.

#### **2. Materials and Methods**

Seismic hazard in the study area was evaluated using two approaches: (a) a Probabilistic Seismic Hazard Assessment (PSHA) and (b) a stochastic approach to worst-case seismic scenarios for known active faults. PSHA remains the primary framework for assessing seismic hazard and is based on prior knowledge of seismicity at a given place. On the other hand, during past decades, finite-fault stochastic ground motion simulations have been proven to be a powerful tool to reliably estimate strong ground motion parameters, such as the Peak Ground Acceleration (PGA) and Peak Ground Velocity (PGV). The combination of these methodologies can provide a holistic evaluation of seismic hazard, which in turn can contribute significantly to the assessment of seismic risk and its mitigation at a given area.

#### *2.1. PSHA*

Concerning PSHA, the approach proposed by Cornell and McGuire [59,60] was followed, utilizing the instrumental earthquake catalogue from Makropoulos 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 Mw ≥ 4 magnitude earthquakes have been recorded in the study area since 2019, the largest being the Mw = 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 Vamvakaris 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 (M0; also considered as the magnitude of completeness (Mc)), (c) the average annual exceedance rate of Mc (λ(Mc)), and (d) the maximum expected magnitude (Mu). The Mc and b-values were estimated for each zone using the maximum curvature (MAXC) method from Wiemer and Wyss [67] and the maximum likelihood function from Aki [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 λ(Mc) was calculated through the analysis of the earthquake catalogue and Mu was estimated via the Robson–Whitlock–Cooke (R-W-C) procedure [72,73]. Regarding the geometry data, the depth of each seismic zone was estimated through depth histograms constructed using the focal depths reported in the earthquake catalogue. It must be noted that ESHM13, ESHM20 and VAM16 provide values concerning the aforementioned parameters, each using a different earthquake catalogue. However, in the present study the b-value, Mc, λ(Mc) and Mu were recomputed for each ASZ, im-plementing the updated earthquake catalogue of Makropoulos et al. [25] to retain ho-mogeneity regarding the seismicity and geometry data. Consequently, only the ASZ were adopted from the three zonation models. The computed seismicity and geometry data are depicted in Table S1.

Ground Motion Prediction Equations (GMPEs) proposed for the Greek territory were adopted for the calculation of the maximum expected ground motion in terms of peak ground acceleration and velocity (PGA, PGV). In particular, for PGA, the GMPEs from [74–78] were employed. Sakkas [77] has not proposed a GMPE for PGV estimation. Furthermore, the model from Skarlatoudis et al. [75] was replaced with an upgraded GMPE for PGV, which was taken into account. As area source zones were herein adopted, it was only feasible to use epicentral distances. As a result, more recent GMPEs [79] that utilize different distance metrics were not considered. In addition to PGA and PGV, Spectral acceleration (Sa) values for different periods (T) were calculated for the three most important urban areas of the WGoC—namely, Patras, Nafpaktos, and Aigion—using the GMPE from Danciu and Tselentis [76]. All GMPEs (except for the one from Margaris et al. [74]) include a term related to the type of focal mechanism, taking a value of 0 or 1 for normal and strikeslip/thrust ruptures, respectively. As a seismic zone could potentially include any type of focal mechanism, it was deemed necessary to compute their respective participation rates. Thus, each GMPE was calibrated with the exact percentage of focal mechanism types in each zone, generating weight-specific GMPEs (hereafter, w-sGMPEs). The final results (PGA, PGV for return periods of 475 and 950 years, and Sa for a return period of 475 years) were computed via an equal-weighted logic tree, where each secondary branch was

a w-sGMPE (except [74]) and each primary branch was a different seismotectonic model (Figure 2). The reason for assigning equal weights to the primary branches was that none of the three zonation models outweighs the others. There are uncertainties that would be increased if a higher weighting factor would be assigned to the ASZ of ESHM13, ESHM20 or VAM16. For example, there are zones including both the western and eastern Gulf of Corinth in ESHM13 and ESHM20 (zones 13 and 1 in Figures S1a and S1b, respectively). In addition, certain small zones of the VAM16 model (zones 5, 13 in Figure S1c), result in an insufficient number of earthquakes, not adequate to reliably compute the seismicity parameters. The software used for PSHA was R-CRISIS and, specifically, its newest version (V20.0) [66].

**Figure 2.** The logic tree constructed to calculate PGA, PGV, and Sa. All branches have the same weighting coefficient. Results were generated for return periods of 475 and 950 years for both PGA and PGV and 475 years for Sa. Abbreviations: VAM16: Vamvakaris et al. [65], CHO18: Chousianitis et al. [78], SAK16: Sakkas [77], DAT07: Danciu and Tselentis [76], SKA07: Skarlatoudis et al. [75], SKA03: Skarlatoudis et al. [75], MAR02: Margaris et al. [74].

#### *2.2. Seismic Scenarios*

We compiled a set of seismic scenarios to assess the seismic hazard in the broader WGoC area in terms of maximum expected strong ground motion from specific seismic sources. These scenarios were generated through stochastic simulations for a predicted maximum magnitude (Mmax) for well-studied known faults in the area. Simulations were performed using a stochastic finite-fault model based on a dynamic corner frequency approach with the EXSIM code [80,81]. The modeling strategy was based upon the discretization of the earthquake fault plane into smaller subfaults, each of which was considered as a potential earthquake source. The point-source stochastic method [82] was employed to generate synthetic time series for each subfault. The summation of the individual contributions of each subfault, along with a suitable time delay, led to the final ground motion parameters at the sites of interest. This approach has been widely implemented worldwide [83], as well as in Greece [84–88].

The anticipated strong ground motion at a given site is a result of a complex physical process that includes the relative contribution of source, path, and local site effects. Source effects describe the characteristics of the accumulated strain release from the fault, which, when released, results in the generation of an earthquake. These include, among others, the source dimensions (length, width), the moment magnitude (Mw), the slip distribution for the causative fault, and the stress parameter (Δσ). After the earthquake nucleation and rupture, seismic waves travel through the Earth's crust; therefore, path properties must also be taken into account. Their propagation is most strongly affected by two path parameters; i.e., the geometrical spreading and the anelastic attenuation, which is controlled

by the quality factor (Q). Ultimately, the local site properties play an important role in the resulting surface ground motions, given that the impact of the shallower layers may lead to significant amplification of the seismic waves. These effects are treated by the EXSIM code using the kappa (κ0) parameter [89] and user-defined soil amplification factors.

As a first step to obtain high-resolution model parameters for the desired earthquake scenarios, we performed a stochastic simulation of the most recent strong (Mw ≥ 6.0) earthquake in the study area; namely, the 1995 Aigion MS = 6.2 (Mw = 6.5 according to the Global Centroid Moment Tensor (CMT) project) mainshock [29]. Modeling of past earthquakes in the study area can significantly help in constraining the path and site components of strong ground motion through comparison with GMPEs, as well as with real strong motion data. In the case of the Aigion 1995 mainshock, 17 recordings from seismic stations up to an epicentral distance of 140 km are available, allowing the calibration of the synthetic results through an iterative procedure to achieve the best fit. The final model parameters are summarized in Table 1. The causative fault's dimensions and geometry (strike and dip) were adopted following [29] and the upper edge of the fault plane was set according to the model from Console et al. [90], which is in compliance with the geometry of the seismogenic layer of the study area. The slip distribution onto the fault plane was modeled through a random slip pattern, given that detailed finite-fault slip inversions are not available for this case. This approach has been also followed by other authors for this mainshock [86,88]. Finite-fault discretization into smaller subfaults was performed using the empirical relationship proposed by Beresnev and Atkinson [91] (Equation (1)):

$$
\log \Delta \text{l} = 0.4 \text{M} - 2,
\tag{1}
$$

where Δl is the length of each subfault and M is the moment magnitude of the mainshock.

The Δσ is a parameter closely related to the actual stress drop and slip velocity due to an earthquake occurrence, but it does not include the natural context of stress drop [92]. It is generally considered to weakly scale with moment magnitude, but in a more regional extent their interconnection appears to be more profound [93]. It was determined by employing an iterative procedure, comparing the synthetic ground motions with recorded ones, as well as with estimated peak values from well-established GMPEs [94,95] appropriate for the study area. The Δσ value of 56 bars, which is routinely used as a mean value for earthquakes in Greece [96,97], was used initially and, afterwards, various other values were tested to find the best fit. Ultimately, Δσ = 30 bars was deemed appropriate and was adopted in the simulations. This deviation from the mean value was foreseen, given that low stress drop has been documented for this earthquake [98], as well as for the study area, especially in comparison with the eastern part of the gulf where stress drop values appear to be larger [99]. Regarding the path effects, the geometrical spreading model, as described by Atkinson [100], was adopted, which divides the slope of the attenuation relation in three distance intervals (<70 km, 70–130 km, and >130 km), reflecting the dominant type of wave in the seismic signal. The anelastic attenuation, on the other hand, is described by the quality factor proposed by Hatzidimitriou [101] for the broader Aegean region (Equation (2)):

$$\mathbf{Q(f)} = 100 \mathbf{f^{0.8}},\tag{2}$$

Lastly, to account for the local site conditions for synthetic ground motions, the average κ<sup>0</sup> value of 0.044 for class C (NEHRP) sites in Greece [102] was adopted, along with the corresponding amplification factors. The selection of class C soil characterization was based on local soil data (e.g., Vs30 profiles) for the major cities in the area and geology data from the locations of the available seismic stations. In order to calibrate and validate our results, the GMPEs of A14 and B14 were employed. Both were derived using strong motion datasets from Italy, Greece, and Turkey.

After validating the basic stochastic simulation parameters for the study area, three well-studied seismogenic faults (Psathopyrgos Fault (PT); Helike Fault (HF); Trizonia Fault (TF)) were selected to assess the worst-case earthquake scenarios (Table 2, Figure 1). Given the limited extent of the area under study and the consistency in the seismotectonic regime, we assumed that path and local site effects were properly constrained from the Aigion 1995 simulations. It was, therefore, only necessary to define the seismic source parameters and, particularly, the fault dimensions and geometry, as well as the Mmax. Fault characteristics were adopted from various studies (see Table 2 caption) and the Mmax was calculated following the approach described by Kourouklas et al. [103], which is a slightly different version of the method proposed by Pace et al. [104]. The adopted approach combined various scaling relationships between magnitude and rupture length with the maximum observed magnitude obtained from historical data. The finally adopted Mmax (Table 2) was weighted and acquired through the combination of all relative Mmax values, along with their standard deviations. In our assessment, the scaling laws from [105–108] were used.

**Table 1.** Modeling parameters used for the stochastic ground motion simulation of the Mw = 6.5 Aigion 1995 mainshock performed with the EXSIM code.


**Table 2.** Stochastic modeling source parameters for the three selected seismic scenarios. Fault orientation and dimensions were adopted from various sources [6,38,90]. Associated past strong events were adopted from Boiselet [38].


Furthermore, the definition of the hypocenter position and the slip distribution onto the fault plane is a crucial step in stochastic simulations. For the purposes of the present study, we divided each fault segment into subfaults using Equation (1), and the hypocenter position was placed randomly among them (minimum of 10 iterations). Different slip rupture patterns (minimum of 5 iterations) were examined, applying random slip weights to each subfault.

#### **3. Results**

#### *3.1. PSHA*

The spatial distribution of the PGA, computed via the logic tree approach (Figure 2) for a return period of 475 years, is illustrated in Figure 3a. The highest values were identified in close proximity to the coastline, near the towns of Nafpaktos and Aigion. The PGA decreased towards the north and reached its minimum of about 250 cm/s2 at the NE edge of the study area. A maximum of ~325 cm/s2 was obtained approximately 5 km SSE of Aigion. The difference between the highest and the lowest values was 75 cm/s2, which indicates that their variation was low to intermediate. The spatial distribution pattern of the PGA for a return period of 950 years was nearly identical (Figure 3b). A maximum of ~400 cm/s<sup>2</sup> was found SSE of Aigion, meaning that there was an increase of 75 cm/s2 compared to RP = 475 yr. The lowest value was about 300 cm/s2, increased by 50 cm/s2 from the respective value obtained for the return period of 475 years.

**Figure 3.** Spatial distribution of PGA values for a return period (RP) of (**a**) 475 yr and (**b**) 950 yr.

Regarding the seismic hazard intensity measure of the PGV, the spatial distribution was similar to the PGA for a return period of 475 years (Figure 4a). A maximum of ~18.5 cm/s was found toward the southern edge of Aigion, where the highest PGA was also determined. The lowest value persistently remained at the NE end of the study area, with an approximate PGV of 14 cm/s. In Nafpaktos, the PGVs were lower than in Aigion but higher than in Patras. The variation in the values was small, as the difference between

the highest and lowest PGV was only 4 cm/s. In Figure 4b, the spatial distribution of PGVs for a return period of 950 years is depicted. The highest PGV was about 24 cm/s, south of Aigion, and the lowest was 18 cm/s at the NE tip. The increase between PGVs for the two return periods was only 5 cm/s for the highest and 4 cm/s for the lowest values.

**Figure 4.** Spatial distribution of PGV values for (**a**) RP = 475 yr and (**b**) RP = 950 yr.

Hazard curves were also generated for the three main WGoC city and towns; namely, Patras, Nafpaktos, and Aigion (Figure 5a). PGA values were calculated for a range of exceedance probabilities in 50 years to observe the increase rate for the PGA as a function of the return period. The hazard curves were proximal to each other, indicating that the level of seismic hazard was similar among the three sites. Patras presented a slightly lower PGA than Nafpaktos and Aigion. It could not be determined whether Aigion or Nafpaktos exhibited higher seismic hazard, as the PGAs were almost identical between the two cities. It is worth noting that the value of 1000 cm/s<sup>2</sup> was exceeded only for very small probabilities of exceedance, which implies that such extreme ground motions are unlikely to be reached. The last output of the PSHA calculations was the computation of Spectral acceleration (Sa) values for periods ranging from 0.1 s to 2.0 s in the three aforementioned sites in order to construct uniform hazard spectra (Figure 5b). The Sa curves for Patras, Nafpaktos, and Aigion were very similar to each other throughout all periods. This was reasonable in light of the close geographical distance among them. As illustrated in both PGA and PGV maps, the city and the two towns were adjacent to regions of comparable maximum expected ground motions. The maximum of the three spectra was about 500 cm/s2 at a period of 0.25 s. The curves seemed to be nearly identical for periods above 0.9 s. Thus, at high periods, there was almost no difference between the three sites. Minor

deviations were detected in the period range [0.3, 0.5] s, for which Nafpaktos seemed to have slightly lower Sa values than Patras and Aigion. The elastic design spectra proposed by the Current National Building Code (EAK2003) [109] and Eurocode 8 (Ec8) [110] were also plotted to investigate their relation with our results. The seismicity was defined as high (type I) and the soil as bedrock (type A) to match the input data of the herein proposed model. The Ec8 spectrum overlay the spectra of the three urban areas for all periods, while that of EAK2003 covered the spectra for the period range [1.4, 2.0] s.

**Figure 5.** Hazard curves in terms of (**a**) PGA for Aigion, Nafpaktos, and Patras and (**b**) UHS for the aforementioned sites, alongside the Ec8 and EAK2003 response spectra. The red rectangle indicates the range of eigenperiods of most buildings in the study area.

#### *3.2. Seismic Scenarios*

#### 3.2.1. Validation of the Aigion 1995 Mainshock Simulation

In order to assess the ground motion variability caused by the Aigion 1995 mainshock, PGA values were computed on a grid enclosing the area under study, with nodes at a spacing of 0.03◦ in latitude and longitude. The spatial distribution of synthetic PGA values is depicted in Figure 6a, along with the surface projection of the activated fault [29] and the available recorded PGA. The highest values estimated onshore exceeded 500 cm/s2, whereas in the city of Aigion, where the highest observed PGA was reported, they were close to 300 cm/s2. The spatial distribution of simulated values was similar to those presented in other studies [88,111], despite using diverse input parameters. The validation of the final parameters used in the simulation is shown in Figure 6b, where simulated PGA values are plotted against those derived from the selected GMPEs and the recorded ones as a function of Joyner and Boore distances (Rjb). The functional form of the GMPEs used corresponded to a normal faulting style and type C (NEHRP) soil condition. As shown in Figure 6b, the PGAs obtained from EXSIM were in good agreement with the GMPE curves throughout the entire Rjb range, except for the points that lay very close to the surface projection of the fault plane (Rjb ≤ 2–3 km). This was, however, anticipated, given that GMPEs are generally not fully capable of reliably reproducing the ground motion in very short distances due to the limited availability of near-fault strong motion datasets, which affects their formulation. Moreover, directivity effects may have a strong influence on near-fault ground motion that is not fully captured by the GMPEs used. In this case, for example, the high PGA value recorded in Aigion (AIGA; Figure 6b) has been attributed to forward rupture directivity, in addition to local soil and topographic characteristics [29,111]. Nevertheless, the overall PGA variability lay inside the ±1σ range for both GMPEs in the entire Rjb range. Regarding the attenuation pattern, PGAs prescribed by the A14 and B14 models appeared to decay slightly faster than the simulated ones; however, taking into account the fact that the same finding also applied to the recorded ones (red triangles; Figure 6b), this highlights the importance of retrieving region-specific parameters that can be incorporated into GMPEs. The finally adopted synthetic PGA values were obtained based on the desire to achieve the right balance between keeping the lowest misfit among synthetic, GMPE-derived, and real values. Consequently, as shown in Figure 6b, our simulations were capable of reproducing the PGAs recorded from the two closest stations (AIGA, AMIA) more closely than the GMPEs, whereas at larger distances the general trend of the recorded PGAs (which presented a relatively high variability in certain cases) was captured by both GMPEs and synthetic values.

**Figure 6.** (**a**) Spatial distribution of simulated PGA values (cm/s2) calculated with the EXSIM code for the case of the Aigion 1995 mainshock. (**b**) Comparison between the recorded (red triangles), GMPE-derived (lines), and simulated (crosses) PGAs, plotted as a function of the Joyner and Boore distance (Rjb). The simulated values on the surface projection of the fault plane (Rjb = 0) were assigned a very small positive Rjb value (~0.01 km) to make them visible on the logarithmic axis.

3.2.2. Spatial Distribution of Simulated PGA for Selected Faults and Site-Specific Analysis

In the finite-fault stochastic scenarios of the present study, a worst-case seismic scenario was adopted in contrast to the PSHA, where peak ground motion parameters were obtained from a large number of possible earthquakes with a certain probability of exceedance. Figure 7 presents the spatial distribution of the simulated PGA values generated by the three selected faults (Table 2). The color scales retain a common, fixed range of PGA values to enable a comparison between each scenario. Simulated acceleration time series for Patras, Nafpaktos, and Aigion are also illustrated along with the resulting PGA values.

In all cases, the maximum PGAs, located at the surface projection of the faults or nearby, appeared to be relatively close, ranging from approximately 350 to 450 cm/s2. The HF may have had larger dimensions among the three cases, but PF produced the highest acceleration values onshore (Figure 7a), even though it was assigned a slightly lower Mmax (6.3) than HF (6.4; Table 2).

Similar to the PSHA framework, hazard response spectra for Aigion, Nafpaktos, and Patras for hypothetical future earthquakes originating from the selected faults (PF, HF, TF) were constructed for 5% damping (Figure 8). In this case, however, it became plausible to gain insight into which seismic sources are more capable of causing damage to each urban area. In all cases, the response spectra of simulated ground motions exhibited a spike behavior at short periods of [0.1–0.3] s, followed by a relatively sharp decline. As shown in Figure 8a, the town of Aigion is mostly threatened by HF and TF, whereas PF poses the greatest threat to Nafpaktos and Patras. The city of Patras (Figure 8b) is exposed to a lower level of seismic hazard compared to the other cases when taking into account the major active faults of the WGoC. Maximum Sa values exceeded 300 cm/s2 for a hypothetical rupture of PF, whereas HF produced the lowest expected values (approximately half of the corresponding ones for PF). Moreover, the town of Nafpaktos, located in between the aforementioned sites, is highly susceptible to high levels of seismic hazard due to the PF, with maximum expected Sa values approaching almost 900 cm/s2. TF and HF, which are located quite far from the town, do not pose significant hazards, given that the predicted maximum Sa values were within the ranges appointed by the National Building Code [109].

**Figure 7.** Spatial distribution of simulated PGA values (cm/s2) calculated with the EXSIM code for the cases of the (**a**) Psathopyrgos, (**b**) Helike, and (**c**) Trizonia Faults. Synthetic acceleration time series for Patras, Nafpaktos, and Aigion are also displayed along with the maximum obtained values.

**Figure 8.** Response spectra of the selected earthquake scenarios for (**a**) Aigion, (**b**) Patras, and (**c**) Nafpaktos corresponding to 5% damping. TF: Trizonia Fault, HF: Helike Fault, PF: Psathopyrgos Fault.

#### **4. Discussion**

The high seismicity of the study area, both onshore and offshore, is due to the high average extension rate of the tectonic rift [5]. It is expressed through both large earthquakes and seismic swarms [7,17,20,39,112–114], which confirm the great importance of seismic hazard assessment as a means to quantify anticipated levels of ground shaking.

When performing PSHA, there are two types of uncertainties, i.e., aleatory uncertainty and epistemic uncertainty [115]. The first one accounts for random variations in PGA, PGV and Sa values due to the implementation of a GMPE, whereas the second one accounts for the accuracy of the values [115]. Epistemic uncertainty is usually handled by the incorporation of a logic tree approach to account for the implementa-tion of more than one GMPE [115–117]. This is the procedure followed in the present study, given that by taking into account alternative models the uncertainties can be reduced [115].

The PGA results for a return period of 475 years, obtained through the Cornell– McGuire approach, can be directly compared with the reference value proposed by the current Greek Building Code [109]. EAK2003 divides Greece into three zones and defines the maximum expected ground acceleration for a return period of 475 years. The Western Gulf of Corinth is within zone II, with a reference value of 0.24 g (~240 cm/s2). However, in the herein proposed model, the entire study area exceeded 240 cm/s2. The computed PGA results in the northeastern part of the area were close to the EAK provisions, even though they still slightly exceeded them. The highest value was about 325 cm/s2, surpassing EAK2003 by 74%. The calculated PGA for a return period of 950 years can be utilized for the construction of buildings of greater importance, such as schools and medical centers. The spatial patterns of PGA and PGV for both return periods were similar; i.e., the highest values were found close to the Gulf, where seismicity is higher, whereas the lowest values were onshore, far from the coastline.

A comparison of the PGA results obtained to those of Banitsiotou et al. [49] and Tselentis et al. [54] was attempted, both of which are for the whole Greek territory. Banitsiotou et al.'s [49] computed PGA values were for certain cities, one of them being Patras, which was assigned a PGA of 0.26 g for a return period of 475 years. In this study, the corresponding PGA was ≈300 cm/s2, which is slightly higher. Tselentis et al. [54] computed PGA and PGV for the same return period for Greece. PGA values varied between 0.40 g and 0.50 g close to Patras and Nafpaktos and between 0.50 g and 0.60 g close to Aigion. Concerning PGV, Tselentis et al. [54] proposed values in the broad range of [20.0–105.0] cm/s. In the present study, PGA and PGV values were in the range of [0.25–0.33] g and [15.0–18.5] cm/s, respectively. The existing deviations can be attributed to the different seismic and geometry data used in each case.

Hazard curves were constructed for the most populated sites of the study area, the city of Patras and the towns of Nafpaktos and Aigion. All three are located close to the coastline and are characterized by high expected ground motion values. The distances between them are relatively small and they belong in the same or neighboring (according to the model) area source zones, parameterized with similar geometry and seismicity data. For this reason, their corresponding hazard curves were quite close to each other, whereas strong ground motions (>1 g) occurred only at low exceedance probabilities.

Uniform hazard spectra were obtained by adopting the GMPE from Danciu and Tselentis [76]. The final results were derived using the same logic tree but without the minor branches, as they were replaced with this GMPE. However, the adopted model was, again, parameterized to take into account all types of focal mechanisms for each zone. The values for the spectral acceleration were almost the same for all three cases, as was also observed in the hazard curves. Differences from the elastic design spectra of EAK2003 [109] and Ec8 [110] were evident. The fundamental periods of interest had ranges of [0.1, 0.5] s (Figure 5) for the majority of buildings in Greece. In this range of values, the spectra for Patras, Nafpaktos, and Aigion significantly surpassed EAK2003. However, they were fully overlaid by the design spectra from Ec8. Moreover, given that the probability of exceeding the respective Sa value was equal for all periods, and since all future earthquake events were taken into account, it is reasonable that higher values were computed. Hence, it is suggested that the Ec8 regulations should be respected for the majority of buildings in the study area.

In addition to the traditional PSHA approach, finite-fault stochastic simulations were performed by generating earthquake scenarios for a predicted maximum magnitude (Mmax), calculated by taking into account past seismicity and fault data. The three modeled active faults were the Psathopyrgos, Helike, and Trizonia Faults. Final results were obtained by testing both random earthquake nucleation points and co-seismic slip distribution on the surface. Prior to the computations, stochastic modeling parameters regarding the path and local site effects were calibrated by performing strong ground motion simulation of

the largest recent event; namely, the Aigion earthquake that struck the study area in 1995. Comparison with real strong motion data and GMPEs proved that this methodology is reliable for the reproduction of the ground motion, regardless of the uncertainties that are involved, which can be ascribed mainly to source complexity and local site condition variability. As a result, this methodology can be treated as a powerful tool for the simulation of the expected strong ground motion of future strong earthquakes, especially for cases where recordings are sparse or not available. It, therefore, provides a unique opportunity to comprehensively evaluate seismic hazard with the synergy between PSHA outcomes and specific scenarios.

Site-specific analysis at the three selected sites was performed through a comparison between response spectra for each seismic scenario. Sa values constitute a good indicator of seismic loading for a variety of structures, as they describe the absolute maximum response of a single degree-of-freedom oscillator to an enforced ground motion. The mean response spectra exceeded the current EAK2003 in some cases, as in Nafpaktos for the hypothetical rupture of the Psathopyrgos Fault. The city of Patras was revealed to be the site less exposed to seismic hazard among those examined for the selected scenarios. However, it should be noted that other faults outside the WGoC may be more dangerous for Patras, such as the Rion–Patras fault zone (Figure 1); the local NW–SE sinistral-slip fault, which was related to the 1993 earthquake [41]; or even the Andravida fault, further southwest, which has caused strong earthquakes in the past [118]. The latter two faults were not taken into account, as the present study assessed seismic hazard in the WGoC area.

Overall, the PGA variability obtained from the specific fault scenarios lay inside the PGA values calculated through the PSHA. Figure 3a,b indicate that the western part of the study area exhibited a similar level of seismic hazard under the PSHA evaluation. The expected peak ground motions in this part, however, were considerably affected by the nearby seismic sources to the west (the Ionian Islands and the westernmost end of the Hellenic Arc), which are among the most active and productive in terms of seismic potential [119]. PSHA provided valuable information reflecting the combined effects of all potential seismic sources, making it possible to assess seismic hazard in terms of statistical likelihood of occurrence. On the other hand, the employment of seismic scenarios made it possible to acquire more realistic results concerning specific strong events through the definition of local site effects and path properties but without taking into account the time frame of occurrence.

#### **5. Conclusions**

A holistic seismic hazard assessment, including both probabilistic and scenario-based methods, was conducted for the highly active Western Gulf of Corinth area in Central Greece. Initially, seismic hazard was evaluated using a probabilistic approach, taking into account all the Mw ≥ 4.0 earthquakes recorded during 1900–2019, in order to determine values for PGA, PGV, and Sa in a dense grid. For this purpose, the Cornell–McGuire method was utilized. Aiming to incorporate a multitude of seismotectonic models in order to qualitatively cope with epistemic uncertainties, three models—namely, ESHM13, ESHM20, and VAM16—were considered, and GMPEs were introduced for each zone. The PGA and PGV results were obtained through an equal-weight logic tree approach in which each major branch was a seismotectonic model and each minor branch was a modified GMPE. The logic tree procedure incorporated with the five GMPEs that were developed using Greek data qualitatively reduced epistemic uncertainties, which simple PSHA approaches may carry. In addition, PGA values were computed for the most populated urban areas (Patras, Nafpaktos, and Aigion) for various exceedance probabilities in 50 years.

The results obtained herein highlight the great significance of seismic hazard assessment using a combined approach that takes into account not only the evaluation of ground parameters in terms of probabilities of occurrence in a given time frame but also the anticipated effects of deterministic worst-case scenarios. Future work could include disaggregation of PSHA results so that the parameter pair of magnitude and epicentral

distance that contributes most to seismic hazard can be identified. In this way, scenarios could be generated based on this result. Furthermore, a study of the Peak Ground Rotational Acceleration (PGRA) and Velocity (PGRV) values could be undertaken, as it has been proven that their results provide aid for engineering purposes [58,77,120]. In addition, the incorporation of a comprehensive microzonation scheme could provide valuable insight into the impact of future strong earthquakes on the major cities of the study area by identifying the possible amplification trends with respect to the structural response. Lastly, the results of the present seismic hazard study can be exploited in the future to assess seismic risk at urban centers in the WGoC area after incorporating structural vulnerability data.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www. mdpi.com/article/10.3390/app122111152/s1, Figure S1. The three zonation models that were implemented for the WGoC, (a) the ESHM13 model [62,63], (b) the ESHM20 model [64] and (c) the model VAM16 from Vamvakaris et al. [65] proposed for the Greek territory; Table S1: The seismicity and geometry data that was computed for each Area Source Zone (ASZ).

**Author Contributions:** Conceptualization, G.K.; Data curation, G.K. and A.Z.; Formal analysis, G.K. and A.Z.; Funding acquisition, G.K.; Investigation, G.K., A.Z., P.B. and V.K.; Methodology, G.K.; Resources, N.V.; Software, A.Z. and P.B.; Supervision, G.K. and N.V.; Validation, G.K. and A.Z.; Visualization, P.B. and V.K.; Writing—original draft, G.K., A.Z., P.B. and V.K.; Writing—review & editing, G.K., A.Z., V.K. and N.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available datasets were analyzed in this study. ESHM13,20 ASZ can be found here: http://hazard.efehr.org/en/home/, accessed on 30 October 2022. In accordance with the earthquake catalogue, the earthquakes for the period 1900-2009 were adopted from Makropoulos et al. [25]. For the temporal expansion (2010-1019) of the catalogue, the reviewed events from International Seismological Centre (http://www.isc.ac.uk/iscbulletin/search/bulletin/, accessed on 30 October 2022) were adopted.

**Acknowledgments:** We express our gratitude to the three anonymous reviewers for their constructive comments that helped improve the manuscript, as well as to the editor and the guest editors for providing us with the opportunity to publish this study. We would like to thank the personnel of the following institutions: (a) the CRL team (CL network, data hosted at RESIF, https://doi.org/ 10.15778/RESIF.CL, accessed on 30 October 2022); (b) the National and Kapodistrian University of Athens (NKUA) (HA network, data hosted at NOA, https://doi.org/10.7914/SN/HA, accessed on 30 October 2022); (c) the University of Patras (HP network, data hosted at NOA, https://doi. org/10.7914/SN/HP, accessed on 30 October 2022), which operates 11 stations jointly with Charles University, Prague; (d) the National Observatory of Athens (NOA) (HL network, data hosted at NOA, https://doi.org/10.7914/SN/HL, accessed on 02 November 2022); and (e) the Institute of Engineering Seismology and Earthquake Engineering (ITSAK) (HI network, data hosted at NOA, https://doi. org/10.7914/SN/HI, accessed on 30 October 2022). Stations operated by the last four institutes are also part of the Hellenic Unified Seismological Network (HUSN, https://www.gein.noa.gr/ en/networks-equipment/hellenic-unified-seismic-network-h-u-s-n/, accessed on 30 October 2022). Figure 1 was drawn using the Generic Mapping Tools (GMT) software [121]. Grapher version 15 (http://www.GoldenSoftware.com, accessed on 30 October 2022) was used for some of the figures.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

