The 2013 earthquake swarm in Helike, Greece: seismic activity at the root of old normal faults

The Corinth Rift in Central Greece has been studied extensively during the past decades, as it is one of the most seismically active regions in Europe. It is characterized by normal faulting and extension rates between 6 and 15 mm yr−1 in an approximately N10E° direction. On 2013 May 21, an earthquake swarm was initiated with a series of small events 4 km southeast of Aigion city. In the next days, the seismic activity became more intense, with outbursts of several stronger events of magnitude between 3.3 and 3.7. The seismicity migrated towards the east during June, followed by a sudden activation of the western part of the swarm on July 15th. More than 1500 events have been detected and manually analysed during the period between 2013 May 21 and August 31, using over 15 local stations in epicentral distances up to 30 km and a local velocity model determined by an error minimization method. Waveform similarity-based analysis was performed, revealing several distinct multiplets within the earthquake swarm. High-resolution relocation was applied using the double-difference algorithm HypoDD, incorporating both catalogue and cross-correlation differential traveltime data, which managed to separate the initial seismic cloud into several smaller, densely concentrated spatial clusters of strongly correlated events. Focal mechanism solutions for over 170 events were determined using P-wave first motion polarities, while regional waveform modelling was applied for the calculation of moment tensors for the 18 largest events of the sequence. Selected events belonging to common spatial groups were considered for the calculation of composite mechanisms to characterize different parts of the swarm. The solutions are mainly in agreement with the regional NNE–SSW extension, representing typical normal faulting on 30–50° north-dipping planes, while a few exhibit slip in an NNE–SSW direction, on a roughly subhorizontal plane. Moment magnitudes were calculated by spectral analysis of S waves, yielding b-values between 1.1 and 1.2 in their frequency–magnitude distribution. The seismic moment release history indicates swarm-like activity during the first phase, which could have acted as a preparatory stage for the second phase (after 12 July) that presented a more typical main-shock–aftershock behaviour. The spatiotemporal analysis reveals that the swarm has occurred in a volume that is likely related with the extension at depth of the NNE-dipping Pirgaki normal fault, outcropping ∼8 km to the south. The slow velocity of eastward migration of the epicentres during June implies triggering by fluids. The situation appears different in the second phase of the sequence, which was probably triggered by a build-up of stress during the first one. The relatively deep hypocentres of the 2013 swarm, compared to the shallower seismic layer within the rift, and their coincidence with the steeply dipping Pirgaki fault, favour an immature rift detachment model. Previous results from instrumental data indicate that approximately the same region had been activated during July–August 1991. The availability of the dense permanent seismological network data thus allowed for a detailed analysis of this crisis, a better understanding of its mechanical context and of the earlier events.


I N T RO D U C T I O N
The Corinth Rift has been identified as one of the fastest expanding (6-15 mm yr −1 of ∼N-S extension) and most seismically active continental rifts around the world (e.g. Jackson et al. 1982;Rigo et al. 1996;Clarke et al. 1997;Papadimitriou et al. 1999Papadimitriou et al. , 2010Briole et al. 2000;Hatzfeld et al. 2000;Lambotte et al. 2014). Detailed maps of the most prominent local tectonic features have been produced by focused studies, estimating their seismic activity through geomorphology and paleoseismology (Koukouvelas et al. 2001;Pantosti et al. 2004), while offshore structures have been mapped with high-resolution bathymetry by the Hellenic Center of Marine Research (Lykousis et al. 2007) and others (Brooks & Ferentinos 1984;Stefatos et al. 2002;Bell et al. 2009).
The presently most active normal faults of the rift are dipping north, resulting in a long-term subsidence of the northern coast and an upward displacement of the main footwalls, forming an asymmetrical rift structure . The study area is dominated by Alpine Carbonate formations of the H1 (External Hellenides Platform) and H2 (Cyclades-Pindus ocean crust) tectonic terranes and Plio-quaternary sediments in the northern shores of Peloponnese. The subsiding Hellenides limestone nappes to the north of the gulf are outcropping in the majority of sites, whereas to the south these formations are mostly covered by a conglomerate layer of several hundreds of meters thickness, which outcrops only on the footwall of the southern active faults Ghisetti & Vezzani 2004). Close to the seaside and offshore, on the hanging walls of the normal faults, the conglomerates are covered by finer Holocene deposits (sands and clay), that reach 150 m thickness in the Aigion harbor (Cornet et al. 2004;Pitilakis et al. 2004).
The study area has experienced several strong earthquakes during both historical and instrumental periods (Makropoulos & Burton 1981;Ambraseys & Jackson 1990;Papazachos & Papazachou 1997;Makropoulos et al. 2012;Stucchi et al. 2013), destroying cities such as Helike in 373 BC and causing severe damage to urban areas around the gulf (e.g. Eratine of Phokida, M w = 6. 3, 1965;Antikyra, M w = 6.2, 1970;Galaxidi, M w = 5.8, 1992;Aigion, M w = 6.2, 1995). Documents about strong earthquakes in the Corinth Rift exist for at least 2000 yr. Nevertheless, seismological catalogues can be considered complete for M w ≥ 6.0 only for the latest 300 yr (Console et al. 2013(Console et al. , 2015. The two most recent events of 1992 and 1995 (M w ≥ 5.8) in the western part of the rift activated shallow northdipping (30-35 • ) offshore faults Bernard et al. 1997) differing significantly from the 45 • to 50 • dipping planes of the 1981 major events (Jackson et al. 1982), at the eastern end of the rift, implying different rheological and/or structural conditions (Bernard et al. 2006).
The seismotectonic environment of the Western Corinth Rift is presented in Fig. 1. East of Psathopyrgos, the faults are oriented on average in an N100 • E direction. The coast of northern Pelopon-nesus at the southern part of the rift is characterized by generally enéchelon, north-dipping faults, with their outcrops separated by a few kilometres, while at the northern part of the rift the mapped south-dipping faults are smaller, more numerous, segmented, irregular and dense. In addition, the latter are antithetic to the main north-dipping faults of the southern side. The topography of the southern shoulder of the gulf and the shape of the coastline are also smoother compared to their northern counterpart, which is steeper and with many peninsulas . The dip angles of the faults at the surface range between 40 • and 70 • (Micarelli et al. 2006). An explanation for shallower dip angles ∼30 • observed in faults at the western part of the Corinth Rift is that they are not newly formed but constitute pre-existing weakness zones which were originally dipping at ∼45 • and have been tilted by ∼15 • on a horizontal axis because of an increase in the extension rate (Hatzfeld et al. 2000). These observations suggest that the Western Corinth Rift is an asymmetric half-graben with the rifting process taking place at its southern shoulder. To the south (outside the bounds of Fig. 1), there are traces of some older north-dipping faults, striking ∼N280 • E, parallel to the newer ones. There is seismological Rigo et al. 1996;Bernard et al. 1997;Lambotte et al. 2014) and geophysical evidence (Le Meur et al. 1997;Pham et al. 2000;Gautier et al. 2006) that the most recently developed faults, bordering the southern shore, are rooting at a depth of 6-8 km beneath the gulf, in a zone which generates continuous background microseismic activity while producing the observed deformation (Bernard et al. 2006). This is supported by several low-angle, north-dipping, normal focal mechanisms which have been determined for events located at depths between 9.5 and 10 km (Rietbrock et al. 1996;Rigo et al. 1996). Hatzfeld et al. (2000) suggested that the seismicity at these depths is related to small faults situated near the transition between brittle and ductile crust and not to a large low-angle north-dipping normal fault. However, recent studies, based on detailed earthquake relocation and fault plane solutions determination, have shown that the seismicity beneath the gulf itself occurs within an ∼1 km thick layer and that a growing shallow dipping detachment exists under the northern coast (Godano et al. 2014;Lambotte et al. 2014). At the southern shore, the Pirgaki-Mamousia fault has been inactive for over 0.4 Ma, probably for 0.7-0.8 Ma (Hemelsdaël & Ford 2014), according to the dating of the sedimentation at the Vouraikos delta (Ford et al. 2007) and the Kerinitis delta (Backert et al. 2010). The Helike fault, on the other hand, which outcrops further north, has been active for 0.7-0.8 Ma (or less, depending on whether there has been acceleration in the uplift rate), according to modelling of footwall uplift using marine terraces (McNeill & Collier 2004;Hemelsdaël & Ford 2014), and is considered as a mature fault with a significant seismic potential.
The background seismicity of the gulf and the episodic swarms or clusters of activity are being recorded by the local Corinth Rift Laboratory (CRL) network since 2000 .  , black open circles represent selected seismicity of the period July-August 1991 , green circles represent the relocated epicentres of the 2001 swarm near Agios Ioannis Kapetanidis et al. 2010) and blue circles represent the relocated epicentres of the 2013 swarm (this study). Local stations are displayed with yellow triangles. The fault lines are after N. Meyer (PhD thesis, in preparation), mostly based on Moretti et al. (2003), Palyvos et al. (2005) and Bell et al. (2008).
In this study we investigate an earthquake swarm that occurred in the period between 2013 May 21 and August, with its epicentres covering a small region roughly located north of the eastern end of the West Helike fault scarp (Fig. 1, blue circles). The same region was activated during the period July-August 1991 (Fig. 1, black circles), partly recorded by a seismological experiment carried out with a network of 51 digital stations . The results were derived using a custom 1-D velocity model and suggested that an earthquake sequence occurred on the Pirgaki fault (referred in that study as Mamousia fault), at depths of ∼7-8 km. Another swarm occurred between May and July 2001 (Fig. 1, green circles) with its epicentres between the outcrops of Pirgaki and Helike faults, near station Agios Ioannis (AIOA) of CRLNET. However, both the calculated focal mechanisms and relocated hypocentres indicated that it activated an NW dipping fault plane and this swarm was associated to the SW-NE trending, NW-dipping, transverse Kerinitis fault, whose extension to the surface is traced in the valley of river Kerinitis (Pacchiani & Lyon-Caen 2010). It is within the scope of this study to gain a better understanding of the seismotectonic regime in the area of the Western Corinth Rift which lies near the junction of the Pirgaki and Kerinitis faults, as well as of the large seismically active layer beneath the gulf.
Earthquake swarms are bursts of seismic activity which lack a distinct, strong, initiating event that would be characterized as a main shock in a typical sequence and would also explain its spatial extent (Mogi 1963;Zobin & Ivanova 1994). Instead, they may contain several major events of comparable magnitude while the strongest one often occurs in the middle or towards the end, rather than the beginning of the sequence (Ferrucci & Patanè 1993;Vidale & Shearer 2006;Lohman & McGuire 2007;Legrand et al. 2011), although this may not always be the case (Benetatos et al. 2004;Lyon-Caen et al. 2004). Swarms are most common in volcanic settings (Hill 1977;Bianco et al. 2004;Smith et al. 2004;Massin et al. 2013;Sigmundsson et al. 2014) but have also been observed in non-volcanic regions (Jenatton et al. 2007;Lohman & McGuire 2007;Ibs-von Seht et al. 2008;Roland & McGuire 2009;Pacchiani & Lyon-Caen 2010;Telesca 2010). Mogi (1963) has suggested that swarms may occur in highly fractured zones or otherwise heterogeneous media. These conditions, combined with a gradual increase of regional stress, enable the local concentration of high stress in a multitude of cracks or small faults. In the case where this volume is seismically activated, the lack of a distinct major active fault limits the maximum magnitude potential, as well as the spatial extent of a generated swarm. Such conditions can be present, for example, in volcanic environments (Hill 1977) or at the intersection of buried tectonic features (Bisrat et al. 2012). Swarms can also be caused by a source of highly localized stress. This mechanism usually implicates pressurized fluids, either of magmatic nature (Hill 1977;Ferrucci & Patanè 1993;Sigmundsson et al. 2014), hydrothermal activity (Bianco et al. 2004;Massin et al. 2013), meteoric origin (Kraft et al. 2006;Jenatton et al. 2007;Leclère et al. 2012;Kassaras et al. 2014b) or mantelic CO 2 degassing (Fischer et al. 2014), which alter the stress field and trigger seismicity (Roland & McGuire 2009;Smith et al. 2004;Hainzl & Ogata 2005). The fluids are diffused through the microcrack network, increasing pore fluid pressure while reducing the effective normal stress, causing seismicity to expand radially away from the initial source of injection (Shapiro et al. 1997;Fournier 1999;Hainzl et al. 2012), usually at low rates of 10-100 m d −1 (Pacchiani & Lyon-Caen 2010;Kassaras et al. 2014b). This phenomenon has been well documented in cases of seismicity induced by hydraulic fracturing (Fehler et al. 1998;Parotidis et al. 2005;Albaric et al. 2014). High migration rates (∼1 km h −1 ) have been observed in swarms caused by dyke propagation (Sigmundsson et al. 2014) or driven by aseismic creep (Lohman & McGuire 2007). Static stress changes caused by a swarm itself may also play a crucial role in its evolution, especially when the activity leaps abruptly to different focal areas (Aoyama et al. 2002). The fluid-driven mechanism of swarms and the possible presence of asperities provide favourable conditions for the generation of repeating earthquakes or multiplets (Hemmann et al. 2003;Bourouis & Bernard 2007;Bisrat et al. 2012;Massin et al. 2013), whose characteristic is the high similarity in their waveforms (Geller & Mueller 1980;Poupinet et al. 1984). The frequency-magnitude distribution of swarms in volcanic areas is often characterized by unusually high b-values, reaching up to 2.5 (Smith et al. 2004;Legrand et al. 2011), while in other regions they may be either slightly above or below unity (Jenatton et al. 2007;Ibs-von Seht et al. 2008;Pacchiani & Lyon-Caen 2010;Telesca 2010).
In the framework of the present study, the earthquake swarm of 2013 May 21 to August 31 in Helike was extensively studied. A large number of events, recorded by the local stations in the Western Corinth Rift, were detected and manually analysed. Their absolute locations were improved using a local velocity model that was determined by an error minimization method. The seismic waveforms were cross-correlated and highly similar events were grouped in multiplets. Double-difference relocation was applied, using both catalogue and cross-correlation differential traveltimes, to reduce relative location uncertainties and distinguish subclusters within the main body of the swarm. Moment magnitudes were calculated by spectral analysis of S waves, in order to study their frequency-magnitude relation and the temporal evolution of their cumulative values. The seismo-tectonic study was completed by the determination of focal mechanism solutions for over 170 events using regional waveform modelling, for the major events, and a combination of P-wave first motion polarities, S/P amplitude ratios and S-wave polarizations for the smaller events. Spatial and temporal characteristics of the sequence were studied in detail to investigate its dynamics. The results are discussed in the context of the seismotectonic regime of the broader area of Western Corinth Rift, with emphasis on the neighbouring major faults.

