Understanding fault behavior in carbonates is critical because they represent loci of earthquake nucleation. Models of fault-slip mode generally assume: (1) seismic sliding and aseismic sliding occur in different fault patches, (2) creep is restricted to lithology-controlled weak domains, and (3) rate-weakening patches are interseismically locked. We studied three carbonate-hosted seismogenic normal faults in central Italy by combining (micro)structural and geochemical analyses of fault rocks integrated with new seismic coupling estimates. The (upper bound) seismic coupling was estimated to be ~0.75, which indicates that at least 25% of the long-term deformation in the study area is released aseismically in the upper crust. Microscopy and electron-backscatter diffraction analyses revealed that whereas the localized principal slip zone records seismic slip (as ultracataclastic material, calcite crystallographic preferred orientation (CPO), and truncated clasts), the bulk fault rock below behaves differently. Cataclasites in massive limestones deform by cataclastic flow, pressure solution, and crystal plasticity, along with CPO development. Foliated tectonites in micritic limestones deform by pressure solution and frictional sliding, with CPO development. We suggest these mechanisms accommodate on-fault interseismic creep. This is consistent with experimental results reporting velocity-strengthening behavior at low slip rates. We present multiscale evidence of coexisting seismic and aseismic slip along the same fault in limestones during the seismic cycle. Our results imply that on-fault aseismic motion must be added to seismic slip to reconcile the long-term deformation rates and that creep is not exclusive to phyllosilicate-bearing units. Our work constitutes a step forward in understanding fault behavior and the seismic cycle in carbonates, and it may profoundly impact future studies on seismogenic potential and earthquake hazard assessment.

Carbonates commonly host seismic faults worldwide and, notably, in the tectonically active Apennines orogenic system, both in extensional (Cello et al., 1997; Galadini and Galli, 2000; Tondi and Cello, 2003; Mirabella et al., 2008; Chiaraluce et al., 2017) and contractional (e.g., Govoni et al., 2014) kinematics. Stress around upper-crustal faults in carbonates is accommodated by a variety of deformation mechanisms, which span from pressure-dependent brittle-frictional processes to temperature- and rate-dependent deformation, such as diffusive mass transfer and low-temperature crystal plasticity (Tondi et al., 2006; Billi, 2010; Molli et al., 2010; Bauer et al., 2018; Antonellini et al., 2020). Increasing confining pressure (amid other factors such as temperature, strain rate, etc.) favors mesoscopically ductile cataclastic flow and, ultimately, crystal-plastic flow over transgranular cracking and brittle faulting (Fredrich et al., 1989; Baud et al., 2000; Paterson and Wong, 2005).

Faults accommodate displacement through a broad spectrum of slip rates ranging from ≥1 m s−1 during earthquakes to ~10−7–10−3 m s−1 during slow-slip events and tremors up to several millimeters or centimeters per year during aseismic stable sliding (Scholz et al., 1969; Peng and Gomberg, 2010). Spatiotemporal variations in slip behaviors along crustal faults have been mainly attributed to: (1) heterogeneous fault structure and related rheology (e.g., Bedford et al., 2022), (2) variations in pore-fluid pressure (e.g., Dal Zilio and Gerya, 2022), and (3) fault roughness (e.g., Harbord et al., 2017).

Mixed-mode fault-slip behavior along carbonate-hosted faults has been mostly attributed to frictional and structural heterogeneities (Collettini et al., 2011; Gratier et al., 2013a; Bullock et al., 2014; Carpenter et al., 2014; Tesei et al., 2014; Chiaraluce et al., 2017). In this view, seismic sliding and aseismic sliding occur in discrete fault patches with different rheology. Fault patches in more competent units (e.g., limestone) are thought to dominantly behave in a velocity-weakening manner (e.g., Carpenter et al., 2014; Tesei et al., 2014). These patches are locked during the interseismic period and build up high differential stress that is partially released by coseismic slip at failure. Fracturing and cataclasis are primary faulting processes and account for breccia, cataclasite, or gouge formation (e.g., Sibson, 1977). Fault-zone cataclastic products are commonly associated with potential seismic markers, such as mirror-like surfaces, pulverized rocks, truncated clasts, and crystallographic preferred orientation (CPO), among others (Smith et al., 2011; Rowe and Griffith, 2015; Delle Piane et al., 2017). Shear experiments on carbonate gouge have shown that coseismic deformation may also be accommodated by crystal-plastic mechanisms (e.g., Pozzi et al., 2019). Features such as CPO and fault mirror surfaces have also been documented in calcite gouges sheared at subseismic rates (Verberne et al., 2014; Delle Piane et al., 2018). Conversely, fault patches in frictionally weak units (e.g., marly limestone, phyllosilicate-rich fault rock) promote a velocity-strengthening behavior and deform aseismically by frictional-viscous creep involving concurrent fluid-assisted pressure solution, frictional sliding, and foliation development (Rutter and Mainprice, 1979; Bos et al., 2000; Tondi, 2007; Collettini et al., 2011; Gratier et al., 2013a; Bullock et al., 2014). The mixed-mode model envisages that “stick-slip patches” in competent units are interseismically locked, and they should cumulate as much elastic strain as the adjacent weak fault patches slide by aseismic creep.

Despite recent progress in understanding the mechanics of carbonate-hosted faults, there are still some unclear aspects regarding strain partitioning and fault-slip behavior, mainly when applied to the Apennine carbonate multilayer: (1) Comparison of geodetic and seismic moment rates shows that seismic release accounts for only a fraction of the measured geodetic deformation (e.g., Carafa et al., 2020), with a significant fraction of energy released in the interseismic and postseismic periods by aseismic deformation (e.g., Avouac, 2015); (2) slow aseismic rupture and fast earthquakes may spatially coexist (e.g., Amoruso et al., 2002; Smeraglia et al., 2017); (3) experimentally sheared carbonate gouges may creep aseismically (e.g., Delle Piane et al., 2018); and (4) juxtaposition of competent and weak units, which should be characterized by mixed-behavior fault patches (seismic vs. aseismic, respectively), is a long-term result of fault evolution and is not representative of the initial “stratigraphic” relationship (e.g., Fig. 1B). Mixed-mode fault-slip behavior does not explain how the current setting was reached. Changes in clay content and porosity, as well as fabric development, that may affect localization and fault frictional properties (e.g., Bos et al., 2000) are also achieved at a late stage during fault evolution (e.g., Chester et al., 1993), especially if the protolith is a phyllosilicate-poor rock. These considerations indicate that the model whereby (1) stable slip and unstable slip are assumed to be separated in space and to occur in different fault patches and (2) aseismic creep is restricted to lithology-controlled weak domains fails to explain the totality of geological and geodetic evidence.

In this study, we assessed the hypothesis that carbonate-hosted upper-crustal normal faults may undergo interseismic creep. We achieved this by studying three exhumed faults, both active and inactive, in the Northern Apennines (Italy). The general objective of this work was to further investigate the partitioning between seismic and aseismic sliding in the overall slip, as well as the conditions that favor each behavior. We based our work on (micro)structural and geochemical evidence and fault mechanics considerations integrated with new estimates on the seismic coupling using geodetic data and the earthquake catalog. Both analyses showed that part of the long-term deformation is released aseismically in the upper crust. Our work provides new insights into fault-slip behavior in carbonate rocks and may profoundly impact current comprehension of the seismic cycle and fault seismogenic potential, and ultimately earthquake hazard evaluations.

The study area was located within the Umbria-Marche mountain belt of the Northern Apennines (Italy). The latter is an E- to NE-verging fold-and-thrust belt that developed since the late Oligocene as a result of the westward subduction of Adria lithosphere beneath the European plate within the framework of Africa-Europe convergence (e.g., Malinverno and Ryan, 1986). Since the late Miocene, the progressive eastward migration of the thrust front has resulted in the extensional exhumation of the internal region (e.g., Elter et al., 1975; Butler et al., 2004).

