Spatio‐Temporal Evolution of Earthquake Static Stress Drop Values in the 2016–2017 Central Italy Seismic Sequence

The static stress drop of an earthquake is an indicator of the stress state of a specific fault before rupture initiation. The stress state is primarily controlled by the ambient stress field, fault strength, fault complexity, and the presence of fluids. This study aims to investigate the spatio‐temporal distribution of static stress drop values of the 2016–2017 multi‐fault rupture seismic sequence in central Italy, which includes three earthquakes with Mw ≥ 5.9 (Amatrice, Visso, and Norcia earthquakes), and over 95,000 aftershocks (M 0.5–6.5). We estimate stress drop values using a circular crack model with corner frequency and seismic moment estimates from single‐spectra fitting, a cluster‐event method, and spectral‐ratio fitting. The temporal distribution of stress drop values shows an apparent increase of stress drop following a large earthquake (Mw ≥ 5.9). The spatial distribution shows comparably high stress drop values for early aftershocks surrounding the mainshock rupture area. High stress drop events correlate with fault complexity, such as fault intersections at depth and reactivated thrust fronts. We observe a constant stress drop for Mw ≥ ∼3, in contrast to previous studies. Instrument response and signal‐to‐noise bandwidth limitations likely govern the observed decrease in stress drop with decreasing magnitude for events with Mw ≤ 3. The spatio‐temporal distribution of stress drop values in a complex seismic sequence could support a more complete understanding of the earthquake rupture process and the evolution of seismic sequences. It could also highlight areas where stress loading is focused, which would have implications for short and intermediate term seismic hazard estimates.

rupture area (e.g., Aki, 1967;Abercrombie, 1995;Prieto et al., 2004). Stress drop -magnitude independence, or self-similar scaling, follows from theoretical considerations (Aki, 1967;Haskell, 1966) and has long been supported by observational considerations, at least for moderate to large earthquakes (Kanamori & Anderson, 1975). However, contrasting research shows the breakdown of self-similar scaling on both the regional and global scale, in particular for smaller-magnitude events (e.g., Bindi et al., 2020;Kanamori et al., 1993;Mayeda et al., 2007;Oth et al., 2011;Trugman et al., 2017). The discrepancy between stress drop estimates obtained with different methods complicates the controversy even further (e.g., Abercrombie, 2014;Abercrombie, 2015;Kaneko & Shearer, 2015;Neely et al., 2020;Shearer et al., 2019). Self-similarity over a given magnitude range would imply that earthquake static stress drop may function as a proxy parameter for physical properties of the fault plane. The properties include, for example, the rupture speed, the dynamic stress level, and stress loading conditions. Estimating stress drop values accurately is important for interpreting real variations in the context of faulting. For example, global-scale studies suggest that thrust-style earthquakes at subduction interfaces have lower stress drop values compared to normal faulting or strike-slip events in intraplate regimes (e.g., Allmann & Shearer, 2009). Many interpret the lower observed values of thrust-style earthquakes to result from the fact that they predominantly occur in the shallow part of the subducting slab with lower rigidity. Stress drop values also seem to vary significantly on the regional scale. For example, stress drop is observed to vary both spatially in southern California (e.g., Shearer et al., 2006), and temporally for repeating earthquakes on the San Andreas Fault at Parkfield in central California (e.g., Abercrombie, 2014). Stress drop values are also commonly found to increase with depth (Goebel et al., 2015;Hardebeck & Aron, 2009), but depth dependence has been interpreted as an observational artifact (Abercrombie et al., 2020), to vary significantly along strike (Moyer et al., 2018;Shearer et al., 2006;Trugman, 2020), and to depend on fault rheological and geometrical complexities (Harrington & Brodsky, 2009;Kirkpatrick et al., 2020). Some work has shown that stress drop values increase after large mainshocks (Chen & Shearer, 2013;Trugman, 2020), where high stress drop events surround the area of maximum slip during a large earthquake (Yamada et al., 2010). Stress drop values are also observed to be lower in regions where fluid pressures and heat flow is high, and where associated fault strength is likely low (e.g., Chen & Shearer, 2011;Hough & Kanamori, 2002;Oth, 2013;Yamada et al., 2015). Conversely, relatively high stress drop estimates correlate with relatively high fault strengths due to inferred low seismic coupling (e.g., Moyer et al., 2018) or high fault heterogeneity (e.g., Goebel et al., 2015).
The spatial and temporal variability in earthquake stress drop estimates and the correlation with certain types of faulting conditions suggest it is a parameter that could be used as a proxy to describe the mechanics and strength of faults. However, a significant body of work shows that stress drop estimates can be highly biased by different methodologies and the limitations of available seismic data, particularly for events smaller than M ∼4 (e.g., Abercrombie, 2021). We focus here on estimating stress drop values of the 2016-2017 Amatrice-Visso-Norcia (AVN) earthquake sequence in central Italy, where a large number of events were recorded with dense station coverage. We use different methodologies to investigate stress drop variation over time, space, and event magnitude using high-quality data (e.g., high signal-to-noise waveform recordings and precise catalog locations) in a region hosting complex rupture behavior with several normal faults, intersecting antithetic faults, and reactivated pre-existing thrust fronts (e.g., Chiaraluce et al., 2017;Improta et al., 2019;Michele et al., 2020). In addition, the AVN seismic sequence activated two major fault systems with different structural inheritance (e.g., Barchi et al., 2021;Porreca et al., 2018;Sebastiani et al., 2019). Therefore the data and seismicity distribution provide an optimal data setting to explore the relation between fault mechanics, geological setting, and stress drop variation. We apply methods previously unused in the study area to broaden the availability of stress drop estimates in the region, validate previous interpretations, and use the spatial and temporal variation in stress drops as a proxy for the seismic response of the fault systems.
The study area is located in the central and northern Italian Apennines and hosts complex patterns of seismicity as a result of a poly-phased tectonic history. In the Neogene, the Apennines were situated in a compressional regime where mainly Miocene flysch deposits were overthrusted by Meso-Cenozoic carbonate rocks (e.g., Lavecchia et al., 1994). The north-northeast-trending Olevano-Antrodoco-Sibillini (OAS) tectonic alignment is one example of an inherited thrust front (Centamore & Rossi, 2009). Because of the movement of the thrust front towards the Adriatic Sea (black barbed line in the inset of Figure 1), the Apennines transformed into a back-arc extensional regime in the Pleistocene (Carminati & Doglioni, 2012). The extension resulted in a complex system of mainly, but not exclusively, NE-SW trending active normal faults (e.g., Pizzi et al., 2009). Pre-existing thrust KEMNA ET AL.

10.1029/2021JB022566
3 of 29 fronts play a complex role in the current faulting system by exhibiting a range of mechanical behavior. They either limit the extension of normal fault rupture segments with observed segment lengths of up to 20 km (e.g., Pizzi et al., 2009), or are reactivated as normal faults with mostly optimal orientation in the present-day NW-SE trending extensional regime (e.g., Chiaraluce et al., 2005;Improta et al., 2019). They also control the depth extent of the seismogenic layer by acting as subhorizontal shear zones (e.g., Chiaraluce et al., 2017;Michele et al., 2020;Vuan et al., 2017). The cumulative regional tectonic extension ranges from 1-3 mm/year (Carafa & Bird, 2016) and has resulted in several multi-fault ruptures in recent times (e.g., MX Irpinia-Basilicata 1980 (Bernard & Zollo, 1989), MIX Colfiorito 1997 (Chiaraluce et al., 2005), and MIX L'Aquila 2009 (Valoroso et al., 2013).  Villani et al., 2018) in the main map and the current approximate location of the thrust front in the top right overview map. The arrows in the overview map outline the direction and amount of extension in central Italy. Blue, yellow, and red hexagons mark the Amatrice, Visso, and Norcia earthquake epicentral locations with associated time domain moment tensors (Scognamiglio et al., 2006). Brown points and diamonds denote all earthquake epicenters and earthquakes with M w ≥ 4 used in this study, respectively. Blue triangles denote a subset of stations operated by the INGV National Seismic Network and real-time temporary stations. The DEM is taken from Tarquini et al. (2012).

of 29
The most recent AVN seismic sequence activated several tectonic structures in the central Apennines. It initiated on 24 August 2016 with the M w 6.0 Amatrice earthquake and activated an 80-km-long NNW-SSE trending system of faults that resulted in several surface ruptures. The sequence produced a substantial rise in seismicity, which reached its peak after the M w 5.9 Visso and M w 6.5 Norcia earthquakes on 26 October and 30 October Tinti et al., 2016). The three largest earthquakes show clear NNW-SSE trending and SW-dipping normal faulting mechanisms along segments of the Laga fault system (LFS) in the south and of the Vettore-Bove fault system (VBFS) in the north (Cheloni et al., 2017;Lavecchia et al., 2016;Pizzi et al., 2017, Figure 1). The AVN seismic sequence further ruptured several antithetic NE-dipping normal faults intersecting the VBFS at depth and reactivated optimally oriented segments of the OAS tectonic alignment, such as the Mount Sibillini Thrust (MST) Improta et al., 2019;Scognamiglio et al., 2018;Michele et al., 2020).
Here we analyze over 95,000 earthquake source parameters in the time window between 1 January 2016 to 17 January 2018. The time window includes nearly 3 months of background and foreshock seismicity, the three largest earthquakes with their specific aftershock sequences, and over a year of aftershock activity following the Norcia earthquake. We estimate static stress drop values of all earthquakes using three methods: single spectra, cluster-event, and spectral-ratio methods. We then investigate and interpret the spatial distribution of stress drop values with fault complexity inferred from high-quality aftershock locations, surface rupture patterns, slip models, and previous interpretations of the complex rupture behavior. In addition, we analyze the spatial and temporal distribution of stress drop values and the correlation with earthquake magnitude. We will show in the following that stress drop variation in space appears to be governed by fault complexity and that it shows a clear correlation with geological conditions. That variation in space likely corresponds to temporal fluctuations related to stress unloading during the sequence.

