The Northern Apennines fold-and-thrust belt have been developing since the Late Cretaceous following the closure of the Ligurian-Piedmontese ocean, intervening between the Europe and Adria plates (Marroni et al., 2017, and references therein). Currently, the Northern Apennines chain is composed of a complex nappe pile in which the uppermost tectonic units are the ocean-derived Ligurian Units (Jurassic-Eocene), overthrusting the Tuscan and Romagna Units (Triassic-Miocene) deposited on the thinned continental passive margin of the Adria Plate (Bonini et al., 2014, and references therein; Figure 1A).

The polyphase compressive tectonics which gave rise to the current fold and thrust belt continued until recent times, as suggested by recent out-of-sequence thrusting and reactivation of older thrusts (Boccaletti and Sani, 1998; Finetti et al., 2005; Bonini et al., 2014), and is still ongoing in the external sector of the chain (i.e., Po Plain) as the compressive seismic activity demonstrates (e.g., Maestrelli et al., 2018, and references therein). However, along the main Apennine divide the stress field is currently extensional, with a roughly NE-SW-oriented minimum stress axis (σ₃), as suggested by focal mechanisms (Pondrelli et al., 2002; Chiarabba et al., 2005; Pondrelli et al., 2006; Chiarabba et al., 2015) and by well developed normal fault systems, superimposed on previous compressive structures (Sani et al., 2009; Bonini et al., 2016; Sani et al., 2016, and reference therein).

The belt of the outermost hinterland intermontane Tuscan basins (Mugello, Casentino, Val Tiberina, Figure 1A), located near the main divide, is crucial for assessing the active tectonics and the related seismic hazard of the Northern Apennines. The Mugello is one of the more active areas of this belt and its seismicity shows a high variability in terms of intensity, depth, and kinematics (Bonini et al., 2016; Figure 1B). The Mugello Basin is inferred to have developed in Early Pleistocene under a compressive regime and was later affected by normal faults as the compressive regime ceased around the Early-Middle Pleistocene transition (Sani et al., 2009).

The study area is characterized by a medium intensity seismicity with historical earthquakes with highest macroseismic magnitude of Mw = 6.02 (1542) and Mw = 6.38 (1919) (Rovida et al., 2020; Rovida et al., 2021). A large SSW-dipping normal fault system, referred to as “Ronta fault system,” (RFS hereinafter) affects the pre-basin substratum mostly made of Miocene sandstones, and delimits the northeastern basin margin (Sani et al., 2009; Figure 2). This fault system is ∼25–30 km long and displays remarkable morphostructural evidence (Sani et al., 2009; Bonini et al., 2016). Fault-slip data indicate mainly dip-slip normal kinematics (average rake = −85) and an average slip-rate of 0.16–0.37 mm/year has been estimated (Sani et al., 2009; Bonini et al., 2016).

The southwestern margin of the basin is affected by a system of NE-dipping antithetic normal faults, whose morphologic expression is not well pronounced (Benvenuti and Papini, 1997; Sani et al., 2009). However, it is believed to control the asymmetric location of the Sieve River on this side of the basin (Benvenuti and Papini, 1997), and previous studies locate the basin master fault along this margin (Martini and Sagri, 1993), which would be connected to the regional NE dipping low-angle Etrurian normal fault system (Boncio et al., 2000). The DISS Working Group (2021) catalog follows this interpretation and subdivides the fault system into a north-western and south-eastern segment, interpreted to have caused the earthquakes of 1542 and 1919, respectively (Basili et al., 2008). Transverse faults at the north-western and south-eastern basin margins might represent other potentially active structures (Delle Donne, 2005).

During the seismic event of 1919, coseismic east–west-trending fractures affecting the ground surface were observed in the central-eastern part of the Mugello Basin (Sani et al., 2009). Moreover, liquefaction phenomena are reported nearby the Sieve River (Galli and Meloni, 1993). Unfortunately, the rough location of these deformation features does not allow identifying the source fault of that earthquake. The distribution of the sites with highest damage may suggest that a transverse fault may have also been activated by the 1919 event (Delle Donne, 2005; Sani et al., 2009). Based on length, lateral continuity, and its clear morphological evidence, the RFS could be considered the most likely source of the biggest earthquakes that hit the Mugello area, including the 1919 event (Bonini et al., 2016).

