The 30 October 2020, MW = 7.0, Samos earthquake: aftershock relocation, slip model, Coulomb stress evolution and estimation of shaking

We study the major MW = 7.0, 30 October 2020, Samos earthquake and its aftershocks, by calculating improved locations using differential travel times and waveform cross-correlations. We image the rupture of the mainshock using local strong motion data, and we examine the Coulomb stress evolution prior to the mainshock, as well as the coseismic stress changes. Lastly, we estimate the produced shaking using all the available information from strong motion data and testimonies. Earthquake relocations reveal the activation of the E-W oriented Kaystrios fault, in the North basin of Samos with a possible extension to the West. The kinematic rupture inversion suggests non-uniform bilateral rupture on a ∼60 km × ∼20 km fault area, with the main rupture propagating towards the West and maximum slip up to approximately 2.5 m. Improved locations of the aftershock sequence are anti-correlated with areas of maximum slip on the fault surface. Similarly, the Coulomb stress change calculations show that only off-fault earthquake clusters are located within lobes of increasing positive static stress changes. This observation is consistent with assuming a fault area of either uniform slip, or variable slip according to the obtained slip model. Both scenarios indicate typical stress patterns for a normal fault with E-W orientation, with stress lobes of positive ∆CFF increments expanding in E-W orientation. In the case of the variable slip model, both negative and positive stress changes show slightly larger values compared to the uniform slip model. Finally, Modified Mercalli Intensities based on the fault model obtained in this study indicate maximum intensity (VII +) along the northern coast of Samos Ιsland. Spectral acceleration values at 0.3 s period also demonstrate the damaging situation at Izmir.


Introduction
On October 30, 2020, at 11:51 UTC, a major M W 7.0 earthquake occurred offshore the north coast of Samos island in the eastern Aegean Sea, in close proximity with Asia Minor coast, causing the death of two people in Vathi (Greece) and 115 people in Izmir (Turkey) due to severe building collapse. The earthquake caused heavy damages which resulted in 19 fatalities and more than 1030 injuries in Samos Island and Turkey (GEER, 2020). A variety of geological effects such as coastal uplifts, ground fractures, ground deformation, were reported after the earthquake and attracted multiple scientific working groups in situ. A moderate tsunami was generated that mainly impacted the northern coast of Samos and the SW coastline of Izmir Province, Turkey . In Samos, runup exceeded 1.8 m in the town of Karlovasi causing minor damages, whereas the low-elevation waterfront of Vathi was impacted by a series of waves with maximum overland flow depth reaching ∼ 1 m (Kalligeris et al. 2021). Along the Aegean coastline of Turkey, a maximum wave runup of 3.8 m was measured in Akarca and flow depth values as high as 1.4 m were recorded in the worst-hit Kaleici region of Sigacik (Dogan et al. 2021). Tsunami warning messages were issued within 11 min after the earthquake by all three Tsunami Service Providers operating in the Eastern Mediterranean under the NorthEastern Atlantic, Mediterranean and connected seas Tsunami Warning System (NEAMTWS) of IOC UNESCO, and were followed by tsunami-ongoing messages following the detection of the tsunami by several tide gauges installed in the Aegean Sea . A few hours after the mainshock the strongest M 5.3 aftershock followed in close epicentral distance from the mainshock (Fig. 1). According to the fault plane solutions reported by various agencies (Table S1 in the Supplementary material) both earthquakes signify normal faulting, which is in agreement with the seismotectonic regime (see also Sect. 2). Historical archives indicate that the broader area of Samos has been occasionally struck by destructive earthquakes, with the first record around 200 BC to be found on ancient inscriptions which describe Samos island suffering from damages due to a strong earthquake (Papazachos and Papazachou, 2003). There was an absence of information for hundreds of years until the 18th century.
Τhere is evidence for nine earthquakes with estimated magnitude M ≥ 6.0 until 1955 when the M=6.9 Samos destructive earthquake occurred (Fig. 1). Furthermore, the island has been affected by strong earthquakes in the northern Aegean and western Turkey which induced damages and losses at Samos and generated tsunamis (Altinok et al. 2005;Melis et al. 2020). In recent times, microseismicity in the area of Samos and Kusadasi has been investigated with the deployment of a temporary network by Tan (2013) where dense earthquake clustering was evident. The 2005, M W =5.8, moderate magnitude seismic sequence South of Karaburun peninsula has also been studied (Benetatos et al. 2006;Melis and Konstantinou 2006). The significance of earthquake interaction through stress transfer and the evolution of stresses due to strong events have been investigated for western Turkey by Paradisopoulou et al. (2010), for the northeastern Aegean Sea by Nalbant et al. (1998), and for the northern Aegean Sea by Papadimitriou and Sykes (2001) and Rhoades et al. (2010). Nevertheless, the Samos 2020 causative fault has not been considered in previous studies. In this study, we attempt to shed some light on the characteristics of the M W =7.0, Samos earthquake, its rupture history, as well as its aftershock sequence and its relation to the tectonics of the study area. Moreover, we investigate the mechanisms of stress transfer prior to the occurrence of the mainshock, and we assess the impact of Coulomb stress changes on the evolution of the aftershock activity.