Data and Spectral Analysis
We take advantage of the large body of research investigating the AVN seismic sequence by combining catalog hypocenter information from several previous studies to build the groundwork of our analysis. We use a combination of four earthquake catalogs, including one with double-difference relocations, to calculate stress drop values of two years of seismicity around the AVN seismic sequence. We begin by estimating the source spectral parameters for a specific earthquake using the single spectra, cluster-event, and spectral-ratio methods with data from all available stations with manual phase arrival picks from National Institute of Geophysics and Volcanology (INGV) described in the following section. In the following subsections, we first describe the collection of earthquake catalogs, waveform recordings, and focal mechanisms. Subsequent sections give an overview of the spectral analysis with details of each method used to estimate a stress drop value for a given earthquake.

Earthquake Catalog and Waveform Data Collection
We analyze approximately 95,000 earthquakes between 1 January 2016 and 17 January 2018. Earthquakes detected before the 2016-2017 Central Italy seismic sequence are taken from the Italian Seismic Bulletin (BSI; http://terremoti.ingv.it/bsi). Following the initiation of the AVN seismic sequence on 24 August with the M w 6.0 Amatrice earthquake, all subsequent earthquakes are then taken from a combination of catalogs, including the BSI, Chiaraluce et al. (2017), referred to as C2017, Improta et al. (2019), referred to as I2019, and Michele et al. (2020), referred to as M2020. We rank locations following the order of M2020, I2019, C2017, BSI for events common to multiple catalogs ( Figure S1 in Supporting Information S1), where priority is based on the reported location uncertainties. We further associate focal mechanisms to events with M w ≥ 3 by merging available focal mechanism solutions from INGV Time-Domain Moment Tensors (TDMT, Scognamiglio et al., 2006, http://cnt.rm.ingv.it/tdmt) and the Moment Tensor Solutions for Italy (Herrmann et al., 2011, http://www.eas.slu. edu/eqc/eqc_mt/MECH.IT). To ensure the quality of the focal mechanism solutions, we take only A and B quality solutions (INGV, see Scognamiglio et al., 2006 for details) and solutions with a goodness-of-fit 0.5 (MECH, see Herrmann et al., 2011 for details), which results in 424 events with associated focal mechanism solutions. For focal mechanism solutions contained in both catalogs, we take the TDMT solution.
We use waveforms from a total of 289 stations from 13 seismic networks in the spectral analysis. We use only three-component stations with a flat displacement instrument response in a magnitude-dependent frequency range (see Supporting Information S1 for further details). All waveforms are publicly accessible from the ( ) = ( ) * ( ) * ( ). (1) The static stress drop released by faulting is a function of the seismic moment and the spectral corner frequency, and is governed solely by the source time function, e(t) in Equation 1. The source time function is commonly described in the frequency domain e(f) with the following analytical expression: where the long-period spectral amplitude Ω 0 is proportional to the seismic moment (or event size), and the corner frequency f c is inversely proportional to the event duration. Two additional parameters influence the shape of the earthquake spectrum, the frequency fall-off rate n, and a spectral shape constant γ. Both parameters are typically empirically related to attenuation conditions expressed by the G(t) term in Equation 1 (e.g., Abercrombie, 1995). The spectral shape constant, which controls the sharpness of the spectral corner, is usually set to γ = 1 or γ = 2, which corresponds to the source models of Brune (1970), referred to as Brune-model and Boatwright (1980), referred to as Boatwright-model, respectively. The path attenuation effects in the frequency domain G(t) can be described with where t is the travel time and Q the frequency-independent seismic quality factor (Abercrombie, 1995). Once the instrument response term I(t) has been removed from s(t) for a specific station, the e(f) in Equation 2 can be modified with Equation 3 to describe the theoretical displacement spectra Ω t (f) of a direct wave (P-or S-wave) in the frequency domain The parameter Ω 0 is proportional to the earthquake size, M 0 . Assuming that the rupture velocity is constant, then f c is also inversely related to the rupture dimension. Once the parameters M 0 and f c are obtained through spectral fitting, they can then be paired with a given source model to estimate the static stress drop Δσ for an earthquake.
The static stress drop describes the average stress released during the faulting process. It can be calculated from the average slip D(M 0 ) and the general rupture dimension L, the latter of which is typically inferred from the corner frequency and assumed rupture velocity (Kanamori & Brodsky, 2004). We calculate the rupture dimension and stress drop from the M 0 and f c assuming a circular crack model with radius r (Eshelby, 1957): The variable β is the shear wave velocity, and k is a source-model-dependent constant. Here, we use k = 0.35 for P and k = 0.26 for S-waves assuming a symmetrical circular model with a rupture velocity of 0.8 β (Kaneko & KEMNA ET AL. 10.1029/2021JB022566 6 of 29 . The depth-dependent value of β comes from the 1D gradient-velocity model from Carannante et al. (2013) at the hypocentral depth to avoid an artificial inferred increase in stress drop value with depth (e.g., Allmann & Shearer, 2007).

Spectral Estimation
To estimate f c , M 0 , and subsequently Δσ, we use the displacement amplitude spectra for each available event waveform from the direct wave (P-or S-wave). In contrast, previous studies used coda wave spectra to estimate source parameters both generally in the Apennines (e.g., Malagnini et al., 2010) and specifically for the AVN sequence Morasca et al., 2019). One useful secondary objective of this study will be to explore the common variation in Δσ estimates that result across different methodologies, which is another reason for using direct waves. We estimate direct wave spectra from time windows cut around the direct phase arrivals (P-or S-wave), estimate the amplitude spectra, and use three fitting methods. The fitting methods include single-spectra fitting, cluster-event method, and spectral-ratio fitting, and are described in the subsections that follow.
The 1.4 million phase picks used for the spectral estimation correspond to a magnitude range of 0.5 ≤ M ≤ 6.5. The window length around the phase arrival can play a crucial role in biasing the corner frequency estimate if it is not appropriate for the event magnitude. For example, a time window that is too short could limit the ability to resolve low corner frequencies and bias estimates toward the reciprocal of the window length. On the other hand, a time window that is too long window could include long-period noise or secondary phase arrivals and could result in low signal-to-noise ratios. One way to avoid such pitfalls is to use a magnitude-based time window length Ruhl et al., 2017).
Selecting appropriate magnitude-based time-window lengths rests upon observations that most stress drop values typically range between 0.1 and 100 MPa (Hanks, 1977). Constant stress drop values (within approximately three orders of magnitude) coupled with theoretical considerations imply that M w and f c follow a predictable scaling relationship that can be used for the window length selection. To avoid potential bias, we use a combined approach that first restricts the minimum and maximum bounds on window length using theoretical stress drop values. The window length is then selected within the bounds based on the radiated energy of each event on a given station. To first select the minimum-maximum bounds, we follow (Ross & Ben-Zion, 2016) and calculate minimum and maximum theoretical stress drop values and their corresponding corner frequencies (Equation 6). We define the upper and lower bounds by the inverse of the corner frequency values corresponding stress drop values of 0.1 and 10 MPa multiplied by a factor of five (in order to avoid being too pre-emptively restrictive). Furthermore, we require a minimum window length of 1 s, and a maximum window length of 30 s for P-waves (governed by the S-P time for the largest events recorded at the most distant stations) and 100 s for S-waves. We also remove P-wave spectra where the P-wave window length is 150% S-P time, or where the S-wave arrival is missing, in order to avoid S-wave energy contamination of P-wave spectra. We then window the waveforms for a specific station using 150% of the maximum window length bound ( Figure S3 in Supporting Information S1) and calculate the cumulative energy as the integral of the velocity waveform starting 0.2 s before the reported arrival time. A window length of 150% of the S-P time could still contain some S-wave energy for the closest stations. However, the tapering performed in the spectral estimation and placing the start of the window 0.2 s before the P-phase maximizes the window length and makes the contribution of S-wave energy negligible for the small percentage of stations affected. We use the vertical and the root sum squared (RSS) of all components for P-and S-waves, respectively. We then set the window length for a given event on a specific station corresponding to the time window containing 90% of the energy on stations recorded within a hypocentral distance of 25 km; time window corresponding to 80% of the energy for distances between 25 and 50 km; and windows corresponding to 70% of the energy for greater distances (Bindi et al., 2018;Pacor et al., 2016). If the window is smaller or larger than the minimum (1 s) or maximum (30/100 s) window length, respectively, the window is adjusted accordingly. We then estimate the displacement spectral amplitude from the windowed waveforms starting 0.2 s before the phase pick using the multi-taper approach implemented by Prieto et al. (2009).
Each spectrum has to pass a signal-to-noise (SNR) criteria before using the spectrum in the different spectral fitting procedures detailed in the subsequent subsections. We estimate a noise spectrum for each estimated signal spectrum using the same window length as the signal spectrum ending 0.4 s before the reported P-wave arrival time (dashed lines Figure S4a in Supporting Information S1). Noise spectra for cases where no P-wave arrival is available are derived from noise windows ending 10 s before the respective S-wave arrival time. We define a magnitude-dependent criteria to define the frequency band of fitting where the SNR must be ≥3, similar to the 7 of 29 time-window length estimate. The corner frequencies corresponding to theoretical stress drop values of 0.01, 1,000 MPa correspond to the lower and upper frequency limits of frequency fitting band. In addition, for lower-magnitude earthquakes (M ∼1.5) we impose limits of 4 and 40 Hz for the frequency corners of the fitting band (gray-shaded areas in Figure 3). The imposed limits are based on the instrument response functions, where all instruments have a flat response function between ∼4 and ∼40 Hz. Each event-station combination noise spectrum is combined across channels as for the signal and then validated against the SNR criteria SNR ≥3 in the respective frequency band.

