Detection of Spatial and Temporal Stress Changes During the 2016 Central Italy Seismic Sequence by Monitoring the Evolution of the Energy Index

We consider approximately 23,000 microearthquakes that occurred between 2005 and 2016 in central Italy to investigate the crustal strength before and after the three largest earthquakes of the 2016 seismic sequence (i.e., the Mw 6.2, 24 August 2016 Amatrice, the Mw 6.1, 26 October 2016 Visso, and the Mw 6.5, 30 October 2016 Norcia earthquakes). We monitor the spatiotemporal deviations of observed radiated energy, ES, with respect to theoretical values, ESt, derived from a scaling model between ES and M0 calibrated for background seismicity in central Italy. These deviations, defined here as Energy Index (EI), allow us to identify in the years following the Mw 6.1, 2009 L’Aquila earthquake a progressive evolution of the dynamic properties of microearthquakes and the existence of high EI patches close to the Amatrice earthquake hypocenter. We show the existence of a crustal volume with high EI even before the Mw 6.5 Norcia earthquake. Our results agree with the previously suggested hypothesis that the Norcia earthquake nucleated at the boundary of a large patch, highly stressed by the two previous mainshocks of the sequence. We highlight the mainshocks interaction both in terms of EI and of the mean loading shear stress associated to microearthquakes occurring within the crustal volumes comprising the mainshock hypocenters. Our study shows that the dynamic characteristics of microearthquakes can be exploited as beacons of stress change in the crust and thus be exploited to monitor the seismic hazard of a region and help to intercept the preparation phase of large earthquakes.