The Umbria-Marche ridge consists of, at outcrop, Mesozoic–Paleogene sedimentary units deposited above the basement of the Adria plate (i.e., Permian–Triassic continental siliciclastic material and Paleozoic crystalline basement; e.g., Pierantoni et al., 2013). The deformed passive-margin succession includes Upper Triassic evaporites (anhydrites and dolostones) and a Jurassic–Oligocene multilayer carbonate consisting of platform (Calcare Massiccio Formation) to calcareous-marly pelagic basin units. This succession is stratigraphically overlain by Neogene foredeep and wedge-top basin deposits (Fig. 1B; e.g., Pierantoni et al., 2013). The Umbria-Marche ridge is bounded to the east by the Umbria-Marche-Sabina thrust zone, extending for ~250 km from the Foglia River valley in the northern part, where it doubles the Umbria-Marche sedimentary succession, to Olevano Romano in the southern section, where it superimposes the Umbria–Sabina Mesozoic–Cenozoic succession onto the Latium-Abruzzi carbonate platform (Fig. 1; e.g., Mazzoli et al., 2005). In the framework of NW-SE–trending extension, the axial zone of the chain has been dissected by Pliocene–Quaternary normal, oblique-slip, and strike-slip faults postdating thrust structures and controlling the seismicity and active tectonics of the area (e.g., Boncio et al., 1996; Cello et al., 1997; Galadini and Galli, 2000; Mazzoli et al., 2021). Such fault systems also delimit coeval small-scale depressions, infilled with continental deposits, which commonly coincide with densely inhabited areas (e.g., Amatrice, Gubbio, Norcia; Fig. 1A). Conversely, the external sectors of the orogen and the adjacent foreland (Adriatic offshore) are dominated by active crustal shortening involving both reverse and strike-slip faulting (Basili and Barba, 2007; Govoni et al., 2014; Mazzoli et al., 2014; Chicco et al., 2019). In both these seismotectonic contexts, moderate to large earthquakes nucleate within the Mesozoic–Eocene carbonates and the stratigraphically underlying Upper Triassic evaporites (Mirabella et al., 2008; Govoni et al., 2014; Chiaraluce et al., 2017). The evaporites are interpreted to occur at depths of 5–8 km, where the mainshocks of the normal-faulting seismic sequences that recently struck central Italy are located: Colfiorito 1997, Mw 6.0; Gualdo Tadino 1998, Mw 5.3; and Central Italy 2016, Mw 6.5 (Mirabella et al., 2008; Chiaraluce et al., 2017).

Monte Vettore Fault

The Monte Vettore fault is part of a NNW-SSE–trending, dominantly WSW-dipping, normal fault system that extends for ~30 km along the Sibillini ridge dissecting the western slope of the Neogene NE-verging mountain range front (Figs. 1A and 2A). This fault system consists of several en-echelon segments, each a few kilometers long, which are associated with various shallow synthetic and antithetic splays. The Monte Vettore fault offsets mostly the Mesozoic limestones of the Umbria-Marche succession in the footwall, and it bounds the late Quaternary Castelluccio continental basin to the east (Fig. 2A; e.g., Cello et al., 1997). Displacement on the normal faults in the area is on the order of a few hundred meters, yet it exceeds 1 km along the Monte Vettore fault, i.e., the fault at the base of Monte Vettore, bounding the basin (Fig. 2A; e.g., Pierantoni et al., 2013). The depth of exhumation is ~2–3 km.