Single-Spectrum Fitting
As a first-order approach to fitting source spectra to recover M 0 and f c , we fit individual spectra for each individual station-event-phase combination passing the SNR criteria. The single-spectra fitting method is based on Equation 4, which we use to fit each individual displacement spectrum with instrument response removed (we refer the reader to the Supporting Information S1 for details on the instrument response removal). We perform each event-station-phase fit individually, and average the results to obtain a single Ω 0 and f c estimate per event ( Figure 2). Several parameters influence the fit, which can be held fixed or allowed to vary during the fitting process. In particular, there is an inherent tradeoff between Q and f c during the fitting procedure. Consequently, several studies have found that attenuation and site effects likely bias the estimation of f c to lower values during single-spectra fitting. Site effects are not considered by the single-spectra fitting approach. We, therefore use it as a first-order approach to be refined by the two method descriptions that follow (CEM and spectral-ratio fitting), where data quality permits (e.g., Abercrombie, 2015;Huang et al., 2016;Shearer et al., 2019).

We fit Equation 4 to each individual spectrum using a trust-region-reflective least-squares minimization algorithm implemented by Newville et al. (2014) inside a water-level defined frequency band (Figures S4b and S4c
in Supporting Information S1). As we are primarily interested in the frequency band surrounding the corner frequency and spectra are typically more unstable at higher frequencies, we define a dynamic spectra-based fitting band. We use the lower SNR frequency limit and define the higher frequency fitting limit based on a water-level of 40 (reduction from the maximum displacement spectral amplitude in dB). If the upper frequency fitting limit is beyond the frequency where the SNR criterion is satisfied, we only fit the band up to the limit where SNR ≥3. The fitting band is therefore dependent on the respective station-event spectrum. We interpolate the displacement spectra amplitude to 100 equally spaced points in log-space inside the fitting band to further decrease the influence of high-frequency instabilities.
We use the Boatwright model (γ = 2), which gives consistently lower χ 2 statistics compared to the Brune model (γ = 1) and fix the long-period spectrum Ω 0 to the median of the highest 10% amplitude values in a frequency range defined by the lower and 1/4 of the upper SNR band frequencies. We let the parameters Q (100-1,000) and n (2-4) vary within the specified ranges due to the large number of events. The large ranges also consider the complex geological and tectonic environment, where we expect substantial heterogeneity in ray-path attenuation.
We then use the individual station-event-phase spectrum estimates of Ω 0 to estimate the seismic moment M 0 . We use the estimated Ω 0 at each component (Z, E, N) of a three-component seismometer with the following relation: where ρ and c are the density and wave velocity at the earthquake hypocenter, R is the hypocentral distance, and U ϕθ is the mean radiation pattern (0.52 and 0.63 for P-and S-waves, respectively Madariaga, 1976). We then convert M 0 to M w following Kanamori (1977): We use only the vertical component for P-wave and the root sum squared (RSS) of all three components for S-wave spectra, which yields the highest number of estimates ( Figure S5 in Supporting Information S1). Nevertheless, the P-wave estimates show similar values compared to the S-wave estimates. As the S-wave estimates are also more numerous and more stable across stations, we use only the S-wave estimates in the two other approaches to refine corner frequency estimates ( Figure S5 in Supporting Information S1). We estimate mean f c , M 0 , and Δσ (Equation 6) values for each event. The 95% confidence intervals result from using the delete-one station KEMNA ET AL.

10.1029/2021JB022566
8 of 29 jackknife-mean in the logarithmic domain of all estimates for a specific event and phase (e.g., Prieto et al., 2006). All reported values result from events with a minimum of 10 individual station estimates and a maximum azimuthal gap of 180°. The presentation of the results refers to single spectra S-wave averaged estimates with the acronym SSSA ( Figure 2). The left column in Figure 2 summarizes the parameter fits (f c , Q, and n) as well as the calculated source parameters based on the fits (M 0 , and Δσ) as detailed above.