From both small-scale laboratory experiments and megathrust studies, it emerges that the spatiotemporal evolution of microearthquakes (called acoustic emission in the former case) provide information on seismic friction and coupling, which in turn may help to predict laboratory earthquake and to intercept the preparation processes of large earthquakes (e.g., Bouchon et al., 2013;Rouet-Leduc et al., 2017;Schurr et al., 2014;Socquet et al., 2017;Tape et al., 2018). Considering the crust as a complex system, microearthquakes can thus help us to learn about the system's condition and its evolution toward instability. Seeber and Armbruster (2000) proposed to use microearthquakes as beacons for the mechanical state of the crust by studying their static Coulomb stress changes. The key issue is that microearthquakes are too small to interact (Stein & Lisowski, 1983) and they contribute little to crust deformation. Nevertheless, their characteristics and distribution in time and space are very sensitive to stress changes.
Seismic sources can be characterized in terms of seismic moment (M 0 ) and radiated energy (E S ), which are parameters of particular interest to monitor rupture processes and the crust condition evolution. These parameters are defined independently from any source rupture model and look at the seismic source from different perspectives (Bormann & Di Giacomo, 2011). M 0 provides a static measure of the earthquake size, while E S provides information on the rupture kinematics and dynamics Picozzi et al., 2017Picozzi et al., , 2018. Therefore, the ratio between E S and M 0 , known as "scaled energy" (Ide & Beroza, 2001) or "slowness parameter" (Newman & Okal, 1998), provides unique information on how efficient the rupture is in radiating seismic energy, which in turn is related to the effective applied stress in the source volume (Walter et al., 2006). McGarr (1999) studied E S and M 0 for earthquakes occurring in different environments and circumstances (i.e., in laboratory, induced and natural earthquakes). He found that the seismic efficiency (i.e., the fraction of released energy in the rupture process i.e., radiated seismically) of all considered earthquakes had a narrow range of variation and an upper bound equal to 0.06. This evidence led McGarr (1999) to conclude that seismic efficiency is controlled by a few fundamental aspects that are common to both shallow earthquakes and laboratory frictional stick-slip events. Implications of McGarr's work (1999) are that laboratory results can provide important insights to interpret earthquake energy budget data, and that measurements of E S and M 0 can be used to estimate shear stress at seismogenic depths.
Interestingly, Hulbert et al. (2019) showed that a machine learning algorithm could decipher the complex temporal pattern of the energy released by acoustic emissions during fast and slow laboratory earthquakes, and it was able to predict both the occurrence time and duration of the main events. Rouet-Leduc et al. (2020) found that regarding the automatic identification and characterization of slow earthquakes in the Cascadia region, the statistical feature identified as most informative was very similar to those found in the laboratory (i.e., those related to signal power and thus energy). Similar observations about the importance of radiated energy were also made by Kano et al. (2018) for the characterization of slow slip transients in the brittle-ductile transition zone along subduction plates in Japan. These authors concluded that by monitoring tremor activity, and especially the associated radiated energy, the strength of tremor patches in transition zones worldwide can be inferred and monitored.
The energy-to-moment ratio is also important for large tectonic earthquakes. High scaled energy values indicate earthquake ruptures that radiate larger energy at high frequencies compared to what would be expected given the earthquake size as measured by its seismic moment, with significant implications for hazard assessment (Choy & Kirby, 2004).
All these recent studies covering a wide range of spatial scales support the idea of McGarr (1999) about the importance of E S , and its scaling with M 0 , for characterizing the crustal stress at which earthquakes occur and for monitoring the evolution of crustal processes leading to earthquakes. Picozzi et al. (2022) investigated the temporal evolution and spatial distribution of the ratio E S to M 0 for microearthquakes that occurred during the preparatory phase of the Mw 6.1, 2009 L'Aquila earthquake, with the aim of testing whether these parameters could capture hints of preparatory processes leading to the mainshock. To this purpose, Picozzi et al. (2022) studied the deviation of E S for foreshocks from a reference model (calibrated during a period characterized by background seismicity) describing the scaling of E S with respect to M 0 .
The 2016 central Italy seismic sequence led to recording more than 500,000 earthquakes through mid-2017 (e.g., Spallarossa, Cattaneo, et al., 2021;Zhang et al., 2019), covering nearly 60 km along the strike of the normal fault system in central Apennines (Chiarabba et al., 2020). Such big data called for a standardized strategy for the characterization of earthquake source parameters. Spallarossa, Picozzi, et al. (2021) implemented an innovative service for the rapid assessment of seismic moment, M 0 , and radiated energy, E S , in central Italy (i.e., RAMONES, http://www.distav.unige.it/rsni/ramones.php). RAMONES exploits continuous data streams stored in free repositories (i.e., ORFEUS-EIDA, IRIS, and DPC). The automatic procedures implemented in RAMONES consist of events detection, location, magnitude estimation, and extraction of features from S-wave waveforms (for more details, see Scafidi et al., 2016;Spallarossa et al., 2014;Spallarossa, Picozzi, et al., 2021). M 0 and E S are then estimated using features calculated for S-waves and applying empirical station correction factors and attenuation models calibrated for a large region of central Italy (i.e., a region bounded by 40.0°N and 44.5°N in latitude and 10.50°E and 16.50°E in longitude). We summarize the procedure implemented in RAMONES for estimating E S and M 0 in Text S1 of Supporting Information S1.
Here, we consider 23,273 earthquakes (M W > 1.8, depth between 2.5 and 17.7 km) to investigate the crustal strength before and after large earthquakes in central Italy during the 2016 seismic sequence. Starting from estimates of E S and M 0 as provided by RAMONES, we study the spatiotemporal evolution of EI to assess if it can help to shed light on the crustal conditions before and after the three mainshocks of the seismic sequence. In other words, as proposed by Seeber and Armbruster (2000), we want to explore whether microearthquakes can act as beacons of stress change.
Several previous works have studied the 2016 central Italy sequence from multiple perspectives. Among them, Improta et al. (2019) provided an overview of the complex interaction of shallow normal faults with the inherited deeper thrust; Kemna et al. (2021) studied the temporal evolution of static stress drop; Chiarabba et al. (2020) showed the existence of precursory velocity changes before the Mw 6.5, Norcia earthquake; Malagnini et al. (2022) observed time variation of seismic attenuation; and Sugan et al. (2022) interpreted the microearthquake temporal and spatial evolution as progressive unlocking of the main faults.
With respect to previous studies, here we provide a novel and complementary overview of the temporal and spatial evolution of the 2016 central Italy seismic sequence in terms of microearthquakes dynamic characteristics. In the first part of the work, we estimate the energy index, and we study its spatial distribution before the three mainshocks. Then, we focus on crustal volumes around the three mainshock hypocenters and we investigate their interaction by studying the temporal evolution of both EI and the mean shear stress associated with the microearthquakes, derived following McGarr (1999).

Data
We use estimates of M 0 and E S provided by the RAMONES service (Spallarossa, Picozzi, et al., 2021). In RAMONES, the P and S onsets are automatically detected, the events located, and their local magnitudes estimated. Then, the squared velocity signals integrated over the S-wave time window (IV2 S ) and the S-wave peak displacement (PD S ) are computed for each waveform, corrected for propagation and site effects, and used to estimate M 0 and E S at each seismic station. The average of the station estimates provide the M 0 and E S associated to a given earthquake.
The attenuation models and site station correction factors used in RAMONES were calibrated by Spallarossa, Picozzi, et al. (2021) considering 6,515 earthquakes recorded in central Italy and detected by 464 stations (∼210,000 waveforms). The high number of events and stations in the calibration data set made the attenuation models and site correction terms very robust for the region under study. Since January 2020 RAMONES has been operating with a configuration set to daily updates (http://www.distav.unige.it/rsni/ramones.php). More details on the RAMONES' procedures are provided in Text S1 of Supporting Information S1. We refer to Spallarossa, Picozzi, et al. (2021) for a detailed description.
Here, we consider the period from January 2005 to December 2016. As we said, the analyzed seismicity covers a wide area, which includes the seismicity occurred on the Paganica Fault, where the Mw 6.1, L'Aquila earthquake occurred in 2009, the Campotosto Fault, where one of the largest aftershocks of the Mw 6.1, L'Aquila earthquake occurred (i.e., the Mw 5.4 earthquake that occurred on 7 April 2009 at 19:47 UTC time), the Laga Fault system (LFS), where the Mw 6.2 Amatrice earthquake occurred, and the Vettore-Bove Fault system (VBFS), where the 26 and 30 October 2016 Mw 6.1 Visso and Mw 6.5 Norcia earthquakes occurred, respectively (Figures 1a  and 1b). We represent faults as boxes projected at the surface (Luzi et al., 2017;Malagnini et al., 2012). Uncertainties in event location are mostly within 1 km both horizontally and vertically ( Figure S1 in Supporting Information S1). The number of stations per event used to estimate M 0 and E S ranges mostly between 10 and 40 ( Figure 1 and Figure S1 in Supporting Information S1). Picozzi et al. (2022) considered a smaller portion of this data set (i.e., from 2005 to 2009, for a total ∼6,000 earthquakes) to investigate the preparatory process of the Mw 6.1, L'Aquila earthquake in terms of the E S versus M 0 scaling, and hence focusing on the radiated energy of foreshocks with respect to that of previous background events. In this study, we extend the observational period to the end of 2016 aiming to include the three mainshocks of the 2016 central Italy sequence and their aftershocks ( Figure 1). The original catalog consists of 41,794 earthquakes recorded at 393 seismic stations ( Figure 1b). However, we select the data considering a minimum threshold magnitude equal to Mw 1.8. Such magnitude threshold was chosen both by Spallarossa, Cattaneo,  Spallarossa, Picozzi, et al. (2021) to calibrate the attenuation models for RAMONES and by Bindi et al. (2020), who evaluated the resolution limits of source parameters for earthquakes occurring in the same region considered in this study. According to their results, for earthquakes with magnitude equal or higher than Mw 1.8, M 0 and E s can be considered free from bias due to unaccounted attenuation effects. We will return on this issue when discussing the results.
Furthermore, we consider events with depths between 2.5 and 17.7 km. We exclude very shallow and deep earthquakes because we consider them not representative of the dynamic conditions at which the main ruptures occurred.
The resultant M 0 and E s for a reference period (see Section 3.1) are used to define a reference E S -to-M 0 scaling model. Then, the energy index, EI, is estimated considering events with magnitude within the range Mw 1.8 and Mw 3.5. For the considered selection in magnitude and hypocentral depth, the EI data set consist of 23,273 earthquakes in total ( Figure 1a). We exclude events with Mw > 3.5 because, as suggested by Kanamori and Heaton (2000), the rupture process of larger magnitude events could be influenced by melting and fluid pressurization, which would alter the E S to M 0 scaling. Furthermore, since we are interested in studying earthquakes as beacons for the mechanical state of the crust, we wish to consider events that are sufficiently small so that they do not interact with each other. Figure 2a shows the temporal evolution of seismicity in terms of magnitude for the 23,273 earthquakes considered in this study (we report the distribution of data in magnitude in Figure S2 of Supporting Information S1). We identify two periods characterized by spatially diffuse seismicity ( Figure 2b) and the absence of events with significant magnitude (Mw > 4.5, Figure 2a). These two periods of time occurred before and after the 2009 L'Aquila and before the 2016 central Italy seismic sequences. We use the seismicity that occurred over these two "interseismic" periods to calibrate a reference scaling model between the base-10 logarithm (indicated as "log") of M 0 , and E S . Period #1 spans from January 2005 to December 2007 for a total of 461 earthquakes, while period #2 includes 1,524 earthquakes from January 2011 to December 2015 ( Figure 2b shows the distribution of earthquakes for these two periods). Figure 2c shows the E S and M 0 estimates for the whole data set and for the two identified periods. The E S -to-M 0 scaling models for periods #1 and #2 are consistent with each other (Figure 2c, model parameters in Table 1). Thus, we combine the E S and M 0 estimates for the two periods and we calibrate an E S -to-M 0 scaling model (i.e., model #1 + #2, Table 1) that we consider as the reference model for the background seismicity in the investigation area. The E S -M 0 scaling for central Italy obtained from the RAMONES estimates agrees with those from other studies and other seismic regions ( Figure S3 in Supporting Information S1). In Figure S4 of Supporting Information S1, we show the standard error associated to the E S -M 0 scaling, with data colored per number of stations used. The high spatial density of the Italian seismic networks allows for dozens of single-station source parameter estimates for most earthquakes, which leads the E S estimates to be particularly robust.

Reference E S -to-M 0 Scaling Model and Energy Index
Following Brown and Hudyma (2017) and Picozzi et al. (2022), we compute the differences EI = log(E S ) − log(E St ) between the experimental E S estimates and theoretical E St values derived from the median E S -to-M 0 reference scaling model. Positive EI values indicate that the earthquakes have radiated more energy per unit-slip and unit-area (i.e., per seismic moment, M 0 ) than expected according to the reference model. On the contrary, negative EI values are associated with earthquakes showing an excess of slip or larger rupture area and lower stress drop with respect to what is expected from the reference model.

Figures 2d and 2e
shows the distribution of EI with respect to M 0 and with data colored per hypocentral depth for the calibration and the complete data sets, respectively. Figures 2f and 2g shows the distribution of EI for the two data sets. As expected, EI values are on average evenly distributed around the zero value. We observe that the EI variability is larger for shallower hypocentral depths than for deeper ones (likely due to a larger crustal heterogeneity at shallower depths).

Energy Index Spatial Variation
To map the spatial variability of EI, we follow the approach of Picozzi et al. (2021). We create a regular three-dimensional grid of size 2 km and select for each grid node the events within a maximum distance from the node equal to one and half times the grid size. Nodes with less than five earthquakes are discarded. Following Nandan et al. (2016), if the ratio of the distance between a grid node and a source (D) with respect to the source event length (L) is more than 10, we can assume that the static stress transfer becomes negligible. We compute L for all events considering the Brune's model (Brune, 1970) and assuming a constant stress drop equal to 1 MPa. For each node, we keep only those events compliant with the above criteria (D/L ≤ 10). To provide an estimate of the uncertainty associated with EI, we follow a bootstrap approach (Efron, 1979), repeating the computation of EI images with 100 random sampling realizations of the original data set with replacement. We then compute the mean of EI, ⟨ ⟩ , and the standard error considering the 100-bootstrap realization. We map ⟨ ⟩ for three distinct periods preceding each of the mainshocks characterizing the 2016 central Italy seismic sequence. Since the period before the 24 August 2016, Mw 6.2 Amatrice earthquake was characterized by a rather low level of background seismicity, we defined the beginning of this period as the beginning of our data set (i.e., 1 January 2005). Similar images for the background data set and for the period following the Mw 6.5 Norcia earthquake (i.e., between the 31 October 2016 and the 31 December 2017) are provided in Figures S5 and S6 of Supporting Information S1, respectively. Figure 3 displays the spatial distribution of ⟨ ⟩ for the years preceding the Mw 6.2 Amatrice earthquake, and it also includes the Mw 6.1, L'Aquila earthquake in the southern part of the investigated areas and its aftershocks. We provide both a map view of ⟨ ⟩ (Figure 3a) and its distribution in depth along a direction corresponding to the mean strike of the main faults in the region (Figure 3c). Figures 3b and 3d shows the standard error associated with each grid node. The spatial distribution of EI for the years preceding the Mw 6.2 Amatrice earthquake shows a few peculiar patterns. We observe a large patch of negative EI values at the Paganica fault. Going northward, the EI progressively increases at the Campotosto faults and further increases, even if with less spatial coherence, at the Laga and Vettore Fault Systems, LFS and VBFS, respectively ( Figure 3a). When we look at the distribution in depth of EI along the mean faults strike (Figure 3c), we observe that at shallow depths the crust including the Paganica and Campotosto faults is characterized by negative EI values. These results seem in good agreement with the pore pressure diffusion process following the Mw 6.1, L'Aquila earthquake and affecting a wide area described by Malagnini et al. (2012). Furthermore, for depths between approximately 10 and 15 km and extending from the base of the Paganica fault to the base of the Laga Fault system (LFS) and the Vettore-Bove Fault system (VBFS) where the Amatrice and the Norcia mainshocks occurred, respectively, EI depicts a seamless wide area of high values. According to these results, the microseismicity suggests a higher stress level along the shear zone (SZ) bounding at depth (i.e., >7 km) the normal faults system where mainshocks occurred. We recall that Sugan et al. (2022) and Vičič et al. (2020) identified this shear zone as the main structure where the progressive unlocking process culminating into the Amatrice mainshock has occurred. Going further north, we also observe large ⟨ ⟩ values below the Vettore-Bove Fault system (VBFS), where the Mw 6.1, Visso and the Mw 6.5, Norcia earthquakes occurred. Figures 3b and 3d shows that the high and low EI main patches are associated to low standard error values. Figure 4 considers the period immediately after the Amatrice earthquake and a few minutes before the Mw 6.1 Visso earthquake. We observe coherent spatial patterns of ⟨ ⟩ in this period too. The south-western part of the LFS is characterized by larger ⟨ ⟩ values with respect to other areas (Figure 4a). We find interesting that the aftershocks occurring in the areas where the Visso and Norcia earthquakes occurred show high ⟨ ⟩ values too,  Table 1 in the main text). (c) Logarithm of radiated energy E S with logarithm of seismic moment M 0 for the whole data set (gray dots), period #1 (red), period #2 (green), and ±1 standard deviation for the estimated parameters (horizontal and vertical gray lines, respectively). Best fit scaling model ±1 standard deviation for period #1 (red lines), period #2 (green lines), and period #1 + #2 (black lines). (d) Energy index (EI) versus logarithm of M 0 (dots colored per hypocentral depth) with associated standard errors for the data belonging to period #1 and period #2. (e) The same as (d) but for the whole data set. (f) Histogram showing the distribution of the EI for period #1 and period #2. (g) The same as (f) but for the whole data set. while the cluster of earthquakes located in the north-eastern sector of the LFS are characterized by low ⟨ ⟩ values. Figure 4c provides further indications about the existence of a stressed patch (i.e., with high EI) at depths between approximately 5 and 10 km at the contour of which the Norcia hypocenter is located. Also in this case, the standard error figures suggest that the main features in the spatial distribution of EI are reliable.
Figures 5a and 5c shows the spatial distribution of EI for the last considered period, which starts just a few minutes after the Visso earthquake and ends right before the Norcia earthquake. The most striking features in these figures, confirmed by the standard error maps (Figures 5b and 5d), are the patch with low ⟨ ⟩ values north of the Visso hypocenter and the large patch with high ⟨ ⟩ values at the VBFS, immediately north-west of the future Norcia earthquake hypocenter. Considering that high EI values can be related to high stress level (Picozzi et al., 2022), this result agrees with the interpretation of Pino et al. (2019), which suggested that the Norcia hypocenter is located at the border of a large patch with a positive Coulomb stress changes after the two previous mainshocks of the sequence.

Temporal Evolution of the Energy Index at Mainshock Hypocenters
An overview of how EI evolves over the whole region and the three periods considered earlier can be seen in Figure 6. In short, during the period preceding the Mw 6.2 Amatrice earthquake (Figure 6a), the seismicity occurring after the Mw 6.1, 2009 L'Aquila earthquake mostly occurs in the northern sector of the Paganica fault (i.e., blue box in Figure 6a), while for the fault southern sector we observe a smaller number of events (i.e., green box). Besides the number of events for fault sectors, we observe that higher EI values are mostly occurring in the region north of the Paganica fault. The progressive increase of EI toward the LFS agrees with the spatial distribution of EI previously discussed (Figure 3). A similar pattern in terms of number of events and EI is also observed for the period lasting from the Mw 6.2 Amatrice and the Mw 6.1 Visso earthquakes (Figure 6b; please note that in this case also, we highlight the fault sector with highest number of seismic events and highest EI by means of a blue box). This result confirms the progressive northward migration of seismicity with higher EI toward the crustal volume including the Mw 6.1 Visso hypocenter. After this last event and before the Mw 6.5 Norcia earthquake, however, the seismicity seems split in two main clusters associated to the northern and southern fault portions. In this case, we observe a progressive southward migration of microseismicity with large EI (blue box in Figure 6c), in agreement with the EI spatial distribution of Figure 5.
Next, we analyze the temporal evolution of EI focusing on cubic crustal volumes with side 5 km centered on the hypocenters of the three mainshocks. We focus on these three events with the aim to verify whether they interacted. As discussed before, for this analysis also, we consider only events with magnitude between Mw 1.8 and Mw 3.5. Figure 7a displays EI for the crustal volume centered on the Amatrice earthquake hypocenter. We observe that right after the Amatrice earthquake, EI shows high and low values rapidly oscillating around zero for a few days (∼4), which highlights the occurrence of events with different dynamic characteristics within the considered crustal volume. However, already after 1 week from the mainshock (i.e., after the 1 September 2016), the seismic activity in this crustal volume becomes modest and the EI values deviate only slightly from zero. This low activity in the crustal volume around the Amatrice hypocenter also continues in coincidence with the Mw 6.1 Visso earthquake. Hence, we can argue that this second mainshock did not significantly alter the stress level around the Amatrice hypocenter (i.e., few events with EI close to zero). Nevertheless, interestingly, we observe a significant change in the EI trend in coincidence with the Mw 6.5 Norcia earthquake. At that time, in fact, we observe a high and protracted seismic activity for the Amatrice crustal volume, with EI values showing even larger amplitude oscillation than during the Amatrice earthquake itself. The crustal volume around the Visso earthquake hypocenter seems not to be affected by the occurrence of the Amatrice earthquake (Figure 7b). Very few sporadic events occur in the period of time between the Amatrice and Visso earthquakes. On the contrary, the seismicity rate is high within this crustal volume in the immediate aftermaths of both Visso and Norcia earthquakes.
Finally, we consider the crustal volume including the Norcia earthquake hypocenter (Figure 7c). Interestingly, already on occasion of the Amatrice earthquake, this volume shows a significant number of events with EI comparable to those occurring in the Amatrice hypocenter volume. Furthermore, we also observe a high number of events in the case of the Visso earthquake and, clearly, in coincidence with the Norcia earthquake itself. To highlight the temporal evolution of the considered crustal volumes, we display in Figure 7d the cumulative of EI in time.
Here, it appears more evident that as the Amatrice earthquake strikes, EI increases within both the Amatrice and Norcia crustal volumes. Moreover, we also observe a similar increase for the Norcia volume on the occasion of the Visso earthquake. This latter earthquake, however, shows a peculiar behavior. Approximately 1 day after the mainshock, the Visso crustal volume is characterized by a decrease in the cumulative EI. This evidence suggests a change in the rupture dynamic for microearthquakes (i.e., larger slip or smaller stress drop). This peculiar behavior also persists after the Norcia earthquake, and only after a few days the cumulative of EI recovers and reaches a rather stable, horizontal trend, which suggests a similar number of events with positive and negative EI values. As previously observed, Figure 7d shows a peculiar increase of EI for the Amatrice volume in coincidence with the Norcia earthquake.

Shear Stress at Mainshock Hypocenters
The radiated energy to seismic moment ratio corresponds to the apparent stress, τ a , over the rigidity modulus, μ (i.e.,  causing earthquake fault slip, , considering data from laboratory stick-slip friction experiments, earthquakes artificially triggered and natural tectonic earthquakes; thus, covering a wide range of seismic moments. The results of McGarr (1999) confirmed the linear relationship between τ a and through the seismic efficiency (i.e., = ) and that varies for crustal earthquakes in a narrow range of values (i.e., from 0.02 to 0.08). Since estimates of are very difficult to obtain because the shear stress is rarely known for earthquakes, seismologists generally estimate the ratio of apparent stress over stress drop, which is also considered a measure of rupture efficiency and known as Savage and Wood efficiency, SW (Savage & Wood, 1971). This is the case, for instance, for central Italy, where Bindi et al. (2020) estimated SW for a set of ∼4,000 earthquakes covering a wide range both in time and M 0 (i.e., earthquakes from 2009 to 2017 and M 0 from 10 11.2 Nm to 10 18.5 Nm). Therefore, first, we calibrate a linear model reconsidering the η and η SW values from McGarr (1999) ( Figure S7a and Table S1 in Supporting Information S1), and then we apply it to estimate η for earthquakes studied by Bindi et al. (2020) and occurred in the same area of our study ( Figure S7b in Supporting Information S1, shows with respect to hypocentral depth and M 0 , and it indicates that primarily scales with M 0 ). We divide the M 0 range in 50 bins equally spaced in logarithmic scale, and we estimate the average for each M 0 -bin and the associated standard deviation. We observe a progressive increase of with M 0 until magnitude Mw ∼3.8 (i.e., which includes the range of magnitude of interest for this work between Mw 1.8 and 3.5), while for larger magnitudes appears more scattered ( Figure S7c in Supporting Information S1). The binned relation between and M 0 obtained considering the data by Bindi et al. (2020) is used to obtain estimates for the earthquakes considered in this study. Then, we assume μ = 30·10 9 Pa to estimate τ a from the E S to M 0 ratios. Therefore, τ a is in turn combined with to obtain estimates of the shear stress (McGarr, 1999) for microearthquakes that have occurred within a limited crustal volume centered on the three mainshocks of the 2016 central Italy seismic sequence. Of course, assuming a uniform μ, while in reality elastic parameters are likely heterogeneous, and estimating from the procedure described above, may lead to large uncertainties. Therefore, we remark that the shear stress values shown in Figure 8 must be considered first order estimates. Nevertheless, both their evolution in time and relative difference when considering different crustal volumes around the mainshock hypocenters provide clues about faults interaction. Daily average of (orange). (b) The same as (a) but for a crustal volume centered on the Mw 6.1 Visso earthquake hypocenter. (c) The same as (a) but for a crustal volume centered on the Mw 6.5 Norcia earthquake hypocenter. Figure 8 shows the temporal evolution of and its daily average for the crustal volumes including the three mainshocks. Of course, both the general trend of each shear stress time series and the interaction among different crustal volumes mimic the trend observed for EI (Figure 7). Here, we find interesting the range of , which varies between 0.15 and 27 MPa for the Amatrice earthquake, 0.18 and 42 MPa for Visso, and finally between 0.12 and 230 MPa for Norcia. We find peculiar that the highest value is observed for the Norcia hypocenter crustal volume but in coincidence of the Visso earthquake, likely resulting from a static stress transfer from Visso earthquake to the already highly stressed patch near the Norcia hypocenter. A further interesting result supporting the idea that the faults interacted with each other is related to the daily average of . We observe that the highest daily average for the Amatrice crustal volume is observed at the time of the Norcia earthquake, while for the crustal volume of this latter event the highest daily average is seen for the Visso earthquake (Figure 8).
The clues for fault interactions are better highlighted considering the daily differences per cubic km (hereinafter, or stress rate, Figure 9). To compute , we use the convex envelope algorithm implemented in MATLAB (i.e., convhull) to estimate the daily volume affected by microseismicity. Since, we wish to associate an estimate of the variability, we apply a Monte Carlo approach. For each earthquake, we randomly and independently sample 5,000 times the associated Gaussian distribution for E S , M 0 , and (i.e., each of the considered parameters are characterized, for each earthquake, by a mean value and an associated standard deviation) and we derive a empirical distribution.
By this set of samples, we obtain 5,000 time series that are in turn compared with the one obtained by using the mean values of E S , M 0 , and for each earthquake (i.e., our best estimate of time series). Following Lomax et al. (2000), the exponential of the likelihood function obtained comparing the sample and the best estimate time series provides the relative probability of the former curves with respect to the latter (Figures 9a-9c). These figures provide a first idea of the variability of when we consider the variability of the input parameters E S , M 0 , and used to derive .
Except for the time series obtained extracting E S , M 0 , and from the tails of their distribution, most of the time series depict a clear narrow trend (i.e., with low variability) around our best estimate . This result supports the idea that our estimates are robust.
Interestingly, the Amatrice crustal volume shows, after the Norcia earthquake, oscillating between ±∼0.5 MPa·day −1 km −3 , which is almost double of the values observed immediately after the Amatrice earthquake itself. The Visso crustal volume shows no activity during and after the Amatrice earthquake and modest oscillations during the Visso and Norcia earthquakes. Finally, the Norcia crustal volume, in agreement with previous observations, sees significant oscillations during the previous mainshocks (i.e., ±∼0.5 MPa·day −1 km −3 ) and a significant one (i.e., +∼1 MPa·day −1 km −3 and −∼0.5 MPa·day −1 km −3 ) for the Norcia earthquake and the following days.
It is also worth noting that for both Amatrice and Norcia earthquakes the oscillations last only a few days after the mainshocks. Figure 9d displays all curves together for facilitating the comparison.
We also repeated the analysis for estimating for a set of control points spread around the three main faults ( Figure S8 in Supporting Information S1). The oscillations for these control points show differences both in their relative maxima amplitude and in the time at which the maximum is reached. Points #1, #2, and #3, which are located at north and east of Visso earthquake, are characterized by oscillations only in coincidence of this earthquake, while for both Amatrice and Norcia only small values are observed ( Figure S9 in Supporting Information S1). Points #4, #5, #7, and #8, which are located to the west of the Vettore-Bove Fault system (VBFS), show high oscillations only on occasion of the Norcia earthquake ( Figure S10 in Supporting Information S1). All points to the south (i.e., #9, #10, #11, and #12) show almost no oscillation, with the interesting exception of a small one that occurred a few days after Norcia earthquake at #12, which is located within the Campotosto fault ( Figure S11 in Supporting Information S1).
These results highlight, in our opinion, a significant activity in the Norcia crustal volume following the Amatrice earthquake and vice versa. EI and derived from microearthquakes provide a confirmation of the interaction between the faults that caused the mainshocks of the 2016 seismic sequence proposed in previous studies (e.g., Improta et al., 2019;Pino et al., 2019).

Discussion and Conclusion
We showed that the Energy Index of microearthquakes is an indicator of crustal stress conditions, and it can help shedding light on the spatiotemporal evolution of stress before and during seismic sequences. EI depicts crustal volumes with anomalously high radiated energy in the proximity of which the Mw 6.2 Amatrice and the Mw 6.5 Norcia earthquakes have nucleated. Focusing on the crustal volumes around the three mainshock hypocenters, we showed that EI clearly indicates interplay between the Mw 6.2 Amatrice and the Mw 6.5 Norcia earthquakes.
Below we summarize the main findings of previous works on the 2016 central Italy seismic sequence. Then, we discuss how our results are related to other studies. Finally, we summarize the distinctive attribute of the microearthquakes EI parameter that led us to believe in its high potential for seismic monitoring.
The 2016 central Italy seismic sequence has been recorded by a dense network of 155 stations, with a mean separation in the epicentral area of 6-8 km (Spallarossa, Cattaneo, et al., 2021). The quality of recordings and the number of earthquakes (i.e., more than 500,000 through mid-2017) made this data set a terrific source of information about the structures and processes driving the seismic sequence. Indeed, even if only a few years have passed, the 2016 central Italy seismic sequence is one of the most studied in the literature (Di Giacomo et al., 2014).
In light of these studies, we wondered, what are the main results that emerged regarding the seismic sequence? According to Chiarabba et al. (2018), the large shocks of the 2016 central Italy seismic sequence have affected adjacent normal faults that were reactivated by high pore pressure at the footwall. High crustal fluid pressure is thought to act as triggering mechanism for the nucleation of large earthquakes in preexisting faults in the Apennines, which in turn often also lead to the activation of other closely spaced fault segments (Di Luccio et al., 2010;Doglioni et al., 2014). In the 2016 seismic sequence, the V P /V S anomalies found by Chiarabba et al. (2018) support the idea of fluid pressurization within the carbonate unit. This process started with the Amatrice earthquake and culminated with the rupture of the main asperity during the Norcia earthquake. Gunatilake and Miller (2022) analyzed the cumulative number of aftershocks of the 2016 central Italy seismic sequence using a nonlinear diffusion model with a source term. These authors provided further evidence about the key role of high-pressure CO 2 in aftershocks occurrence. Furthermore, they suggest that, besides the existence of fluids of deep origin, high-pressure CO 2 might have been generated by thermal decomposition from carbonates. Pino et al. (2019) and Convertito et al. (2020), who studied Coulomb failure function changes, provided evidence of faults interaction. Following the Mw 6.2 Amatrice earthquake, an asperity close to the future Mw 6.5 Norcia hypocenter was affected at its contours by stress corrosion processes, which were likely related to pore pressure increase caused by fluid flow in the upper crust. This stress corrosion mechanism was considered to lead to a clock advance for the Norcia earthquake (Pino et al., 2019).
Besides crustal fluids, the tectonic complexity of the area also played an important role. The Monti Sibillini Thrust (MST) represents a lateral ramp that separates the Laga Fault system (LFS) from the Vettore-Bove Fault system (VBFS) (i.e., the faults where the Amatrice and Norcia earthquakes occurred), and it is a key factor determining the complex geometry of these latter faults. Such geometric complexity of crustal structures has been called into question for explaining the complexity of both coseismic ruptures and aftershocks distribution . The spatiotemporal analysis of Improta et al. (2019) well illustrates the complex interaction occurred among the normal faults where mainshocks took place and the reactivation of the inherited Monti Sibillini Thrust (MST). Vuan et al. (2017) and Vičič et al. (2020) highlighted a further important aspect. These authors hypothesized that the shear zone (SZ) at the base of the seismogenic layer, which is observed throughout the whole central Apennines, might play a crucial role in the progressive unlocking of the overlying normal faults. The SZ is considered by these authors and by Sugan et al. (2022) as one of the main drivers of an 8-year progressive crustal weakening process that culminated with the Amatrice mainshock. According to Sugan et al. (2022), the spatiotemporal evolution of microseismicity in the years before the Amatrice earthquake well agrees with the progressive deformation localization model proposed by Ben-Zion and Zaliapin (2020). Besides the seismic observations, it is also worth mentioning that geodetic transients migrating northward support the existence of an expanding and propagating slow slip pulse in the shear zone along the strike (Vičič et al., 2020). Similar transients associated to slow slip events were also observed in the months preceding the Mw 6.1, 2009 L'Aquila earthquake (Borghi et al., 2016;Sugan et al., 2014).
How do our results relate to this scientific framework of the 2016 central Italy seismic sequence?
To better highlight the spatiotemporal evolution of EI, we prepared time lapse movies of EI values using both a map and a section along strike views (see Movie S1 and Movie S2, respectively). Since the months following the 2009 L'Aquila earthquake, we can clearly see a progressive spread northward of microseismicity characterized by low EI values (i.e., it can be seen both in Movie S1-map and Movie S2-section views and for a depth range between ∼5 and ∼15 km). These results agree with the diffusion of crustal fluids described by Malagnini et al. (2012) and more recently also by Vičič et al. (2020). Indeed, low EI values are consistent with fluids triggering microearthquakes at reduced stress levels. This spreading phase of low EI microseismicity lasts about 5 months and is compatible with the slow slip along strike proposed by Vičič et al. (2020).
The crustal volume surrounding the Amatrice earthquake is affected by high EI values for the first time in 2011, but for a short time only. Then, starting from 2015, we observe a patch located north of the Amatrice hypocenter that is characterized by high EI values. If, again, we take the relation between EI and stress as valid, and consider the relative position of the patch with high EI and the Amatrice hypocenter, our results support, as for the Norcia earthquake, the existence of a weakening stress corrosion process at the patch boundaries that culminated with the Amatrice earthquake.
Immediately after the Amatrice earthquake, the distribution of EI is characterized by high values near the hypocenters of the next future large earthquakes north of Amatrice, while fewer events with low EI are observed southward (Movie S2, section view). This pattern agrees with the along-strike rupture directivity observed by Calderoni et al. (2017).
It is important to highlight that after the Amatrice earthquake the crustal volume surrounding the Norcia hypocenter was characterized by high EI. Interestingly, the activity around the future Norcia hypocenter also persists several days after the Amatrice event, when there is no activity in the Amatrice hypocentral volume. This evidence support the interpretation of the preparatory processes for the Mw 6.5 Norcia earthquake proposed by Pino et al. (2019). The aftershock area of the Amatrice earthquake extended northward considerably, but the seismicity remained substantially confined and persistently active around the Norcia hypocenter (which is compatible with our results shown in Figures S8-S11 of Supporting Information S1). Pino et al. (2019) suggested the existence of an asperity characterized by positive Coulomb stress changes after the previous mainshocks of the sequence, which likely was progressively eroded by pore pressure increase due to fluid flow in the upper crust. This chain of events advanced the Norcia earthquake occurrence.
We observe that the crustal volume around the Visso earthquake shows low activity, whereas the increase of microseismicity around its hypocenter with increasing EI occurs only a few days before the mainshock. After this second mainshock, the crustal volume around its hypocenter seems mostly characterized by moderate and low EI. Likely, these results are related to the existence of crustal fluids that favored larger slip per unit M 0 like in the L'Aquila earthquake.
Our results show that the spatiotemporal evolution of the energy index, EI, starting from the 2009 L'Aquila earthquake to the 2016 seismic sequence well complements the previous information and interpretations, as for instance the results by Kemna et al. (2021) in terms of stress-drop spatial distribution.
To better highlight the usefulness of EI from microearthquakes with respect to the interaction of faults, we studied the temporal evolution of both EI and of the derived first order estimates of the mean loading shear stress ( ) for crustal volumes around the three mainshock hypocenters. Focusing on these specific volumes, we show that the Amatrice and Norcia hypocenter crustal volumes interacted with each other very rapidly. Thanks to EI from microearthquakes, we can thus track the stress increase in the Norcia crustal volume starting already from the Amatrice earthquake. Since high EI values correspond to high energy radiation, our results support the conceptual scheme according to which the crustal volume nearby Norcia earthquake hypocenter was highly stressed and was activated by a static stress transfer and/or dynamic triggering from the Amatrice earthquake. Moreover, we show that the stress for the Norcia hypocenter crustal volume was further enhanced by the Visso earthquake.
What is the potential of the microearthquakes EI parameter for seismic monitoring?
Our results support the McGarr's (1999) vision on the importance of studying the microseismicity radiated energy to seismic moment scaling (which is communicated as scaled energy or apparent stress), as this parameter helps to characterize the crustal stress levels and its spatiotemporal evolution. We think that more efforts from the scientific community are necessary for renewing and increasing the set of seismic efficiency measures for different lithologies and tectonic contexts. In this study, we showed the potential of McGarr's (1999) vision, with the caveat that we assumed that the seismic efficiency values of that study are transferable to the tectonic context of central Italy. Combining the results by Bindi et al. (2020) with that of McGarr (1999), we observed the seismic efficiency of microseismicity scales with seismic moment. The scaling of seismic efficiency with seismic moment is certainly a further interesting line of research, which, if confirmed, would provide important indications about the different physics of small and large ruptures. Picozzi et al. (2022) suggested that by monitoring the deviation of the radiated energy from a reference model calibrated studying the background seismicity (i.e., the energy index, EI), it might be possible to intercept anomalous trends related to the preparatory phase of large earthquakes; or at least, this is what these authors observed for the 2009 L'Aquila earthquake. In this study, we confirmed that by studying the spatiotemporal evolution of microearthquakes EI, it is possible to obtain important information about the evolution of crustal stress and interaction of faults.
Our results renew the idea of Seeber and Armbruster (2000) that microearthquakes can be considered as beacons of stress change in the crust. These authors suggested to focus on Coulomb stress changes. Here, we highlight that EI could complement Coulomb stress changes providing indications about the mechanical state of the crust.
In our opinion, the simplicity of the EI approach, which relies on parameters free from source model assumptions, makes it particularly suitable for both long-term and near real-time monitoring processes leading to fault failure. In the long term, EI can thus be useful to highlight the existence and evolution of progressive deformation localization processes responsible for the nucleation of large earthquakes, as proposed by Ben-Zion and Zaliapin (2020). In near real-time monitoring applications, estimating systematically the radiated energy and the seismic moment for small and moderate magnitude earthquakes, as done for central Italy by the RAMONES service (Spallarossa, Picozzi, et al., 2021), might become a key ingredient for short-term hazard estimates. For instance, the integration of the EI with seismicity rate, b-value, V P /V S , and other parameters in a procedure like the one proposed by Gulia and Wiemer (2019) might help improving the real-time discrimination of foreshocks and aftershocks, and therefore gather hints regarding the preparation phase of large earthquakes.
As discussed by Vičič et al. (2020), among the faults that were active in central Italy from the L'Aquila 2009 to the 2016 central Italy sequence, the Campotosto fault is considered still mostly unbroken and thus capable of producing large earthquakes (e.g., among many, Cheloni et al., 2019;Chiaraluce et al., 2011;Valoroso et al., 2013). In agreement with these latter studies, Figure S6 in Supporting Information S1 shows that the region between the Laga Fault system and the Campotosto fault presents higher EI values than before the 2016 central Italy sequence. Therefore, future efforts will be directed to a long-term testing of our EI approach on specific areas, as for instance the Campotosto fault.
Our study focused on central Italy, but our approach is applicable to any tectonic context. To this aim, it is important to export strategies to estimate E S and M 0 at local scales for moderate and small magnitude earthquakes.
For Italy, we have shown that the RAMONES procedure allows us to collect useful information, and we believe that in the near future we will be able to extend the EI monitoring also to other areas. However, it is worth noting that for extending the EI monitoring to other regions, it is necessary to calibrate ad hoc attenuation relationships for the considered area. The large amount of data available for Italy makes us confident about the possibility of extending the EI to the whole country.