Influence of Fault System Geometry and Slip Rates on the Relative Role of Coseismic and Interseismic Stresses on Earthquake Triggering and Recurrence Variability

We model Coulomb stress transfer (CST) due to 30 strong earthquakes occurring on normal faults since 1509 CE in Calabria, Italy, including the influence of interseismic loading, and compare the results to existing studies of stress interaction from the Central and Southern Apennines, Italy. The three normal fault systems have different geometries and long‐term slip‐rates. We investigate the extent to which stress transfer can influence the occurrence of future earthquakes and what factors may govern the variability in earthquake recurrence in different fault systems. The Calabrian, Central Apennines, and Southern Apennines fault systems have 91%, 73%, and 70% of faults with mean positive cumulative CST in the time considered; this is due to fewer faults across strike, more across strike stress reductions, and greater along‐strike spacing in the three regions respectively. In regions with close along strike spacing or few faults across strike, such as Calabria and Southern Apennines, the stress loading history is mostly dominated by interseismic loading and most faults are positively stressed before an earthquake occur on them (96% of all faults that ruptured in Calabria; 94% of faults in Southern Apennines), and some of the strongest earthquakes occur on faults with the highest mean cumulative stress of all faults prior to the earthquake. In the Central Apennines, where across strike interactions are the predominant process, 79% of earthquakes occur on faults positively stressed. The results highlight that fault system geometry plays a central role in characterizing the stress evolution associated with earthquake recurrence.

At the scale of the fault system and over the long-term (60 ka to multiple Ma), earthquake recurrence intervals are mainly controlled by the regional strain-rate (Mouslopoulou et al., 2009;Nicol et al., 2005Nicol et al., , 2016)).For individual faults, the variability in recurrence intervals on shorter timescales (<20 ka) is primarily influenced by fault length and slip-rates (Nicol et al., 2010(Nicol et al., , 2016)), but may be modulated by other factors, such as stress interactions between neighboring faults (Cowie et al., 2012;Mouslopoulou et al., 2009;Nicol et al., 2005).Calculations of cumulative Coulomb stress changes that include coseismic and interseismic stress can help identify differences in the way stress builds up across a fault system, depending on the geometry of its constituent faults and shear zones, with a fluctuating stress loading history for faults that are interacting coseismically across strike, and an almost linear increase in stress, dominated by interseismic loading, for faults that are isolated, that is, do not have other faults across strike (Sgambato, Faure Walker, Mildon, & Roberts, 2020).This suggests that earthquake recurrence might be different for regions where faulting is localized across a narrow fault system, compared to those where the deformation is distributed across multiple faults, and this has been previously demonstrated using numerical simulations (Cowie et al., 2012;Visini & Pace, 2014).Understanding the impact of this behavior on the recurrence time is important for seismic hazard assessment because, although earthquake recurrence is directly related to slip-rate, and thus to the regional strain-rate (Mouslopoulou et al., 2009;Nicol et al., 2005Nicol et al., , 2010Nicol et al., , 2016)), fault interaction will affect the variability in earthquake recurrence intervals (Cowie et al., 2012;Nicol et al., 2005Nicol et al., , 2010Nicol et al., , 2016)), defined with the coefficient of variation (CV), and the slip-rate variability (SRV) (Cowie et al., 2012).CV is calculated as the ratio between the standard deviation of the recurrence time and the long-term average earthquake recurrence.SRV is calculated as the ratio of standard deviation of short-term slip rates over a sliding time window of fixed length and the long-term average slip-rate (Cowie et al., 2012).Within a fault system, both the mean recurrence interval and CV vary with fault length and displacement rates, with a higher variability in recurrence interval, thus higher CV, for smaller faults with lower slip-rate (Mouslopoulou et al., 2009;Nicol et al., 2010Nicol et al., , 2016)).Thus, both CV and SRV vary spatially in an extensional fault system, depending on the length and deformation rates, and on the position and orientation of the faults, due to fault interaction (Cowie et al., 2012;Robinson et al., 2009).Specifically, in areas of the fault system where the average slip rates are higher and deformation is localized on one or few major faults, the dominant effect of nearby ruptures is the along-strike stress increase, leading to short average recurrence intervals and lower CV and SRV values.Where the deformation is distributed across multiple faults, the combination of increases and decreases in stress through time leads to higher CV and SRV values, and the distribution of earthquake recurrence intervals becomes exponential (Cowie et al., 2012).This suggests that fault system geometry and stress transfer control the earthquake recurrence intervals (Cowie et al., 2012;Visini & Pace, 2014), yet many existing models for probabilistic seismic hazard assessment (PSHA) do not include variations induced by stress changes, and constant values of CV are imposed along entire fault systems (e.g., Pace et al., 2006;Toda et al., 1998;Valentini et al., 2017).Where the effect of fault interaction is included (e.g., Akinci et al., 2017;Pace et al., 2014;Xiong et al., 2017), often only the effect of coseismic CST is considered to recalculate probability of occurrence of the next earthquake(s).
In this paper, we focus on investigating the effect of the fault system geometric complexity, indicative of the number of faults in a system and their relative location (i.e., whether they are arranged predominantly across-strike or along strike) on the distribution of stress and its influence on earthquake recurrence, and provide a natural example of the variations in fault system geometry and stress loading of the type proposed by Cowie et al. (2012).We model the sequence of earthquakes with Mw > 5.5 that occurred on the active normal faults in Calabria and Northern Sicily since 1509 CE (Rovida et al., 2020(Rovida et al., , 2022)), including interseismic loading from their underlying shear zones.The earthquake sequence in Calabria is compelling because many large earthquakes were concentrated in clusters in the years between 1638 and 1908 CE, and separated by a few months or even days (Galli & Bosi, 2002).Some of these events were the strongest recorded in Italy, such as the devastating Mw 7.1 Messina earthquake in 1908 with a death toll of ∼70,000 people (Baratta, 1910), and the earthquake sequence in 1783 CE, which included 5 mainshocks with Mw > 6 within 2 months and numerous smaller events that lasted for about 3 years (Rovida et al., 2020(Rovida et al., , 2022)).Thus, we investigate the stress loading history and the spatial distribution of stress for the faults that ruptured in these strong events.Finally, we compare, over similar timescales, the coseismic and cumulative Coulomb stress changes for three extensional fault systems in Italy with different geometry, the Calabria fault system (this work), and the Central and Southern Apennines, previously studied by Mildon et al. (2019) and Sgambato, Faure Walker, Mildon, and Roberts (2020) respectively, calculated in a comparable timescale.We investigate how different fault geometry and stress loading from underlying shear zones affect the stress distribution through time across the three studied fault systems.We discuss the results in terms of their implications for fault behavior and their relative influence on seismic hazard.

Geological Background
The Calabrian Arc is one of the most seismically active regions in Italy, linking the NW-SE oriented Southern Apennines chain with the E-W oriented Maghrebian chain in northern Sicily (L.Tortorici et al., 1995) (Figure 1).The Quaternary tectonic evolution of Calabria is characterized by extension with direction WNW-ESE, related to the opening of the Tyrrhenian Sea, and to the ongoing subduction of the oceanic crust of the Ionian Basin beneath Calabria.These processes are associated with the south-eastward migration of the Calabrian Arc (Malinverno & Ryan, 1986;Monaco & Tortorici, 1995;L. Tortorici et al., 1995) and the uplift and deformation of Quaternary marine deposits (Monaco & Tortorici, 2000;Westaway, 1993).Since the Middle Pleistocene, ESE-WNW directed extension prevails (L.Tortorici et al., 1995), creating the Siculo-Calabrian rift zone (Monaco et al., 1997;Monaco & Tortorici, 2000), a normal fault system extending 370 km along the coast of Sicily and the western side of the Calabria Arc (G.Tortorici et al., 2003).At the same time (∼0.7 Ma) strong regional uplift of the Calabrian Arc started, possibly as a response to the isostatic rebound due to the detachment of the Ionian subducting slab (Westaway, 1993).As a result of the uplift and Quaternary sea-level changes along the Tyrrhenian coast of Calabria, flights of Quaternary and Holocene marine terraces and wavecut platforms are presently exposed up to hundreds of meters above the present sea-level (Westaway, 1993).The preservation of marine terraces allowed the evaluation of an uplift rate of about ∼1 mm/year (Westaway, 1993).However, it is debated if the uplift is regional, and controlled by the subduction process, or local and related to the individual faults (Roberts et al., 2013).The other prominent structures characterizing the area are normal fault segments between 10 and 45 km long, strongly deforming Upper Pleistocene deposits (L.Tortorici et al., 1995) and associated to a young morphology, with escarpments hundreds of meters high, corresponding to a total throw accumulated in the last 700 ka, and often with well-developed triangular facets (Monaco & Tortorici, 2000).The slip-rates on these faults, measured from Holocene throws at the surface and displacement of raised marine terraces are rela tively well-constrained (Akinci et al., 2017;Aloisi et al., 2013;Catalano et al., 2008;Galli & Bosi, 2002;Galli et al., 2010;Meschis et al., 2018Meschis et al., , 2019;;Monaco et al., 1997;Monaco & Tortorici, 2000;Roberts et al., 2013;Roda-Boluda & Whittaker, 2017;L. Tortorici et al., 1995;G. Tortorici et al., 2003), although the crystalline rocks appear not to support scarp preservation, thus the lack of well-exposed fault scarps and hence fault dip measurements hinders precision in throw measurements (Roberts et al., 2013) (see Table S1 in Supporting Information S2 for slip-rate values and their uncertainties).The seismicity of the Calabrian Arc includes large crustal earthquakes distributed along the Tyrrhenian side of Calabria and the Ionian coast of Sicily (Monaco & Tortorici, 1995), and subcrustal earthquakes (depth 200-350 km), located along the inner side of the arc (L.Tortorici et al., 1995).The majority of the historical events recorded in Calabria occurred between 1638 and 1908 CE (Galli & Bosi, 2002).