Cluster-Event Method
One approach to refining f c estimates while still retaining a large number of events for analysis is to use clusters of similar events with similar source-station paths to reduce the number of free parameters during fitting. As shown by several studies (e.g., Ko et al., 2012;Stachnik et al., 2004;Sonley & Abercrombie, 2006), individual f c estimates from single-spectra fitting are commonly underestimated due to the trade-off between f c , Q, and, if . To estimate path attenuation effects more rigorously and to obtain a more robust f c estimate, we follow the cluster-event method (CEM) proposed by Ko et al. (2012).
The CEM assumes that earthquakes in close spatial proximity (referred to as clusters) share similar ray paths and, as a result, identical integrated Q-values along those paths ( Figure 2). Based on common assumed Q and n values for individual cluster-station pairs, all individual events can be fitted with a single respective f c and Ω 0 value common to all stations to reduce the degrees of freedom. Assuming K events recorded at J stations, the CEM decreases the number of degrees of freedom to K + J from (K ⋅ J) ⋅ 2 for single spectrum fitting. The assumption of identical Q-values depends primarily on the spatial extension of the cluster (Ko et al., 2012).
The present complex geological and tectonic setting of the AVN seismic sequence likely causes highly variable attenuation paths for distinct source-station pairs. The setting and the large number of earthquakes that potentially generate dense clusters make this data set ideally suited to using a CEM approach. Nevertheless, the seismicity in the AVN seismic sequence is mainly clustered around elongated tectonic structures, and singular clusters are challenging to define. We therefore use all events with a single spectrum fitting estimate to define a cluster containing the closest (hypocentral) neighbors of a given event (maximum distance 4 km). We require at least 10 neighbors with a maximum of 40 events to define a given cluster. One drawback of the approach rests on the assumption of a single corner frequency across several stations, which neglects the presence of high-frequency directivity (which is observed in the study region Pacor et al., 2016) and radiation pattern effects, where the f c is Figure 3. Estimated moment magnitude vs. corner frequency vs. for all events from S-wave spectra. All abbreviations are defined in Figure 2 and Section 2.2. The numbers after each abbreviation in the figure legend denote the number of estimates for each method. The SSSA ( Figure S9 in Supporting Information S1) and CEMA estimates are colored by the density of points around a given estimate (see colorbar taken from the legend). Reviewed spectral ratio estimates for eGf events (RREA, turquoise-colored rectangles) that lie beyond the bandwidth resolution limit of the seismic instruments used in this study (∼40 Hz), are fixed to a corner frequency f c = 40 Hz and marked with a rightward-pointing arrow ( Figure S14 in Supporting Information S1). Dashed diagonal lines denote stress drop isolines for a shear-wave velocity of 3,200 m/s. The lower-upper Signal-to-Noise (SNR) boundary is based on a theoretical stress drop value of 0.01-1,000 MPa with minimum and maximum corner frequency values of 4 and 40 Hz, respectively (Section 2.2.1).
10 of 29 predicted to vary over the focal sphere (Madariaga, 1976). However, Ko et al. (2012) showed that potential biases due to radiation and directivity effects could be neglected when using a dense seismic array. To ensure the density of the estimates for a given cluster, we solve only for clusters that have at least 10 stations with at least five overlapping station-spectra pairs across the events in the cluster. In addition, we require clusters to have a maximum azimuthal gap of 180° between stations.
To solve individual clusters with K events recorded at J stations to estimate K ⋅ f c values, we follow the approach described in Ko et al. (2012) that formulates a misfit function Θ in terms of the L 2 norm: D kj , S kj are the observed spectra and calculated spectral fits (in our case Equation 4) for the k th event at the j th station, respectively. The weighting factor Ψ k is event-dependent, in contrast to the frequency range-dependent weighting factor used by Ko et al. (2012). It normalizes the individual event-misfit to the number of stations with a spectral fit before summation. We use the weighting factor due to the fact that each event theoretically has a different number of stations that meet the SNR criteria. Station-event spectra not passing the SNR criteria are treated as zero in the summation, and do not influence the event misfit. We use the spectra from the single spectrum fitting and the free parameters f c (0.01-40; for each event), Q (100-2000; for each station), n (2-4; for each cluster), where γ is fixed to 2 as in the single spectrum fitting. The middle column in Figure 2 summarizes the parameter fits for the CEM approach. Because the f c value is held fixed for a given event within a cluster, it can be considered as an average of the variation that would be expected from station observations at a range of azimuths.
As the inverse problem is multi-dimensional ( Figure S6 in Supporting Information S1) and highly non-linear, several local minima are expected (Ko et al., 2012). Here, we use a Self-adaptive Differential Evolution algorithm (Biscani & Izzo, 2020 to solve each cluster for f ck , Q j , and n (see Supporting Information S1 for more details). Each cluster produces K ⋅ f c , J ⋅ Q values, and one n value. We use the M 0 estimates from single spectra, as they are not typically biased.
After solving each individual cluster, we use all f c estimates for a given event from different clusters and use the delete-one jackknife-mean in the logarithmic domain (as done for the single spectrum fitting) to estimate mean values and confidence intervals. We only report values with an f c error ≤10%. Cluster averaged cluster-event method estimates are referred to as CEMA ( Figure 2). We then use the CEMA f c and with the SSSA estimated M 0 values to calculate the CEMA Δσ values.

Spectral-Ratio Fitting
As mentioned before, earthquake spectra can be affected significantly by non-source-related terms, where the influence may be reduced using the CEM as opposed to single-spectra fitting. Nevertheless, both methods are expected to be biased by non-source-related terms, particularly for smaller magnitude events (M ≤ 3) as the SNR decreases and the expected frequency range of the source increases. To get a better constraint on the f c values and estimate uncertainties in the single-spectra fitting and CEM results, we use a spectral ratio approach (Bakun & Bufe, 1975).
In general, spectral ratio estimates allow better constraints on f c values by minimizing the imprint of non-source related terms on the spectral shape. The approach is based on the fact that the G(t) and I(t) terms in Equation 1 cancel in the ratio of two event spectra with a common source location recorded at a given station (i.e., with identical travel paths), leaving only relative source information between the two events (Hartzell, 1978;Hough, 1997). In general, the approach uses co-located event pairs with a sufficient magnitude difference to resolve their respective corner frequency values, or the corner frequency value of the large event if the smaller event corner frequency lies outside the bandwidth of resolution. The event similarity, magnitude difference, and SNR requirements over a broad frequency band generally limit the number of events available for analysis. However, the resulting refined f c estimates are arguably subject to less bias compared to single spectrum and CEM approaches (assuming station coverage is comparable).
We calculate spectral ratios between pairs of highly similar events (Figure 2) selected using a full-waveform cross-correlation criterion (we refer the reader to the supplemental material for details on event pair selection).
10.1029/2021JB022566 11 of 29 Dividing the two event spectra cancels the non-source related terms, and the remaining displacement spectral ratio Ω r (f) is then written as: where f c1 and f c2 are the corner frequencies of the larger magnitude target and the smaller magnitude empirical Green's function events (eGf), respectively. As noted above, if the corner frequency of the eGf is outside the bandwidth of resolution (defined by SNR ≥3), then the smaller event source time function is effectively considered as a delta function and only the source time function of the target event remains in the ratio. Alternatively, suppose one is interested in constraining the f c value of the target event only. In that case, the spectral ratio can be fit up to frequency values before the inflection point related to the f c of the eGf. In such cases, one may consider only the influence of the target event on the spectral ratio, in which case one would fit the target event source in the fitting frequency band: The parameters γ and n are the same as described above. Similarly, the γ value is fixed during the fitting process and n is fit along with the corner frequency values, as summarized in the right column of Figure 2. We detail below how we first use Equation 11 to obtain stable f c estimates of the target events in an automated approach, and then use a semi-automated approach to obtain f c estimates for eGfs where bandwidth limitations permit.
We calculate the spectral ratio from the spectra of the target and eGf using the procedure described for the single-spectra fitting without removing the instrument response. Instead, we apply a high-pass filter with corner at 0.5 ⋅ f c , where f c is the CEMA estimate of a given target event, if available, or with the SSSA estimate if not. The high-pass filter reduces the impact of long-period noise on the ratios. Typically the long-period noise is Gaussian distributed and would cancel between the target and eGf event. However, as we use comparably small time windows to enhance the signal around the phase arrival, we observe a non-Gaussian, non-canceling long-period noise distribution in particular windows. Manual revision of spectral ratio pairs shows that the high-pass filter stabilizes the spectral ratio, and helps preserve the SNR requirements. We also use the same window length for both the target and eGf events to ensure the same frequency resolution. After applying the same SNR criteria as in the signal spectrum fitting, we cut all station ratios of a specific event pair to the length of the shortest frequency vector across station spectra. We then divide the spectra at each station to get station ratios and stack all available station ratios for a given target-eGf pair using the geometric mean. We then use an automated approach to estimate target event f c values and a semi-automated approach to estimate spectral ratio values. In both approaches, we use a variable n (2-4) and fixed γ (2; Boatwright-model). We further require an amplitude ratio ≥2 of the stacked spectral ratio.
In the automated workflow, we fit only the source model (Equation 11) to estimate f c of the target event (f c1 in the following). The benefit of the source model lies in its independence of the visibility of f c2 , which could be beyond the upper frequency limit. Therefore, the source model could be used in an automatic way to obtain a higher quantity of estimates compared to a workflow that needs manual revision of the spectral ratios. We consider all pairs with a cross-correlation value CC ≥ 0.7 and require at least five station ratios to pass the SNR criteria for a particular event pair (we refer the reader to the supplemental material for details on event pair selection). We then fit Equation 11 to the station-stacked ratio using a water-level of 5 (similar to single-spectra fitting, Figure  S7 in Supporting Information S1) to avoid fitting the high-frequency asymptote above eGf corner frequency. We interpolate the spectral ratio to 100 equally spaced points in logarithmic space to further reduce the impact of high-frequency uncertainties (red points in Figure S7 in Supporting Information S1). We fix Ω 0r to the mean of the 10% largest values and use the same algorithm as for single spectra to obtain f c1 . In addition, we consider only target events with a minimum of two individual eGf events. Individual f c1 estimates and uncertainties result from the average of all estimates of a given target event using the delete-one jackknife mean in the logarithmic domain (as done for the single spectrum fitting). We report only target events with f c1 errors ≤20% to ensure that the f c1 estimate is consistent across different eGf events. We refer to the f c1 estimates from averaged target events using the automated source model fitting workflow as RSTA (Figure 2).

of 29
The semi-automatic workflow then further evaluates event pairs with stable target f c1 estimates to obtain estimates of both the target (refined) and the eGf (f c2 ) events. Because the eGf events are smaller in magnitude and subject to frequency bandwidth limitations, we impose stricter criteria for the procedure. We first pre-calculate the spectral ratios for pairs with CC ≥ 0.8 and a minimum of 10 station ratios and manually review the spectral ratio fits. We perform the fits in the same way as described for the source model fitting described above, but the SNR frequency band defines the upper fit band frequency instead of a fitting band defined by a water-level. In the manual revision, we adjust the frequency fitting band. Depending on whether the inflection point in the spectral ratio near f c2 is clearly visible or not, we fit the ratio using either the source (Equation 11) or the spectral ratio (Equation 10). We combine all f c1 , f c2 estimates for individual target, eGf events from different pair combinations using the delete-one jackknife mean as described before. We refer to the f c1 and f c2 estimates from averaged target and eGf events using the semi-automated workflow as RRTA, RREA, respectively ( Figure 2).
For each of the three spectral ratio f c estimates described above, we calculate corresponding Δσ values using the M 0 estimates from single spectra, and spectral ratio corner frequency estimates corresponding to automated (RSTA), and semi-automated (RRTA, RREA) estimates.

Results
We recover seismic moment, corner frequency, and stress drop estimates for 16,613 events using single spectrum fitting ( Figure S8 in Supporting Information S1). We note that our calculated M w values for M ≤ 3 are slightly higher than the M w vs. M L relationship from previous studies in the Apennines (Munafò et al., 2016;Malagnini & Munafò, 2018, Figure S8 in Supporting Information S1).
The corner frequency estimates from the single-spectra fitting and CEM show similar results. However, single-spectrum estimates deviate to lower values for M ≤ ∼ 2.5 ( Figure 3 and Figure S9 in Supporting Information S1). Furthermore, CEM estimates are largely missing for M ≥ 5 earthquakes (as would be expected due to the smaller number of larger earthquakes). The corner frequency estimates from CEM and the spectral ratio approach show similar values for overlapping event estimates (Figure 3 and Figure S10 in Supporting Information S1), where spectral ratio estimates show a generally higher trend below M ∼3. Overall, the Δσ estimates show a log-normal distribution with a median value of 3.32 MPa for all events, and of 5.21 MPa for events with M w ≥ M cspec (Figure 4c). M cspec is the magnitude of completeness of all spectral estimates, which we calculate following Aki (1965) and Wiemer and Wyss (2000) (Figure 4b). We use M cspec to estimate the magnitude threshold under which we expect spectral estimates are likely to be missing, which is in our case M w ≤ 2.2. The stress drop values also show generally decreasing values with magnitude for M w ≤ ∼3. Due to the scaling of stress drop with magnitude and the fact that most events are smaller than M cspec , the stress drop results are not directly comparable in an absolute sense over the full magnitude range considered here.
To overcome biases related to the apparent (natural or systematic) scaling of stress drop with magnitude, and because most events are smaller than M cspec ∼2.2, we consider the relative differences of stress drop values in the following discussion and interpretation. By using relative differences, we transfer the stress drop estimates into mutually comparable values over a wide magnitude range as we eliminate the influence of stress drop scaling with magnitude. Doing so helps avoid bias in any interpretation of the spatial and temporal evolution of stress drop from the more abundant small magnitude events. To generate a catalog with relative stress drop estimates, we first combine all stress drop estimates and use one estimate for each event, where the order of priority is: RREA, RRTA, RSTA, CEMA, SSSA (only M w ≥ 3 for SSSA). We then calculate the delete-one jackknife stress KEMNA ET AL.

10.1029/2021JB022566
13 of 29 drop mean in the logarithmic domain for several magnitude bins (M w ≤ 1.5, 0.5 steps for 1.5 < ≤ 4, M w > 4; BJM; red points in Figure 4a), and normalize each stress drop estimate by the estimate from its respective magnitude bin to calculate what we refer to in the discussion as the magnitude-normalized stress drop Δ using the logarithmic change Δ = 10(Δ ∕Δ ) . Computing the magnitude normalized values allows us to discuss the spatio-temporal distribution of relative stress drop values without any bias related to stress drop -magnitude scaling. Synthetic tests show that the magnitude-normalized stress drop values of the CEM approach are able to reliably recover relative differences in stress drop values for all magnitudes. However, the test indicates that the correlation is marginally less robust for events with M w < 2 ( Figure S11, see Supporting Information S1 for further details). We therefore only consider events with M w ≥ 2 in the discussion of spatio-temporal patterns in magnitude-normalized stress drop. The color bars in Figure 4a show the theoretical range of Δ for a specific magnitude bin. We then look at spatio-temporal patterns in magnitude-normalized stress drop Δ in the AVN seismic sequence as a proxy for seismic response and the stress state in a specific location at a specific time during the seismic sequence. Figure 5a shows the temporal variation of Δ along the strike of the Laga Fault System (LFS) and Vettore-Bove Fault System (VBFS). Figures 5b and 5c show the moving median Δ over time for the whole dataset (b), for events in the VBFS (c), and events in the LFS (d). Events in the VBFS exhibit generally higher Δ compared to events in the LFS (Figures 5 and 6). The Δ values show a median increase in the time period directly following the three largest earthquakes, with its highest peak after the Amatrice earthquake. The (c) and (d) panels exhibit the different responses in both fault systems.
In the VBFS, the median Δ remains high for several weeks during pronounced aftershock activity at the northwestern edge of the VBFS. The median Δ tends to increase again after several months of below-average Δ , starting from the middle of 2017 until the end of the study period (Figure 5c). In the LFS, the median Δ increases after a few weeks of aftershock activity following the Amatrice earthquake and persists at the same level until several weeks after the Norcia earthquake. The median Δ then decreases drastically following the Norcia mainshock. A short and sharp increase of the median Δ occurs during aftershock activity at the south-eastern edge of the LFS in concert with several M w ≥ 4 earthquakes, followed by a steady increase until the end of the study period (Figure 5d).
We will describe in the following the spatio-temporal patterns of Δ from the early aftershocks (within 4 days) of each mainshock, between the Amatrice early aftershocks and the Visso mainshock, and after the Norcia early aftershocks. Figure 6 provides a map view of all earthquakes, and Figure 7 shows a cross-sectional view of the double-difference relocated earthquakes projected on the fault model of Lavecchia et al. (2016). Each subplot shows the spatial distribution of Δ in the five time windows indicated on the top right. Figure 8 shows cross-sections in greater detail indicated by the cross-sections marked in the respective time windows of Figure 6.

Amatrice Earthquake
The M w 6.0 Amatrice earthquake (24 August 2016) activated two distinct fault segments in a bilateral rupture: the LFS to the south and the VBFS to the north (Calderoni et al., 2017;Lavecchia et al., 2016). The median Δ increases directly after the Amatrice earthquake, with higher values in the VBFS compared to the LFS. The Δ values then decrease rapidly in the following days of aftershock activity ( Figure 5). There is a dearth of aftershocks in the LFS surrounding the Amatrice hypocenter corresponding to the area of maximum slip from finite-fault inversions (Lavecchia et al., 2016;Walters et al., 2018, Figure 7a). The aftershocks of the Amatrice earthquake show a binary pattern of elevated and decreased Δ values to the north and south of the maximum slip area, respectively ( Figure 7a). The areas of maximum slip are broadly outlined in Figure 7 for each mainshock and are based on a composite of various slip models (Lavecchia et al., 2016;Pizzi et al., 2017;Scognamiglio et al., 2018;Walters et al., 2018). To the north, the Amatrice earthquake activated several synthetic, antithetic, and subhorizontal structures in the southern part of the VBFS (Michele et al., 2020, Figures 8a and 8b). In addition, the lateral ramp of the MST was reactivated in at least one thrust-faulting event (Improta et al., 2019;Michele et al., 2020, Figure 8b and Figure S12a in Supporting Information S1). The majority of northerly, early Amatrice aftershocks have an elevated Δ , where the highest values are located in zones near fault intersections (Figures 8a and 8b).
In the two months following the Amatrice earthquake, the median Δ decreases relative to the early part of the aftershock sequence. In the VBFS, the median Δ stays elevated relative to the median of all cataloged earthquakes and the median of earthquakes preceding the Amatrice earthquake (Figures 5b and 5c). In the LFS, the median Δ decreases rapidly after the early aftershocks to values well below the mean. However, the median value recovers to median-levels after a few weeks of aftershock activity. Cross-sectional views that are roughly perpendicular to strike show that earthquakes are mostly located in areas of the early aftershocks of the Amatrice earthquake (Figures 7b, 8c and 8d) but are more concentrated in the northern part of the Amatrice rupture zone. The aftershocks in the VBFS show higher Δ than events in the LFS (Figure 6a). One notable feature is that events with elevated Δ in the VBFS are located at the intersection zone of fault zones as described above (Figure 8b). Furthermore, events with low Δ occur in areas with few early aftershocks and at the southern edge of the early aftershock cloud of the Amatrice earthquake.

Visso Earthquake
The M w 5.9 Visso earthquake (26 October 2016) ruptured the northern and central segments of the VBFS, occurring as a doublet together with its M w 5.4 foreshock (Improta et al., 2019;Michele et al., 2020). The median Δ of the early aftershocks increases slightly but stays below the increase after the Amatrice earthquake ( Figure 5). Earthquakes are sparse in areas with maximum slip inferred from slip-inversion, and show a similar pattern to the aftershocks of the Amatrice earthquake (Pizzi et al., 2017;Walters et al., 2018, Figures 7c and 8e). They also reflect similar patterns to areas with previous aftershock activity, such as in the vicinity of the MST (Figure 7c). A number of earthquakes with both elevated and decreased Δ occur around the maximum slip area with a lack of an obvious spatial pattern (Figure 7). Many events with a slightly higher Δ occur in deeper parts of the fault zone at the intersection of the subhorizontal shear zone (∼8 km depth; Figure 8e).

Norcia Earthquake
Four days after the Visso earthquake, the M w 6.5 Norcia earthquake (30 October 2016) ruptured multiple faults segments of the VBFS, the northern part of the LFS, and several secondary and antithetic faults. The earthquake ruptured the surface a length of over 35 km along the VBFS with a maximum surface displacement of up to two meters ( Figure 6, red lines) (e.g., Improta et al., 2019;Michele et al., 2020). The median Δ increases drastically  Figure 1). Each grid rectangle (1.5 × 1.5 km) is colored by the mean Δ of the earthquakes contained inside. Only double-difference relocated earthquake hypocenters from Michele et al. (2020) are plotted. Light gray, black outlined, and slate gray circles denote earthquakes before, in (without a spectral estimate), and after the specific time window, respectively. Mainshock locations (colors and symbol are the same as in Figure 1) are shown in all subfigures independent of the timing of occurrence relative to a specific time window. Areas of maximum slip are roughly outlined for each mainshock with the dashed circular shape in the corresponding color from Figure 1 (after Lavecchia et al., 2016;Pizzi et al., 2017;Scognamiglio et al., 2018;Walters et al., 2018). Earthquakes potentially related to fluid-diffusion are broadly marked with the cyan dashed ellipse in (b).

