Long-Term Recurrence Pattern and Stress Transfer along the Kefalonia Transform Fault Zone (KTFZ), Greece: Implications in Seismic Hazard Evaluation

Abstract

An effort is exerted to investigate the recurrence pattern of large earthquakes (M_w ≥ 6.0) in the Kefalonia Transform Fault Zone (KTFZ), Greece, by considering the incorporation of the 74-year (1948–2022) evolving stress field. Four earthquake occurrence models—the Poisson, Poisson with the incorporation of the static stress changes (Poisson + ΔCFF), Brownian passage time (BPT) and Brownian passage time with the incorporation of the static stress changes (BPT + ΔCFF)—have been applied to estimate the occurrence probabilities of nearly characteristic earthquakes for the seven fault segments of the study area. The mean recurrence time, T_r, is estimated using the physics-based seismic moment rate conservation method. The results show large variability depending upon fault parameters. Incorporating the state of stress into T_r results in both advanced and delayed recurrence patterns. The occurrence probability estimates for the next 10, 20 and 30 years indicate that the fault segment most likely to be ruptured is the Paliki North fault segment in all models. Overall, the occurrence probabilities, combined with the state of stress along the fault segments, emphasize the high seismic moment rate of the study area. The application of time-dependent models (BPT, BPT + ΔCFF) resulted in significant increases or decreases in the associated seismic hazard.

1. Introduction

Large earthquake (e.g., M ≥ 6.0) recurrence behavior on specific faults or fault segments is the primary input for developing long-term earthquake rupture forecast (ERF) models [1]. These models combine a series of parameters (maximum observed magnitude, fault dimensions, long-term slip rates, available recurrence times) for the estimation of the likelihood of the occurrence of nearly characteristic magnitude earthquakes [2] in a specific time span. This occurrence is neither purely periodic, as described by the elastic rebound theory [3], nor completely random in time, but exhibits short- and long-term clustering (e.g., [4,5]). This complex temporal behavior could be explained by physical processes, such as fault heterogeneity [6] and the occurrence of lower-magnitude earthquakes that take up a sufficient part of the accumulated strain. One important factor influencing the recurrence time on certain fault segments is the interaction among them due to the permanent and transient static stress changes caused by the coseismic slip of large earthquakes [7,8,9]. Specifically, static stress changes among the interacting causative and receiving fault segments are capable of affecting the mean recurrence time, T_r, by moving a given segment towards (promoting) or away (delaying) from failure. This implies that the recurrence time should be addressed using probabilistic models rather than deterministic ones (e.g., [10]).

The time-independent Poisson model and a series of renewal models, which consider the elapsed time since the last large earthquake on a certain fault, have been proposed to investigate the recurrence pattern. The key input parameters are the mean recurrence time, T_r, and its aperiodicity, α. Aperiodicity can be considered analogous to the coefficient of variation of the normal distribution and represents the level of randomness, with values ranging between 0 to 1 (0 ≤ α ≤ 1), addressing the physical meaning of large earthquake recurrence [11]. A precise and robust estimation of T_r requires the inclusion of as many large earthquakes as possible, including both historical and instrumental data associated with individual fault segments, along with the selection of the appropriate statistical model. However, observational data of T_r per fault segment are often limited, with only a few cases having 3 to 10 observations [12,13] due to the long duration of the stress rebuild and the shortness of available large earthquake records. Additional data and methods, such as paleoseismic records (e.g., [14]), slip-per-event constraints (e.g., [15]) and Bayesian methods (e.g., [16]), have been proposed to refine the parameter space of the recurrence models.

An alternative method for addressing the difficulties and limitations associated with the absence of multiple large earthquakes associated with specific fault segments is the estimation of T_r through the application of the seismic moment rate conservation method [17]. T_r could be defined as the ratio of the maximum expected seismic moment corresponding to a large earthquake with the maximum observed magnitude, M_{max_obs}, occurring on a given fault segment to the seismic moment that can be released by the respective fault segment due to the tectonic loading, assuming a nearly characteristic earthquake model. This method can provide precise T_r estimates because the overall rate and size distribution of earthquakes should reflect the tectonic loading in the brittle part of the crust, and this function is typically constrained by invoking conservation of seismic moment rate. The method is widely used for relevant applications worldwide, including the Corinth Gulf Fault Zone in Greece [18], Calabria [19] and Central Apennines [20,21] Fault Zones in Italy and the Sea of Marmara region in Turkey [22,23]. It is employed not only for the single estimation of the mean recurrence time but also for combining these estimations along with the effect of the stress interactions between adjacent segments, considering the advanced or delayed time shifts due to the static stress changes.

Kefalonia Transform Fault Zone (KTFZ; red rectangle in Figure 1) is recognized as the active boundary connecting the Hellenic Arc, formed by the subduction of the oceanic lithosphere of the Eastern Mediterranean under the continental Aegean microplate (Figure 1; [24]), to the south and the continental collision of the Adriatic microplate and Eurasia to the north parallel both onshore and offshore to the western coastal areas of Greece and Albania [25]. KTFZ exhibits the highest crustal deformation rates within the Aegean area with slip rates up to 25–30 mm/yr [26,27,28], resulting in the frequent occurrence of earthquakes with M_w ≥ 6.0 during both historical and instrumental periods [29] (Figure 2). The temporal occurrence pattern of these large earthquakes can be explained by triggering due to the stress transfer among the adjacent fault segments [30]. These facts provoke the detailed study of the T_r of earthquakes with M_w ≥ 6.0 along the KTFZ fault segments, embodied with the effects of the evolutionary stress field caused by both long-term tectonic loading and coseismic displacements. The main objectives of this study are the determination of the T_r of large (M_w ≥ 6.0) earthquakes associated with the major fault segments of the study area and the calculation of the evolutionary stress field. The results of these approaches could form the basis for the statistical analysis by both the Poisson model and the Brownian passage time (BPT) renewal model, allowing the estimation of the occurrence probabilities of nearly characteristic earthquakes for the fault segments of the study area.

2. Fault Segmentation Model of Kefalonia Transform Fault Zone and Slip Rate Constraints

Kefalonia Transform Fault Zone is characterized by right-lateral strike-slip motion with a minor thrust component [31,32,33,34] and is the most active seismic zone in the Aegean area [35]. It can be distinguished into two main fault branches, namely the Lefkada branch to the north covering its northern part, striking NNE–SSW and the Kefalonia branch to the south with a slightly different NE–SW strike (Figure 2). The mean thickness of the seismogenic layer along the KTFZ is equal to 12 km, with earthquake focal depths ranging between 3 and 15 km [36]. The large earthquakes that occurred in the study area since 2003 (Figure 2; 2003 with M_w = 6.2 in the northern part of Lefkada; 2014 with M_w = 6.1 and 6.0 in the Paliki Peninsula in Kefalonia; 2015 with M_w = 6.5 in the southern part of Lefkada), along with the one that occurred offshore southeast of the Island of Kefalonia in 1983 with M_w = 7.0, have attracted research interest ([37,38,39,40,41,42,43], among others), and thus a detailed segmentation model for KTFZ is available. Each fault segment considered in this study is modeled using the standard source characterization approach, promoting specific fault segments with fixed endpoints [44,45]. The workflow of the proposed fault segment network model is primarily based on the elaboration of the precise relocated aftershocks’ spatial 3D distribution of large earthquakes, along with their available focal mechanisms (GCMT solutions, [33]).

Figure 2. Historical and instrumental seismicity of central Ionian Islands. Magenta stars depict all known M_w ≥ 6.0 historical earthquakes until 1948. Epicentral distribution of after 1948 M_w ≥ 6.0 earthquakes is denoted with the yellow stars. Small, moderate and large white circles depict the 3.0 ≤ M_w < 4.0, 4.0 ≤ M_w < 5.0 and 5.0≤ M_w < 6.0 earthquakes that occurred in the study area since 2003. Fault plane solutions are shown as equal area lower hemisphere projections with the compressional quadrants depicted in red and green colors for the strike-slip and thrust mechanisms, respectively. Red solid lines depict the major fault segments, and the red arrows represent the right-lateral strike-slip motion (S1: Lefkada North fault segment; S2: Lefkada South fault segment; S3: Paliki North fault segment; S4: Paliki South fault segment; S5: offshore Kefalonia fault segment; S6: Argostoli thrust fault segment; S7: Ainos thrust fault segment).

