The stratigraphic setting of the area containing 2016 Norcia seismic sequence is mostly represented by the successions belonging to the Umbria-Marche and (only marginally) Gran Sasso domains (UM and GS, respectively in Figure 1D). They reflect the evolution, starting from the late Triassic, of a passive continental margin characterized by platform-to-deep-water environments (Pierantoni et al., 2005; Cosentino et al., 2010; Barchi et al., 2012; Pierantoni et al., 2013 and references therein).

The units cropping out in the UM domain (Figure 1D) are represented, from the bottom to the top, by: Upper Triassic-Early Jurassic platform dolostones and limestones evolving to middle-late Jurassic basinal limestones and cherty limestones; Lower Cretaceous-Eocene (p.p) pelagic units marked by a significant increase of the pelitic content; early-middle Miocene hemipelagic successions (for a more detailed description we relate the reader to the focus provided in the next subsection).

In the Gran Sasso domain (GS in Figure 1D) the pre-orogenic stratigraphy is characterized by deposits settled along the transitional zone connecting the UM basinal domain and the carbonatic platform realm (Vezzani and Ghisetti, 1998; Cosentino et al., 2010; Adamoli et al., 2012; and references therein). The pelagic succession bears frequent calcareous-clastic turbidites interbeds. Levels of macroforaminifera-bearing breccias are also frequent in the middle Eocene p.p.-lower Miocene p.p. terms. Glauconitic limestones, calcareous marls alternating with calcareous-clastic gravity flows and by a late Miocene (p.p) clay interval rich in planktonic foraminifera predate the onset of the siliciclastic foredeep deposits. Typical Jurassic-Cretaceous carbonate platform succession crops out only in the GS westernmost sector.

The uppermost Miocene stratigraphy is mostly represented by the foreland basin deposits corresponding to the ‘Laga Flysch’ auctorum. The formation shows a complex internal architecture characterized by several vertical and horizontal transitions marked by changes in sandstone/pelite ratio and turbidite facies, by the presence of gypsum–arenite horizons (middle Messinian) and a volcaniclastic layer settled during the uppermost (5.5 Ma) Messinian (Scarsella, 1953; Crescenti, 1966; Ricci Lucchi, 1973; Ricci Lucchi, 1975; Centamore et al., 1991; Centamore et al., 1992; Odin et al., 1997).

The structural evolution of the area is the result of a late Miocene-early Pliocene compressional phase which originated the two arcuate regional thrusts (Figure 1D) of the Sibillini Mts (Koopman, 1983; Lavecchia, 1985; Barchi et al., 1988; Calamita and Deiana, 1988; Lavecchia et al., 1988) and the Gran Sasso Mt (Ghisetti and Vezzani, 1991). This phase was also responsible for doubling the UM sedimentary sequence, to the north, and for overthrusting of the UM succession on the GS, to the south (Olevano-Antrodoco line - Salvini and Vittori, 1982; Cipollari et al., 1997; Deiana et al., 2003). Folds and thrusts consequent to this phase define structural trends which are roughly N-S and E-W trending, in the Sibillini Mts- and Gran Sasso sectors, respectively.

The last tectonic phase recognizable in the area relates to the late Pliocene-Quaternary extensional tectonics (Brozzetti and Lavecchia, 1994; Lavecchia et al., 1994; Cavinato and and De Celles, 1999; Galadini and Galli, 2000).

Normal faults cross-cut the Mio-Pliocene folds and thrusts along NNW-SSE (to the North) to E-W (to the South) trending alignments with a prevalent western and southern dip, respectively. Their activity started in Lower Pleistocene (Calamita et al., 1994; Lavecchia et al., 1994; Cavinato and and De Celles, 1999) and promoted the formation of continental basins (e.g., the Norcia and Castelluccio basins) filled by hundreds of meters of lacustrine and fluvial deposits (Blumetti et al., 1993; Boncio et al., 1998; Boncio et al., 2004; Villani et al., 2019). Middle-Upper Pleistocene and Holocene successions (slope debris, alluvial fans interlayered with fluvial deposits, glacial- and fluvioglacial deposits) also occur within the continental basins and along the slopes of the main reliefs (Figure 1D).

In this section, we briefly describe (Figure 1D) the stratigraphic succession cropping out in the study area (see location in Figure 1D) and reported in the geological sections discussed later (Methods and Materials and Discussion and Conclusions).