Following that destructive earthquake, the area has been affected by weak seismicity, except for short-duration seismicity bursts characterized by a typical mainshock-aftershock sequence, and maximum magnitudes between 4 and 4.5. The most recent of these sequences occurred in 2008 (Mw = 4.5), 2009 (Mw = 4.2), 2015 (Mw = 4.3). This latter sequence, which occurred to the NW tip of the area of study and externally to the basin, was then followed by a minor sequence in 2017 located further NW of the 2015 cluster. The most recent sequence occurred in 2019 (Mw = 4.5), at the NW margin of the basin. In our research all these seismic sequences have been relocated, allowing a better definition of orientation and dip of the causative fault segments. Particular attention is dedicated to the 2019 sequence, which was also recorded by the dense, temporary, local seismic network.

Phase arrival times from the integrated catalog were inverted using NonLinLoc (Lomax et al., 2009), a non-linear, probabilistic procedure based on the original work of Tarantola and Valette, (1981). Theoretical travel times were calculated in reciprocal geometry using a finite-difference scheme, applied to the minimum-misfit 1D seismic velocity model derived by Piccinini et al. (2014) for the same area from inversion of local earthquakes. The resulting catalog amounts to 11,234 locations, obtained by inverting at least six wave arrival times, with a minimum of 2 S-wave phases. Within the study area, most hypocenters are distributed according to four main spatial and temporal clusters, which occurred in March 2008 September 2009, January 2015 and December 2019 (Figure 3). For all these sequences, the mainshocks had similar magnitudes, spanning the range Mw = (4.2–4.5).

Earthquakes are then relocated using HypoDD, a computer code which implements the double-difference (DD) algorithm of Waldhauser and Ellsworth (2000). The DD technique exploits the similarity of ray paths between the source and a common station when the hypocentral separation between two earthquakes is small with respect to the event-station distance and the wavelength of velocity heterogeneity. This approach permits to reduce systematic errors related to the prediction of travel times in poorly known velocity structures, and hence to achieve more focused images of hypocenter clustering. We only consider earthquakes included in a box of corners [43.95 N, 11.05 E] [44.30 N, 11.80 E], and located with a minimum of 12 phases.

Source mechanisms for the most relevant earthquakes are taken from the Time Domain Moment Tensor (TDMT) catalog (Scognamiglio et al., 2006), which reports the solutions obtained using the full-waveform moment-tensor inversion of Dreger and Helmberger (1993), as described in Scognamiglio et al. (2009).

Stress relieved during an earthquake can influence the stress state on faults surrounding the ruptured fault (e.g., Stein, 1999; Kilb et al., 2002; Kilb, 2003; Hill and Prejean, 2007; King and Devès, 2015). We thus evaluate the potential role of dynamic and static stress changes on the sequence of small to moderate earthquakes that struck the northern margin of the Mugello basin in 2008–2019.

Static stresses can permanently modify the stress in the crust around the fault rupture, and attenuate away from it as 1/R ³, R being the epicentral distance (Aki and Richards, 1980). Consequently, stress acting on faults nearby the earthquake can be changed depending on the position, geometry and rake of the considered receiver fault. Coulomb stress changes imparted by a mainshock were computed in a homogeneous elastic half-space and resolved on the fault rupture associated with the subsequent mainshock using available source fault models (Supplementary Tables S1, S2) and the Coulomb 3.3 software (Lin and Stein, 2004; Toda et al., 2005). Coulomb stress changes (ΔCFF) are defined as ΔCFF = Δτ + μ′ Δσ_n, where Δτ is the shear stress change (positive in the direction of slippage), Δσ_n is the fault-normal stress change (positive for fault unclamping), and μ′ is the apparent friction, which is an approximation of fault properties (e.g., static friction, pore-pressure conditions). Failure is encouraged if ΔCFF is positive, whereby a fault will fail if the applied stress increments exceed a stress threshold specific for the fault. Increased shear stress and fault unclamping are normally considered favorable conditions for failure. Fault failure is instead hindered if ΔCFF is negative, which occurs for decreased shear stress on faults and if the fault is clamped.

Dynamic stresses, σ _d , produced by passing seismic waves decay much slower with epicentral distance (as ∼1/R ² or ∼1/R ^1.5; Hill and Prejean, 2007), and can thus trigger earthquakes distantly (Hill et al., 1993; Kilb et al., 2002; Brodsky et al., 2003). Most observations indicate that dynamic triggering occurs in concomitance of the transit of the largest-amplitude seismic arrivals, such as the Rayleigh waves (e.g., Hill and Prejean, 2007). However, other studies suggest that such transient stress changes may alter the properties of a pre-existing fault, accelerating its evolution toward failure (Parsons, 2005). In that case, the triggering may be delayed by up to several years (e.g., Kilb, 2003).

