Spatiotemporal characteristics and earthquake statistics of the 2020 and 2022 adjacent earthquake sequences in North Aegean Sea (Greece)

ThetwomoderateearthquakesthatoccurredcloseandtothenorthoftheNorthAegeanTrough (NAT) on 26 September 2020 (M w 5.3) and 16 January 2022 (M w 5.4), both followed by aftershock activity, are examined. Seismic activity along the NAT and its parallel branches is continuous and remarkable, with numerous strong instrumental (M≥6.0) earthquakes. Yet, the frequency of moderate (5.0≤M<6.0) earthquakes outside these major fault branches is rather rare and therefore their investigation provides the optimal means to decipher the seismotectonic properties of the broader area. The temporal and spatial proximity of the two seismic excitations from late September of 2020 through early 2022, intrigues for exhaustive investigation of seismic activity with the employment of earthquake relocation techniques, moment tensor solutions and statistical analysis. Our research revealed that this seismic activity purely falls inside the Mainshock – Aftershock type, with fast aftershock decay rates and moderate productivity. According to our findings, the two seismic sequences, despitetheircloseproximity, exhibitdistinctivefeaturesasaresultoftheintricatestressfieldgen-erated at the western termination of the NAF system in an extensional domain


Non-technical summary
On 26 September 2020 and 16 January 2022 two moderate earthquakes (Mw5.3 and Mw5.4,respectively) occurred at the North Aegean Sea, southern of Chalkidiki peninsula.Their close proximity in space and time and the rare manifestation of such moderate events in the area promotes their analysis in order to better understand the faults and the state of stress in the broader area.We used the available seismological stations in the area to enhance the quality of the earthquake catalog and thoroughly investigated the properties of the two seismic sequences.We found that this activity is not directly related to the prevailing seismotectonic feature of the area, namely the North Aegean Trough, but it is derived from secondary features, typically found in the vicinity of such large complex systems.

