The 2020 Samos (Aegean Sea) M7 earthquake: a normal fault with rupture directivity and near surface slip explaining the tsunami generation and coastal uplift

The 30 October 2020 M7 Samos earthquake occurred offshore the Greece-Turkey cross border region, and will be recalled as among the deadliest (118 fatalities) that affected both countries. It generated a strong tsunami and caused coseismic uplift of 20 to 35 cm of the NW part of the Samos Island. It ruptured a ~60 km long, north-dipping normal fault, related to the back-arc extension of the Aegean Sea area. Using picks from regional strong motion and broad-band waveforms we relocated the mainshock and the aftershocks, applying suitable velocity models. The closest strong-motion recordings, constrained the finite fault slip model, suggesting up-dip and westward propagation of the rupture. The westward rupture propagation is independently confirmed by the apparent source time functions inferred using the empirical Green’s function method from nearregional broad-band and strong-motion waveforms. Static displacements measured by GNSS stations constrain near surface slip of ~1 m, explaining the tsunami and the island uplift. The 2020 Samos event showed that normal faults bounding the basins in the back-arc Aegean region can host M7 earthquakes and when combined with tsunami generation, constitute a constant threat for the nearby coastal areas of both Greece and Turkey.


Introduction
On 30 October 2020 11:51UTC (13:51local time) an earthquake of moment magnitude M7.0 (GCMT) occurred ~ 9 km offshore the northern coast of Samos Island (Greece) in the 690 m deep Samos Basin, at eastern Aegean Sea (Fig. 1). This event is among the strongest events to occur during instrumental times in the Aegean Sea. It will be recalled in the collective memory as the event that caused severe damage ~70 km away from the epicenter, namely to the metropolitan area of Izmir (Turkey), a city of ~4M inhabitants. It caused 2 fatalities and 19 minor injuries at Samos Island, but 116 fatalities and over 1,030 injuries in Izmir. The damage in Izmir Bay area (Bayrakli district) was mainly attributed to the amplification of ground shaking at site frequencies in the range 0.7 to 1.6 Hz for both stiff and soft soil sites (Cetin et al., 2020).
This earthquake is yet another event along the Turkey-Greece cross border region, severely affecting both nations. Detailed information regarding the seismic history of the region is provided in Cetin et al. (2020). Dominant E-W striking normal faults control the northern Samos coastline, the marginal fault of the Samos Basin, is a structure of ~35 km length and of 20-30 km width (Nomikou et al., 2021), inferred from sea topography. It is this fault, e.g. the Samos Basin Fault (Nomikou et al., 2021) or else catalogued as Kaystrios Fault (Caputo and Pavlides, 2013) or as North Samos Fault (Chatzipetros et al., 2013) that the sequence was associated to. The Samos earthquake predominantly ruptured the western segment of this fault, that had not ruptured in documented historical times (since ~1700) while a number of strong events (1873, M~6.5) may be associated with the eastern segment of this fault (Cetin et al., 2020; see chapter 1).
The dislocation of the seafloor produced a tsunami, with maximum inundation (e.g. max horizontal intrusion) and maximum runup equal to 2.31m and 3.82m, respectively, which affected nearby Samos Island and cities along the Aegean coast of Turkey, resulting in substantial property losses. Field survey on Samos Island indicated that sea recession was the leading motion (Triantafyllou et al., 2021). More specifically, at a location 5.8 km east of Karlovasi (site K4 in Triantafylllou et al., 2021), sea recession was followed by waves of ~3.35m height which occurred 2 min and 4 min after the earthquake's origin time, respectively. Moreover, at this site the authors suggest that coupling of tsunami with landslide generated tsunami might have taken place. In general, the tsunami arrived within ~10 min to the coast of NW Samos Island and within 20 min to the coast of Turkey. This short time window halved any immediate protection measures, however a tsunami warning alert was disseminated to the residents of the Greek Aegean Sea Islands at 12:15 UTC (Triantafyllou et al., 2021). Field reconnaissance measured co-seismic shorelines uplift of 20±5 to 35±5 cm at the western coastline of northern Samos (Evelpidou et al., 2021), and a tectonic uplift of 10 cm was detected from geodetic data on Samos Island (Ganas et al., 2020).
At a broader scale the tectonic setting is governed by i) the rapid (~24 mm/y) westwards escape of the Anatolia block towards the Aegean (escape-tectonics), and ii) the even faster (~35 mm/yr) southward migration (trench retreat) of the Hellenic subduction zone, due to the slab rollback (Jolivet et al., 2013, Faccenna et al., 2014. The westward extrusion of Anatolia relative to Eurasia is partly driven by the Arabia-Eurasia convergence and partly by gravity, e.g. the forces originating from gravitational potential energy differences between the high topography of eastern Anatolia and the low elevations of the Aegean Sea (Karabulut et al., 2019), but it is still unresolved which driving force prevails. From a kinematic point of view, the western motion of Anatolia is facilitated by the dextral strike-slip North Anatolian Fault (NAF), in conjunction with the operation of the sinistral strike-slip East Anatolian Fault (EAF). The trench retreat drives the widespread ~N-S extension and crustal thinning, as evidenced by many active grabens in the back-arc Aegean region and western Anatolia.
The above briefly described tectonic regime is reflected in the earthquake focal mechanisms.
Along the cross-border Greece-Turkey region, a mixed mode of normal and strike-slip faulting exists, due to the back-arc extension and the shear imposed from the westward escape of the Anatolian block (Chatzipetros et al., 2013;Kiratzi, 2002Kiratzi, , 2014Tan et al., 2014). The North Anatolian Fault as it enters into the Aegean Sea, broadens the shear zones, as it splays into enechelon NE-SW dextral strike-slip faults and associated grabens (pull-apart structures), The Samos sequence occurred in a location with good station coverage from regional networks of Greece and Turkey. We harness abundant broadband (BB) and strong motion (SM) digital recordings and GNSS static displacements to: i) relocate hypocenters of the 2020 Samos sequence and benchmark the efficacy of available velocity models, ii) seek for source complexity by studying the rupture history, and its kinematics, and iii) discuss the results within the regional seismotectonic context. We provide a consistent finite -fault model, explaining both strong motion and GNSS data. We document a shallow slip localized close to the northern-most coast of Samos and a westward rupture directivity towards the Greek mainland. The along-strike rupture directivity for normal faults has been documented for the normal-faulting earthquakes along the Apennines in Italy (Ren et al., 2017;Calderoni et al., 2015Calderoni et al., , 2017Pizzi et al., 2017;Tinti et al., 2016;Ameri et al., 2012), and our results come as another example for directivity related to normal faulting within the Eastern Mediterranean.