Figure 2. Historical and instrumental seismicity of central Ionian Islands. Magenta stars depict all known M_w ≥ 6.0 historical earthquakes until 1948. Epicentral distribution of after 1948 M_w ≥ 6.0 earthquakes is denoted with the yellow stars. Small, moderate and large white circles depict the 3.0 ≤ M_w < 4.0, 4.0 ≤ M_w < 5.0 and 5.0≤ M_w < 6.0 earthquakes that occurred in the study area since 2003. Fault plane solutions are shown as equal area lower hemisphere projections with the compressional quadrants depicted in red and green colors for the strike-slip and thrust mechanisms, respectively. Red solid lines depict the major fault segments, and the red arrows represent the right-lateral strike-slip motion (S1: Lefkada North fault segment; S2: Lefkada South fault segment; S3: Paliki North fault segment; S4: Paliki South fault segment; S5: offshore Kefalonia fault segment; S6: Argostoli thrust fault segment; S7: Ainos thrust fault segment).

In detail, the KTFZ comprises five major dextral fault segments with strikes ranging from 12° to 40°, lengths of 12–40 km and typical rake values for right-lateral strike-slip faulting (Figure 2 and Table 1). Starting from the western coasts of Lefkada Island, the Lefkada North and Lefkada South segments (S1 and S2 in Figure 2, respectively) are defined. The Lefkada North segment (S1) is defined by considering the precise relocated aftershock seismicity of the 2003 M_w = 6.2 main shock, as analyzed by [37,46] and the respective GCMT solution (Figure 1 and Table 2; https://www.globalcmt.org). The Lefkada North fault segment strikes at 18° (φ = 18°), dips east-southeast with an angle of 60° (Δ = 60°) and has a rake angle equal to −175° (λ = 175°). Its dimensions equal 16 and 10 km in length and width, respectively [37,46]. The Lefkada South segment (S2) is considered the causative fault of the 2015 M_w = 6.5 main shock [40], taking into account again the relocated aftershock seismicity and the GCMT solution (Figure 2) as well. The S2 fault segment also strikes NNE-SSW NEN–SWS (φ = 22°) and dips to the east-southeast at a slightly larger angle equal to 64° and rake angle equal to 179°. Its length is equal to 20 km (L = 20 km), and its width equals 12 km (w = 12 km) [40]. The thickness of the seismogenic layer in this part of the fault zone is equal to 9 km (ranging from 5 to 14 km), according to [37,40,46].

In the western part of Kefalonia Island, along the Paliki Peninsula, the Paliki North and South fault segments (S3 and S4, respectively, in Figure 2) are extended. S3 and S4 have almost the same dimensions and are associated with the 2014 Kefalonia doublet (the 26 January and the 3 February 2014 earthquakes with M_w = 6.1 and 6.0, respectively). The geometric and kinematic parameters of these two fault segments are defined by combining the relocated aftershock seismicity of the doublet [38] and the available GCMT solutions (Figure 2 and Table 2). Specifically, their length and width are both 12 km and 10 km, respectively [38]. The S3 fault segment strikes NNE–SSW (φ = 20°) and dips to E-SE at a steep angle of 65° (Table 1), whereas the S4 fault segment strikes NS (φ = 20°) with a dip angle of 45°. The fifth and largest fault segment of KTFZ, namely the offshore Kefalonia fault segment (S5), is associated with the 1983 M_w = 7.0 earthquake (Figure 2 and Table 1). The relocated aftershocks indicate a NE–SW direction, in agreement with one of the nodal planes of the GCMT solution of the main shock [31]. The S5 fault segment has a length equal to 33 km and width equal to 20 km, NE–SW strike (φ = 40°), dip angle equal to 45° (Δ = 45°) and slip angle equal to 168° (λ = 168°) [31]. The seismogenic layer in Kefalonia Island is slightly thicker than in Lefkada and equals 15 km (from 3 km to 18 km) [38].

Onshore Kefalonia thrust faulting dominates, as the available geological studies reveal [47,48,49], among others. Geological mapping identifies two major thrust fault segments in the southeastern part of Kefalonia Island, namely the Argostoli (S6) and Ainos (S7) segments. The Argostoli fault segment has a length of 15 km, a width of 15 km (w = 15 km), strikes at 299°, dips at a low angle of 30° and exhibits pure dip-slip thrust faulting (λ = 90°) [47]. Further east, the Ainos fault segment is located along the southwestern slopes of Mount Ainos (S7; Figure 2 and Table 1). S7 has a strike of 300°, a dip of 30°, a rake angle of 100°, a length of 35 km [48,49] and a width of 24 km, as calculated from its dip angle and the thickness of the seismogenic layer (12 km).

Slip rates gradually increase from north to south, ranging from 10 mm/yr for Lefkada Island to 25 mm/yr for the southwestern margin of the KTFZ [27,28]. Additionally, [50] suggested slip rate values equal to 19.5 mm/yr along the Paliki Peninsula in the western part of Kefalonia Island. More recently, [51] confirmed this gradual increase in deformation rates by proposing slip rate values of 9.4 mm/yr in the northern part of KTFZ (Lefkada Island) and 14.3 mm/yr in the southern part of the Paliki Peninsula. They also concluded that the KTFZ is the only area within the Aegean region where the tectonic deformation is fully coupled. This latter conclusion agrees with the results of [26,52], who compared the seismic moment rates with the long-term geodetic moment rates and showed that the total amount of the geodetic deformation is translated into seismic activity. In our study, we assigned slip rates of 10 mm/yr to both Lefkada fault segments (Lefkada North and Lefkada South) and slip rates of 19.5 mm/yr to the Kefalonia Island fault segments (Paliki North and South and offshore Kefalonia fault segments), as proposed by [27,50,51]. Slip rate values for the Argostoli and Ainos thrust fault segments are considered equal to 4.9 mm/yr [53] (Table 1).

Table 1. Geometric and kinematic parameters of the fault network model of the Kefalonia Transform Fault Zone considered in the present study.

Fault Segment Name	Fault ID	Upper Left Edge of the Fault Segment		Strike (deg)	Dip (deg)	Rake (deg)	L (Km)	W (Km)	Slip Rate (mm/yr)	Stressing Rate $\dot{\tau }$ (bar/yr)	Ref.
		Lat	Lon
Lefkada North	S1	38.69	20.56	18	60	−175	16	10	10 ± 0.5	0.8459	[27,37,46]
Lefkada South	S2	38.55	20.49	22	64	179	20	12	10 ± 0.5	0.6907	[27,40]
Paliki North	S3	38.26	20.0	20	65	177	12	10	19.5 ± 0.5	1.9046	[38,50]
Paliki South	S4	38.15	20.35	12	45	156	12	10	19.5 ± 0.5	1.9046	[38,50]
Offshore Kefalonia	S5	37.91	20.08	40	45	168	33	20	19.5 ± 0.5	0.8518	[31,50]
Argostoli	S6	38.08	20.55	299	30	90	15	15	4.9 ± 1.0	0.3495	[47,53]
Ainos	S7	38.03	20.79	300	30	100	35	24	4.9 ± 1.0	0.1809	[48,53]

Table 1. Geometric and kinematic parameters of the fault network model of the Kefalonia Transform Fault Zone considered in the present study.

Fault Segment Name	Fault ID	Upper Left Edge of the Fault Segment		Strike (deg)	Dip (deg)	Rake (deg)	L (Km)	W (Km)	Slip Rate (mm/yr)	Stressing Rate $\dot{\tau }$ (bar/yr)	Ref.
		Lat	Lon
Lefkada North	S1	38.69	20.56	18	60	−175	16	10	10 ± 0.5	0.8459	[27,37,46]
Lefkada South	S2	38.55	20.49	22	64	179	20	12	10 ± 0.5	0.6907	[27,40]
Paliki North	S3	38.26	20.0	20	65	177	12	10	19.5 ± 0.5	1.9046	[38,50]
Paliki South	S4	38.15	20.35	12	45	156	12	10	19.5 ± 0.5	1.9046	[38,50]
Offshore Kefalonia	S5	37.91	20.08	40	45	168	33	20	19.5 ± 0.5	0.8518	[31,50]
Argostoli	S6	38.08	20.55	299	30	90	15	15	4.9 ± 1.0	0.3495	[47,53]
Ainos	S7	38.03	20.79	300	30	100	35	24	4.9 ± 1.0	0.1809	[48,53]

Once the geometric and kinematic parameters of the seven fault segments are defined, the stressing rate, $\dot{\tau }$, due to the long-term tectonic loading, is calculated. The stressing rate is assumed to be unchanged by the stress perturbations, being identical for each fault segment. It is calculated based on the fault segment dimensions (L, w) and long-term slip rate, V, using the equation proposed by [54]:

$$\dot{\tau }=\frac{32\mu V}{{\pi }^2\sqrt{Lw}}$$

(1)

The highest stressing rates are observed on the Paliki North and South fault segments, with values equal to $\dot{\tau }$ = 1.9046 bar/yr. The lowest stressing rate values are found for the Argostoli and Ainos thrust faults, with values equal to $\dot{\tau }$ = 0.3495 bar/yr and $\dot{\tau }$ = 0.1809 bar/yr, respectively, due to their significantly lower slip rate values (Table 1). An intermediate $\dot{\tau }$ value is calculated for the largest strike-slip fault segment of KTFZ, the offshore Kefalonia segment, which is equal to $\dot{\tau }$ = 0.8518 bar/yr. Intermediate stressing rate values characterize the Lefkada Island fault segments, with values equal to $\dot{\tau }$ = 0.8459 bar/yr and $\dot{\tau }$ = 0.6907 bar/yr for the Lefkada North and South fault segments, respectively. The higher value of the North Lefkada fault segment is due to its smaller fault area despite both segments being assigned with the same slip rate.

3. Large Earthquake Occurrence (M_w ≥ 6.0) along Kefalonia Transform Fault Zone

Large earthquake focal parameters are equally important to the detailed fault segmentation model for both the mean recurrence time, T_r, and the evolutionary stress field estimations. In both approaches, it is necessary to define the association between large earthquakes and specific fault segments. Given that the seismic moment rate conservation method refers to characteristic or nearly characteristic earthquakes that correspond to the rupture of the entire fault segment, the M_w ≥ 6.0 threshold is considered. This threshold is considered because M_w ≥ 6.0 earthquakes are capable of rupturing the entire seismogenic thickness of the upper crust [6,55]. In this study, the focal parameters of the M_w ≥ 6.0 earthquakes since 1948 are compiled from the parametric earthquake catalog of [56] and the regional earthquake catalog of [57] (http://geophysics.geo.auth.gr/ss/catalogs_en.html; last accessed on 10 June 2023). This dataset comprises 11 earthquakes associated with all major fault segments described in the previous section (Table 2 and Figure 2).

Both catalogs are homogeneous with respect to the earthquake magnitude scale; that is, all magnitudes are equivalent moment magnitudes (either directly computed or appropriately converted from other magnitude scales [58]). Location uncertainties concerning the pre-1980 earthquakes included in the parametric catalog of [56] are of the order of 10–20 km, whereas those of the earthquakes that occurred after 1980 are even smaller, equal to a few km [29]. It is also worth noting that the focal parameters of all the M ≥ 6.0 earthquakes considered in the present study were relocated after [59] with epicentral uncertainty less than 3 km. Additionally, supporting information about early instrumental earthquakes is compiled from the parametric catalog of [60].

The association between the M_w ≥ 6.0 earthquakes after 1983 and their causative faults is primarily determined by the combination of their focal mechanisms from the GCMT database and the precise 3D spatial distribution of their aftershock sequences, as already discussed in Section 2. Likewise, the 1972 earthquake, for which a fault plane solution is also available from [33], is associated with the Paliki North fault segment by considering both its epicentral location and the good agreement between its focal mechanism and the geometric and kinematic parameters of the respective fault segment (Figure 2; Table 1 and Table 2). The pre-1972 M_w ≥ 6.0 earthquakes are associated with the respective ruptured fault segments by considering all the available macroseismic data, including maximum macroseismic intensities, I_max, the spatial distribution of the macroseismic intensity assignments and descriptions of the damages caused by these earthquakes, as reported by [29], aiming to identify evidence of the localized coseismic deformation. These criteria were already applied by [42] for the 22 April 1948 M_w = 6.5 and the 30 June 1948 M_w = 6.4 earthquakes, suggesting that these earthquakes are associated with the Lefkada South and Lefkada North fault segments, respectively. Following this suggestion, the focal mechanisms and depths of these earthquakes are similar to the ones of the instrumental earthquakes (the fault plane solutions of the 2015 and 2003 earthquakes for the case of 22 April and 30 June 1948, respectively, denoted with asterisks in Table 2) associated with the respective fault segments since the seismogenic layer of the study area is well constrained by [59].

For the association of the 9 August 1953 M_w = 6.4 and the 12 August 1953 M_w = 7.2 earthquakes that occurred in Kefalonia Island with the Argostoli and Ainos thrust fault segments, respectively, the same criteria are considered. For this reason, all the available macroseismic information provided by [29] and [61] are compiled. According to the macroseismic data, Argostoli town was severely affected by the 9 August 1953 M_w = 6.4 earthquake, with most of the buildings damaged, whereas less damage was observed in other towns and villages of Kefalonia Island. Taking this information into account, the 9 August earthquake is associated with the Argostoli thrust fault. The 12 August 1953 M_w = 7.2 totally devasted the southeastern part of Kefalonia Island, whereas the islands of Ithaki and Zakynthos were severely damaged as well [29,61]. Further, [62] and the references therein reported anomalies in the spatial distribution of the macroseismic intensity values, with the largest ones (IX–X) to be observed in the southeastern part of Kefalonia and the Ithaki Islands, while in the northern part of Kefalonia, macroseismic intensities are significantly lower (of the order of VI), indicating the localization of the coseismic deformation in the footwall of the mapped Ainos thrust fault segment. These facts, along with the available focal mechanism of the 12 August earthquake [63] and the length of the Ainos fault segment, which is capable of hosting an M_w = 7.2 earthquake, led to its association with the Ainos fault segment.

Table 2. Focal parameters of the with M_w ≥ 6.0 occurred within KTFZ since 1948. Focal mechanism parameters are taken from GCMT database, and related references are given in the last column of table unless denoted with asterisks.

Date	Time	Epicenter		Depth (km)	Mw	Mechanism (deg)			Ref.
		Lat. (°N)	Lon. (°E)			Strike	Dip	Rake
22 April 1948	10:42:45	38.620	20.570	13.9 *	6.5	22 *	64 *	179 *	[56]
30 June 1948	12:21:13	38.800	20.600	11.0 *	6.4	18 *	60 *	−175 *	[56]
9 August 1953	07:41:07	38.430	20.500	11.0	6.4	299	30	90	[47,56]
12 August 1953	19:23:52	38.100	20.350	11.0	7.2	300	30	100	[56,60,63]
17 September 1972	14:07:15	38.300	20.300	8.0	6.3	46	66	−174	[33,56]
17 January 1983	12:41:30	38.100	20.200	9.0	7.0	39	45	175	[56]
23 March 1983	23:15:05	38.200	20.300	7.0	6.2	31	69	174	[56]
14 August 2003	05:14:55	38.815	20.606	11.0	6.2	18	60	−175	[56]
26 January 2014	13:55:41	38.199	20.434	13.5	6.1	20	65	177	[57]
3 February 2014	03:08:44	38.269	20.410	9.4	6.0	12	45	156	[57]
17 November 2015	07:10:07	38.677	20.577	13.9	6.5	22	64	179	[57]

Table 2. Focal parameters of the with M_w ≥ 6.0 occurred within KTFZ since 1948. Focal mechanism parameters are taken from GCMT database, and related references are given in the last column of table unless denoted with asterisks.

Date	Time	Epicenter		Depth (km)	Mw	Mechanism (deg)			Ref.
		Lat. (°N)	Lon. (°E)			Strike	Dip	Rake
22 April 1948	10:42:45	38.620	20.570	13.9 *	6.5	22 *	64 *	179 *	[56]
30 June 1948	12:21:13	38.800	20.600	11.0 *	6.4	18 *	60 *	−175 *	[56]
9 August 1953	07:41:07	38.430	20.500	11.0	6.4	299	30	90	[47,56]
12 August 1953	19:23:52	38.100	20.350	11.0	7.2	300	30	100	[56,60,63]
17 September 1972	14:07:15	38.300	20.300	8.0	6.3	46	66	−174	[33,56]
17 January 1983	12:41:30	38.100	20.200	9.0	7.0	39	45	175	[56]
23 March 1983	23:15:05	38.200	20.300	7.0	6.2	31	69	174	[56]
14 August 2003	05:14:55	38.815	20.606	11.0	6.2	18	60	−175	[56]
26 January 2014	13:55:41	38.199	20.434	13.5	6.1	20	65	177	[57]
3 February 2014	03:08:44	38.269	20.410	9.4	6.0	12	45	156	[57]
17 November 2015	07:10:07	38.677	20.577	13.9	6.5	22	64	179	[57]

4. Methods

4.1. Mean Recurrence Time Estimation

We estimated the mean recurrence time, T_r, by applying the seismic moment rate conservation technique [17] in the absence of a sufficient number of recurrence intervals for each of the 5 strike-slip fault segments of KTZF and the two thrust fault segments of Kefalonia Island. This technique assumes that the total amount of accumulated seismic moment on the fault can be released by a large earthquake with nearly characteristic magnitude. It considers the maximum observed magnitude (M_{max_obs}) and its corresponding uncertainty (ΔM) and the maximum seismic moment corresponding to accumulated strain on this specific fault segment due to the tectonic loading. The mean recurrence time is calculated as the ratio of the seismic moment rate that can be released by the maximum magnitude earthquake divided by the seismic moment due to strain accumulation on the fault:

$$T_r=\frac{M_o{}_{max}}{\mu LwV}$$

(2)

where $M_o{}_{max}$ is the maximum seismic moment that can be released by a large earthquake with a magnitude within the range M_{max_obs} ± ΔM, μ is the shear modulus, whose typical value for faults in the Earth’s crust is equal to 3.3 × 10⁵ bar (μ = 3.3 × 10⁵ bar), L and w the length and width of the fault segment (in km), respectively, and V the long-term slip rate. Maximum seismic moment is calculated with the equation of [64]:

$$M_o{}_{max}={10}^{\left(1.5M_{max\_obs}+9.1\right)}$$

(3)

in SI units. Uncertainties in the T_r computation are caused by the estimates of the maximum observed magnitude and the values of the long-term slip rate. To account for the effect of these uncertainties, two approaches for computing its corresponding confidence interval are followed.

The first approach is the Monte Carlo simulations, in which T_r is iteratively calculated, considering the uncertainties of maximum observed magnitude and the long-term slip rate ([19,23], among others). The procedure is implemented by the sequential calculation of T_r by generating N random values corresponding to the ranges M_{max_obs} ± ΔM and V ± σ_{_V} for the maximum observed magnitude and the long-term slip rate, respectively, assuming a certain statistical distribution. By applying this iterative process, a set of N values of T_r is obtained, whose median value is the final estimated T_r. The corresponding 67% confidence interval (corresponding to ±1σ) of this estimate can be easily defined as the range between the 16.5th and 83.4th percentiles of the simulated values, respectively.

The second approach [65] estimates T_r by direct application of Equation (2), and its corresponding uncertainty is calculated by applying an error propagation technique [66]. According to this technique, the general relationship describing the error propagation of a function can be described by a Taylor series:

$${\sigma }_{T_r}^2\approx {\left(\frac{{10}^{\left(1.5M_{\max 0bs}+9.1\right)}1.5\log \left(10\right)}{\mu LwV}\right)}^{2}d{M}_{max}^2+{\left(\frac{{10}^{\left(1.5M_{\max obs}+9.1\right)}}{\mu Lw{V}^2}\right)}^{2}d{V}^2$$

(4)

where $d{M}_{\max \_obs}$ and $dV$ represent the uncertainties of ΔM and σ_{_V} of the maximum observed magnitude and slip rate, respectively. Knowing the variability of mean recurrence time, ${\sigma }_{T_r}^2$, the standard deviation, σ, can be easily calculated. The advantage of this approach is that having estimated T_r and σ, one can compute the corresponding aperiodicity values, α, as the coefficient of variation, C_v, of T_r, which is the second appropriate parameter of renewal model applications.

We applied the Monte Carlo simulation approach by creating N = 1000 random samples of maximum observed magnitude and long-term slip rate between the ranges M_{max_obs}± ΔM and V ± σ_{_V}, respectively, assuming the uniform distribution as [19,23] proposed. This assumption for both maximum observed magnitude and long-term slip rate is made because M_{max_obs} mainly corresponds to historical large earthquakes that occurred in the study area, and therefore, their uncertainties are rather large. Further, there is no evidence that the slip rate follows a certain known statistical distribution to qualify it over the uniform one. In the case where the maximum observed magnitude refers to an earthquake that occurred after 1970, the magnitude uncertainty is considered equal to 0.2 (ΔM = ±0.2); otherwise, ΔM is considered equal to 0.3 (ΔM = ±0.3) [56].

4.2. Static Stress Change Calculation

The seismic moment rate conservation method assumes each fault segment is an isolated structure that does not interact with its neighboring segments and is not affected by the permanent changes in the stress field. It is well accepted that large earthquake occurrence is controlled by fault interactions through redistribution of stress, with a particularly clear example being the Kefalonia Transform Fault Zone [30]. For this reason, the cumulative changes in stress that are assumed to arise from tectonic loading on the major regional fault segments and coseismic displacements associated with large (M_w > 6.0) earthquakes are considered. This approach is widely applied in the context of identifying possible future rupture zones worldwide [67,68,69,70].

Interseismic stress accumulation is modeled by introducing “virtual negative displacements” across these faults, using the best available information on long-term slip rates [71]. Hence, tectonically induced stress builds up in the vicinity of faults during the time intervals between earthquakes. All computed interseismic stress accumulation is associated with the deformation caused by the time-dependent virtual displacement on major faults extending from the free surface to the seismogenic depth, i.e., the depth at which earthquakes and brittle behavior cease (~18 km depth). Stress build-up is released entirely or in part during the next large to great earthquake, with real displacements on given fault segments, which are considered positive in the model.

Earthquakes occur when the stress exceeds the strength of the fault. The closeness to failure was quantified using the change in the Coulomb failure function (ΔCFF) (modified from [6,72] and references therein). It depends on both changes in shear stress $\Delta \tau$ and normal stress $\Delta \sigma$, and in the presence of pore fluids, it takes the form of the following:

$$\Delta CFF=\Delta \tau +\mu \left(\Delta \sigma +\Delta p\right)$$

(5)

where $\Delta \tau$ is the shear stress change (computed in the slip direction), $\Delta \sigma$ is the fault-normal stress change (positive for extension), $\Delta p$ is the pore pressure change within the fault and μ is the friction coefficient, which ranges between 0.6 and 0.8 [72] and references therein). Both $\Delta \tau$ and $\Delta \sigma$ are calculated from the stress tensor and for a fault plane at the observing (field) point. When compressional normal stress on a fault plane decreases, the static friction across the fault plane also decreases. Both positive $\Delta \tau$ and $\Delta \sigma$ move a fault toward failure; negative $\Delta \tau$ and $\Delta \sigma$ move it away from failure. A positive value of ΔCFF for a particular fault denotes the movement of that fault toward failure (that is, the likelihood that it will rupture in an earthquake is increased).

Throughout this study, we ignore the time-dependent changes in pore fluid pressure and consider only the undrained case [73], meaning that $\Delta p$ depends on the fault-normal stress, whereas the fluid mass content per unit volume remains constant. Induced changes in pore pressure resulting from a change in stress under undrained conditions, according to [74], are calculated from the following:

$$\Delta p=-B\frac{\Delta {\sigma }_{kk}}{3}$$

(6)

where B is Skempton’s coefficient (0 ≤ B < 1) and Δσ_kk indicates summation over the diagonal elements of the stress tensor. If the air fills the pores, then B is nearly zero, whereas if water fills the pores, it is typically between 0.5 and 1.0 for fluid-saturated rock and close to 1.0 for fluid-saturated soil. Sparse experimental determinations of B for rocks indicate a range from 0.5 to 0.9 for granites, sandstones and marbles [74]. We assume a B = 0.5 and μ = 0.75 (as in [75], among others). These values are also proposed by [30] for the investigation of the evolutionary stress field in our study area by testing the influence of a range of apparent (effective) coefficient of friction, ${\mu }^{\prime }$ (from 0.2 up to 0.6), concluding in the value ${\mu }^{\prime }$= 0.4 as the optimal.

If in the fault zone $\Delta {\sigma }_{11}=\Delta {\sigma }_{22}=\Delta {\sigma }_{33}$, so that $\frac{\Delta {\sigma }_{kk}}{3}=\Delta \sigma$, then the apparent coefficient of friction is defined as ${\mu }^{\prime }=\mu \left(1-{\rm B}\right)$. The above-selected values for B and μ result in a value of the apparent coefficient of friction close to 0.4, which is widely used in studies of Coulomb stress modeling. In Εquation (6), Δσ_kk is the summation of the stress normal components, which, along with $\Delta \tau$, are calculated according to the fault plane solution of the next earthquake in the sequence of events, whose triggering is inspected. $\Delta \tau$ is positive for increasing shear stress in the direction of the relative slip on the observing fault, while $\Delta \sigma$ is positive for tensional normal stress. For the calculations of ΔCFF, the values of shear modulus and Poisson’s ratio are selected to be fixed at 3.3 × 10⁵ bar and 0.25, respectively.

Once the cumulative static stress changes are computed by adding the contributions from all the 7 fault segments of KTFZ, the incorporation of their respective values to the advance or delay in the earthquake cycle is estimated. The estimation is carried out by calculating the sudden time shift, $\Delta t$, to the next large earthquake, considering the state of stress since 31 December 2022 over the fault plane of the 7 fault segments and their corresponding stressing rate, $\dot{\tau }$, values (Equation (4)):

$$\Delta t=\frac{\Delta CFF}{\dot{\tau }}$$

(7)

Positive ΔCFF values result in the promotion of the next rupture and an advanced time shift of the mean recurrence time, whereas negative values can cause a negative time shift, resulting in the delay of the next rupture.

4.3. Large Earthquake (M_w ≥ 6.0) Recurrence Models and Occurrence Probabilities

The large earthquake mean recurrence time, T_r, is modeled by both a time-independent Poisson model and a renewal model, aiming at the estimation of the occurrence probabilities of M_w ≥ 6.0 earthquakes on each segment of KTFZ for the next 10, 20 and 30 years since 1 January 2023. Poisson process can be expressed by the exponential distribution with probability density function (PDF) given by the following:

$$f\left(t|T_r\right)=\frac{1}{T_r}\exp \left\{-\frac{t}{T_r}\right\}$$

(8)

where T_r is the mean recurrence time of large earthquakes associated with a specific fault segment. For modeling the large earthquake occurrence as a renewal process, the Brownian passage time (BPT) distribution [76] is applied. The PDF of the BPT model is given by the following:

$$f\left(t|T_r,\alpha \right)={\left(\frac{T_r}{2\pi {\alpha }^2{t}^3}\right)}^{1/2}\exp \left\{-\frac{{\left(t-T_r\right)}^2}{2T_r\alpha t}\right\}$$

(9)

where T_r is also the mean recurrence time and α is the aperiodicity, which can be considered analogous to the coefficient of variation, C_v, of the normal distribution.

The occurrence probabilities of the next large earthquake on a certain fault in a specific time span can be computed by applying the simple Poisson probability model given by the following:

$$P\left(t\le T\le t+\delta t\right)=1-e^{-\delta t/T_r}$$

(10)

where $\delta t$ is the forecast duration and the conditional probability corresponding to the BPT renewal given by the following:

$$P\left(t\le T\le t+\delta t\right)=\frac{{{\displaystyle \int }}_t^{t+\delta t}f\left(t\right)dt}{{{\displaystyle \int }}_t^{\infty }f\left(t\right)dt}$$

(11)

where t is relative to the date of the last earthquake, conditioned by the fact that it has been $t+\delta t$ years since the last one.

Alternatively, occurrence probabilities of the next large earthquake could be evaluated by the estimation of the hazard function, H(t), of both exponential and BPT models. Such an analysis is very useful for concluding future rupture scenarios because the values of the hazard function, or, in other words, the hazard rate, are equivalent to the conditional probability estimate in a specific time span. The hazard function of a given distribution can be easily defined using its corresponding probability density, f(t), and cumulative density, F(t), functions as follows:

$$H\left(t\right)=\frac{f\left(t\right)}{S\left(t\right)}=\frac{f\left(t\right)}{1-F\left(t\right)}$$

(12)

where S(t) is the survival function of the distribution.

For both models, the estimates of T_r and C_v derived by the application of the [65] technique are used. The results of the time shift estimates caused by the accumulated static stress changes are also inserted in the calculations. Specifically, the modified mean recurrence time, $T_r^{\prime }$, is computed by considering the time shift, $\Delta t$, (Equation (7); [7] and references therein):

$$T_r^{\prime }=T_r+\Delta t.$$

(13)

5. Application and Results

5.1. Mean Recurrence Time of Large Earthquake (M_w ≥ 6.0) Estimation

For the estimation of the large earthquake mean recurrence time, T_r, along the five right-lateral strike-slip KTFZ fault segments and the two thrust faults located in the southeastern part of Kefalonia Island, the appropriate parameters for the application of the seismic moment rate conservation method are defined. These parameters include the fault dimensions (L, w), the maximum observed magnitude, M_{_max_obs}, their long-term slip rates (V) and their corresponding uncertainties (ΔM and σ_{_V}). In this respect, the geometric and kinematic parameters of each fault segment given in Table 1 are used for the estimation of T_r. Regarding M_{max_obs}, the most recent large earthquake associated with the respective causative fault is considered if it is associated with more than one large earthquake (Lefkada North and South and Paliki North and South fault segments). This consideration is made due to the progressively higher accuracy of the magnitude estimation over time. Specifically, for the Lefkada North fault segment (S1), which is associated with the 30 June 1948 M_w = 6.4 and the 14 August 2003 M_w = 6.2 earthquakes (Table 2), the magnitude of the 2003 earthquake is selected as the M_{max_obs} (M_{max_obs} = 6.2 ± 0.2). For the Lefkada South fault segment (S2), the maximum observed magnitude is defined as equal to M_{max_obs} = 6.5 ± 0.2 since the two large earthquakes associated with this fault segment (the 24 April 1948 and the 17 November 2015; Table 2) were assigned the same magnitude.

The M_{max_obs} values of the Paliki fault segments in Kefalonia Island (Paliki North and South; S3 and S4, respectively) are considered those of the most recent earthquakes (the 2014 doublet; Table 2), which are equal to M_{max_obs} = 6.0 ± 0.2 and M_{max_obs} = 6.1 ± 0.2, for the Paliki North and South fault segments, respectively. The offshore Kefalonia fault segment (S5) is associated with only one large earthquake, the 17 January 1983, and consequently, its maximum observed magnitude is considered equal to M_{max_obs} = 7.0 ± 0.2. The M_{max_obs} for Argostoli and Ainos thrust fault segments (S6 and S7, respectively) are taken equal to M_{max_obs} = 6.5 ± 0.3 and M_{max_obs} = 7.2 ± 0.3, respectively.

Table 3. Estimated mean recurrence time, T_r, of M_w ≥ 6.0 earthquakes of the 5 strike-slip fault segments of KTFZ and the Argostoli and Ainos thrust faults, along with their respective maximum observed magnitudes, M_{max_obs} and elapsed time, T_e, (since 31 December 2022).

Fault Segment Name	Mmax_obs	Elapsed Time, Te (Years)	Monte Carlo Method		[65] Method
			Tr (Years)	67% C. I. (Years)	Tr (Years)	σ (Years)	Cv
Lefkada North	6.2 ± 0.2	19.4	60.8	[34.4, 109.9]	59.9	34.3	0.6
Lefkada South	6.5 ± 0.2	7.1	76.1	[39.1, 153.2]	79.6	32.3	0.4
Paliki North	6.0 ± 0.2	8.9	29.4	[23.0, 37.3]	33.5	17.1	0.6
Paliki South	6.1 ± 0.2	8.9	28.3	[22.7, 36.4]	33.5	16.9	0.6
Offshore Kefalonia	7.0 ± 0.2	39.9	197.2	[118.9, 301.1]	210.1	84.1	0.4
Argostoli	6.4 ± 0.3	69.4	201.5	[81.2, 487.8]	194.6	94.2	0.7
Ainos	7.2 ± 0.3	69.4	309.7	[152.4, 624.2]	298.5	148.2	0.5

summarizes the values of M_{max_obs} that are used in the calculation of mean recurrence times with both techniques and the respective T_r estimates, along with the elapsed time, T_e, since the last earthquake in each fault segment (since 31 December 2022).

The estimated T_r values (Table 3 and Figure 3a) from both approaches are in good agreement, with no systematic overestimation or underestimation of one over the other, but occasional fluctuations reported between them. Nevertheless, significant differences are observed between their corresponding confidence intervals. The confidence intervals of the Monte Carlo simulation approach (vertical red solid lines in Figure 3a) are considerably larger, encompassing both extremely smaller and larger values compared to the [65] method (Table 3). The significant variability of the Monte Carlo simulation confidence intervals is due to the larger sensitivity of this approach in the maximum magnitude and long-term slip rate uncertainties. This sensitivity is becoming even larger due to the exponential scaling of conversion of the maximum observed magnitude into the seismic moment and the selection of the uniform distribution for these simulations.

Focusing on the estimated T_r values, it is derived that they range from 33.5 years for the North and South Paliki fault segments up to almost 300 years for the Ainos thrust fault segment (Table 3 and Figure 3a). These variations are obviously related to the dimensions and the slip rate of each fault segment, resulting in different stressing rate values (Table 2) and their maximum observed magnitude. As already said, the results of the [57] approach are considered in the next steps of the current study due to their ability to provide insights about the aperiodicity of T_r, considering their C_v estimates.

Starting from the estimates for the Lefkada North and South fault segments, the [65] approach resulted in mean recurrence times equal to T_r = 59.9 years and T_r = 79.6 years, respectively. The shorter mean recurrence time of the Lefkada North segment is influenced by its larger stressing rate due to its smaller fault dimensions (Table 3). Regarding the respective C_v values, it is reported that the T_r of earthquakes with M_{_max_obs} = 6.2 on the Lefkada North fault segment exhibit slightly higher aperiodicity (C_v = 0.6) than the ones with M_{_max_obs} = 6.5 on the Lefkada South segment (C_v = 0.4). The T_r of the North and South Paliki fault is found to be equal to T_r = 33.5 years for both because they have equal dimensions and slip rates and almost the same M_{max_obs}. Their coefficient of variation values indicated a quasi-periodic to slightly aperiodic recurrence behavior (C_v = 0.6). The mean recurrence time of M_{max_obs} = 7.0 of the offshore Kefalonia fault is found to be equal to T_r = 210.1 years, exhibiting an almost quasi-periodic behavior, with a coefficient of variation value of C_v = 0.4. For the Argostoli and Ainos thrust fault segments, the mean recurrence time of M_{max_obs} = 6.5 and M_{max_obs} = 7.2 large earthquakes are found to be equal to T_r = 194.6 years and T_r= 298.5 years, respectively. The coefficient of variation for the Argostoli fault segment is the largest one, indicating high aperiodic behavior (C_v = 0.7), whereas the Ainos fault exhibits intermediate periodicity (C_v = 0.5).

One more important factor for the statistical modeling of T_r, especially for time-dependent approaches, is the elapsed time since the last earthquake. More specifically, the ratio of the T_e over T_r could be considered as an initial indicator of the stage of the earthquake cycle for a certain fault. Values of the ratio approaching zero indicate the early stage of a new earthquake cycle, while increasing values indicate that the fault is closer to the next rupture. Figure 3b shows the ratio of the elapsed time and the mean recurrence time of the seven fault segments. The smallest value is reported for the Lefkada South fault segment, which is the most recently ruptured one (T_e = 7.1 years). For the other six fault segments, the ratio ranges between 0.2 and 0.35, also showing that the elapsed time is considerably shorter than the mean recurrence time, indicating that all faults are at an early stage of their seismic cycle.

5.2. Stress Evolution of Kefalonia Transform Fault Zone

We calculated the Coulomb stress changes caused by the 11 M_w ≥ 6.0 earthquakes since 1948. We aim to examine whether later earthquakes were triggered by the earlier ones and, finally, to identify the current status of stress onto each fault segment. We are primarily interested in whether stress change on each fault, at successive stages of our examination, exceeds the commonly accepted triggering threshold of ~0.1 bars (e.g., [77,78]).

The coseismic slip, u, of the 11 earthquakes included in our model is estimated from the definition of the seismic moment, M_o. For the earthquakes that occurred after 1972, the seismic moment, which is estimated through waveform inversion, is used, whereas, for the pre-1972 earthquakes, M_o is calculated via [64] formulation. Coseismic static stress changes (ΔCFF) are computed according to the type of faulting of the target fault. Figure 4a shows the ΔCFF due to the coseismic slip of the 22 April 1948 M_w = 6.5 earthquake. A bright zone (positive ΔCFF) appears at the northeast location of the Lefkada North fault segment, most probably triggering the second 1948 earthquake (30 June 1948 with M_w = 6.4). Coseismic ΔCFF of both 1948 earthquakes, calculated according to the type of faulting of the Lefkada North fault segment (Figure 4b), developed a wide shadow zone covering the majority of Lefkada Island and a large area to NW and SE at the western and central part of Kefalonia Island. Figure 4c depicts the state of stress before the occurrence of the two 1953 earthquakes (9 August M_w = 6.4 and 12 August M_w = 7.2, calculated for the thrust faulting type. Both causative faults are inside areas of positive ΔCFF due to the tectonic stress accumulation. Cumulative static stress changes after the 12 August 1953 M_w = 7.2 earthquake (Figure 4d) created a broad stress shadow zone to NNE and SSW of the study area.

The 17 September 1972 M_w = 6.3 earthquake, associated with the Paliki North fault segment (white solid line in Figure 4e), occurred within an area with high positive ΔCFF values, as derived from the evolutionary stress field, calculated for strike-slip faulting type (Figure 4e). The state of stress before the 1983 January M_w = 7.0, associated with the offshore Kefalonia dextral fault segment, is shown in Figure 4f. This rupture occupies an area with a cumulative increase in ΔCFF since the reference year of the calculations (1948; Figure 4a–e). Cumulative ΔCFF before the 23 March 1983 M_w = 6.2 earthquake (Figure 4g) created a broad shadow zone to the east, covering the central and south parts of Kefalonia Island.

Stress accumulation before the 14 August 2003 M_w = 6.2 earthquake, resolved for the type of faulting of the Lefkada North fault segment, is shown in Figure 4h. The evidence that the causative fault (white solid line in Figure 4h) is located in a positive ΔCFF region, especially its southwestern fault tip. Cumulative static stress changes before 26 January 2014 (Figure 4i) produced an increase in ΔCFF on the Paliki Peninsula. The cumulative ΔCFF pattern before the second 2014 earthquake (3 February 2014 with M_w = 6.0; Figure 4j) created a positive stress zone in the northern part of the study area. Considering both contributions of tectonic loading and the coseismic stress changes in all previous earthquakes, Figure 4k shows the stress evolution just before the 17 November 2015 M_w = 6.5 earthquake in the southern part of Lefkada Island. The northern tip of the causative fault (Lefkada South; white solid line in Figure 4k) is in a positive static stress change area. Figure 4l shows the cumulative ΔCFF pattern after the 2015 M_w = 6.5, the coseismic slip of which created an extended shadow zone along both Lefkada fault segments, while it enhances the stresses in the northern and central part of Kefalonia island.

The evolutionary stress field calculation is further extended until 2022 (74 years in total), even though in 2015, the last earthquake occurred, aiming after determining the current state of stress to calculate the occurrence probabilities incorporating the current ΔCFF in each fault segment. The evolved stress field until 2022 is calculated after considering both tectonic loading and static stress changes calculated at a depth of 10 km. The calculation is implemented by considering three representative faulting types for the Lefkada branch of KTFZ (strike: 20°, dip: 62°, rake: 180°), the Kefalonia branch of KTFZ (strike: 30°, dip: 50°, rake: 180°) and the thrust faulting (strike: 300°, dip: 30°, rake: 95°).

Figure 5a shows the state of stress on 31 December 2022 for the Lefkada branch faulting type. The Lefkada fault segments are located in a stress shadow region caused by the 2015 earthquake. The Kefalonia branch fault segments (Paliki North and South and offshore Kefalonia; Figure 5b) are located inside stress-enhanced regions, with the Paliki South being inside an area with the highest positive ΔCFF values. Figure 5c depicts the state of stress calculated for thrust faulting (Argostoli and Ainos thrust fault segments), where it is observed that both segments are inside either negative (to the west) or positive (to the east) static stress change values.

The accumulated static stress changes are also resolved onto the rupture plane of each fault segment (Figure 6). From these calculations, which are summarized in

Table 4. Minimum, average and maximum Coulomb stress change values calculated onto the rupture fault plane of the five strike-slip fault segments of KTFZ and the Argostoli and Ainos thrust faults, along with the respective time shift, Δt, values calculated by averaging ΔCFF values.

Fault Segment Name	ΔCFF (Bar)			Time Shift, Δt (Years)
	Min	Aver	Max	Aver
Lefkada North	−637.02	−34.03	124.88	−40.23
Lefkada South	−641.06	−60.18	490.64	−87.14
Paliki North	−355.29	19.91	669.70	10.45
Paliki South	−8.44	3.48	28.62	1.83
Offshore Kefalonia	−44.69	10.17	53.48	11.95
Argostoli	−656.51	22.07	715.62	63.17
Ainos	−615.85	5.05	535.77	27.96

, it is clearly derived that both Lefkada fault segments (Lefkada North and South, Figure 6a and 6b, respectively) attain negative average values equal to ΔCFF = −34.03 bar and ΔCFF = −60.18 bar, for the Lefkada North and South fault segments, respectively. Limited areas with high positive static stress change values are also observed (Figure 6a,b; Table 4). The rupture planes of the Paliki North and South segments of Kefalonia Island (Figure 6c and 6d, respectively) are mainly covered by high positive ΔCFF values, with an average equal to ΔCFF = 19.91 bar and ΔCFF = 3.48 bar, respectively. The considerably larger average ΔCFF values of the Paliki North fault segment are caused because it incorporates most of the enhanced stress area created from the 2015 earthquake. The offshore Kefalonia fault segment (Figure 6e) is also covered by intermediate to large positive static stress change values with an average equal to ΔCFF = 10.17 bar. Static stress changes onto the Argostoli and Ainos fault segments (Figure 6f and 6g, respectively) demonstrate a more complex pattern of both high negative and positive values, especially onto the Ainos segment, exhibiting large variations (Table 4). The average values of ΔCFF are found to be equal to ΔCFF = 22.07 bar and ΔCFF = 5.05 bar for Argostoli and Ainos fault segments, respectively, highlighting that they are affected mainly by the positive static stress changes.

The effect of accumulated static stress changes is then considered for calculating the time shift, Δt (Equation (7)). Although it is widely recognized that a higher percentage of earthquakes occur in regions of enhanced static stress changes rather than in stress shadows (e.g., [23,79]), in the present study, the average of both negative and positive ΔCFF, as shown in Table 4, are considered. This choice is made because it is not known where nucleation will take place, and thus, the stress conditions onto fault planes need to be more representative. The results of these calculations are shown in the last column of Table 4, indicating that in the Lefkada North and South fault segments, the stress accumulation results in a next rupture delaying time shift, whereas in all other cases, a promoting time shift is reported.

5.3. Recurrence Models and Occurrence Probabilities

We applied the exponential and BPT distributions, representing the Poisson and renewal models, with the mean recurrence time, T_r, and C_v values obtained by a [65] method (Table 3). Additionally, the incorporation of the state of stress until 2022 is taken into account by calculating the modified mean recurrence time, $T_r^{\prime }$, as defined from Equation (13) by considering the respective time shift for each fault segment. In this respect, the statistical analysis focuses on four different models: the simple Poisson and BPT models only based on the T_r estimates, and the Poisson and BPT models with the ΔCFF incorporation and the $T_r^{\prime }$.

Figure 7 shows the hazard functions of all four statistical models applied, in which their significant differences in modeling the T_r are highlighted. Starting from the Lefkada North fault segment (Figure 7a), the constant hazard rate of the simple Poisson model (blue solid line) is considerably larger than the hazard rate of the BPT model (red solid line), which exhibits an increasing trend, at the time corresponding to T_e. The effect of the ΔCFF on both recurrence models (Poisson+ ΔCFF and BPT+ ΔCFF) significantly influences the associated results, whereas the hazard rates of both statistical distributions are reduced, indicating a delay in the occurrence of the next large earthquake. This influence is clearer in the case of the Lefkada South fault segment, which is the last rupture within the 74 years of our analysis. Both Poisson and Poisson+ ΔCFF models give rather low hazard rates, whereas the hazard functions of both renewal models (BPT and BPT+ ΔCFF; red solid and magenta dashed lines, respectively) show that the fault segment is at an early stage of a new earthquake cycle at time t corresponding to the elapsed time (vertical black solid line).

On the other hand, the BPT model application on the Paliki North fault segment (Figure 7c) reports hazard rates almost equal to the hazard rate of the simple Poisson model. Hazard rate values of both models significantly increase if the incorporation of ΔCFF is considered since this is the fault segment for which a large advance time shift is calculated (Table 4). Both renewal models (BPT and BPT+ ΔCFF; red solid and magenta dashed lines, respectively) indicate that the Paliki North fault segment is at an advanced stage of a new earthquake cycle, even though its elapsed time is short. This Is due to the low values of T_r, its large stressing rate (the largest, along with the Paliki South, among the fault segments of KTFZ) and the significant influence of the static stress changes. The Paliki South fault segment is affected less than the previous one from the state of stress until 31 December 2022. Comparing even the time-independent or the renewal models (Figure 7d), it is observed that no significant advancement in T_r is reported. Same as in the latter case, the offshore Kefalonia fault segment (Figure 7e) does not exhibit a significant influence of its positive time shift (Table 4) because its T_r is considerably larger than both Δt and the corresponding elapsed time. Focusing on the hazard function curves of both Poisson and Poisson+ ΔCFF and BPT and BPT+ ΔCFF models, one can observe that they almost coincide.

The application of the four models on the T_r of Argostoli fault segment shows that hazard rate values of the simple Poisson and BPT models are almost equal at a time corresponding to the elapsed time, resulting in the same occurrence rates. The positive static stress changes onto the Argostoli fault segment, resulting in an advanced time shift, moving the new earthquake cycle earlier (dashed green and magenta lines in Figure 7f for the Poisson+ ΔCFF and BPT+ ΔCFF models, respectively). This fact is also affected by the aperiodicity of large earthquakes, T_r, of this fault segment, which is found to get the highest value (α = 0.7). The Ainos fault segment exhibits the same recurrence behavior as the offshore Kefalonia segment. Specifically, its positive time shift does not significantly influence the hazard function of the applied statistical models (Figure 7g) since both Poisson and Poisson+ ΔCFF models and BPT and BPT+ ΔCFF are slightly different.

The aforementioned results are quantitatively confirmed by the calculation of the occurrence probabilities based on all four models (Poisson, Poisson+ ΔCFF, BPT and BPT+ ΔCFF) by applying Equations (10) and (11) for the next 10, 20 and 30 years.

Table 5. Occurrence probabilities of M_w ≥ 6.0 earthquakes for the next 10, 20 and 30 years since 31 December 2022 for the five strike-slip fault segments of KTFZ and the Argostoli and Ainos thrust faults, according to Poisson, Poisson+ ΔCFF, BPT and BPT+ ΔCFF models.

Fault Segment Name	Poisson			Poisson + ΔCFF			BPT			BPT + ΔCFF
	P10	P20	P30	P10	P20	P30	P10	P20	P30	P10	P20	P30
Lefkada North	0.15	0.27	0.38	0.10	0.18	0.26	0.12	0.29	0.44	0.02	0.08	0.16
Lefkada South	0.12	0.22	0.32	0.06	0.11	0.16	10 × 10−6	3 × 10−³	0.03	2 × 10−12	2 × 10−7	3 × 10−5
Paliki North	0.26	0.45	0.59	0.40	0.62	0.70	0.19	0.47	0.57	0.50	0.78	0.89
Paliki South	0.26	0.45	0.59	0.27	0.46	0.60	0.19	0.47	0.57	0.24	0.52	0.72
Offshore Kefalonia	0.05	0.09	0.13	0.05	0.09	0.14	6 × 10−5	5 × 10−4	3 × 10−3	2 × 10−4	10 × 10−4	6 × 10−3
Argostoli	0.05	0.09	0.14	0.06	0.12	0.17	0.04	0.09	0.14	0.07	0.14	0.21
Ainos	0.03	0.06	0.09	0.04	0.07	0.10	2 × 10−3	6 × 10−3	0.01	3 × 10−3	9 × 10−3	0.02

and Figure 8 summarize the results. According to the simple Poisson model, the probability estimations vary from low to intermediate for all segments for the next 10 years, whereas they will increase up to 40% for large earthquakes for the Paliki North and South fault segments for the next 20 and 30 years (Table 5 and Figure 8b,c). The influence of the time shift caused by the state of stress until 2022, as previously discussed, is represented by the occurrence probabilities of the Poisson + ΔCFF model. Specifically, a half-value decrease in the 10-year probability (P₁₀ = 0.06 instead of P₁₀ = 0.12 for the Poisson and Poisson + ΔCFF models, respectively; Table 5) is observed for the Lefkada South fault segment (S2 in Figure 8d). Further, an increase is reported for the Paliki North fault segment with the 10-year probability of the Poisson + ΔCFF model to be equal to P₁₀ = 0.40 instead of the one of the simple Poisson (P₁₀ = 0.26). The 20- and 30-year occurrence probabilities do not exhibit significant differences between the two time-independent models.

Focusing on the time-dependent probabilities, namely the ones obtained from the BPT and BPT + ΔCFF models, it is derived that the 10-year occurrence probabilities range from low to intermediate values for all the fault segments (Table 5 and Figure 8g). The Paliki North and South fault segments reported probability values up to 45% for the next 20 and 30 years for both time-dependent models (P₂₀ = 47% and P₃₀ = 57% for both segments according to the simple BPT model; P₂₀ = 78% and P₂₀ = 52% and P₃₀ = 89% and P₃₀ = 72% for the Paliki North and South, respectively, calculated with the BPT + ΔCFF model). It should be mentioned that for the Paliki South fault segment, the 20-end 30-year probabilities are slightly decreased if the influence of static stress is considered. A slight increase from low (p < 10%) to intermediate (10% ≤ p < 40%) values is observed for the Agrostoli fault segment from the application of the BPT + ΔCFF model. On the other hand, occurrence probabilities of the Lefkada North segment decrease if the BPT + ΔCFF model (Figure 8g–l) is considered for all cases of the next 10, 20 and 30 years, whereas the ones of the Lefkada South fault segment are found to be lower than 10% in both time-dependent models (Table 5).

In summary, the Paliki North fault segment gets the highest occurrence probability values in all models because of its relatively short mean recurrence time, its intermediate aperiodicity, its large stressing rate and the great influence of stress transfer onto the fault plane. On the contrary, the lowest values are observed in the Lefkada South fault segment, which reports no greater than 32% probability (the highest reported value for the next 30 years according to the simple Poisson model; Table 5), due to its relatively large mean recurrence time, its relatively low stressing rate, along with the great decrease in the static stress caused by the 2015 earthquake. The offshore Kefalonia and Ainos fault segments, which are associated with the largest M_{max_obs} earthquakes during the reference period (M_w = 7.0 and M_w = 7.0, respectively), report low to intermediate probability values for all models due to their large recurrence times in respect with the elapsed time and their relatively low stressing rate values.

6. Concluding Remarks

Four earthquake occurrence models (Poisson, Poisson+ ΔCFF, BPT, BPT+ ΔCFF) are applied, aiming at the calculation of occurrence probabilities of nearly characteristic (M_w ≥ 6.0) earthquakes in the central Ionian Islands. We used a detailed segmentation model for the estimation of the mean recurrence time, T_r, by applying the physics-based seismic moment rate conservation method. The results show that T_r ranges from about 30 years up to almost 300 years. The Paliki peninsula fault segments (Paliki North and South) exhibit the shortest T_r values due to their high slip rates and their relatively small dimensions. Intermediate T_r values are observed in the two segments along the western coasts of Lefkada Island (Lefkada North and South fault segments) and for the Argostoli thrust fault segment. The largest mean recurrence times are observed for the offshore Kefalonia and Ainos thrust fault segments due to both their significantly larger dimensions and their M_{max_obs} (M_w = 7.0 and M_w = 7.0, respectively).

The incorporation of the 74-year (1948–2022) evolutionary stress field into the recurrence pattern modeling indicates that the mainly affected fault segments are the Lefkada South and Paliki North ones. On the one hand, the stress accumulation of the Lefkada South fault segment reveals a significant time delay, as it is the last ruptured fault segment. On the other hand, the Paliki North fault segment was affected by high static stress changes caused by both the static stress changes in the 2015 earthquake and its high stressing rate. For the other five fault segments, either intermediate negative (Lefkada North fault segment) or positive (Paliki South, offshore Kefalonia, Argostoli and Ainos fault segments), static stress changes are calculated. The least influence of the evolutionary stress is observed in the offshore Kefalonia and Ainos fault segments due to their large fault area and the high variable distribution of ΔCFF onto their fault planes. Overall, it could be stated that the influence of the evolutionary stress field is affected by both the geometric and kinematic parameters of the fault segments.

Occurrence probability estimations for the next large earthquake during the next 10, 20 and 30 years, associated with the seven KTFZ fault segments, reveal that the most probable fault segment to be ruptured is the Paliki North fault segment for an M_w = 6.0 earthquake, according to all four models. For example, the 20-year occurrence probabilities are equal to P₂₀ = 45% and P₂₀ = 47%, according to the Poisson and the BPT models. Taking into account the incorporation of static stress changes, these probabilities significantly increase to 62% and 78% for the Poisson + ΔCFF and the BPT + ΔCFF models, respectively. Although both time-independent and time-dependent modeling results cannot be given any preference for seismic hazard studies, since there is no statistical evidence of one performing better than the other, their results could equally be considered as the basis of fault-based large earthquake occurrence models.