Spatiotemporal Properties of Seismicity and Variations of Shear-Wave Splitting Parameters in the Western Gulf of Corinth (Greece)

: The Western Gulf of Corinth (WGoC) exhibits signiﬁcant seismicity patterns, combining intense microseismic background activity with both seismic swarms and short-lived aftershock sequences. Herein, we present a catalogue of ~9000 events, derived by manual analysis and double-difference relocation, for the seismicity of the WGoC during 2013–2014. The high spatial resolution of the hypocentral distribution permitted the delineation of the activated structures and their relation to major mapped faults on the surface. The spatiotemporal analysis of seismicity revealed a 32-km-long earthquake migration pattern, related to pore-pressure diffusion, triggering moderate mainshock-aftershock sequences, as ﬂuids propagated eastwards in the course of ~15 months. The anisotropic properties of the upper crust were examined through automatic shear-wave splitting (SWS) analysis, with over 2000 SWS measurements at local stations. An average fast shear-wave polarization direction of N98.8 ◦ E ± 2.8 ◦ was determined, consistent with the direction of the maximum horizontal regional stress. Temporal variations of normalized time-delays between fast and slow shear-waves imply alterations in the level of stress or microcrack ﬂuid saturation during the long-lasting pore-pressure diffusion episode, particularly before major events. The present study provides novel insights regarding seismicity patterns, active fault structures, anisotropic properties of the upper crust and triggering mechanisms of seismicity in the WGoC. seismic swarms and mainshock-aftershock sequences, presenting complex dynamics. We apply double-difference relocation to reduce the relative location uncertainties between clustered events and to enhance the spatial resolution of the catalogue. We examine the spatiotemporal evolution of seismicity, focusing on earthquake migration patterns. We examine phenomena related to the diffusion of ﬂuids or aseismic creep by means of earthquake migration speeds, hydraulic diffusivities and diffusion exponents, with respect to potential origin points of ﬂuid injection. The above are studied in conjunction with the temporal variations of the SWS parameters,


Introduction
The Western Gulf of Corinth (WGoC), located in Central Greece, is one of the bestmonitored seismogenic sites in Europe. Its high deformation rate (14-15 mm·yr −1 extension in a N9 • E direction; [1,2]) is the driving force of the intense microseismic activity that has been recorded by the Corinth Rift Laboratory (CRL) local seismological network since 2000 [3]. The background microseismicity is interspersed by seismic swarms (e.g., [4][5][6][7]) and earthquake sequences [8][9][10]. The WGoC has been the site of major events, with the most recent one being the 15 June 1995 M s = 6.2 Aigion earthquake [11]. The Gulf of Corinth is an active continental rift, formed by conjugate E-W-trending normal faults ( Figure 1; [12,13]). Major structures of the southern flank of the WGoC include the north-dipping Aigion, Helike and Pirgaki faults, as well as the Psathopyrgos fault [14], linked with the Aigion fault via the en echelon Lambiri, Selianitika and Fasouleika NW-SE to WNW-ESE-trending fault segments (Kamarai fault system [11]). In the Greek Database of Seismogenic Sources [14], the Aigion fault is represented by a crustal-scale individual seismogenic source (GRIS503; Figure S20), whereas the Fasouleika and Selianitika faults are minor branches. Opposite to these structures, at the northern shore, the SW-NE-trending, SE-dipping Platanitis and Nafpaktos faults are mapped. Towards the east, the south-dipping Marathias and Kallithea faults are aligned E-W and WNW-ESE, respectively. Another south-dipping structure to the north of Marathias is the Trikorfo fault [14], extending eastwards to the Sergoula fault [15]. Offshore faults have also been mapped, such as the south-dipping Trizonia fault, the West Channel fault and the North and South Eratini faults [16,17].
The background seismicity in the WGoC is scattered with epicenters, mainly located offshore between Aigion and Trizonia Island, but also north of the Psathopyrgos fault. Background seismicity rates are of the order of 10 events/day, automatically detected by the local CRL network, reaching hundreds of events/day during earthquake swarms or aftershocks [5,8]. This seismic activity, driven by high strain-rates, is associated with a fragmented seismogenic low-angle, north-dipping layer below the gulf, where the major normal faults of the region are rooting. Within the background seismic activity in the WGoC, long-term repeaters have been identified, associated with recurring slip on the same fault patch over the course of several years [9]. Such observations indicate the existence of small asperities surrounded by a partially creeping fault surface. Onshore seismicity is mostly projected towards the northern shores, near Efpalio and Eratini, mainly associated with activity at the deeper part of a north-dipping growing detachment (e.g., [8,9]).
Herein, we investigate the spatiotemporal evolution of the earthquake activity in the WGoC and examine patterns that could be associated with the involvement of aseismic factors driving the seismicity, such as propagating fluids or creeping of the surrounding fault surface. For this purpose, we also examine the anisotropic properties of the crust in terms of their statistics and temporal variations, which could be related with alterations in the level of stress or pore space fluid saturation. Seismic anisotropy in the upper crust is usually attributed to the presence of microcracks with a specific orientation [18]. As shear-waves propagate through an anisotropic medium, the Shear-Wave Splitting (SWS) phenomenon occurs, producing a fast shear-wave (S fast ) that becomes polarized to the direction favored by the medium's anisotropic properties, whereas a slower one (S slow ) is polarized to a perpendicular direction. The Anisotropic Poro-Elasticity (APE) model considers the presence of vertical, fluid-saturated microcracks in the upper crust, taking into account the influence of the physical properties of the fluids [18]. The APE model attributes the SWS properties to the balance between the pore-fluid pressure and the local stress field, particularly its maximum horizontal component, σ Hmax [19,20]. It also ascribes variations of the SWS parameters to stress field changes, enabling their correlation with fluid processes, such as diffusion. Changes in the microcrack orientations or fluid saturation have been related to pore-pressure diffusion mechanisms associated with precursory stressload on faults or post-seismic fluid release [21]. In Greece, a significant number of SWS studies have been performed. In the Eastern Gulf of Corinth [22][23][24], the Santorini Volcanic Complex [25] and Marathon [26], upper-crust anisotropy was attributed to the presence of stress-oriented fluid-filled microcracks, according to the APE model. However, there have also been observations where SWS has been related to CO 2 degassing processes [27]. Since 1990, studies in the WGoC have revealed a relation between σ Hmax and the polarization direction of the fast S-waves (e.g., [11,28]). Since 2000, with the installation of the CRL network [3], the large amount of recorded microseismic events has rendered the region an ideal site for the study of variations of the SWS parameters.
In this work, we present a high-resolution earthquake catalogue of~9000 events for the WGoC during 2013-2014, along with a large dataset of shear-wave splitting parameters. This period is characterized by several seismic swarms and mainshock-aftershock sequences, presenting complex dynamics. We apply double-difference relocation to reduce the relative location uncertainties between clustered events and to enhance the spatial resolution of the catalogue. We examine the spatiotemporal evolution of seismicity, focusing on earthquake migration patterns. We examine phenomena related to the diffusion of fluids or aseismic creep by means of earthquake migration speeds, hydraulic diffusivities and diffusion exponents, with respect to potential origin points of fluid injection. The above are studied in conjunction with the temporal variations of the SWS parameters, focusing on periods characterized by the occurrence of seismic swarms or major events producing aftershock sequences. We determined more than 2000 fast shear-wave polarization direction and time-delay measurements over a period of two years; a number significantly larger from what is included in most upper-crust SWS studies. The wealth of measurements ensures the validity and the reliability of the herein presented results. We examine the spatial distribution of relocated hypocenters through detailed cross-sections and provide a 3D interactive visualization to showcase the relation between the inferred activated structures at depth and the mapped faults on the surface. The present study aims to provide a high-resolution view of the active fault geometries in the WGoC and their interplay, as seismicity migrates from one structure to another. We further investigate and discuss the possible triggering mechanisms of seismicity associated with the observed migration patterns, providing further insights from the SWS analysis of the upper crust in the area of the WGoC. . Beachballs depict the focal mechanisms of significant earthquakes: red after NKUA-SL, blue after GI-NOA, green for the 2010 Efpalio earthquakes after [29], yellow for the 1995 M s = 6.2 Aigion earthquake [11] and gray for the 2001 Agios Ioannis M w = 4.2 event [3,4]. Triangles depict seismological or accelerometric stations that operated during the study period (2013-2014), with colours corresponding to different network codes (legend at the bottom-left). Fault lines after [4,12,[15][16][17]30]. Pt.f.: Panachaikon transfer-fault.