Data and methods
For the subsequent analysis, we used waveform data, either strong motion records or broadband records or a combination of both, from networks operated in Greece and Turkey. For each section of our analysis, we have included in the Supplement additional sections and figures with detailed description of the methods adopted, the data used and their application-specific processing. In the following, we focus on the key elements and results of our work.

Mainshock hypocenter location
We manually picked P-and S-wave arrivals from strong motion (SM) and broadband (BB) stations, 81 in total, at epicentral distances from 20 to 290 km (Fig. S1). Several 1-D velocity models (VM) applicable to the region are available which allowed us to benchmark their efficacy to provide consistent results regarding the mainshock hypocenter, and get an estimate of the uncertainties. We examined 8 velocity models and located the mainshock using the probabilistic NonLinLoc code of Lomax et al., 2000 (see Section 1A, and Figs S1, S2 in the Supplement for details). The location of the epicenter is very stable, at 37.900N and 26.817E (±2.5 km) within all the models tested ( Fig. 2 and Table 1). As a rule, depth is the most difficult parameter to constrain, however all models prefer depths in the upper crust, at ~12 km and shallower. The models of Konstantinou (2018) and Novotný et al (2001) provide comparable P-and S-residuals ±1s at all epicentral distances (Fig. S2). The model of Novotný et al (2001), which was derived from surface-waves dispersion, has proven to be very efficient also in describing wave propagation in the Aegean area, especially at low frequencies used in the finite-fault modeling.