DATA A N D L O C AT I O N P RO C E D U R E
The area around the West Helike fault (Fig. 2) is located within the bounds of the CRL local seismological network. Since its in-stallation in mid-2000, the CRLNET has recorded many tens of thousands of earthquakes in the Western Corinth Rift, with magnitudes as low as M ≈ 1 Lambotte et al. 2014). In collaboration with the Seismological Laboratory of the University of Athens and the Seismological Laboratory of Patras University, real-time telemetry has also been established in the recent years, enabling the automatic detection and location of current seismicity (http://www.geophysics.geol.uoa.gr, http://crlab.eu, last accessed 25 June 2015).
In this study, data from more than 15 seismological stations at epicentral distances up to 30 km have been used, including stations of the CRLNET, the Hellenic Unified Seismological Network (HUSN) and a temporary station at Helike (HELI), installed on 2013 May 23. More than 1500 events belonging to the earthquake swarm have been detected during the period between 2013 May 21 and August 31. Their P-and S-wave arrival times were manually picked to obtain well constrained initial locations using the algorithm HYPOINVERSE (Klein 1989), with a local velocity model for the area of Western Corinth Rift .
An average traveltime residuals and location uncertainties minimization technique (Kissling et al. 1994;Kaviris et al. 2007;Papadimitriou et al. 2010) was applied for the determination of a custom velocity model (CM) for this particular focal region and network geometry. This was done in order to acquire better estimates of the initial hypocentral locations and to take advantage of the available local network. The procedure was performed on a subset of 300 events located with uncertainty of less than 1 km and having RMS <0.15 s, after estimating the optimal distance weighting parameters, limiting station data up to epicentral distances of ∼30 km. The resulting optimal velocity model (CM) is presented in Table 1 and Fig. 3(a), along with the initial model (RM) by Rigo et al. (1996). The latter was derived from seismicity that was recorded during the period July-August 1991 by a dense temporary seismological network which expanded on both coasts of the Western Corinth Rift and the Rio graben, with the area of Helike being at its SE end. Seismic activity in approximately the same area, between the W. Helike and Aigion fault scarps, was present during the 1991 experiment; however the data also included hundreds of earthquakes with epicentres offshore and in several of the surrounding onshore regions.
In the model determined in the present study (CM), the V p /V s ratio was found equal to 1.80, as also in the RM, using a modification of the Wadati (1933) method, as formulated by Chatelain (1978), which estimates the V p /V s ratio by the slope of differential S-wave arrival times (T Si -T Sj ) plotted against the corresponding P-wave arrival-time differences (T Pi -T Pj ) for all combinations (i, j) of station pairs (Fig. 3b). The same value (V p /V s = 1.80) also minimizes location uncertainties. However, the P-wave velocities have been found larger by 0.2-0.6 km s −1 , compared to the initial model (Fig. 3a). Higher velocities compared to the RM have also been observed in another seismicity study in the Aigion region (Novotny et al. 2007), more specifically in the upper layers, including a higher ratio V p /V s = 1.83. It is noteworthy that the estimated ceiling depths are very similar to the RM. The homogeneous half-space starts at a shallower depth, as for this particular dataset and with the stations situated in epicentral distances up to ∼30 km most seismic rays are upgoing.
The results obtained using both velocity models are presented in Fig. 4. The epicentres have become less diffuse and the depths are better constrained using the velocity model determined by the present study. This is particularly evident in an apparently southdipping lobe at the northern part of the spatial distribution of the   (Table 4), FP2 fault plane related to events with oblique-normal, NW-dipping focal mechanisms. Table 1. Comparison between the local velocity model by Rigo et al. (1996) and the custom velocity model calculated for the 2013 swarm in Helike (this study). Rigo et al. (1996) (Fig. 4b), which has been embedded in the main body of the swarm with the CM (Fig. 4d). The RMS errors have been improved by about 17 per cent (median value 0.10 s), while the horizontal and vertical location errors were reduced by 6-7 per cent (median values 0.48 and 0.76 km, respectively). The depth range is similar in both models, roughly between 8 and 10 km. The custom velocity model improved the results by minimizing the traveltime residuals for the specific ray paths travelling from the foci of the swarm's events to all the stations of the local network. However, it is only an average approximation and, being 1-D, it does not contain information concerning lateral velocity variations, geology or topography characteristics at each station. In order to further improve the initial locations we have calculated and applied station delays in both P and S waves, using the mean traveltime residuals of well-located events. These static station corrections account for unmodelled lateral variations in the velocity structure and station elevation which may cause a bias in the traveltime residuals (Richards-Dinger & Shearer 2000). The station corrections improve the horizontal concentration and enhance the dip towards the north (Fig. S1). Numerically, there is only slight reduction in the location and RMS errors. Initially, an average timing error ERR = 0.3 s was considered, which is more than double compared to its true average value, as will be estimated in a later section. This includes errors in the calculated traveltimes caused by unmodelled velocity  Rigo et al. (1996) and the one determined in the present study (CM), (b) V p /V s ratio estimated using a modification of the Wadati method (Chatelain 1978). structure, mostly reflected by the average RMS error, which is also the suggested value of the ERR parameter for HYPOINVERSE. Source-specific station corrections (Richards-Dinger & Shearer 2000) were not deemed necessary, as the events are confined in a small region compared to (most) epicentral distances of the stations, while any remaining uncertainties were minimized by the double-difference relocation technique that follows.

M U LT I P L E T C L A S S I F I C AT I O N -R E L O C AT I O N P RO C E D U R E
The absolute locations in the initial catalogue have been improved by the use of a custom velocity model and the application of station corrections. There are still, however, significant uncertainties in the relative locations of the hypocentres, which make it difficult to distinguish details within the swarm, such as the existence of subclusters. In order to further constrain the hypocentral locations, we have applied a double-difference relocation procedure using the HypoDD algorithm (Waldhauser & Ellsworth 2000), incorporating both catalogue and cross-correlation differential traveltime data. This method works on the assumption that the ray paths of two earthquakes are practically similar when the distance between their foci is much smaller compared to their (common) hypocentral distance from a station. Hence, the differences in their traveltimes can be attributed to the distance between their foci, measured along the initial ray direction. The relocation works by minimizing the double-difference between observed and calculated traveltimes of correlated pairs of events, both for P-and S-waves. As a result, relative location errors caused by unmodelled velocity structure are minimized (Waldhauser & Ellsworth 2000), as long as the chosen velocity model is appropriate (Michelini & Lomax 2004). Direct differential traveltime measurements can also be used by incorporating waveform cross-correlation data to account for arrivaltime picking inconsistencies within multiplets (Kapetanidis & Papadimitriou 2011).
It is expected from the theory of elastic wave propagation that earthquake events with similar focal parameters produce similar waveform recordings (Geller & Mueller 1980;Poupinet et al. 1984). Such events are characterized as multiplets or repeating earthquakes and are of particular interest in studies of hypocentral relocation (Hauksson & Shearer 2005;Bohnhoff et al. 2006;Waldhauser & Schaff 2008;Statz-Boyer et al. 2009), propagation velocity variations (Poupinet et al. 1984;Nishimura et al. 2000;Pandolfi et al. 2006;Cociani et al. 2010), empirical Green's functions (Ichinose et al. 1997) and others. In this study, we measure the similarity between pairs of waveform recordings of equal sampling rate (timeseries x(n) and y(n)) using the normalized cross-correlation function, XC(d) (eq. 1): where N is the number of samples of the longer time-series (the shorter one is padded with trailing zero samples to reach N) and d is the time-lag measured in samples, with 1 − N ≤ d ≤ N − 1. Both waveforms are bandpass filtered prior to the cross-correlation procedure to remove long-period and high-frequency noise. For local microearthquakes, the frequency band between 2.5 and 23 Hz is usually adequate in most cases (Fig. 5), as described in more detail in Appendix A. The global maximum XC max = XC(t m ) is considered as a measure of the similarity between the two waveforms. The cross-correlation data can provide direct measurements of arrivaltime picking inconsistencies (Seggern 2009). Their incorporation into the relocation procedure can lead to reduction of location uncertainties due to arrival-time reading errors (Waldhauser & Ellsworth 2000). We applied full signal cross-correlation on events recorded on the three components of stations TEME, situated within the epicentral region, and LAKA, in an average distance of ∼10 km from the centre of the swarm. A combined cross-correlation matrix was constructed by averaging the XC max from the six measurements. This was achieved using a vector-length approach (Kapetanidis √ N c , where N c ≤ 6 (depending on the availability of the common component recordings) is the number of XC max measurements for each event pair. The hierarchical clustering was calculated using the nearest-neighbour linkage algorithm on the matrix data.
Events were grouped in multiplets according to their waveform similarity, using a correlation threshold C th = 0.80, an optimal value which maximized the difference between the size of the largest multiplet and the sum of correlated events (Kapetanidis et al. 2010). This resulted in 113 multiplets containing a total of 1189 events (∼80 per cent of the catalogue). Afterwards, cross-correlation measurements (XC max and t m ) were made for all combinations of pairs of events within each multiplet for P or S waves of all stations within a 30 km radius from the epicentres in all working components. The t m values were calculated with subsample precision (1 ms) by fitting a resampled parabola on the global maximum of the cross-correlation function (Deichmann & Garcia-Fernandez 1992;Schaff et al. 2004). At this stage, a lower correlation threshold (0.5) was used to allow for indirect relatives (Kapetanidis & Papadimitriou 2011), but measurements with XC max < 0.5 or |t m | > 0.5 s were rejected.
The distribution of time-lags, t m , for the P or S waves of a station can provide an estimate of the corresponding arrival-time reading error (Seggern 2009) for a certain multiplet, supposing that, on average, the picks are scattered around the correct arrival time (for the observed similar waveform onset). Reading errors are expected to be smaller for closer stations than for farther ones, and larger for S than for P waves, due to the superposition of the S-wave onset with the P-wave coda or the waveforms of secondary phases. The mean absolute time-lags (averaged over all events at all stations) were found to be about 0.05 and 0.09 s, for P and S waves, respectively. This justifies the choice of the ERR parameter for HYPOINVERSE equal to 0.15 s, including other timing errors or traveltime residuals (given that RMS ≈ 0.07 s for the final model with station corrections). The reduced error parameter yields median location uncertainties ERH = 0.25 km and ERZ = 0.40 km.
The relocation procedure was divided in two main stages, the first with strong a priori weights on catalogue data and the latter on cross-correlation data, resulting in a stable solution with minimal origin shift. A second pass was also considered, using the relocated foci of the first one as initial positions and repeating the procedure on each of the individual multiplets separately, using the Singular Value Decomposition method for the smaller ones and the weighted least-squares (LSQR) method for the larger ones.
Consequently, the dispersion of the spatial distribution was significantly reduced, especially in the vertical axis, and multiplets were compacted in dense spatial clusters. A total of 1489 events (99 per cent) were successfully relocated (Figs 6 and 7), including those which were only linked through catalogue data (they did not belong to multiplets). Both horizontal and vertical scattering were decreased, with the epicentres becoming concentrated in a 3 km × 2 km patch, oriented approximately E-W. A 3-D visualization of the relocated spatial distribution suggests that the events have occurred at a depth range between 8 and 10 km, on a curved surface with an uncertainty of about ±250 m (the width of the distribution can be observed in the perpendicular cross-sections of Figs 7b and d), striking approximately N270 • -280 • E at its north-dipping middle and eastern part and N260 • E at its NNW-dipping western side. Its estimated horizontal extent is ∼3.7 km in an N100 • E direction. The surface is not evenly covered with hypocentres but the activity is focused in strong spatial clusters, separated by gaps, some of which may correspond to asperities which have hosted the stronger events of the sequence.  A spatial segmentation of the distribution in groups was performed to provide visual aid for better understanding of the hypocentral shifts achieved by the relocation procedure, but also to enable the spatially targeted analysis and descriptions in the next sections of this study. The inter-event hypocentral (3-D) distance matrix of the relocated foci was calculated and Ward's linkage (Ward 1963) was applied to create a clustering hierarchy. A threshold value was imposed, dividing the distribution in exactly 11 clusters (Figs 6b and 7b), as this number was considered optimal for the adequate separation of spatial clusters which could also be discriminated by visual inspection. Although the clustering is 3-D, the horizontal projections of the clusters are also well-separated. The clusters appear Downloaded from https://academic.oup.com/gji/article/202/3/2044/610654 by CNRS -ISTO user on 17 August 2020 to be divided in two rows (both horizontally and vertically), with the exception of the middle region of the swarm which is shared by clusters #3, 8 and 9. Their dynamics will be examined in detail in the following section.

S PAT I O T E M P O R A L / M U LT I P L E T A N A LY S I S
The characteristics of the spatiotemporal distribution, which can shed light on the dynamics, reveal patterns in its behaviour and provide clues for the interpretation of the phenomenon, are being thoroughly examined in this section. The Helike seismic swarm was initiated on 2013 May 21 by a series of small earthquakes, escalating up to an M w = 3.1 event on 23 May 05:31 UTC, spreading through the polygon area marked in Fig. 8(a). A second outbreak (Fig. 8b) began in 2013 May 27, culminating to an M w = 3.6 earthquake on 28 May 01:13 UTC. The latter event occurred at the westernmost point of the main cluster of activity of period a and triggered a patch of the fault surface further to the west ( Fig. 9, 1st star). This was followed by another burst of seismicity (period c) including a major M w = 3.7 event on 31 May 08:57 UTC (2nd star in Fig. 9, cluster #3), which occurred in a spatial gap of the previous activity, approximately in the middle of the swarm's epicentral distribution. By the end of period c, since 10 June, the seismicity started to migrate gradually towards the eastern portion, while at the western and middle part the activity diminished. The eastwards migration continued during period d with small bursts every 2 d, including one with three earthquakes of magnitude between 3.0 and 3.5, followed by two more clusters, the latter of which involved an M w = 3.6 event on 27 June 17:24 (the fourth event of the sequence with M w ≥ 3.5) that generated smaller events at the easternmost edge of the distribution.
A sudden change in the spatiotemporal pattern occurred in mid-July ( Fig. 8(e), with a dashed circle marking the epicentral area in all panels for spatial reference), when the western portion of the swarm was activated, with several M w ≥ 3.0 events, including an M w = 3.7 earthquake on 15 July, 20:07:55 UTC, one of the major events of the sequence. It is worth noting that the latter occurred in a gap between spatial clusters #1, 2 and 6 which were initiated 1 d earlier with some smaller events. The last period (f) is characterized by a small spatiotemporal cluster at the middle-northern part of the swarm and a gradual decline of activity towards the end of August.
Several differences exist between the first (periods a-d) and the second (periods e-f) phase of the sequence. The 'middle region' (in the range of Y between 0 and 0.5 km in Fig. 9) is mostly active during the first 30 d. The eastern region (Y > 0.5 km), active during the first phase, exhibits a gradual migration of seismicity from the middle of the swarm towards the east, with a ∼400-500 m long patch of activity propagating at a rate of ∼30-40 m d −1 , mostly during June. In the second phase, the activity is almost exclusively occupying the westernmost portion of the swarm, for a range of Y between −1.5 and −0.5 km (Fig. 9), with the exception of a small spatiotemporal cluster in the middle of the swarm during period f. The focal depths during phase #1 are constrained between 8.5 and 9.5 km, gradually reaching ∼10 km by the end of the second week (Fig. 9, bottom panel). In phase #2, the activity is more constrained between 9 and 10 km, mostly involving spatial clusters #6 and 1 with persistent activity in a constant volume. Cluster #2 was also triggered at the beginning of period e and cluster #3 by period f, both concentrated at about 9.5 km.
The differences between the two phases were also evident in both the cross-correlation matrix and the consequent hierarchical clustering. A history of the occurrence of repeating events, with increasing IDs reflecting the origin time of the first event in each multiplet, can be observed in Fig. 10. Sudden bursts of seismic activity are usually accompanied by the generation of new multiplets in regions previously unbroken. With the initiation of phase #2, about 20 new multiplets were formed during the first 4 d, some of which containing over 10 events. There are also a few groups of repeating earthquakes, mostly doublets with the first event of the pair having already occurred in phase #1 (among IDs 4-25 and 35-55). The same is true for one of the larger ones (ID #44), with all events but one in phase #2. This shows that there is some overlap between the activated regions in the two phases. The largest multiplet, by far, is the first that was generated (ID #1), containing 158 events, and occurring between 2013 21 May and 9 June. This corresponds to the spatial cluster #9, which lies in the middle of the swarm. The largest 18 multiplets contain over 20 events and have been generated in regions of strongly spatially (and often temporally) clustered activity.
The hypocentral relocation uncertainties were assessed using bootstrap techniques described in detail in Appendix B. Three different sources of errors were examined: (1) the influence of stations, (2) the influence of initial locations and (3) the influence of each event to the relocation of the rest. The results, which are presented in Table 2, reflect the well-defined spatial clusters which have been formed, with uncertainties of the order of 100 m or even lower for events with strong cross-correlation links.

F O C A L M E C H A N I S M S
The highly detailed spatial distribution provides some clues on the geometry of the activated fault surface. In order, however, to gain a more complete understanding of the seismotectonic characteristics of the 2013 seismic crisis in Helike, Focal Mechanisms (FMs) were initially estimated by P-wave first-motion polarities (FMP) for a wide range of events with an average of 17 polarity measurements, using azimuths and take-off angles as calculated for the relocated foci. A grid-search was performed to identify all combinations of strike, dip and rake which give solutions that are in agreement with the observed FMP. While a significant number of FMs could be safely derived by this method alone, S-Wave Polarizations (SWPs) and S to P amplitude Ratios (SPRs) were also measured and incorporated, in a second stage, to better constrain the solutions. The differences between the observed and theoretical directions of SWP, as well as between the logarithm of observed and the logarithm of theoretical SPR were included in a weighting scheme. This was applied in the calculation of an average solution for each FM, thus acquiring a strike, dip and rake combination which confirms all FMP and is also in favour of minimizing the differences between theoretical and observed SWP and SPR (for more details on this procedure refer to Appendix C).
The seismic moment tensors of major events were calculated by synthetic waveform modelling at regional distances using the frequency-wavenumber integration method (Bouchon 2003). Synthetic seismograms were generated using the discrete wavenumber method and a modification of the Axitra code (Bouchon 1981) for five elementary types of faulting. The obtained synthetic seismograms were compared to the observed ones, to which instrumental response correction and bandpass filtering in the range 0.03-0.125 Hz was applied. Then, a grid-search procedure was performed June-11 July, (e) 12 July-31 July and (f) 1-31 August. The colours correspond to different spatial clusters. The dashed black polygon and dashed red circle correspond roughly to the epicentral area activated during periods a and e, respectively, and are used for spatial reference.
to determine the values of strike, dip, rake and focal depth that correspond to the minimum normalized RMS misfit between observed and synthetic waveforms (Fig. 11). This methodology has been used previously in Western Greece (Papadimitriou et al. 2012), providing stable solutions even for events of moderate magnitudes (M w < 4.0). Waveform modelling with ISOLA code (Sokos & Zahradník 2008) was also used in this study to calculate the moment tensor solutions of major events, providing similar results.
The FMs which were derived by waveform modelling were compared to the ones derived independently by FMP, weighted by SWP and SPR, and were found to be highly compatible, with the one method confirming the results of the other. The only difference in the modelling approach concerns the focal depth for some of the events (being larger than 10 km instead of between 8 and 10 km). In these cases, the depths of the relocated foci were considered and the inversion was repeated for the strike, dip and rake parameters only. The fixed depth restriction at 9 km was found to not affect the solutions in a significant degree, although the waveform misfit was slightly larger. These differences are mainly due to the different (regional) velocity model that was used for the generation of the synthetic waveforms and to the incorporation of data from stations in epicentral distances larger than 30 km. The complete list of results for FMs derived by waveform modelling is presented in Table 3. The normalized RMS misfit, which is used as a measure of the solutions' uncertainties, is computed using the formula of eq. (2):  Fig. 6a). Colours correspond to different spatial clusters (similar to Figs 6b and 7f). The circle size is proportional to the event magnitude. A set of 11 major events M w > 3.2 is displayed with stars (see events marked with asterisk in Table 3 where 'obs' represents the observed waveforms, 'syn' the synthetic waveforms and 'VR' the Variance Reduction, which can be alternatively used as a measure of the goodness of fit. A composite solution per spatial cluster was calculated for the events with FMs determined either by FMP or waveform modelling ( Fig. 12 and Table 4), by averaging their normalized moment tensors, weighted by the number of FMP available for each individual solution. A similar procedure has also been applied in Trichonis Lake, Western Greece (Kassaras et al. 2014a), and in the Santorini Volcanic Complex . The analysis re-vealed that most clusters exhibit a consistent composite solution that is more or less similar to the mean mechanism, averaged over all 179 individual observations. This corresponds to a normal fault, striking ∼258 • and dipping 39 • to the north, which is favoured against its auxiliary, south-dipping plane, as it is in agreement with the plane indicated by the relocated spatial distribution (Fig. 13).
However, some events have been found with significantly different P-waveform onsets. While the FMP in most observations are dilatational in the middle and compressive at the top and bottom of the stereo-net, group 'X', displayed on the right hand side of   Fig. 12, has compressive polarities in the SSW half of the stereo-net while dilatations are limited to the NNE half. These solutions were mostly found in spatial clusters #2 and #6, at average depths of ∼9.5 and ∼9.2 km, respectively, located almost exclusively at the western part of the swarm, with occurrence times mostly on 28-29 May or after 14 July. However, they were grouped separately for the calculation of a composite mechanism ( Fig. 12 and Table 4, group 'X'). Its rake is marginally positive (considering the angle wrapping at 180 • ), primarily due to the influence of the stations TEME, ALIK (compressive first motion) and LAKA (dilatational motion) which constrain the subvertical auxiliary nodal plane. The subhorizontal rupture plane dips slightly WNW, however the slip direction of the footwall is towards SSW. This is in agreement with the motion of the low-angle detachment which lies underneath the northern coast of the Western Corinth Rift at similar depth, where evidence of similar focal mechanisms has been found previously (Rietbrock et al. 1996). The co-existence of purely normal faulting events in the same volume with events that have occurred on a subhorizontal rupture plane with horizontal slip and a small reverse component increases the complexity in the western part of the swarm. It is noteworthy that the major event of 14 July also belongs to this group and similar FMs have been derived through waveform modelling. FMs belonging to the spatial clusters with CLIDs 2, 5, 6, 7, 9 and 11 are mostly suggestive of a north-dipping normal fault. FMs for CLID 8 have an NW-dipping nodal plane and a steeper SSW-dipping counterpart. In this case, the latter could well be a small antithetic rupture plane (Fig. 13b, fault line labelled 'F8' in cross-section f 1 -f 2 ). The same could be true for CLID 10, or at least for some of its individual events. Both clusters #8 and #10 are located in relatively shallower average depths of ∼8.8 and 8.5 km, respectively. CLID #1, located at the westernmost and deepest (∼9.7 km) part of the swarm, dips roughly N-NNW, in agreement with the change in the curved surface that can be fitted to the spatial distribution. Clusters #4 and #3, the latter of which includes the oblique-normal, major event of 31 May, although they are located in the middle of the swarm, also suggest an NW-dipping or a steeper SE-dipping rupture plane. In general, the dip angles of the majority of the composite solutions vary near 45 • , while the corresponding rake is around −90 • , indicating mostly normal faulting with a small lateral component in some groups. Both pure-and oblique-normal solutions are obtained throughout the whole sequence, both spatially and temporally. Table 3. Focal mechanism solutions for 18 major events of the 2013 sequence in Helike, derived by synthetic regional waveforms modelling, where CLID is the ID number of the spatial cluster where each event belongs. The hypocentral co-ordinates are from the results of the relocation procedure. The moment magnitude has been calculated using the spectral fitting method described in this study. Asterisk on the left marks the major 11 events presented with stars in several figures. For most groups, the P-axes are generally subvertical, oriented in various directions but with their plunge angle around or above 80 • . The T-axes, on the other hand, are mostly horizontal, with plunge angles below 10 • , in a general NNE-SSW direction. That is with the exception of clusters #8 and #3, where the T-axis is less horizontal and, more importantly, group 'X', which is the obvious exception, as it represents rupture on a subhorizontal plane. Even in the latter case, both P-and T-axes for the composite solution of group 'X' are directed roughly NNE-SSW, with trends 29 • and 194 • , respectively. These observations are compatible with both the known tectonic regime of the Western Corinth Rift and the regional extension field.

S E I S M I C M O M E N T
In the framework of the present study, seismic moment magnitudes (M w ) have been computed using a spectral fitting method which estimates the source parameters in the frequency domain from the inversion of earthquake S-wave displacement spectra. We assume a Brune-type (Brune 1970) source, with ω −2 spectral decay and invert for the seismic moment M o , the corner frequency f c and attenuation t * , simultaneously. The procedure introduces signal/noise ratio as a weighting factor to reduce the influence of noise that can mask the f c and the spectral level of M o at low frequencies (Fig. 14d).
Data processing and quality control play a crucial role in the determination of M w : it is particularly important for small earthquakes, as seismic noise contaminates the low-frequency spectral amplitudes. For this reason, we selected records characterized by a signal-to-noise ratio larger than 2 and for which an accurate manual picking of the P-wave arrival was available. S waves are analysed by combining the two horizontal components and selecting a 5 s time window starting at 1 s before the S-wave onset (Fig. 14a). The use of a 5 s time window around the direct S-wave arrivals is a compromise between taking a window of signal large enough to include the most important part of the energy but not too large in order to reduce the influence of propagation effects that can affect the estimation of the parameters. When an S-wave arrival has not been manually picked, a V p /V s ratio of 1.8 is used to estimate its arrival time from the P-wave picks. A pre-P-wave noise window is selected with duration as long as possible to enable the correct recovery of any long-period noise. Optimally, a 10 s noise window that ends 10 s before the P-wave arrival time is chosen.
In order to reduce distortions due to the windowing of the signals, a cosine taper function is applied to S-wave time-series before computing the amplitude spectrum. A bandpass filter between 1 and 50 Hz is used for the short-period recordings and between 0.5 and 50 Hz for the broad-band recordings. Finally, the S-wave displacement spectra from velocity/acceleration time-series are computed from integration/double-integration of velocity/acceleration in the frequency domain.
We applied a nonlinear least squares algorithm to fit the observed S-wave displacement spectra in the frequency range of 0.5-40 Hz (Fig. 14c). In this approach, the source models are determined using progressive correction of the spectra for the effects of site amplification and de-amplification. The spectral analysis was, in fact, performed in three steps: (1) Preliminary spectra inversion to obtain M o , f c and t * .
(2) Computation of non-parametric Site Response Function (SRF) to take into account characteristic amplification or deamplification for certain frequencies at each station (Fig. 14b). SRF is computed by averaging over all the events the difference between the observed spectra and the theoretical one.
(3) Computation of M o for each record from the displacement spectra corrected for the SRF. We invert separately the spectra for each available record passing the noise criteria. For each event, the source parameters are computed by averaging the M o that has been determined in all available stations and the standard deviation is assigned as the uncertainty. Fig. 15(a) shows a histogram of computed M w values for the events analysed in this study. The Gutenberg-Richter (GR) law for the whole set is displayed in Fig. 15(b). The GR law suggests a b-value of ∼1.1, calculated by linear regression between the completeness magnitude M c = 1.7 and 3.5, but with a deficit in the larger magnitude events that can be interpreted as a variation of the b-value between smaller magnitudes (1.7-2.8) and larger ones (2.8-3.7). The cumulative number of events is displayed in Fig. 16(a), while Fig. 16(b) shows the cumulative seismic moment. The dashed lines in these panels mark the origin times of 11 major events of the swarm (M w ≥ 3.2, marked with asterisk in Table 3).

C H A R A C T E R I S T I C S A N D I N T E R A C T I O N S O F T H E 2 0 1 3 S WA R M
The earthquake swarm of 2013 in Helike was one of the most intense sequences that have occurred in the vicinity of Aigion city during the last decade. The interaction between the north-dipping normal faults in the southern part of the rift and the seismogenic layer on which they are rooting has created a complex network of smallscale faults which were possibly involved in this swarm (Fig. 13). The geometrical and temporal characteristics of the 2013 swarm and their relation to previously activated neighbouring areas are discussed in the following subsections.

Geometrical characteristics
The coarse geometry of the spatial distribution in combination with the normal FMs is in agreement with the known tectonic regime, suggesting a north-dipping normal fault. A possible interpretation is that the main fault that hosted the swarm is the Pirgaki fault, which outcrops about 8 km to the south (Fig. 2). Its strike direction ranges N095 • -100 • while the dip angle varies between 40 • and 70 • to the north (Micarelli et al. 2006). This is consistent with the mideastern portion of the spatial distribution which is aligned in an E-W direction dipping roughly 40 • towards the north. The highresolution analysis which was undertaken in the present study has also provided details on secondary faults which were implicated at the western part of the swarm that was activated at a later temporal phase. There are a couple of spatial clusters (CLIDs: 1 and 6) which deviate from the N100 • E oriented plane and make the spatial distribution appear curved at its western portion, dipping ∼45 •  Table 4. Results of composite focal mechanisms per spatial cluster identified by CLID. Column '# evt'. refers to the total number of events in each spatial cluster, while column '# f.m'. is the number of focal mechanisms that corresponds to the cluster. a For the spatial clusters the depth column refers to median depth. b Some outliers were removed before the spatial clustering procedure was applied. c Excluding group 'X'. , which coincides with the centre of the middle cross-section g 1 -g 2 . The colours and numerical labels correspond to the spatial clusters defined in Figs 6(b) and 7(F). X 2 and X 6 mark the subhorizontal focal mechanisms of group 'X' (Fig. 12) belonging to clusters #2 and #6, respectively. Dashed lines represent possible faults involved in the swarm: P.f., Pirgaki fault; F8, possible antithetic fault corresponding to cluster #8; Fx2 and Fx6: possible subhorizontal faults corresponding to subclusters X 2 and X 6 , respectively.  NNW. The activated fault segment is located about 2 km deeper than the seismically active weak layer of diffuse deformation  which is situated at 6-8 km beneath the gulf (Fig. 17). This result is independent of the velocity model that has been used for the hypocentral locations (Figs 4b and d).
The extension of the fault plane, as indicated by the mean focal mechanism (strike = 261 • , dip = 41 • , Table 4) at 9 km focal depth, meets the surface approximately at the traces marked with FP1 and FP2 in Fig. 2. This includes Kerinitis fault, which links the eń echelon fault segments of Pirgaki and Mamousia. The Kerinitis river runs between the latter as well as the Western and Eastern Helike fault segments. An SW-NE trending fault trace has been mapped, dipping 70 • NW (fault segment 13d described by Ghisetti & Vezzani 2005), while its NE extension is drawn with dashed line in Fig. 1. Its existence is supported by seismological evidence (Pacchiani & Lyon-Caen 2010) and field observations (Ford et al. 2007;Backert et al. 2010). The intense swarm of 2001 near station AIOA (Agios Ioannis village) has been attributed to this fault, which is roughly parallel to the orientation of the Hellenic napes and is characterized as a structure that is cross-cutting the currently active N100 • E trending faults Lambotte et al. 2014). This hypothesis has been based on the geometry of the best-fit plane indicated by the relocated spatial distribution of that swarm (Kapetanidis et al. 2010;Pacchiani & Lyon-Caen 2010), as well as the FM solution of its major event (Zahradnik et al. 2004) and was also suggested to have been triggered by fluids, originating approximately at Kerinitis river (Pacchiani & Lyon-Caen 2010).
Clusters such as 11 and 5 suggest a more E-W trending, northdipping fault at a bit steeper angle, which best fits the easternmost part of Pirgaki fault, although the difference is of the order of 5 • . Subhorizontal FMs have also been observed at the same western region, at depths of about 9.2-9.7 km (subgroups X 2 and X 6 ,  Fig. 13a, profile e 1 -e 2 ). It is safe to assume that there are at least two parallel subhorizontal discontinuities at the westernmost segment, which are depicted by the dashed lines Fx2 and Fx6 in Fig. 13(a), as they are grouped in clusters 2 and 6, respectively. The small-scale subhorizontal fault Fx2, which also hosted the major event of 14 July, likely extends to the next segment, shown in Fig. 13(b) (profile f 1 -f 2 ). However, the focal mechanisms of group X 2 in that crosssection are more disperse, possibly due to location errors (which is why the Fx2 line is drawn to pass through their mean depth), or would, otherwise, suggest the existence of additional subfaults at different depths. The FMs of group 'X' have an average dip of 8 • towards W-WNW, likely occurring on subhorizontal weakenings in the crust, while providing an SSW-NNE slip direction consistent with the regional extension. The FM of the major event of the 2001 swarm is similar (one subvertical WNW-ESE plane), but with larger dip angle on the NW-dipping plane. Regarding a subset of oblique-normal FM solutions which were observed during the 2013 crisis, these could be attributed either to a small antithetic southdipping fault or to an NW-dipping one. This suggests that while the swarm most likely took place on the easternmost, deepest part of Pirgaki fault, the focal zone is probably affected by an interaction either with Kerinitis fault or by other tectonic features. However, the best-fit plane of the 2001 swarm passes above the 2013 hypocentres, suggesting no direct mechanical influence of the Kerinitis fault to the 2013 swarm.
In addition to the 2001 sequence, important seismic activity had occurred earlier, in 1991, in the vicinity of the 2013 swarm. That was an aftershock sequence which followed an M L = 4.5 event on 1991 July 3, described as the Aigion cluster 'Cl1' by Rigo et al. (1996). A relocated subset of its spatial distribution ) is displayed in Fig. 17. Its hypocentres are concentrated at a depth of around 8 km, shifted slightly towards the north (Fig. 2). Rigo et al. (1996) observed that they are diffused through the whole block bounded by Pirgaki and Helike faults. The probability of bias in the original locations caused by the RM velocity model, as one of its discontinuities lies near 8 km (Table 1), is considered to be low, taking into account that similar results were obtained with the CM and, also, after performing double-difference relocation using either model. The divergence from the hypocentral region of the 2013 swarm could also be due to the difference in network geometry or data quality. Another possible explanation could be that the 1991 sequence has occurred on a different (possibly blind) normal fault, parallel to the Pirgaki one. This scenario might be supported by the observation of a small bulge in the topography (Fig. 17). No superficial major fault trace is shown on the map of Fig. 1, however some small fault segments have been reported in the literature (e.g. segment 7 described by Ghisetti & Vezzani 2005). A third possibility would be that the 1991 sequence has occurred on a fault segment with its superficial trace a bit more southerly of Pirgaki fault, with a smaller dip angle. This fault segment intersects with the main Pirgaki fault at about 1 km depth (segment 13b described by Ghisetti & Vezzani 2005), retaining its lower dip angle, which means it should reach the depth of ∼8 km more northerly than Pirgaki fault. The location uncertainties are such that the spatial distribution of the 1991 sequence does not permit the definition of a fault plane dipping north. According to a different scenario, it would have occurred on a small antithetic, south-dipping fault. It is noteworthy that the mean FM solution for the 1991 sequence  is quite similar to the one estimated for the 2013 swarm.