As highlighted in the available literature maps (e.g., Barchi et al., 2012; Pierantoni et al., 2013), here the basinal Umbria Marche succession includes, from the bottom: 1) Early Jurassic carbonate platform deposits (Calcare Massiccio Fm), 2) Early to-Late Jurassic pelagic limestones, marls and cherty limestones (Corniola Fm, Rosso Ammonitico Fm/Marne del Serrone Fm - Calcari e Marne a Posidonia Fm, Calcari Diasprigni Fm) and 3) Lower Cretaceous-Lower Miocene p.p. limestones, marly limestones and marls (Maiolica Fm, Marne a Fucoidi Fm, and the FMs of the “Scaglia group”).

The aforementioned stratigraphy refers to the so-called “complete” succession, which was settled in most of the subsiding sectors of the Umbria-Marche basin. Conversely, the Jurassic palaeo-highs, which suffered a very slow drowning, are characterized by “condensed” or “reduced” successions, in which the Bugarone Fm replaces the whole -or part of- Jurassic pelagites.

Upward, the early-middle Miocene hemipelagic succession (Bisciaro Fm, Schlier Fm) consisting of marly siliceous limestones and clayey-silty marls with calcarenitic intercalations, predate the siliciclastic sin-orogenic deposition which is locally represented by the Arenarie di Camerino Fm (Serravallian p.p.-Messinian p.p.; Calamita et al., 1979).

Several logs of deep wells for oil exploration (e.g. Trevi 1, Antrodoco 1, Perugia 2, ViDEPI, 2020) highlight that the succession described above lies above Upper Triassic evaporites (Anidriti di Burano) which pass below to continental metapelites and quartzite (Verrucano Fm., Late Triassic- Permian (?)).

As recorded by instrumental seismicity, the Mt Vettore-Mt Gorzano area was characterized, by small earthquakes with magnitude M_L≤4.0 and very low seismicity rate (1981–2016) (Supplementary Figure S1) while neighboring active fault systems generated intense and important seismic sequences as the Norcia, 1979, M_W5.9, Colfiorito 1997, M_W6.0, and L’Aquila 2009, M_W6.3. The modern observations of such seismic sequences highlighted their multi-phase spatial-temporal evolution which progressively activated adjacent master faults, synthetic and antithetic segments and illuminates complex normal fault systems.

The Norcia 2016 seismic sequence showed similar characteristics. In fact, starting from EQ1, which occurred nearly 8 km NNW from the town of Amatrice, the seismicity migrated toward Cupi and Pievebovigliana towns, to the north, and to about 20 km south to Amatrice (see sect. Introduction).

The time-space analysis of the seismic sequence is shown in Figure 2. Here, we compare the seismicity, mainly associated with the VBF and GF (Figure 2A) with the ones of the study area (black inset in Figures 2A,B). The seismic sequence mostly occurred from August 2016 to the ending of June 2017 (75% of events, see Figure 2). An accelerated seismic release has been also later observed in April-May 2018 (about 10% of the total number of analyzed earthquakes).

The 2016–2018 seismicity vs time highlights that four significant periods (Figure 2D), in terms of number of earthquakes per day and seismic moment release, can be detected: 2016/08/24–2016/10/26, 2016/10/26-/10/30, 2016/10/30–2017/01/17 and 2017/01/18–2018/12/31, hereinafter called periods 1–4.