Relocation of aftershocks
From all the velocity models examined regarding their location uncertainties and their hypocenter distribution (see Section 1B and Figs S3 to S6 in the Supplement for details), we selected the Özer et al. (2018) as our preferred result. This model was derived from the most recent seismic experiment in the study area, and provided the best travel-time data fit and better hypocenter distributions without artificial linear depth concentrations. Figure 3 summarizes the final HypoDD (Waldhauser 2001) relocated dataset of ~1500 events in total. The hypocenters have -mean-square residuals of 0.03s. Formal estimates of the uncertainties in NS, EW, depth and of the origin time are 0.07km, 0.05km, 0.1km and 11ms, respectively. From the aftershock's spatial distribution it is evident that the area west of the hypocenter is depleted in aftershock productivity compared to the eastern area (see also Figs. S6 in the supplement). This observation is a first proxy for the inferred locus of the major slip, later confirmed by our slip model. The cross sections (inset in Fig. 3) clearly depict that: i) the sequence operated in the upper crust (h<15 km) and the aftershocks are mostly updip from the hypocenter ii) the causative fault dips to the north, and iii) secondary structures were activated. The main cluster comprising the hypocenter, borders the Samos northern coastline, and all aftershock mechanisms associated with it clearly depict pure normal faulting along E-W striking planes (Karakostas et al., 2021).
Whereas, the westernmost cluster with a NE-SW alignment following the Ikaria basin topography, (Nomikou et al., 2021)

Multiple point source (MPS) modeling of the mainshock
To parameterize the finite-fault model, we initially investigated the source process of the 2020 mainshock using the tools available in the ISOLA software package (Zahradník and Sokos 2018; Liu and Zahradník, 2020 and references therein), adopting the velocity structure of Novotný et al (2001) to calculate full-waveform synthetics. We specifically seek to approximately describe the rupture process as a sequence of points of moment release episodes, known as subevents. Our goal is to identify the number of subevents involved in the rupture process, their focal mechanism parameters, and their relative space and time separation.
We started with usage of waveforms from BB, stations in the epicentral distance range of 263 to 456 km and their centroid moment tensor (CMT) inversion in the frequency range between 0.005 Hz and Fmax, varying Fmax from 0.02, 0.03 and to 0.04 Hz. For the lowest Fmax value, we found a poor spatial resolution of centroid and an overestimated moment magnitude (Mw = 7.04, compared to Mw = 7 of GCMT). For the other values of Fmax we obtained correct Mw (7.0 and 6.9, respectively), satisfactory waveform fit (variance reduction VR=0.8 and 0.6) and a consistent position of the centroid, shifted 20-km westward of the epicenter, and with the centroid time ~10s relative to origin time. The latter can be taken as a proxy of the half-duration, as later confirmed by our more detailed modeling. The optimal C-depth was between 6 and 8 km, representing an improvement compared to the fixed 12-km depth of GCMT. Besides, we found the C-depth almost independent on the used velocity model, being significantly more stable than the hypocenter depth ( Figure 2). Regarding the focal mechanism we found a stable, high percentage double-couple deviatoric source (DC>85%), with strike/dip/rake (s/d/r) angles equal to 270/50/-100, slip vector azimuth/plunge= 15/49, close to the GCMT solution (Table 2).
Increasing the frequency range from 0.04 to 0.09 Hz and using SM stations in the distance range 30 to 155 km, the earthquake appears as a multiple point source. To stabilize the inversion, we kept the depth fixed at 6 km and searched for best-fitting positions of point sources along a W-E striking horizontal line. Using double-couple constrained inversion, we identified three relatively stable subevents, shown in Figure 4, alongside our subsequently discussed slip model. Their focal mechanisms denote normal faulting, featuring some variation among them. Their tensor sum yields s/d/r equal to 246/44/-125, slip vector azimuth/plunge = 20/35, being close to the GCMT solution (Kagan angle 21), Mw=7.0 and VR=0.63. Further tests (not shown here) proved a fourth subevent as unstable, providing only negligible improvement to the waveform fit.
The first subevent is situated 4 km east of the epicenter, indicating an initial eastward propagation. The other two subevents, located at 12 and 32 km west from the epicenter, support the predominant westward propagation, further confirmed by our modeling. The subevents occur at 6, 10 and 15s after origin time, providing ~20s of source duration in agreement with the source time function provided by Geoscope, using the SCARDEC method (Vallée et al., 2011). Using the subevents' source mechanisms and the non-negative-least squares (NNLS) method of Lawson and Hanson (1974), we calculated the moment rate function (MRF); see appendix of Zahradník and Sokos (2014). This MRF is compared with the one obtained from the slip inversion model (Fig. 4). Both functions depict a total source duration of ~ 20s.