Introduction
The rollback of the oceanic lithospheric plate of eastern Mediterranean that subducts beneath the continental crust of the Aegean microplate (Pichon and Angelier, 1979;Papazachos and Comninakis, 1971) is the driving mechanism for the fast extensional deformation in an almost N-S direction of the back arc region (Konstantinou et al., 2017;Kapetanidis and Kassaras, 2019).The westward prolongation of the North Anatolian Fault (NAF) into the Aegean forms the North Aegean Trough (NAT) that constitutes the active boundary between the Aegean microplate and the Eurasian lithospheric plate (inset of Figure 1).McKenzie (1972) showed that the northward motion of the Arabian plate pushes the smaller Anatolian plate that is moving westerly relative to the Eurasian plate along the NAF, with an average velocity of about ~24 mm/yr.An additional N-S deformation of ~11 mm/yr in the Aegean further enhances this motion, resulting in a total SW motion of ~41 mm/yr of the south Aegean relative to Eurasia.More recent studies (e.g.England et al., 2016;Bitharis et al., 2023) comply with this overall pattern.A large part of the deformation occurs seismically, as the Aegean shows a total seismic slip rate of the order of 20 mm/yr relative to Eurasia (Papazachos and Kiratzi, 1996).Seismicity is intense and continuous in the back arc region forming specific seismic zones striking almost E-W in the normal faulting environment and ENE-WSW to NE-SW for the dextral strike-slip and oblique f).aulting seismic zones (Papazachos et al., 1998).Conjugate faulting sinistral strike-slip faults are also present (Karakostas et al., 2003) accommodating less frequent M≥6.0 main shocks as well as moderate (M≥5.0)earthquakes.
Several strong (M≥6.0)historical and instrumental earthquakes have struck the NAT (Figure 1) since 360 BC (Papazachos and Papazachou, 2003).Their temporal distribution evidences catalog incompleteness at least until 1300 AD, and afterward missing events are the ones of M<6.5 (Kourouklas et al., 2018).Most recently, the 2014 M6.9 main shock ruptured the middle part of the NAT (Kiratzi et al., 2016), taking place in a previously identified seismic gap (Karakostas et al., 1986).The aftershock activity was astonishingly weak with a maximum magnitude aftershock of M4.9.However, the western termination of NAT against the Greek mainland appears to be relatively quiescent (Kourouklas et al., 2018(Kourouklas et al., , 2022)).The strike-slip faulting of the NAT is confirmed by the fault plane solutions of the stronger earthquakes in the most recent part of the instrumental era, when seismological networks were capable to provide adequate data, such as in the case of the 1983 seismic sequence (Rocca et al., 1985).The results from waveform modeling performed for the 1968 (M7.5;Pacheco and Sykes (1992)), 1981(M6.8) and 1983 (M6.6) earthquakes established the dextral strike-slip nature of the major faults in North Aegean with the axis of maximum extension, T, striking N-S and being nearly horizontal (Kiratzi et al., 1991).
Because it is difficult to determine the causative faults of the moderate (M≥5.0)earthquakes, our understanding of the seismogenic setting of the secondary faults in our study area is rather limited.An M5.8 main shock occurred in 2013 on one of the ENE-WSW trending parallel dextral strike-slip fault branches in Northern Aegean, in the continuation of the 1968 large (M=7.5)rupture (Drakopoulos and Ekonomides, 1972).Its rich aftershock sequence was scrutinized and the off-fault seismicity was perfectly explained with Coulomb stress changes when the parameters of friction and Skempton's coefficients attained values of µ>0.5 and B=0.0, respectively, implying high fault friction (Karakostas et al., 2014) Thus, clarifying the properties of the causative faults of moderate earthquakes is crucial to reveal the seismogenic structures and their relation to the major fault zones and evaluate the subsequent seismic hazard.
Microseismicity studies are not common in the study area, given that active faults are off-shore, and nearfault instruments are not in operation.The western termination of NAT was investigated by Hatzfeld et al. (1999) who confirmed the strike-slip motion which has probably been active since the Pliocene.This strike-slip motion is transferred into normal faulting (with the direction of principal extension being constant) in continental Greece.Accurate relocation of seismicity during 2011-2016 along the NAT yielded a pick of the focal depth distribution at 8 km diminishing down to 20 km (Konstantinou, 2017).The relocated seismicity defines distinctive clusters, one of which is located inside our study area, but the spatial distribution is quite disperse and thus not adequate to identify a specific activated structure.
On 26 September 2020 (at 22:50:24 UTC), an M w 5.3 earthquake occurred (right magenta circle in Figure 1), to the north of the NE-SW trending western portion of NAT and very close to its inferred trace, in a way that it would be considered as associated with failure on a fault patch along the NAT interface.It occurred in the middle of the night and was widely felt, alarming plenty of citizens even at distances greater than 200 km (EMSC felt reports).The GCMT (Global Centroid Moment Tensor; https://www.globalcmt.org/)solution suggests strike-slip faulting with one of the nodal planes striking NE-SW that complies with the sense of motion along the NAT (focal mechanisms configured with magenta compressional quadrants in Figure 1).The appreciable aftershock activity, however, persistently aligned in a NW-SE direction in agreement with the strike of the second nodal plane.These first observations provoked our interest in investigating the characteristics of the activated structure.Then, fifteen months later, a second moderate M w 5.4 earthquake occurred on 16 January 2022 (at 11:48:06 UTC) very close to the epicenter of the first 2020 shock (left magenta circle in Figure 1).The 2022 aftershock activity exhibited the same pattern as in 2020, and the second main shock's GCMT solution implies oblique faulting with the one nodal plane also striking NW-SE.In this study, we relocate the aftershocks and determine the fault plane solutions of the strongest among them.We further study the spatiotemporal evolution of the activity and perform statistical analysis through the application of the Epidemic Type Aftershock Sequence (ETAS) model.(after Leptokaropoulos et al., 2012), as compiled from the regional parametric earthquake catalog of the Seismological Station of Geophysics Department of the Aristotle University of Thessaloniki (http://geophysics.geo.auth.gr/ss/cata-logs_en.html; Aristotle University of Thessaloniki, 1981).The small light green and moderate orange circles depict the epicenters of 4.1≤M<5.0and 5.0≤M<6.0earthquakes, respectively, and magenta circles the epicenters of for the earthquakes investigated in this study.White stars display the epicenters of the historical M≥6.5 earthquakes between 1845-1900.Yellow stars depict epicenters of the M≥6.0 earthquakes occurred during the instrumental era.The available fault plane solutions of the M≥5.0 earthquakes as taken from Global Centroid Moment Tensor (GCMT; https://www.globalcmt.org/) are plotted as equal area lower hemisphere projections with the compressional quadrants colored according to their magnitude as before.Solid thick red line represents the main branch of the North Aegean Trough Fault Zone, whereas the parallel red arrows imply the right-lateral strike-slip motion.Inset map shows the main seismotectonic features of the Aegean region (solid red lines; KTFZ: Kefalonia Transform Fault Zone; RTF: Rodos Transform Fault; NAF: North Anatolian Fault).

Earthquake Relocation
An initial earthquake catalog was compiled by retrieving the recordings of the Hellenic Unified Seismological Network (HUSN) after routine analysis accomplished in the central Seismological Station of the Geophysics Department of the Aristotle University of Thessaloniki (http://geophysics.geo.auth.gr/ss/)for earthquakes with M≥2.5.Additional phase picking and initial location were performed by the authors of the present study to achieve the inclusion of lower magnitude events (even lower than M=1.0) that were detected by stations quite close to the activated area.In total, 841 earthquakes were initially located with magnitudes 0.6<M<5.4.For the earthquake relocation procedure, we used data from 14 seismological stations at epicentral distances up to approximately 150 km (Figures S1 & S2).The median azimuthal gap for all initial locations equals to 108 ο indicating a sufficient station azimuthal coverage, considering that the activated structures are offshore.
The first step towards improved focal coordinate estimation is relocation with the HYPOINVERSE code (Klein, 2002).This step requires the manually picked Pand S-phases, a local velocity model, the ratio of compressional to shear wave velocity (V p /V s ) and the corresponding station time delays.At first, we estimated the V p /V s through the Wadati method (Wadati, 1933) by taking the minimum number of phase arrival pairs (P-and S-) equal to 8 for each event, and the process resulted in a ratio equal to 1.74 ± 0.007.Similar values were also estimated for the broader Aegean region in previous studies (e.g.Karamanos et al., 2007;Mesimeri et al., 2018;Andinisari et al., 2020;Karakostas et al., 2021).Next, we proceeded to the determination of a crustal model using the VELEST algorithm (Kissling et al., 1994) by testing multiple published models as reference ones (e.g.Akyol et al., 2006;Karabulut et al., 2006;Konstantinou, 2018).All resulting models after multiple iterations consistently converged to a similar model (apart from the first few kilometers where differences could be found) very close to the model of Karabulut et al. (2006).We thus decided to adopt this model as a reference one (Figure S3) and use our derived model (Table 1) for the relocation procedure.We also added the appropriate station corrections calculated with the VELEST application to incorporate lateral inhomogeneities in the 1-D crustal model.
We then run the hypoDD program (Waldhauser, 2001) that employs the double difference algorithm (Waldhauser and Ellsworth, 2000), to improve the accuracy of the focal coordinates as they were obtained from HY-POINVERSE.Travel time differences between manually picked phases in the earthquake catalog were calculated, and then we kept event pairs with at least eight (8) observations, resulting to 32,956 P-phase pairs and 25,870 S-phase pairs (13 links per pair on average) in total.The differential times data set was analyzed using five sets with five iterations on each set gradually decreasing the residual threshold (secs) and the maximum distance (km) between catalog linked pairs.The hori-Depth (km) Vp (km/s) 0.0-2.03.56 2.0-4.0 4.32 4.0-7.05.43 7.0-14.06.08 14.0-20.06.22 20.0-26.0 6.55 26.0-29.06.8 ≥30 7.29 Table 1 P-wave velocity (V p ) model adjusted for the study area zontal and vertical errors between the initial and the relocated catalog indicate a substantial reduction of location uncertainty (Figure S4).The final catalog contains 817 relocated earthquakes out of 841 (~97%) that constituted the initial data set.The 477 earthquakes occurred before the 2022 main shock and the 340 are the aftershocks of the second sequence.

Fault plane solutions
Fault plane solutions of 12 earthquakes with M≥4.0, which occurred during the whole study period, were calculated.One of them occurred slightly later, in April 2022, but in the same area.The moment tensor inversions were conducted using the Grond software (Heimann et al., 2018) which operates under the Pyrocko toolbox framework (Heimann et al., 2017).The methodology applied within Grond aims to minimize the misfit between synthetic and observed data by implementing a Bayesian bootstrapping inversion in parallel bootstrap chains.In the present case, the bootstrap algorithm was applied to minimize the L2 norm misfit for 22000 iterations, 3000 to uniformly sample the solution space and 19000 for the direct sampling.The number of parallel bootstrap chains was set equal to 200.We used the recordings of the regional broadband seismological stations of HUSN in distances ranging between 50 and 300 km (Figure S1).Green's functions were calculated for the local velocity model (Table 1) using the QSEIS program (Wang, 1999).A point source model was considered to perform the Bayesian optimization for a deviatoric moment tensor.Lastly, bandpass filters were applied, using frequencies between 0.04 and 0.11 Hz (the corresponding window for each earthquake is presented in Table 2) and the data fitting was carried out in the time domain.The NW-SE striking nodal planes of the focal mechanisms (Table 2; Table S1) exhibit predominantly sinistral strike-slip sense of motion, ranging from almost pure strike-slip to oblique faulting.Table 2 Fault plane solutions calculated for eleven earthquakes of the relocated catalog and one that occurred one month later (No 12).The focal parameters provided correspond to the relocated ones.In the last column, the square of the misfit of each solution is reported.

The Epidemic-Type Aftershock Sequence (ETAS) model
Short-term temporal properties of both the 2020 and 2022 seismic excitations were studied through the ETAS stochastic model (Ogata, 1988(Ogata, , 1998)).The temporal ETAS model (Ogata, 1988) expresses the seismicity occurrence rate density as the summation of two parts, the constant background seismicity rate, µ, assumed timeindependent, and the occurrence density rate of triggered earthquakes, λ i (t), and is given by a conditional intensity function, λ(t): where K and α are the aftershock productivity parameters, c and p are the parameters of the temporal aftershock decay rate of the modified Omori law (Utsu, 1961), m i is the magnitude of each earthquake occurred at time t i and m c is the completeness magnitude.The parameter K represents the intensity of aftershocks generation above m c triggered by an earthquake with m = m c , whereas α parameter describes the efficiency of earthquakes in triggering their own aftershocks.Large α values [1.3-3.1]indicate that large magnitude earthquakes trigger a large number of aftershocks, whereas small α [0.35-0.85]implies relatively higher triggering capabilities of small earthquakes.This means that Mainshock-Aftershock sequences typically tend to have larger α values, dominated by the main shock magnitude, while a swarm-like activity is characterized by small α values (Hainzl and Ogata, 2005;Ogata, 1992).
The exponent p of the modified Omori law controls the aftershocks decay rate, where the decay is becoming faster as the value of p is increasing.The parameter c is linked with the short-term aftershock incompleteness soon after the occurrence of a main shock (t = 0), aiming to avoid singularities at occurrence times very close to t = 0.The five parameters of the model were estimated via the Maximum Likelihood Estimation (MLE) method using the simulated annealing technique proposed by Lombardi (2015) and implemented through an algorithm which is part of the SEDA package (Lombardi, 2017).
For the model evaluation, a residual analysis was performed Ogata (1988Ogata ( , 1998)), which offers a qualitative model evaluation through visual display.Specifically, the occurrence times, t i , were converted into transformed times, τ i , according to the function: (2) Transformed times, τ i , express the number of earthquakes that are expected to occur in the time interval [0, t i ].If the estimated model adequately describes the temporal seismicity evolution, then transformed data, so-called residuals, behave like a stationary unit rate Poisson process.Otherwise, the transformed process will show some systematic departure from the linear Poisson process.Positive and negative departures in-dicate that the estimated model is under-and overpredicting the observed seismicity, respectively.In order to evaluate whether the transformed times, τ i , are described by the Poisson process, the one-sample Kolmogorov-Smirnov test (KS1 test; Massey (1951)) was applied.More specifically, if τ i are modelled by the unit rate Poisson process, then the transformed earthquake interevent times, t i+1 − t i , should be independent and identically drawn from an exponential distribution (Llenos and Michael, 2013).The KS1 test is implemented under the null hypothesis that the earthquake interevent times are following the exponential distribution based on the p-value returned by the test, compared with the 0.05 significance level.If p-value is greater (or lower) than 0.05 then the null hypothesis can either be rejected or accepted.

Spatiotemporal evolution of seismicity
The number of precisely relocated earthquakes equals to 817, with 477 of them belonging to the period starting from the initiation of the seismic excitation (26 September 2020) until the start of the second seismic sequence (16 January 2022) and the rest 340 constituting the second aftershock sequence.Regarding pre-seismic activity, only two (2) foreshocks were detected prior to the first main shock with the largest one occurring 4 hours in advance (Table 2).Our interpretations on aftershock activity and the associated active faults can be ensured by the quality criteria that the relocation procedure has fulfilled concerning the uncertainty estimation.These criteria do not prevent from systematic bias, but with our efforts to approach the velocity structure as best as possible we are confident that this bias could be considered small.We are interested in how seismicity evolved spatially or temporally in the activated area on the temporal scale of several months, for which that this activity persisted.We started our relocated data set since 1 January 2020, almost 10 months before the first main shock occurrence.We may observe that seismicity is quite sparse with a couple of 2.0<M<2.9shocks located close to the first main shock epicenter (green dots since the beginning of the space-time plot up to 26 September 2020, in Figure 2a) but not close in time for being considered foreshock activity as part of the nucleation phase.An M w 4.3 earthquake, instead, occurred about 4 hours before and in close distance with the main shock (Table 2).Intense aftershock activity followed the main shock on 26 September 2020, extended in an area more than 15 km long, much larger than the rupture length of an M w 5.3 main shock, (roughly 6 km) as prescribed by empirical scaling laws (Wells and Coppersmith, 1994;Thingbaijam et al., 2017).As Figures 2a & 3a show, seismicity is mainly concentrated on a 20 km length prolonged zone, with a general strike of 145°(thick dashed line in Figure 3a).The epicentral distribution of the two distinct clusters present a spatial overlap of slightly over 50%.The determined fault plane solutions unveil the dominance of strike-slip faulting in the study area given that pure strike-slip (e.g., No 1 & 6, Table 2; Figure 3a) and strike-slip with a normal component (2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12) focal mechanisms are present.The complex faulting patterns coming into play are revealed by the presence of either pure normal (No 5, Table 2; Figure 3a) or oblique (normal with a strike-slip component) focal mechanisms (No 3 & 7, Table 2; Figure 3a).
In Figures 4 & 5 the relocated seismicity of the two activated structures is shown separately aiming to detail the properties of each seismic sequence.Regarding the 2020 M w 5.3 sequence (Figure 4a), the epicentral alignment (NW-SE) of the aftershocks is in good agreement with the strike of the one nodal plane of the main shock focal mechanism.The total length of the activated area is approximately 15 km, which significantly exceeds the estimations of frequently used empirical relationships between main shock magnitude and the causative fault length (Wells and Coppersmith, 1994;Thingbaijam et al., 2017).The duration of the intense aftershock activity was short, with over half of the aftershocks occurring within the first week after the main shock (Figure 3).Moreover, all M≥4.0 aftershocks took place in less than 3 days (Table 2).An interesting feature of this seismic sequence is also the occurrence of a few foreshocks with the strongest among them having a magnitude of M w 4.4 (Table 2), placed very close to the main shock epicenter (Figure 2) although in deeper parts of the seismogenic layer.
The aftershock activity during the first 2 days defines a seismogenic zone of 9 km (13 -22 km) (Figure 4b).As time advances, however, the inclusion of all earthquakes (Figure 4c) evidences an expansion of the aftershock zone in more shallow depths (~5km).Normal to the strike cross section encompassing events within the larger part of the aftershock distribution (3 km of either side of the cross section) provides a complete picture of the activated fault plane (Figure 4d).The total depth extent of the aftershock zone ranges from 6 to 24 km (Figure 2b).
Shifting our focus to the second seismic sequence that emerged from the strongest event included in our catalogue (M w 5.4) at the beginning of 2022 we might observe that the epicentral alignment of the aftershocks is analogous to the previous case (NW-SE) however their spatial extent is developed further to the north.The total length of the activated area, as inferred by the aftershock epicentral distribution is approximately equal to 10 km (Figures 2a, 5a).The aftershock activity of the first 2 days occupies a depth range of 9 to 21 km, inside a 12 km extent (Figure 5b), with only minimal lowmagnitude events outlying.Figure 5c, which includes all events at a distance of 1.5 km either side of the normal cross section, appears to be very similar to the one of Figure 4b, highlighting the fast decay of the aftershock rate after the first few days.Moreover, the wider  normal cross section (Figure 5d) reveals the close proximity of the strongest aftershocks (11-13 km in depth).

Temporal features from the ETAS model fitting
The temporal ETAS parameters (µ, K, α, c and p) were estimated through the MLE method after the determination of the magnitude of completeness, m c , for the entire initial earthquake catalog covering the period 2019-2023.The m c was identified through the Goodnessof-Fit method (GFT; Wiemer and Wyss (2000)), considering the 95% confidence level of residuals (Figure S5) and was found equal to m c =1.9 (residuals value equal to 3.66%) resulting to a data set of 470 earthquakes with m ≥ m c and a b-value equal to 0.78 (b=0.78).The model was first applied in the entire initial earthquake cata-  Table 3 Temporal ETAS parameters estimates (µ, K, α, c and p) for the periods 2019-2022, January 2020-June 2021 and July 2021-December 2022, along with the respective number of observations and the p-values of the one sample Kolmogorov-Smirnov goodness-of-fit test (KS1) between the earthquake interevent times and the exponential distribution.
log of the period 2019-2022 and then in two additional and distinctive sub-periods, namely from January 2020 to June 2021 and from July 2021 to December 2022, aiming to compare the temporal properties of the two seismic sequences.The selection of each sub-period was based on the adequacy of the number of earthquakes before the occurrence of each excitation, in which the parameters were estimated, representing the learning phase of the ETAS model application.The calculated parameter values are given in Table 3.
The ETAS application for the entire period (2019-2022) ascertained that the model adjusts well the observed earthquake rate, with slight discrepancies soon after the occurrence of the 2020 M w 5.3 and 2022 M w =5.4 main shocks (Figure 6a).The good fit of the estimated model in respect to the observations is highlighted by the residuals analysis application (Figure 6b), since the lines that depict the expected and the observed number of events against the transformed times (red and black lines in Figure 6b, respectively) almost coincide.This latter fact is confirmed by the result of the KS1 goodness of fit test shown in Table 3. Specifically, the calculated p-value of the test is equal to 0.62, much larger than the 0.05 confidence level.The estimated parameters are acquiring values typical for Mainshock-Aftershock sequences (Ogata, 1992).In more detail, the observed seismicity during the period from 2019 to 2022 is characterized by a very low background rate equal to µ=0.046 event/day.This means that almost 88% of the total number of earthquakes are offsprings of the 2020 M w 5.3 and 2022 M w =5.4 sequences, and only 68 out of the 470 are assumed to be independent.Furthermore, the productivity parameter, α, was estimated equal to α = 1.78, indicating also Mainshock-Aftershock type of sequence, since Mainshock-Aftershock activity is typically characterized by α values ranging between [1.3-3.1],whereas swarm-like activity is ascribed to values ranging between [0.35-0.85](Hainzl and Ogata, 2005;Ogata, 1992).The parameters expressing temporal characteristics, c and p, attain expected values, in comparison with those reported in previous studies ranging from 0.03 to 0.07 and from 1.16 to 1.25, respectively (Chu et al., 2011;Kourouklas et al., 2020).
Placing emphasis on the two sub-periods (January 2020-June 2021 and July 2021-December 2022), when the M w 5.3 and M w 5.4 main shocks occurred, the small discrepancies between the observed earthquake rates and the modeled ones are becoming more explicit (Figures 6c and 6e, for the January 2020-June 2021 and July 2021-December 2022, respectively).Specifically, it is observed that the estimated models for both periods slightly underestimate the earthquake rates soon after the occurrence of the M w 5.3 and M w 5.4 earthquakes.This result is also visible from the residual analysis plots (Figures 6d & 6f), in which the observed transformed times (red lines in Figures 6d & 6f) slightly deviate from the unit rate Poisson process (black lines in Figures 6d &  6f).However, the p-values of the KS1 test (p-value=0.60 and 0.57 for the periods January 2020-June 2021 and July 2021-December 2022, respectively; Table 3) are again much larger than the significance level (0.05), suggesting a good performance of the temporal ETAS model to the data from a statistical point of view.
Comparison of the parameter estimates of the two sub-periods (Table 3) shows that throughout the first one (January 2020-June 2021), during which the M w 5.3 earthquake occurred (September 26 th of 2020), the background rate is estimated equal to µ = 0.047 event/day, a value almost equal to the entire period estimate).On the contrary, background rate is found quite smaller (µ = 0.037 event/day) during the period from July 2021 to December 2022, in which the M w 5.4 earthquake occurred (January 16 th of 2022).Both estimated values are again indicating a clear Mainshock-Aftershock activity as well as the estimated values of the productivity parameter, α, (α = 1.75 and 1.86 for the 1 st and the 2 nd sub-periods, respectively; Table 3) with the one attached to the second sub-period being higher.Additionally, the estimated value of p parameter of the modified Omori law is larger in the second sub-period's application (p = 1.17 instead of p = 1.25 for the 1 st period) indicating the smaller duration of the second sequence.As the westward prolongation of the North Anatolian Fault (NAF) system into the Aegean Sea, the North Aegean Trough (NAT) is characterized by dextral strikeslip faulting.The intense N-S extension of the back-arc area due to the roll back of the subducted slab in southern Aegean, resulted to a complex transtensional basin dated back to early Pliocene or Pleistocene.The conjunction of major strike-slip systems with a wide variety of secondary structures, especially at fault zone tips or between linked structures, the damage zones (e.g Petit and Barquins, 1988;Kim et al., 2004;Kim and Sanderson, 2006) is very commonly met.The reason behind their development can be attributed to stress concentrations (e.g Cox and Scholz, 1988) or to host displacement alterations along the fault zones (Kim et al., 2000).The North Aegean area is characterized by a broad spectrum of strike-slip damage patterns (Figure 7), with the most important and relevant to our study area being: a) Horsetail structures: The western endpoint of the NAT is characterized by a right-stepping horsetail structure (red patch in Figure 7) expressed through numerous oblique splays stemming from the main strike-slip zone forming a number of transtensive basins (Sakellariou et al., 2017).b) Positive or negative flower structures: Negative flower structures are found inside most of the subbasins bounded by the oblique splays.They consist of opposing-dip normal fault structures and merge into a single strike-slip fault at deeper parts (upper left part rectangle in Figure 7).Positive ones are also present but outside the area under study, along the south Marmara and the Ganos fault segments, expressing transpressive structures owing their formation to strike-slip faults with a reverse component (Rodriguez et al., 2023).c) Conjugate faults: They are faults that intersect with the main structure at a high angle (typically more than 60°) and exhibit a sense of displacement opposite to the dominant one, with which they are often not connected.In our case, they represent the sinistral counterparts of the prevalent dextral strike-slip faults that dominate the northern Aegean area.Conjugate faults may be combined with branch faults (see below) and create block rotation (e.g.Nicholson et al., 1986;Kim et al., 2003).d) Branch faults: They represent shear fractures having similar sense of motion as the main strike-slip zone (dextral strike-slip in our case, blue patch in Figure 7).They can act together with the main zone to create splays (e.g.Kim et al., 2004) or with other structures and form more complex patterns.
The diversity characterizing the termination of strikeslip fault systems, especially on such a large scale such as the NAT makes it very challenging to uncover its fine details, due to their complex nature and high intercorrelation.Typical examples are the fault segments associated with the two main shocks of this study which are not included in the simplified tectonic map of the North Aegean based on the interpretation of the air gun lithospheric profiles by Papanikolaou et al. (2006)   Based on our analysis, the activated structures hosting the 2020 M w 5.3 and 2022 M w 5.4 main shocks can be attributed to conjugate faults at the termination of the main strike-slip zone, having a considerable vertical component consistent with the extensional character of the NAF termination (Figure 7).The primary indications leading to this suggestion are the determined focal mechanisms of the main shocks and the stronger aftershocks (Table 2), which indicate an intricate stress field consisting of strike-slip and normal faulting style.Adding to this, the strike-slip moment tensor solutions exhibit a left-lateral displacement, as opposed to the dextral strike-slip motion characterizing the NAF.Moreover, the spatiotemporal distribution of aftershocks (Figures 3-5) further demonstrates the high angle (almost perpendicular) in which the activated structures lie in comparison to the NAF (Figure 1).The scenario of a transtensional flower structure cannot be ruled out given that such systems are closely related and common to conjugate fault systems.A valid reason for this argument can be ascribed to the relative differences of the focal depths of the two sequences (Figure 2).Even though the spatial distribution of both aftershock sequences appears to extend in more or less the same area in map view (Figure 3), the mean focal depths from the aftershocks originated from the 2020 M w 5.3 sinistral strike-slip main shock are found in deeper parts compared to those from the 2022 M w 5.4 oblique (λ=-43°; normal with sinistral strike-slip component) one.In this case, a depressed area, or even a pull-apart basin (on a larger scale) is most commonly formed.Bathymetric maps (Papanikolaou et al. (2006) and references therein) cannot provide us with a definite answer however our study area is located at the northern margins of the main North Aegean basin and more specifically in an area where a bulge is disrupting the almost linear NE-SW oriented edges of the basin.
Statistical analysis of the short-term clustering features of the 2020-2022 excitations highlights their Mainshock-Aftershock nature.The estimated temporal ETAS parameters show non-significant background seismicity rate changes for both the average model (2019-2022) and the two sub-period models, which are referring to the 2020 and 2022 sequences, indicating that the excitation is driven by the regional tectonic loading.An interesting remark arises from the estimated parameters of the second sequence (2022) modeling, in which the temporal parameters (c and p) are larger.These values indicate a faster decay rate compared to the first sequence's aftershocks, even though the magnitude of the 2022 main shock is larger than the 2020 one.This remark could be likely explained by the fact that the study area was recently activated and consequently rather stress relaxed.
Moreover, the low background rate, the estimated productivity parameters and the fast aftershock decay rates in all evaluated models provide an indirect insight into the properties of the activated fault segments, qualifying them as rather weak.This property characterizing fault segments in the broader area (including the main branch of NAT fault zone) has already been suggested by physics-based earthquake simulator results (Kourouklas et al., 2021).

Conclusions
We thoroughly investigated the seismic sequences of two moderate main shocks occurred in the region north of NAT on 26 September 2020 (M w 5.3) and 16 January 2022 (M w 5.4), by means of their spatiotemporal distribution, moment tensor solutions and statistical analysis.We constrained the dimensions of the rupture areas (10x10 km 2 for the 2020 main shock and 11x10 km 2 for the 2022 main shock) based on the relocated early aftershocks.Fault plane solutions highlight the complexity of the faulting patterns, with most of them exhibiting mainly sinistral strike-slip faulting, with some oblique cases (strike-slip with a normal component) being also present.
All things considered, we attribute the occurrence of the investigated seismic sequences as a result of secondary faulting related to the termination of an active strike-slip plate boundary, expressed as transtensional deformation.We may consider these main shocks as a beneficial example of better recognition of the complex behavior characterizing the NAT, but through the lens of seismic hazard, we weigh this kind of events much less prominent than the overall seismic hazard of the broader area.

Figure 1
Figure 1 Epicentral distribution of M≥4.1 crustal (0≤h≤20 km) earthquakes occurred in northern Aegean Sea from 1975 to 2022 (after Leptokaropoulos et al., 2012), as compiled from the regional parametric earthquake catalog of the Seismological Station of Geophysics Department of the Aristotle University of Thessaloniki (http://geophysics.geo.auth.gr/ss/cata-logs_en.html; Aristotle University of Thessaloniki, 1981).The small light green and moderate orange circles depict the epicenters of 4.1≤M<5.0and 5.0≤M<6.0earthquakes, respectively, and magenta circles the epicenters of for the earthquakes investigated in this study.White stars display the epicenters of the historical M≥6.5 earthquakes between 1845-1900.Yellow stars depict epicenters of the M≥6.0 earthquakes occurred during the instrumental era.The available fault plane solutions of the M≥5.0 earthquakes as taken from Global Centroid Moment Tensor (GCMT; https://www.globalcmt.org/) are plotted as equal area lower hemisphere projections with the compressional quadrants colored according to their magnitude as before.Solid thick red line represents the main branch of the North Aegean Trough Fault Zone, whereas the parallel red arrows imply the right-lateral strike-slip motion.Inset map shows the main seismotectonic features of the Aegean region (solid red lines; KTFZ: Kefalonia Transform Fault Zone; RTF: Rodos Transform Fault; NAF: North Anatolian Fault).

Figure 2
Figure2(a) Spatiotemporal distribution of the relocated seismicity since 1 January 2019 up to 1 November 2022.Yellow stars represent the two main shocks, red circles the 4.0≤M w ≤4.9, orange circles the 3.0≤M≤3.9,green circles the 2.0≤M≤2.9,and small black dots the M≤2.0 earthquakes.(b) Histogram of the focal depths with different colors denoting the two aftershock sets.

Figure 3
Figure 3 Seismicity of the study area, since the beginning of 2019 until the end of 2022.(a) Stars indicate the relocated epicenters of the main shocks, whereas circles denote either the relocated aftershock locations or the background seismicity.All epicenters are scaled in compliance with the corresponding earthquake magnitudes and colorcoded according to the temporal scale at the bottom of the figure.Fault plane solutions determined in this study (Table 2) are also shown as lower hemisphere equal area projections, with the compressional quadrants colored in blue.(b) Temporal evolution of seismicity for the 2019-2022 period versus magnitude.Blue, yellow and cyan circles and shaded areas indicate the pre-2020, 2020 and 2022 seismic activity.Red solid line depicts the cumulative number of earthquakes for the same period.

Figure 4
Figure 4 a) Spatial distribution of the relocated aftershocks of the M w 5.3 main shock on 26 September 2020 (red star).The epicenters of the aftershocks are shown as circles of different sizes and colors according to their magnitude.The black dashed line shows the general strike of the epicentral distribution, whereas the constant one signifies the normal to the strike cross sections including events b) within 1.5 km either side of the section, during the first 2 days after the main shock, c) within 1.5 km for the entire duration and d) 3.0 km either side of the section, again for the entire relocated catalog.

Figure 5
Figure 5 a) Spatial distribution of the relocated aftershocks of the M w 5.4 main shock on 16 th January 2022 (red star).The epicenters of the aftershocks are plotted with circles of different sizes and colors according to their magnitude.The black dashed line shows the general strike of the epicentral distribution, whereas the constant one signifies the normal to the strike cross sections including events b) within 1.5 km either side of the section, during the first 2 days after the main shock, c) within 1.5 km for the entire duration and d) 3.0 km either side of the section, again for the entire relocated catalog.

Figure 6
Figure 6 Observed (black lines) and expected by the estimated temporal ETAS model (red lines) cumulative earthquake number against ordinary (a,c,e) and transformed (b,d,f) time for the periods 2019 -2022 (a & b, respectively), January 2020 -June 2021 (c & d, respectively) and July 2021 -December 2022 (e & f, respectively).Magnitudes of earthquakes with m ≥ m c (right y-axis) against time are shown as pink dots in all subplots.

Figure 7
Figure 7 Strike-slip damage patterns (schematic illustrations) and their connection to the study area (simplified map).The green and yellow stars indicate the epicenters of the main shocks of this study.