of 29
for the early aftershocks in the VBFS, but stays below the Amatrice earthquake maximum increase ( Figure 5). The majority of high Δ events are located west and south of the mainshock at depths exceeding five kilometers, and similar to the Amatrice aftershocks, are mostly absent in the zone of maximum slip (Scognamiglio et al., 2018, Figure 7d). In the LFS, high Δ values are located in areas with high aftershock production following the Amatrice aftershock sequence. However, aftershocks are fewer in number compared to the Amatrice aftershock sequence (Figures 7d and 8f). Furthermore, the depth distribution of early Norcia aftershocks resembles the  Figure 7. The sketches (black-dashed lines) show interpretations of fault geometry at depth inferred from previous studies Improta et al., 2019;Michele et al., 2020, the sketch of [b] applies to [d,f]). Focal mechanism solutions are taken from Scognamiglio et al. (2006). Focal mechanisms corresponding to the Visso (e) and Norcia (f) mainshocks follow the same color scheme as in Figure 1, and the same symbol denotes the hypocentral location. The dashed hexagons with a white-to-red color gradient correspond to the future hypocenter of the Norcia earthquake (b, d).

10.1029/2021JB022566
19 of 29 depth distribution of the early Visso and Amatrice aftershocks. The earthquakes south-east of the Norcia rupture zone show a very complex interaction between synthetic, antithetic, and subhorizontal structures with elevated Δ values near the junction of fault surfaces, and contain several M w ≥ 4 events (Figure 8f).
In contrast to the latter part of the Amatrice aftershock sequence, the late aftershock sequence of the Norcia earthquake shows significantly higher Δ values persisting for several weeks. High Δ events mostly locate in the subhorizontal shear zone at the bottom of the seismogenic zone and around the intersection with the MST fault zone (Figures 6e and 7e). Low Δ events mostly fill areas with previous low aftershock activity (Figure 7e). In addition, a large cluster of high Δ events rupture the northern tip of the VBFS. The northern tip contains intersections of synthetic and antithetic structures and several strike-slip faults linking the dominant normal fault structures (Figures 7e, 8g and Figure S12e in Supporting Information S1). Following the cluster of events at the northern tip of the VBFS, the median Δ decreases drastically to below-average values for several months. A second large cluster of events ruptures the southern tip of the LFS with several M w ≥ 4 events, leading to a sharp increase Δ values (Figures 5d and 8h). The median Δ stays below-average at values with small fluctuations close to the median before the Amatrice earthquake, where the highest Δ events consistently occur in the southern part of the VBFS, that is, the vicinity of the Mt. Vettore fault (Figure 7e). In the second half of 2017, the median Δ starts to increase, however, event density also thins out with increasing time ( Figure 5).