Finite -Fault kinematic rupture model
We used the Linear Slip Inversion (LSI) method of Gallovič et al. (2015) to infer kinematic finite-fault description of the rupture process. The method has been applied to various earthquakes in Greece (Sokos et al., , 2016(Sokos et al., , 2020. Table 3 lists the quantities describing the setting of the calculation. In LSI, slip rate functions, spanning the full rupture duration, are discretized in time and space. Synthetic Green's functions are calculated by the discrete wave number method adopting the Novotný et al. (2001) crustal model (as in ISOLA analysis) in frequency range 0.02 to 0.15 Hz. At higher frequencies, the details in source and propagation medium (e.g., wave scattering due to local heterogeneities) could not be adequately captured by our 1D velocity model. The data are displacement waveforms acquired from local strong motion stations ( Fig. 4) filtered in the same way as the Green's functions, using 4th order causal (singlepass) Butterworth filter. We also employ static GNSS data adopted from Ganas et al. (2021).
We stabilize the inverse problem by i) assuming prior covariance function with k -2 decay at large wavenumbers k, ii) prescribing seismic moment inferred by the CMT inversion using ISOLA modules, and iii) positivity of the slip rates. Regarding the latter, we use the NNLS approach of Lawson and Hanson (1974). We point out that the source description in the LSI method is very general with no a-priori constraints on the position of the nucleation point, rupture speed, and shape of slip-rate functions. A drawback of this loose parameterization is that the result is sensitive to artifacts and biases imposed by the imperfect station distribution and smoothing. To this end, the result must be carefully interpreted considering lessons learned from previous synthetic tests and real-data applications (Gallovič, 2016;Gallovič et al., 2015;Gallovič and Zahradník, 2011).
The fault is modelled as a rectangle, 100 km × 24 km along strike and dip, respectively. We gridsearched its position and mechanism: strike (240°, 250°, 260°, 265°, 270°), dip (35°, 40°, 45°) and rake (-100°, -110°, -120°, -130°), position in the north-south direction was varied by ±5 km with respect to the centroid. The best waveform fit (waveform variance reduction VR=0.618 and GNSS VR=0.978, Fig. 5) and the least occurrence of artifacts, was attained for the initial fault position (with its center fixed in the centroid) and for strike/dip/rake=265°/40°/-110°. The preferred slip model (Fig. 4) shows that the main slip episode occurred west of the epicenter. The slip is located both at depth, but also close to the surface. We point out that the shallow slip is illuminated almost exclusively by the GNSS data. Indeed, as we show in the supplementary Fig. S7, the surface slip is not revealed and the closest GNSS is underestimated by about 50% if the GNSS data are neglected in the inversion. Noting that the waveform improves only slightly when GNSS are omitted (VR=0.623), the seismic data proves insensitive to the temporal evolution of the shallow slip, at least in our frequency range and with our station coverage.
Snapshots of the inferred slip rates (Fig. 5) suggest that the rupture started close to the epicenter (not a priori prescribed in the inversion). On the dipping fault, the hypocenter lies at 6 km depth. and propagated bilaterally mainly up-dip for the initial ~3 to 4 s (see the first subevent from the MPS inversion delayed by 6 s after the origin in Fig. 4a, and the slip-rate peaks at 4 to 8 s in the LSI snapshots of Fig. 4d). After that, the rupture continued to the west, i.e. towards the mainland Greece, creating the major slip within 8 to 16 s after the origin time, extending ~40 km westwards of the hypocenter and to shallow depths (see the second and third subevent in the MPS model and the dominant slip patches in the EGF and LSI models). We note that since the shallow slip is constrained just by the GNSS data with almost no effect on the waveform fit ( Fig.   S7), we consider the temporal evolution of the surface rupture as poorly constrained. We hypothesize that at shallow depths the rupture had relatively long slip rate duration due to long rise times and/or slow rupture propagation, with small slip-rate peaks and thus weak radiation of seismic waves, especially at periods dominating the displacement waveforms (~10 to 20s).
For this model parameterization, the rupture speed can be only roughly estimated, keeping in mind that only very smoothed image has been revealed . For the dominant westward faulting, the speed can be estimated from the slip-rate peak position at 45 km along strike at 8s and the termination of the rupture 8s later occurring at about 70 km (see the snapshots in Fig. 5). This suggests relatively high rupture speed of about 3 km/s. Obviously, the possibility that during each moment release episode the rupture could have propagated slower/faster or more episodically cannot be ruled out.