Temporal characteristics
The spatiotemporal distribution of the 2013 swarm revealed that the earthquake activity during the first phase was developed bilaterally in several sudden bursts, with a partial, weak gradual migration towards the east. The parabolic curves drawn in Fig. 9 correspond to the theoretical triggering front r = r(t) (eq. 3) caused by fluids infiltration and pore-pressure migration (Shapiro et al. 1997): where y o and t o represent the source position and time of the occurrence of the fluid injection and D is the hydraulic diffusivity of the medium. The theoretical r(t) lines have been drawn bilaterally for different values of D with a pair of (y o , t o ) parameters estimated by visual inspection of the spatiotemporal diagram for each of the two phases. The y o corresponds approximately to the projection of the major events of 31 May and 14 July, for the first and second phase, respectively, while t o has been chosen ∼1 d before the initiation of the corresponding subsequence. All events triggered by fluid intrusion should occur within the envelope defined by the rupture front, r(t), with the farthest events following close. The eastwards migration of the first phase nearly follows this pattern for D = 0.1 m 2 s −1 , while the western branch is better enveloped by the D = 0.2 m 2 s −1 line. Although the second phase is within the bounds of the triggering front of the first one, it should most likely be considered independently. A fast, bilateral migration during the latter phase, lasting for about 7 d, is implied by the D = 0.4 m 2 s −1 curve in Fig. 9, combined with persistent activity near its initial source.
An alternate view is provided in Fig. 18, which shows the radial distance from the starting point of each phase. The D = 0.4 m 2 s −1 curve fits well on the beginning of the first phase, while D = 0.1 m 2 s −1 fits better its later part, after day 15, at least in terms of migration rate. In the second phase, the D = 0.4 m 2 s −1 curve does not manage to envelope the events of cluster #2 (y ≥ 0.5 km) as they occurred almost simultaneously to those at cluster #1. It does, however, match well the major events of 14 and 15 July which are displayed by stars in Fig. 18. A similar pattern can be observed for the first phase if the D = 0.1 m 2 s −1 line is shifted to the right by ∼5 d. What is common in both phases is that, as in other cases of swarms in the literature (Aoyama et al. 2002;Vidale & Shearer 2006;Roland & McGuire 2009;Legrand et al. 2011), the stronger events appear later in the sequence, rather than the beginning. In many instances they are preceded by a number of foreshocks, while they behave as part of the sequence without triggering large series of aftershocks.
There are significant similarities between the histogram of repeating events per day (Fig. 10, top panel) and the total number of events per day (Fig. 9, top). This was expected, as their majority (80 per cent) has been found to belong to multiplets, although not necessarily large ones. However, the most striking similarity is between the diagram of the cumulative seismic moment (Fig. 16b) and the multiplet history dendrogram (Fig. 10, main panel). This shows that periods of outbursts in the seismic activity were followed by generation of new multiplets, which is interpreted as a spread of the activity in new areas rather than repeated slip in places that were ruptured earlier, even though repeating events tend to occur on a fault patch soon after the generation of a new multiplet. The major event of 2013 May 31 did not seem to cause an adequate number of aftershocks reflecting its magnitude (no sudden increase in Fig. 16a). It was, however, soon followed by events that belong to multiplets which were briefly re-activated during the second phase of the sequence, such as the event of 1 June, 04:36 UTC, which belongs to multiplet #44 with more than 45 repeating earthquakes in phase #2 (Fig. 10). These events belong to spatial cluster #2 and are located within the bounds of the circle marking the epicentral region of period e, at its northern part (Fig. 8e). This implies influence of the seismicity during the first phase in the preparation of the region that was activated in the second one.
The behaviour of Phase #1 is characterized by energy release in several intermittent bursts. Furthermore, it begins smoothly rather than abruptly, and proceeds with a few significant events accompanied by a large number of smaller events which do not contribute significantly to the released seismic moment. The second phase, beginning roughly around 12 July, is more typical of a main-shockaftershock sequence, starting with a large number of events, including some with M w ≥ 3.5, accompanied by strong energy release, then slowly decaying with relatively smaller shocks. The released energy can be described in terms of cumulative moment magnitude (CMM) (Fig. 19). The CMM for the whole sequence reaches 4.55, while the corresponding values are 4.39 and 4.30 for the first and second phase, respectively. Although the total energy is comparable for both phases, the latter released it more abruptly than the former (dashed red line compared to the blue line in Fig. 19). This is mainly explained by the fact that most of the major events of phase #2 occurred in a very short interval of 1-2 d (14-15 July, stars in Fig. 18), while the corresponding events of phase #1 were widespread in a time-span of ∼30 d, which is more typical of swarm-like behaviour. This observation indicates that the first phase, possibly triggered by fluids intrusion, has contributed to the preparation of the more abrupt second phase, which resembles a typical mainshock-aftershock behaviour in terms of temporal energy release, without, however, starting with a major event that could be characterized as a main shock of significantly higher magnitude than the rest of the events in its sequence.