Discussion
The AVN seismic sequence shows significant spatio-temporal variations in stress drop values. The median Δ increases after the three largest earthquakes (Amatrice, Visso, Norcia earthquake). Stress drop patterns also change in time within different fault systems (LFS, VBFS). The values tend to increase at fault intersections following large earthquakes and near the terminus of maximum slip regions inferred from finite fault inversions. Below we discuss the general stress drop scaling with magnitude, the spatial and temporal variation of the sequence as a whole, and the relation to geological structures. The last subsection discusses the general implications, including the relation to other studies.

Stress Drop Scaling
The breakdown of self-similar scaling suggested by Figures 3 and 4 becomes apparent for earthquakes with M w ≤ 2-3. Wang et al. (2019) interpreted this breakdown in scaling at small magnitudes to result from an insufficient increase of slip with rupture dimension due to a strongly heterogeneous stress field (Ben-Zion, 2008). Complex fault geometry may also play a factor in limiting rupture, either in dimension and/or slip (Chiaraluce et al., 2005;Improta et al., 2019;Michele et al., 2020). In the presence of complexity, an expected break in scaling might be interpreted as resulting from a greater capacity of larger earthquakes to rupture through and break potential rupture barriers, whereas smaller earthquakes are less likely to do so. An alternative explanation for the break in the observed scaling may be the unexceptional use of surface stations in this and in previous studies in the region. Ground motion measurements from surface stations would be depleted in high frequency energy above ∼25-30 Hz, depending on the station installation site. The depletion of high frequency energy would bias the measured and estimated f c to lower values, which would disproportionally affect smaller magnitude earthquakes with higher theoretically expected f c values. Based on an assumed constant stress drop scaling and the observed SNR for a given station, one can infer where the observational bandwidth limits are for a given data set, and the theoretical f c and magnitude values that would be affected.
For example, Figure 4a shows the stress drop estimates ranging over six orders of magnitude and a theoretical measurement limit of = 25 Hz. Frequencies above 25 Hz are in a range where the depletion of high frequencies is expected for surface stations; 25 Hz is also 1/2 of the Nyquist frequency of the majority of the instruments (100 Hz sampling rate). Therefore, 25 Hz likely represents the upper limit of the observable bandwidth to estimate f c values, even if SNR values exceed 3 for larger events at frequencies up to the Nyquist value. Abercrombie (2014) suggests that the frequency bandwidth limit for estimating is ∼2-3 ⋅ f c using spectral ratio methods. Therefore, 25 Hz could be the initiation point at which the instrument response bandwidth limit is approached, even if SNR is high up to 40/50 Hz. The decreasing trend in stress drop values measured here tends to follow the 25 Hz isoline (Figure 4a), suggesting it indeed represents the upper frequency limit of fitting. Further, KEMNA ET AL.

10.1029/2021JB022566
20 of 29 the synthetic test shows an apparent boundary at ∼25 Hz (Figure S11a in Supporting Information S1). The fact that the upper boundary of stress drop estimates so closely follows the 25 Hz isoline points to a systematic bias related to the use of surface stations and instrument response bandwidth limitations that shape the observation of decreasing Δσ with M w . The systematic bias would suggest that high stress drop estimates (≥50 MPa) would become unresolvable for M w < ∼3 due to the bandwidth limit. The spectral ratio eGf f c estimates that are beyond the instrument response bandwidth suggest that high Δσ are possible, but would not be resolved by the current seismic network of surface stations (see symbols for RREA with arrows suggesting a lower bound estimate, Figure 4 and Figure S14 in Supporting Information S1). Therefore, we interpret the observation of decreasing Δσ with M w to more likely result from observational limits rather than geological factors. We want to emphasize that we are not able to refute any observed stress drop scaling for events below M w < ∼2,3. However, as other studies are likely to be affected by the same instrument response and SNR limitations, the available data do not allow establishing or refuting definitive scaling of stress drop with M w < 3. The results presented here suggest that stress drop values may show no discernable trend with M w < 3. It does, however, show (relative) changes in space and time, as suggested by the magnitude-normalized values.