Apparent Source Time Functions (ASTFs)
Rupture directivity is a key element of the physics of earthquakes and here we seek to explore this feature employing the empirical Green's functions (EGF) approach. Several methods were developed to obtain ASTF (e.g., Mueller, 1985;Mori and Hartzell, 1990;Bertero et al., 1997;Courboulex et. al., 1999;McGuire, 2004;Vallée, 2004;Roumelioti et al., 2009;López-Comino and Cesca, 2018). Most of them are based on spectral deconvolution, requiring careful stabilization. Here we suggest a simple alternative technique based on NNLS technique, fully operating in time domain, assuming the ASTF is implicitly positive and seismic moment is constant across the stations (see supporting Section 3 and Figs. S8-S12).
We calculate apparent source time functions (ASTFs) from regional waveform data by the NNLS technique (see supporting Section 3 and Figs. S8-S12) and investigate their duration and amplitude variation with azimuth. This method only requires to find an aftershock to serve as an empirical EGF, originating at similar depth and location with the mainshock and having a similar focal mechanism. No further assumptions are made (no source or velocity model is needed). The aftershock sequence was depleted in strong aftershocks, specifically in the magnitude range M5 to M6. This narrowed significantly the number of suitable EGF's. Exploiting the available data we finally selected the M5 of 31 October 2020 05:31 UTC. At each station, we invert the full seismogram, including P and S waves and all three components. We have also calculated ASTFs using the P or S waves groups only, but the results were similar to the full records inversion. The frequency band of inversion was chosen at 0.05 to 0.5 Hz but higher frequencies can be used also. We have inverted the original acceleration, or velocity or integrated records, and the results were stable, that is to say the basic character of ASTFs was not sensitive to the input data. The ASTFs are searched in a time interval from -5s to 35s relative to origin time, i.e. 40s in total.
The inferred ASTFs (Fig. 6) from stations located orthogonal to the fault strike (EFSA, PRK, SOMA, TVSB in the north and KLNA, ASTA, ARG in the south) depict longer pulse duration and lower amplitudes, compared to those located along strike (KARY, VLY, TNSA in Greece and NAZL in Turkey), which supports westward rupture propagation. More specifically, NAZL lies in the backward direction of rupture propagation, whereas KARY, VLY and TNSA in the forward direction, exhibiting narrow, high-amplitude pulses.
Assuming a horizontal rupture propagation featuring a unilateral rupture propagation on a part of the fault, apparent duration τ ( f ) as a function of station azimuth f can be described by (1) Here T D =T 1 +T 2 is the total rupture duration, T 1 is the rupture duration corresponding to nondirective part of the fault, and T 2 =L 2 /V R is the rupture duration of the fault portion L 2 with assumed unilateral rupture propagation at speed V R . TheV P ,S is the velocity of P or S waves, and α is the rupture directivity azimuth.
We have tested several combinations of T D and L 2 V P , S to find the optimum ones that provide the best match with the observed durations of the ASTFs (Figure 6b). The direction of rupture α = N265 ˚ is fixed, obtained from our fault slip model. The best fit to the data is for values:T D = 22s ±2s, L 2 V P , S =7 s(green curve in Figure 6b). If we consider the S waves velocity in the source depth (e. g., V s = 3.5km/s in Novotny's model), the length of the directive zone is L 2 = 24.5km, which corresponds to the estimate from the kinematic finite-fault modeling. We note that the rupture velocity V R cannot be inferred from the durations of the ASTFs (Eq. (1)).

Conclusions and Discussion
The 30 October 2020 M7 Samos earthquake ruptured a north-dipping normal fault bounding the northern coast of Samos island. The event caused 116 fatalities in the city of Izmir, 70 km away from the epicenter, and 2 fatalities on Samos Island. It will be recalled in the collective memory as among the deadliest events to affect the cross-border Greece -Turkey region in modern times.
It generated a strong tsunami, the largest in the Aegean Sea since the one associated with the Amorgos 1956 M7.5 earthquake (Triantafyllou et al. 2021). The prevalence of E-W striking normal faulting declares that this earthquake was the immediate consequence of the N-S stretching of the upper Aegean Sea plate from mechanisms that combine the effects of the Hellenic subduction, e.g. the slab roll-back and trench retreat. At a broader scale this prevailing extension is accommodated by the opening of parallel oriented grabens in Aegean Sea and western Anatolia as for example the Gediz, Izmir, Samos and Menderes grabens. Evidently, these basins spread within the Aegean Sea and can encompass ~M7 earthquakes, increasing the hazard to nearby islands and coastal areas. As a concluding remark, we want to stress the hazard imposed to the urban regions across the cross-border Greece -Turkey region, as documented from the on-going time-cluster of strong events within the last 6 years or so, encompassing two ~M7 earthquakes (e.g. the Samothrace 2014 and the Samos 2020). Additionally to strong events, many swarms along the cross-border region have also occurred pointing to an increased level of localized stresses across the region. To this end, both countries should join efforts to better understand the connection and continuation of the structures, especially in the offshore areas, that are less explored.

Data and Resources
Digital seismic waveforms were retrieved from the ORFEUS Eida-nodes ( orfeus-eu.org ) , and  Tables   Table 1 Mainshock hypocenter parameters as obtained from P-and S-wave arrivals benchmarking different velocity models within NonLinLoc code. 'Depth' is the formal, best-fitting value. For uncertainty, see Figure 2.The location provided by EMSC is listed for comparison.

Velocity Model
Origin Time   Table 3 Geometry     We relocated the hypocenter of the mainshock, using phase arrivals from broadband and strong motion stations from regional networks (see Data and Resources). We handpicked 81 P-and 40 S-wave arrival times at 81 stations, at epicentral distances from 20 to 290 km (Fig. S1). Using a subset up to 160 km, (circle in Fig. S1) we obtained similar results. We excluded any S-phase arrivals with large residuals or unclear arrivals. We did care though to include S-phases from the nearest stations to constrain the depth. We selected eight velocity models (VM) as most suitable for this application (Fig. 2a,    of Turkey (Fig. S3). The manually picked events recorded during the first 40 days were jointly used for the initial location. The Vp/Vs ratio was set equal to 1.69 based on the Wadati diagram results (Fig. S4). As previously stated we examined eight velocity models: Akyol et al., 2006, Kalafat, 1987, Kayapack&Gokkaya, 2012, Konstantinou, 2018, Novotný et al., 2001, Pasyanos, 2004and Özer et al., 2018 (   October to 4 December 2020. The star symbol represents the relocated mainshock.
Section 2. Finite -Fault kinematic rupture model Here, we show the results of the LinSlipInv kinematic slip inversion without including the GNSS geodetic data. Note the deficiency of the slip model in the shallow slop patch that reached the surface and the misfit on the predicted synthetic displacement on the Samos inland station SAMO, located in the area where the maximum uplift was measured and where the tsunami initiated. The localized shallow slip in our preferred model (Fig. 4) is crucial in fitting the observed 35.7 cm static displacement in SAMO. Figure S7. a) Map view of slip distribution when the geodetic displacements are not included in the inversion. Note the misfit between the GNSS daily solution (blue) and the predicted synthetic displacements (red) at SAMO. b) Moment rate functions. c) Slip distribution (with slip rate functions superimposed (maximum scaled to 0.5m/s). d) Waveform fit of displacement recordings. On the synthetics we have superimposed those (grey) of the inversion when geodetic data were included to show that no significant change is observed to the seismic data fit. d) Slip rate snapshots of the model (Fig. S7a,c). All other notations as in Figure 4 and 5 of the main text.
Section 3. New Empirical Green's function (EGF) method to calculate Apparent