During period-1 (green circles), after the occurrence of EQ1, aftershocks migrated bilaterally. The study area (inset in Figures 1D, 2A) was not affected by seismicity (Figures 2A,B,D,E). For two months, the intense seismic activity (Figure 2D) remained confined between Mt Gorzano and the southern segment of VBF (Vettoretto-Redentore segment in Brozzetti et al. (2019)). After the occurrence of EQ2, the seismicity migrated to the northern segment of the VBF system and mainly concentrated in between Norcia and Pievebovigliana towns for about 30 km (Figure 2A). As shown in Figures 2B,E, only after EQ2, the Pievebovigliana sector became involved in the seismic sequence (period-2, orange circles). The major event (M_W4.1) of this period occurred on 2016/10/27 in the southern part of the study area (focal mechanism n°1 in Figure 2C) and produced the first increase of the cumulative moment release (red curve in Figure 2E). The number of events per day was over 200 (Figure 2E), about one third of the total events (Figure 2D), with a noteworthy increase of seismicity with respect to the seismic activity that had previously characterized the study area (maximum number of events per day = 10, Supplementary Figure S1 - lower inset). After EQ3 (M_W6.5), the VBF system was involved for about 65 km (period-3, yellow circles). During period 3, the study area was affected by more than 100 events/day and the seismicity extended to about 4 km north and west to Pievebovigliana (Figure 2B). This period was characterized by a significant increase of seismic moment release (red line in Figure 2E) with the occurrence of five earthquakes with M_W≥3.9 (yellow beachballs in Figure 2C), never observed before the 2016 Norcia seismic sequence. In fact, in the instrumental period (1985–2016) the Pievebovigliana sector was mainly struck by earthquakes with M_L≤3.5 (Supplementary Figure S1). The major event of the period-3 had a magnitude of M_W4.8 and occurred on November 1, 2016. During the period-2 and -3, the study area was uniformly affected by aftershocks and the map distribution of earthquakes well correlate with the average trend highlighted by focal solutions (average focal mechanism in the inset of Figure 2C).

On January 18, 2017, period-4 began with the occurrence of EQ4 that involved the southern portion of Gorzano fault. The increase of seismic activity of the study area started on 2017/01/27, 9 days after the occurrence of EQ4 (Figure 2A). The number of events/day and cumulative seismic moment release highlight two main sub-periods, some minor peaks (black arrows in Figure 2E) and the average number of earthquakes that remained ten times over the one observed before the occurrence of the Norcia seismic sequence. It is very interesting to note that, in this phase, the seismicity was characterized by additional westward and northward migration from Pievebovigliana. During the sub-Period-4/2, the seismicity was mainly concentrated in the study area (compare Figures 2D,E). The most significant earthquake (M_W4.6) of this period occurred on April 10, 2018 during the sub-period 4/2 (Figure 2E). This event represents the last significant increase in term of moment release.

All the focal mechanisms (TDMT solutions in Scognamiglio et al., 2006) in Figure 2C, show normal-normal oblique sense of motion with a SW-NE trending T axis, apart from a strike-slip cluster of small events M_W≤3.5 located between Cupi and Fiastra with T axis compatible with the average one (inset in Figure 2C).

We tested the effectiveness of the topographic derivative first along the western slope of the Mt Vettore fault system (inset a in Figure 1D) where we explored the coherence between the curvature values with (already) mapped fault segments. In this study, we used those reported in Brozzetti et al. (2019). The workflow includes the following steps:

1) Apply a low-pass filter on DEM, to remove roughness in the surface (and/or local anthropic artifacts). We used the Focal Statistics tool and a Rectangular cell (50×50m) Neighborhood–Mean statistic setting;

The interpreted lineaments are reported in Figure 3C. For ease of understanding, we refer to each fault segment to its number (from n°1 to 11)

Following the test along the Mt Vettore fault, we computed a curvature map also in the sector immediately ahead of the VBF northern tip, i.e., between the Pievebovigliana village and the Mt Val di Fibbia area (insets b and c in Figures 1D, 4A–D), where the seismicity did not correlate to any relevant discontinuity. We applied the same methodology and we interpreted the map according to the following criteria:

− comparing continuous (well evident) pseudo-linear tonal edges with the rock bedding and lithological boundaries as reported in available maps at scale 1:100,000 (Regio Ufficio Geologico, 1941; Servizio Geologico d’Italia, 1967) and 1:25,000 (Pierantoni et al., 2013), to ease the distinction between the latter and possible tectonic lineaments (Figure 4E–G);

We mapped lineaments (1–7 in Figures 4B,D,F,G) and we also provided a preliminary dip direction of the inferred faults exploiting, locally, the ArcMap ‘Profile Graph’ tool to indicate the inferred raised and lowered blocks. Based on the trend and dip of the VBF system, our first attempt was to highlight possible W-dipping normal faults (NW-SE striking) considering the prevailing geometry of the active fault system in the epicentral area. Nevertheless, the existence of possible antithetic (E-dipping) structures was not ruled out from our map interpretation.

In the sector between the Cupi village and Mt Val di Fibbia (inset b in Figure 1D), resistant rock types crop out and flat surfaces on top of the Meso-Cenozoic carbonates appear to be offset (Figure 5A). The only mapped structural element is represented by the north-western tip of the VBF (Cupi-Ussita fault section–CUS in Brozzetti et al. (2019)) which affects partially the carbonates, and only toward the west.