Clustered Earthquake Sequences in Calabria
On the 5 February 1783 CE, a Mw 7.1 earthquake struck the southern part of Calabria and started a sequence that lasted for more than 3 years, including four other large events recorded between February and March 1783 CE.Some of these earthquakes were characterized by multiple faults rupturing at the same time; the whole sequence caused catastrophic damage and ∼35,000 deaths (Guidoboni et al., 2018(Guidoboni et al., , 2019) ) (Figure 3).Previous studies have identified the faults responsible for the 1783 CE events based on damage distribution (Jacques et al., 2001;L. Tortorici et al., 1995).The first and largest event occurred on 5 February 1783 CE (Mw 7.1), and it is hypothesized to have ruptured the Cittanova and Sant'Eufemia faults, with a total rupture length of ∼45 km (Figure 3a).The second event (Mw 6.3), on 6 February, struck the western coast of Calabria, rupturing 18 km of the Scilla fault (Figure 3b).We have modeled two scenarios for this earthquake, rupturing either the northern or the southern segment of the Scilla fault, Figure 3b shows the northern section rupture, while the alternative scenario is shown in Figure S1 in Supporting Information S1.On 7 February, a third mainshock of Mw 6.7 occurred further north, most probably rupturing the southern section of the Serre fault for ∼30 km (Figure 3c).The northernmost section of the Serre fault ruptured on 1 March 1783 CE with a Mw 5.9 (Figure 3d).The event that occurred on the 28 March was the last major event of the sequence (Mw 7.0), and it has been tentatively assigned to the Serre fault (Figure 3e).However, there is uncertainty on the source responsible for this event, due to the cumulative damage caused by the previous earthquakes.Other examples exist in Calabria of earthquakes separated by only a few days or months.The 27 March 1638 CE earthquake is hypothesized to consist of a series of two or three smaller events occurring on the Crati West, Lamezia and Serre faults (Galli & Bosi, 2003;Galli & Scionti, 2006), generating a cumulative magnitude Mw 7.0.A few months later, an event occurred on the Lakes fault (Galli & Bosi, 2003) (Figure 1).All these faults are in close proximity, suggesting the possibility of a triggered sequence.However, we are aware that this sequence occurs too early in the time considered to draw further conclusions.Other earthquakes smaller in magnitude occurred on the Crati West and Crati East faults in a possible example of a migrating sequence(s) in 1767, 1835, 1854, 1870CE (Galli & Bosi, 2003) and later in 1886, 1887 CE.

Modeling Coseismic CST From Historical Earthquakes
We modeled 23 faults in Calabria, including their mapped strike-variable geometry, and modeled the combined coseismic and interseismic loading to calculate the cumulative CST.Calculations are made in Coulomb 3.4 (Lin & Stein, 2004;Toda et al., 2005) using an elastic half-space model.
Following the approach of Mildon et al. (2016Mildon et al. ( , 2019) ) and Sgambato, Faure Walker, Mildon, & Roberts (2020), we calculate the changes in Coulomb stress transfer (∆CST) on faults with realistic strike-variable geometry.The ∆CST is defined as: where ∆τ is the shear stress change, ∆σ the normal stress change, ∆P is the pore pressure change and μ is the friction coefficient (Harris & Simpson, 1992).We use a value of 0.4 for the friction coefficient.Other works have investigated the effect of pore fluid pressure changes on Coulomb stress changes that are hypothesized to trigger aftershock sequences (e.g., Albano et al., 2017Albano et al., , 2019;;Miao et al., 2021).We have not considered this effect in our model (a) due to the lack of knowledge of the permeability properties around the faults, and thus lack of information on the magnitude and rate of pore pressure changes, and (b) because we are interested in seeing if we can explain the sequence without modeling the pore pressure changes.This is also discussed in Section 3.5.
The CST is calculated for 30 events with Mw > 5.5 occurring on 16 faults since 1509 CE using Coulomb 3.4 software (Lin & Stein, 2004;Toda et al., 2005) and a code developed by Mildon et al. (2016Mildon et al. ( , 2019) ) (Movie S1).The Coulomb stress changes in this paper are indicated in bar (1 bar = 0.1 MPa).To build 3D faults, the trace at the surface is divided into 1 km segments and faults are generated by discretizing the fault plane in ∼1 × 1 km elements.The surface fault traces are based on observation on geological maps, review of the literature and analysis of Google Earth satellite imagery and a DEM (Digital Elevation Model).The traces at the surface are used to project the fault to 15 km depth where the brittle-ductile transition is assumed to be.A single dip value is assigned along the fault and at depth, which is calculated as the mean value of the measurements of fault dip collected in the field, or if such data are unavailable it is taken from the literature.Where no information is available, a dip of 65° is assigned, which represents the average dip of the faults in the region.We use a rake value conventionally set to −90° for normal faults (Aki & Richards, 1980).
To calculate the CST from historical earthquakes occurring on the faults in Calabria, earthquakes with Mw > 5.5 have been selected from the Italian earthquake catalog CPTI15 (Rovida et al., 2020(Rovida et al., , 2022) ) and assigned to an active fault (Figure 1).Only events with Mw > 5.5 were selected, since the Italian historical catalog is thought in Supporting Information S1).The last event of the sequence (e) has been tentatively attributed to the Serre fault.Although the highest damage of the last event was localized south of the city of Catanzaro, this event was felt in most part of southern Italy.Due to the difficulties in the comparison and homogenization of the data, and the reliability of the sources, the epicenter location, like for most other events in this sequence, carries some uncertainties (CFTI15Med, Guidoboni et al., 2018Guidoboni et al., , 2019)).

10.1029/2023JB026496
7 of 26 to be complete for events above this magnitude after 1349 CE (Michetti et al., 1996).Each earthquake has been attributed to one or more faults.The criteria adopted to assign an event to a seismic source are based on published information; where palaeoseismological trenches or other sources of information were available, these have been used to assign the events to a seismogenic source (see Table S2 in Supporting Information S2 for details on references used to identify seismic sources or whether the source was assigned).In some cases (27% of the earthquakes) the distribution of the macroseismic intensity points (CPTI15, Rovida et al., 2020Rovida et al., , 2022) ) has been used to choose a source responsible for lesser-known events (see Table S2 in Supporting Information S2).In this case, shaking points are used as an indication of the fault being the source of the event considered, specifically when the highest damage is concentrated in the hangingwall of the normal faults.The modeled slip is consistent with the regional strain, since the normal active faults in Calabria are perpendicular to the NNW-SSE direction of the principal horizontal strain-rate axis calculated from geodetic data (Serpelloni et al., 2005), and consistent with the regional stress field (Montone & Mariucci, 2016).Note that the CST can also be modeled on optimally oriented faults using the regional stress orientation, a solution preferred by other authors (compare e.g., to Verdecchia & Carena, 2015).
Finally, a slip distribution is generated to model the earthquakes.The rupture length of each historical earthquake is calculated from the magnitude using the Wells and Coppersmith (1994) relationship that relates subsurface rupture length (RLD) to earthquake magnitude (M).The relationship is: a and b are coefficients and values are given for strike-slip, reverse and normal faults or for all faults (Wells & Coppersmith, 1994); the values used herein are for normal faults, with a = −1.88 and b = 0.50.Where the calculated rupture length for a given magnitude is too small for the whole fault to have ruptured, it is hypothesized that a partial rupture occurred along the fault.In the case of partial ruptures, we used information from the literature or the distribution of the shaking points as an indication of which section of the fault was more likely to have ruptured.Although the general location of each earthquake is well-known from macro-seismic records (Rovida et al., 2020(Rovida et al., , 2022)), and this constrains the assignment of each earthquake to a small number of scenarios, the lack of seismometer readings for most events means there is no information to create detailed slip distributions.Thus, we created a symmetrical concentric slip distribution, with the maximum slip in the center of the fault, 10% of the slip maximum at depth reaching the surface, and zero slip at the base and edges of the fault surface.This resembles the relationship between observed surface slip of 10 cm (EMERGEO Working Group, 2010;Rossetto et al., 2011;Vittori et al., 2011) and maximum slip values at depth of up to 0.8 m (Walters et al., 2009;Wilkinson et al., 2015) for the Mw 6.3 2009 L'Aquila Earthquake, Italy.This is consistent with previous approach (Mildon et al., 2016(Mildon et al., , 2019;;Robertson et al., 2020;Sgambato, Faure Walker, Mildon, & Roberts, 2020).In the case of a partial rupture, a concentric slip distribution is applied to the portion of fault that slips.For each of the 30 earthquakes modeled in Calabria, the maximum slip is iterated until the calculated magnitude matches the historical magnitude indicated for that event.Once the correct magnitude is reached and the slip distribution is modeled, the Coulomb stress changes are calculated in Coulomb 3.4 (Lin & Stein, 2004;Toda et al., 2005), onto receiver faults with realistic 3D geometry.Moreover, we assume that, following an earthquake, the Coulomb stress is reset to zero on the portions of faults that slip (Mildon et al., 2016(Mildon et al., , 2019;;Smith-Konter & Sandwell, 2009;Verdecchia et al., 2018).This is because, although the stress drops for the large earthquakes are typically between 10 and 100 bar (Nostro et al., 2001), due to lack of information of how stress drop in these historical earthquakes relates to finite stress, the actual stress drop on a fault following an earthquake is unknown.This will provide a minimum estimate of stress accumulated at a given time (Smith-Konter & Sandwell, 2009).