The effects of the dynamic stresses are evaluated by estimating, for each mainshock, the peak ground velocity (PGV) at the faults which ruptured subsequently. PGV is derived from the ground motion predictive equation of Bindi et al. (2011), which is routinely used for the production of shake maps (Michelini et al., 2020; Supplementary Figure S3). This analysis admittedly involves some degree of uncertainty related to the simplifying approximation that PGV predicted at the earth surface is similar to that at hypocentral depth. Nevertheless, the decrease with depth of amplitude and strain of surface seismic waves will be balanced at some extent by the tendency of elastic moduli to increase with depth (Hill and Prejean, 2007). Peak dynamic stresses (σ _pd ) are derived from PGV following the relationship (Hill et al., 1993; Hill and Prejean, 2007; Fan et al., 2021): σ p d = P G V   G β − 1 (1) where G is the shear modulus and β is the average S-wave velocity, here taken equal to 32 GPa and 3,200 m/s, respectively.

The high-resolution relocation of seismic data clearly indicates that the four main sequences that affected the study area during the past 13 years were associated with the subsequent activation of adjacent segments of a larger fault system bordering the NE sector of the Mugello Basin and adjoining areas. In particular, the 2019 and 2009 sequences are clearly correlated with two segments of the RFS. Both focal mechanisms and the SW-dipping fault planes imaged by the hypocenter clusters correlate remarkably well with the fault geometry delineated by the surface geology (Sani et al., 2009; Bonini et al., 2016).

The 2008 and 2015–2017 sequences are located externally to the basin, but in continuity with the RFS. The 2008 sequence is related to a NE-dipping plane, in agreement also with what reported by Piccinini et al. (2014), and Bonini et al. (2016). The fault segment activated during the 2015 Castiglion de’ Pepoli sequence is steeply dipping toward NE, and it is likely related to a system of normal faults affecting the Castiglion de’ Pepoli anticline (Figures 2, 5). The geometrical relationships between the RFS and the NE-dipping fault planes are unclear, even though the differences in fault geometry and morphological characteristics (with the RFS being characterized by marked morphological features, and the NE-dipping fault planes being mostly blind) may suggest that they belong to distinct fault systems.

The activation of the 2009 and 2019 fault ruptures fits well with the overall distribution of historical earthquakes, which are arranged along a NW-SE-trending belt, which runs sub-parallel to the RFS, near the North-Eastern margin of the basin.

These findings have obvious implications for the regional seismotectonic setting of the study region. In particular, both geometry and depth extent of the ruptured faults are not consistent with the activity of a low-angle, NE-dipping normal fault (i.e., the Etrurian Fault System; Boncio et al., 2000; DISS Working Group, 2021) which is hypothesized to be emerging at the South-Western margin of the basin. As a matter of fact, simple geometrical considerations suggest that the rupturing depths (which are of the order of 10 km) are deeper than the prosecution at depth of such a hypothetical low-angle normal fault. Finally, this observation also rules out the possibility that the SW-dipping normal faults represent antithetic structures to the Etrurian Fault System. These findings thus challenge the ongoing activity, or even the existence of the Etrurian Fault System in the Mugello region.

The four main shocks occurred at relatively close distances and during a relatively short time span, and therefore a role of seismic triggering could potentially be involved. Nevertheless, the triggering mechanisms should explain the relatively long time intervals between the mainshocks, which occurred with delays varying between 1.5 and 10 years. Actually, identifying the processes that control delayed activation is intricate: on the one hand, permanent static stresses can be predicted as a possible mechanism, but they are small. On the other hand, dynamic stresses are relatively large, but transient, with the triggering more than a year after the passage of the seismic waves being difficult to explain.

A triggering threshold for static stress changes has been generally set at 10 kPa (Stein, 1999), while dynamic stresses as small as 0.1–5 kPa are reported to have triggered seismicity (e.g., Brodsky and Prejean 2005; van der Elst and Brodsky 2010; Saccorotti et al., 2013). If such small stress thresholds can be also valid for static stress changes, then the consideration that some faults (before rupturing) experienced similar small stress magnitudes raises the possibility that static stresses may have played a role in the triggering (especially for the 2009 Acqua Panna-Marcoiano earthquake, for which the highest static stress changes are estimated). Notwithstanding the modest magnitude, Coulomb stress changes have loaded all the considered fault ruptures, and if Coulomb stress changes produced by any earthquake are cumulative, even small stress increases may “advance the clock” of the next earthquake, although it is currently unknown whether the considered faults exhibit such a sensitivity.

