Seismicity supports the theory of incipient rifting in the western Ionian sea, central Mediterranean

The present work focuses on earthquake locations and seismogenic stress in the eastern offshore of Sicily, a sector of the central Mediterranean region where geophysical information available is not good enough, yet, for proper geodynamic modeling. I have applied to an updated seismic database of the study area a Bayesian non-linear hypocenter location method already proven to be more effective than linear methods when the recording network geometry is poor, like in the present case. Then, I have selected from literature and official catalogs the local earthquake focal mechanisms computed by waveform inversion, and inverted them for stress tensor orientations. The results confirm the main finding of the previous investigations, i.e. that NW-trending convergence between Africa and Eurasia is a main source of tectonic stress in this area, however they also furnish evidence of additional tectonic factors locally acting together with convergence. In particular, extensional dynamics are detected inside the convergence-related compressional domain: these are characterized by a minimum compressive stress oriented SW-NE (perpendicular to convergence) and can be related to the rifting process (opening SW-NE) detected by previous investigators at the southwestern edge of the Ionian subduction slab. The findings of the present study may also concur to answer several open questions left by previous investigators.


Earthquake space distributions
The Figure 3a shows the earthquakes of local magnitude over 2.5 occurring between 1985 and 2018 at depths less than 70 km in the calabro-sicilian region, according to the Italian national seismic database (http://istituto.ingv.it/lingv/archivi-e-banche-dati/). Hypocenter locations of the used database are estimated with a 1-D velocity structure of all Italy, which is a reasonable approach for bulletin activity. I use these locations here (i) to furnish an introductory overall view of seismicity in the study region and (ii) as starting data for more accurate locations to be performed in the sector of our greatest interest in the present study (shadowed internal rectangle).
I have focused my attention on the western Ionian because, as explained in the previous Sections, the geophysical knowledge in this offshore sector is still relatively poor and its exploration can be decisive for answering several open geodynamic questions in the region. I have relocated the hypocenters of earthquakes of the shadowed sector of Figure 3a, after integration of P-and S-wave readings of the Italian database with those available from databases of local seismic networks operating in Sicily and Calabria during 1985. For hypocenter relocations, I have selected the events for which a minimum of 12 P+S arrival times were available, and used the Bayesian non-linear location algorithm named Bayloc [Presti et al., 2004[Presti et al., , 2008. As well known from the literature [Lomax et al. 1998[Lomax et al. , 2000Lomax and Michelini, 2001;Husen and Smith, 2004;Presti et al. 2004Presti et al. , 2008Lippitsch et al., 2005; among others], the non-linear probabilistic location methods furnish more accurate estimates of hypocenter locations and relative errors compared to linearized methods when the network geometry is not optimal: this is the situation of our offshore study sector (Figure 3b). Starting from seismic phase arrival times at the recording stations, Bayloc computes for an individual earthquake a probability cloud marking the hypocenter location uncertainty. Then, Bayloc estimates the spatial distribution of probability relative to a set of earthquakes by summing the probability densities of the individual events. This method has been shown to help detection of seismogenic structures through better hypocenter location and more accurate estimation of location errors compared to linearized methods [Presti et al., 2008].

5
Seismicity supports rifting in the West Ionian Bayloc's locations have been performed in a 3D velocity structure estimated for the study region by Orecchio et al. (paper in preparation) who applied the method used by Orecchio et al. [2011] and previously proposed by Waldhauser et al. [1998Waldhauser et al. [ , 2002. This method is based on: (i) integration of different types of velocity data available from the literature (seismic profiles, earthquake tomography, surface wave inversion, Moho depth maps, etc.); (ii) LET inversion where this is allowed by the available P-and S-wave arrival data. Velocity data from literature were taken from Kennet et al. [1995], Tesauro et al. [2008], Orecchio et al. [2011, and references therein], Neri et al. [2012], Scarfì et al. [2018]. The 3D velocity structure used here for hypocenter locations covers the depth range 0-300 km in the whole area concerned by seismic rays travelling from hypocenters of the study sector to recording stations ( Figure 3b).
Epicenter distributions relative to different hypocentral depth ranges and a SW-NE vertical section of hypocenters computed by Bayloc for earthquakes occurring during 1985-2018 in the study area are displayed in

Seismogenic stress inversion
With the purpose of analyzing seismic faulting styles and related tectonic stress in the area of interest of the present study, I have collected from literature and international catalogs the waveform-inversion focal mechanisms of the earthquakes occurring in the area of Figure 4. I have limited the selection to seismic events of magnitude over 2.5 occurring in the period 1977-2018 at depths less than 70 km. The map of these focal mechanisms is shown in Figure 5, the list of their parameters is furnished in Table 1 In the dataset of Figure 5 and Table 1, the focal mechanisms computed by the CAP method [Li et al., 2007;D'Amico et al., 2010;Presti et al., 2013;Orecchio et al., 2014;Polonia et al., 2016;Totaro et al., 2016] are affected by errors of the order of 8-10 degrees. The literature and the bibliographic sources of the other focal mechanisms of Figure 5 and Table 1 indicate that these are typically characterized by fault parameter errors of the order of 10-15 degrees [see, e.g., Helffrich, 1997;Frohlich and Davis, 1999;Pondrelli et al., 2006;Hjörleifsdóttir and Ekstrom, 2010]. Then, the fault parameter errors of the solutions displayed in Figure 5 are, in general, smaller than those of focal mechanisms computed by inversion of P-onset polarities in areas of critical network geometry like ours [D'Amico et al., 2011;Presti et al., 2013;Scarfì et al., 2013;Musumeci et al., 2014]. The overall level of uncertainty of focal mechanisms of the dataset of Figure 5 makes it suitable for application of the method by Gephart and Forsyth [1984] for calculating seismogenic stress directions in the study region (Figure 5a-d). This method searches for the stress tensor showing the best agreement with the available focal mechanisms (FMs). Four stress parameters are calculated: three of them define the orientations of the main stress axes; the other is a measure of relative stress magnitudes, R = ( 2 -1 )/( 3 -1 ), where 1 , 2 and 3 are the values of the maximum, intermediate and minimum compressive stresses, respectively. In order to define discrepancies between the stress tensor and observations (FMs), a misfit variable is introduced: for a given stress model, the misfit of a single focal mechanism is defined as the minimum rotation about any arbitrary axis that brings one of the nodal planes, and its slip direction and sense of slip, into an orientation that is consistent with the stress model. Searching through all orientations in space by a grid technique operating in the whole space of stress parameters, the minimum sum of the misfits of all FMs available is found. The confidence limits of the solution are computed by a statistical procedure described in the papers by Parker and Mc Nutt [1980] and Gephart and Forsyth [1984]. The size of the average misfit corresponding to the best stress model provides a guide as to how well the assumption of stress homogeneity is fulfilled [Michael, 1987]. In the light of results from a series of tests carried out by Wyss et al. [1992] and Gillard et al. [1996] to identify the relationship between FM uncertainties and average misfit in the case of uniform stress, I will make the following assumptions. I assume that the condition of homogeneous stress distribution is fulfilled if the misfit, F, is smaller than 6°, and that it is not fulfilled if F>9°. In the range 6°<F<9°, the solution is considered as acceptable, but may reflect some heterogeneity.
The advantage of using Gephart and Forsyth's [1984] method instead of other more recent stress inversion methods [such as, for example, Arnold and Townend, 2007;Vavrycuk, 2014;Karagianni et al., 2015] is that the former is more conservative concerning the relative orientation of seismogenic stress and seismic dislocation surface. In this connection, caution is appropriate in the present study because I do not make any assumption here concerning the date of formation of the faults which generated the earthquakes of the dataset. The tectonic stress orientation may have changed since the date of formation of the fault producing the study earthquake, therefore I must consider a relatively wide range of possible angles between the today-acting stress and the fault surface. This more conservative approach is better represented by the method of Gephart and Forsyth [1984]. The effectiveness of this method in several conditions, also compared to more recent methods, is well documented in the literature [see Hardebeck and Hauksson, 2001;Maury et al., 2013;Karagianni et al., 2015; among others].
I report in Figure 6  i.e. same sector of plot c, depth range 30-70km. See Table 2 for numerical values of stress inversion results. Table 2. Stress tensor inversion of earthquake focal mechanisms performed for the earthquake sets indicated in Figure 5 and described in the text. N is the number of earthquakes (= focal mechanisms) belonging to the inversion set.  Plots (c) to (e) display the stress inversion results obtained for the dataset of plot (b) by using the stress inversion methods by Gephart and Forsyth [1984], Vavrikuk [2014] and Arnold and Townend [2007], respectively.

Discussion
In the most recent analysis of regional seismicity and stress fields, Totaro et al. [2016] proposed a geodynamic scheme highlighting extensional processes of the Apennine-Maghrebian chain occurring inside the overall compressional domain due to Africa-Eurasia convergence (see, in particular, the Figure  domain [Totaro et al., 2016, among others]. Figure 5c and Table 2 (set c) display the results of stress inversion performed in the area of incipient rifting proposed by Polonia et al. [2016Polonia et al. [ , 2017, approximately located between the Alfeo-Etna and Ionian Fault Systems. In this case, the F-value drops to 6.7° that is quite smaller than the values of the larger datasets of Figures 5a and 5b, but still larger than the value of 6° assumed as approximate upper bound of values corresponding to stress homogeneity (see previous Section). Therefore, stress is moderately heterogeneous in the area between the Alfeo-Etna and Ionian Fault Systems. Here, the best model of stress coming from inversion is characterized by a NW-trending, horizontal 1 matching well with the direction of convergence of Africa and Eurasia in this part of the plate margin. However, the 95% confidence limits of stress orientations reveal that 1 orientation is practically unconstrained from horizontal NW-SE to vertical. This may suggest that some extensional process opening SW-NE acts together with NW-trending plate convergence in this sector.
A close analogy can be noted between the stress confidence limits of Figure 5c (real earthquakes) and those obtained by inversion of synthetic focal mechanisms (Figure 6c). This suggests that the focal mechanisms available between the Alfeo-Etna and Ionian Fault Systems (Figure 5c (ii) an extensional stress with SW-NE opening direction which can be plausibly related to the rifting process hypothesized by Polonia et al. [2017]. I have reported in Figure 7 the individual misfits of the earthquakes of Figure 5c with respect to the corresponding stress solution, as a function of focal depth. The best model of stress in 5c is compressional with SE-NW 1 , which recalls regional stress related to Africa-Eurasia convergence. According to the plot of Figure 7 plate convergence seems to be the only geodynamic process active at depths between 30 and 45 km, approximately. In fact, at these depths, nearly all the individual misfits are lower than 10° (average error of focal mechanisms in the dataset). Conversely, the individual misfits of the earthquakes shallower than 30 km (in half cases exceeding 10°; Figure 7) suggest some degree of stress heterogeneity in the depth range 0-30 km. Based on the information furnished by the plot of Figure 7, I have performed a stress inversion run excluding from the dataset of Figure 5c all the earthquakes shallower than 30 km. The results of this additional run are shown in Figure 5d and Table 2 set d. The value of F (3.3°) and the small confidence limits of stress orientations indicate that stress is homogeneous and well constrained at depths between 30 and 45 km in the sector of Figure 5d. In this framework, the best model of stress characterized by a NW-SE sub-horizontal 1 , indicates that the action of plate convergence is dominant in this depth range in the area of incipient rifting identified by Polonia et al. [2017]. According to the results of Figure 5c-d and to the plot of Figure 7, the extensional processes associated to rifting are to be confined to the upper ca. 30 km in this area. Unfortunately, the low number of focal mechanisms available above the depth level 30 km (only 10) does not allow a separate stress inversion and this operation has to be postponed to a moment when additional data will be available. It is also worth mentioning that the analysis of the individual misfits of the earthquakes of Figure 5c as a function of magnitude (not reported graphically for conciseness) evidences that all the events of magnitude over 4.0 (maximum magnitude of 4.5) can be imputed to the convergencerelated regional stress (misfits less than 10°). On the other hand, the weaker events (magnitude less than 4.0) show in several cases misfits larger than 10°, what means that they in part reflect stress heterogeneities due to more localized processes, such as rifting. In this case, too, the relatively low number of focal mechanisms available does not allow stress inversion of proper subsets partitioned according to earthquake magnitude. I look at a future availability of additional data in the study area to explore in a greater detail the local space variations of stress and the contributions by the different tectonic factors. On the other hand, I evidence the new contribution to knowledge of regional geodynamic processes given by the present study in comparison

Debora Presti
to the most recent stress inversion analyses carried out in the same area [Totaro et al., 2016]. In fact, although the primary seismogenic role attributed to plate convergence by Totaro et al. [2016] is confirmed here, the use of a more conservative stress inversion method allows me to better detect in the present study the heterogeneity of stress in the study area (see, e.g., Figure 5c). In particular, signatures of a rifting-related extensional component have, here, been detected into the compressional domain produced by plate convergence (Figure 8)  (approximate average error of focal mechanisms in the dataset). This indicates that plate convergence is likely the only geodynamic process active at these depth levels. A certain degree of stress heterogeneity is revealed by individual misfits at depth shallower than 30 km.

Conclusion
Earthquake relocations and stress inversion of focal mechanisms in the still poorly resolved geodynamic domain of western Ionian lead to the following conclusions: (i) the main finding of previous investigations indicating the primary tectonic action of Africa-Eurasia convergence in this part of the Mediterranean region [Montone et al., 2012, Montone andMariucci 2016;Totaro et al., 2016] is confirmed; (ii) the level and style of stress heterogeneity detected in the study area furnish, however, evidence of additional tectonic factors acting together with convergence.
In this regard, an additional factor can be recognized in the rifting process at the southwestern edge of the subducting slab recently hypothesized by Polonia et al. [2017]. SW-NE opening in the NW-trending belt comprised between the Alfeo-Etna and Ionian Fault Systems would add an extensional stress component to convergencerelated compression in the offshore of Eastern Sicily: this is really observed in the plot of stress orientations obtained by inversion of focal mechanisms (Figure 5c). Also, the hypocenter relocations evidence that the dislocation processes between the subducting slab (located northeast of the Ionian Fault System) and the adjacent lithosphere (southwest) are distributed over a relatively wide zone, probably because the subduction kinematics are slow and do not mimic fast STEP dynamics [Gallais et al., 2013;Orecchio et al., 2014]. This wide zone located between the Alfeo-Etna and Ionian Fault Systems (Figures 3d and 4d) corresponds to the NW-trending highly fractured zone of rifting where serpentinite diapirs rise from deeper depths [Polonia et al., 2017]. The stress inversion results of Figure   5c-d and the analysis of earthquake individual misfits in Figure 7 suggest that the seismic effects of rifting are confined to the upper 30 km of the area between the Alfeo-Etna and Ionian Fault Systems (a final sketch of the results into the geodynamic context is given in Figure 8c).
Data and sharing resources. Data used in the present study were collected from the databases of Istituto Nazionale di Geofisica e Vulcanologia (http://istituto.ingv.it/it/archivi-e-banche-dati) and from catalogs and bibliographic sources indicated in detail in the article.