Modeling Shear Zones and Interseismic Loading
Shear zones are built from 15 to 24 km depth, as a continuation of the brittle portion of the faults (Figure 4).The lower bound is chosen at 24 km since it has been shown that deformation on normal faults in the Apennines is driven by viscous slip occurring at depth on localized shear-zones, with a lower limit for viscously deforming material at 710 ± 120 K, a temperature reached at ∼24 km depth (Cowie et al., 2013).The interseismic loading is modeled as annual Coulomb stress transferring from the shear zones slipping at depth to the brittle portions of the faults (Mildon et al., 2017(Mildon et al., , 2019(Mildon et al., , 2022;;Sgambato, Faure Walker, Mildon, & Roberts, 2020;Wedmore et al., 2017).The annual rate of slip is taken to be the slip-rate measured at the surface from fault scarps aged 15 ± 3 Kyr, as it is assumed that these are representative of the long-term slip-rate (Cowie et al., 2013), thus the 10.1029/2023JB026496 8 of 26 resulting stress transferred to the brittle fault reproduces the spatially variable slip along the faults.The cumulative slip we measure includes any postseismic afterslip localized on the fault, as it is resulting from the summed displacement of each of the earthquakes that would have occurred in the last 15 ± 3 Kyr.
For this work, due to the variable preservation of Holocene fault scarps, we use a combination of published slip rates from Holocene fault scarps and uplift of marine terraces, thus they represent the slip-rates over 15-700 Kyr (Akinci et al., 2017;Aloisi et al., 2013;Catalano et al., 2008;Galli & Bosi, 2002;Galli et al., 2010;Meschis et al., 2018Meschis et al., , 2019;;Monaco et al., 1997;Monaco & Tortorici, 2000;Roberts et al., 2013;Roda-Boluda & Whittaker, 2017;L. Tortorici et al., 1995;G. Tortorici et al., 2003).Where no information from the literature is available, we assume a maximum slip-rate of 0.2 mm/year, since higher slip-rates would produce fault scarps that would be recognized if present (Morewood & Roberts, 2000).Indeed, considering that a slip-rate of 0.2 mm/year would produce a 25 m high fault scarp for the 125 ka marine terraces, it would be surprising if we did not recognize such scarps.However, for a sliprate of 0.1 mm/year, it might not be possible to spot a 12.5 m fault scarp on a 125 ka marine terrace on a DEM.The dip of the shear zones is maintained equal to the dip of the brittle portion of the fault (Cowie et al., 2013;Mildon et al., 2017Mildon et al., , 2019Mildon et al., , 2022;;Sgambato, Faure Walker, Mildon, & Roberts, 2020;Wedmore et al., 2017).This is an approximation because dip may vary with depth, but we choose to model the shear zones with the same dip due to lack of detailed information at depth.Further, in support of our assumption, where high resolution deep seismic reflection profiles are available, these show localized shear zones beneath the upper crust, with dip similar to those in the upper crust (e.g., DRUM profile north of Scotland, Reston, 1990).To build the interseismic slip distribution, a 10.1029/2023JB026496 9 of 26 quasi-triangular slip profile is assumed, by linearly interpolating between all measurements of slip-rate measured at the surface and zero at the tips.We assume that short faults (length < 15 km) maintain an aspect ratio of 1 for strike dimension/dip dimension, a typical aspect ratio value associated to laterally restricted faults (Nicol et al., 1996).Thus, faults that are <15 km in length are not assigned shear zones extending beyond the depth of the seismogenic layer (Figure 4) (see also Mildon et al., 2017Mildon et al., , 2019;;Sgambato, Faure Walker, Mildon, & Roberts, 2020).

Cumulative Stress-Combining Coseismic and Interseismic CST
Once the coseismic stresses associated with each historical earthquake in the sequence and the interseismic stress for 1 year of loading have been calculated, these are combined to calculate the cumulative CST, also called pre-stress (Mildon et al., 2019) (Movie S2).In other words, the cumulative CST is the sum of coseismic and interseismic stress accumulated on each individual fault prior to an earthquake.Since the Italian catalog of historical seismicity is considered complete for earthquakes with Mw > 5.8 since 1349 CE (Michetti et al., 1996) and hence information on events prior to this are likely to be incomplete, we start our model of CST on the Calabria faults with zero stress at 1509 CE, the date of the first well-constrained earthquake in the region since 1349 CE, and an approach similar to previous works (Mildon et al., 2019;Sgambato, Faure Walker, Mildon, & Roberts, 2020).However, as it is recognized there would not be zero stress at 1509 CE, the start of the analysis is set after 1743 CE, to give the model some burn-in time away from the initial stress state, so the data have been analyzed after this time (Figure 5).The interseismic CST prior to each earthquake is calculated by multiplying the annual rate of slip and stress transfer by the number of years between two subsequent earthquakes.Then the cumulative CST before each earthquake is calculated by summing the coseismic effects due to earlier earthquakes on neighboring faults and the interseismic loading due to slip on the shear zone underlying the fault in question and other neighboring shear zones, prior to each historical event for each 1 × 1 km element constituting a fault plane.The stress accumulated on the part of a fault that ruptures is reduced to zero following each earthquake.

Analysis of Cumulative Coulomb Stress Changes
To compare the modeled histories of stress accumulation from 1509 CE to the present day produced by either summed coseismic and interseismic stress changes (cumulative CST), or by solely the coseismic stress changes, we calculate the mean cumulative CST and the percentage of faults with mean stress that is positive for all faults through time (Figure 5), and then specifically for the faults that ruptured in the 1783 CE sequence (Figure 6).Moreover, for some of the faults that ruptured in the largest earthquakes in the area, we show the mean cumulative CST prior to each historical earthquake (Figure 7), and the percentage of positively stressed elements prior to each event for the whole fault surface and for a sample of elements of the fault surface, located at 7 and 15 km depth (Figure 8).Following the approach of Sgambato, Faure Walker, Mildon, & Roberts, 2020, we compare the stress history for faults that ruptured across strike from each other, and faults that did not have any earthquake across strike, to investigate if different geometries influence the distribution and accumulation of coseismic and interseismic stresses.We show the mean cumulative CST and percentage of positively stressed element for faults in two areas of the Calabria fault system, one that includes the Cittanova (CIT), Serre (SER), Scilla (SCI), Sant'Eufemia (SE) faults, and one that includes the Messina-Taormina (MT) fault (Figures 7 and 8).The Cittanova, Scilla and Sant'Eufemia faults are located in the central section of the fault system, and are arranged across-strike from each other, while the Serre fault is located along strike to the Cittanova fault (Figure 1); these faults ruptured during the 5 major events of the 1783 CE sequence (Figure 3).The MT fault is located in the distal southern section of the fault system and is across strike from other faults that did not rupture in the time considered (Figure 1).

Model Limitations
Our modeling includes a number of assumptions, relating to the use of a simplified modeling approach for the slip distribution of historical earthquakes, the uncertainties in average displacement rates, the effect of transient postseismic stress changes, and to the choice of physical parameters that have been suggested to affect Coulomb stress changes.
One of the main assumptions of our model is that, since we do not have detailed information on the slip distribution from historical earthquakes, we assume a simple concentric slip distribution.The precise location and slip distribution of earthquakes can affect the pattern of the stresses across fault surfaces, especially close to the source fault.A heterogeneous distribution of stresses may be decisive in the occurrence of future earthquakes, and this should be considered when investigating interaction among multiple closely spaced faults (e.g., those in the Central Apennines fault system).For example, Mildon et al. (2019) show that although most of the earthquakes that occurred in the Central Apennines within the last few centuries were on faults that had a positive mean cumulative CST, the faults that had an earthquake with a negative mean cumulative CST were characterized by a heterogeneous pattern of stress, with some large positively stressed areas, where the earthquakes may have nucleated.We acknowledge that only taking into consideration the mean value of cumulative CST is limiting in terms of understanding the complicated patterns of stress that bring a fault to rupture, but we consider this a reasonable assumption in the absence of detailed slip distributions for historical events.
Another parameter to consider is the coefficient of friction (μ), which has been suggested to range between 0 and 0.85 (Verdecchia & Carena, 2015).We assume an average value of 0.4, which has been used in many other studies (e.g., King et al., 1994;Mildon et al., 2019Mildon et al., , 2022;;Stein et al., 1992;Verdecchia & Carena, 2015, 2016).Studies have shown that changes in the coefficient of friction have little impact on the magnitude of CST, with higher values of CST obtained for higher values of μ, and no impact on the spatial pattern of stress (King et al., 1994;Verdecchia & Carena, 2015, 2016).Thus, since we are interested in the relative changes in CST and not absolute stress values, we used an average value of 0.4.
We have not calculated the effects of pore fluid pressure changes on CST, although other studies have hypothesized a correlation of pore pressure changes with temporal and spatial distribution of aftershocks (e.g., Albano et al., 2017Albano et al., , 2019;;Miao et al., 2021), and sequences of earthquakes (e.g., Bosl & Nur, 2002;Noir et al., 1997;Vidale & Shearer, 2006).We think our approach is conservative, because we have no data on fluid pressure changes for historical earthquakes, and hence we investigate whether these earthquakes can be explained to an extent without invoking fluid pressure changes (e.g., Mildon et al., 2019).Moreover, the effects of pore pressure changes may modify the patterns of CST, but this effect is likely to be more pronounced in the region close to the source fault (Albano et al., 2017(Albano et al., , 2019)).Pore pressure changes could be important when evaluating sequences of triggered earthquakes, as they cause CST increases in the area surrounding the source fault immediately after a shock, however fluid overpressure tends to dissipate in a short time (60 days to 2 years) (e.g., Albano et al., 2017Albano et al., , 2019;;Miao et al., 2021).Nevertheless, our results show that the sequence of earthquakes in Calabria can be explained without the effects of fluid pressure changes.
As stated in Section 3.2, we assume that, if present, any postseismic afterslip would be included in the value of slip measured at the surface, so the interseismic creep we have modeled includes also postseismic afterslip.In other words, this means that contrary to other works (e.g., Bagge et al., 2019;Bagge & Hampel, 2017;Verdecchia & Carena, 2015, 2016;Verdecchia et al., 2018) we have not needed to separately model the stresses that occur due to postseismic viscoelastic relaxation after earthquakes.The transient postseismic viscoelastic stress changes are likely to be a second order effect, particularly in the few years following an earthquake, when postseismic afterslip is the dominant effect (e.g., Diao et al., 2014).However, postseismic viscoelastic relaxation may affect triggering of earthquakes over longer time periods (hundreds of years) and across large distances (tens of kilometers) (Verdecchia & Carena, 2015;Verdecchia et al., 2018).Postseismic CST may be comparable in magnitude to the interseismic CST calculated in the same time period (compare postseismic CST (Verdecchia et al., 2018) and interseismic CST (Mildon et al., 2017) calculated over 94 years before the 2009 L'Aquila earthquake, Italy), suggesting that interseismic and postseismic CST may be equally important for tectonic loading.
Fault displacement rates are sensitive to the time-period of observation, and this should be considered in our modeling, since the throw-rates we use to calculate yearly rates of interseismic stress loading for the Calabria faults are measured over a wide range of time periods (15-700 ka, see Section 3.2).Although displacement rates are relatively uniform on long temporal scales (>300 ka), over short time periods (<20 ka) they can vary greatly, by up to three orders of magnitude (Mouslopoulou et al., 2009;Nicol et al., 2009Nicol et al., , 2020)).However, for the study area, the agreement between geological rates and those measured with GPS geodesy and seismic moment release rate (e.g., Cowie et al., 2013;D'Agostino et al., 2011;Faure Walker et al., 2010, 2012;Meschis et al., 2022;Serpelloni et al., 2005Serpelloni et al., , 2010)), suggests that the geological observation periods may well be long enough to be compatible with the long-term average for elastic deformation rates and hence stress loading rates.Thus, we have relatively robust constraints on the rates needed to perform interseismic stress loading calculations.Information on the temporal window over which throw-rates are calculated are included in Table S1 in Supporting Information S2.
Finally, the interseismic CST might be influenced by the depth of the brittle portion of the faults, since it has been shown that stress accumulation rates are inversely proportional to the fault locking depth, therefore along-strike variations in stress may be caused by variations in locking depth (Smith-Konter & Sandwell, 2009).However, observations on hypocentral depths in the Apennines and Calabria suggest a relatively consistent locking depth at ∼15 km (Chiarabba et al., 2005;Chiaraluce et al., 2011), thus variations in interseismic stress due to depth variations would be minimal.

Interseismic Loading Rate in Calabria
Figure 4 shows the result of modeling 1 year of interseismic loading from slip on the shear-zones.Most of the faults in Calabria appear to have their entire brittle surface that is positively loaded by the stress transferred from shear zones, with a higher stress at the base of the fault.However, there are few exceptions: the Scilla fault (SCI), Sant'Eufemia fault (SE), Reggio Calabria fault (RC), and Southern Calabria fault (SC), have the upper portion of their brittle faults subjected to a stress reduction (blue areas across the fault surface), since they lie above shear zones of other neighboring faults.Despite this effect, all but one of the faults have a mean loading rate, calculated across the entire brittle fault surface, that is positive, with values ranging between 0.0007 bar (Coccorino fault) and 0.02 bar (Cittanova fault) per year; the only fault with a negative value for the mean interseismic loading rate is the Reggio Calabria fault (−0.001 bar), due to the presence of the MT shear zone underneath the brittle portion of the Reggio Calabria fault (Figure 4).

Effects of Cumulative (Coseismic Plus Interseismic) Versus Coseismic Coulomb Stress Changes in Calabria
Considering the cumulative CST (including the effects of both coseismic and interseismic CST), we find that all but one of the earthquakes after 1743 CE (and indeed since 1609 CE) occur on a fault with a positive mean cumulative CST (Figures 5a and 5b).In detail, the mean cumulative CST on the faults that ruptured historically was >0.5 bar in 91% of the cases since 1743 CE, and this is significant because the value of 0.5 bar is exceeding the typical hypothesized triggering threshold values (e.g., 0.5 bar King et al., 1994;0.1 bar Lorenzo-Martín et al., 2006; 0.2-1.0bar Toda et al., 1998;0.2 bar Verdecchia et al., 2018).For the one example where a negative cumulative stress is modeled, note that two alternative scenarios have been proposed and hence modeled for the 1905, Mw 6.9 earthquake: one that proposes a rupture of the Sant'Eufemia Gulf fault (SEG) with a cumulative mean CST of 0.98 bar on the whole fault (Loreto et al., 2013, see scenario iii in Figure 5a), and another which proposes the Vibo fault (VI) as the source responsible (Piantanesi & Tinti, 2002, see scenario iv in Figure 5a), which has a negative cumulative mean CST (−5.7 bar).This suggests that it is more likely that the 1905 CE earthquake occurred on the Sant'Eufemia Gulf fault.Thus, it may be that all earthquakes after 1743 CE occurred on faults with positive mean cumulative CST.Some of the strongest earthquakes in the area occurred on faults that had the highest mean cumulative CST of all studied faults prior to the event (red squares in Figure 5a).For example, before the first event in 1783 CE (Mw 7.1), the mean cumulative CST on the Cittanova fault was 5.41 bar; the 1908 CE event (Mw 7.1) occurred on the MT fault when this had a mean cumulative CST of 6.56 bar (Figure 5a; Tables S6 and S7 in Supporting Information S2).Other strong events with Mw > 6.5 either (a) occurred on faults whose geometry and location are not well-constrained, such as the 1905 CE, Mw 6.9 event, which may have occurred on the Sant'Eufemia Gulf fault (SEG), located offshore, or on the VI, and the 1832 CE, Mw 6.6 event, occurring on the Lakes South fault (Figure 1), or (b) they occurred on faults that ruptured multiple times in a short period, such as the Serre fault, which ruptured three times during the 1783 CE sequence.The percentage of faults that have a mean cumulative CST that is positive is 91% since 1743 CE (Figure 5b, Table S8 in Supporting Information S2).In contrast, the results obtained considering coseismic stresses only, produce a scenario where there is a less clear pattern on stress change and earthquake rupture.Only 48% of the faults have a mean coseismic CST >0.5 bar prior to an earthquake since 1743 CE (Figure 5c, Tables S3 and S4 in Supporting Information S2), and four earthquakes occurred on faults with negative values for the mean coseismic CST.The percentage of faults with positive mean coseismic CST through time is only 73 ± 5% (Figure 5d, Table S5 in Supporting Information S2).In this scenario some of the largest earthquakes in the area, such as the Mw 7.1 event on 5 February 1783, the Mw 6.9 1905, and the Mw 7.1 1908, occurred on faults with low mean coseismic CST, that is (a) 0.14 and 0.05 bar on the Cittanova and Sant'Eufemia faults respectively in 1783 (see also Figure 6 These results suggest that coseismic CST alone might not be dominant enough to explain the occurrence of the events in Calabria.In this fault system, characterized by faults oriented mainly along strike, few ruptures have occurred that caused stress reductions on across-strike faults, thus these faults should not be predominantly negatively stressed.Yet the faults that ruptured were often characterized by low values for the mean coseismic CST and were never the faults with the highest mean coseismic CST prior to an earthquake within the time period of the analysis (Figure 5c).
On the other hand, including interseismic loading in the calculations implies that most faults are positively stressed throughout the analyzed period (91% of faults have a positive cumulative CST since 1743 CE, Table S8 in Supporting Information S2), and most faults that ruptured in the same period have a mean cumulative CST above the threshold value (91% of faults with mean cumulative CST >0.5 bar).This reinforces the idea that the study of solely coseismic CST may not be adequate to study the full stress history of active fault systems when studied over multiple earthquake sequences (Mildon et al., 2019;Sgambato, Faure Walker, Mildon, & Roberts, 2020).

Coseismic and Cumulative CST During the 1783 CE Sequence in Calabria
Details of the coseismic CST and cumulative CST for the faults that ruptured in the 1783 CE sequence are shown in Figures 6a-6c; this allows investigation of whether the occurrence of this earthquake sequence is likely to be the result of triggering by coseismic stresses only or resulting from a combination of coseismic and interseismic stress loading.The first event of the sequence is suggested to have occurred on the Cittanova and Sant'Eufemia faults, followed by the Scilla fault in the second event, and the Serre fault which ruptured the southern section, northern section and possibly the whole surface during the third, fourth and fifth events respectively (Jacques et al., 2001;Sgambato, 2022) (Figure 6b).Some of these faults had already ruptured since the starting date of the analysis, while others had not experienced any earthquake before the 1783 CE sequence in the time period considered.This, along with the geometry and orientation of these faults, affects the state of stress across the fault surfaces.In this scenario, we consider that the Cittanova, Sant'Eufemia and Scilla faults presumably ruptured before the time considered in this analysis, but did not rupture since 1509 CE, and these faults have no other faults across strike, thus they accumulated stress undisturbed until 1509 CE.The Serre fault, on the contrary, ruptured four times before 1783 CE in the time considered: events associated to this fault are the 1626, 1638, 1659 and 1743 CE.In terms of cumulative CST the Cittanova fault had the highest mean cumulative CST value at depth before the first mainshock in 1783 CE (Figure 6a), this is likely because the Cittanova fault has the highest slip-rate among these faults (1.6 mm/year, compared to 0.7-0.8mm/year for the other faults; Table S1 in Supporting Information S2) and therefore the rate of interseismic loading is higher.This is shown in Figure 6a, where the Cittanova fault has the highest mean cumulative CST among all faults prior to the first earthquake in the sequence.All faults that ruptured in 1783 CE had a positive mean cumulative CST >1 bar before the first earthquake: 5.41 bar for the Cittanova fault, 1.93 bar for the Serre fault, 1.68 bar for the Scilla fault, and 1.36 bar for Sant'Eufemia fault (Figure 6a).Earthquakes 1 and 2 increase the mean cumulative CST on the Serre fault to a value of 3.01 bar before the third earthquake (Figure 6a).The Serre fault then re-ruptured with a smaller magnitude (Mw 5.9) during earthquake 4, when its mean cumulative CST value was 1.2 bar.The last earthquake, which is here hypothesized to have ruptured the Serre fault, occurred when the fault had a mean cumulative CST value of 1.00 bar (Figure 6a).
In contrast, if only coseismic stress changes are considered, the values are <0.5 bar before the first event of the sequence for all faults, except the Serre fault: that is 0.92 bar for the Serre fault, 0.14 bar for the Cittanova fault, 0.06 bar for the Scilla faults, 0.05 bar for the Sant'Eufemia fault (Figure 6b).The Scilla fault had a negative mean coseismic CST value (−1.24 bar) before rupturing in the second event of the sequence (Figure 6b).The Serre fault ruptured with the third event when it had accumulated a mean coseismic CST of ∼2 bar (Figure 6b).The Serre fault had a mean coseismic CST value of 1.07 and 0.95 bar before the fourth and fifth event respectively (Figure 6b), a value similar to that obtained for the mean cumulative CST (Figure 6a).More information is provided by the pattern of CST across faults that ruptured in 1783 CE: (a) the Cittanova fault is highly stressed before earthquake 1 (Figure 6c (i)); (b) although the first two earthquakes increased the stress on the Serre fault which subsequently ruptured in earthquake 3, the Serre fault was already highly stressed before earthquake 1 (Figure 6c (i, ii)); (c) the third and fourth earthquakes occurring on the Serre fault left a patch ∼5 km long in the central section of the fault that is highly stressed (>∼3 bar) (Figure 6c (iii-v)).Since there is no detailed information on rupture length and stress drop for these historical earthquakes, they are modeled to rupture across the whole depth of the brittle fault surface; therefore, it is not possible to know in detail the pattern of positive stress that might have been still present across the Serre fault surface after earthquake 3 and 4 of the sequence (Figure 6c (v)).A plausible hypothesis is that a patch of high stress like the one on the Serre fault before earthquake 5 might have contributed to starting the last rupture, since there are no regions of negative stress across the fault that might represent an impediment in the rupture propagation across the fault surface (e.g., Mildon et al., 2017;Steacy & McCloskey, 1998), allowing the occurrence of a large magnitude earthquake (Figure 6c (v)).
In summary, it appears that if solely coseismic stresses are considered, then earthquakes appear to have occurred on faults with very low or even negative stress loading values in the 1783 CE sequence.In contrast, if interseismic and coseismic stresses are combined to calculate the full CST then earthquakes occurred on faults with positive stresses, with values above postulated triggering thresholds.These results suggest that even for sequences of earthquakes occurring in a short period of time, the interseismic loading should not be ignored.

Fault Interaction and Stress Heterogeneity in Calabria
The analysis of the mean cumulative CST shows that some faults in Calabria that ruptured in the largest earthquakes also had the highest value for the mean cumulative CST in the region at the time of their rupture (Figure 5a).
To further investigate the stress loading patterns and their link to fault geometry, we analyzed the stress loading history of these faults, and analyzed the stress changes occurring at the center of the brittle portion of the fault (7 km depth) and at the boundary between brittle faults and their underlying shear-zones (15 km depth) (Figure 7, Table S9 in Supporting Information S2).We focused on two different sub-sections of the Calabria fault system, that is the section including the Cittanova, Sant'Eufemia, Scilla, Serre faults, which ruptured in the 1783 CE sequence, and a section that includes the MT fault, which ruptured with the highest magnitude earthquake in the time considered (Mw 7.1 Messina earthquake in 1908 CE).
The overall feature that emerges is that before an earthquake occurs, the mean CST averaged across the whole fault surface increases over time on all the faults listed above (Figure 7a).The Cittanova, Sant'Eufemia, Scilla, and Serre faults show decreases or increases in stress, when an earthquake occurs across strike, interrupting the long-term loading (Figures 7a-d).Superimposed on this steady accumulation of cumulative CST for the whole fault surface we observe both positive and negative stress changes produced by interactions with earthquakes on neighboring faults.For example, superimposed on the steady-increase on the MT fault, we see an increase of ∼0.4 bar due to the events in 1783 CE, which occurred along strike of the fault (Figures 7a-e).Other examples are clear when analyzing a sample of elements at 7 km depth, where the effects of across-strike interaction exist for the Cittanova, Sant'Eufemia, Scilla, Serre faults, and the increase in stress due to along-strike interaction is visible at 7 km depth on the MT fault (Figures 7b-e).At 15 km depth, the effect of the coseismic CST is overwhelmed by the interseismic loading in the vicinity of the shear zones beneath the faults (Figure 7c), the exception being the The Cittanova, Sant'Eufemia, Scilla, and Serre faults show a complicated pattern of stress loading, characterized by stress changes due to events occurring across strike.
No changes due to earthquakes across strike are observed for the MT fault.See Figure S3 in Supporting Information S1 for mean cumulative CST across Sant'Eufemia fault for the scenario in which the Scilla fault ruptures its southern section.
Scilla fault, where the stress drops at 15 km depth due to the 1783 Mw 7.1 event on the Cittanova and Sant'Eufemia faults (Figure 7c).No changes in stress due to earthquakes along or across strike are observed at 15 km depth on the MT fault because the pattern is dominated by interseismic stresses induced by the underlying shear zone (Figures 7c-e).
To study the stress spatial distribution across the same set of faults, we calculated the percentage of positively stressed elements before each earthquake for the whole fault surface and at specific depths (7 km, 15 km, Figure 8, Table S9 in Supporting Information S2).There is significant variability in values for the mean and standard deviation for the percentage of positively stressed elements on individual faults calculated over the historical sequence.We note that there is a wider range for these values where multiple faults are across strike, for example, for the Cittanova, Sant'Eufemia and Scilla faults compared to the MT fault (Figure 8a).The value is 83 ± 18% (mean and standard deviation) for the Cittanova fault, 78 ± 8% for the Sant'Eufemia fault, 76 ± 13% for the Scilla fault and 93 ± 15% for the Serre fault (Figure 8a).Whereas, since there are no ruptures across strike in the time considered, the MT fault shows a constant value of 99% of elements that have positive stress averaged over 1609-1928 CE, where averaged across the whole fault (Figure 8d).Similar patterns are observed at 7 km depth for the Cittanova (80 ± 21%), Sant'Eufemia (85 ± 13%), Scilla (85 ± 19%), Serre faults (92 ± 17%), but the value is 100% for the MT fault (Figure 8b).At 15 km depth the percentage is generally constant through time for the Cittanova, Sant'Eufemia, Scilla faults, because the coseismic stress changes are, in general, overwhelmed by stress changes due to loading from the underlying shear zone; with the exception of the events occurring during the 1783 CE sequence, the three faults show a percentage of positive elements that is 100% over time (Figure 8c).The percentage at 15 km is 100% through time for the MT fault.The Serre fault shows a higher variability at 15 km due to across strike interaction, such as the effect of the stress decrease caused by the 1905 CE earthquake which in the model is shown as occurring on the Sant'Eufemia Gulf fault (SEG, see Figure 1), across-strike from the Serre fault (Figures 8c-d).These results indicate that faults that rupture across-strike cause more variability in the stress loading, compared to faults that did not experience any stress reduction due to ruptures across strike, such as the MT fault.This is also shown as a variable percentage of fault surface that is positively stressed before an earthquake, since the faults investigated ruptured when between 30% and 100% of faults elements comprising the entire surface were positively stressed (Figure 8a).The earthquake on the Cittanova fault occurred when 100% of its surface was positively stressed (Figures 8a-a).The values are 76%-77% for Sant'Eufemia fault, when it ruptured in 1783 CE and 1894 CE (Figures 8a and b).The Scilla fault had 49% of its whole surface that was positively stressed before the second event of the 1783 sequence (Figures 8a-c).A wider range of values is observed for the Serre fault, with 30%-100% of the faults element that had a positive CST at the time of historical earthquakes (Figures 8a-d).This wide range created because this fault ruptured three times in a row and we assume that the CST on ruptured fault sections reduces to zero, which would not be counted as positive CST in our analysis.The value for the MT fault, with no faults across strike was 99% for the 1908 CE earthquake (Figures 8a-e).

Comparison of Coulomb Stress Models in the Central Apennines, Southern Apennines and Calabria Fault Systems
The results presented above show that for the interseismic periods of individual faults, rupture of neighboring faults produce stress reductions or increases that are superimposed on the steady-rate increase in stress produced by its underlying shear zone.
To further investigate the influence of fault and shear zone system geometry on the stress loading, we compare the study of CST across faults in Calabria with previous studies of the Central and Southern Italian Apennines, which have different fault system geometries and similarly relatively long and complete records of past earthquakes.Figure 9 shows a comparison between the mean cumulative CST (Figure 9a) and the mean coseismic CST (Figure 9b) calculated on all faults before each historical event in each region.The percentage of faults that have a positive mean cumulative CST (Figure 9c) and positive mean coseismic CST (Figure 9d) through time for each region are also shown.The values of mean cumulative CST on all faults are generally higher in the Central Apennines and Calabria compared to the Southern Apennines (up to ∼12 and ∼11 bar, and ∼3 bar respectively) due to the relatively higher 15 ± 3 Kyr slip-rates measured during fieldwork and hence input into our modeling (Figures 9a-c).After the initial burn-in phase for the Central Apennines model, all earthquakes occurred on faults that are not characterized by the highest value of mean cumulative CST among all faults in the system, except in one case (1915 CE, Fucino fault, Mw 7.0; Figure 9a); 79% of the faults that ruptured had a positive mean cumulative CST (Mildon et al., 2019), so some events occurred on faults that had a mean CST value that is negative (Figure 9a), although the authors noticed that 5 out of 6 of the faults with negative mean cumulative CST had some parts of their surface that was positively stressed before the rupture (Mildon et al., 2019).In the Southern Apennines, some events occurred on faults with the highest mean cumulative CST in the region at the time of the earthquake (e.g., 1694 and 1980 CE, Irpinia fault, Mw 6.8 and 6.9); 94% of the faults that ruptured had a positive mean cumulative CST (Figure 9a and b; Sgambato, Faure Walker, Mildon, & Roberts, 2020).The fault system in Calabria is characterized by many earthquakes in a short time (22 large magnitude events over 185 years since the starting date of the analysis), all occurring on faults that are positively stressed, except for the 1905 CE event in the alternative scenario ii in Figures 9a-c, which would bring the percentage of faults rupturing with a positive mean cumulative CST to 96% if this was the chosen scenario.In Calabria, like the Southern Apennines, some earthquakes occurred on faults with the highest cumulative CST value before an earthquake, such as the Mw 7.1 1783 CE, and the Mw 7.1 1908 CE (Figures 9a-c).Furthermore, the average value for the mean cumulative CST on faults that ruptured is highest in Calabria, with an average of 2.16 bar, with 1.09 bar in the Central Apennines, and 0.90 bar in the Southern Apennines.This is due to the effect of stress decrease due to ruptures occurring across strike, an effect that is less pronounced in Calabria than in the Central Apennines.
A comparison of the percentages of faults with a positive mean cumulative CST through time for the three areas shows that, after the initial burn-in period, most of the faults in Calabria are constantly positively stressed (91%), compared to 73% of faults in the Central Apennines (Mildon et al., 2019) and 70% of faults in the Southern Apennines (Figure 9c).This difference cannot be explained by rate of interseismic loading because the Central Apennines faults also have generally high throw rates (up to 1.56 mm/year), therefore the same behavior would be expected.The difference can be explained by the presence of multiple faults interacting across-strike, reducing the stress on some of the faults in the system, and thus reducing the percentage of faults that are positively stressed through time (Figure 9c).
When considering coseismic stress only (Figure 9b) it is shown that although all earthquakes occurred on faults with low mean coseismic CST, relative to the maximum coseismic CST on faults in the region at each time step, there are some differences between the three regions.The percentage of studied earthquakes that occurred on faults with mean coseismic CST >0.5 bar is 0%, 35%, and 42%, for the Central Apennines, Southern Apennines, and Calabria respectively (Figures 9a-c), suggesting a large influence of across strike interaction in the Central Apennines faults, producing stress reduction, compared to the other studied regions.This can be considered as evidence that in regions with complicated geometry and many faults across strike there is a greater need to investigate the cumulative CST to understand which faults are close to rupture.In Figure 9d, the percentage of faults with mean coseismic CST that is positive is plotted through time and compared for the three regions.While for the Southern Apennines and Calabria the number of faults that have a positive mean coseismic CST is variable, and generally lower than the number of faults with positive mean cumulative CST, this is not the case for the Central Apennines faults, which show at times a higher percentage of faults with positive mean coseismic CST than that calculated with mean cumulative CST (compare Figures 9c  and 9d).This is due to the closely spaced faults across strike in the Central Apennines, so that when shear zones are modeled most of the faults show positive stress transferred to the base of the shear zone, but some faults experience negative transferred stress, since they lie just above the shear zone of another across-strike fault.Thus, when calculating the percentage of faults that have a positive mean cumulative CST through time, some faults will have a portion of their brittle fault that is negatively stressed, but this effect is not present when calculating the coseismic stress only, since the shear zones are not modeled.This explains the difference between the percentages in Figures 9c and 9d.In the Southern Apennines and Calabria the interseismic loading generally increases the stress on all faults through time, and since few earthquakes across strike occur, most of the faults remain constantly positively stressed through time.In the Central Apennines many faults are across strike from each other, and therefore it is inevitable that some faults will experience stress increase, and some will experience stress reduction, therefore producing a larger range.
Finally, an analysis of the percentage of positively stressed elements across the surface of faults that ruptured in the three regions allows us to investigate how heterogeneous the stress is across the faults that are close to rupture, when the cumulative CST is calculated (Figure 10).For the time period after the burn-in times for the three areas, it appears that faults in the Southern Apennines tend to rupture when most of their surface is positively stressed (Figure 10b).Indeed, the range of values for the percentage of positively stressed elements are more clustered at higher percentages in the Southern Apennines (Figure 10b), whereas there is a greater spread of data for the Central Apennines and Calabria (Figures 10a and 10c).Faults that ruptured with the highest cumulative CST value were also characterized by a high percentage of positive elements across the fault surface, such as the 1694 and 1980 CE in the Southern Apennines, with a value of 95% (Figure 10b), and the 1783 and 1908 CE in Calabria with 100% and 99% respectively (Figure 10c).These results suggest that, where deformation is more localized, as is in the Southern Apennines and Calabria, the geometry influences the distribution of the stresses across the faults, and in turn can possibly determine a higher magnitude of CST on faults prior to a rupture.

Discussion
Whether fault interactions influence variability in earthquake recurrence is a key question for seismic hazard assessment.We investigate this question by studying three active fault systems, which share similar tectonics but have different geometry, and where there is a good knowledge of past seismicity, fault geometry and kinematics, allowing us to analyze fault interactions over the last few hundred years.Our analyses highlight the combination of effects of coseismic loading and interseismic loading with fault system geometry that influence the stress loading history of the earthquake sequences in Calabria and in the Central and Southern Apennines.
We suggest that the geometry of the fault system in Calabria, specifically a lack of multiple faults across strike from each other that have ruptured during the analyzed period, provides an environment in which high interseismic loading rates dominate the stress loading, overwhelming the effect of along-and across-strike stress changes.One striking example is represented by the MT fault, where the main control on the stress loading is the interseismic stress from the underlying shear zone, which led the fault to rupture in 1908 CE with a Mw 7.1 (Figure 7e).Nonetheless, some faults in Calabria are influenced by coseismic across strike interactions, and these show a higher complexity in their stress loading history (Figures 7a-7d) and a distribution of CST values that is heterogeneous across the fault surface (Figures 8a-8d), in agreement with previous work (Sgambato, Faure Walker, Mildon, & Roberts, 2020).These results show that there is a difference in the way stress is accumulated through time among the studied faults, with evidence of across-strike interactions on multiple parallel faults creating complex heterogeneous stress patterns, and negligible along-strike interaction when this is compared to the influence of regional stress loading.Overall, stress loading history and spatial distribution of the stress across a fault surface are dependent on the geometry of the overall fault system because this determines whether a particular fault has fault/shear-zone structures located across strike or not.Thus, variability in stress build up at the scale of individual faults is related to their location and orientation within a fault system and the overall fault system geometry, in agreement with what previously suggested by Cowie et al. (2012) for numerical models.
At the scale of the fault system, the number of faults across strike influences the mean values of cumulative CST when faults rupture and whether the values are positive or negative.This is shown as a high number of faults that ruptured in Calabria with a positive value for the cumulative CST before an earthquake occurred on them in the time considered, that is 91% since the starting date of the analysis (Figures 3b and 9c).In contrast, this value is 73% for the Central Apennines (Mildon et al., 2019) and 70% for the Southern Apennines (Figure 9c).Although high values of slip-rates are found for both faults in Calabria and the Central Apennines, our modeling suggests that the geometry of the fault systems differentiates their stress loading patterns.The difference is that in Calabria, the high slip-rates of faults in the area and the fault system geometry constitute favorable conditions to the accumulation of a uniformly high level of stress across more faults, and this is associated with the occurrence of several destructive earthquakes in Calabria in a short time (Figure 9a), and some propagating triggered sequences (see Section 2.1).While in the Central Apennines the regional loading rate is distributed in places between ∼10 faults across strike, in Calabria this is accommodated by only 2-3 faults across strike (Figure 8).Thus, faults in Calabria have a high proportion of the regional rate accommodated on each fault, and the loading is not reduced by across-strike interactions to the same extent as the Central Apennines.In the Central Apennines the stress loading of faults is more spatially heterogenous due to the preponderance of across-strike coseismic interactions.This suggests that, with similar throw rates and slip rates, the geometry of the fault system is a major influence on the stress loading of faults.
Our results suggest that the propagating earthquake sequence in 1783 CE may have occurred in relation to the cumulative CST across the faults, rather than being an effect of coseismic stresses only, as previously suggested (Jacques et al., 2001).The 1783 CE earthquake sequence started with a rupture along the Cittanova fault, when this was characterized by the highest mean cumulative CST of all faults, 5.41 bar (Figure 6a), as opposed to a mean value of 0.14 bar if only coseismic stresses are considered (Figure 6b), and all faults that ruptured during this sequence were all positively stressed with a mean cumulative CST across the entire fault surfaces between 0.5-6.0bar before the first event of the sequence in 1783 CE (Figure 6).For this reason, calculations of the cumulative CST should not be ignored when studying triggered sequences of earthquakes that occur in a short period of time.Moreover, large events may occur where stress is relatively uniform across the fault, since rupture propagation is not prevented by low stress heterogeneities across the fault surface, similar to the ideas previously proposed by Steacy and McCloskey (1998).Perhaps the occurrence of many events in a short time in Calabria is the result of this behavior that cause the faults to be highly stressed and with less heterogeneity in stress, so that when an earthquake occurs it cascades along the system.Although there are examples of large earthquakes occurring within areas of negative stress across the fault (e.g., the 30 October 2016 Mt Vettore earthquake; Mildon et al., 2017Mildon et al., , 2019)), and this has been observed in other tectonic regimes (e.g., strike-slip fault model, Duan & Oglesby, 2005, 2006), where heterogeneous stress is generated at the regions of fault bends or fault overlaps and has a significant effect on rupture initiation and propagation (Duan & Oglesby, 2005, 2006).However, other processes that have not been considered in this work may promote triggering of earthquakes, such as dynamic stress triggering (e.g., Belardinelli et al., 2003;Gonzalez-Huizar & Velasco, 2011;Velasco et al., 2008), and the role of fluids (e.g., Bosl & Nur, 2002;Miao et al., 2021;Noir et al., 1997).
Overall, these results show that regions of high geometric complexity are characterized by a higher variability in stress loading both spatially along individual fault surfaces and through time.This has potential implications for different aspects of seismic hazard studies.Alterations to the probability of occurrence for future events can be induced by CST that are thought to delay or advance the next earthquakes (Gomberg et al., 1998;Stein et al., 1997;Toda et al., 1998).The effects of static CST are considered in some probabilistic seismic hazard models by changing the time-dependent conditional probability of earthquake occurrence (e.g., Iacoletti et al., 2021;Pace et al., 2014;Stein, 1999;Stein et al., 1997;Toda et al., 1998;Verdecchia et al., 2019).In these models the interseismic loading is indirectly included in the calculation of the earthquake probability, by increasing the elapsed time or reducing the mean recurrence interval by the value of calculated CST divided by the tectonic loading rate (e.g., Stein et al., 1997;Toda et al., 1998;Verdecchia et al., 2019).Many studies have included coseismic, interseismic and postseismic stress in calculation of CST (e.g., Albano et al., 2017Albano et al., , 2019;;Bagge et al., 2019;Bagge & Hampel, 2017;Freed & Lin, 2001;Luo & Liu, 2010;Verdecchia & Carena, 2015, 2016;Verdecchia et al., 2018Verdecchia et al., , 2019) ) and demonstrated the influence of Coulomb stress changes on time dependent and time-independent probability calculations (Verdecchia et al., 2019), and we follow their lead.However, the spatial variability of coseismic, interseismic and postseismic stress changes due to fault system geometry is often not considered.The results in this paper show that fault and shear zone geometry, in other words the number of faults and their shear zones across strike, controls the style of cumulative CST, and therefore might control recurrence intervals.Thus, we suggest that spatial variability of cumulative stress changes, and thus the variability in recurrence intervals due to fault/shear zone geometry, should be considered in time-dependent PSHA.
Knowledge of the current state of stress on faults is fundamental to assess whether future ruptures will occur (Cowie et al., 2012;Marzocchi et al., 2009;Mildon et al., 2019;Sgambato, Faure Walker, Mildon, & Roberts, 2020), since differences in fault system geometry may be associated with different recurrence time distributions, depending on the complexity of the fault system (Cowie et al., 2012;Marzocchi et al., 2009;Visini & Pace, 2014).In general, faults with high stress rates have shorter recurrence intervals, and vice versa, faults with low stress accumulation rates have longer recurrence intervals (e.g., Smith-Konter & Sandwell, 2009).Our results, however, show that for faults with similar stress loading rates, the geometry of the fault system may have a decisive role in governing the variability of earthquake recurrence.Moreover, our results suggest that the fault and shear zone system geometry is responsible for the spatial and temporal distribution and magnitude of the Coulomb stress changes across faults.This could influence where nucleation occurs and determine rupture size and continuity of rupture along a fault or multiple faults.This is particularly relevant for evaluating lengths and probabilities of ruptures, which are a fundamental input for seismic hazard analyses (Biasi & Wesnousky, 2016).Thus, we suggest that it could be beneficial to the calculation of the probability of occurrence of earthquakes to incorporate region-specific constraints that aim to characterize the complexity of fault and shear zone system geometry since, as shown here, this represents a major influence on fault interactions through CST.Geometric complexity of fault and shear zone systems is one of the main controls on stress loading and thus may control how earthquakes are distributed in space and time; more work is needed to further explore the application of these results to seismic hazard models, possibly involving determination of a quantitative relationship between CST, fault system geometry and earthquake probability.

Conclusions
Calculations of cumulative CST from 30 earthquakes that occurred in Calabria between 1509-1928 CE show that all but one of the events occurred on faults that had a positive mean cumulative CST value.In two cases, the faults that ruptured had the highest cumulative stress value of the faults in the region prior to the earthquakes, and these were among the strongest events recorded in the area (Cittanova fault, Mw 7.1 in 1783 CE; Messina-Taormina fault, Mw 7.1 in 1908 CE).The faults that ruptured were highly stressed by interseismic loading from the underlying shear zone, with little influence from across-and along-strike stress transfer.Results show that the effects of along and across-strike interactions can be overwhelmed by high interseismic loading from high slip-rates, as observed for the Messina-Taormina fault.The mean stress history for faults that ruptured during the 1783 CE sequence highlight the importance of considering all known historical earthquakes and the interseismic stress, in order to fully analyze a triggered seismic sequence.Only by calculating the cumulative CST can we explain the occurrence of strong earthquakes and identify highly stressed faults that might be close to failure.A comparison of CST changes over the last few centuries in Calabria, Central and Southern Apennines, shows that in Calabria a larger number of faults were positively stressed over time in comparison to the other areas, a high number of strong earthquakes occurred in a relatively short time, and large magnitude earthquakes appear to have occurred on highly stressed faults.This behavior differentiates faults in Calabria from the Central Apennines and Southern Apennines faults, and this can be attributed to their different fault system geometry.Overall, these results highlight the predominance of interseismic stress driven by regional tectonics in Calabria, and the importance of the fault system geometry and long-term slip-rates in controlling the relative effects of coseismic and interseismic stresses, which in turn might affect the occurrence of earthquakes.

Figure 2 .
Figure 2. Map of the Italian active fault regions and their associated historical earthquakes.(a) Locations of regions within Italy; (b) Calabria (this work); (c) Southern Apennines (modified from Sgambato, Faure Walker, Mildon, & Roberts, 2020); (d) Central Apennines (modified from Mildon et al., 2019).Dashed lines in (c) and (d) indicate an area of overlap between the Central Apennines and Southern Apennines.Faults and associated earthquakes in this area were included in both analyses of Coulomb stress transfer in the respective published works (see Figure 9); the analyses are unchanged in this work.The maps show the difference in geometry and slip-rates in the three fault systems, with many faults across strike in the Central Apennines, few faults across strike in the Southern Apennines and Calabria.

Figure 3 .
Figure 3. Panels showing the 1783 CE earthquake sequence in Calabria.Red lines are active faults; CIT, Cittanova fault; SE, Sant'Eufemia fault; SCI, Scilla fault; SER, Serre fault.The sections of faults that ruptured in each event are highlighted in yellow, the sections of faults that rupture next are in blue.The epicenters of the historical earthquakes (CFTI15Med, Guidoboni et al., 2018, 2019) are indicated with stars.Events (a)-(d) have been assigned to the faults following Jacques et al. (2001).Two scenarios have been modeled for the 6 February 1783 earthquake, rupturing either the northern segment of the Scilla fault (b) or its southern segment (see FiguresS1 and S2in Supporting Information S1).The last event of the sequence (e) has been tentatively attributed to the Serre fault.Although the highest damage of the last event was localized south of the city of Catanzaro, this event was felt in most part of southern Italy.Due to the difficulties in the comparison and homogenization of the data, and the reliability of the sources, the epicenter location, like for most other events in this sequence, carries some uncertainties (CFTI15Med,Guidoboni et al., 2018Guidoboni et al., , 2019)).

Figure 4 .
Figure 4. Annual interseismic loading rates calculated on 3D strike-variable active faults in Calabria.Coulomb stress transfer (CST) represents 1 year of interseismic loading rate transferred from shear zones to the brittle portion; stress-loading onto the brittle portion of the fault and annual slip distribution on the underlying shear zones are shown.Most of the faults are wholly positively loaded by the stress transferred from shear zones, except the Scilla fault (SCI), Sant'Eufemia fault (SE), Reggio Calabria fault (RC), and Southern Calabria fault (SC), since they lie above shear zones of other neighboring faults, resulting in a stress reduction on portions of their brittle fault surfaces.
); (b) −1.38 bar or −8.91 bar for one of the two possible scenarios for faults that ruptured in 1905, the VI or Sant'Eufemia Gulf faults respectively; and (c) 0.35 bar on the MT fault before the 1908 earthquake.

Figure 5 .
Figure 5. Analysis of cumulative and coseismic Coulomb stress transfer (CST).Magnitudes for faults that ruptured are indicated with color coded symbols.Values of CST for all other faults are in gray.(a) Mean cumulative CST (coseismic plus interseismic stresses) calculated for all faults in Calabria prior to each historical earthquake.The faults that ruptured show positive mean cumulative CST since 1609 CE, except for one of the alternative scenarios for the 1905 CE earthquake, that possibly occurred on the Vibo fault (VI) (red square indicated with d).Magnitudes are indicated with color coded symbols.(i) and (ii) represent the mean cumulative CST on the Sant'Eufemia fault before the 1894 CE earthquake, in the case of the Scilla fault rupturing either its northern (i) or southern segment (ii) on 6 February 1783.(iii) and (iv) are alternative scenarios for the 1905 earthquake occurring either on the Sant'Eufemia Gulf fault (iii) or VI (iv).(b) Percentage of faults that have mean cumulative CST that is positive through time.(c) Mean coseismic CST calculated for all the faults prior to each historical earthquake.(d) Percentage of faults in Calabria that have mean coseismic Coulomb stress that is positive.The cumulative and coseismic CST on faults during the 1783 CE sequence is shown in detail in Figure 6.

Figure 6 .
Figure 6.Analysis of the 1783 CE seismic sequence.(a) Mean cumulative Coulomb stress transfer (CST) on faults that ruptured in the 1783 CE sequence.All the faults that ruptured in the sequence have a positive mean cumulative CST prior to the earthquake.The sequence started on the Cittanova fault (CIT) that had the highest cumulative mean stress among all faults in the system at the time of the event, and this was one of the strongest earthquakes recorded in the study region.(b) Mean coseismic CST on faults that ruptured in the 1783 CE sequence.(c) Spatial distribution of mean cumulative CST before each earthquake in the 1783 CE sequence.Green lines are fault surface traces.Faults that rupture next are shown in purple and indicated with black arrows.The Cittanova fault ruptured contemporaneously with the Sant'Eufemia fault in the first event (i); the second earthquake occurred across strike, on the Scilla fault (ii), here it is shown the northern segment of the fault rupturing, an alternative scenario is shown in FiguresS1 and S2in Supporting Information S1.Then the sequence continued along strike on the Serre fault, which ruptured first its southern section (iii), then its northern section (iv), and then possibly its whole length (v).

Figure 7 .
Figure 7. Cumulative Coulomb stress transfer (CST) over time calculated for the Cittanova, Sant'Eufemia, Scilla, Serre and Messina-Taormina (MT) faults.(a) Mean cumulative CST across the whole fault surface; (b) mean cumulative CST calculated at 7 km depth; (c) mean cumulative CST calculated at 15 km depth.Orange circles indicate partial ruptures, red circles are total ruptures of the whole fault surface; blue circles represent the imposed condition of zero stress on fault, immediately after a total rupture.Events marked (i)-(v) indicate events occurring in the same year.Dashed boxes highlight stress changes due to rupture occurring across or along strike.The Cittanova, Sant'Eufemia, Scilla, and Serre faults show a complicated pattern of stress loading, characterized by stress changes due to events occurring across strike.No changes due to earthquakes across strike are observed for the MT fault.See FigureS3in Supporting Information S1 for mean cumulative CST across Sant'Eufemia fault for the scenario in which the Scilla fault ruptures its southern section.

Figure 8 .
Figure 8. Analysis of cumulative Coulomb stress transfer heterogeneity on the Cittanova, Sant'Eufemia, Scilla, Serre and Messina-Taormina (MT) faults.(a) Percentage of positively stressed elements across the whole fault surface, calculated prior to each historical earthquake.(b) Percentage of positively stressed elements calculated at 7 km depth.(c) Percentage of positively stressed elements calculated at 15 km depth.Black lines represent the mean percentage over time; the gray area is one standard deviation.While the Cittanova, Sant'Eufemia, Scilla, Serre faults, the range of positively stressed elements is variable, for the MT fault, the percentage of positively stressed elements is constant through time.See FigureS3in Supporting Information S1 for percentage of positively stressed elements across the Sant'Eufemia fault for the scenario in which the Scilla fault ruptures its southern section.

Figure 9 .
Figure 9. Comparing cumulative and coseismic Coulomb stress transfer (CST) for Central Apennines, Southern Apennines, and Calabria.(a) Mean cumulative CST calculated on all faults that did and did not rupture before each earthquake in the three regions.(b) Mean coseismic CST on each fault before each earthquake in the three regions, calculated as a sum of coseismic CST through time.(c) Percentage of faults with positive mean cumulative CST through time.(d) Percentage of faults with positive mean coseismic CST.Data points in gray are for faults that did not rupture; data points in red are for faults that ruptured.The starting date of analysis for each region is indicated with a dashed line.The arrows highlight strong earthquakes occurring on faults with the highest cumulative CST at that time of the event.

Figure 10 .
Figure10.Comparing the percentage of positively stressed elements across the entire surface of faults that ruptured in historical earthquakes in the Central Apennines (a), Southern Apennines (b), and Calabria (c).In the Southern Apennines earthquakes tend to occur on faults that have a high percentage of positively stressed elements; there is more variability in Calabria and in the Central Apennines.