Materials and Methods
The primary data used in this study are the recordings of the seismological stations of the CRL network and of the local stations of the Hellenic Unified Seismological Network (HUSN). These were acquired through the web services of the international Federation of Digital Seismograph Networks (FDSN); specifically, from the Réseau sismologique et géodésique français (RESIF; http://seismology.resif.fr/) (accessed on 9 July 2021), for stations of the CRL network, and from the National Observatory of Athens data center for the European Integrated Data Archive (EIDA@NOA; http://eida.gein.noa.gr/) (accessed on 9 July 2021), for stations of the HUSN [31]. In the following subsections, we present the various methodologies that were followed to process the available data and produce the results which are further described in Section 3.

Location and Relocation of the Catalogue
We employed a dataset of over 9000 events that occurred during 2013-2014 in the WGoC. Preliminary event detections and arrival-time data for these events were available from the websites of the Geodynamics Institute of the National Observatory of Athens (GI-NOA; https://bbnet.gein.noa.gr/HL/) (accessed on 9 July 2021) and CRL http:// ephesite.ens.fr/~eworm/ (accessed on 9 July 2021). Additional P-and S-wave arrival time picks at local stations of CRL and HUSN were manually added to enrich the dataset, which herein is referred to as "Catalogue 1", or briefly, Cat1.
Initial locations were obtained with the HypoInverse code [32], using the local velocity model of [33] for the WGoC. The only exception is the seismic cluster of the 2013 Helike swarm, for which an optimal velocity model [6] was employed. We applied a first approximation of station corrections, by removing the average P-or S-wave travel-time residual on each station and recalculating locations. At this point, we separated the dataset into spatial groups based on the hypocentral distribution. Initial locations indicated that seismicity related to the shallow structures of the WGoC reaches a depth of~15 km. Seismicity at depths Z < 15 km was divided into clusters with the application of Ward's [34] linkage to the matrix of interevent hypocentral distances. Smaller clusters were merged by visual inspection to produce a total number of 9 spatial groups ( Figure 2). The seismicity around the margins of the study area or deeper than 15 km was placed in a 10th "miscellaneous" group.
All the following procedures were performed independently for each spatial group. At least one reference station, in close proximity to the respective epicenters, was selected to perform waveform cross-correlation. Each recording was band-pass-filtered in the range 2-15 Hz (e.g., [6,25]) and then cropped in windows containing P, S and coda waves. A correlation matrix was constructed by registering the global maximum (XC max ) of the cross-correlation function between all event-pairs with available recordings at the reference stations. Then, a combined correlation matrix was calculated, taking into account the maximum value of the individual matrices for each event-pair. Nearest-neighbor linkage was performed and each group was divided into multiplets by selecting an optimal threshold, C th . The latter is determined as the value that maximizes the difference between the size of the largest cluster and the sum of clustered events [35], with a minimum of 0.60. Finally, for event-pairs within each multiplet, cross-correlation was performed on the recordings of all stations with available phase-pick data, in windows containing only P-or S-waves, after alignment of the waveforms on their respective P or S arrival-times. The XC max and its corresponding time-lag were registered on each case.
These data, along with catalogue travel-times, were used as input for the doubledifference relocation code, HypoDD [36]. This algorithm performs iterative inversion to minimize the double-difference between observed and calculated travel-times in pairs of neighboring events. The HypoDD code relocates hypocenters by attributing their traveltime differences to their inter-event distance, reducing uncertainties caused by unmodelled velocity structure and arrival-time reading errors. As a result, spatially correlated events become more concentrated and aligned, increasing resolution and highlighting the activated structures. Location statistics for Catalogue 1 are presented in the Supplementary Section S1. In addition, magnitudes were calculated using a spectral fitting method, similar to the one applied by [6] for the 2013 Helike swarm. Details on this procedure are described in the Supplementary Section S2. Compared to the initial locations (Figure 2a), the map of relocated epicenters ( Figure 2b) presents a clearer, enhanced image of the spatial distribution of the seismicity, with better defined sub-clusters within each of the 9 groups and fewer sparse events in the gaps among clusters. A robustness test for the relocation procedure is described in Section S1 of the Supplementary Material. Different colours and numbers represent the 9 spatial groups to which the dataset has been divided, in addition to the "miscellaneous" group (#10), presented in gray.
Besides Cat1, which contains a selection of manually revised events, the location of additional events was sought to provide more candidate events for SWS measurements. For this purpose, the automatic catalogue of CRL http://ephesite.ens.fr/~eworm/ (accessed on 9 July 2021) was employed, containing over 43,000 events that occurred during 2013-2014. After removing events that were already available in Cat1, a quick examination by visual inspection was performed. Events with low signal-to-noise ratio (SNR), or iden-tified as regional/out-of-network events being erroneously located inside the gulf, were discarded. Afterwards, automatic picking was performed by applying auto-regression [35] and correlation detector methods (e.g., [37]), using the manually picked P-and S-waves as templates (see Section S3 in the Supplementary Material). In total, 10,493 events were automatically resolved and located with HypoInverse. As the solutions of these additional events are considered to be of lower quality, they were placed in a separate "Catalogue 2", or Cat2. More details on Cat2 are presented in Supplementary Section S3 ( Figure S7). The combination of the two datasets will be referred to as "Catalogue 3", or Cat3.

Shear-Wave Splitting (SWS) Analysis
A preliminary set of 14,825 events from both the manual (Cat1) and the automatic (Cat2) catalogues were selected as candidate events for an upper-crust SWS analysis using local stations. The first selection criterion is the "shear-wave window" [38], which requires the S-wave seismic rays to arrive at an incidence angle smaller than 45 • , to avoid the influence of converted pS phases. The Signal-to-Noise Ratio (SNR) is measured for a signal window containing the S-wave first pulses. This is performed for a pre-defined range of frequency bands used for the application of a band-pass Butterworth filter, and the filter yielding the highest SNR is selected [39]. Waveform recordings with SNR < 1.5 are discarded.
Seismic anisotropy at the upper crust is studied via SWS, which is quantified by two parameters: a) the polarization direction, ϕ, of the fast shear-wave and b) the timedelay, t d , between the arrivals of the two split shear-waves at a receiver. To process this large dataset, we employed the open-source software Pytheas [40], which offers a variety of popular SWS analysis methods and has successfully been applied in recent crustal anisotropy studies (e.g., [26,41]). We adopt the Eigenvalue (EV) [42], rather than the Rotation-Correlation (RC) method [43], as the former is considered to be more reliable [44]. The EV method performs a grid-search at a pre-defined range of ϕ, t d value pairs. For each grid point, the horizontal waveforms are rotated to all candidate S fast and S slow coordinates, according to the examined ϕ, and then S slow is shifted by -t d for all time-delay values in the pre-defined range. The covariance matrix of the corrected waveforms is calculated in a window containing the S-wave pulse, and the major and minor eigenvalues λ 1 and λ 2 are determined. The solution that yields the minimum λ 2 is considered as the optimal pair of ϕ and t d . A 95% confidence interval in the range of λ 2 values around the optimal grid point is considered to define the measurement uncertainties δϕ and δt d ( Figure 3).
An important advantage of the Pytheas software is that it permits fully automatic SWS processing by defining the optimum S-wave signal window. This is performed with the application of a Cluster Analysis procedure [45]. A series of different window configurations are examined, within the bounds of the middle of S-P interval and the dominant period of the S-wave pulse multiplied by a user-defined factor [39,40]. For each window, ϕ and t d are measured using the EV method and the solution pairs are sorted into clusters [45]. Then, the window yielding the (ϕ, t d ) solution pair with the smallest errors from the most constrained cluster is selected. The quality of the SWS parameter measurements is assessed by automatically calculating a score value, taking into account the ratio between δϕ or δt d and a respective user-defined maximum threshold. A quality grade between A (best) and E (worst) is then attributed, depending on their resulting average score. The E grade is also reserved for cases of low SNR < 1.5. In the following, only measurements with grades A, B or C were taken into consideration. Furthermore, measurements are marked as "null" when their ϕ value is within ±10 • from p or p + 90 • , i.e., sub-parallel or perpendicular to the polarization direction, p, of the corrected (for anisotropy) shear-wave. In those cases, it is not possible to distinguish whether the medium is isotropic or the initial polarization coincides with that of the S fast or the S slow , which is why "null" measurements are discarded [46]. Example of the automatic determination of the SWS parameters with the Pytheas software, using the Eigenvalue (EV) and Cluster Analysis (CA) methods. Horizontal recordings (a) in their initial N-S/E-W co-ordinate system (N-E), after the application of band-pass Butterworth filter between 0.5 and 5.0 Hz, (b) rotated into the S fast /S slow co-ordinate system (F-S), having measured ϕ = N94 • E, then shifted S slow leftwards by t d = 60 ms and (c) rotated back to the N-S co-ordinate system, after being corrected for anisotropy, where the p = N58.8 • E polarization direction is measured. The S-wave particle motion in the N-S co-ordinate system is also presented in the hodograms of panels (d) for the initial waveforms and (e) after correction for anisotropy. (f) Contours of the distribution of λ 2 values during the grid-search of (ϕ, t d ) pairs performed for the application of the EV method, with the polarization direction displayed in the form of ϕ-baz, i.e., degrees relative to the back-azimuth, wrapped in a range of ±90 • with respect to the R-axis direction. The selected (ϕ, t d ) pair, corresponding to the global minimum of λ 2 , is marked by the dashed cross-hair, pointing to the middle of a dark-shaded region that indicates the 95% confidence interval for λ 2 [42], from which δϕ = 0.3 • and δt d = 0.03 ms are measured. (g-k) Diagrams of the Cluster Analysis (CA) method for the determination of the best window presenting: (g) the (ϕ, t d ) measurement pairs, (h) the performed clustering, (i) the selection of the most constrained cluster and measurement with the smallest errors, (j) the series of ϕ and (k) the series of t d measurements per trial window index, with their respective uncertainties. Figure 3 shows an example of the automatic analysis with the combination of the EV and CA methods for an event recorded at station PSAR. Initially, the automatic filtering algorithm decided the optimal band-pass Butterworth filter range between 0.5 and 5.0 Hz. Then, SWS measurements were performed in 102 trial signal windows. For each window, the EV method was applied to determine the ϕ and t d values. The horizontal components were first rotated into the Radial-Transverse (R-T). Then, grid-search was performed, rotating the waveforms to the various ϕ values between R−90 • T and R90 • T, and shifting the transverse component leftwards by a time-delay t d between 0.0 and 250.0 ms. In this example, the solution pair that yielded the minimum λ 2 value was ϕ = N94.0 • E ± 0.3 • and t d = (60.00 ± 0.03) ms in a 95% c.i. with grade "A", indicating a measurement of the best quality. Figure 3g presents the selection of the optimal solution. It is noteworthy that despite the ϕ values being spread in a range roughly between 80 • and 120 • , the majority are clustered near the selected value of 94 • , with the respective t d at 60 ms (the preferred value). The diagram of Figure 3f shows the global minimum of λ 2 for the selected signal window, marked by the dashed cross-hair. The validity of the ϕ, t d measurements is verified in Figure 3e, where the resulting particle motion becomes nearly linear after the removal of anisotropy.
We verified the instrument orientation during the study period (2013-2014) using the method of [47], to ensure the validity of the polarization direction measurements. This method compares the polarization direction of teleseismic P-waves with the respective back-azimuth of the epicenters to determine the orientation of the positive direction of the instrument's N-S component, which should be pointing towards the North. The results were compared with the orientations documented in the metadata of the stations supplied by the FDSN web services. Significant differences were identified at stations MG06 (N-S azimuth = N83 • E) and the temporary station HELI (N-S azimuth = N340 • E). The method of [47] was not applied to station ALIK, as it could not properly record the teleseismic waveforms. In that case, an orientation correction of N299 • E was adopted after [48], who followed a similar approach, but using the P-wave polarization determined by local earthquakes [35]. For other stations, especially for borehole instruments, the nominal orientations supplied in the stations' metadata were verified.
Out of the 14,825 candidate events from both catalogues, a total of 37,263 initial SWS results (15,083 from Cat1 and 22,180 from Cat2) were acquired, with 10,456 events providing a result from at least one station. Other events were discarded due to low SNR or errors during processing, in particular at the clustering stage. In the combined Cat3, 25,266 measurements with grade E and 4223 with grade D were rejected due to low quality. In addition, 3298 measurements were considered as null, having ϕ values sub-parallel or sub-perpendicular to the corrected source polarization, p. Finally, in Cat3, 192 valid measurements are with grade A, 1299 with grade B and 2985 with grade C, summing up to 4476 (Table S3 in the Supplementary Material).

Spatial Distribution of Seismicity
A total of 8944 events of Cat1 were successfully relocated. As presented in the crosssections of Figure 4, seismicity is mainly concentrated on a low-angle seismogenic layer at depths between 6.5 km and 10 km below the gulf (see also Figure S1d), reaching~12 km down-dip the activated structures. Several clusters, mainly related to the occurrence of larger events, indicate the activation of higher-angle faults, rooting at the weak layer. The latter is further supported by the reported focal mechanisms ( Table 1). The magnitudes of major events mentioned in the present study refer to their respective M w values calculated herein by spectral fitting (see Section S2 in the Supplementary Material). Table 1. Focal mechanisms of significant events during 2013-2014 in the WGoC. Sources: NKUA 1 (Seismological Laboratory of the National and Kapodistrian University of Athens), NOA 2 (Geodynamics Institute-National Observatory of Athens) and K2015 [6]. Latitude, longitude, (centroid) depth and moment magnitude M w * values are those reported from the respective reference. M w is the moment magnitude determined by spectral fitting in this study and CLID is the spatial group ID (see Figure 4). Three important events of the dataset, related to specific sequences, are marked with bold.

Date
Hour The westernmost part of the Nafpaktos-Psathopyrgos group #4 (green), activated after the 21 September 2014 M w = 4.6 event, could possibly be associated with the Rion-Patras SW-NE-trending dextral oblique-normal fault zone, or with related sub-parallel structures dipping towards NW. The N10 • E cross-sections of Figure S13 imply gradual deepening of seismicity towards the western profiles, which is also highlighted in Figure 4a 1 -a 2 , as well as in the N315 • E-oriented cross-sections of Figure S14 and the N100 • E-oriented cross-section of Figure S15. However, the focal mechanisms associated with the larger offshore events of group #4 (e.g., the 21 September 2014 M w = 4.6 earthquake) suggest an E-W to WSW-ENE-striking fault. The latter, along with the hypocentral distribution, could indicate an association with the Ano Kastritsi fault, further south. The two available focal mechanisms for events below Platanitis support SW-NE faulting, with one being normal and the other dextral, both compatible with the kinematics of the Rion-Patras fault zone.
The N10 • E cross-section b 1 -b 2 of Figure 4 contains seismicity of groups #7 (yellow) and #8 (brown), the latter associated with the June 2014 outburst in the gulf. Although the seismic activity is mostly compacted in the depth range 8-10 km, distinct clusters can be associated with individual structures, some being located at shallower depths (6-8 km) towards the north. The 8 June 2014 M w = 4.3 event was likely hosted on a north-dipping structure, with a~45 • dip according to the focal mechanism reported by NKUA (Table 1). This makes it partly compatible with the Fasouleika fault segment, if a small degree of listricity is considered. Likewise, the different strike (NW-SE) of the 24/01/2014 M w = 3.9 event could support an association with the Selianitika or Lambiri faults. Seismicity at the northern part of the group #8 is distributed on a gently dipping layer, with a dip of~23 • . This is sub-parallel to the seismicity of group #7, further north, which seems to spread deeper and northwards along the same trend (Figure 4b 1 -b 2 , dashed yellow line with a northwards dip of~27 • ). Shallower seismic activity in group #7 could be related with antithetic, south-dipping structures, such as the Sergoula fault which is associated with the eastern part of the Trikorfo seismogenic source (GRCS441 in [14]; Figure S20). Colours and numbers correspond to the 9 spatial groups into which the dataset was divided, as well as the "miscellaneous" group #10 (gray). Beachballs depict the focal mechanisms of major events (Table 1) Another important profile is presented in the cross-section c 1 -c 2 of Figure 4, which includes groups #6 and #9. The 2013 Helike swarm (group #9, blue) occurred mainly in the depth range of 8-10 km. A detailed analysis of this swarm [6] has revealed focal mechanisms associated with both north-dipping normal faulting and oblique-normal faulting. Such kinematics have also been observed for the 2001 Agios Ioannis swarm [4,8,49], associated with the northwest-dipping Kerinitis fault. The extension of the north-dipping fault plane of the 2013 Helike swarm to the surface is traced near the junction between Pirgaki and Kerinitis faults [6]. The relation of these faults to the structures activated during the 2013 Helike swarm is presented in Figure S25 (see also the provided interactive Matlab figure; File S1). It should be noted that the Pirgaki fault is likely a branch of the Mamousia seismogenic source (GRCS524 in [14]; Figure S20).
The southern seismicity in group #6 (red) is mostly part of the aftershock sequence of the 7 November 2014 M w = 4.8 event and appears to be associated with the west Helike fault. The focal mechanism of the mainshock suggests a northwards dip of 30 • , while the hypocentral distribution of its aftershocks (south of its epicenter) corresponds to a roughly~40 • dip. The latter could indicate some degree of listricity, or the existence of multiple smaller activated structures within the seismogenic layer, each dipping at a higher angle (e.g., 60 • ). Seismicity of group #6 extends further north, with the focal mechanism of a smaller event that occurred on 28 January 2013 having a dip of 28 • . This low dip angle is similar to that of the northern part of group #8 (brown), as well as of the deeper part of group #7 (yellow). It is noteworthy, however, that when examined in a N100 • E cross-section ( Figure S15, along the Kamarai fault zone), there appears to be a vertical offset between these shallow dipping segments of groups #4, #8 and #6, tending to deepen from west to east (see also File S1). This is likely related to the en echelon faulting configuration at the southern shoulder of the WGoC.
Other parts of seismicity ( Figure S13) include the northern groups #5 (orange, between Nafpaktos and Pyrgos, i.e., station PYRG) and #2 (magenta, near Eratini; Figure 1). Both are associated with sparse seismicity on a shallow-angle, north-dipping plane, reaching depths larger than 10 km below the northern coast. Group #1 (black), likely related with the deeper and almost flat sector of the South Corinth Gulf composite seismogenic source (GRCS500 in [14]; Figure S20), follows the trend of group #2 southwards at depths of 8-12 km. At the eastern margin of the study area, the isolated cluster #3 (cyan) is associated with a nearly dip-slip faulting type ( Figure 4, Table 1) on a sub-vertical WNW-ESE-trending fault, slightly tilted northwards. Seismicity in the miscellaneous group #10 (gray) includes activity around the margins of the study area or deeper than 15 km. The latter is the case for events southwest of group #4, in the form of small clusters located at depths of~17-22 km. Some deeper events are also located below Aigion, including an M w = 4.0 earthquake that occurred on 28 April 2013 at a depth of 47 km with its epicenter offshore, east of group #9.

Spatio-Temporal Analysis of Seismicity
The temporal evolution of seismicity in the WGoC during 2013-2014 is examined in terms of seismicity rates, spatiotemporal seismicity migration speeds and pore-pressure diffusion. Figure 5 (top) shows the cumulative number of events per spatial group, for selected seismicity above the magnitude of completeness, which is estimated at M c = 1.8 for the manual dataset ( Figure S19). The application of this threshold enhances the occurrence of seismicity bursts, usually triggered by events of larger magnitude, related to specific sequences whose initiation is marked by vertical gray lines in Figure 5. Some spatial groups are characterized by generally constant background seismicity rates, i.e.,~0.1 events/day in group #1, 0.2 events/day in group #2 and 0.4-0.5 events/day in group #5. The miscellaneous group #10 (gray) presents a background seismicity rate of 0.4-0.7 events/day, increased to~1.1 events/day in the second half of 2014. In other groups, their low seismicity rates, excluding periods of bursts, range between ≤0.1 and 0.6 events/day. A general sense of the spatial distribution of background seismicity in the WGoC can be observed during January-August 2013 ( Figure S22a). During that period, the more significant cluster was the 2013 Helike swarm in the south, whereas only sparse seismicity and small clusters are distributed in the rest of the study area. Figure 5 (bottom) shows the projection of epicenters along the N100 • E profile line A-B (Figure 4, top) with respect to the origin time, without filtering out events with M w < M c . In early 2013, most seismic activity was concentrated mid-gulf at group #6 (red), following an M w = 3.8 event on 28 January 2013, and at group #8 (brown) with a few minor bursts. Groups #2 (magenta), #4 (green), #5 (orange) and #7 (yellow) exhibited relatively constant seismicity rates of 0.2-0.4 events/day. Then, on 21 May 2013, the first phase of the Helike swarm (group #9, blue) was initiated with a series of events reaching M w 3.6-3.7, mainly activating the eastern half of the cluster. The activity in Helike culminated with a second phase during July, involving an M w = 3.8 event on 15 July 2013, which triggered seismicity at the western part of the group. The first two phases of the swarm, until August 2013, are described in detail by [6]. A third phase began on 26 October 2013, lasted for about 20 days and was concentrated at the mid-northern part of the cluster. Group #7, near Marathias, presented a small burst of activity in early September 2013.
An increase in the seismicity rate to~1.1 events/day began in October 2013 at the offshore group #4, in the region between Nafpaktos and Psathopyrgos, presenting bursts associated with the occurrence of larger events (M w ≥ 3.0). This first wave of seismic activity in group #4 lasted up to February 2014, with its seismicity rate returning to 0.6 events/day. Seismicity of group #4 was initially concentrated in the Y-range of 10-20 km ( Figure 5, bottom), but gradually spread bilaterally towards Y = 25 km in January and Y = 5 km in February 2014. At that time, the seismic activity was also increased, temporarily, in the neighboring groups #7 (yellow) and #8 (brown). The isolated cluster #3 (cyan) near Itea presented swarm-like activity between 8 March and 14 September 2014, following a series of events, the largest being an M w = 4.0 on 10 May 2014.
On 8 June 2014, 15:10:51 UTC, an M w = 4.3 event initiated an aftershock sequence in group #8, at the offshore region SW of the Trizonia Island, lasting for about one month. Increase in the seismicity rate of group #6 (red), including several M w ≥ 3.0 events, began in July 2014 ( Figure A1). Part of that group was triggered by the 8 June 2014 event in group #8 ( Figure A2), likely due to the partial overlap between the involved structures, as these groups cannot be fully distinguished with spatial margins. More importantly, a significant pattern of bilateral seismicity migration began to emerge in mid-July 2014 in group #4, originating at Y ≈ 15 km ( Figure 5, bottom), spreading towards Y = 25 km and Y = 10 km in a period of two months. Seismic activity at the western part of group #4 was abruptly increased in the form of an aftershock sequence after the occurrence of an M w = 4.6 event on 21 September 2014 [48]. The last major seismicity outbreak at the WGoC during the study period began on 7 November 2014, with an M w = 4.8 event in group #6 (red), offshore Aigion [48]. The seismic activity spread rapidly in the form of an aftershock sequence, bilaterally along a~7-km-long patch in an E-W direction ( Figure A3), but mostly south of the mainshock's epicenter ( Figure S22c).
Regarding the physical mechanisms behind the spatiotemporal seismicity migration patterns, stress-transfer and pore-pressure diffusion are considered. Slip occurs when the applied shear stress on a fault plane, aided by the pore-pressure in fluid-filled microcracks, surpasses the sum of the cohesive strength of the material and the internal friction. Changes in the level of stress in a specific location may result from Coulomb stress transfer due to the occurrence of an earthquake in a nearby region, which is the main mechanism behind the triggering of aftershocks. Aseismic slip, e.g., due to creeping, also causes Coulomb stress changes that may load a neighboring area, accelerating earthquake occurrence, or relaxing another region. Finally, the diffusion of pressurized fluids can generate earthquakes in the volume behind a triggering front that spreads through the fracture network, increasing pore-pressure and lowering the effective normal stress on a fault plane, thus facilitating slip. When fluids are injected in a homogeneous and isotropic, poro-elastic, fluid-saturated medium, the expanding spherical triggering front of induced seismicity follows Equation (1) [50]: where r(t) is the triggering front's radius (distance from the injection point), t is the time since the beginning of the pore-pressure perturbation at time t o and D is the hydraulic diffusivity. The same principle may apply as a first approximation to a fluid-saturated crust and an initial pore-pressure perturbation. The determined spatiotemporal seismicity migration rates ( Figure 5, bottom) provide a first indication regarding the possible physical triggering mechanism of seismicity in the WGoC. Initiating in September 2013, offshore Nafpaktos, seismicity started to migrate eastwards at a speed of~0.13 km/d and then slowed down to 0.04-0.05 km/d after February 2014. Figure 6 shows the spatiotemporal evolution of the offshore groups #4, #8 and #6 (from west to east), mostly confined below the bold parabolic envelope drawn after Equation (1) with D = 2.2 m 2 /s. This triggering front starts in September 2013 at the Y ≈ 12 km position along the A-B profile line of Figure 4 (top). As the fluids propagate eastwards, earthquake sequences are triggered in bursts near the expanding front. This suggests a combination between pore-fluid pressure diffusion through the fracture network and seismic swarms occurring along the way. The latter evolve with additional stress transfer due to tectonic slip. In natural fluid injections, the time, location and geometry of the injection source can only be deduced by the observed seismicity migration patterns. In that sense, secondary triggering front envelopes can be inferred for each sub-sequence, e.g., at Y ≈ 15 km in mid-December 2013 or near Y ≈ 25 km in January-February 2014, with D = 0.5 m 2 /s. An eastwards-propagating patch with a width of~10 km seems to be active behind the D = 2.2 m 2 /s triggering front. However, a deficit of seismicity is observed west of that patch ( Figure 6). The latter indicates that the initial fluid injection was likely terminated and a back-front r b (t) is formed according to the relation [51]: where t d is the duration of the injection and D b the hydraulic diffusivity corresponding to the back-front. In Figure 6, a back-front curve of D b = 2.2 m 2 /s is shown with a bold red line, beginning~120days after the initiation of the triggering front at point Y ≈ 12 km on 1 September 2013. In that case, the post-injection hydraulic diffusivity is the same, considering pore-pressure diffusion in three dimensions (Equation (2); [51]). Meanwhile, in July 2014 the Nafpaktos-Psathopyrgos swarm began to evolve, starting very close to the supposed injection site at Y ≈ 15 km (dashed red ellipse in Figure S22b,c).
However, its evolution cannot be fit with a parabolic envelope, but is better approximated by a roughly linear triggering front, up to a point when it even tends to accelerate as it spreads ( Figure A4). Its seismicity migration speed was 0.13 km/d eastwards and 0.03-0.07 km/d westwards between 15 July and~20 August 2014 and then the eastwards migration speed increased to 0.17 km/d. The 21 September 2014 event caused seismicity to spread eastwards at over 0.35 km/d, while the western branch began to expand rapidly at over 1 km/d. The latter speed is likely too fast to be associated with pressurized fluids in a tectonic environment and is more akin to triggering by aseismic slip (e.g., [52]).
We also estimated the mean squared distance (msd) x 2 (t) of seismicity with time t, which provides a proxy for the rate of triggered earthquake migration [53]. In normal (Gaussian) diffusion, the msd scales linearly with t. Instead, the hallmark of anomalous diffusion is the non-linear scaling of the msd with t, that frequently takes the form of a power-law (e.g., [54]): The power-law, or diffusion exponent a, characterizes the various domains of the diffusive process, designating sub-diffusion for 0 < a < 1 and super-diffusion for a > 1, while for a = 1 normal diffusion is recovered.
Taking as an approximate reference (x 0 = 0, t 0 = 0) the Y ≈ 12.3 km position along the A-B profile line of Figure 4 (top) on 1 September 2013 (similar to the reference point of the triggering front shown in Figure 6), we estimate the msd of seismicity with time. This is performed for the events with M w ≥ M c , enclosed within the parabolic envelope defined by the triggering-and back-fronts of seismicity for D = 2.2 m 2 /s ( Figure 6). The msd is calculated as ( [49]): where x n (t) is the distance of the n th event from the reference point at time t and N is the total number of events. The result is shown in Figure 7 in double-logarithmic axes and logarithmically spaced bins. The overall growth of the msd with t signifies the diffusive process. This growth is not continuous, but rather evolves in sequential "steps", as observed after the passage of~50, 130, 272 and 418 days (Figure 7). These days mark the activation of the sequential seismic clusters in groups #4, #6 and #8 towards the east, within the parabolic envelope defined in Figure 6. On the other hand, during the in-between time periods, the msd is almost constant or grows slowly with t. The overall scaling behavior of the msd with t can be approximated with Equation (3), yielding a diffusion exponent a = 1.81 ± 0.17 (R 2 = 0.98). In terms of the Continuous Time Random Walk (CTRW) model (e.g., [55]), this value of the power-law exponent indicates super-diffusion. Compared to the range of exponents reported in the literature that mainly correspond to slow earthquake sub-diffusion [49,[56][57][58][59][60], the diffusion exponent in this case is significantly larger and comparable to that reported by [61]. Hence, the large diffusion exponent likely corresponds to large rates of earthquake migration triggered by aseismic factors, such as pore-pressure diffusion or aseismic creep.

Shear-Wave Splitting Parameters
The distribution of the fast shear-wave polarization direction (ϕ) measurements from Cat1 for selected stations, representing sites with different average ϕ directions, is presented in Figure 8 in the form of rose diagrams. The respective diagrams for the rest of the stations are available in Figure S16. Rose diagrams, widely used in shear-wave splitting studies (e.g., [22]), are circular histogram plots which display directional data (in the current case, fast shear-wave polarization direction data) and the frequency of each class, (i.e., the number of measurements for each ϕ bin). The SWS parameter statistics for Cat1 are presented in Table 2. The average ϕ value is calculated with circular statistics [62] to avoid wrapping between 0 • and 180 • . In most stations, the mean ϕ value coincides with the direction that corresponds to the larger amount of valid measurements. There are, however, a few exceptions, such as station HELI ( Figure S16), which presents two main fast-shearwave polarization directions, one at~N30 • E and another at~N100 • E. The mean ϕ in that case was calculated at 77.6 • ± 11.6 • . However, this value does not correspond to either of the two maxima of the distribution. There are also some stations with significant ϕ scatter, such as VVK, and certain stations of the Magoula array, i.e., MG03, MG05 and MG07, for which the mean value differs from the azimuth that corresponds to the maximum of measurements. The ϕ value of valid measurements for all stations (Cat1) is 98.8 • ± 2.8 • . Important deviations from this direction are observed for stations LAKA (ϕ = 134.7 • ± 5.9 • ) and TEME (ϕ = 145.1 • ±16.3 • ), in the south, and stations SERG (ϕ = 58.4 • ±7.5 • ) and MALA (ϕ = 68.1 • ± 11.0 • ), in the north. Stations with significant ϕ scattering, indicated by the large uncertainty, include MG03 (ϕ = 67.1 • ± 30.8 • ), PANR (ϕ = 132.8 • ± 36.6 • ) and PYRG (ϕ = 142.2 • ± 41.6 • ). The time-delays, t d , are partially depended on the spatial distribution of the seismic sources relative to the receiver. Since the time-delay is accumulated along the path-length of the seismic ray as it travels through the anisotropic medium, longer ray-paths will result in larger t d values. The time-delays are normalized (t n ) by dividing t d with the hypocentral distance to remove the effect of the ray-path length. The mean t n values are overall (9.2 ± 0.1) ms/km, while average t n per station ranges from 5.0 ms/km (AGRP) to 12.5 ms/km (ALIK). The raw t d measurements have a mean of (112.3 ± 1.4) ms.  Figure S16 for rose diagrams of the other stations. Table 2. Mean values of SWS parameters, i.e., mean polarization direction of the fast shear-wave, ϕ, mean time-delay between the two split shear-waves, t d , and mean normalized time-delay, t n , along with their respective uncertainties per station for valid measurements (grade A-C) of Cat1. Errors δϕ correspond to the 95% c.i. of the mean value, whereas δt d and δt n are the respective typical errors of the mean. N is the number of valid measurements per station. The first line in bold corresponds to all valid measurements of Cat1. The last column, "OC", indicates the orientation correction applied on some stations, given in the form of the true azimuth of the instrument's N-S component. The distribution of ϕ values for null measurements, with ϕ being within ±10 • from p or p + 90 • , i.e., sub-parallel or perpendicular to the polarization direction of the corrected (for anisotropy) shear-wave, p, is presented in Figure S17. The distribution for all stations indicates a ϕ null = 5.4 • ± 6.1 • , i.e., in an approximate N-S direction, that is perpendicular to the respective ϕ of the valid measurements. This suggests that the majority of null ϕ values likely correspond to the polarization direction of the S slow wave. However, there are (fewer) cases where ϕ null measurements of certain stations are similar with their ϕ, e.g., SERG (ϕ null = 64.9 • ± 12.1 • ). In most stations, the average t n for null measurements is larger than the respective for valid ones.
Regarding the temporal variations of the SWS parameters, we take advantage of the large number of measurements from the combination of the manual and automatic catalogues to investigate patterns of change with time. The time-delay between the S fast and S slow waves, normalized by the ray-path length, can be used as an indicator for the degree of anisotropy in the medium. However, factors besides the traveled distance also affect t d . These include the difficulties in the precise measurement of ϕ or t d related to the radiation pattern (focal mechanism), the distortion of polarization due to non-perpendicularity of the incident seismic rays to the surface, variations in the frequency content between S fast and S slow and larger attenuation of the S slow [44,63], effectively contributing to a significant degree of scatter. To exclude measurements that are less affected by changes in the geometry of the fluid-filled microcracks, we select event-station pairs whose ray-paths belong to "Band-1". This is defined as the double-leafed solid angle between 15 • and 45 • on either side of the vertical crack-plane (whose azimuth is defined by ϕ) passing through the station (see Figure S21 for a schematic view of Band-1). Ray-paths of event-station pairs travelling through the ±15 • solid angle from the crack-plane are classified as "Band-2" measurements and are more affected by the microcrack density [64]. Optional filters are applied to remove part of the t n measurements, such as those affected by a higher angle of incidence at the surface, having strong spatial dependence or having a lower quality grade. In addition, temporal variations of t n are also examined after discarding results related to ϕ measurements that significantly differ from the ϕ of that station measured from Cat1, which may be potentially erroneous. For events of Cat2 (automatic), we applied a focal depth correction, as their lower quality hypocenters tended to be scattered at much shallower depths than their neighboring events of Cat1, causing inconsistencies in the t n dataset of Cat3. This was performed by applying, for each event of Cat2, the average focal depth of the events of Cat1 with neighboring epicenters (excluding those with Z > 15 km). This modification reduces some outliers in Cat2 that yielded erroneously high t n values due to the shorter distance (shallower depth). Small fluctuations of ±1 km or even ±2 km from that average depth don't seem to affect the measured t n significantly, as the resulting changes are much smaller than the observed measurement scatter.
We focus on patterns of significant linear increase or decrease of t n values observed at more than one station during specific periods, as this may suggest a change that could be related to more widespread stress variations. To assess this task, we performed a series of linear regressions on sliding windows along the t n data (with or without the application of the above-mentioned filters) versus origin time. The regressions were applied using a range of different window-lengths, measured either in number of successive data points or days, to allow for various time scales and data densities to be taken into account. For each window, least-squares inversion is performed to calculate the best-fit line, along with the respective correlation coefficient, R 2 . On each case, the Student's t-test is performed to examine whether the null hypothesis of zero slope can be rejected with a 95% confidence interval, and cases with p > 0.05 are discarded. We then examined the results and identified patterns of increase or decrease, which either persist or emerge, as stricter filters are applied. Furthermore, we examined periods during which a group of stations presented the same behavior regarding the t n change, whereas another group exhibited the opposite change behavior. The 3-point moving average is also calculated, to account for non-linear temporal variations or weak correlations. Figure 9 shows the t n temporal variations for selected stations. In all cases, Band-1 was selected, while filtering out measurements with ϕ outside the ϕ ± 30 • range (with ϕ as determined from Cat1) usually manages to further reduce scattering. Certain common characteristics are identified between different stations. In January-April 2013, stations PSAR, KALE and ZIRI exhibit a moderate decrease of t n values. Weak decrease is also observed at stations TRIZ and PYRG, whereas station SERG has increasing t n values during February-May 2013. That is the period before the occurrence of the 2013 Helike swarm. Station TRIZ presents a mild increase of t n in the period of May-July 2013, while timedelays at PYRG and ALIK (with events from groups #6 and #8 only; Figure 9) also follow a weakly increasing trend that lasts until December 2014. On the other hand, in October 2013-January 2014, several stations in the north present a significant decrease of t n values, including PSAR, KALE, TRIZ and most of the stations of the Magoula array, e.g., MG04 (Figure 9). One of the trend lines of station PYRG also indicates a strong decrease of t n between November 2013 and March 2014. This period could be of particular significance, as it coincides with the time when the proposed seismicity triggering front of Figure 6 passes near the middle of profile A-B (Y ≈ 25 km). On the contrary, some stations in the south, such as AIOA and ALIK (the latter for events of groups #6 and #8), exhibit an increasing pattern in t n values during that period. February 2014 marks the beginning of new trends for many stations, with increasing t n values for PYRG and ZIRI and decreasing ones for AIOA and TRIZ. Lastly, a decreasing t n trend is exhibited at several stations, including ZIRI, KALE, AGRP and SERG, towards the last quarter of the study period; a time-span that begins with the M w = 4.2 event of 8 June 2014 and is marked by the seismic swarm in the Nafpaktos-Psathopyrgos area. A shorter decrease during August-September 2014 is observed at station ALIK (groups #6 and #8), following a longer period of moderate increase, while longer decreasing trends are observed at stations AIOA, TRIZ, beginning around February 2014. The latter trends last until September-November 2014, when the 21 September and the 7 November events occurred at the western and central parts of the study region, respectively. This decrease of t n values might be related to crack coalescence that takes place before the stress reaches the critical point where the M w = 4.8 earthquake occurs.

Spatial Distribution of Seismicity and Possible Active Faults
The high-resolution hypocentral distribution of the present study reveals that the seismic activity is mostly concentrated in a~2-km-thick, low-angle-dipping layer below the gulf, mainly between depths of 7 and 10 km. This is composed of distinct clusters associated with planar structures protruding towards shallower depths (see also Figure A5 and the interactive figure of File S1 in the Supplementary Material). Focal mechanisms ( Figure 4 and Table 1) indicate that these clusters are related with normal faults of moderate dip angle. A connection of structures at depth to mapped faults at the surface is not always straightforward, either due to the activation of blind faults, rooting at the fragmented seismogenic layer, or to focal depth uncertainties and biases.
The latter case is important, as some systematic differences at hypocentral depths can be observed when compared to other studies. For example, [9] presented relocated seismicity of the period 2000-2015 in the WGoC using data from the CRL network, with focal depths mainly in the range of 6-8 km. Although in the northern part the focal depths of Cat1 are similar to those of [9], a deviation is observed in the southern and western parts of the study area, with absolute locations of our catalogue being 2-4 km deeper. This is likely a bias introduced by the inclusion of HUSN stations at epicentral distances further than 30 km, to better constrain epicenters near the margins of the local network, particularly towards the western part. The use of stations at distances up to 30 km tends to pull hypocenters towards shallower depths, especially at the southern part below Aigion, but also increases scattering in the distribution of epicenters. As a result, the foci in the present study are more clustered, but likely at the cost of a 1-2 km offset in absolute depths. The 2013 Helike swarm is an exception to this, as it has been separately processed using similar data, stations and velocity model as in [6], but considering it as a whole spatial group rather than on a sub-cluster level.
Concerning the larger events of the dataset, the 8 June 2014 event exhibits nearly E-W normal faulting (Table 1). [9] related this event with the south-dipping Trizonia fault. However, as observed in cross-section b 1 -b 2 (Figure 4), this cluster of aftershocks seems to be consistent with a north-dipping structure that is likely associated with the Aigion-Selianitika fault system (GRIS503 in [14]; Figure S20) in group #8. As highlighted by [9], this earthquake seems to have also triggered seismicity to antithetic structures, sub-parallel to the Trizonia fault. Furthermore, its activity likely spreads eastwards to the Fasouleika segment, in group #6, particularly by the end of July 2014, in an opposite direction to what was observed in 2003-2004 [5].
Another significant event during the study period is the 21 September 2014 earthquake, which has been estimated at M w = 4.6 (M w = 4.9 from MT at NKUA and GI-NOA; Table 1). This event triggered aftershocks at the westernmost part of the 2014 Nafpaktos-Psathopyrgos swarm (group #4). Activity in this part of the WGoC has been related with the junction between the Rion-Patras fault zone, dipping NW and having a dextral strike-slip component, and the Psathopyrgos north-dipping normal fault [9]. The focal mechanism for this event suggests E-W to WSW-ENE normal faulting (strike = 254 • , dip = 40 • ). Although the cross-section a 1 -a 2 in Figure 4 is oblique to the focal mechanism, the hypocentral distribution shows a better fit with a NW-dipping plane. The 21 September 2014 event seems to coincide with the down-dip extension of the Ano Kastritsi fault. The triggered aftershocks, however, could have been spread to another structure, such as the Panachaikon transferfault (Pt.f., Figures 1 and 4; [30]), which is parallel to the Rion-Patras fault but further east. When surface faults and deeper structures associated with the seismicity of group #4 are viewed in the 3D space ( Figure S24, File S1 in the Supplementary Material), it is observed that the down-dip extension of the Rion-Patras fault and Panachaikon transfer fault, with a typical 60 • dip, intersect the hypocentral distribution. If the seismicity cloud is interpreted as a combination of two parallel NW-dipping structures at a high angle (as in Figure S24), it is plausible to associate it with those two faults on the surface. The 7 November 2014 event was previously referred to as having M w = 5.0 (e.g., [48]), following the MT calculated at NKUA. However, spectral fitting provided M w = 4.8 for this event, which is closer to the value reported by GI-NOA MT inversion (M w = 4.8) and [9] (M w = 4.7). This event is relocated within group #6, with a spatial gap to its south, followed by a cluster of aftershocks. The hypocentral distribution of this group suggests activation of a north-dipping structure at a relatively low angle (25 • for the north-dipping plane of the focal mechanism). The southern sub-clusters of group #6, taken as a whole, have an apparent dip of~40 • . Allowing for a~2 km offset at depth and some degree of listricity, the activated structure is compatible with the west Helike fault. [9] related the activity at the southern part of group #6 with the Aigion fault, as their hypocenters reach depths shallower than 5 km. Considering that group #6 is composed of high-angle structures, when viewed in the 3D space ( Figure S23, File S1 in the Supplementary Material) it seems plausible to partly relate the seismic activity of this group with the down-dip extensions of the w. Helike fault, the Aigion fault and the Fasouleika segment.

Physical Mechanisms for Earthquake Migration
The temporal evolution of seismicity in the WGoC during 2013-2014 presents a series of bursts with significant earthquake migration patterns. The most characteristic one was observed in the offshore region between Nafpaktos and Psathopyrgos (group #4), mainly during July-October 2014, with migration speed in the range of 0.13-0.18 km/day eastwards and~0.07-0.12 km/day westwards ( Figure S22c). The occurrence of the 21 September 2014 earthquake in the west produced an aftershock sequence that severely accelerated earthquake migration rates to over 1 km/day in its vicinity. Seismicity migration speeds of the order of 100 m/day are associated with fluid-driven seismicity, as observed in several cases in the WGoC (e.g., [4][5][6]). The speed of the triggering front, due to pore-pressure diffusion, depends on the distance from the fluid injection source (Equation (1)) and from the permeability k of the rock, embedded in D = Nk/η, with other factors being the poro-elastic modulus, N, and the pore-fluid dynamic viscosity, η [50]. On the other hand, earthquake migration rates of over 1 km/day (e.g., 0.1-1 km/h) are likely attributed to aseismic creep episodes (e.g., [52]). [65] have documented the case of a seismic swarm in Malamata, i.e., below station MALA (Figure 1), which exhibited a pattern of expanding seismicity at a rate of~125 m/day. However, it was manifested as a series of intermittent bursts, with seismicity migration speeds of 2-10 km/day over a short duration and distance. This example indicates how fluid-driven seismicity can be combined with short-lived aseismic slip episodes.
We presented a long-range earthquake migration pattern, covering a~32 km distance from Nafpaktos to Aigion over the course of~15 months ( Figure 6). Although this sequence evolved eastwards at an average speed of~70 m/day, it progressed in a series of bursts involving events of moderate magnitude that triggered their own sub-sequences. This pattern is consistent with pore-fluid diffusion, evidenced by an expanding triggering front of hydraulic diffusivity D = 2.2 m 2 /s through an active patch of 11-12 km, delimited by a back-front of seismicity deficit ( Figure 6). A similar pattern has been observed in 2011 in Oichalia, South Peloponnese, where, likewise, a swarm evolved southwards in a series of episodes associated with major events of NNW-SSE-oriented normal faults, within an active patch following a parabolic triggering front [66]. In that case, the hydraulic diffusivity was D = 1.1 m 2 /s, the active fault patch was~8-km-long and it traveled a distance of~10 km over 160 days, i.e., at an average speed of 63 m/day.
In addition, we estimated the mean squared distance (msd) of seismicity with time for the eastwards expanding seismic activity, enclosed within the active patch defined by the triggering-and back-fronts. The msd shows a step-like increase with time that can be approximated, overall, by a power-law function, with a diffusion exponent of a = 1.81 ± 0.17 that indicates super-diffusion. [60] studied separately the distinct seismic clusters occurring in the spatial groups #4, #6 and #8 (e.g., Figures A1, A3 and A4 in Appendix A), yielding diffusion exponents lower than 0.2 and sub-diffusion. However, herein we observe that if the seismic clusters are jointly studied on a larger scale, a much faster diffusive process is obtained that spans from the Rion strait in the west to offshore Aigion in the east (Figure 4, top). In terms of the CTRW theory, the overall super-diffusive process resembles a Lévy walk (e.g., [67]), with intermediate asymptotics having different exponents. High diffusion exponents, close to or greater than unity, are found in anomalous fluid transport processes in heterogeneous and fractured media [68]. Numerical simulations of fluid transport in heterogeneous fault zones also show anomalous diffusion with exponents greater than 1, as fluid flow is channeled in fractal fault zones with varying permeabilities [69]. In contrast, stress transfer effects due to subcritical crack growth or to cascading effects yield quite low diffusion exponents and sub-diffusion [57,58], most commonly observed in earthquake diffusion [49,60]. The super-diffusive behavior of seismicity reported herein may, thus, provide an additional indication of fluid overpressures.
Consequently, a candidate physical mechanism for the observed seismicity migration pattern in the offshore region of the WGoC during the~15 months period may well correspond to an initial intrusion of overpressurized fluids from the area of Nafpaktos. In a fluid-saturated and critically stressed crust, fluid overpressures can reduce the effective normal stress along the active fault structures, triggering the observed seismicity offshore Nafpaktos during September-December 2013 ( Figure S22b). Fluid diffusion along an E-Woriented, permeable pathway, defined by the intersection of the active fault structures with the shallow-dipping seismogenic layer beneath the WGoC, may have caused the eastward migration of seismicity. This eventually activated the easternmost part of group #4 and the westernmost part of group #8 during January-February 2014, which included certain M w > 3.8 events ( Figure S22b, yellow colour). Increased stresses in this area trigger additional events towards the east during March-May 2014, until 8 June 2014, when the M w = 4.3 earthquake occurred in group #8. The latter initiated an aftershock sequence, including several M w ≥ 3.0 events in group #8 and eastwards in group #6 during July-October 2014. Seismic activity then culminated with an M w = 4.8 event on 7 November 2014 in group #6, offshore Aigion, initiating an aftershock sequence that spread rapidly in an E-W direction and mostly south of the mainshock's epicenter until the end of 2014 ( Figure S22c).
The overall seismicity migration pattern is, thus, likely produced by a combination of fluid diffusion and stress transfer effects caused by the larger magnitude events. Enhanced permeabilities due to fluid flow in the area of group #4, followed by healing/sealing processes in the active faults (e.g., [70]), may have caused the deficit of the seismic activity, delimited by the back-front observed during March-June 2014 ( Figure 6). In mid-July 2014, the seismic activity at group #4 increases again in the same area, as during September-December 2013, possibly triggered by another intrusion of overpressurized fluids. This resembles a conventional fault-valve behavior [71], in which the low-permeability Phyllade nappe below the shallow-dipping seismogenic layer acts as a cap for a deeper fluid reservoir evidenced by V P /V S anomalies [5,65,72,73]. A combination of fluid diffusion, aseismic creep, supported by fast earthquake migration speeds, and stress transfer effects caused by the M w > 4 events during September 2014, may have caused the bilateral migration of seismicity that was observed offshore Nafpaktos during July-October 2014.

Seismic Anisotropy in the WGoC
The comprehensive dataset of shear-wave splitting parameters of the present study may provide evidence regarding changes in the local stress field, fluid diffusion or saturation that can affect the anisotropic properties of the medium. Previous upper-crust anisotropy studies in the WGoC have focused in the 2010 Efpalio sequence [74], the 2013 Helike swarm [75] and the November 2014 sequence [48]. Their findings indicate that the fast shear-wave polarization directions are consistent with the maximum horizontal stress component, σ Hmax (e.g., [76]), which is in a WNW-ESE direction, parallel to the strike of major normal faults. Herein, the obtained mean ϕ directions are collectively presented in Figure 10. Notably, the results of the automatic analysis are, on average, consistent with those of [48,75], who performed manual SWS measurements through visual inspection of polarigrams and hodograms. In the following, the presented average values of the SWS parameters are from measurements of Cat1, unless specified otherwise. These measurements concern the volume between the sources and the receiver as sampled by the seismic ray-paths, whose length is used for the normalization of the time-delays. In addition, we provide SWS measurements for stations in the WGoC that were not analyzed previously, such as AIOA, HELI and TEME in the south, as well as PANR and PSAM (Cat3 only) in the north. Station TEME presents a NW-SE ϕ direction (145.1 • ± 16.3 • ), which does not appear to be consistent with σ Hmax or the strike of local structures, but is similar to that of LAKA (ϕ = 134.7 • ± 5.9 • ), to the west. On the other hand, HELI, a temporary station installed to record the 2013 Helike sequence [6], exhibits a WSW-ENE ϕ direction (77.6 • ± 11.6 • ) that is similar to that of the neighboring AIOA station to the south (ϕ = 79.4 • ± 7.7 • ). However, the ϕ value at station HELI deviates from the direction indicated by the largest percentage of ϕ measurements, which is~E-W, as can be observed from the respective rose diagram ( Figure S16), with a secondary SSW-NNE ϕ direction. Stations AIOA and HELI are related with the epicentral regions of the 2001 Agios Ioannis [4,49] and 2013 Helike [6,48] swarms, both being either directly or partially related with the re-activation of older structures, which are not optimally oriented to the present-day stress-field. Specifically, the 2001 Agios Ioannis swarm has been attributed to the Kerinitis fault (Figure 1), with a strike that is SW-NE to WSW-ENE, consistent with ϕ measured at AIOA and HELI. The ϕ direction measured at TEME, however, is almost perpendicular to that of its neighboring station HELI, although the stations are located in close vicinity (~2 km apart).
In the north, the ϕ directions at stations TRIZ, PSAR and PSAM are nearly E-W, consistent with both σ Hmax and fault orientations. An important difference was previously observed at stations PYRG and SERG, with ϕ oriented SW-NE, a deviation that was attributed to the strike of local structures [48,75], particularly to mapped offshore faults with significant strike-slip component (e.g., [17]). In the present study, however, station PYRG exhibits a much different NW-SE ϕ direction in Cat1 (142.2 • ), albeit with a very significant error (±41.6 • ). On the other hand, in Cat3, where more measurements are available, station PYRG shows a more E-W-oriented ϕ, although with a still large uncertainty (102.4 • ± 29.6 • ). This is due to the significant scatter in ϕ measurements (Figure S16), while a secondary ϕ maximum at a SW-NE direction, that would be in agreement with previous studies, is also observed. The existence of two groups of ϕ measurements at stations PYRG and HELI indicates either a strong spatial dependence or that shear-waves propagate through two distinct media with different anisotropic properties [21,77,78].
The measurements at the neighboring station SERG are very consistent, with ϕ = 58.4 • ± 7.5 • . Further west, stations EFP and MALA present an E-W to WSW-ENE ϕ direction, consistent with previous findings [48,75]. Regarding the Magoula seismic array, station MG00, with ϕ = 82.2 • ± 7.5 • , is in agreement with the results of [48]. Other stations of the array present variability in their ϕ values, ranging between 60.7 • and 106.3 • (Figure S16). Station MG06, however, with more consistent measurements and smaller error, yields ϕ = 85.6 • ± 8.9 • , in general agreement with MG00. The WSW-ENE ϕ orientation at station MALA (68.1 • ± 11.0 • ) appears as a smooth transition in relation to the respective ϕ at the neighboring stations EFP and MG00, and is also consistent with the strike of the Nafpaktos fault. Observations of ϕ at stations MALA and SERG favor the existence of structurally controlled anisotropy [41,79,80]. Further west, stations AGRP and VVK present an E-W to WNW-ESE ϕ direction, more consistent with σ Hmax [48,75]. Notably, in the WGoC the mapped faults' strike changes from E-W in the north to WNW-ESE in the south, the latter being in agreement with ϕ measured at stations ROD3, ZIRI and ALIK. Stations TEME and LAKA show a more pronounced NW-SE direction, consistent with the results of [81], regarding the latter station.
In addition, we investigated patterns of changes in the normalized time-delays, t n , between the fast and slow split shear-waves, possibly associated with the observed evolution of seismicity in the WGoC. We point out three cases of linear decrease in t n values, which could be related to stress-relaxation due to microcrack coalescence. The first is associated with the period before the occurrence of the 2013 Helike swarm, between January and May 2013. The second was observed between November 2013 and January 2014, i.e., near the beginning of the fluid-driven, long-range seismic sequence. In particular, a cluster within group #4 occurred in late November 2013, while mid-to late December marks the beginning of the back-front, or the end of the "injection" period (~120 days). The pore-pressure pulse then travels at a speed of~0.13 km/day until January-February 2014 at the eastern end of group #4 and part of groups #6 and #8 (Y = 23-30 km in Figure 6). This epicentral area north of the Psathopyrgos fault has been previously proposed as the source of a strain transient event of equivalent M w = 5.3 that occurred on 3 December 2002 and was detected by a dilatometer at Trizonia Island [82]. To the east of this likely creeping fault patch, there is an area where larger events have occurred, including the recent M w = 5.3 event of 17 February 2021 (Figure 1; [83]). The last decreasing t n pattern was observed during the months before the 7 November 2014 event. A similar decrease pattern has been reported by [48] at stations ZIRI and SERG during that period, while at ROD3 station both an increase and a decrease period were identified; a pattern that could possibly be considered as precursory behavior.
Regarding the upper-crust anisotropy in the WGoC, a future goal would be the identification of the depth of the anisotropic layer. However, this would require recordings at the same site both on the surface and at various depths. Unfortunately, there are currently no such stations in the study region. Thus, we consider the whole medium between the hypocenter and the station as anisotropic. Another prospective step is the computation of the spatial averaging of the fast shear-wave polarization direction at the surface [84], in order to determine changes of ϕ at areas where no seismological stations are available, but this is beyond the scope of the present study.

Conclusions
To summarize, the presented results highlight the complex seismicity patterns that were observed during a two-year period in the Western Gulf of Corinth. The high resolution earthquake catalogue enabled the recognition of the activated structures at depth and their possible relation to mapped faults on the surface. The detailed spatiotemporal evolution of seismicity revealed interactions between neighboring structures via a migrating active patch, likely associated with increased pore-pressure. Measurements of shear-wave parameters were performed on a larger scale than ever before in the WGoC through an automated procedure, validating previous results from manual analysis and providing new information at different sites or study periods. The seismic monitoring of the WGoC over the past 20 years has been providing valuable information on the manners with which seismicity is evolved, driven by fluids that propagate through the microfracture network, enabling aseismic creep, triggering swarms and mainshock-aftershock sequences on their way. The high deformation rate in this part of Central Greece ensures that such patterns will continue to emerge and may eventually lead to significant and destructive events. The understanding of the dynamics of seismicity may therefore play an important role for the assessment of seismic hazard in the area.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/app11146573/s1, Section S1: Location/Relocation statistics for Catalogue 1, Section S2: Seismic moment measurements, Section S3: Automatic picking for Catalogue 2, Section S4: Additional tables and figures, File S1: 3D Interactive Matlab Figure, Dataset D1: Catalogue of seismicity, Dataset D2: Anisotropy measurements, Table S1: Absolute location error and other statistics per spatial group (see Figure 2) for the seismicity of 2013-2014 (Catalogue 1) in the Western Gulf of Corinth (WGoC). Table after [60], Table S2: Statistics of the station jack-knife test for the relocation robustness per spatial group of the 2013-2014 seismicity in the WGoC, with dx, dy and dz values representing perturbations from the true relocated hypocenter. The 95% confidence interval ellipse major semi-axis length refers to the horizontal ellipses of Figure S2, as estimated for each whole spatial group, Table S3: Mean values of SWS parameters and their respective errors per station for valid measurements (grade A-C) of Catalogue 3. Stations with additional data from Catalogue 2 are marked with a star (*) next to their station code. Errors δϕ correspond to the 95% c.i, whereas δt d and δt n are the typical error of mean. N is the number of valid measurements per station. The first line in bold corresponds to all valid measurements of Cat3, Figure S1: Location statistics for Catalogue 1. (a) RMS travel-time error, (b) horizontal (ERX, ERY) and vertical (ERZ) location uncertainties, taking into account an average arrival-time reading error of 0.13sec, (c) number of phases with non-zero true weight, (d) focal depth (including relocated hypocenters with HypoDD), (e) azimuthal gap and (f) epicentral distance of the nearest station with arrival-time data. Absolute location errors and RMS travel-time errors are nominal values reported by HypoInverse, Figure S2: Distribution of 95% c.i. ellipses for epicentral perturbations derived by the station jack-knife test for the robustness of the relocation procedure, Figure S3: Cross-sections a 1 -a 2 to c 1 -c 2 ( Figure S2) with vertical and horizontal bars representing the vertical and horizontal (square root of dx 2 + dy 2 ) average perturbation of each hypocenter, respectively, derived by the application of the station jack-knife test for the robustness of the relocation procedure, Figure S4. Example of seismic moment calculation by spectral fitting for the 15 July 2013, M w = 3.8 event on the short-period station TEME (top) before residual correction and (bottom) after residual correction. The blue and dashed red lines correspond to signal and noise spectra, respectively. The S-wave window is also presented in the inset panel at the bottomleft, Figure S5. Examples of spectral moment fit residual distribution (left) for the short-period station TEME and (right) for the broad-band station KALE, Figure S6. Same as Figure S4 Figure S9. Same as Figure S8, but for an automatic pick on an S-wave, Figure S10. Example of a rejected automatic pick due to excess noise resulting in scattered AIC minima, Figure S11. Example of automatic picks on P-and S-waves using the template matching method. The top waveform (black) belongs to the target event (to be picked). Red and green waveforms correspond to the best matching P-and S-wave templates, respectively. Arrival-time picks shown as vertical lines on the templates, with their colour corresponding to their correlation coefficient (CC). The final P-and S-wave pick is marked on the target waveform with a down-pointing arrow at the top of a red and green thick line, respectively. Available picks from more templates (with lower CC) are also shown on the target waveform, Figure S12. Same as Figure S11, but for a different target waveform, Figure S13. (Top) Map of relocated seismicity for 2013-2014 in the WGoC. Colours correspond to the 9 spatial groups in which the dataset was divided, along with the "miscellaneous" group #10 (gray). Beachballs depict the focal mechanisms of major events (Table 1). (Bottom panels) Vertical cross-sections a 1 -a 2 to p 1 -p 2 , drawn in a N10 • E direction, corresponding to the dashed rectangles shown on the map (top). Beachballs depict the far hemisphere projections of focal mechanisms for major events. Black dashed lines at the top of the cross-sections indicate the position of mapped active faults with inferred extension at a 60 • dip. For mapped faults' names refer to Figure 1, Figure S14. (Top) Map of relocated seismicity for 2013-2014 in the WGoC. Colours correspond to the 9 spatial groups in which the dataset was divided, along with the "miscellaneous" group #10 (gray). Beachballs depict the focal mechanisms of major events (Table 1) Figure S16: Rose diagrams of the distribution of the fast shear-wave polarization angle ϕ for valid measurements with grade A (best) to C (acceptable) from Catalogue 1 (only manually located/relocated events) per station (station code shown on top of each diagram). Colours correspond to the respective normalized time-delay t n . The red arrows indicate the average ϕ direction, Figure S17: Rose diagrams of the distribution of the polarization angle ϕ for null measurements (where ϕ is in the range of p ± 10 • or (p + 90 • )±10 • ) from Catalogue 1 (only manually located/relocated events) per station (station code shown on top of each diagram). Colours correspond to the respective normalized time-delay t n . The red arrows indicate the average ϕ direction of null measurements, Figure S18. Map of candidate events for SWS measurements from Cat3 (manual and automatic) in the WGoC during 2013-2014. Blue lines represent the average fast S-wave polarization direction, ϕ (orientation), and time-delay, t d (length), at the respective station (triangle, coloured according to network code). Coloured circles correspond to events with valid SWS measurements (grade A to C), whereas gray circles are events without valid measurements (only D, E or null), Figure S19. Gutenberg-richter law for the seismicity of Catalogue 1. M c = 1.8 is the magnitude of completeness. The b-value is calculated as 1.235±0.015 from linear regression (correlation coefficient R = −0.998) or b MLE = 1.185 from the maximum likelihood method, Figure S20. Composite Seismogenic Sources (orange) and Individual Seismogenic Sources (yellow) from the Greek Database of Seismogenic Sources (GreDaSS, [14]) in the WGoC area, mentioned in the main text. GRCS441: Trikorfo source, associated with Trikorfo and Sergoula faults, GRCS500: South Corinth Gulf source, associated with Sela, Lakka, w. Helike and e. Helike faults, GRCS524: Mamoussia source, associated with Pirgaki and Mamoussia faults, GRIS502: west Helike fault source, also associated with Lakka fault, GRIS503: Aigion fault source, also associated with Fasouleika and Selianitika faults, Figure S21. Schematic representation of Band-1 and Band-2 zones, modified after [85]. The receiver at the surface is marked with S. The ABCD plane of vertically aligned microcracks is oriented towards the mean fast shear-wave polarization direction, Φ. Band-1 contains seismic rays traveling towards the seismic station at directions contained within the solid angle EFGH-to-S (blue),at angles between 15 • and 45 • (shear-wave window) on either side of the crack-plane. Band-2 contains seismic rays that travel within the solid angle ADEHG-to-S (red) towards the seismic station, i.e., within ±15 • from the crack-plane, Figure S22. Maps of the temporal evolution of the 2013-2014 sequence in the WGoC (Catalogue 1, all spatial groups except #3 and #10). The colours represent relative origin time (see labels on the colour-scale of each panel). (a) January -August 2013: the map shows sporadic activity inside the gulf, with the most significant cluster being the 2013 Helike sequence. (b) September 2013 -14 July 2014: several clusters were observed at the western part of group #4, with activity going back and forth as shown with the red arrows. The dashed red ellipse shows a possible source of fluid injection. The black arrows show the direction (N100 • E) and speed of seismicity migration or inferred fluid diffusion (see Figure 5). Seismicity gradually migrated eastwards at~0.13km/day, activating the clusters with yellow colour (days 380-400). Fluid diffusion at a rate of 0.04km/day towards N100 • E precedes the triggering of the June 2014 sequence (dark blue), during which seismicity migration at 0.75 km/day was observed within the cluster (see also Figure A2). (c) 15 July -December 2014: bilateral seismicity migration observed in the spatial group #4, starting from the region marked with the red dashed ellipse (same as in panel b). Seismicity migration at a speed of 0.05km/day towards N100 • E between groups #8 and #6 is marked with a black arrow. The 7 November 2014 M w = 4.8 event (green star) triggers seismicity mostly southwards, Figure S23. Annotated snapshot of the interactive 3D figure (File S1), showing the spatial group #6 and associated structure planes, along with selected surface fault planes (gray) extended downwards to the depth of 8 km at a 60 • dip. The view is towards west. No vertical exaggeration has been applied, Figure S24. Annotated snapshot of the interactive 3D figure (File S1), showing the spatial groups #4 (green) and #5 (orange) and proposed associated structure planes. Selected surface fault planes (gray) are extended downwards to the depth of 8 km at a 60 • dip. The view is towards north-east. No vertical exaggeration has been applied, Figure S25. Annotated snapshot of the interactive 3D figure (File S1), showing the 2013 Helike swarm (spatial group #9) and associated structure planes that were activated (blue). The surface fault planes (gray) of Pirgaki and Kerinitis faults are shown in gray. The view is downwards and towards NNE. Although the swarm is mainly associated with the north-dipping Pirgaki fault, it also contained events with oblique-normal focal mechanisms, likely related to NW-dipping structures sub-parallel to Kerinitis fault [6].