Source Time Functions (ASTFs)
The waveforms s(t) and S(t) of the weak event (EGF) and the mainshock, respectively, are defined by equations (1) and (2): Green's function g(t) is the same for both events and need not be known. The m(t) and M(t) are the moment rate functions. We assume a frequency range (detailed below) in which m(t) can be approximated as an isosceles triangle, centered at time t=0,whose duration is shorter than the duration of M(t). Function M(t) is expressed as a set of equidistantly shifted functions m(t); see equation (3) The ratio of the scalar moments of the mainshock and the EGF event (the relative moment) provides a constraint for the weights.
Generalizing for a three-component station (total number of time samples M), equation (4) with real data S, and equation (5)   ) .
In practice, the last row of the matrix and the the last value of the data column must be multiplied by a constant, i.e. const. w i = const. M o /m o . The constant does not alter equation (5), but guarantees its proper balance with eq. (4). The numerical value of the constant depends on relative values of S with respect to M o /m o , depending on whether data S are in counts or in meter per second. In this paper we use const. ~ 1.0e7.
Further assuming that w i ≥ 0 for each i, we solve the system (6) by non-negative leastsquares inversion (NNLS) after Lawson & Hanson (1974). The quality of the inversion is measured by the fit between the real mainshock record and the synthetic record (eq. (4)), quantified by variance reduction. Equation (3) then provides M(t), the desired nonnegative ASTF for a given station.
The procedure is performed on the s(t) and S(t) records equally filtered with a band-pass filter (Harris, 1990). Assuming frequency band (Fmin, Fmax), the duration of the triangle m(t) is defined as 1/Fmax. As such, the shortest temporal variation of the ASTF that can be resolved is 1/Fmax.
A Fortran code and Gnuplot graphics scripts have been developed to perform the inversion and to automatically visualize the results. In the case of quality data, the software provides an ASTF, which is: 1. Non-negative (by definition).
2. Causal, i.e. starting generally at origin time (t=0); for discussion of possible small signals before t=0, see below.
3. Stable, i.e. having generally only small artifacts beyond the major ASTF part. For details about the artIfacts, see below.
4. The area of ASTF is proportional to the relative moment (Mo/mo) at each station.
The program reads the three-component ASCII waveforms (time, NS, EW, Z) of the mainshock and EGF, recorded by the same instrument. No instrumental correction is needed. Before the first code run, both seismograms are aligned to have the same Pwave arrival times. Since the locations of the mainshock and EGF are not exactly identical, the P-wave alignment does not guarantee the S-wave alignment. Therefore, if inverting the whole record, or only S waves, we must allow for a possible start of the resulting ASTF before t=0, mentioned above as the small acausal effect. Regarding compactness of the ASTF, code is executed repeatedly, using either the whole set of the calculated weights w 1 ,… N , or just w J ,… K , where J≥1 and K≤N. If the fit between real and synthetic seismograms is similar for the <J, K> interval of the weights, the weights outside of this interval are considered to be noise (artifact).
The entire inversion process is controlled by a single configuration file. The user can set up several parameters: a) Using the whole seismogram, or define a (smoothly tapered) time window for inversion that contains, e.g., P or S waves only. b) Selection of the station components to be used in the inversion. c) Time interval T = <t 1 , t 2 >, where weights can be inverted. It must be greater than the largest expected ASTF duration. A small time interval before origin time (i.e. t 1 <0) is advisable. d) Time shifts  of the weights . The number of unknowns N in equation (6) is given by the time interval T and time shift , N = T/ e) Scalar seismic moments, 0 , 0 .
f) Option whether to invert the original or integrated records.
g) The parameters of the band-passfilter .
Graphical outputs: Two gnuplot scripts (seismo.gpl and rstf.gpl,) display the results. The examples of the output are at Figs S9 -S12.
The three-component normalized waveforms and amplitude spectra of the mainshock and EGF event are displayed by seismo.gplgnuplot script (Figs S9,S10). The numbers at the waveforms panel are the true amplitudes. The waveforms are filtered in the same frequency range as that used in the inversion. The frequency range is marked by green zone in the amplitude-spectra plot.
The output of rstf.gpl script (Figs S11, S12) contains the header which summarized the general parameters (station name, original or integrated input records, frequency band in Hz and seismic moment ratio) and shows the legend of the waveform panel.
The main result of inversion, the weights are shown in the top left panel in time interval T = <t1, t2>. Figure below is a Moment rate function constructed according to equation (3) using the elementary triangles of width =1/Fmax, which is shown in the legend together with the maximum moment rate value. As a next plot is the Cumulative sum of weights.
Synthetic seismograms at the right panel are calculated according to equation (4) for all calculated weights (red) and for their subset w JK ,(green), mentioned above. The numbers are the true amplitudes and variance reductions for each component. The user can compare how the subset of the weights fits the observed mainshock (blue) relative to the fit employing all weights. In this way, redundant weights are identified and the corresponding 'tail' of ASTF is removed as a noisy artefact.
Amplitude spectra of the NS components for EGF (black) event, observed (blue) and synthetic (red) mainshock are placed at the left bottom part of the figure. Figure S10. As in fig. S9 for an example strong motion station (ASTA) Figure S11. Summary plot of the EGF method for the example station ARG. See text for detailed explanation. Figure S12. As in Figure S11 for station KARY.