Spatio-Temporal Stress Drop Evolution
Larger earthquakes with longer fault lengths may no longer reflect local physical properties of specific fault patches, but rather, averaged properties over fault planes with multiple rupture patches. We may therefore assume that stress drop values for earthquakes with M ≤ 4 (corresponding to rupture patches ≤1 km) may reflect the stress release and the physical properties on a specific fault patch (also assuming no significant overshoot or undershoot) (Yamada et al., 2017). Because static stress drop is an averaged parameter, stress drop estimates of small earthquakes have the potential to illuminate finer details related to spatial (and temporal) changes in stress release or fault rock properties compared to estimates for larger earthquakes. Furthermore, they offer insights into the heterogeneity in the ambient stress field and its evolution in time. For example, larger earthquakes (M w ≥ 5) during a given seismic sequence release a considerable amount of stress rapidly, resulting in a significant reduction of the overall stress level. Nevertheless, the stress release is not uniform and may be distributed over several fault patches. This non-uniform distribution then governs the evolving stress field surrounding the main slip areas of successive events in a sequence, as well as the corresponding aftershock nucleation (Das & Henry, 2003;Hainzl & Marsan, 2008). In that context, the regional stress field of the AVN seismic sequence would be expected to evolve in a transient, highly heterogeneous fashion.
In order to lend support to the interpretations of the spatio-temporal variation of stress drop estimates, we perform a first-order investigation of the implied stress field evolution. We use focal mechanism solutions passing the quality control (as described in Section 2.1, Figure S12 in Supporting Information S1) to estimate the temporal evolution of the principal stress directions σ 1 , σ 2 , σ 3 , and the shape ratio R = (σ 1 − σ 2 )/(σ 1 − σ 3 ) through an iterative joint focal mechanism inversion (Vavryčuk, 2014). To quantify temporal and spatial changes in the stress tensor, we first associate the focal mechanisms to the VBFS and LFS (separated by the MST) and group the two subsets into moving windows of 20 events overlapping by 15 events (i.e., we move the the window of events in increments of five). We refer the reader to the Supporting Information S1 for further details of the standard stress inversion procedure. Figure 9 shows the temporal evolution of the shape ratio (a) and the plunge of the stress axis (b) in the two fault systems for each window of events, as well as polar stress plots for six specific inversion windows indicated in panels (a) and (b). Windows 2, 4, and 5 all have a nearly vertical σ 1 in an undisturbed stress field in the LFS and VBFS (Figure 9), which is then repeatedly perturbed by the occurrence of large earthquakes (M w ≥ 4). The shape ratio and the stress axis orientations are particularly variable in the VBFS, and changes correlate with the late aftershocks of the Amatrice earthquake and both the early and late Norcia aftershocks (Figures 9-1,3). In the LFS, the stress field changes are more gradual and less extreme in magnitude during the course of the seismic sequence, which is also due to fewer earthquakes with available focal mechanism solutions. Specifically, the stresses appear perturbed during aftershock activity at the southern tip of the LFS associated with the occurrence of several M w ≥ 4 events (blue vertical lines). We note that similar patterns emerge when testing different sliding windows (25 and 30 with 2025 overlapping events). The larger windows produce patterns that are less pronounced, which might be an indicator of short-lived changes related to local stress rotations near the large earthquakes. If so, stress field variations may not be representative of the regional stress field over longer periods of time. Nevertheless, the changes are consistent with the spatio-temporal distribution of magnitude-normalized stress drop, as discussed below. While the stress inversion results shown in Figure  signal related to the large mainshocks, the correlation with the spatial distribution of stress drop results suggest that at least some of the signal is related to the geological conditions. Thus, further detailed study devoted to the evolution stress orientation during the sequence is likely warranted.
The observed patterns of stress orientation and shape ratios correlate with the observed patterns in magnitude-normalized stress drop Δ both in space and time (Figures 5 and 6). In space, the VBFS exhibits higher magnitude-normalized Δ values compared to the LFS and correlates with higher variability of the stress field in the VBFS. Large earthquakes also likely cause temporal perturbations of the stress field and may even rotate the principal stress directions. The temporal variations observed here are not unique. Hardebeck and Okada (2018) presented a series of examples of stress rotations following large earthquakes for a range of tectonic settings. They infer that the observed stress rotations could suggest a low background differential stress level that may also highlight spatial stress heterogeneity. Similarly, the spatial stress heterogeneity inferred here is consistent with the observed heterogeneity in Δ and the stress field inversion. The evolving stress field could then lead to the nucleation of aftershocks on faults/fault patches with a high fault strength as suggested by an elevated Δ , which would typically not otherwise be activated under the ambient loading conditions. Therefore, Δ changes are likely highly influenced by changes in the stress field and could be used as a proxy to highlight local changes in stress loading.
Furthermore, the Δ values are higher in the VBFS relative to the LFS (Figures 5c, 5d and 6). In addition to previously described differences, the VBFS is characterized by increased fault complexity, which we infer from mapped faults (east-and west-dipping normal faults, Figure 1), the stress inversion, and observed fault junctions at depth (Figure 8). The higher magnitude-normalized stress drop values observed here are consistent with previous studies in southern California, which also associated higher absolute stress drop values to increased fault complexity (Goebel et al., 2015) and stress field heterogeneity (Goebel et al., 2017). Fault complexity may also lead to heterogeneous fault properties, including strong asperities, which could in some cases arrest growing ruptures before they reach large sizes (Sammonds & Ohnaka, 1998). Both the stress field heterogeneity and fault complexity are likely to be mutually supportive in generating a heterogenous distribution and evolution of differences in magnitude-normalized (and absolute) stress drop values. For example, an increased fault complexity with different potential asperities would facilitate the transient change of the stress field. The high Δ values in the vicinity of the MST that persist over the entire study period may also highlight high local fault strength in general (Figures 5 and 7). If the MST represents a locally strong fault system, it may also have limited the Amatrice earthquake rupture as well as provided the region of stress concentration necessary to nucleate the Norcia earthquake (Improta et al., 2019).
The median magnitude-normalized Δ values also change significantly over time, which may highlight temporal changes in fault frictional strength. For instance, high Δ earthquakes occur adjacent to maximum slip patches inferred from slip-inversion models, which could reflect the release of stress concentrations at areas with high frictional strength that did not rupture during the mainshock (Chen & Shearer, 2011;Oth & Kaiser, 2014;Yamada et al., 2010, Figure 7). The median Δ values gradually decrease to values below the median following the early aftershocks. The latter, lower Δ could be associated with a temporary reduction in fault strength within the mainshock rupture zone (e.g., Abercrombie, 2014;Kim et al., 2016;Trugman et al., 2017) and/or highlight fault segments that have not fully healed from the large earthquake ground motion (e.g., Moyer et al., 2018;Shaw et al., 2015). Consistent interpretations by Yamada et al. (2021) associated high stress drop values with a gradual weakening of the shear strength due to stress concentration during the coseismic slip of the mainshocks. Following the temporal decrease, the median Δ values then increase in the subsequent months, consistent with fault zone healing and an increase in fault strength (e.g., Goebel et al., 2015;Moyer et al., 2018;Vidale et al., 1994). In our case, the median Δ value increases above the median of earthquakes prior to the Amatrice earthquake Hexagons denote the occurence of the three largest earthquakes with the same color scheme used in Figure 1. The vertical teal-colored lines denote the occurence of M w ≥ 4 earthquakes. Black filled diamonds denote the specific inversions plotted in (c) and are numbered accordingly (1-6). (b) Plunge of the principal stress axis (principal stress colors are defined in the legend) vs. time for the VBFS (top) and LFS (bottom). Colors and sybmols of (a) apply. Errorbars are noted on the symbols, and are commonly within the boundaries of the symbol. (c) Polar plots of specific stress inversions denoted by the numbers (1-6) in (a, b) (see Figure S13 in Supporting Information S1 for the inverted FM solutions). Colors for the principal stresses from (b) apply. The diamonds in (c) denote the stress inversion values using all focal mechanisms and the filled circles resulting from 80 random noise realizations.

10.1029/2021JB022566
23 of 29 ( Figure 5). However, as there are only a few estimates between 1 January 2016 and the Amatrice earthquake, we cannot robustly determine the latent regional Δ median value.