Triggering power of the 2013 swarm
The rather large cumulative seismic moment magnitude of the swarm, close to M = 4.5, can be used to estimate the Coulomb stress change on Helike fault, which is regarded as the most threatening active fault in the area. We discarded the Pirgaki fault, as it is considered inactive for over 0.4 Ma, according to geological studies (Hemelsdaël & Ford 2014). Two alternative fault geometries have been proposed by Lambotte et al. (2014) for the Helike fault, about 5 km away at its closest distance from the 2013 swarm hypocentres. A listric fault, with dip angles changing from 50 • near the surface to 30 • beneath 3 km, reaching the root of the Aigion fault at around 6 km in depth; and a straight fault plane reaching larger depths (10-12 km) at a constant dip 50 • (Fig. 17). Mechanical considerations, involving the fluctuations of seismicity rate and the stability of opening strain rate from GPS, have led Lambotte et al. (2014) to favour the steep, deep Helike fault model. The Coulomb stress change was evaluated for both fault models on 3 km × 3 km subfault patches with the Coulomb 3.3 software (Toda et al. 2011). The result is a Coulomb stress change by less than 0.01 MPa (increase or decrease, depending on the position of the subfaults). The positive changes thus represent a clock advance of less than 1 yr, supposing a recurrence time of 300 yr on the West Helike fault and a 3 MPa stress drop (Briole et al. 2000).

