https://orcid.org/0000-0001-6273-514X
SUMMARY
Extensional faults in Southern Calabria (Italy) have been widely studied for their capability of generating high magnitude earthquakes (Mw 7–7.2). An example is the historical seismic sequence occurred in 1783, which caused numerous fatalities near the villages located along the longest faults of this region: the Cittanova and the Serre faults. In this work, we estimated the seismic potential of these two faults by a kinematic block modelling approach using GNSS data of both campaign points and permanent stations. Our results indicate that both faults are accommodating the recognized extensional velocity gradient (∼ 1 mm yr−1) by long-term slip rates (∼ 2 mm yr−1). To estimate the back slip distribution and the interseismic coupling degree of the Cittanova and Serre faults, we discretized these by a triangular dislocation elements mesh. This approach has allowed us to distinguish the fault areas where elastic seismic rupture is more likely to happen from those affected by aseismic creeping behaviour. The obtained results show that the highest values of coupling are located near the shallow portion of the fault planes and near the southern tip of the Cittanova fault. We therefore estimated a set of possible rupture scenarios finding that the Southern Calabria domain is accumulating an interseismic moment rate at most equal to 2.16 × 1016 Nm yr−1, the equivalent of an earthquake of Mw 4.86 for each year.
1. INTRODUCTION
Analysis of ground deformation in active tectonic domains is possible using geodetic techniques since these are capable to estimate millimetre scale velocity variations (Segall & Davis 1997; Burgmann et al. 2000). Geodetic measurements, such as Global Navigation Satellite System (GNSS) data, are useful to identify structural lineaments (e.g. faults), define their kinematics and ability of accommodating the observed velocity gradient. This approach is adopted by several authors to explain deformation processes related to active faults observed at the interseismic timescale (e.g. Savage & Burford 1973; Meade & Hager 2005). In this work, we focus our attention on the Southern Calabrian Arc (Southern Italy, CA hereinafter), an active tectonic domain located in the central Mediterranean (Figs 1a and b). The CA is affected by large-scale tectonic processes, such as the relative convergence of the African and Eurasian plates. The CA belongs to a structural domain defined by the presence of a large accretionary wedge, located in the Ionian Sea, formed by the subduction of the Ionian oceanic crust beneath the CA (Malinverno & Ryan 1986; Jolivet & Faccenna 2000) (Fig. 1). In this complex tectonic framework, normal faults accommodate the extensional stress regime which is still active since the end of the Miocene (Westaway 1993; Monaco et al. 1996; Wortel & Spakman 2000; Goes et al. 2004; Faccenna et al. 2011). The regional crustal seismicity in the CA domain is mainly related to the activity of extensional faults, these are capable of generating M ∼ 7 events (e.g. the 1783 seismic sequence, Mw 6.8–7, the 1908 Strait of Messina earthquake, Mw ∼ 7.1) (Fig. 1b). Other remarkable seismic events are evidenced by historical and instrumental seismic records (Jacques et al. 2001; Barbano et al. 2005; Locati et al. 2022; Rovida et al. 2022; Andrenacci et al. 2023): The 1905 September 8 earthquake (Mw 7.5; Presti et al. 2017), the 1638 March 27 event (Bosi & Galli 2004) and the Crati earthquake, occurred in 1870 (Baratta 1901), (blue dots in Fig. 1a). In this paper we study the Cittanova and Serre faults (CF and SRF hereinafter), the longest normal faults slicing across the CA. Several authors (Jacques et al. 2001; Galli & Bosi 2002), agree in defining the CF and SRF as two West-dipping, ∼40 km long faults responsible of the 1783 seismic sequence. For instance, in Jacques et al (2001), the main shock of the 1783 seismic sequence is associated to the activity of the CF according to the historical reports of de Dolomieu (1784), who observed a 20 km long scarplet at the base of the Aspromonte Mountain, interpreted as a coseismic rupture along the CF. The CF has been the object of numerous studies such as paleoseismological trenches (Galli & Bosi 2002; Bosi & Galli 2004; Galli & Peronace 2015), which ascribed the 1783 February 5 event to the activity of this fault and defined return times for earthquakes similar to that occurred in 1783 (ranging from 1300 to 1800 yr). With regard to the other events of the same sequence, the February 7 and March 1 events, several authors ascribed these to the activity of the SRF. This hypothesis is supported by geological field observations, by the distribution of the epicentral areas of the 1783 sequence (Baratta 1901), by revisited macroseismic data (Andrenacci et al. 2023) and paleoseismological studies (Galli & Bosi 2002). Given the complex seismotectonic history of the CA, we decided to investigate the seismic potential of the CF and SRF. Indeed, even if these faults are widely studied, their role in determining possible future seismicity is still unclear, nevertheless the high seismic moment release against the low geodetic strain rate that characterize the CA (∼ 10–20 nanostrain yr−1) is still doubtful (see also Carafa et al. 2018 and reference therein).
(a) Seismotectonic framework of Calabrian Arc in the context of Europe-Africa relative convergence. (b) Location of the studied area represented by a dashed line, red line is the actual front of the collisional belt (Giuffrida et al. 2023 Modified, and reference therein). Dots represent the location of the main historical earthquakes occurred in the area, arrows represent the actual extensional dynamic. LCFS Lamezia Catanzaro Fault System; FC: Fagnano Castello Fault; SM: S. Marco Fili Fault; MR: Montalto Rende Fault; SRF: Serre Fault; CNFS: Coccorino-Nicotera Fault System; CF: Cittanova Fault; ARF: Armo Fault; RCF: Reggio Calabria Fault; W-F: West-Fault; TF: Taormina Fault; TFS: Tindari Fault System; AFS: Alfeo.Etna Fault System. Base map was made using DEM and Bathymetric data from: https://emodnet.ec.europa.eu/.
To study the seismic potential of the CF and SRF we investigated the interseismic deformation related to these by applying a kinematic Block Modelling approach (Mccaffrey 2002; Meade & Hager 2005; Meade & Loveless 2009), using campaign and continuous GNSS data acquired from 1995 to 2021 (e.g. Pirrotta et al. 2021; Giuffrida et al. 2023) (Fig. 2). We modelled the kinematics and the geodetic slip rates of the faults and investigated whether their presence and characteristics justify the recognized velocity gradient. We therefore estimated the back-slip distribution on fault planes (Savage 1983), and computed the interseismic coupling coefficient as the ratio of the back slip and the long-term slip derived from the blocks motions that border these faults. This step of analysis allowed us to distinguish the fault portions where the elastic contributions should occur from areas affected by creeping behaviour. Indeed, according to this approach, high coupling fault portions (coupling close to 1) are those that contribute elastically to reproduce the observed velocity gradients, whereas areas of low coupling (close to 0) creep. Areas of fault with medium coupling (∼0.5), are affected by complex behaviour (i.e. seismic ruptures and creep). The obtained results allowed us to propose a set of possible rupture scenarios by estimating the moment magnitude (Mw) and the recurrence interval (R) for each coupling interval on the SRF–CF fault plane.
Velocity field of the Southern-Central Calabria and north-eastern Sicily, computed considered the Hyblean Plateau fixed (Pirrotta et al. 2021 and Carnemolla 2021, PhD thesis); Dots indicate the IGM95 network; red triangles indicate the stations belongs to the RING/RDN network; yellow triangles indicate the TopNetLive stations and blue triangles indicate the permanent stations computed by the Nevada Geodetic Laboratory. Dashed lines are the section traces of the velocity profiles. Thick lines are the SRF and CF traces. Base map was made using DEM and Bathymetric data from: https://emodnet.ec.europa.eu/.
1.1 Seismotectonic framework
The Southern Calabria is an arc shaped terrain formed in the context of the relative convergence between the European and African plates (Malinverno & Ryan 1986; Faccenna et al. 2004; Rosembaum & Lister 2004) (Fig. 1). The CA is part of a large accretionary wedge formed as a consequence of the Miocene to Quaternary Subduction process of the Ionian crust (Deway et al. 1989; Faccenna et al. 2001). It is delimited to the North and to the South, by the Palinuro and the Tindari Fault Systems, respectively (Finetti and Del Ben 1986; Del Ben et al. 2008; de Guidi et al. 2013). The last tectonic phase of the CA is characterized by an WNW-ESE oriented extensional stress regime, which is still ongoing and has contributed to the formation of normal and oblique faults (Monaco & Tortorici 2000; Galli & Bosi 2002; Ferranti et al. 2008; Tripodi et al. 2018). The actual configuration of the Southern Calabrian Arc is defined by Upper Pliocene to Quaternary extensional basins arranged parallel to the main fault systems (Monaco et al. 1996; Jacques et al. 2001). The latters border the Aspromonte and the Serre Mountains, which are constituted by basement metamorphic rocks ascribed to the European margin (Ogniben 1969; Amodio Morelli et al. 1976).
As attested by several authors (e.g. Serpelloni et al. 2010; Palano et al. 2012), extensional faults with a prevalent dip-slip kinematic accommodate the ∼NW-SE stretching by generating high magnitude earthquakes (Mw ∼ 7). Examples are the CF and SRF, these almost 60° west-dipping faults have been regarded as the main seismic sources of catastrophic events, such as the ones occurred in 1783 (Jacques et al. 2001). The dimensions, strike and dip angle of the CF and SRF are well constrained by recent studies (see also Giuffrida et al. 2023 and references therein), widely documented through geological and geostructural observations (Tortorici et al. 1995; Monaco & Tortorici 2000; de Guidi et al. 2013; Scudero et al. 2020), paleoseismological studies (Galli & Bosi 2002; Galli & Peronace 2015) and seismological data (Scarfì et al. 2018, 2021). The CF and SRF delimit the extensional basins of Mesima and Gioia Tauro to the west from the Serre and Aspromonte Mountains to the east, respectively. These two faults exhibit clear morphological evidence such as cumulative scarps up to 400 m high, triangular and trapezoidal facets and v shaped valleys (Monaco & Tortorici 2000).
The complex seismotectonic framework of the Southern Calabrian Arc is characterized by the presence of other fault systems (Fig. 1). To the south, the 25 km long, NNE-SSW oriented, west-dipping Reggio Calabria Fault (RCF) was active during the Quaternary and various authors inferred its reactivation during the catastrophic 1908 earthquake (Monaco & Tortorici 2000). In the same area other faults have been recognized and considered capable of generating high magnitude events: The west-dipping, high angle Armo normal fault (ARF in Fig. 1; Aloisi et al. 2012) and the W-fault, a low angle (∼ 45°) east-dipping normal fault (W-F in Fig. 1), located in the Messina Strait. The latter was recently unlighted by high resolution reflection profiles (see Barreca et al. 2021 for further information). An active oblique belt, the SSW-dipping Coccorino-Nicotera Fault System (CNFS), separates the Gioia Tauro Basin from the Mesima Basin (Monaco & Tortorici 2007). To the North, the Lamezia Catanzaro Fault System (LCFS) separates the Southern Calabrian Arc from the Northern Calabrian Arc, bordering the Catanzaro basin. This fault system is composed by WNW-ESE oriented left stepping oblique faults with a prevalent normal kinematic and a subordinate left lateral component of motion (Dijk et al. 2000; Tansi et al. 2007; Pirrotta et al. 2021). Historical and instrumental seismicity of this sector testify the recent activity of the LCFS (Pirotta et al. 2021). In the Northern Calabrian Arc, the Crati Valley is bordered to the west by normal faults: the Fagnano-Castello (FC), the S. Marco Fili (SM) and the Montalto-Rende (MR) east-dipping faults (Fig. 1a), together forming an almost 45 km long, NNW-SSE oriented fault belt (Monaco & Tortorici 2000). These faults are considered responsible of the 1870 earthquake according to historical reports ( Baratta 1901) (Fig. 1a).
2. GEODETIC DATA
We computed a horizontal velocity field using a dense data set of about 50 GNSS stations belonging to various geodetic networks (Fig. 2). We occupied the IGM 95 (Istituto Geografico Militare) GPS points between July and September 2020, our survey consisted of at least 10 hr of occupation time for each geodetic benchmark, using Topcon Hiper SR and Hiper V receivers (L1/L2 double frequency). Moreover, we improved the robustness of the data set using the permanent stations of RING/RDN (https://webring.gm.ingv.it:44324/rinex/RING), the TopNet Live Italy network (https://rtk.topnetlive.com/italy/networks/topnet-live-italy) and the time-series of stations processed by the Nevada Geodetic Laboratory (https://geodesy.unr.edu/NGLStationPages/gpsnetmap/GPSNetMap_MAG.html). Permanent stations covered a large period, the older stations are those which belong to the RDN/RING networks (2005–2009), stations belonging to the TopNet Live networks were installed between 2012 and 2014 (see supplementary materials). The campaign points of the IGM-95 have time-series starting from 1994.
GNSS data were processed using GipsyX 1.5 (Bertiger et al. 2020), following the same methodology described in Pirrotta et al. (2021) and Giuffrida et al. (2023). We obtained the coordinates of the selected geodetic benchmarks for each DOY (Day Of the Year). These were used to compute time-series in ITRF 2014 reference frame (Altamimi et al. 2016). Subsequently, we referred the obtained velocities to the eastern Sicily block (Hyblean Plateau fixed) with the aim of visualizing the velocity field of the study area with respect to a favourable local reference system (see also Pirrotta et al. 2021 and Carnemolla 2021, PhD Thesis, for further information). The resulting velocity field and the related error ellipses (95 per cent level of confidence) are plotted in Fig. 2.
2.1 Velocity gradient across the Serre and Cittanova faults
Both campaign and permanent GNSS stations exhibit horizontal velocities (Fig. 2) of the same order of magnitude, ranging from a minimum value of about 1.5 mm yr−1 to a maximum value of about 4 mm yr−1. The obtained velocity field (referred with respect to the Hyblean Plateau fixed reference frame, Pirrotta et al. 2021 and Carnemolla 2021, PhD Thesis) shows a common motion towards the SE direction that is almost orthogonal to the traces of the main fault systems, that is, SRF (N40°E) and CF (N35°E), and exhibit a gradual increase moving away from West to East. Fig. 3 shows various sections orthogonal (numbered dashed lines in Fig. 2) to the traces of SRF and CF where is plotted the parallel component of the velocity rates of both permanent and campaign stations (blue and red dots, respectively, in Fig. 3), sections 1 to 3 are 10 km width, sections 4 to 6 are 15 km width. The velocity profiles allow us to visualize a predominantly extensional gradient. In particular, sections 1 to 3, orthogonal to the CF, exhibit a differential change in velocity from 2–2.5 mm yr−1 to almost 4 mm yr−1 moving towards the southeastern sectors. Similarly, sections 4 to 6, orthogonal to the SRF, depict a differential velocity change moving from NW to SE showing an increase in velocity from ∼ 2–2.5 to 3.5 mm yr−1 (see Fig. 3 caption for accurate minimum and maximum values of velocity for each section based on the weighted regression lines and the related R2 coefficients). The MNST station shows a slightly higher parallel velocity component (∼5 mm yr−1) than those of other nearby stations (e.g. ARSE and MOS3) indicating some possible errors and/or local processes not taken into account.
Parallel component of velocities across the CA. sections 1 to 3 are orthogonal to CF, whereas sections 4 to 6 are orthogonal to SRF. Section 1): min. vel. 2.10 ± 0.74 mm yr−1; max. vel. 3.48 ± 1.35 mm yr−1; R2 0.5. Section 2): min. vel. 2.42 ± 0.55 mm yr−1; max. vel. 3.11 ± 1.13 mm yr−1; R2 0.11. Section 3): min. vel. 2.68 ± 0.70 mm yr−1; max. vel. 3.0 ± 1.18 mm yr−1; R2 0.06. Section 4): min. vel. 2.51 ± 0.71 mm yr−1; max. vel. 3.99 ± 1.16 mm yr−1; R2 0.45. Section 5): min. vel. 2.52 ± 0.82 mm yr−1, max. vel. 3.88 ± 1.33 mm yr−1; R2 0.38. Section 6): min. vel. 2.18 ± 1.1 mm yr−1; max. vel. 3.5 ± 1.7 mm yr−1; R2 0.85. Blue dots are velocities obtained by using permanent stations (RING/RDN; TOpNetLive Italy), red dots are velocities of campaign stations (IGM-95), black dashed lines are the weighted regression lines for each section. Red vertical lines indicate the location of the CF and SRF for each section.
The velocity gradient recognized across the studied faults is also consistent with the orientation of the predominant extensional strain documented for these areas by structural analysis (de Guidi et al. 2013; Scudero et al. 2020), as well as with the seismotectonic framework of the Southern-Central Calabria (Tortorici et al. 1995; Monaco & Tortorici 2000; Ferranti et al. 2008; Serpelloni et al. 2010; Palano et al. 2012).
3. BLOCK MODELLING APPROACH
To reproduce the measured GNSS velocity field and estimate kinematically consistent long-term slip rates of the main faults in the study area, we applied a block modelling approach (Mccaffrey 2002) exploiting all the available geological and seismological information and using the BLOCKS MATLAB code (Meade & Loveless 2009). According to the Block modelling approach (McCaffrey 2002), the velocity gradient observed at surface by means of geodetic techniques is modelled considering the lithospheric layer divided into rigid blocks, bordered by rectangular fault planes embedded in an elastic, homogeneous and isotropic half-space (Okada 1992). At a more local scale of observations, blocks rotate with respect to a reference one considered fixed, allowing to study more in detail the measured velocity gradient. The block modelling theory follows the back-slip concept (Savage 1983), according to which the long-term blocks motions can be defined as the sum of the interseismic and coseismic deformations. Changing the order of the terms, it is possible model the interseismic deformations as the sum of the long-term steady velocity, related to the blocks motion and a long-term fault slip rate with a negative sign (i.e. acting in opposite directions with respect to the coseismic slip sense), called back-slip rate or slip deficit rate.
Accordingly, we can consider the interseismic velocity field (Vi) observed at the surface as the sum of two contributions (eq. 1): Vr, that represents the long-term rotational velocity of a point belonging to a rigid plate, and Ve, representing the elastic contributions to the velocity field due to the back slip rates of the modelled fault planes. Moreover, block model theory can consider a third parameter, Vs, as an additional deformation source for a possible internal strain of blocks. Following the above considerations, we can write Vi as:
One of the greatest strengths of the block modelling approach is the use of the rotation poles of the bounding blocks to compute the slip rate on faults, this ensures that the model is kinematically self-consistent (McCAffrey 2002). On the other hand, one limitation of this approach is that it must be supported by a good knowledge of fault geometry because the fault parameters (e.g. locking depth, dip-angle, strikes, faults length) are a priori fixed in the modelling.
3.1. Blocks configuration of the Calabrian Arc
Due to the complexity of the study area, we tried to evaluate the role of the numerous active faults within the occurring tectonic deformation. To explain the recognized geodetic velocity gradient and verify the role of the SRF and CF within the extensional regional framework, we tested three different model settings (Fig. 4), see also supplementary materials (Section 1; Tables S1, S2 and S3) for a summary of the input parameters used (i.e. blocks and faults geometries and locking depth). For each model, we constrained the position, extension, locking depth and dip angle of the faults from several seismotectonic studies and available data bases. As a first step, we considered the simplest model in which the velocity gradient is only accommodated by extensional faults, separating the two main blocks of the model (Model 1 hereafter, see Fig. 4). For both the SRF and CF we constrained the locking depth and dip angle (equal to 20 km and 60°, respectively) according to previous studies (Giuffrida et al. 2023). For the faults that border the Crati valley (Fig. 1), we considered a 10 km locking depth and 60° dip angle using the available information from the DISS data base (Data base of Individual Seismogenic Sources; https://diss.ingv.it/download-diss-3-3-0). To the south, we opted to consider the RCF as an active fault, due to its intermediate position with respect to the other faults slicing across this area (i.e. the ARF and the W-F). In this case we used a dip angle of 60° and a locking depth of 15 km, according to available studies (Ghisetti 1992; Monaco & Tortorici 2000). Moreover, we used two other fault segments taken from the DISS data base: the Piano Laga Savuto Fault, probable associated to the 1638 March 27 seismic event (Bosi & Galli 2004) (PLS in Fig. 1, located to the south of the Crati), with a locking depth of 10 km and dip angle of 60°, and the Taormina Fault, located in the eastern coast of Sicily, characterized by a locking depth of 15 km and dip angle of 60° (TF in Fig. 1).
Tested models for the given geodetic velocity gradient: Model 1 (2 blocks + 1 fixed), Model 2 (3 blocks + 1 fixed), Model 3 (4 blocks + 1 fixed).
The following model configurations included the oblique faults, the LCFS for Model 2, and the CNFS for Model 3 (Fig. 4). For both faults we considered high angle planes with 15 km locking depths according to available literature studies and DISS information. The more the number of blocks increases, the better the model reproduces the data, and it is important to evaluate whether the fit improvement is statistically significant. As the number of blocks was increased, we performed the F-test according to Stein & Gordon (1984) (Table 1) using the χ2 value (chi squared value of residuals), the number of stations used for the inversion process, the rotational parameters of the blocks involved (lat. long. and angular velocity of each block) and the number of blocks (n + 1, because 1 block is fixed). We find that the F values of the simplest model versus the complex ones are lower than the minimum F value of the statistical acceptance reported in Stein & Gordon (1984), with a 95 per cent level of confidence for the used data set (F 95 per cent conf. in Table 1). This means that our number of stations cannot solve for model block scenarios more complex than model 1 (2 blocks + 1 fixed) with an acceptable level of statistics.
Statistic F-Test, for each model were reported the computed residual chi2 (χ2), number of blocks and the F values of minimum statistical acceptance from Stein & Gordon (1984).
| N. blocks | χ2 | F | F (95 per cent conf.) | |
|---|---|---|---|---|
| MODEL 1 | 2 | 45.317 | Model 1→Model 2: 0.76 | 2.76–2.78 |
| MODEL 2 | 3 | 44.153 | Model 2→Model 3: 1.69 | 2.76–2.79 |
| MODEL 3 | 4 | 41.632 | – |
| N. blocks | χ2 | F | F (95 per cent conf.) | |
|---|---|---|---|---|
| MODEL 1 | 2 | 45.317 | Model 1→Model 2: 0.76 | 2.76–2.78 |
| MODEL 2 | 3 | 44.153 | Model 2→Model 3: 1.69 | 2.76–2.79 |
| MODEL 3 | 4 | 41.632 | – |
Statistic F-Test, for each model were reported the computed residual chi2 (χ2), number of blocks and the F values of minimum statistical acceptance from Stein & Gordon (1984).
| N. blocks | χ2 | F | F (95 per cent conf.) | |
|---|---|---|---|---|
| MODEL 1 | 2 | 45.317 | Model 1→Model 2: 0.76 | 2.76–2.78 |
| MODEL 2 | 3 | 44.153 | Model 2→Model 3: 1.69 | 2.76–2.79 |
| MODEL 3 | 4 | 41.632 | – |
| N. blocks | χ2 | F | F (95 per cent conf.) | |
|---|---|---|---|---|
| MODEL 1 | 2 | 45.317 | Model 1→Model 2: 0.76 | 2.76–2.78 |
| MODEL 2 | 3 | 44.153 | Model 2→Model 3: 1.69 | 2.76–2.79 |
| MODEL 3 | 4 | 41.632 | – |
Thus, we chose the simplest model with two blocks (Model 1) as the best to reproduce the observed geodetic velocity gradient. Model 1 allowed us to infer that the recognized extension of the region is mainly accommodated by long-term slip rate on the SRF and CF (uniformed distributed on fault planes), ranging from ∼1.3 to ∼ 2.5 mm yr−1, (Fig. 5). We also estimated slip rates for the RCF and TF, respectively ∼ 2.7 and ∼ 2.8 mm yr−1.
Long-term slip rate values of CF and SRF (red lines). Arrows are the residuals of velocities. The mean residual magnitude is 0.65 mm yr−1.
We conducted a sensitivity analysis to verify if the estimated slip rates are particularly dependent on locking depth variations. Accordingly, we ran the optimal block model setup (Model 1) by changing the locking depth of the SRF and CF from 5 km to 30 km with a constant step of 1 km. We also investigated the variability of the Chi2 value for each model, as well as the change of the mean residual magnitude of velocity, (see section S2 and Fig. S1 in supplementary material, for further information). The sensitivity tests indicate at most 0.4 mm yr−1 of variation of our estimated SRF and CF slip deficit rates, moreover, the minimum value of Chi2 is denoted for the model that assumes a 10 km of locking for the SRF and CF (see Fig. S1c in supplementary materials).
4. DISTRIBUTED BACK-SLIP OF SRF-CF SURFACE
To better evaluate the seismic potential of SRF and CF (which are located in the area with the highest data density), we computed the back-slip distribution on a 3-D fault surface starting from the best block model setting that had reproduced the recognized velocity gradient (Model 1, see section 3.1 and Fig. 4). In particular, the fault surfaces of SRF and CF have been divided into triangular patches (TDEs), using Gmsh Software (Geuzaine & Remacle 2009) and including the known faults attitude, dip direction and a maximum depth of 30 km (Fig. 6a). Thus, we made a new model run considering the same locking depth of Model 1 for all faults except for the SRF and CF, which we discretized by a triangular mesh on which we solved for back-slip. Estimating the back-slip distribution on fault planes is useful for discriminating fault portions where seismically coupled asperities and aseismic creeping areas are located (Meade & Loveless 2009; Anderlini et al. 2016). Coupled regions represent portions of the fault that entrain the elastic component of the velocity gradient recognized on the surface, whereas creeping portions accommodate the aseismic motion of the blocks in the long term (Sammis & Rice 2001). To reproduce this concept and following the approach of Anderlini et al. (2016), we constrained the back-slip values to be at most equal to the maximum values of the long-term slip rate obtained from the Model 1 by applying a linear constrained algorithm during the inversion (Coleman & Verma 2001). This inversion approach is useful to study areas where low interseismic velocity occurs (Anderlini et al. 2016), considering that for this kind of domain the original BLOCKS code (Meade & Loveless 2009) enables to obtain a good solution for the rotational block parameters, but not for the evaluation of spatial variation of the slip deficit on fault planes. Indeed, without any kinematic constraints implemented for the estimate of the slip deficit distribution, the high number of degrees of freedom given by the slip on each patch can allow the code to converge on non-kinematically consistent solutions. We constrain the back-slip rate distribution to taper to zero at the bottom of the fault planes in order to consider a possible brittle–ductile transition in depth (Scarfi et al. 2018, 2021). The back-slip rate distribution is regularized using a Laplacian operator (as in Meade & Loveless 2009), weighted by a smoothing factor (β), whose optimal value has been estimated by a trade-off reiterative approach (Harris & Segall 1987). This approach consists of a series of runs by varying the smoothing value in a range of values with constant steps: the more important the weight of the Laplacian operator that smooths the back-slip distribution solution is (i.e. low roughness of the solution), the more the residuals between data and model increase. At the end of the iterations, the trade-off plot compares the misfit of each solution with the slip roughness, and a ‘knee point’ is chosen in order to ensure the minimization of the residuals in terms of WRSS (Weighted Residual Sum of Square), and, at the same time, avoid an oversmoothed solution (see Fig. S2 in supplementary material). This approach allows to retrieve the main spatial features of the slip distribution, corresponding to the areal extent of possible fault asperities.
(a) Back slip rate distribution on fault plane (Dip slip component). Blue arrows represent the observed horizontal velocity field, red arrows represent the modelled horizontal velocity field resulting from the block modelling inversion. (b) Velocity sections AA’ and BB’ (30 km width) across SRF and CF respectively, showing the profile-parallel velocity components of both observed and modelled velocity field. Black lines indicate the polynomial best-fitting curves of the forwarded velocities. Relocated seismicity and faults trace from Giuffrida et al. (2023).
The back-slip rate pattern, for the modelled mesh, exhibits a distribution of negative values (corresponding to extensional rates), reaching zero moving down dip in the central portion of the surface, and maximum magnitude of about 3 mm yr−1. Variable values of back slip rates for the entire fault plane indicate that the surface is potentially coupled to some degree along the entire fault plane. Moreover, we noted a good fit (WRMS = 0.509 mm yr−1) between the observed horizontal velocity field and the modelled one (Fig. 6, blue and red arrows, respectively).
To compare the results of the modelled versus observed extensional velocity pattern and investigate the role of the SRF and CF to accommodate it, we built a regular grid of experimental points (EPs), with resolution of 1/60°, covering the entire southern Calabrian area. For each EP we computed the corresponding velocity derived from a direct computation of our block model solution. Our results are shown in Fig. 6(b) along two sections (AA’ and BB’) where we plotted the parallel component of the modelled velocities (red circles) obtained for each EP. For each section we also plotted the geometry of faults and the relocated seismicity of the area (see Giuffrida et al. 2023). Section AA’ shows a velocity gradient across the SRF increasing from west to east. The observed velocity ranges from 2–2.5 mm yr−1, in the hanging wall domain, to almost 3.5 mm yr−1 in the footwall sectors. In order to estimate the velocity gradient across the studied faults as accurately as possible, we fit a polynomial function to the modelled velocity field for each section (black lines in Fig. 6b). Then we computed the difference between the minimum and maximum values of the velocities and the associated propagated errors. The obtained parallel components of the GNSS modelled velocity suggest an extensional gradient of ∼ 0.9 ± 0.12 mm yr−1 along the overall section AA’ and a gradient of ∼ 0.39 ± 0.13 mm yr−1 through the SRF fault trace. Section BB’ shows the velocity pattern across the CF. Here, the stations near the footwall (MASP, MTSI, BOVA) suggest that the velocity increases from west to east from ∼2.5 to ∼3.5 mm yr−1. The forward modelled velocity pattern shows an overall extensional gradient of ∼0. 62 ± 0. 15 mm yr−1 and of ∼0. 31 ± 0. 18 mm yr−1 through the fault.
4.1 Interseismic coupling of SRF-CF
The slip deficit rates resulting from block modelling are useful to determine the interseismic coupling degree (IC). This parameter is helpful to assess the spatial diffusion of coupled asperities that are potentially capable of generating earthquakes (Mccaffrey 2002; Meade & Loveless 2009; Graham et al. 2021). The IC is defined as the ratio of the elastic back-slip distribution and the interseismic long-term slip rates derived from the long-term motion of the blocks (see also Ader et al. 2012). The IC varies from 0, where no back slip on the fault should occur (areas of fault with no coupling), to 1 where back slip occurs at the same rate of the long-term block motions (areas of fault with full coupling). These high coupling areas represent the spatial extent of asperities that may generate earthquakes. According to our model, the IC distribution (Fig. 7a) suggests the presence of two possible coupled asperities along the fault surface (IC > 0.5). The first asperity is almost 35–40 km long and is located in the shallowest portion of the fault extending down to 10 km of depth. The second asperity is located at the southern termination of the fault, covering almost the entire depth (IC close to 1 from 0 km to ∼ 20 km and IC > 0.5 from ∼ 20 km to ∼ 30 km). Areas of low coupling (IC < 0.5) are detected in the deepest portions of the surface and the shallowest sector of the northern tip, close to which most of the recorded microseismicity is located (Figs 7a and b), that is probably ascribed to creeping behaviour according to Rubin et al. (1999). Areas of intermediate coupling (IC close to 0.5) may suggest complex fault portions where the recognized extensional gradient can be accommodated by both elastic and creeping processes. To compare the IC distribution with the microseismicity, we selected an overall ∼ 1000 seismic events, falling within the SRF-CF plane, from the original data set of ∼ 8500 relocated events in Giuffrida et al. (2023). The selected earthquakes have a magnitude between 1 and 4.1 and covering a period from 1984 to 2021. Fig. 7(c) shows the frequency of the coupling value of the closest TDE associated with each seismic event, accordingly, a high frequency is visible for coupling values of 0.4 and 1. Normally, the microseismicity should fit the areas of low coupling, this is in line with the high frequency of events at 0.4 of IC. The high frequency of events close to IC 1 can be explained by the overlap of the microseismicity with the coupled areas located in the southern tip of the CF (Figs 7a and b), this could be indicative of an IC pattern that is not well detected in this fault region. Fig. 7(d) shows the frequency of distances of each seismic event to the closest TDE of the model surface with an average value of ∼ 5 km.
(a) 3-D view of the triangular patches fault plane (dipping towards the West sectors) coloured as function of the interseismic coupling distribution, black circles represent the relocated seismicity across the fault plane; black lines are the coupling contour. (b) Map view of the CF-SRF surface and the relocated seismicity belongs to the fault area, black lines are the coupling contour. (c) Histogram shows the frequency of the coupling values associated with the seismic events closest to the modelled surface. (d) Histogram shows the frequency of distances related to the seismic events closest to the modelled surface. (e) Resolution length distribution expressed in km. Black lines are the resolution contours. (f) Map view of the resolution contour lines and distribution of the stations used for the modeling approach. Red triangles are GNSS of RING network, yellow and blue are stations of TopNetLive and Nevada Geodetic Laboratory respectively, black circles are IGM95 campaign stations. Figure made with MOVE software (Petex, academic grant).
To verify the reliability of the obtained IC distribution and test how well the spatial location of our data can resolve the spatial extent of the recognized coupled asperities, we performed a resolution length analysis (Ader et al. 2012; Anderlini et al. 2016) (Fig. 7e). An asperity can be considered well resolved if the resolution length has the same order of magnitude or is smaller than the asperity size. In Fig. 7(e) is possible visualize the obtained resolution length on the fault plane in km, whereas Fig. 7(f) shows the map view of the resolution contour lines and the distribution of the stations used. The lowest resolution is recognized in the deepest portion of the surface, where it reaches the value of 50 km. On the contrary, the highest resolution is detected in the shallowest portion of the surface (up to 15 km of depth) and ranges from 6 km to almost 20 km. Thus, according to the resolution analysis, the deepest part of the surface, where low coupling occurs, corresponds to a low resolution, as in this sector our data cannot resolve for small asperities. Conversely, the 35–40 km long asperity, recognized in the shallowest portion of the fault, seems to be well resolved by the data since the resolution length is smaller than its size. This coupled portion is well resolved also in its lateral portions, as here we detected low coupling values accompanied by low resolution length. Finally, if we consider the microseismicity distribution, it seems to be clustered below the asperity, where low coupling values occur.
A more complex scenario can be considered for the asperity located at southern tip of the fault. In this sector the resolution length reaches 30 km, this value is higher than the size of the detected coupled portion, accordingly, the IC distribution is not well resolved by the data, especially in the deepest and lateral edges of the asperity. We can also consider the distribution of the microseismicity: it seems to overlap with the coupled area and it would not justify the elastic behaviour of this fault portion, even though the horizontal velocity gradient is mostly accommodated by the high back-slip rate. This also explains the high frequency of events close to IC 1, detected in Fig. 7(c).
Moreover, we conducted further resolution tests to investigate the resolvable patches size using a checkerboard approach. These tests are useful to evaluate the reliability of the obtained coupling pattern and the ability of the GNSS data to recover the spatial distribution of the asperities on fault plane (see Section 4 and Figs S3 and S4 in supplementary materials). The checkerboard test provided further validation of the good resolution for the shallow part of the fault.
5. DISCUSSION
We used a dense data set of GNSS stations belonging to permanent and discrete networks to estimate the interseismic velocity field of the southern central Calabria and verify the role of the SRF and CF faults within their geodynamic setting. We found a prevalent ∼WNW-ESE oriented extension with low velocity gradient orthogonal with respect to the SRF and CF. The low velocity gradient and strain rate associated with the high seismic moment release of this area, point out several unresolved questions on the seismic assessment of the Calabrian Arc and testify that this domain represents a challenging region to study (Carafa et al. 2018). The recognized velocity gradient across the SRF and CF is well expressed by the forward model of the velocities (Fig. 6b), suggesting that the activity of these faults is probably accommodating most of the extension. However, we must also consider the presence of other faults, such as the CNFS and the LCFS, that may contribute to reproducing the observed velocity pattern. Additionally, the complex dynamics affecting the Messina Strait, constituted by a long extensional faults belt (see also Barreca et al. 2021) should be considered.
To better understand the seismic potential of the study area, we applied a self-consistent kinematic block modelling approach by testing different block settings and tectonic scenarios where normal and oblique faults were considered to variably accommodate the observed geodetic velocity gradient. Then, we performed a statistical F-Test analysis to evaluate the model that best reproduces the geodetic velocity gradient with the minimum necessary number of degrees of freedom, finding that the simplest model considering the activity of the SRF and CF faults is the best candidates, while estimating the corresponding long-term faults slip rates (assuming uniform slip on fault planes). Then we modelled the SRF and CF as a single 3-D west-dipping 30 km deep fault surface subdivided into TDEs, considering the geometric parameters provided by previous studies (Giuffrida et al. 2023 and reference therein), to compute the slip-rate deficit distribution on fault. A reliable locking depth estimate for these faults is challenging because of the absence of onshore seismic reflection profiles in these areas. For this reason, the modelled 3-D mesh extended in depth up to 30 km, that allows us to consider the overall crustal layer for which we solved for slip deficit, according to the known depth of the Moho from tomographic Vp model of Scarfì et al. (2018). Moreover, the above step of analysis was helpful to investigate the crustal portion that is likely affected by a prevalent elastic behaviour, neglecting the viscous-elastic layer which is not contemplated in the block modelling approach.
The obtained coupling degree distribution shows two main asperities. The first one is located in the shallowest portion of the plane, within the first 10 km of depth, and has a total length of about 35–40 km. The resolution analysis allowed us to infer that this asperity is well resolved, both laterally and in depth. The highest resolution is consistent with the location of the closest GNSS stations (i.e. MASP, CELL, OPMA, MTLM SSBR and FILA), roughly covering the central portion of the fault plane. We can do the same consideration for the uncoupled portions of the fault where high resolution occurs, for instance in the lateral portion of the asperity previously described and below it; this is consistent with the presence of stations PZCL, VIBO, FILA, SPAR and GTRO.
This coupled shallow fault area has been compared with the available information on the subcrustal structure of the study area (Fig. 8) as well as with the distribution of the relocated microseismicity from Giuffrida et al. (2023). In Fig. 8 the red dashed line represents the top of a low velocity anomaly layer, well defined through the Vp model of Scarfì et al. (2018). The latter is interpreted as a 15-km-thick weak material made up by the sediments accreted in the internal portion of the Calabrian wedge. This could be an interesting feature if compared with the estimated IC distribution, that seems consistent with the presence of creeping portion where the weak material is located, as well as with the distribution of the microseismicity.
3-D view of the SRF-CF surface coloured as function of coupling degree (blue: coupled; light red: uncoupled), the black dotted line represents the top of the Moho (from Scarfì et al. 2018) that intersects the modelled 3-D mesh. The red dashed line represents the top of the low Vp anomaly layer (Scarfì et al. 2018), its elevation spanning from −17 km in the northern sector to −15 km in the southern sector. Figure made with MOVE software (Petex, academic grant).
The second asperity concerns the southern tip of the fault; it is 10–15 km long arranged for almost the total depth of the surface. From the resolution analysis we found that this asperity is not well resolved, especially for its deepest portion, where a high distribution of microseismicity also occurs. In this area, the top of the Moho discontinuity (black dotted line in Fig. 8) from the Vp velocity model of Scarfi et al. (2018), (considering a Vp velocity of 7.6–7.8 km s−1 according to Rabbel et al. 2013), rises up to 22 km of depth, interacting with the modelled 3-D mesh and suggesting a possible visco-elastic behaviour of this region. Thus, we suggest that this asperity can be reliable only for its shallowest portions (i.e. 5–10 km of depth), where high resolution occurs, consistent with the presence of the stations MASP and CELL.
5.1 Rupture scenarios
To investigate possible consequences of the proposed model in terms of seismic potential of the Southern Calabria area, we conducted a detailed study based on both evaluation of the area of the coupled asperity and the corresponding moment magnitude (Mw). We implemented the code used in Graham et al. (2021), as an extension of the original BLOCKS code, to evaluate the interseismic moment accumulation rate of the detected central asperity for each coupling increment, and the time needed to release this energy based on the geodetically constrained slip deficit rate derived from the block modelling. According to Graham et al. (2021), we used known empirical law scaling (in this case Wells & Coppersmith 1994, table 2a, all slip type) to estimate the corresponding Mw (and the related uncertainties ranging from ± 0.12 to ± 0.14), furthermore, we estimated the recurrence time interval using Aki (1972) inverse formulation (R = M0/μAs), where s is the slip deficit rate of each coupled triangles within the mesh, μ is the shear modulus (assumed to be 30 GPa), A is the area of the considered fault surface and M0 is the seismic moment previously computed with Hanks and Kanamori relation (Hanks & Kanamori 1979), (Fig. 9).
(a) Rupture scenarios in Southern Calabria based on the proposed interseismic coupling for the SRF-CF faults surface. Each subplot in figure indicates a different coupling fraction (from 0.1 to 0.9) and the moment magnitude (Mw) given by clusters of triangular elements coupled at or above this coupling increment. Each cluster is coloured based on Mw; R indicates the corresponding recurrence interval (years) for each coupled cluster. (b) plot describing the relation between Mw, R and coupled areas of fault (in km2), for each coupling increment.
Following the described approach, a possible rupture scenario is estimated for each coupling increment. In particular, coupling ranges from ≥ 0.1 (which include triangles with almost free slip) and coupling ≥ 0.9 (including triangles with slip deficit rates almost equal to the long-term block motion), with a step of 0.1. The area of the coupled triangles decreases as the coupling threshold increases. The estimated Mw follows this trend as does the computed recurrence interval (R). In general, rupture scenarios including high values of coupling can be considered more hazardous than the weak ones because they are accompanied by short recurrence time intervals with still high Mw. Furthermore, using the 0.5 coupling threshold scenario to estimate an average interseismic moment accumulation rate of the studied area, we achieved a value of 1.62×1016 N m yr−1. Whereas considering the 0.1 coupling fraction scenario, we achieved a maximum value of 2.28×1016 N m yr−1. We also compared the obtained interseismic moment accumulation rate, derived from our slip deficit calculation, with the average coseismic moment release rate derived from the relocated seismicity in Giuffrida et al. (2023), for the first 25 km of depth in the period 1990–2018. We found that the southern Calabria is releasing 1.19×1015 N m yr−1 of seismic moment. By calculating the difference between the interseismic moment accumulation rate and the coseismic moment release rate, we suggest that the southern Calabria is at most accumulating energy ranging from an equivalent Mw 4.86 earthquake (if we would consider the 0.1 scenario), to an equivalent Mw 4.75 (if we would consider the average 0.5 scenario), for each year. Moreover, we should also consider that the total of the energy can be released not only by elastic contribution but also through aseismic behaviours, thus our scenarios represent an end-member scenario, where the total amount of energy is released elastically. Finally, we used the maximum area of the coupled triangles (see rupture scenario for 0.1 of coupling threshold) to calculate the time needed to accumulate energy capable to generate an earthquake of Magnitude 7.2, similar to that happened during the 1783 seismic sequence as attested by several studies (i.e. Jacques et al. 2001 and reference therein). Considering our slip deficit rate calculation, we found that the time required to accumulate such energy (corresponding to a M0 ∼2.76 ×1019 N m) is equal to 1 300 ± 230 yr, consistent to the recurrence interval estimated in Galli & Bosi (2002), (1.3 to 1.8 kyr) by paleoseismological trenches.
6. CONCLUSION
Using both permanent and campaign GNSS stations, we estimated the velocity field of the southern-central Calabrian domain finding a general trend of motion towards the SE Ionian sector with respect to the Hyblean Plateau (Pirrotta et al. 2021 and Carnemolla 2021, PhD Thesis). The velocity rate is between 3 and 4 mm yr−1, which is consistent with previous studies (Palano et al. 2012; Giuffrida et al. 2023 and references therein). Moreover, the obtained velocities show a general extension across the Calabrian Arc, arranged orthogonal with respect to the main fault systems.
By applying a kinematic block modelling approach, we found that the recognized extensional pattern is mainly accommodated by long-term slips rate on the SRF and CF, two large West-dipping normal faults bounding the extensional basins of Mesima and Gioia Tauro Plains, respectively. These faults are responsible for several historical earthquakes, for instance the 1783 seismic sequence; in particular, the Mw ∼7 main shock (which occurred on 5 February) is one of the most destructive earthquakes to have hit the Italian Peninsula.
Thus, to investigate the seismic potential of these faults, we modelled the fault surface by TDEs; then we computed the back-slip rate for each TDE starting from the long-term slip rate and blocks velocities previously obtained. Finally, to detect potential asperities capable of generating earthquakes and areas affected by aseismic behaviour, we estimated the interseismic coupling degree as the ratio of the back-slip rate variability and the long-term slip rate.
Our results show the presence of at least one asperity for the first 10 km of depth, located in the central portion of the modelled surface and extending for almost 40 km in length. It is well resolved by the spatial data distribution since the resolution length analysis exhibits low values for the corresponding high coupled area. A second asperity was recognized at the southern tip of the modelled fault plane. Considering the resolution analysis, we suppose that this asperity is not well resolved in its deepest portion, but it is reliable only for its shallowest parts. Moreover, the low coupled area of the plane apparently fit with the microseismicity distribution, but we cannot rule out aseismic creeping processes in these areas, since the data resolution seems to be insufficient to determine this, especially regarding to the southern asperity.
Finally, using Wells & Coppersmith (1994) empirical relation, we developed a set of potential rupture scenarios and estimated the corresponding recurrence time intervals related to the fault area for which the IC is greater than a specific values, varying from 0.1 to 0.9. Then we inferred that the southern Calabria is actively accumulating a total seismic moment rate between 1.5 ×1016 N m yr−1 (considering IC ≥ 0.5) to 2.16 ×1016 N m yr−1 (considering IC ≥ 0.1), respectively equivalent to an earthquake of 4.7–4.8 Mw each year. We also estimated the time needed to accumulate the energy capable of generating an earthquake of similar magnitude to the one that occurred in February 1783 (Mw ∼ 7), which is about 1 300 ± 230 yr according to our slip deficit calculation, consistent to previous estimates of the recurrence interval for the studied faults.
ACKNOWLEDGMENTS
This research is part of SG’ s doctoral project at the University of Catania and was funded by the MUSE 4D project—Overtime tectonic, dynamic and rheologic control on destructive multiple seismic events—Special Italian Faults & Earthquakes: from real 4D cases to models in the frame of PRIN 201, grant number ‘2017KT2MKE’; and has been supported by the PIAno di inCEntivi per la RIcerca di Ateneo (PIACERI 2020/2022) (Fund manager Giorgio De Guidi). This research was also funded by the AdriaLab project- On the dynamics of the Adria plate: a natural laboratory for multiple subduction systems in the frame of PRIN 2022, under grant number 2022XZ3W22 (fund manager Carmelo Monaco).
AUTHOR CONTRIBUTION
S. Giuffrida (Conceptualization [equal], Investigation [equal], Methodology [equal], Writing – original draft [equal]), L. Anderlini (Conceptualization [equal], Methodology [Equal], Software [equal], Writing – review & editing [equal]), F. Carnemolla (Data curation [equal], Writing – review & editing [equal]), F. Brighenti (Visualization, Writing – review & editing [equal]), G. De Guidi (Conceptualization [equal], Supervision [equal], Writing – review & editing [equal]), F. Cannavò (Visualization [equal], Writing – review & editing [equal]), S. Graham (Conceptualization [equal], Software [equal], Writing – review & editing [equal]), and C. Monaco (Supervision [equal], Writing – review & editing [equal])
CONFLICT OF INTEREST
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The authors also acknowledge the use of MOVE Software Suite granted by Petroleum Experts Limited (www.petex.com).
DATA AVAILABILITY
Data used to achieve the presented results are available on the following links: Permanent stations from NGL (Nevada Geodetic Laboratory, https://geodesy.unr.edu/NGLStationPages/gpsnetmap/GPSNetMap_MAG.html), RING Network ( https://webring.gm.ingv.it:44324/rinex/RING) and TopNET live Italy Network (https://rtk.topnetlive.com/italy/networks/topnet-live-italy). Campaigns points belong to IGM Network (instituto Geografico Militare, https://www.igmi.org/it/Home). GNSS data solutions are processed using GipsyX licensed under JPL (Jet Propulsion Laboratory). GypsyX is available at: https://gipsyx.jpl.nasa.gov/, (Bertiger et al. 2020). The software used in this paper for blocks modeling (BLOCKS code) works on MathLab environment (https://it.mathworks.com/help/install/ug/install-products-with-internet-connection.html, MathWorks Inc. license, academic grant) and is licensed under MIT (Massachusetts Institute of Technology) and published on GitHub https://github.com/jploveless/Blocks (Meade & Loveless. 2009). Images are made by MOVE Software Suite, granted by Petroleum Experts Limited, (https://www.petex.com/pe-geology/move-suite/, academic grant) and ArcMap granted by Esri (https://www.esri.com/en-us/arcgis/products/arcgis-desktop/resources, academic grant).