According to the prevalent dip of the normal fault system responsible for the entire seismic sequence, we investigated the role that west-dipping (unmapped) normal faults may have in displacing possible PaS remnants in the sector. We integrated the previous analysis with the computation of combined slope and aspect maps to enhance the morphostructural setting. We computed the derivatives (in ArcMap environment) according to the following steps:

1) slope map - to isolate (flat) sub-areas dipping less than 15°;

Moreover, we drew five NE-SW oriented topographic profiles (labels Tp1 to Tp5 in Figure 1D). The profiles were drawn perpendicularly to the trend of the neighbor VBF system assuming, by analogy, that faults active in the late Quaternary should form within the same stress field acting along the central Italy extensional alignment. We computed a possible vertical offset (Figure 5C) and discuss them also by the light of the stratigraphic setting of the limestones.

A new geological field survey at scale 1:10,000 was carried out in the sectors neighboring the village of Pievebovigliana and the Mt Val di Fibbia (see insets b, c in Figure 1D). The survey was directed in the area where most of the interpreted lineaments concentrated and where topographic anomalies suggested the possible interference of active tectonic structures with the landscape evolution. The survey allowed us to differentiate between possible tectonic- or lithologic control on the curvature value outcomes.

We compared the survey outcomes with the long-term deformation reconstructed for the sector. In this perspective, we drew four 20 km-long interpretative geological cross-sections along transects perpendicular (N60°E) to the trend of the main tectonic structures (traces in Figure 1D, n°1 to 4). The cross-sections were drawn starting from available geological maps at scale 1:40,000 (Pierantoni et al., 2013) and 1:25,000 (Barchi et al., 2012). Where the coverage was not attempted by the previous cartography, available maps at scale 1:100,000 (Regio Ufficio Geologico, 1941; Servizio Geologico d’Italia, 1967) were used. Evidence of normal faulting as deduced from the previous approaches were also incorporated into the new geological cross-sections (Figure 7). For the in-depth interpretation of the cross-sections, we relied on the deformation style proposed in specific papers on the topic (Lavecchia, 1985; Barchi, 1991; Barchi et al., 1998; Lavecchia et al., 2016).

To validate the seismogenic nature of the PBF and to find evidence of an association between the PBF fault (at the surface) and the seismicity (at depth) during the evolution of the seismic sequence, we examined the relocated earthquake dataset available in Chiaraluce et al. (2017). We investigated the in depth-distribution of the hypocenters available for the study area in the time interval 2016/08/24–2016/11/30 and we compared the seismicity clustering with the collected geological and geomorphological evidence. We performed the analysis collecting the events with depth 10 km, within a half-width of 2.5 km starting from four transects (coinciding with the geological cross-sections reported in Figure 7) and moving northwestward. We imported all the layers (e.g., surface geology, geological cross-sections and earthquakes) in the Move Suite (see sect. Methods and Materials) to compare the data in a 3D environment (see the 3D block model in Figure 8).

We investigate the location and geometry of the highest stress changes within the study area and assess if they were compatible with the PBF activation, computing the stress transfer induced by the largest subsequent earthquakes of Norcia 2016 seismic sequence EQ1, EQ2 and EQ3 on the study area. EQ4 was not considered because the interactions acting on volumes are proportional to the rupture dimension (e.g., Hardebeck et al., 1998). Hence, EQ4 being the least energetic (M_W5.5) and the most distant event with respect to the study area, we considered the possible interaction negligible. The imparted stress was computed both considering the single events and their cumulative effects. We used the Coulomb code 3.4 (Lin and Stein 2004; Toda et al., 2005) and we analyzed the results considering the stress changes along a specific plane having the main geometric characteristic of the Mt Vettore-Gorzano fault System (method 1- M1), inserting into the calculations the influence of the regional stress (method 2- M2), and projecting the stress changes on the reconstructed PBF (method 3- M3).

To achieve our aims, we considered the EQ1, EQ2 and EQ3 source models, with variable slip, proposed by Chiaraluce et al. (2017), and a friction coefficient (μ) of 0.4 generally reported in the literature (Sibson, 2000) and previously used for central-southern Apennines studies (Troise et al., 1998; De Natale et al., 2011).