D I S C U S S I O N A N D C O N C L U S I O N S
The seismic crisis of 2013 in Helike has offered the opportunity to investigate the fault geometry and mechanics at depth in a poorly known region of the Western Corinth Rift which has been relatively inactive since 1991. The high-resolution relocation and a multitude of focal mechanism solutions determined in the framework of this study have enabled the detailed description of both spatial and temporal characteristics of this swarm. Its geometry is a slightly curved surface which dips about 50 • N-NW, a result which is also confirmed by the focal mechanisms, consistent with the eastern, downdip continuation of the presently inactive Pirgaki fault. The thickness of the swarm's spatial distribution (∼500 m), which may be considered large relative to the relocation uncertainties, combined with the observation of subhorizontal planes for some of the events in the western part at depths between 9.0 and 9.5 km and, possibly, some antithetic south-dipping faults at 8.0-8.5 km, are indicative of fault branching and, likely, fault edge. We thus propose that this swarm marks the downdip limit of the seismogenic part of the Pirgaki fault.
Clustering and spatio-temporal analysis revealed evidence for triggering by diffusion of pressurized fluids, especially during the first phase of the sequence, which lasted for ∼50 d and was characterized by a series of outbursts with a weak migration to the east. The western portion of the swarm was then activated abruptly, releasing roughly the same amount of energy in a much shorter time. The sudden nature of the activity during the second phase is also reflected by larger migration rate, which is described by a hydraulic diffusivity parameter D = 0.4 m 2 s −1 compared to the D = 0.1 m 2 s −1 value that better describes the eastern migration of seismicity during the first phase. The fast, radial migration during the first day of the swarm's activity (Fig. 18), with a roughly constant rate of 1 km d −1 , suggests a significant contribution of creep propagation and faulting interaction for this early stage. This implies intricate mechanical processes, involving coupling between creep and pore-pressure diffusion, similar to what was revealed by Bourouis & Bernard (2007) in the context of microseismicity induced by fluid injection. During the rest of the first phase, the mean migration rates, roughly 20-100 m d −1 , and the hydraulic diffusivity values, D, ranging from 0.1 to 0.2 m 2 s −1 , are all characteristic of porepressure diffusion. The faster stress relaxation for the second phase may be due to the transfer of coseismic stresses and possibly of the high pore-pressure from the first phase to the western fault system: the latter would thus have been loaded closer to its failure strength, making it more reactive when failure started diffusing within it. This combination of triggering mechanisms (fluid diffusion and tectonic stress transfer) has been previously observed in seismic swarms, such as the 2000 swarm in Vogtland/NW Bohemia (Hainzl & Ogata 2005). The limited extension of the swarm upwards and northwards may be due to the interaction of the Pirgaki fault with others that are intersecting it, like the NW-dipping Kerinitis fault, which may act as barriers for stress migration through creep and/or pore-pressure diffusion. The narrow depth range of the hypocentres can be mainly attributed to their strong clustering and relatively small magnitudes. Similarly narrow ranges, however, are also observed in cases of induced seismicity by man-made hydraulic stimulations (Ake et al. 2005;Baisch et al. 2010;Bachmann et al. 2012;Albaric et al. 2014), which makes the fluid intrusion hypothesis even more likely. The 2013 swarm possibly took place at the deepest edge of the Pirgaki fault where brittle fracture is feasible. At larger depths, the deformation of the crust becomes gradually ductile and, even if the Pirgaki fault extends deeper, it is not expected to be seismogenic. At shallower depths, the Pirgaki fault is likely intersected by the transverse Kerinitis fault, which decouples this seismic activity from the upper layers, either by lateral creep transfer or via the control of pore-pressure migration by impermeable barriers (upwards) or drainage (sidewards). Above this intersection it is probable that the Pirgaki fault becomes simpler, with a stable dip angle, thin gouge and narrow damage zone, and is completely sealed and locked towards the surface, consistent with the absence of activity for more than 0.4 Ma, as reported by geologists (Hemelsdaël & Ford 2014). This is in contrast with the case of the 2001 swarm, during which the microseismicity went a longer way updip on the Kerinitis fault (4-5 km), due to the activation of a rather thick system of subparallel faults, which probably favoured the upward migration of fluids (Pacchiani & Lyon-Caen 2010).
The low value of the calculated static stress perturbation on the nearby Helike fault is independently supported by the absence of any triggered microseismicity in that region during the 2013 swarm. It should be noted that the equivalent M = 4.5 cumulative moment magnitude may be an underestimation of the total released seismic moment, as pore-pressure effects and additional aseismic creep may have possibly increased the global strain at the seismogenic volume. However, as no significant microseismicity has been detected in the central part of the rift, despite its high sensitivity to triggering (Godano et al. 2014), one can exclude a cumulative magnitude of the order of M ∼ = 5. This indicates that even if significant aseismic deformation has occurred it could not have been much larger than the seismic one.
Finally, taking a broader view of the rift mechanics into account, the 8-10 km depth range of this swarm clearly demonstrates the existence of high shear stresses in a necessarily elastic-brittle crust much deeper than the 6-8 km depth range of the main seismically active layer in the centre of the rift. This makes unlikely a dominant role of a shallow north-dipping detachment controlling the western rift opening, favouring a rifting model, as proposed by Lambotte et al. (2014), with a rather immature character of the detachment structure, and a deeper, more centered anelastic strain source.