Study area and its seismotectonic setting
The study area is situated in the eastern edge of the Aegean arc within the transition zone between the fast moving Aegean and the Anatolia microplate where deformation is transferred into the fast moving Aegean Sea, as deduced from GPS and seismological data (Papazachos, 1999). The broader area has been repeatedly struck by destructive historical earthquakes which are related to active seismogenic faults, built in a complex seismotectonic setting. The broader Aegean area undergoes widespread NNW-SSE extension orthogonal to the subduction of the eastern Mediterranean plate under the Aegean microplate (Papazachos and Comninakis, 1971). The westward extrusion of the Anatolia plate and the prolongation of the North Anatolian fault into the Aegean Sea, which started 5 Ma ago, further reinforced the existing extensional forces (McKenzie 1972). The Aegean microplate Fig. 1 a General map of the study area (enclosed by the red rectangle) within the Aegean regime with respect to the main seismotectonic features of the North Anatolian Trough (NAT) and the Hellenic Trough; b: Seismotectonic map with the approximate locations of the most significant historical earthquakes (M ≥ 6.5) since 1850 obtained from the historical catalogue of Papazachos and Papazachou (2003), along with their inferred fault plane solutions (details are given in Table 3). Plotted with the black star is the 2020 Samos mainshock along with the Global CMT best fitting double-couple fault plane solution, whereas the circles in red outline refer to the strongest earthquakes which occurred within four months of aftershock activity. Main active faults as reported in the GEM active fault database (Styron and Pagani 2020) are also depicted 1 3 accommodates a southwestward movement relative to the stable Eurasia with a velocity rate of 32-35 mm/yr (LePichon et al. 1995;McClusky et al. 2000). Dextral strike slip faulting is dominant in the Northern Aegean as revealed by tectonics and fault plane solutions (Taymaz et al. 1991;Kiratzi, 2003). The trend of the extensional axis has been gradually rotated from NE-SW to NNE-SSW allowing the formation of new structures and causing the older NW-SE and NE-SW trending faults to acquire a strike slip component (Kissel and Laj, 1988). Currently the extensional axis strikes in an almost N-S direction, according to geodetic measurements (Armijo et al. 1996) and fault plane solutions for strong earthquakes (Papazachos et al. 1998).
Back arc tectonics and Tertiary volcanism are the dominant characteristics of the Aegean and the coast of western Asia Minor, which have given rise to the formation of several neotectonic basins (Ikaria and Samos basins). In the Western Anatolian Extension Province, which is dominated by N-S extension, a significant number of elongated E-W grabens like Gediz, Kucuk Menderes and Buyuk Menderes have been developed (Sengör et al. 1984), along with offshore and onshore N-S to NE-SW steeply dipping oblique slip faults, especially in Kusadasi peninsula and Izmir gulf showing a transpressional character (Ocakŏglu et al. 2005). According to the same authors, E-W compression in this area causes the N-S trending reverse faults, NE-SW dextral and NW-SE left-lateral strike slip faults, like the Karaburun fault. Onshore seismotectonic research in Samos Island by Mountrakis et al. (2003) highlighted the existence of active normal faults which shape the northern and southern coasts of the horst-like tectonics of the island and bound the Quaternary basins. NNW-SSE basins were initially formed by low angle detachment zones, but NE-SW extension in the Miocene imposed high angle faulting. The successive rotation of the stress field from NE-SW to NNE-SSW resulted in the formation of new E-W normal high-angle faults, with the NW-SE and NE-SW ones being reactivated by acquiring a strike slip component. The active fault databases of GreDaSS (Caputo and Pavlides, 2013) and GEM (Styron and Pagani, 2020) present the Kaystrios normal fault dipping offshore to the north of Samos Island with a slip rate of 1.0 mm/yr (Pavlides et al. 2009). Regarding offshore seismotectonic investigation (Lykousis et al. (1995) suggest that strike slip deformation is active at the eastern part of the asymmetric Ikarian basin between Ikaria and Samos Islands as also proposed by Stiros et al. (2000). The bathymetry analysis by Nomikou et al. (2021) suggested the existence of an E-W normal fault bounding Samos basin from the south with an average dip 45º and total throw 650 m since early Pleistocene, whereas slopes get steeper to the western part of the island mostly related to the NE-SW Ikaria margin.

Data
For the study of the mainshock as well as its aftershock sequence, we combined parametric phase arrivals and waveform data. Details on the availability, temporal and space distribution of the data being used are provided below.

Parametric phase arrival data
We downloaded hypocentral parameters and all the available phase arrival times of P and S phases up to 250 km in epicentral distance for all the earthquakes within our study area with magnitude M ≥ 2.0, using the web services of the European-Mediterranean Seismological Centre (EMSC, https:// www. seism icpor tal. eu/ fdsn-wseve nt. html, database last accessed March 2021), as it combines data from different providers within hours after an earthquake has occurred. Details on the search parameters are providedshown in Table 1. This search yielded parametric data for 2122 earthquakes with 67,775 associated P and S phase arrivals. Figure 2a shows the spatial distribution of the seismicity and seismic stations being used in the current study. The associated phase arrivals are mainly manually picked from the permanent seismic stations of the Hellenic Unified Seismic Network (HL, HT, HA, HP, HC, National Observatory of Athens, 1997; Aristotle University of Thessaloniki Seismological Network, 1981;University of Athens, 2008;University of Patras, 2000;Technological Educational Institute of Crete, 2006), coordinated by the Institute of Geodynamics, National Observatory of Athens (NOA-IG), Greece, the Kandilli Observatory and Earthquake Research Institute, Bŏgaziçi University (KO, Kandilli Observatory and Earthquake

Fig. 2 a
Map showing the spatial distribution of the seismic stations used in the current study with respect to the mainshock (orange star) and its aftershock sequence (white circles). Red inverse triangles show stations which have been used in earthquake relative locations, whilst green inverse triangles represent stations equipped with accelerographs whose recordings have been used to determine the mainshock's slip model; b: comparison of 1D velocity models used for the earthquake relative locations (Akyol et al., 2006), and Earth's structure used for the determination of the mainshock's slip model, which is extracted from the CRUST 2.0, 3D velocity model (Bassin et al., 2000), with respect to the 2º x 2º cell, with reference to the mainshock's epicentre Research Institute, Bŏgaziçi University, 1971) and the Turkish National Seismic Network (TU, Disaster and Emergency Management Authority, 1990), operated by the Disaster and Emergency Management Authority (AFAD).

Waveform data
Along with parametric phase arrivals, we used three-component continuous recordings from the seismic networks mentioned in Sect. 3.1 in order to calculate differential travel times. Each waveform's start time was defined as the earthquake's origin time and the end time was set to 20 s past the theoretical S arrival time based on the ak135 velocity model (Kennett et al., 1995).

Macroseismic observations and input data for shaking estimation
The "Did You Feel It" (DYFI) testimonies reported to EMSC (https:// seismicportal.eu/ event-details.html?unid = 20201030_0000082), along with Peak Ground Acceleration (PGA) values obtained from strong motion data (see also, Askan et al., 2021) were used for the estimation of Modified Mercalli Intensity (MMI) maps. This data is being used at NOA in a routine manner, under SeisComP3 monitoring software (Helmholtz-Centre Potsdam-GFZ German Research Centre For Geosciences and GEMPA GmbH, 2008). A module (scwfparam, https:// www. seisc omp. de/ doc/ apps/ scwfp aram. html? highl ight = scwfparam) is in operation to measure peak ground acceleration (PGA), peak ground velocity (PGV) and the pseudo absolute acceleration elastic response spectrum (PSA) at periods 0.3 s, 1.0 s and 3.0 s (see Sect. 4.4 for more details).

Earthquakes relocation
With the majority of seismicity located offshore at the north of Samos, only a few stations are located in close proximity to the epicentres (< 10 km), which is key to accurate earthquake location (e.g., Bondar and McLaughlin, 2009). Moreover, the fact that most of the aftershocks following the mainshock are of small magnitude (M < 3.0), large azimuthal and secondary azimuthal gaps are often closely associated with the picking of very few phase arrival observations (see Fig. 3 and Figure S1 in the Supplementary material).
In order to improve the locations of the existing catalogue (see Sect. 3.1), we applied a double-difference location technique (Waldhauser and Ellsworth, 2000) by calculating differential travel times (tt) obtained from both catalogue data and waveform cross-correlations. Considering two neighbouring seismic events i and j with hypocentral parameters i i h i t i T and j j h j t j T , respectively, the double-difference problem for any n phase observation relative to the two events i and j can be defined as: where Δm ij = Δd ij Δd ij Δdh ij Δdt ij T represents perturbations in the model space (m) which is defined by the relative hypocentral parameters between the two events i and j, and the second part of Eq. 1 represents the differential travel time residuals. Equation 1 can be used to form a system of linear equations for each station pair and can be solved by means of least squares in an iterative way. By linking as many neighbouring events together as possible (typically thousands) within small distances of a few kilometers, high resolution relative hypocentre locations can be achieved. In our case, we set the maximum separation distance to 10 km and the minimum number of links per pair to six, when at least four phase pairs are available. This setup yielded 18,096 event pairs with 153,301 P phase and 55,303 S phase differential travel time pairs. We then calculated cross-correlation differential travel times based on 694,001 seismograms, after removing the mean and applying a zero phase bandpass filter from 1 to 10 Hz to each waveform. Only phase pairs with a correlation coefficient above Using the 1D velocity model of Akyol et al. (2006), and adopting a V P /V S ratio of 1.73 ( Fig. 2b), both phase picks and cross-correlation differential times were combined in a dynamically weighted double-difference inversion giving more weight to catalogue phase data in the first stages of the inversion, whereas control is passed to cross-correlation differential times in the last stages, allowing the cross-correlation data to constrain only event pairs with separations smaller than 5 km.

Kinematic slip model
In order to calculate a kinematic rupture model for the mainshock, we applied the technique developed by Gallovič et al. (2015). The displacement wavefield u in space (r) and time (t) is described by the representation theorem: where G is the Green's function which contains the responses of point sources from subfaults distributed along the fault surface S, as described by the model. The term Δu( , ) represents the slip rate function in space ( ) and time ( ) , which is parameterised by overlapping Dirac functions distributed along the fault surface. This type of formulation imposes no constraints on the nucleation point, the rupture velocity, or the shape of the slip rate function. As a result, the inverse problem consists of a large number of model parameters, namely, samples of the slip velocity as a function of spatial coordinates and time which are linearly related to the wavefield. The inverse problem is then solved by applying smoothing and a non-negativity constraint on the slip rates as regularisation.
In the current study, we used three-component strong motion data (Fig. 2a, see also Sect. 3.2), removed the mean, filtered the accelerograms between 0.05 Hz and 0.5 Hz and converted it to displacement. Moreover, for the station SAMA in Samos Island, where the accelerograph is installed in a public building and oriented in parallel to the walls of the building, we rotated the two horizontal recordings to N-S and E-W directions prior to the processing mentioned above. Finally, based on the results from the relative locations (see Sect. 5.1) we set the fault surface measuring 60 km along strike and 20 km along dip and we calculated the Green's functions based on the Global Centroid Moment Tensor model (GCMT, https:// www. globa lcmt. org/, database last accessed March 2021, Dziewonski et al., 1981;Ekström et al., 2012) using the Earth's structure extracted from the CRUST2.0 velocity model (Bassin et al., 2000) with reference to the location of the main earthquake ( Fig. 2b), as this velocity model was built by combining both travel times and waveform data (surface wave dispersion measurements and normal modes).

Coulomb stress transfer
Changes in the stress field arise due to the coseismic stress changes induced by the occurrence of strong earthquakes along with the interseismic long-term stress accumulation which is accommodated on major faults and is driven from relative plate motions. In the case of Samos sequence, the stress state was investigated prior to the mainshock by investigating the successive stress changes imposed by the occurrence of strong earthquakes in the surrounding area (M ≥ 6.5) and the changes during the occurrence of the M W =7.0 mainshock. The methodology that was followed relies on the elastic rebound theory, according to which the stress released in an area existed prior to the event and the applied technique follows Deng and Sykes (1997). For the calculations of the interseismic strain accumulation, the "virtual dislocation" concept was introduced according to which the released coseismic stress pre-exists in the brittle part of the crust and is determined by assuming backward fault slip. The Coulomb Failure Function criterion examines the conditions under which failure occurs on rocks when shear stress exceeds rock strength (Scholz, 2002). Changes in Coulomb Failure Function (∆CFF) depend on changes in shear stress, ∆τ, and normal stress, ∆σ, resolved onto the earthquake fault plane according to: where µ′ is the apparent coefficient of friction. Positive ∆CFF values denote a high likelihood for future failure, therefore locations with advanced stress changes indicate areas close to rupture, whereas negative values indicate that fault failure is prevented. Subsequent earthquakes preferentially occur on locations with positive increments whereas negative values are considered areas of seismic quiescence described as shadow zones. The term µ′ describes the effect of the pore pressure change due to pore fluid and for dry materials ranges between 0.5 and 0.8 (Harris, 1998). For calculating ∆CFF, source models for large earthquakes are constrained to approximate the rupture geometry. Despite their heterogeneity, faults can be simply approximated as planar rectangular geometric structures, which dip into the brittle part of the crust. The geometrical parameters, such as the fault length and width, the co-seismic horizontal (u SS ) and along dip (u DS ) slip and the fault plane solution (strike, dip and rake) adequately describe the rectangular rupture models and are used as input for stress change evaluation. The selection of these parameters is crucial for the definition of the stress field since the variation of these parameters strongly affects the final stress pattern. The geometrical parameters are usually deduced using the local seismotectonic information. If this information is not available absent, empirical equations are applied: i.e. Wells and Coppersmith (1994) (2004) fits the average slip deduced from the slip analysis on Samos earthquake (1.64 m, see also Sect. 5.2), while fault lengths show no significant variations for different calculations (~ 2 km). Therefore, the set of equations given by Papazachos et al. (2004), regarding length, L (Eq. 4) and average coseismic displacement, u (Eq. 5) for dip slip continental faults was applied in place of the missing seismotectonic information: In cases of instrumental recordings where the seismic moment (M o ) is known, the average coseismic displacement, was directly calculated from Eq. 6: where G is the shear modulus in the seismic source, estimated approximately equal to 3.3 × 10 11 dyn × cm −2 (Stein et al., 1997), u is the average displacement and S corresponds to the fault area. Fault width, w, was obtained from the equation w = h∕sin( ) , where δ corresponds to the fault dip and h is the downdip distance from the upper to the lower edge of the fault, always taking into consideration the value of the aspect ratio (fault length over width, L/w). The width of the seismogenic part of the crust for which ∆CFF were calculated, was determined from the distribution of the relocated earthquake foci and ranges between 3 and 14 km.

Shakemap calculations
ShakeMap calculations were derived by combining PGA, PGV, PSA (at 0.3 s, 1.0 s and 3.0 s) values computed using the swfparam module in SeisComP3 and the DYFI testimonies collected by EMSC. We attempted to estimate and map the shaking of the M W 7.0, Samos earthquake, using the United States Geological Survey (USGS) ShakeMap4 (Worden, 2016) standard procedure. PGA values were converted to MMI following Worden et al. (2012), which satisfies shallow crustal events in Greece . Site effects were taken into account using the V s 30 gridded layer produced by Stewart et al. (2014) and made available at the USGS ShakeMap repository and github (https:// usgs. maps. arcgis. com/ apps/ webap pview er/ index. html? id= 8ac19 bc334 f747e 48655-0f328 37578 e1, https:// github. com/ usgs/ earth quake-global_ vs30/ tree/ master/ Greece), combined with the Ground Motion Model (GMM) of Boore et al. (2021). The location and dimensions of the fault plane were based on the slip model obtained in the current study using the technique presented in Sect. 4.2, which corresponds to the best fitting double-couple model from GCMT.

Relative locations
Relative locations were obtained following the procedures described in Sect. 4.1 after experimenting with different setup parameters and weighting schemes. In fact, in some cases high damping was needed to stabilise the inversion, as expressed by the ratio of the largest to smallest eigenvalue of the system (condition number). This might be an indication of weak links between events and/or the presence of data outliers, possibly due to the fact that the majority of the seismic events in our set are located offshore where the density of the seismic stations in the close proximity of the epicentres is not ideal. A way to overcome this requires more links between event pairs in order to form continuous clusters, changes in the weighting of the catalogue and cross-correlating differential travel time data, or even the generation of differential time data sets that allows for more neighbours for each event by taking into account more distant events. Nevertheless, this can be a balancing act between the resolution of the relocations and the number of relocated events.  and those that were relocated using differential travel times from catalogue and waveform cross-correlation data (orange circles); b: map showing the spatial distribution of earthquake density before the relocation on a 5 km × 5 km grid; and c: same as in b but using the relocated epicentres We achieved acceptable condition numbers and we obtained relative locations for 1357 out of 2122 seismic events by combining catalogue phase picks and cross-correlation differential travel times (Fig. 5). Sparse earthquakes initially located off the main cluster which lies offshore the north shore of Samos island were rejected (Fig. 6a) as these events did not meet the requirements set regarding the minimum links per pair and/or maximum separation (see also Sect. 4.1). Moreover, event pairs associated with similar wave trains showing high correlation coefficients (> 0.7) yielded a higher density of earthquakes located in a narrower zone along the Kaystrios fault in the Samos basin ( Fig. 6b and c). Figure S2 in the Supplementary material shows the evolution of the relocation process every time a different weighting scheme was applied in the inversion, carrying out 25 iterations in total. When the catalogue data and cross-correlation differential travel times were equally weighted, distinct earthquake clusters formed as early as in the 10 th iteration. Essentially, the inversion converged approximately at the 20 th iteration when catalogue phase data are down-weighted and cross-correlation differential travel times almost entirely control the inversion. Figure S3 in the Supplementary material provides information regarding the quality of the event pairs constrained by cross-correlation data and how it relates to event separation. In general, travel time residuals are centred around zero seconds following a normal distribution and the RMS residuals have the tendency to increase with increasing event offset, as expected. This increase with offset distance may be explained as being a result of scattering along the ray path and/or discrepancies in the source mechanisms for each event pair.
The onset of the Samos earthquake sequence was marked with the occurrence of the M W =7.0 mainshock, which was followed by the strongest aftershock M=5.2 a few hours later and in very close distance. The stem diagram of Fig. 7a shows the temporal evolution of the seismic sequence for the initial earthquake catalogue reported from EMSC for earthquakes with magnitude M ≥ 2.0. Dense earthquake occurrence was observed within the first 15 days of November 2020, while by the end of the study period only 10 moderate events with magnitude M ≥ 4.5 occurred. For a clear investigation of space-time earthquake relations, distinct colours are assigned to the most significant spatial clusters (Figs. 7b and c). Seismic activity expanded equally to the West and to the East of the main shock, along an elongated zone offshore, North of Samos Island. Synchronous seismic activity appeared to the West by forming two clusters. The most numerous cluster (in red) shows a quite diffused distribution of epicentres with an increased activity within the first 10 days, which coincides with the activity and duration of the smaller blue cluster. Some activity of small magnitude earthquakes bursts to the East simultaneously, but it is not E-W aligned. Some hours after the mainshock, seismic activity migrated to the eastern coast of the island (in green) and two seismic bursts also appeared to the northeast (yellow and magenta), in the proximity of the Turkish coast. Seismic rate is considerably decreased all over the area, approximately 15 days after the mainshock.
In order to investigate the seismic faults activated by the mainshock and its aftershocks, we examine the cross-sections of Fig. 5. Cross-sections A and B are oriented along strike to the Kaystrios fault, whilst the rest are oriented perpendicular to it. All cross-sections are 5 km in width and those that are parallel slightly overlap to each other.
Along strike (E-W) cross-sections (Fig. 8) cover almost the entire seismicity in the study area and show earthquake clusters at a total length of 60 km approximately, with the mainshock's epicentre located 40 km from the West (A1) and 20 km from the East (A2). The vast majority of the relocated earthquakes are observed at depths ranging from 3 to 15 km. Four distinct clusters are formed along this direction, with the largest observed at close proximity to the mainshock and the largest aftershock (M=5.2) hypocentres. Cross-section B, which is located slightly to the South, offers a clearer view to these distinct earthquake clusters, suggesting the activation of possible parallel seismic faults to the main Kaystrios fault. The cluster located to the western part covered by the cross-sections A and B, could be considered as a result of a possible extension of the Kaystrios fault to the West, as suggested by fault mechanisms ( Figure S4 in the Supplementary material). North-South oriented cross-sections (Fig. 8) indicate that the main cluster observed at the Samos basin is associated with the Kaystrios fault, which is rather shallow dipping for a normal fault (∼40°), in agreement with the earthquake fault plane solutions in the area. Some minor activity may be associated with the presence of other smaller parallel seismic faults, suggested by the similarity of source mechanisms ( Figure S4 in the Supplementary material).
In general, sharp images of relocated seismicity, especially as shown in C and D crosssections, revealed a possible system of listric faults in the Samos basin, dipping to the North. The earthquake cluster observed at the East coast of Samos is not associated with the North-dipping Kaystrios fault and may be explained by the activation of other seismic faults in the area with steep dipping angles (see cross-sections B and F in Fig. 8). Since the seismicity of this cluster is of rather low magnitude (∼3.0), there is a lack of source mechanisms which could reveal the characteristics of their associated sources, and hence, it is not safe to draw any conclusions at this point regarding this earthquake cluster prior to further investigation.

Slip model
Using the relocated hypocentral solution obtained for the mainshock in Sect. 5.1, we carried out a kinematic slip inversion described in Sect. 4.2 in order to determine its slip model. Based on the distribution of the relocated aftershock sequence, we assumed a planar fault for simplicity, with a fault rupture area of 60 km × 20 km. The nucleation point was placed 20 km from the West boundary along strike, and 8 km in the up-dip direction (see also cross-section A at Fig. 8). Several source models determined for this earthquake (Table S1 in the Supplementary material) show a rotation angle (Kagan, 1991) up to ∼30º which is typical among source models obtained using different data and techniques (e.g., Lentas et al. 2019). Despite these variations within source parameters based on the nucleation point or centroid based solutions, we used the best fitting double couple solution from the Global CMT (Dziewonski et al. 1981;Ekström et al. 2012) as an average representation of the source. We used the seismic moment determined by the GCMT model and we set the source duration twice the GCMT half-duration (7.6 s), as being a good approximation of the total rupture time (Duputel et al. 2013;Lentas et al. 2013). The GCMT solution assumes a triangular source time function with half duration determined by a constant stress drop scaling relation, proportional to the seismic moment. Table 2 summarises the input parameters used in the kinematic slip inversion.
Unlike synthetic tests, real case applications usually suffer from non-uniform station azimuthal coverage or data quality issues. Our case is no exception regarding the station coverage due to the topography of the study area, where stations to the west are sparse and far from the epicentre (green triangles in Fig. 2a). Hence, we tried to overcome this issue by slightly Fig. 7 a Temporal distribution of the reported seismicity for the study area, used as input to the relocation analysis from 30-10-2020 to 01-03-2021 along with the cumulative number of earthquakes in red colour; b: Map view of the relocated seismicity; and c: Spatial-temporal diagram of the relocated seismicity projected in E-W orientation. Different colouring is used to highlight distinctive earthquake clusters both in the map view and the space-time plot 1 3  Figure 9 compares the data against synthetic waveforms with respect to the calculated slip model. Synthetic displacement seismograms show a very good fit against observed data in most cases, especially for the stations closest to the epicentre (up to 100 km). Nevertheless, the station located at Tinos island (TNSA) may suffer from timing errors, thus, it was strongly down-weighted. Some discrepancies between the amplitude fit of waveform data and CRUST 2.0 synthetics may be improved by the use of a more accurate velocity model for the area. Figure 10 shows the rupture evolution and the composite slip model obtained for the mainshock. Within the first three seconds the rupture propagated asymmetrically in the up-dip direction, mainly towards the East, whereas, in the next four seconds (4-8 s) it showed signs of near-simultaneous failure of two asperities both up-dip and downdip. Since our kinematic inversion is based on a single source, this could potentially be an artefact and less likely an indication of two sub-events. Even though multiple sub-events are not very common, nevertheless, there are cases of strong earthquakes (M W ≥ 6.5), where multiple sub-events have been identified. For example, Sokos et al. (2016) carried out kinematic slip inversions based on multiple point source modelling for the 2015, M W = 6.5, Lefkada earthquake, and showed that the rupture must have involved at least two sub-events with a time gap of approximately 4 s. Their findings were further supported by independent studies based on geodetic data (Bie et al. 2017).
Next, the rupture is characterised by westward propagation, with the maximum slip observed up-dip (Fig. 10a). Individual slip rate functions (Fig. 10b) indicate a more prominent slip patch in the up-dip direction, whereas down-dip propagation was slightly shorter in time. Slip towards the East is almost negligible, whilst, slip in the West fades out at the 14th second, with just a very short episode of rupture (pulse-like) in the last second. However, this might be just an artefact, since no matter how we set the total rupture time in the inversion, the ending is always ambiguous, possibly due to the fact that the station coverage to the West is sparse. Nonetheless, based on the assumed source duration (15 s) and fault dimensions, the obtained seismic moment (4.01 × 10 19 Nm) from our slip model agrees well with that of the GCMT.

Stress state prior to Samos 2020 earthquake
Stress accumulation on the causative Samos fault before the 2020 earthquake was investigated by incorporating the interseismic deformation along the fault according to its slip rate and the coseismic ∆CFF due to the occurrence of the known historical earthquakes in the vicinity of the study area. Uncertainties related to estimated magnitude and location of historical earthquakes along with lack of data alter the real stress state prior to the earthquake. For the purposes of the study, strong earthquakes reported after 1881 when the devastating Chios-Cesme earthquake occurred are considered reliable to be involved in the stress field reconstruction. Six earthquakes with magnitude M ≥ 6.5 have struck the study area from 1881 until present, most of which regard normal fault plane solutions, shown in Fig. 1 and further described in Papazachos and Papazachou (2003). Information on the parameters used for approximating the source models and the determination of their coseismic   (2003) stress pattern is given in Table 3. Coseismic stress changes caused by the occurrence of each earthquake were computed according to the source properties described in Table 3 but resolved onto the dipping plane of Samos fault according to the GCMT solution (see Table S1 in the Supplementary material), which is the receiver fault (planar calculations at 8 km depth are shown in Figure S5 in the Supplementary material). Shear modulus and Poisson's ratio are fixed at 3.3 × 10 4 MPa and 0.25, respectively. Figure S6 in the Supplementary material subplots exhibit the successive evolutionary state results for the receiver fault. Results incorporate the cumulative effect due to the tectonic loading according to the slip rate and the progressive coseismic changes. Before 1881, stress was presumed to be zero. Figure  The cumulative effect of the coseismic stress changes for the given rupture models along with the estimated aseismic deformation shows a progressive stress built-up along the Samos causative fault, which indicates a future failure promotion. A planar view at 8 km depth presents the stress state before the occurrence of the 2020 earthquake according to which, bright zones are formed along the fault area north of Samos, as well as to the east and west of Samos Island ( Figure S7 in the Supplementary material) mainly due to the cumulative contributions of the 1881, the 1902 and the 1955 strong earthquakes.

Samos 2020 earthquake coseismic coulomb stress changes
Coulomb stress changes due to the Samos mainshock were calculated with the use of the uniform and finite fault models determined in this study. The GCMT solution which signifies a pure normal fault (-90º rake) striking 276º with 34º dip was employed for constructing the rupture model. In the first case, the source fault zone was approximated with a rectangular plane 35 km long and 20 km wide, which coincides with the area of maximum slip and the length of the aftershock zone, and is slightly smaller than the fault rupture length estimated using the empirical scaling relation (Eq. 4, L = 43 km). Average slip was defined according to Eq. 6 for M o = 4.01 × 10 19 Nm determined by the slip inversion analysis with an along dip component (u DS ) equal to 1.64 m (u SS = 0). ΔCFF was calculated for three different horizontal layers of the seismogenic zone at 5 km (Fig. 11a), 8 km (Fig. 11c) and 10 km depth (Fig. 11e) where the majority of hypocentres are found. The stress pattern is typical for a normal fault shedding a broad stress shadow in a N-S direction where potential rupture is prohibited for a similar faulting type. Stress lobes with positive ∆CFF increments expand in an E-W trend, and enhance stress changes over the Ikaria-Samos basin and the coasts of Turkey and Menderes basin to the West. Most of the aftershocks, including the strongest aftershock (Fig. 11c) lie along the main dislocation plane and are concentrated in the central and the western parts of the shadow zone. Visual inspection ascertains a good correlation between seismicity and spatial distribution of the positive stress values only for the off-fault seismicity clusters. The western cluster is entirely located within the bright zones, whereas the cluster to the southeast is not fully explained by the distribution of stress resolved onto the north dipping plane and is probably attributed to the activation of a secondary fault. The cluster is effectively explained by the Coulomb stress evolutionary model, with the majority of the recent seismicity located where increasing positive static stress changes are calculated.
In the case of the finite fault model, a 60 km long and 20 km wide and 1 km × 1 km gridded fault plane was introduced in order to investigate the determined slip patches across the entire rupture zone. Stress was calculated for the depth of 5 km, 8 km and 10 km shown in Figures 11b, d and f along with the seismicity for the corresponding depths. The shadow zone is expanded as expected, with patches of positive stress changes observed to the east of the aftershocks at 8 km and at 10 km depth. Both negative and positive stress changes experience larger values compared to the uniform slip model. For investigating the correlation of the earthquakes with the given rupture model, six cross-sections were plotted following the cross-section definitions provided in Fig. 5, where ∆CFF changes are projected onto vertical planes (Fig. 12). Earthquake hypocentres in the 10 km range are additionally plotted. A1-A2 and B1-B2 profiles run the entire zone in an ENE-WSW direction. The southern profile shows that the hypocentres are embedded into the shadow zone, whereas the two zone terminations, especially the eastern part where seismicity is dense, exhibit positive stress increments. To the North of the zone (B1-B2), the shadow zone is thinner but also coincides with the hypocentral distribution. N-S sections indicate a positive correlation between hypocentres and stress distribution to the eastern fault edge where increased ∆CFF is found especially for seismicity between 8 and 10 km depth, as shown in sections D1-D2 and E1-E2. Section F1-F2 encompasses the southeastern cluster which seems not to be favoured by a stress pattern resolved on the north-dipping dip slip fault.

Estimation of shaking
Modified Mercalli Intensities (MMI) with emphasis given to the epicentral area are shown in Fig. 13, where the maximum expected intensity (VII +) is observed in the North of the Samos island, following the location of the rupture area. In fact, the NW part of the island suffered greater damage, as many rural villages were devastated, where old and poor quality stone masonry structures did not withstand the shaking. In general, this coincides with the results of the expected intensity calculations. Moreover, we note high values of ∼25% g PSA at 0.3 s period (spectral acceleration at 0.3 s) in the city of Izmir, which are expected to affect the high rise buildings in the greater vicinity. The latter combined with the amplification of the ground motion due to the basin soil conditions could explain the collapse of the buildings in this area, where the greater number of fatalities was observed. The above trend is also evident in static maps of macroseismic intensity (MMI), peak ground acceleration (PGA), peak ground velocity (PGV) and pseudo-spectral acceleration (PSA) at 0.3 s, 1.0 s and 3.0 s periods, which are shown in Figures S8-S13 in the Supplementary material. Diagrams of the fitted Ground Motion model  also demonstrate excellent fit to PGA, PGV and PSA from strong motion data, whilst converted values from testimonies are somewhat explained by the model (see Figures S14-S16 in the Supplementary material).

Discussion and conclusions
In the current study we analysed the major M W 7.0, October 30, 2020, Samos earthquake, and its aftershock sequence. We relocated the mainshock relative to its aftershocks, using differential travel times from phase catalogue and waveform cross-correlation data (Waldhauser and Ellsworth, 2000), covering a broad time period of four months. We carried out a kinematic slip inversion  for the mainshock and we examined Coulomb stress transfer (Deng and Sykes, 1997) as a result of the main earthquake, using a uniform slip model based on a simple planar fault, as well as using the fault slip distribution that we obtained from our kinematic slip inversion. Finally, based on the rupture determined from our slip model, we estimated the shaking in the epicentral area by combining strong motion data and DYFI EMSC-testimonies.
Improved relative locations of the mainshock and aftershocks were obtained by combining catalogue and waveform cross-correlation data based on the hypocentres downloaded from the EMSC database, and the use of a 1D velocity model (Akyol et al., 2006) for the broad area of the East Aegean Sea. The EMSC is expected to offer an acceptable basis of initial hypocentral solutions, as it combines phase arrival data from various sources, which are being used in the relocated EMSC hypocentres. Other studies also suggest the determination of absolute locations prior to the relocation using differential travel times and the calculation of a more appropriate 1D velocity model for the study area when working with local stations (i.e., Matrullo et al., 2013;Konstantinou, 2018;Konstantinou et al., 2020). This also offers the advantage of calculating station corrections which reflect unmodelled 3D structure and can lead to lower station time residuals. This is evident as phase arrival data in local epicentral distances is very sensitive to lateral heterogeneities as the rays are mostly up-going. Nevertheless, synthetic tests based on velocity models with substantial differences, have shown that relative locations based on linearised double-difference equations can still converge close to the true locations (Waldhauser and Ellsworth, 2000). This is also enhanced by the use of cross-correlation differential travel times which are very Fig. 13 Estimated intensity MMI contours using the resulting rupture model, strong motion maximum PGA values at recorded stations converted to MMI and EMSC collected testimonies. Filled triangles correspond to recorded stations and filled circles to EMSC testimonies. Colours are following USGS MMI ShakeMap legend effective in relocating repeating earthquakes. In fact, hypocentral separation is one of the key factors in accurate relocations of seismic events using double-difference and cross-correlation methods (Waldhauser and Schaff, 2008).
The relative locations of the earthquakes in our dataset revealed sharper images of the main tectonic characteristics of the fault area, which extends 60 km along strike (∼270º), and approximately 20 km along dip. Fault plane solutions of the mainshock and aftershocks are in excellent agreement with the orientation of the Kaystrios fault in the Samos basin. Moreover, the available source models suggest dip angles that range from 30 to 55º and indicate a possible system of parallel listric faults in Samos basin. This is also in agreement with morphotectonic studies (i.e., Nomikou et al., 2021) which support the existence of an E-W normal fault in the Samos basin with an average dip of 45º.
Using the improved hypocentral solution of the mainshock obtained from the doubledifference relocation, and the best fitting double-couple solution from the Global CMT, we parameterised the fault area defined by the relative locations of the aftershock sequence, accordingly, and we attempted to calculate a kinematic slip model for the mainshock using displacement data from local strong motion stations.
The model space of linear slip inversion problems is characterised by a large number of model parameters, which typically leads to non-unique solutions of finite-fault source models (Zhang et al. 2014). As such, slip inversions are ill-posed problems that require some kind of regularisation in order to stabilise the inversion. The latter is usually implemented by applying nonlinear constraints, such as positivity and/or spatiοtemporal smoothing in order to damp any arising artefacts (Gallovič and Zahradník, 2011). Since spatial smoothing is strongly affected by station coverage and weighting, it is important to be implemented in a way that equalises the effect of nearby and distant stations (Gallovič and Zahradník, 2011;Gallovič et al. 2015). For example, nearest stations which typically have larger amplitudes compared to distant stations, can introduce biases to the obtained slip model, in a similar manner that non-uniform station azimuthal coverage can lead to false unilateral rupture propagation towards the azimuth with the highest station density. Moreover, data noise and/or inappropriate velocity models can further give rise to spurious effects and small-scale heterogeneities in the obtained slip model.
After careful weighing of the waveform data in our kinematic slip inversion, we imaged the bilateral rupture of the mainshock within the first four seconds, as well as the main westward propagation towards the ending of the rupture, suggesting two main episodes of rupture which are also evident in the obtained moment rate function. Patches of maximum slip in our model are anticorrelated with the spatial distribution of the relocated aftershocks as expected (e.g., Konca et al., 2007;Kim and Dreger, 2008;Sladen et al., 2010), and the main cluster of the aftershock sequence around the mainshock, seems to lie on the fault area with the minimum slip (Fig. 14). Only a small cluster to the West is associated with moderate slip, which might be an indication of a false slip patch due to unmodelled 3D structure in the inversion, non-uniform station azimuthal coverage, or a more complex rupture that cannot be fully explained by the assumption of a planar seismic fault. Joint inversions of seismic and geodetic data could possibly provide a more robust slip model (i.e., Sladen et al., 2010).
Coulomb stress analysis was carried out in order to investigate the stress state on the fault before the 2020 Samos earthquake as well as the coseismic stress disturbances and its impact on the spatial evolution of the aftershock activity. The evolution of the stress field included the cumulative impact of coseismic stress changes of six strong earthquakes with magnitude M ≥ 6.5 beginning with the 1881 Chios-Cesme earthquake and the aseismic deformation on the fault. We assumed that loading on Kaystrios The GCMT best fitting double couple solution used for the determination of the slip model is also shown on the map. The relative locations of the aftershock sequence are shown as black dots; b: moment rate function determined from the kinematic slip inversion for the mainshock; c: same as in a, but in 3D plot. The surface trace of the planar fault assumed in the slip model inversion is shown in green for reference fault which was simply approached as a rectangular fault source, 36 km long and 20 km wide, is accommodated by a moderate rate equal to 1 mm/yr obtained from seismotectonic information (Pavlides et al., 2009). Before the 1881 earthquake occurred, stress on the fault was zero, which means that all the energy along the fault is released and starts building up from the beginning, gradually affected by nearby strong earthquakes. In all cases the stress field is resolved onto the properties of the Samos receiver fault. The fault plane solutions regard pure normal and oblique normal faults dipping to the North and to the South. The shadow zone that was cast after the 1883 event was recovered by the 1904 earthquake at South Samos which increased the ∆CFF values for more than 1 bar. After this point, positive stress changes are reinforced on the fault, moving the fault closer to failure even if the earthquakes responsible for the corresponding stress changes have occurred on antithetic south-dipping faults like the 1928 or the 1955 earthquakes (see also Table 3). After the last strong earthquake in 1955 the entire fault zone is characterised by positive stress increments with patches of stress with more than 2 bars as shown at the cross-section and the surface fault projection.
Static stress changes influence the location and timing of subsequent strong earthquakes or aftershock activity (Scholz, 2002). Positive correlations between Coulomb stress changes induced by a mainshock and the distribution of aftershock hypocentres has been generally affirmed (King et al., 1994;Parsons, 2002). In our case, the investigation between the main shock coseismic stress changes and the locations of the aftershock activity was performed by modelling coseismic stress changes both for a uniform slip 36 km long fault and a variable slip 60 km long fault zone. In both cases the vast majority of the aftershock locations is spatially distributed over the elongated shadow zone with high negative stress changes which is shed along fault strike and denotes the rupture area. Due to the pattern of the stress field, which is typical for a pure normal fault, the clusters which evolved at the terminations of the faults are located within the stress bright lobes. Especially the off fault seismicity at the western edge of the Kaystrios fault (red colour) is totally embedded within the bright zone and favoured by stress enhancement, as also the calculation of ∆CFF at the foci of the earthquakes show. The stress field analysis based on the variable slip model is significant because it sufficiently explains the existence of the majority of aftershocks at the eastern part of the rupture zone, since it coincides with the existence of a patch with positive stress increments, which probably have triggered microseismicity. The cluster in green which is located at the southern part of the seismogenic fault cannot be fully explained by the stress pattern, probably because its activation is attributed to a fault with different fault properties. Calculations for a stress field resolved for a N-S fault with dextral strike slip component as seismotectonic investigation indicates for the broader area, explain more clearly the existence of the activity.
Our analysis of the major M W 7.0, Samos earthquake, and its aftershock sequence, revealed the activation of the E-W oriented, Kaystrios normal fault in the North basin of the Samos Island. The obtained slip model, as well as the Coulomb stress changes due to the mainshock, are in agreement with the improved relative locations of the aftershocks, suggesting that the main rupture propagated mostly to the West. The estimated shaking in the epicentral area explained the observed severe damage to the NW part of Samos, as well as the collapse of buildings in the city of Izmir, further enhanced by the local soil conditions.