In M1, we choose to project stress changes on generical receiver faults having a 156° strike, 50° dip and –90° rake, that is the average geometry of the main active fault alignments constrained from the EQ1–3 source models (Chiaraluce et al., 2017).

To compute the stress imparted on the surrounding crust volume with optimally oriented stress calculations (M2) we considered the background regional stress field evaluated for central Italy by Ferrarini et al. (2015) (σ1 = 292/85, σ2 = 139/04, σ3 = 048/02).

We iteratively ran the stress changes considering the temporal evolution of seismic sequence and comparing, in map and cross-sections, the redistribution of stress after each significant event Specifically, the map and cross-sections represent the maximum value of Coulomb stress changes over the considered depth range (0–10 km, with an incremental step of 1 km). Figures 10, 11 show the map of cumulative stress changes (EQ1-EQ2 and EQ1-EQ3 respectively), computed at 8 km of depth, based on the hypocentral depth of the EQ1 and EQ2 main events, and the related in-depth sections spaced 5 km, distributed along the study area in correspondence of the four geological ones.

We overlaid to such stress cross-sections the aftershock locations made available in Chiaraluce et al. (2017) to verify the correspondence between the occurrence of earthquakes and the positive Coulomb stress lobes. In the Supplementary Figures S3–S9 are reported all the maps and related cross-sections obtained considering the consecutive occurrence of EQ1–3 and the single events.

We, then, focus on the stress imparted on Pievebovigliana geological structure defining a 3D fault model. To this aim, we used the Mildon et al. (2017) code, which allows variable strike faults, and we used the results obtained combining geological and seismological data (field data, geomorphological analysis, seismicity depth distribution and focal mechanisms). We used the same parameters of the previous calculations) and ran several tests varying the dip of receiver source model (50°–60°) and analyzed the stress imparted on PBF by the single-source models and by their cumulative effects. Figure 12 and the Supplementary Figure S9 synthesize the obtained results.

Along the western slope of the Mt Vettore fault system (Figures 1, 3), the analysis of the curvature derivative successfully identified most of the active faults recently mapped in the literature. With the exclusion of the mountain topographic crests, the comparison between the fault segments reported in Figure 3A and the curvature value distribution pointed out that most of the outcropping fault traces correspond to (Figure 3B):

- intermediate curvature values between two pseudo-linear- and coupled maximum (red) and minimum (blue) ones. These latter have been found related to the concave and convex shapes (on the hanging- and footwall blocks, respectively) resulting from the intersection of the outcropping fault plane and the slope. Evident examples (transparent gray lines in Figure 3C) relate to the points n°2, 4 and 9;

- breaks in the minimum (blue) values (points n°1, the northern tip of the n°5, 6, 7 and 8), related to the fault offset affecting locally the pseudo-linear ridges of the mountain crests;

- breaks in the maximum (red) values widely outcropping in the sector even if with a weak “derivative signal” (points n°3, 10 and 11). These latter have been found related to fault planes often located along the mountain slopes and bearing slope deposits on the hanging wall blocks; and/or in the Castelluccio plain, thus often totally buried (see for detail the geological map provided in Figure 4 Brozzetti et al. (2019)). Most of the fault segments belonging to the VBF system with a length >1.5–2 km showed a good match with the anomalies observed in the profile curvature map. The results supported the effectiveness of this DEM derivative as a useful tool at least in the perspective to direct field survey where no other information is available from the literature background.

In the interpretative geological cross-sections (Figure 7), it is possible to recognize the different deformation phases related to both the Neogene and Quaternary tectonics.

The pre-Quaternary deformation is overall recognizable in the compressional structures rooted in the Upper Triassic evaporites (see subsect. Structural-Geological Setting) and in the deeper basement (Lavecchia et al., 1988; Coward et al., 1999; Speranza and Chiappini, 2002; Porreca et al., 2018). They deform mostly the basinal succession giving rise to both cylindrical- and box anticlines, the latter localized in correspondence to the main thrust planes.

The Plio-Quaternary change of the tectonic regime is evident from the superposition of the normal faults which offset the pre-existing structures and led to the formation of the intra-mountain basins. Along the geological cross-sections, the Quaternary extensional system is recognizable (from SW to NE) in the Colfiorito (CoF, cross-sections n°1 and n°3), Norcia (NF in cross-sections n°4) and Cupi-Ussita (CUS) faults (cross-sections n°3 and n°4). The latter bears the northernmost evidence of coseismic displacements (Figure 9A, according to Brozzetti et al., 2019).