A C K N O W L E D G E M E N T S
The present study was co-funded through the FP7-ENV-2011/Collaborative Project 'REAKT: Strategies and tools for Real Time EArthquake RisK ReducTion'. We wish to thank the scientists and personnel who participated in the installation or maintenance of the stations belonging to the HUSN and CRL networks. We would also like to thank the two anonymous reviewers for their constructive comments which helped improve this paper. Maps and cross-sections were created using Generic Mapping Tools (Wessel & Smith 1998).

S U P P O RT I N G I N F O R M AT I O N
Additional Supporting Information may be found in the online version of this paper: Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

A P P E N D I X A : O N WAV E F O R M S I M I L A R I T Y A N D F R E Q U E N C Y B A N D S
The degree of waveform similarity within multiplets depends, among other parameters, on the selection of the frequency filter that is applied. For the case of microearthquakes recorded by stations in local distances, Kapetanidis & Papadimitriou (2011) have performed extensive tests with correlation detectors in a wide range of frequency bands. The tests have shown that a frequency range Downloaded from https://academic.oup.com/gji/article/202/3/2044/610654 by CNRS -ISTO user on 17 August 2020 Figure A1. The same waveform recordings as in Fig. 5, but recorded in station TEME, (a) filtered between 2.5 and 23 Hz, (b) filtered between 2.5 and 10 Hz. The waveforms are sorted by the same order as in Fig. 5, regardless of their similarity to the first one at the top. Colours refer to the XC max value for the cross-correlation of each waveform with the one at the top of the stack (indicated as XC). between 2.5 and 23 Hz is usually adequate for a default filter. The 2.5 Hz corner frequency of the high-pass manages to remove the long-period trend, which can be significant in broad-band stations, while the low-pass filter at 23 Hz preserves enough signal complexity in the waveforms to enhance the distinction between neighbouring multiplets. Similar frequency bands, permitting frequencies higher than 10 Hz, have also been used for cross-correlation by other authors in the literature (Scarfì et al. 2003;Carmona et al. 2009;Häge & Joswig 2009;Chun et al. 2010;Yukutake et al. 2010). It is important that the filter is focused in the frequency range where the signals are coherent (Kita et al. 2010), include their dominant frequency (Massin et al. 2013) and also aim to improve or retain adequate signal-to-noise ratio (Hemmann et al. 2003;Massa 2010). A wider frequency range is also useful to prevent cycle-skipping (Akuhara & Mochizuki 2014). However, frequencies below 10 Hz are preferred in other studies (Shearer 1998;Massa 2010;Myhill et al. 2011;Massin et al. 2013;Pirli et al. 2013;Vavryčuk et al. 2013), which may depend on factors such as the type of earthquakes (e.g. volcanotectonic), the longer epicentral distances or, simply, the requirement for less waveform complexity to achieve higher correlation coefficients.
In the framework of the present study, multiplets were constructed by applying nearest neighbour linkage on a cross-correlation matrix, which combined measurements of three-component recordings from two stations (LAKA and TEME), with a threshold C th = 0.8. An example of similar waveforms recorded at the vertical component of station LAKA, situated at an epicentral distance of ∼10 km, is presented in Fig. 5. While, in this case, the lowest XC max value is 0.795 (the bottom waveform compared against the top one), there can be many events in the same multiplet with their waveforms having an XC max value well below the C th , when compared against the top waveform of Fig. 5. This is due to the selection of the linkage method, which only implies that there is always a chain of intermediate events linking an event A with another event B, in which successive pairs satisfy the condition XC max > C th , while this may not be true when A is directly compared to B. Even so, this applies to the combined cross-correlation matrix which was used for the linkage and not to a single component of a certain station.
The same events as in Fig. 5, as recorded by station TEME, which is situated within the epicentral region, are presented in Fig.  A1(a), in the same order and using the same filter (2.5-23 Hz). It is worth noting that, in this case, the degree of similarity between individual waveforms and the top one is significantly different than in Fig. 5. Although the station is very close to the epicentres, thus the signal-to-noise ratio is expected to be larger than in LAKA, many XC max values are below 0.7. There are, also, significant variations in the amplitudes of the S waves which in some cases are largely responsible for the reduction of XC max . The low correlation values are not due to noise but can be mainly attributed to changes in the radiation pattern caused by variations either of focal mechanisms (or the relative position of the station on the stereo-net) or in the hypocentral locations, especially in focal depth. The former is more significant for stations such as TEME, within the epicentral region and most likely near the P-axis (Fig. C1), than for stations in longer distances. The latter may introduce extra phases in the waveform train due to reflections/refractions on small-scale discontinuities in the crustal structure. This situation is partially improved, similarity-wise, when a stronger filter is applied between 2.5 and 10 Hz, removing the higher frequencies in the 10-23 Hz band (Fig.  A1b). While the radiation pattern effects still remain (polarities and S/P amplitude ratios) and the time-lags, t m , are unchanged, the waveform complexity is reduced, leading to increased XC max values. This increases the tolerance and enables more events to be grouped in the same multiplet, which may be necessary when similar events Downloaded from https://academic.oup.com/gji/article/202/3/2044/610654 by CNRS -ISTO user on 17 August 2020 Figure B1. 95 per cent c.i. error ellipses for the relocation procedure, (a) for the influence of stations, (b) for the influence of the initial locations (relocation robustness test) and (c) for the influence of each event on the others. are scarce. However, the wider 2.5-23 Hz band is useful when the enhancement of dissimilarities between neighbouring multiplets is required, as in the present study. This allows for the distinction of several subclusters within the swarm, which would have been merged if a narrower band was used.