The Monte Vettore fault has been active at least since the late Pleistocene (e.g., Galli et al., 2019) and in recent times, during the 2016 central Italy seismic sequence (Mw 6.5), with focal mechanisms indicating NE-SW extension along NW-SE–trending planes, with most of the seismicity compartmentalized between ~3 km and 10 km depth (http://cnt.rm.ingv.it/tdmt). During this sequence, the fault has recorded centimeter-to-meter-scale coseismic surface rupture (Civico et al., 2018), and maximum coseismic slip at hypocentral depths (3–7 km) was larger than 2 m (Cheloni et al., 2017). Paleoseismological studies suggest return times of around 1.8 k.y. for both the Norcia and Monte Vettore fault systems, the latter of which has a slip rate of 1.3 mm yr−1 (Galli et al., 2019).

Gola della Rossa Fault

The Gola della Rossa fault is a NNE-SSW normal fault located in the homonymous gorge, and it crops out in a quarry ~1.5 km to the SW of Serra San Quirico village. It cuts through the carbonate sequence of the Umbria-Marche succession, along the back of the Neogene NE-verging mountain front (Figs. 1A and 2B). The total fault length is unknown, but it has been documented for no less than 3 km (Fig. 2B). This fault is (at least?) post-Aquitanian in age, since it cut through the Scaglia Cinerea Formation (Bartonian–Aquitanian pro parte), and it is buried by Pliocene deposits. However, no definite evidence has been found to determine the tectonic context in which this fault formed. For instance, Chicco et al. (2019, see NFS1 in their work) described an ~50-km-long, NE-dipping, normal-to-oblique-slip fault system that passes through Serra San Quirico–Urbino dissecting parts of the eastern front of the Apennine ridge and nucleates in the Lower Pliocene section pointing to synorogenic stage development. Geological evidence indicates that this fault is inactive.

Gubbio Fault

The Gubbio fault is an ~22-km-long, WSW-dipping normal fault with inferred maximum displacement of 3–3.2 km (Mirabella et al., 2004), and the outcropping extensional fabric was formed at depths <2.5–3 km (Bussolotto et al., 2007). The Gubbio fault cuts through a NE-verging Miocene fold-and-thrust system; it exposes the Jurassic–Oligocene carbonate strata of the Umbria-Marche succession in its footwall and delimits the Quaternary Gubbio continental basin to the east (Figs. 1A and 2C; Menichetti and Minelli, 1991). This fault records multiple deformation stages, including early Miocene pre-orogenic extension and late Miocene synorogenic thrusting; however, most of its displacement is attributed to early Pleistocene extensional slip (Mirabella et al., 2004).

Both instrumental and historical records of seismicity in the area are limited and characterized by weak to moderate earthquakes (http://cnt.rm.ingv.it; Rovida et al., 2022). Focal mechanisms are mainly normal with NW-SE–oriented nodal planes and hypocenters localized along steeply SW- and NE-dipping zones down to ~10 km of depth. The strongest event was the 1984 Gubbio sequence (Mw 5.1), with hypocenters located at 3.5–7 km depth (Collettini et al., 2003). Pucci et al. (2003) argued that the low-angle, deeper portion of the southernmost section of the Gubbio fault was activated. Collettini et al. (2003) suggested, instead, that the sequence involved a region located in the hanging wall of the Gubbio fault. Most of the seismicity related to the 2013–2014 Gubbio seismic sequence (Mw 3.9) instead involved a region in the footwall of the Gubbio fault. In the southern sector, earthquakes occur along and around the low-angle deeper section of the Gubbio fault, and no events have been recorded along its shallower high-angle portion. This inspires questions as to whether it is a locked active fault or an inactive one (Valoroso et al., 2017). These authors suggested a mixed-mode (seismic/aseismic) slip behavior. Active fault databases (DISS Working Group, 2021) list this fault as seismogenic.

(Micro)Structural Analysis

The attitude of structural elements, such as faults and foliation in fault rocks, was recorded in the field. Rock slabs were cut perpendicular to faults and foliation surfaces and parallel to the slip direction from oriented samples collected in the field for microanalytical investigations. Samples were collected within the slip zone (i.e., fault-rock materials embedding the slip surface, which altogether constitute the fault core; e.g., Smith et al., 2011) of each fault, at well-known outcrop exposures. Samples may include the localized principal slip zone and adjacent fault-related rocks, which were described following Sibson (1977). In total, 27 thin sections were prepared by thinning the slabs down, after epoxy impregnation, to ~30 μm in thickness using an Al-based abrasive paste and a diamond-bearing, 0.25-μm-thick powder for polishing. All analyzed samples are listed in Item S1 of the Supplemental Material1.

Optical and scanning electron microscopy (SEM) and cathodoluminescence (CL) were used for the microstructural and microchemical analysis of fault rocks on thin sections, and to investigate the fabric, microdeformation mechanisms, and crosscutting relationships between different microstructures. Optical imaging was performed in transmitted light (TL) by means of a ZEISS Axioscope 5. A ZEISS SIGMA 300 field emission (FE) SEM operating at 15 kV with a working distance of 8.1–8.5 mm, coupled with energy-dispersive X-ray spectroscopy (EDS), was used for microchemical analysis and imaging. The samples were sputtered with carbon before SEM analysis. CL analysis was performed using a Cathodyne® (NewTec Scientific; operated at 8 kV beam energy and 200 μA beam current) equipped with a standard petrographic microscope (Leica DM2700 M). CL imaging has been mainly used to overcome the difficulties in detecting microstructural details via backscattered SEM images, given the chemical homogeneity of the fault-rock material (i.e., calcite). Electron-backscatter diffraction (EBSD) analysis was performed using a Tescan Solaris focused ion beam (FIB) FE-SEM coupled with an Oxford Simmetry EBSD detector, operated at 15 kV accelerating voltage and 3 nA probe current. EBSD data were acquired with a 0.4–1.49 μm step size at a working distance between 7.7 mm and 19.3 mm on samples tilted 70°. The EBSD analysis was conducted both on the bulk gouge (at distances ≫500 μm from the localized principal slip zone; e.g., see Demurtas et al., 2019) and along the principal slip zone. Thin sections for EBSD analysis were Syton-polished and coated with a thin (~5 nm) layer of carbon. Map processing was performed with the MATLAB toolbox MTEX 5.8.1 (Bachmann et al., 2010). Grains were calculated considering a critical misorientation of 10° and a minimum of 7 pixels per grain. For grain-subgrain discrimination, a misorientation angle of 2° was used. Pole figures were obtained from the grains’ orientation distribution function (computed with a de la Vallée Poussin kernel with 10° half width and accounting for one point per grain) and are equal-area, lower-hemisphere projections, oriented according to the fault zone’s kinematic reference frame (where z is the vorticity axis, and x is the slip direction).

X-ray Powder Diffraction

The mineralogical composition of fault-rock samples was characterized by means of X-ray powder diffraction (XRPD) using an automated Philips Bragg–Brentano diffractometer equipped with graphite monochromator. The long fine focus Cu tube was operated at 40 kV and 25 mA. Spectra were recorded in the 2θ range 3°–75° with a 0.02° step and counting time of 2 s. The structural parameters and scale factors of calcite and quartz present in each sample were refined with the program GSAS (Larson and Von Dreele, 2004) to obtain the unit cell parameters of calcite and the quantitative analysis of the phases in each sample. The refinement was carried out in the space group R-3c (calcite) and P3221 (quartz). Cell parameters, scale factors, and the background polynomial functions were free variables during refinement. The intensity cutoff for the calculation of the profile step intensity was initially set at 1.0% of the peak maximum and was lowered to 0.1% of the peak maximum in the final stages of the refinement. Final convergence was assumed to be reached when the parameters shifts were <1% of their estimated standard deviation.

We compared the set of structural parameters obtained using different refinement strategies on the same diffraction data to get an alternative estimate of the accuracy of the refined structural data. These comparisons showed that realistic estimates of the error bars are ±0.0005 Å for the cell parameters of calcite, whereas quartz cell parameters could not be refined due to the very low amount of quartz. The error in the estimation of phase content was ±1% by weight. No attempt was made to determine the amounts and type of clay minerals present by XRPD (e.g., illite, smectite, chlorites) because of the absence of peaks related to clays in the samples examined.

Seismic Coupling in the Central-Northern Apennines

Earthquakes occur when tectonic and resistive forces equalize along the fault plane. Then, the elastic deformation of adjacent rock volumes stored in the interseismic period is suddenly released by rupture propagation and fault slip. Over several earthquake cycles and a sufficiently long time, the elastic strain rate in the lithosphere is marginal to the seismic strain rate. Hence, the long-term tectonic moment rate graphic can be expressed as the sum of the long-term average of seismic moment rates (graphic, resulting from earthquakes) and aseismic moment rates graphic. We can relate these three long-term estimates by:




where c is the long-term seismic coupling. Calculations at the global scale for different tectonic regimes and regions determine stable estimates of c (Bird and Kagan, 2004), whereas quantifying it in small regions is often inconclusive due to the paucity of recorded earthquakes, which does not allow for a stable estimation of the long-term average graphic. On the other hand, using frequency-magnitude distributions allows the contributions of very large (but very rare) earthquakes and small (but very frequent) events to be included simultaneously. This approach dismisses the statistically inconsistent determination of graphic as the summation of moment rates of earthquakes that occurred in a small region. Among different frequency-magnitude distributions, the tapered Gutenberg-Richter (TGR) is preferred because it is computationally simple and manageable and is based on two main parameters: the asymptotic spectral slope at small moments (β), and the corner moment (Mc). For a statistically sound earthquake sample (see Kagan, 2014), the TGR describes the fraction (G) of earthquakes with a moment exceeding M as:


where Mt is the threshold moment for completeness of the catalog. Considering the total number of earthquakes α0 with M ≥ Mt in a region, the seismic moment rate graphic can be calculated as:


where Γ is the gamma function. For the Apennines, Carafa et al. (2017) determined graphic and graphic. The Mc upper bound (defined as the 95% confidence range) cannot be determined in the likelihood map surface, but the maximum length of active faults in the area (~45–50 km) suggests that mc = 7.17 is the maximum admissible corner magnitude. If we consider the CPTI5 v4 (Rovida et al., 2022), a threshold magnitude mt of 5.570 ± 0.115 can be set for the study area according to Meletti et al. (2021). Following a statistical approach, the catalog is complete at the threshold magnitude mt from A.D. 1750, whereas it goes back to A.D. 1650 using a historical approach (Meletti et al., 2021). It follows that α0 can vary in the range of 18–33 earthquakes (in 270 yr) following the statistical approach and 22–38 (in 370 yr) for the historical approach.

For graphic, it has recently been proven that geodetic horizontal velocities adequately describe long-term kinematics at both global (Bird et al., 2008; Bird and Kreemer, 2015) and regional scales (Barka and Reilinger, 1997; Reilinger et al., 2006; Zhu and Shi, 2011). Indeed, geodetic strain rates may be elastic, permanent, or mixed because they are based on time series for which the character is uncertain. The elastic rebound theory, however, indicates that the spatial separation between sites experiencing elastic strain accumulation and nonelastic permanent strain may be as little as a few kilometers. Also, correcting observed strain rate fields to an idealized long-term equivalent requires a complete fault database, which is often not available (see Carafa et al., 2020). Therefore, we can assume that modeled geodetic strain rates approximate long-term permanent strain rates (Bird and Kreemer, 2015; Carafa et al., 2017). This assumption smooths out the graphic local maxima but not its regional value, and the principal components of the geodetic strain rate tensor can be used to determine the long-term seismic moment rate graphic of a crustal volume. In detail, the diagonalization of the 2 × 2 geodetic horizontal-plane strain rate tensor gives the principal horizontal values graphic (more negative or compressional) and graphic (more positive or extensional) and their associated principal axis directions. The vertical (principal) strain rate tensor graphic is inferred under the incompressibility assumption graphic. The three principal values of the geodetic strain rate tensor are now defined and can be identified as graphic, respectively. In such kinematic settings, the active faults have their slip-rate vectors in the plane containing the principal axes graphic and graphic. In addition, a less active conjugate fault set is likely to exist with slip-rate vectors in the plane containing both the principal axes graphic and graphic (if graphic), or graphic (if graphic). Now, if we define the principal strain rates as:





we can determine the long-term tectonic moment rate graphic of a crustal volume with ground surface (A) and depth of the seismicity cutoff (z) as:


where θ1 and θ2 are the average angles dividing the fault planes from graphic in the graphic plane and in the graphic plane, respectively (Carafa et al., 2017). The shear modulus (μ) can be set to globally averaged values, or it can be regionally determined if S (or P) velocities and densities are available. We opted for this second alternative using the tomographic model of Di Stefano and Ciaccio (2014) and setting the averaged upper-crustal densities to ρ = 2670 ± 70 kg m−3 (Improta et al., 2003; Tiberti et al., 2005). The resulting values of μmin 3.11 × 1010 Pa and μmax = 3.90 × 1010 Pa represent the shear modulus bounds propagating the error in densities and considering the lateral variability of P velocities.

The long-term seismic and aseismic moment rates can be written as:




where we refer to the product c · z as the “coupled thickness of seismogenic lithosphere” (Bird et al., 2002; Bird and Kagan, 2004), while graphic is the diffusivity (or potency rate per unit depth) of the deforming volume.

Now, we can take advantage of seismicity and geodesy to calculate the coupled thickness (and infer the seismic coupling) as:


In this work, we relied on the recently published velocity solution of Serpelloni et al. (2022) to determine graphic and investigate the 〈cz〉 variability for our study area by considering the uncertainties in α0, μ, and Mt and lateral smoothing of the geodetic strain rate.

(Micro)Structural Fault Characterization

Mount Vettore Fault

For this study, we selected a prominent synthetic splay of the Monte Vettore fault system. The sampled outcrop (42°50′9.92″N, 13°13′46.37″E) was located in the fault footwall in the vicinity of Colli Alti e Bassi (Fig. 2A). Here, the fault surface is oriented, on average, N150°/65° (strike azimuth/dip) with striae pitch >70°, indicating normal slip with a left-lateral component. The fault juxtaposes the massive limestones of the Calcare Massiccio Formation in the footwall with Maiolica Formation calcilutites in the hanging wall (Figs. 2A and 3A). This splay has a total length of ~2.5 km. The fault recorded coseismic surface faulting during the 2016 seismic sequence, and surface coseismic throw is on the order of 0.2–1.0 m (Fig. 3A). The freshly exposed fault surface (i.e., slip surface) is characterized by smooth, whitish cataclastic slickenlines, grooves, and metric-scale undulations with the long axes oriented parallel to the slip direction (Fig. 3A).

Microstructural observations of the slip zone (i.e., fault-rock material) refer to samples collected from the footwall block (Fig. 3). The slip zone and slip surface constitute the fault core. The slip surface displays a sharp boundary that locally truncates calcite clasts, as observed both directly on the fault surface in the field and in thin section (Fig. 3C). Moving away from the slip surface, the slip zone is characterized by an irregular layer of ultracataclasite a few millimeters thick (Fig. 3C), which was observable also at the sample scale (Fig. 3B). Larger clasts dispersed within this layer are submillimeter in scale, whereas the fine-grained matrix is made up of particles in the order of tens of microns in diameter or smaller (Fig. 3C). Below this layer, there is a well-developed cataclasite up to ~80 mm thick (Figs. 3B and 3C). Within the bulk cataclasite layer, there are shear bands of more comminuted material with respect to the surroundings, and these are subparallel or oriented at an acute angle (dipping toward the SW) to the slip surface (Figs. 3B and 3D). The arrangement of fragments around and within these zones suggests their movement was synthetic to that of the principal slip surface. The cataclasite consists of angular to subrounded particles, generally below 1 mm in size, surrounded by a calcite-rich, fine-grained matrix. Clasts are mainly calcite derived from the parent limestone, either pristine, preserving the typical grainstone texture, or partially recrystallized, or from secondary calcite cement (e.g., reworked veins or pore-filling cement) and cataclasite fragments. In places, the cataclasite is matrix supported. In other places, the clasts are in contact with each other, locally forming trains of aligned grains (Fig. 3E), and they show predominantly dissolution and microcracks at contact points (Figs. 3C and 3E). Pressure solution is mostly intergranular, and it is manifested by the presence of packed grains with sutured contacts (Figs. 3C and 3F). Poorly developed solution seams truncate the cataclastic material and impart a weak, nonrandom fabric to the fault rock (Fig. 3C). Shear fractures at high angle to the slip surface cut through the solution seams. Solution surfaces, “chains” of aligned grains in contact, and shear fractures show a kinematic consistency with slip along the slip zone (i.e., solution seams are normal to the component of maximum stress, whereas shear fractures and trains of aligned grains are parallel to the component of maximum stress) (Ramsay, 1967). Some calcite clasts and crystals show twinning and undulose extinction (Figs. 3D, 3E, and 3G). Finally, there is a zone of fault breccia, which extends in the footwall block, beyond the size of the analyzed samples (Fig. 3B). This sector is primarily characterized by fracturing, calcite veins, and pressure solution (see Supplemental Item S2). Its boundary with the cataclasite is not abrupt, with some fine-grained cataclastic material filling fractures and spaces among coarse-grained clasts within the fault breccia (Fig. 3B).

The EBSD analysis was performed at two sites: (1) site v1, within the ultracataclasite layer immediately beneath the principal slip surface (i.e., localized principal slip zone; Figs. 3B and 4A), and (2) site v2, within the underlying bulk cataclasite layer (Figs. 3B, 4B, and 4C). Pole figures show that calcite develops, in both sites, a weak crystallographic preferred orientation (CPO; Figs. 4D, 4E, and 4F). At site v1, the CPO (multiple of uniform distributions (m.u.d.) up to 1.6) is defined by a clustering of the calcite c axes (0001) parallel or at a small angle (clockwise) to the y direction (i.e., normal to the principal slip surface, which is parallel to the x-z plane; Fig. 4D). At site v2, the CPO (m.u.d. up to 1.5) is produced by clustering of c axes at small angles to the x direction (i.e., parallel to the shear direction of the principal slip surface), and it is inclined in the direction of shear (Figs. 4E and 4F). The a-axis <1120> directions are more dispersed, especially at site v1 (Fig. 4D), and arranged perpendicular to the c axes at site v2 (Figs. 4E and 4F). At both sites, misorientation angle distribution plots show a maximum between 75° and 80° (Figs. 4D, 4E, and 4F) related to e-twin (e{01–18}) development in calcite, with the highest frequency found at site v2 (Fig. 4F). Correlated distributions are close to uniform for the three sites, with only site v2 in Figure 4E showing an increase in frequencies at low misorientation angles (<10°), indicating the presence of low-angle grain boundaries (Fig. 4E; see Supplemental Item S3).

Gola della Rossa Fault

At the sampled outcrop (43°26′3.33″N, 13°0′31.31″E), the fault surface is oriented N27°/66° with striae pitch around 55°–60°, indicating oblique slip (dominant normal slip with a right-lateral component). Here, the fault juxtaposes the massive limestones of the Calcare Massiccio Formation in the footwall with the limestones and marly limestones of the Scaglia series (Scaglia Variegata and Scaglia Rossa Formations) in the hanging wall (Figs. 2B and 5A). The pristine rock fabric and sedimentary bedding attitude are preserved in the undeformed footwall (Fig. 5A). The stratigraphy indicates that this part of the sedimentary succession is condensed (Guerrera et al., 2014). A dip separation of 1600–1700 m was obtained by considering the stratigraphic juxtaposition between the footwall and the hanging wall, resulting in a net slip of 1900–2000 m considering the pitch of striae (Fig. 5A). If we consider that the fault is buried under Pliocene deposits, and by excluding significant erosion of the pre-orogenic succession above the exposed part of the fault prior to faulting (up to the Messinian deposits), we can infer that the fault has been exhumed from ~2300 m depth.

Microstructural observations refer to fault-rock samples collected from the footwall block (Fig. 5). The exposed fault plane is characterized by a polished slip surface showing mechanical striae and truncated clasts (Figs. 5A and 5B). At the sample scale, it is possible to observe an ~20–30-mm-thick, dark-gray layer that includes the principal slip surface. An irregular convoluted boundary running subparallel to the slip surface separates this layer from a lighter-gray one below (Fig. 5B). However, no apparent microstructural or geochemical difference was observed between the two layers (see Supplemental Item S4). The cataclasite layer comprises the whole thickness of the studied sample; i.e., it is at least ~90 mm thick (Fig. 5A). It is made of sub-centimetric comminuted grains with angular to smoothed boundaries, dispersed in a fine-grained matrix composed of sub-millimetric particles. Some larger clasts are truncated by the slip surface (Figs. 5C and 5D). Clasts and matrix show similar features to the fault rocks of the Monte Vettore fault (Figs. 5C and 5D). Several micrometer-scale (transgranular) fractures filled with calcite are subparallel to low angle with respect to the slip surface, and they occur as single fractures or in swarms, with no or limited shear. An ~2-µm-thick, fault-parallel vein is located beneath the slip surface (Fig. 5C). There is no clear evidence to establish the relative timing of these fractures (veins). Intragranular fractures occasionally exploit cleavage planes as preferential breakage surfaces, and, accordingly, the outline of some calcite clasts follows the cleavage orientation (Fig. 5C). When clasts interact, they commonly display sutured contacts and local grain indentation reflecting pressure solution. Grain in contacts sometimes form chains with preferential disposition, and intragranular microcracks are disposed with the long axes oriented parallel to these trains of grains (Fig. 5C). Pressure solution is manifest also through seams that markedly truncate the cataclastic material, and they occur either at high or at small angle (20°–30° toward the WNW) with respect to the slip surface (Fig. 5E; see Supplemental Item S4). The latter seams are located near the slip surface and are characterized by flat-to S-shaped irregular surfaces. The solution seams are consistent with the shear couple acting on the fault plane (Figs. 5B and 5E). Deformation twinning is widespread (Figs. 5C and 5E), and sometimes twins are dissolved at grain boundaries, or they end before the boundary is reached (Fig. 5F). Other clasts, instead, display undulose extinction and lobate grain boundaries (Fig. 5F).

Gubbio Fault

The two sampling sites used for this study were both situated within the footwall block. The Cava Filippi outcrop (43°23′22.55″N, 12°29′51.53″E) is located near the NW fault termination (Fig. 2C). Here, the fault surface is oriented, on average, N141/63°, and it cuts through layered micritic limestones and marly limestones of the Scaglia Rossa Formation (Fig. 6A). Protolith fabric and sedimentary bedding attitude are preserved in the undeformed footwall. Moving basinward, the fault zone exposes distinct structural domains separated by major SW-dipping slip planes (Fig. 6A). Farther southwestward, a major morphologic fault escarpment juxtaposes the footwall domain with the hanging-wall basin-fill deposits (Figs. 2C and 6A). Fault-zone domains in between localized slip surfaces are characterized by both foliated (S-C fabric; sensu Berthé et al., 1979) and random fabric fault rocks. Penetrative solution cleavage (S) and secondary shear planes (C) define sigmoidal lithons and crosscut sets of mostly subvertical, calcite-filled veins (Figs. 6A, 6B, and 6C; see also sketch in Fig. 6D). At the mesoscale, foliation intensity increases toward the marly limestone protolith; still, fault rocks in the micritic limestones show foliated domains at multiple scales, e.g., at outcrop scale (Fig. 6A), at sample scale (Fig. 6D), and at the microscale (Fig. 6E). Slip surfaces are decorated with slickenlines and slickenfibers (Figs. 6B, 6C, and 6D). At times, these accessory structures occur also along C-plane surfaces. All the structural elements trend roughly NW-SE. Principal slip planes are arranged at moderate to high dip angle (50°–70°), and C planes lie approximately parallel to slip planes, whereas S planes gently dip (10°–30°) to the SW (Fig. 6A).

The Cemetery outcrop (43°20′41.96″N, 12°35′47.54″E) is located along the central part of the fault (Fig. 2C). Here, the fault is oriented, on average, N127/67°, and it juxtaposes marls and marly limestones of the Marne a Fucoidi Formation against the micritic limestones of the Scaglia Bianca Formation (Fig. 6F). The general structure of the fault zone is similar to that exposed in the northern outcrop, and major SW-dipping fault planes separate domains of foliated gouges and cataclasites. Kinematic indicators imply normal-slip fault kinematics (pitch >85°; Figs. 6A and 6F). The fault-zone structure at the two sites has been described in more detail by Bussolotto et al. (2007) and Bullock et al. (2014).

For the purposes of this study, our observations and sampling in both sites were focused on the micritic limestones of the Scaglia Rossa and Scaglia Bianca Formations, in the footwall block. Fault rocks in the Scaglia Rossa Formation (Fig. 6; see Supplemental Item S5) are characterized by millimeter- to centimeter-thick, foliated, and less frequently massive, calcite cemented-breccia and cataclasites. Thin (<2 mm) ultracataclasite layers are observed just below the slip surface (Fig. 6E). In Figure 6E, the ultracataclasite is a random fabric, and it is made of calcite clasts in the order of 10 µm and below, dispersed in an ultrafine-sized matrix indistinguishable with a light microscope (hence <1 μm in size). In other cases, the ultracataclasite layer is constituted by calcite particles dispersed in brownish material that has an appearance of a foliation, and it is separated from the calcite below by multiple sharp surfaces, parallel to the slip surface, and characterized by fine-sized compacted calcite clasts (see Supplemental Item S5). In some cases, a calcite vein is found beneath the slip surface. Microcrystalline calcite slickenfibers decorate slip surfaces (Fig. 6D). Foliation is quite ubiquitous but is more pervasive in some domains than in others, and it is defined by thin and irregular surfaces filled with a brownish fine-grained material and where calcite clasts have been truncated (Fig. 6E). Foliation surfaces show evidence of dissolution, shear, and cataclasis (Fig. 6E). Foliations separate domains of cement precipitation (veins) and domains made of brecciated/cataclastic material, which is calcite detrital grains of the parent rock (brownish micritic limestone) and crystalline cement (e.g., reworked veins or intragranular cement; Fig. 6E). Sometimes, detrital grains record multiple phases of fracturing (and veining) and pressure solution, and such features are confined within the grain by the foliation; the most prominent set of veins is approximately normal to the solution cleavage plane (S). In fact, veins filled by coarse calcite that cut through parent-rock detrital grains are bounded by pressure solution surfaces (see Supplemental Item S5). The latter have mutual crosscutting relationships with the intragranular fractures that affect the brecciated material (i.e., some fractures are confined within the detrital carbonate grains, and they are bounded by solution surfaces, whereas others cut through the foliated fabric). Clasts show well-developed twins (Fig. 6E) and undulose extinction and lobate contacts (see Supplemental Item S5).

Fault rocks in the Scaglia Bianca Formation (Fig. 6; see Supplemental Item S5) are characterized by pervasive and anastomosed surfaces with S-C-C′ fabric (Fig. 6G), although C′ planes (secondary shear surfaces dipping toward SW) are not entirely developed. Sigmoids are bounded by pressure solution surfaces that define the foliation, and a thin film of brownish fine-grained material is present within these surfaces (Fig. 6H). Sigmoids are mostly made of coarse-sized crystalline calcite, most likely vein-derived, and brecciated detrital grains of the parent rock (brownish micritic limestone with fossils), without significant cataclasis (Fig. 6G). These grains have been affected by multiple phases of fracturing and pressure solution, and these features are confined within the grain by irregular and indented solution cleavage planes manifesting both dissolution and shear (Fig. 6H). Grains that are near or in between solution seams have undergone partial recrystallization and have generally a fine size (see Supplemental Item S5). Fracturing exploits the cleavage planes, and this is more intense near the foliation surfaces (Figs. 6I and 6J). Calcite clasts are typically intensely twinned; however, twins do not always reach the grain boundary. They also show undulose extinction and lobate contacts (Figs. 6I and 6J). At both outcrops, the fault fabric is compatible with the stress couple acting parallel to the slip of the shear zone (see sketches in Figs. 6D and 6G).

The EBSD analysis was performed on a fault-rock sample of the Scaglia Rossa Formation (Fig. 6D) from three sites: (1) site g1, in a slickenfiber filled with blocky calcite crystals along a slip surface (Figs. 7A, 7B, and 7D); (2) site g2, within a region made up of fine-grained and elongated calcite crystals with respect to the surroundings, and which lies immediately beneath the same slip surface (Figs. 7A, 7B, and 7E); and (3) site g3, within a rhomboidal region made up of both elongated and equant calcite crystals, with finer grain size than the surroundings, along an incipient foliation surface with shear (Figs. 7A, 7C, and 7F). Pole figures of site g1 (Fig. 7D) do not show a clear CPO within the mineral aggregate. In site g2, the clustering of calcite c-axes (0001) parallel or at small angle (counterclockwise) to the x direction (i.e., parallel to the principal slip surface, which corresponds to the x-z plane) produces a strong CPO (m.u.d. up to 6.9; Fig. 7E). The orientation map shows that elongated calcite crystals are also oriented with their long axis (sub)parallel to the x direction (Fig. 7E). The a-axis <1120> directions define a girdle perpendicular to the c axes with maxima close to the y direction (Fig. 7E). In site g3, a relatively weaker CPO (m.u.d. up to 2.6) is produced by the clustering of the calcite c axes at small angle (clockwise) or parallel to the x direction (Fig. 7F). The orientation map shows that elongated calcite crystals are also oriented with their long axis ~20° clockwise to the x direction (Fig. 7F). The a-axis directions define a wide girdle normal to the c axes with maxima close to the periphery of the pole figure (Fig. 7F). At each site, misorientation angle distribution graphs show a strong peak at low angles (<10°; Figs. 7E and 7F; see Supplemental Item S2), indicating the presence of low-angle grain boundaries, with the lowest frequency found at site g1 (Fig. 7D; see Supplemental Item S2). The maxima at low misorientation angles at site g2 and site g3 gradually decrease toward higher angles, and the correlated distribution clearly deviates from the random distribution up to ~20°. This trend is apparent at site g3 (Fig. 7F), and it is consistent with the presence of the CPO at these sites (Figs. 7E and 7F). Except for site g2, misorientation angle distributions show a peak at 75°–80° (Figs. 7D and 7F) due to e-twinning of calcite, with the highest frequency found at site g1 (Fig. 7D).

Geochemical Analysis

XRPD and SEM microchemical analyses showed that fault rocks (cataclasite) developing in the massive limestones (Calcare Massiccio Formation) are constituted essentially by calcite (Figs. 8A and 8B; Table 1; see Supplemental Item S4), which constitutes fragments, fine-size matrix, and fracture filling in veins (Figs. 3 and 5). The amount of clay minerals was not detectable by XRPD data.

XRPD patterns of fault rocks (foliated cataclasite, foliated tectonite) developing in the micritic limestones (Scaglia Rossa and Scaglia Bianca Formation) showed two peaks of quartz; however, calcite still represented the dominant mineral (Figs. 8C, 8D, and 8E; Table 1; see Supplemental Item S5), since it constitutes fragments, matrix, and fracture filling in veins (Fig. 6). Quartz was present both in grains with sizes up to a few hundreds of micrometers and grains finely dispersed in the sample (see Supplemental Item S5). Smectite (here, being dominantly constituted by montmorillonite) and illite dominated over other clay minerals (e.g., kaolinite, chlorites; e.g., Bullock et al., 2014), the amount of which was not detectable by XRPD data (Figs. 8C and 8E). SEM analysis showed that Si and Al content dominated over Ca along stylolitic surfaces (Fig. 8D) and revealed the presence of small amounts of goethite and fluorapatite (see Supplemental Item S5).

Seismic Coupling

In Figure 9, we show the 〈cz〉 variability due to the uncertainties in α0, μ, and Mt. and the lateral smoothing of the geodetic strain rate. To determine this variability, we first defined the long-term seismic moment rate graphic by applying Equation 4 considering the α0 range for both historical and statistical completeness periods. The resulting graphic range is 3.7–4.7 × 1017N m yr−1.

Then, we determined the tectonic diffusivity (and related uncertainty) in the study area by using the global navigation satellite system (GNSS) horizontal velocity from Serpelloni et al. (2022) and the stress orientations from the Italian Present-day Stress Indicators (IPSI) database (Mariucci and Montone, 2020) with the NeoKinema code (Bird and Carafa, 2016; Carafa and Bird, 2016). In Figure 10, we show the fit between the data and the tectonic model. The resulting tectonic diffusivity range is 1.15–1.44 × 103 m2 yr−1. For the variability of the shear modulus (μ), as in the Methods subsection, we set the value in the 3.11–3.90 × 1010 Pa range. The final range of admissible 〈cz〉 was graphic km.

The definition of seismic coupling (c) for the study region from the resulting graphic km requires a quantification of the seismicity cutoff (z). We extracted the study area seismicity from the Italian Seismological Instrumental and Parametric Data-Base (ISIDE; available at http://iside.rm.ingv.it/, last accessed 31 July 2022) and found that 90% of the earthquakes occurred in the first 12–13 km of the upper crust. This range agrees with the seismogenic thickness determined by Chiarabba and De Gori (2016). If we set z = 12 km, then graphic, where the indices define the upper and lower ranges. If we select z = 13 km, then graphic. Because we used only the maximum admissible corner magnitude (mc = 7.17) in our graphic quantification, hence not varying it to lower values, the obtained seismic coupling ranges can be equally seen as upper ranges. Thus, we can reasonably assume that (at least?) 25% of the long-term deformation in the study area is released aseismically in the upper crust.

In this work, we studied fault rocks belonging to carbonate-hosted seismogenic normal faults in central Italy through detailed microstructural and geochemical analyses integrated with new estimates on the seismic coupling. Hereafter, our results will be discussed in terms of: (1) the frictional behavior of faults in carbonate rocks upon reactivation, (2) the deformation mechanisms that may accommodate interseismic fault creep, (3) the partitioning between seismic and aseismic slip, and (4) the implications for the seismic cycle.

Frictional Behavior of Carbonate-Hosted Faults

The studied faults represent exposed analogues of carbonate-hosted normal faults and fault rocks that are representative of those formed at seismogenic depths. At low (subseismic) slip rates and low temperature, both intact calcite samples and gouges show coefficients of friction (μ) in the range 0.6–0.8 (Verberne et al., 2014; Tesei et al., 2014; Delle Piane et al., 2016). This coefficient drops to sub-Byerlee values when samples are sheared at seismic velocities (e.g., Han et al., 2007; Tisato et al., 2012). As gouge thickens, particle angularity decreases, and surfaces smooth in equilibrium with the slip mode, the peak friction can evolve to a lower steady-state value, and slip may stabilize (Marone et al., 1990; Anthony and Marone, 2005). Stick-slip behavior is generally promoted by slip localization and increasing sliding velocity, bulk temperature, and confining pressure (e.g., Verberne et al., 2014; Passelègue et al., 2019), provided the rock remains brittle. Experimental tests, however, show that faults in carbonates do not exhibit exclusively a velocity-weakening behavior (i.e., decrease in strength upon increasing slip rate), which is a prerequisite for earthquake rupture nucleation (e.g., Scholz, 2019). In fact, carbonate gouges sheared at a range of pressure (50–100 MPa normal stress), temperature (~25 °C to >100 °C), and saturation (water saturated, dry) conditions and at low slip rates (0.1–100 µm s−1; i.e., subseismic) can show velocity-strengthening behavior (i.e., increase in strength upon increasing slip rate; Tesei et al., 2014; Verberne et al., 2014; Mercuri et al., 2018). This may result in stable sliding and aseismic creep. Similarly, faulting and subsequent reactivation tests conducted at room temperature on saturated, calcite-rich (≥98%) samples under subseismic slip rates indicated stable slip without slip weakening (Delle Piane et al., 2016, 2018). Based on these data and assuming a typical geothermal gradient of ~20–25 °C km−1 and long-term strain rates in the order of a few millimeters per year (e.g., Carafa et al., 2020), carbonate-hosted faults can be expected to slide aseismically until depths of ~3–6 km. Velocity-strengthening friction at shallow depths (<3–5 km) is consistent with other laboratory tests and seismological observations (Marone et al., 1990; Marone and Saffer, 2015). These depths comprise the upper limit of the Apennines seismogenic zone, which, in some cases, may be as shallow as ~2 km (http://cnt.rm.ingv.it/tdmt). To a first approximation, this evidence may help to explain the (in some way enigmatic) coexistence of slow aseismic phenomena and destructive seismic rupture in central Italy (e.g., Amoruso et al., 2002). In addition, carbonate fault mirrors, owing to extremely low frictional healing rates, are likely to be frictionally stable (rate-and-state friction parameter (a-b) > 0) and to creep aseismically (Park et al., 2021).

This evidence suggests that (1) slip behavior may not be simply governed by lithology or structural heterogeneities and (2) faults in competent carbonate units may undergo aseismic slip. Hence, new physical mechanisms need to be invoked that may accommodate interseismic creep.

Deformation Mechanisms for Interseismic Creep in Carbonate-Hosted Faults

Microstructural observations revealed that cohesive fault rocks in the massive limestone dominantly underwent pervasive intergranular microfracturing, chipping, and resultant extensive grain-size reduction and particle rearrangement (Figs. 3 and 5). These microscale mechanisms define a process commonly referred to as cataclasis (Sibson, 1977; Billi, 2010). Localized cataclasis is a dominant process during brittle faulting (e.g., Sibson, 1977). Cataclastic flow, however, is a mesoscopically ductile, microscopically brittle, granular deformation mechanism involving stable and homogeneously distributed microcracking, frictional sliding, and rolling among fragments, and grain plasticity (Rutter, 1986; Paterson and Wong, 2005). This mechanism is strongly pressure-sensitive, and it might occur at low differential stresses, comparable with upper-crustal conditions (Edmond and Paterson, 1972; Renner and Rummel, 1996). Cataclastic flow appears to be independent of rate, so it may produce similar fault-rock textures during seismic slip as during creep, making it difficult to discriminate the regime from which the cataclastic products resulted (Sibson, 1977; Stünitz et al., 2010).

In the studied faults, the presence of truncated clasts along the slip surface (Figs. 3F and 5D) is indicative of seismic slip (Fondriest et al., 2013). A thin ultracataclasite layer lying immediately beneath the principal slip surface (Figs. 3B and 3C) may also suggest fast slip, defining a localized principal slip zone where most of the coseismic slip is accommodated (e.g., Smith et al., 2011). Other potential markers of seismic slip rates include polished slip planes with frictional wear striae (Figs. 3A and 5A; e.g., Fondriest et al., 2013). However, calcite-rich fault patches composed of tightly packed nanograins similar to mirror-like surfaces have been observed in experiments at stable slip and slow rupture conditions (Verberne et al., 2014; Passelègue et al., 2019). The CPO documented within site v1 (Figs. 4A and 4D) can be compared to that observed in localized slip zones in carbonate gouges sheared at seismic slip rates (i.e., (0001) planes subparallel to the shear plane; <1120> axes subparallel to the shear direction), and this comparison suggests that the CPO may result from grain-size–sensitive and grain-size–insensitive plastic creep mechanisms promoted by frictional heating during coseismic slip (e.g., Pozzi et al., 2019), except that, in our case, low-angle grain boundaries are really close to the random distribution (Fig. 4D). Conversely, deformation in the bulk cataclasite layer below is dominated by distributed cataclasis associated with pressure solution, which is mostly intergranular, and crystal plasticity, i.e., twinning and undulose extinction (Figs. 3 and 5), indicative of dislocation glide. Lobate contacts and sutured boundaries (Figs. 3G and 5F) may suggest incipient grain boundary migration. We suggest that along the studied faults in massive limestone units, once the fault gouge is formed during faulting and the accompanying dynamic dilatancy, interseismic creep can progress via cataclastic flow aided by pressure solution and intracrystalline plasticity. The CPO pattern observed in site v2 (i.e., (0001) axes inclined in the direction of shear to subparallel to the shear plane; <1120> axes approximately perpendicular to the c axes; Figs. 4E and 4F) is different from that observed in the highly deformed, ultrafine-grained region underneath the principal slip surface (Fig. 4D), and in most CPOs documented in slip zones in carbonate rocks (e.g., Delle Piane et al., 2017, and references therein; Pozzi et al., 2019). This may suggest they were generated by different mechanisms. The CPO was also accompanied by development of subgrains in the bigger clasts (as evidenced by low-angle misorientations in Fig. 4E; see Supplemental Item S3) and e-twin production. The development of calcite CPO in the bulk cataclasite layer, which we suggest has been formed during slow interseismic creep, may be the result of different concurring deformation mechanisms: intergranular pressure solution coupled with simultaneous e-twin generation and/or subgrain rotation on the favorably oriented clasts. For instance, Smith et al. (2011) proposed that intergranular pressure solution may steer calcite CPO development, producing a pattern similar to that found in the bulk cataclasite layer of the Monte Vettore fault. Such a CPO pattern can be attributed to a crystallographically dependent anisotropy in dissolution and precipitation rates of grains during deformation; i.e., they would dissolve and grow faster in the direction parallel to the c axis (e.g., Power and Tullis, 1989; Bons and den Brok, 2000). According to this process, grains will grow quicker than others if their c axis is parallel to the direction of minimum principal stress, while they will dissolve faster if their c axis is parallel to the direction of maximum principal stress. Calcite c axes are broadly parallel to the orientation of solution seams and sutured contacts among grains (cf. Figs. 4E and 3C), which form preferentially at high angles to the σ1 direction. Then, the observed calcite CPO pattern (Fig. 4E) may result from preferential crystal growth parallel to σ3. Our interpretation is in accordance with microstructural observations of experimental carbonate gouges deformed at subseismic slip rates and low temperature, which propose the concurrent operation of granular/cataclastic flow and grain-scale creep mechanisms (i.e., intergranular pressure solution and/or crystal plasticity) and CPO development (Verberne et al., 2014; Delle Piane et al., 2018; Demurtas et al., 2019). Intragranular microcracking and cataclastic flow are aided by intergrain weaknesses produced by mechanical twinning (Fig. 5C), and dislocation slip and increasing strain during cataclastic flow may increase twinning density (e.g., Fredrich et al., 1989; Schubnel et al., 2005). Also, brittle creep, if operative, would reduce the stress conditions necessary for frictional wear processes accompanying stable sliding in fluid-present fault zones. For instance, Oncken et al. (2022) suggested that cataclastic flow aided by brittle creep govern slow-slip transients in the seismogenic zone. Hence, our observations on natural fault zones are in general agreement with experimental observations. Simulated (Carrara marble) fault gouges deformed at subseismic velocities (0.3–10 µm s−1) and 100 MPa normal stress display distributed deformation with evidence of cataclasis, pressure solution, and intracrystalline plasticity, and velocity-strengthening behavior (Mercuri et al., 2018). Although the authors of these studies proposed that the activation of plastic deformation was responsible for strain weakening at the slow sliding slip rate (0.3 µm s−1), which in turn would promote unstable fault slip, the weakening they observed could actually facilitate stable aseismic slip (e.g., see Hadizadeh, 1994).

In the micritic limestone, the occurrence of sharp slip surfaces decorated with slickenlines (Figs. 6A, 6B, and 6C), cataclastic grain-size reduction, and thin massive ultracataclasite layers (Fig. 6E) within the fault zone may testify to unstable frictional behavior. Slickenfibers coating slip surfaces (Fig. 6D) suggest fluid-assisted nonfrictional slow aseismic creep processes involving pressure solution of asperities, diffusive mass transfer, and mineral fibers crystallization (e.g., Gratier et al., 2013b). Healed cracks with euhedral crystals in fibers may indicate seismic events interrupting their growth. Bullock et al. (2014) also documented recrystallized calcite textures and clay transformation that have been associated with coseismic thermal decomposition. Deformation in the bulk volume (tectonite), however, is dominated by pressure solution, foliation development, and frictional sliding on discrete surfaces (Fig. 6). The role of brecciation/cataclasis is limited when compared to fault rocks that developed in massive limestones. Deformation twinning, cuspate grain boundaries, and undulose extinction (Figs. 6I and 6J) are microstructural evidence for crystal plasticity. We propose that permanent creep processes, including concurrent pressure solution and slow slip along discrete foliation surfaces, accommodate interseismic stable sliding in micritic limestone. Such coupled frictional-viscous processes, in fact, have been suggested to allow fault creep at low-stress conditions (e.g., Rutter and Mainprice, 1979; Bos et al., 2000), and cataclastic grain-size reduction may promote diffusive mass transfer (e.g., Gratier et al., 2013b). Similar observations have been made by Bullock et al. (2014) for the Gubbio fault. However, they suggested that post- or interseismic creep was dominantly accommodated in the marly limestone domain, while faults in micritic limestone remained locked and accommodated coseismic slip (see Fig. 6A). In these fault rocks, the development of calcite CPO, which we suggest is formed during slow interseismic creep, is likely related to oriented crystal growth during diffusion-precipitation creep (e.g., Giuntoli et al., 2018), as suggested by the consistency among crystal elongation directions and the CPO orientations (Figs. 7E and 7F). Otherwise, it could be attributed to crystallographically controlled mineral dissolution and growth combined with passive grain rotation (e.g., Bons and den Brok, 2000). The strong maximum at low values in the misorientation angle distribution graph (Figs. 7E and 7F; see Supplemental Item S3) suggests that dynamic recrystallization by subgrain rotation was active during deformation. These results confirm the dominant role of dissolution-precipitation creep and crystal plasticity in controlling slow aseismic deformation in these rocks.

Foliated fault rocks are frequently indicative of slow aseismic creep processes, and foliation development is most often associated with phyllosilicate-rich domains (e.g., Collettini et al., 2011). However, foliation also developed in the micritic limestone fault rocks at Gubbio, which are phyllosilicate-poor rocks (Fig. 8; Table 1). By comparing the silica (Si) content estimated via SEM microchemical analysis (Fig. 8D) with the composition of the montmorillonite and illite (http://webmineral.com), in agreement with the XRPD analysis (Figs. 8C and 8E; Table 1), we can conclude that clay content is <1%. This result is overestimated if we consider that part of the Si is quartz (Fig. 8D; Table 1; see Supplemental Item S5). Similarly, Smeraglia et al. (2017) have reported the occurrence of 100%-calcite foliated cataclasites, meaning that foliation is not exclusive of phyllosilicate-bearing rocks, and this calls for other controlling factors in its development. Such fabric development could be related to preexisting anisotropies, rather than just to clay content. For instance, in the Scaglia Rossa Formation, the basic structures constituting a mechanical layering and leading to foliated fault rocks with S-C fabrics could be (1) the original bedding (few tens of centimeters in thickness; Fig. 6A), (2) burial-related, bedding-parallel pressure solution seams, and (3) pressure-solution cleavage produced by initial layer-parallel shortening (e.g., Pierantoni et al., 2013). On the other hand, widespread foliation development in massive limestone fault rocks could (1) be prevented by the lack of an apparent inherited protolith texture and/or (2) need higher confining pressure conditions. For instance, Mercuri et al. (2018) observed pervasive and anastomosed (S-shaped) foliation in calcite-rich (>98%) fault gouges sheared at high normal stress (100–120 MPa) and subseismic rates, whereas fault gouges deformed at low normal stress (10 MPa) and at the same velocities showed localized B-R1 shear zones and no foliation.

We showed that fault rocks indicative of seismic and aseismic slip coexist, suggesting that the same fault may have a bimodal slip mode. We posit the possibility for carbonate-hosted normal faults to accommodate aseismic inelastic strain during the interseismic period. This deformation state is necessarily transient, and eventually the fault regains strength, and elastic strain accumulates with time until the subsequent seismic event occurs.

Transition between Aseismic and Seismic Slip

The conditions that determine the transition from fault creep to seismic rupture are still unclear. There has been some effort in addressing this issue for carbonate-hosted faults (e.g., Verberne et al., 2014; Passelègue et al., 2019). These experiments essentially considered variations in temperature and confining pressure to reproduce different modes of slip, implying that the transition occurs downdip. Interestingly, several lines of evidence, including our work, suggest slow aseismic phenomena and destructive seismic rupture occurring on the same fault patch (e.g., Hadizadeh, 1994; Di Toro and Pennacchioni, 2005; Kato et al., 2012; Smeraglia et al., 2017). This indication indicates that the frictional behavior of upper-crustal faults cannot be suitably described by a simple downdip layering of rate-strengthening and rate-weakening behavior regions.

If carbonate-hosted faults undergo interseismic creep in between seismic rupture events, there must be some mechanism that eventually accounts for fault strength recovery before the subsequent earthquake. Common processes include frictional healing, gouge compaction, pressure solution, cementation, plastic deformation at grain contacts, and resultant mechanical/chemical porosity reduction and growth of contact area (Bos et al., 2000; Renard et al., 2000; Marone and Saffer, 2015; Chen and Spiers, 2016). In our perspective, restrengthening would occur at low rates (e.g., Delle Piane et al., 2016; Park et al., 2021) and not in the immediate postslip period (e.g., Tesei et al., 2014). For instance, triaxial tests of fault reactivation conducted on simulated calcite-rich (~98%) faults at confining pressure up to 30 MPa, at room temperature, and under water-saturated conditions and subseismic slip rates revealed little to no time-dependent friction or cohesion strengthening (Delle Piane et al., 2016). Even lower healing rates can be expected if higher confining pressure and lower water-saturation conditions are considered. A sluggish interseismic healing could favor aseismic sliding of the studied faults. Eventually, aseismic release may be insufficient, and healing would become more significant, and elastic strain would accumulate over time (Fig. 11). The fault-slip mode is greatly affected by the competition between healing and stress relaxation mechanisms (Marone and Saffer, 2015). Aseismic sliding may be then suppressed and transition to unstable failure when the fault restrengthening rate overcomes the rate of stress relaxation resulting from creep and plastic deformation. When considering the faults hosted in micritic limestones, it is likely that fault healing and compaction are achieved by time-dependent and/or temperature-activated deformation mechanisms such as pressure solution, crack sealing, and subordinate gouge compaction (Fig. 11A; e.g., Bos et al., 2000; Renard et al., 2000). A fault hosted in competent massive limestones could experience interseismic creep if it is “ductile” at low stress, and it experiences strain hardening (fault strength ∝ fault slip) prior to failure, until reaching a failure threshold stress (e.g., Beeler et al., 2001). A reactivated carbonate fault allows for failure at relatively low shear stress (e.g., Delle Piane et al., 2016). These conditions, i.e., aseismic slip and strain hardening, can be met along faults with a developed gouge (i.e., after earthquake faulting). Such materials easily undergo cataclastic/granular flow (e.g., Renner and Rummel, 1996), which can accommodate slip aseismically (e.g., Sibson, 1977; Delle Piane et al., 2018) and which may cause strain hardening of the fault rock due to decreasing grain size, particle rearrangement, and inelastic compaction of pore space during increasing strain (e.g., Edmond and Paterson, 1972; Baud et al., 2000). That is especially the case when these changes are coupled with (compacting) pressure solution or crystal-plastic processes (e.g., Hadizadeh, 1994; Fredrich et al., 1989). Cataclastic flow may involve a complex interplay between dilatancy and compaction (Baud et al., 2000). Some experimental data show that cataclastic deformation is associated with strain hardening whether the sample is dilating or compacting (e.g., Schubnel et al., 2005). Aseismic cataclastic flow is presumably transient and may eventually lead to (localized) seismic failure at higher stress, when load-bearing force chains within the gouge layer catastrophically fail (Fig. 11A; e.g., Anthony and Marone, 2005). In addition, pore compaction and/or cement precipitation (cohesion strengthening) produce porosity and permeability reduction and a resulting increase in pore pressure within the fault gouge, which may eventually result in dynamic failure (e.g., Renard et al., 2000; Di Toro and Pennacchioni, 2005; Dal Zilio and Gerya, 2022). We cannot rule out the possibility that fluids play a role in fault-slip mode, since they may either induce aseismic creep or trigger seismicity (e.g., Avouac, 2015; Dal Zilio and Gerya, 2022). Recent experiments simulating induced seismicity have shown that fluid overpressure can stabilize fault slip, thereby promoting aseismic creep instead of triggering dynamic rupture (Scuderi and Collettini, 2016).

Implications for the Seismic Cycle

The key issue in our understanding of the seismic potential of faults and the earthquake cycle is the contribution of seismic versus aseismic fault slip to the tectonic loading rates (Avouac, 2015). Several lines of evidence exist of faults in the uppermost seismogenic crust accumulating little elastic strain (e.g., slip deficit) during the interseismic period. This may be explained by faults experiencing long-term creep (e.g., Bürgmann et al., 2000) or by distributed inelastic deformation (e.g., Fialko et al., 2005). Along the Apennines, a comparison of geodetic and seismic moment rates showed that seismic release accounts for only a fraction of the measured geodetic deformation, and the deficit is probably being released aseismically by long-term and transient fault creep and/or off-fault diffuse deformation (D’Agostino et al., 2009; Carafa et al., 2017, 2020). In this context, the interseismic behavior of an active fault remains largely enigmatic: Is it locked or creeping aseismically? Is it possible for a fault patch to switch from the former to the latter during the seismic cycle?

Several aseismic deformation mechanisms have been recognized that may release a significant fraction of energy and fault-slip budget in postseismic and interseismic periods (e.g., Perfettini et al., 2010; Avouac, 2015). For instance, interseismic creep has been documented along the ENE-dipping Alto-Tiberina fault, both along the fault plane and within its footwall block (e.g., Valoroso et al., 2017), an ~70-km-long, low-angle extensional detachment, to which the Gubbio fault connects antithetically at ~4–6 km depth (e.g., Mirabella et al., 2004). Cheloni et al. (2019) analyzed the rupture process and slip distribution associated with the 30 October MW 6.5 earthquake in central Italy. By applying their three-fault model, they documented that the total seismic moment released was greater than the estimates based only on seismic data, suggesting that the discrepancy might be related to either slow (or aseismic) or postseismic slip along the whole fault network. The absence of large aftershocks around secondary structures supports this hypothesis (Cheloni et al., 2017). In our study area, we found a seismic coupling ≈ 0.75 (Fig. 9), meaning that (at least) 25% of the long-term deformation is released aseismically in the upper crust. This result implies that aseismic slip must be added to seismic slip along fault traces to reconcile the long-term deformation rates.

A comparison of geodetic data, background strain rate, and the earthquake catalog for the Italian peninsula revealed that elastic energy accumulates in those areas where faults are locked and the strain rate is lower, so larger earthquakes occur with higher probability in zones of lower (interseismic) strain rate (Doglioni et al., 2011; Riguzzi et al., 2012). Intuitively, this could mean that regions characterized by relative higher strain rates are areas where no or little elastic strain is accumulated, which may result from an earthquake (e.g., after slip or viscoelastic relaxation of the lower crust) rather than preceding it (Hammond et al., 2010), and/or where faults are creeping. This supports our idea that faults may creep interseismically.

Our data suggest a phase in which faults undergo aseismic motion, but we cannot yet discriminate whether (1) fault creep (i.e., inelastic strain; Fig. 11B, green bar) initially dominates and slows down over time when it is suppressed by fault-healing processes (blue bar), and eventually elastic strain accumulates (Fig. 11B, blue bar), until the subsequent seismic event occurs (Fig. 11B, red bar) or if (2) fault creep (Fig. 11C, green bar) takes place along with elastic strain accumulation (Fig. 11C, blue bar), since it is unable to accommodate all the deformation. In the first case, the time intervals between consecutive earthquakes may be irregular (Fig. 11B) because the elastic strain (and stress) accumulation depends on the efficiency of fault-restrengthening processes and/or fluid pressure conditions suppressing aseismic sliding. This scenario poses additional questions about earthquake recurrence time and the regularity of the seismic cycle. In the second case, elastic strain accumulation, and hence the time intervals, may be nearly regular (Fig. 11C), unless there is interaction between faults (e.g., Tondi et al., 2021).

According to our model, faults in competent carbonate units undergo aseismic creep, and faults therein are not locked for the entire interseismic period, but only for a limited time between repeated events of short-term seismic rupture, postseismic slip, and long-term interseismic creep (Figs. 11B and 11C). The occurrence of aseismic slip on faults in competent limestone units, which are at present considered to be interseismically locked, requires a reevaluation of the seismic cycle and should be considered for more robust earthquake hazard assessment, since it may significantly contribute to the budget of moment release.

A seismic coupling of ~0.75 (upper bound) indicates that (at least) 25% of the long-term deformation in the study area is released aseismically in the upper crust. This implies that on-fault aseismic motion must be added to seismic slip along fault traces to reconcile the long-term deformation rates.

We propose that seismogenic normal faults hosted in competent carbonate units partially creep during the interseismic stage due to their velocity-strengthening behavior at low slip rates. Aseismic slip appears as (1) cataclastic flow accommodated by pressure solution and crystal plasticity, and CPO development in fault gouges developed in massive limestone and (2) pressure solution, frictional sliding on foliation surfaces, CPO development, and minor cataclasis in faults hosted in micritic limestone. During this creeping phase, minor elastic strain is accumulated in the adjacent rock volumes. We suggest that seismic and aseismic sliding can occur along the same carbonate-hosted fault and at the same structural level at different times during the seismic cycle, posing doubt about whether (1) seismic versus aseismic behavior takes place in distinct patches along/across the fault and (2) the lithology largely controls slip behavior. Other mechanisms must be invoked to explain coexisting and overlapping deformation structures related to both aseismic and seismic slip in the same lithology. In addition, we observed that creep is not exclusive to phyllosilicate-rich units. We envisage that different deformation mechanisms and resulting fabrics (foliated versus nonfoliated fault rocks) likely also depend on the inherited fabric (e.g., bedding) rather than just on lithology (e.g., clay content).

We cannot exclude the possibility that part of the aseismic deformation is accommodated by off-fault diffuse strain rather than on-fault creep. Eventually, fault locking and loading could result, according to the mechanism active during creep, through a combination of sluggish frictional healing and (slow) physical/chemical strengthening with progressive slip, i.e., slip-hardening cataclastic flow, pressure solution, plastic deformation at grain contacts, and cementation, which all result in pore reduction (although we cannot exclude the role of fluid pressure cycles). Our work provides new insights into the fault-slip modes in carbonate rocks and may profoundly impact current comprehension of the seismic cycle, fault seismogenic potential, and earthquake hazard assessment of faults in these terrains.

1Supplemental Material. Item S1: Analytical techniques used for fault-rock characterization. Item S2: Additional micrographs of fault rocks of the Monte Vettore fault. Item S3: EBSD data useful for the evaluation of subgrains and plastic deformation. Item S4: XRPD patterns and additional micrographs of fault rocks of the Gola della Rossa fault. Item S5: Additional micrographs and SEM microanalysis data of fault rocks of the Gubbio fault. Please visit https://doi.org/10.1130/GSAB.S.23638359 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Wenjiao Xiao
Associate Editor: Dongfang Song

Funding was partly provided by the University of Camerino FAR Project “Novel Approach for Seismic Hazard Analysis—NoHard” (E. Tondi). We thank two anonymous reviewers and the associate editor for constructive criticisms. We thank A.V. Brovarone for support with cathodoluminescence analysis.

Gold Open Access: This paper is published under the terms of the CC-BY license.