We also introduced the evidence of the PBF normal faulting as we deduced from the field survey and morphotectonic analysis. Our findings suggest that the deformation does not manifest at the surface as a single lineament but, rather, as distributed segments evident in the field as far as to the Pievebovigliana village (segments a to g Figure 9A and Table 1).

In fact, the earthquake-fault association along cross-section n°1 suggest the activation of a structure up to a depth of ∼5 km (Figure 8A), even where no surface evidence was observed. The source was illuminated starting from October 26 and the subsequent seismic activity, in the time window covered by the catalog (up to November 30, 2016) well aligned, also in map view, with the PBF strike and its suggested position (Figures 2A, 9A). No significant events fall within the buffer of cross-section n°1 (only one M_L3.6 event on October 31, 2016, Figure 8A) and considering the time window explored. Nevertheless, if we also project the most energetic event (April 10, 2018, M_W4.6) reported in Scognamiglio et al. (2006) (star in Figure 9A and Table 1) which felt in the study area and the buffer of the cross-section (see also Figures 1C,D, 2A) we could infer its association with the PBF at depth (Figure 8) and in correspondence of the inferred segment h, at the surface. For this reason, we suggest that the seismicity projected along cross-section n°1 intercepts the PBF northwestern tip.

Along cross-section n°2, we reported the evidence from stops discussed in sects. 4.2.2 (Figure 6B,E) while, along the section n°3, we considered the offsets affecting PaS remnants (see subsection Sector ‘Cupi-Mt Val Di Fibbia’ and Figures 5C,D). The deformation we depicted at the surface, along these cross-sections, found also confirmation in the 3D analysis of the hypocenter distribution (Figure 8). The seismicity projected shows clustering in correspondence of the PBF along both cross-sections n°2 and n°3.

In detail, along cross-sections n°2, the seismicity illuminates the suggested fault plane up to ∼7 km depth. The October 30 M_W4.0 earthquake (event 2 in Figures 8A,B) locates along the PBF down-dip prolongation crossing the cross-section (segments d, e, f, g in Figure 9A and Table 1), while the November 3 (M_W4.7) locates at greater depth (event 4 in Figure 8, star in Figure 9 and Table 1). This event could be associated to a deeper PBF patch (Figure 8B), even considering the low (24°) dip angle of the west-dipping plane, or being the result of the activation (at the PBF hanging wall) of minor discontinuities (e.g., ENE-dipping faults) whose existence and attitude are suggested also along the adjacent section n°3. The activation of antithetic structures in the study area has been recently highlighted also by the dataset Michele et al., 2020 (provided in Spallarossa et al. (2021)) and elsewhere in the epicentral area of the Norcia 2016 seismic sequence. Even if not outcropping, they have been associated also with energetic events (e.g., August 24, 2016, M_W5.4), at the hanging wall of the VBF system (Porreca et al., 2018).

Along cross-section n°3, two clusters of hypocentres are evident in the first 5 km depth: the easternmost ones coincide with the segments associated to the PBF (b, c, d in Figure 9A), and the other with the CUS. Both PBF and CUS host M3.9 events (e.g., November 12 and 27) (Figures 2C, 8A). Most of the seismicity clusters around the CUS even if several M>3.5 earthquakes fall around the PBF. Minor antithetic discontinuities are also illuminated by the contours (Figure 8B) as well as an ENE-dipping low-angle (∼15°) plane possibly coinciding with a regional low-angle normal fault (see Chiaraluce, et al., 2017 and references therein).