Another important issue regards the role of dynamic stresses in delayed triggering. Although dynamic stresses are transient, they can occasionally be converted into permanent deformation, for instance when seismic waves breach hydraulic barriers separating different fluid compartments (e.g., Manga et al., 2009). In addition, rapid fluid flow induced by seismic waves can remove temporary barriers within fractures and faults leading to pore pressure redistribution (Brodsky et al., 2003). Accordingly, Kilb (2003) envisaged that peak dynamic stresses induced by the 1992 M7.3 Landers earthquake would have played a role in triggering the next 1999 M7.1 Hector Mine earthquake. In particular, such a delayed triggering has been related to variations in the properties of fault zones (e.g., changes in friction and pore fluid pressures) induced by dynamic stress change oscillations (see also Parsons, 2005), whereby faults become weaker and thus more prone to rupture.

The role of dynamic shaking in delayed triggering could be evaluated considering the tensile strength of rock anisotropies that are conceivably common in fault zones. In this regard, the peak dynamic stresses on the order of 500 kPa estimated at some epicentral locations of the 2008–2019 Mugello seismic sequences (Figure 7), are included into the range of tensile strengths of rock anisotropies (i.e., mineral veins, schistosity planes, joints and faults), as experimentally determined by Shang et al. (2016). Therefore, suitably oriented dynamic stresses could potentially breach some of these discontinuities, opening routes for the escape of fluids from pressurized compartments or altering the scale length of the frictional contacts.

Applying these lessons to the north-eastern margin of the Mugello Basin, we raise the possibility that seismic perturbations would have played some role in triggering some of the mainshocks of the 2008–2019 sequence, although with modalities that may differ depending on the case (see Figures 7, 8, and Supplementary Tables S2, S3).

On the basis of the magnitude of static and dynamic stress changes, the 2009 Acqua Panna-Marcoiano earthquake was seemingly triggered by a combination of both types of stress changes. Static stresses are exceedingly small for the 2015 earthquake, and thus dynamic stress changes can be conceivably envisaged as the only potential driver. More ambiguous is the attribution of the 2019 Galliano-Sant’Agata earthquake, which occurred 10 years after the closest mainshock of the sequence (i.e., the 2009 Acqua Panna-Marcoiano earthquake). Also in this case dynamic stresses are much larger than static stresses (cf. Figure 7 and Supplementary Tables S2–S4), and so the actual role of the very small positive Coulomb stress change (i.e., 1 kPa) is uncertain.

Admittedly, the evidence for the seismic triggering is not solid. Nevertheless, the 4 M > 4 earthquakes clustered over a relatively short time interval (approx. 11 years), marking a significant increase of the seismicity rate with respect to the previous decades (see Supplementary Figure S4). In this respect, it is also worth noting that for that same magnitude range, the seismic catalog of the area is assumed to be complete since at least the second half of the 19th century (Visini et al., 2022). This spatial and temporal clustering may thus suggest that the ruptured faults were relatively close to failure. If so, this condition would be consistent with the hypothesis that earthquake-related stresses have played some role in the triggering of some 2009–2019 mainshocks.

The north-eastern margin of the Mugello Basin has experienced large historical earthquakes, particularly the 1542, Mw = 6 Scarperia earthquake and the Mw ≈ 6.4 1919 Vicchio earthquake (Rovida et al., 2021; Figure 1). Such large earthquakes have been correlated to the activity of the RFS (Bonini et al., 2016; this study). The overall framework resulting from the present analysis clearly points to the activation of distinct segments of the RFS. In particular, earthquake ruptures of 2009 and 2019 have been correlated to two RFS segments located in proximity of the macroseismic epicenter of 1542. Fault scaling relationships indicate that the 1542 rupture had a length of ca. 12 km and a width of ca. 8.5 km (Wells and Coppersmith, 1994). Notwithstanding the uncertainties of macroseismic epicentral locations, the fault rupture area caused by the 1542 earthquake conceivably includes the RFS segment activated during the 2019 sequence. This hypothesis is also supported by the observation that the highest macroseismic intensity of the 1542 earthquake (MCS = 9) has been documented at localities sited over the ruptured area of 2019 (e.g., Sant’Agata and Galliano; Figure 8). The evidence points towards the 1542 earthquake rupturing at least part, if not all the, the fault segment that ruptured in 2019. Less clear is the spatial correlation between the 1542 earthquake and the RFS fault segment that ruptured in 2009.

The four seismic sequences that we have relocated define similar fault segment geometries (Figure 8), with typical lengths and widths spanning the (2, 6)km range. In some cases, especially the 2009 and 2019 sequences, the subsurface fault segment matches in length, position and orientation to the associated surface expression of the RFS. This, coupled with the interpretation that the western end of the 1542 earthquake likely included the fault segment that ruptured in 2019, has the implication that the 1542 event involved the failure of multiple (likely two or 3) segments of the RFS.