Geological, Tectonic, and Hydrological Influence on Earthquake Stress Drop
The earthquakes in the VBFS and LFS, which are cross-cut by the MST, show a clear difference in magnitude-normalized Δ values, where the VBFS hosts consistently higher Δ events (Figures 5a, 5c,, and 6e). The primary differences may be associated with the different geological settings. The MST separates both fault systems and two geological domains, the Laga Foredeep domain to the south and the Umbria-Marche domain to the north. The Umbria-Marche domain consists of stacked thrust sheets of carbonates and evaporites (Umbria-Marche succession), and turbidites of siliciclastic deposits characterize the Laga Foredeep domain, covering the Umbria-Marche succession (Barchi et al., 2021;Porreca et al., 2018;Sebastiani et al., 2019). A second factor dictating strength may be the complexity of the fault system, which is significantly higher in the VBFS compared to the LFS, with multiple mapped synthetic and antithetic normal faults overprinting the pre-existing thrust faults on the surface   Figure 1) and at depth (Barchi et al., 2021;Porreca et al., 2018). The higher fault complexity is consistent with the observed elevated Δ values as well as the stress inversion results (Section 4.2, Figure S12 in Supporting Information S1). The influence of the 1D velocity model used in the stress drop calculation is negligible, as the changes of stress drop values correspond to non-realistic changes in the velocity (e.g., a ∼35%-25% velocity increase-decrease would be required to explain an observed Δ of −0.4 -+0.4 in the magnitude bin M w 3-4, respectively). Furthermore, the geological cross-sections by Barchi et al. (2021) suggest that the earthquakes in the LFS nucleate in the same Umbria-Marche succession, evaporites and carbonates, as in the VBFS or in the transition zone containing the turbidites. Therefore, the lower Δ values in the LFS suggest that the data are able to resolve the effects of rock strength, even within the resolution of the 1D velocity model.
Low magnitude-normalized Δ events could also be related to fluid diffusion processes (e.g., Chen & Shearer, 2011). As observed for the 2009 L'Aquila seismic sequence, fluids seem to be a significant driver of aftershock evolution and migration (Lucente et al., 2010). The presence of fluids is also observed to produce fault fracture meshes at shallow depths (Ross et al., 2017;Sibson, 1996), which could limit the extension of earthquake rupture (Wei et al., 2011) or channel pressure-driven fluid flow (Walters et al., 2018). The latter observed migrating aftershocks in the VBFS preceding the Visso earthquake at fluid diffusion speeds (Figure 7b, cyan-dashed ellipse). The area where the migrating cluster occurs is also associated with several zones with low Δ values, which could suggest a reduction of frictional strength, possibly due to the presence of fluids (e.g., Goertz-Allmann et al., 2011;Yamada et al., 2015). While approaching the Visso earthquake hypocenter, the earthquakes begin to show elevated Δ values. The migrating events are interpreted to play a role in the triggering of the Visso earthquake, bypassing the hypocenter of the Norcia earthquake at shallow depths (Walters et al., 2018). The presence of fluids and possible pressure-driven flow is highly dependent on the geological and hydrological setting. The setting is different to the north and south of the MST, as described above. In the VBFS, alternating impermeable (evaporites) and permeable (calcareous) units are more likely to channel fluid flow or overpressured fluids in the upper layers (Improta et al., 2014). Miller et al. (2004) also proposed that aftershocks of the 1997 Colfiorito Central Italy seismic sequence may be driven by the coseismic release of trapped, high-pressure fluids propagating through zones of damaged rock created by the mainshock. The lithological change could explain the observed seismic migration in the VBFS and the absence of such patterns in the LFS, dominated by the Flysch della Laga Formation's turbidites at the surface. Furthermore, Michele et al. (2020) hypothesized that the missing shallow aftershocks in the LFS (≤5 km) are related to lithology, as the turbidites could be velocity strengthening materials and respond only with aseismic deformation (e.g., stable creeping) when locally stressed. The lithology change could also explain the low Δ events in the transition zone between velocity strengthening and velocity weakening materials. The lithological changes also likely invalidate our assumption of a constant rupture velocity, as the shear velocity would be lower for the weaker turbidites. However, it would not invalidate the common assumption of the rupture velocity being some constant proportion of the shear velocity (0.8 β, Section 2.2). If the turbidites have a lower shear, and thus lower actual rupture velocity compared to the more competent rock units to the north, it would lead to lower stress drop estimates in our calculation (Allmann & Shearer, 2009;Kaneko & Shearer, 2015). Therefore, the lower Δ values in the LFS could also be interpreted as earthquakes with lower rupture velocities and a similar magnitude-normalized stress drop compared to the VBFS.

Comparison to Previous Studies
Previous studies in the Apennines found a breakdown of the self-similar scaling of the seismic moment with corner frequency for small (M w ≤ ∼3 or 4) earthquakes (Bindi et al., 2020;Morasca et al., 2019;Wang et al., 2019). We observe a similar pattern, where the stress drop increases with increasing magnitude for earthquakes with M w ≤ ∼ 3 (Figure 4). For earthquakes with M w ≥ 3 we observe self-similar scaling with a constant stress drop ∼6.5 MPa, in contrast with previous findings that observe increasing stress drop with magnitude over the entire magnitude range studied. Reasons for the discrepancies may rest in the approaches used to estimate f c . For example, Wang et al. (2019), comparable study period and Bindi et al. (2020), longer time frame spanning a wider area used the nonparametric generalized inversion technique (GIT). Both studies found lower average stress drop values for events M ≤ 4, and increasing stress drop values with magnitude. In contrast, we find no apparent scaling of stress drop with magnitude for M ≥ ∼ 3, suggesting discrepancies between methodologies. For example, averaging of results associated with groups of events during the GIT inversion could suppress differences of source properties at nearby locations, where the spectral ratio method is able to resolve such differences (e.g., Ross & Ben-Zion, 2016). A similar suppression would also influence the CEM estimates and could explain the difference between CEM and spectral ratio estimates within this study as well. We try to account for lateral changing source properties and attenuation effects by limiting the cluster spatial extension (≤4 km) and resolving each event in multiple clusters. CEM estimates mostly overlap with the estimates from spectral ratios ( Figure S10), suggesting the suppression, if present, is minimal here. Another set of studies used coda wave spectra (Malagnini & Munafò, 2018;Morasca et al., 2019) and also found scaling of stress drop with magnitude. The discrepancy with our results is in agreement with previously described differences between direct and coda wave based spectral analysis and warrants further study to be solved in a direct comparison of the two approaches (e.g., Abercrombie, 2013). Furthermore, there is a general lack of CEM and spectral ratio estimated f c values for M w ≥ 5 events here due to the data assumption requirements for the two methods, which further complicates the comparison with previous studies. The smaller events generally have SNR at the stations at which events with M w ≥ 5 have estimated spectra. In addition, the recordings of the larger events exceed the dynamic range at closer seismic stations, where the majority of spectra for smaller events are observed. The missing estimates could be explained by missing station spectra for events with M w ≥ 5 compared to smaller magnitude earthquakes, which are necessary to form clusters or serve as eGf events.

Conclusions
The spatial-temporal evolution of absolute and magnitude normalized stress drop values for over 16,000 earthquakes from the AVN seismic sequence shows a correlation to the spatio-temporal distribution of large earthquakes, and to geological conditions. For example, high magnitude-normalized stress drop values correlate with fault junctions and stronger geological host rock. As large earthquakes relieve large amounts of accumulated shear stress, the fault slip over a large area introduces new stress concentration at the boundaries of the the rupture path. We observe temporal variations in magnitude-normalized stress drop corresponding with the stress field variance expected from earthquake magnitude, resulting in an elevated median stress drop in time periods of a few days (Amatrice earthquake) -∼1 month (Visso and Norcia earthquakes) after large earthquakes. High-stress drop events are observed to cluster in the vicinity of high fault slip areas and zones with increasing fault complexity. Zones of high fault complexity include fault intersection zones of the VBFS and the MST thrust front. We observe immediate increases in magnitude-normalized stress drop estimates directly following the three largest mainshocks, followed by decreases below the median and then recovery during continuing aftershock activity. We interpret the recovery observed late in the study period as being associated with fault healing processes.
Similar to previous studies, we also observe a breakdown of self-similar stress drop scaling with magnitude for M ∼ 2.5. In contrast to previous studies, we observe no discernable trend with earthquake size for M ≥ 3, and interpret the scaling for M ≤ 3 to be governed by observational limitations. We note that if there is any natural scaling of stress drop with magnitude in the study area, the available seismic network of surface stations would not be able to resolve it. Finally, our methodology reiterates the commonly observed discrepancy of stress drop estimations between different studies, methods, and portions of the seismic signal used (coda vs. direct waves). We want to emphasize the urgency of a comprehensive comparison between different approaches to calculate KEMNA ET AL.

10.1029/2021JB022566
25 of 29 stress drop values from seismic wave spectra on common dataset(s) to identify the cause(s) of measurement discrepancies.