Finally, along cross-section n°4, some scattered seismicity is observable around the CUS. Two moderately energetic events are reported (October 27 M_W4.4 earthquake (event 1) and the November 1 M_W4.8 one (event 3) (Figure 8, stars in Figure 9 and Table 1). Nevertheless, the location of the event 1 does not fit well the fault plane position and the event 3 has a very shallow location (Figures 8B), tentatively associable with the CUS inner splay (Figures 8, 9A) or suggesting possible linkage (at depth) of the PBF with CUS. In this sector, we did not find other evidence of normal faulting aside from an anomalous syncline at the hanging wall of the thrust T3. The back-limb fold reconstructed from the strata attitudes (Figure 7) could be explained with an ‘incipient’ (and/or distributed deformation) lowering of the PBF hanging wall, along the fault south-eastern tip. On the other side, the scattered seismicity located east of the CUS, including some M>3.5 earthquakes at a depth between 2 and 3 km, could be related to minor discontinuities activated during the sequence. We are not able to discriminate among these different hypotheses and we associated the southern tip of localized deformation to the last evidence found along segment a (Figure 9A and Table 1).

Considering all the collected results and imaging a single master structure responsible for observed distributed deformation, we propose as reliable geometry for the PBF that of at least ∼13 km-long master normal fault which strikes ∼ N155°E, dips SW and is arranged in right-lateral en echelon setting with respect to the VBF system (blue line in Figures 9, 10).

The average west-dipping focal plane we computed in the study area, where the PBF strikes within (Figure 2A), has a dip of 43°. On the other side, the fault attitude measurements reported in Figure 6C show an average dip of 78°. Besides, data collected in the surroundings of the CUS (Testa et al., 2019) show average dip of 65° while all along the central- and southern sections of VBF the structural data in Brozzetti et al., 2019 highlighted average dip values from 65° to 71° (associated to long-term- and coseismic planes, respectively). Considering the heterogeneity of the data source we used a PBF theoretical dip value of 60° to model this fault as input data for the Coulomb stress computation.

The model calculations of the elastic stress change generated by EQ1 suggest that the stress imparted on the study area is negligible (<0.1 bar; see map and cross-sections 1–4, Supplementary Figures S3, S4). This result is also consistent with the lack of seismicity in the area during the period-1 (Figure 2).

On the contrary, the occurrence of EQ2 was important for the redistribution of stress. The largest Coulomb stress changes along a generic receiver fault (M1, subsect. Coulomb Stress Change) identify a well-defined volume of stress changes in the area where the PBF has been suggested (Figure 10). Along the cross-sections 2 to 4, it is possible to identify positive maximum Coulomb stress changes (0.5 bar) greater than the minimum (0.1 bar) commonly required to contribute to the triggering process (Reasenberg and Simpson 1992; Hardebeck et al., 1998). This result is also confirmed by M2 analysis performed along optimally oriented normal faults (Supplementary Figure S5). The negative stress changes on the cross-sections 3-4 are probably due to the high slip values of the causative source EQ2 (see also Figure 12). Moreover, selecting the seismicity of the relocated dataset (Chiaraluce et al., 2017) in the time interval 2016/10/26–2016/10/30, we note the high correlation of the aftershocks and the volumes undergoing increasing Coulomb stress.

The cumulative Coulomb stress changes imparted by EQ1, EQ2, EQ3 are shown in Figure 11 (Supplementary Figure S6) and depict a similar pattern of positive and negative lobes although with higher stress magnitudes. Additionally investigating the contribution of each large event on stress redistribution and modeling only the EQ3 occurrence, it is possible to observe (Supplementary Figures S7, S8) that the resulting cross-sections show evenly distributed stress changes and it is not possible to identify well defined positive stress volumes (cross-section 4 apart). Hence, we can argue that absolute values stress changes induced by EQ3 are added to the previous one enhancing and minimally perturbing the stress pattern induced by EQ2. Overlaying the seismicity, it is possible to observe that a high percentage of events occur on positive stress changes (Hardebeck et al., 1998).

Moreover, comparing the cross-section of Figures 10, 11, it is interesting to point out the westward enhanced stress, along a volume antithetic to the main one (black arrows in the cross-sections 2-3 of Figure 11) where the seismicity migrated and accelerated seismic release was observed in April-May 2018 (period-4 in Figure 2).

The results of stress changes modeling induced by EQ1, EQ2, EQ3 on the reconstructed Pievebovigliana fault (Figure 12) are fully consistent with the ones obtained with M1 and M2. The analysis indicates that the positive stress changes were transferred to the active Pievebovigliana receiver fault, after EQ1-EQ3, with different levels. The EQ1 minimally perturbed the fault while EQ2 and EQ3 imparted significant stress values greater than 0.5 bar along the southern portion of the fault. EQ2 was the crucial event in favouring the activation of PBF. EQ3 induced high-stress values around the southern tip of the PBF, and this result is particularly important because the tips of faults are typically the weakest zones. The cumulative effect enhanced the positive stress perturbation and made PBF a good candidate to accommodate this stress as demonstrated by the significant occurrence of seismicity in the study area until the end of 2018.