A P P E N D I X B : E VA L UAT I O N O F R E L O C AT I O N U N C E RTA I N T I E S
Statistical resampling methods were applied for the calculation of approximate location uncertainties of the relocated events, as those reported by HypoDD for the conjugate gradients method (LSQR) tend to be underestimated (Waldhauser 2001). A jack-knife technique was used to investigate the influence of individual stations to the relocation procedure. This works by removing differential traveltime data from one station at a time and repeating the relocation procedure using, otherwise, the same parameters and starting model (initial hypocentres). This was applied only for stations with data covering at least 5 per cent of the whole set. For each case, the distance between the hypocentre derived by the original relocation (with full data) and the one with modified data were calculated, generating values dx, dy and dz for the two horizontal dimensions and the vertical one, respectively. The 95 per cent confidence interval (c.i.) 2-D horizontal ellipses were calculated for each event to define the orientation and length of the major and minor semi-axes (Fig. B1a). Likewise, the 95 per cent c.i. was also calculated for the vertical dimension as measure of the influence in the focal depth.
A bootstrap technique was used for the measurement of the robustness of the relocation procedure, by estimating the influence of the starting model (hypocentres) to the results. This was performed by generating random permutations to the initial hypocentre of each event in all three dimensions and re-running the relocation procedure. The randomness is constrained within the ERX, ERY and ERZ bounds for each event with respect to the original location and the procedure was repeated 1000 times. The 95 per cent c.i. error ellipses are presented in Fig. B1(b). The influence of each event in the relocation of others was also estimated by a similar jackknife method, removing one-event-at-a-time and repeating the relocation procedure. The 95 per cent c.i. error ellipses were also calculated for this case and are displayed in Fig. B1(c).
Error statistics averaged over all solutions are displayed in Table 2. The 95 per cent c.i. values for dx, dy and dz are derived from the distribution of the offset in each direction (1 degree of freedom) while the length of the error ellipse's major semi-axis is calculated from the distribution of the [dx, dy] horizontal vectors (2 degrees of freedom).
The mean uncertainties for the relocated hypocentres are around 35 and 55 m, for the horizontal and vertical axis, respectively, while the corresponding median values are smaller by ∼35 per cent. The major semi-axis length is an upper limit to the relocation uncertainty, largely overestimated due to outliers, especially in the robustness test. It is noteworthy that the events which are concentrated in spatial clusters tend to have the smaller relative location errors (Figs B1a  and b).
The results show that the most influential parameter is the starting model (initial locations of events). This indicates that it is of utmost importance that the original solutions be of adequate quality. The error ellipses of the robustness test strongly depend on the level of noise generated by the random permutations, which in this case are large enough to make the clustering harder to achieve. However, the relocation procedure could also work with noisy initial data, but it would probably require a different configuration to produce the best results.
The stations' influence is slightly less significant (Fig. B1a). This suggests that geometry of the network is adequate enough for the relocation procedure to be stable even if some stations are missing. An important result from the error statistics is that, on average, the errors are uniformly distributed in various directions, which indicates that the alignment of the relocated foci is not biased towards a certain azimuth. The fact that in practice not all events have the same number of arrival-time data in all the stations of the network creates pseudo-random scattering caused by variation in the effective network geometry (Richards-Dinger & Shearer 2000). The doubledifference relocation that has been applied in the present study is proven robust against this type of noise.
The influence of a single event to the others is almost undetectable in all but a few weakly correlated events. This is expected, as the large number of events provides a significant amount of neighbours to each of them, weakening the influence of a single event to the others. In any case, the strongly linked events tend to concentrate into spatial clusters despite the perturbations in the starting model or the available station/event data. The general image remains the same and it is mainly some events with few or without cross-correlation data which are more susceptible to noise.

A P P E N D I X C : S -WAV E P O L A R I Z AT I O N A N D S T O P A M P L I T U D E R AT I O A S W E I G H T I N G FA C T O R S
It is common practice to use first motion polarities (FMPs) to determine focal mechanisms (FMs) of small magnitude events when there are an adequate number of available local stations (Fig. C1a). However, the network geometry may not always be proper to provide constrained solutions with FMPs exclusively. This is the situation in the Western Corinth Rift, where stations are located in opposite sides of the gulf with relatively large gaps east and west for epicentres inside the gulf. In the present study, besides these gaps in azimuthal coverage there is also a lack of stations in 'middle' distances, while there are enough stations within the epicentral region Figure C1. (a) A well-constrained focal mechanism calculated by averaging individual solutions derived from a grid search using only FMPs; (b) a poorly constrained solution with FMPs alone; (c) indication of SWP (straight lines on stations) and SPR (circles centred on stations) for the average solution of panel (b), where the blue and red straight lines correspond to observed and theoretical SWP directions, respectively, while solid and dashed circles correspond to observed and theoretical SPR, respectively. The radius of each circle is proportional to the absolute value of the logarithm of the corresponding SPR, while blue and red indicate positive and negative logarithms, respectively, and (d) modified average solution, weighted by SWP and SPR differences. and far from the epicentres, at the opposite coast. In this case, as most events exhibit a normal FM, it is difficult to constrain the strike and dip of the nodal planes with FMPs alone.
For this reason, instead of simply considering the average solution of FMs which are in agreement with the FMPs, the S-wave polarizations (SWP) and S to P amplitude ratios (SPR) were also taken into account and used as weights during the averaging procedure for each event. SWP of the source were measured for several events within the shear-wave window of the 4 closest stations with three working components (AIOA, TEME, ALIK and HELI) after applying corrections for anisotropy (procedure described in Kaviris et al. 2015). The results were then used for the calibration of an automatic procedure which calculates polarization directions after applying an appropriate bandpass filter (ranging between 0.5 and 5 Hz) which was determined individually for each station, to minimize the effect of shear-wave splitting in the estimation of the polarization direction. The SWP directions were calculated by singular value decomposition on the horizontal particle motion of the S-wave onset in the four selected stations for the same events in which FMPs were also determined. In addition, S and P amplitudes were also determined automatically by measuring the maximum vector length of the 3-D particle motion in windows containing the first three semi-periods of P and S waves in several stations. The logarithm of the observed SPR was considered in each measurement. The theoretical SWP and SPR were calculated using the radiation pattern of S and P waves (Aki & Richards 1980) for each individual solution which was in agreement with the FMPs for each event. The differences between theoretical and measured SWP and SPR were then incorporated in a weighting scheme of the following form (eq. C1): where the index i refers to an individual fault plane solution which is in agreement with the FMPs,¯ SWPi ∈ [0, 90 • ] is the average of the minimum angular differences between the measured and calculated directions of S-wave polarization, SWP , in all selected stations, SPRi > 0 is the absolute difference between the logarithms of observed and calculated S to P amplitude ratio, SPR , averaged in all selected stations, and C o , C 1 , C 2 , P 1 , P 2 , P 3 are positive constants which can be calibrated to control the effect of SWP and SPR to the averaging of individual solutions. In the present application of this method, all constants were set to unity except for P 3 = 1.5. In practice, the applied scheme has been set to allow stronger weights on the SPR than SWP, as there are larger uncertainties in measuring S-wave polarization directions (partially due to anisotropy effects) than S to P amplitude ratios. After a first run of the procedure, the mean differences between theoretical (for the temporary average solution) and observed SPR were calculated in log scale for each station and a correction factor was determined to account for site effects. Then the correction was applied to the observed SPR measurements in a final run of the algorithm.
An example is presented in Fig. C1(b). A 5 • -step grid search for strike (φ), dip (δ) and rake (λ) was performed to identify ∼910 individual fault plane solutions which satisfy the observed FMPs. The average solution (φ = 256 • , δ = 44 • , λ = −101 • ) is calculated by summing the individual normalized moment tensors. In this case, the mean values of SWP and SPR are 40 • and 0.45, respectively (Fig. C1c). The quality of an average solution derived by the gridsearch on FMPs is defined by the number of valid individual solutions, which in this case is very large, permitting a wide range of strike and dip values for both nodal planes, thus the quality is relatively poor.
By taking into account the proposed weighting scheme, the normalized moment tensor of each individual solution is multiplied by the corresponding weight, W i , before summing up. This resulted in a different average solution (φ = 292 • , δ = 54 • , λ = −59 • ), which is also in agreement with the FMPs but is biased towards reducing the mean SWP and SPR differences to 7 • and 0.41, respectively (Fig. C1d). This is particularly useful when the focal mechanisms cannot be constrained adequately by FMPs alone. In this study, at least 1 and an average between 2 and 3 SWPs were used for each event while SPR were calculated in 10 selected stations.