Models for hydraulic fracturing–induced earthquakes in shales typically ascribe fault activation to elevated pore pressure or increased shear stress; however, these mechanisms are incompatible with experiments and rate-state frictional models, which predict stable sliding (aseismic slip) on faults that penetrate rocks with high clay or total organic carbon. Recent studies further indicate that the earthquakes tend to nucleate over relatively short injection time scales and sufficiently far from the injection zone that triggering by either poroelastic stress changes or pore pressure diffusion is unlikely. Here, we invoke an alternative model based on recent laboratory and in situ experiments, wherein distal, unstable regions of a fault are progressively loaded by aseismic slip on proximal, stable regions stimulated by hydraulic fracturing. This model predicts that dynamic rupture initiates when the creep front impinges on a fault region where rock composition favors dynamic and slip rate weakening behavior.

( A ) Stratigraphic column from the study area (blue, limestone and dolostone; gray, shale and mudstone; orange, sandstone; pink, evaporites). ( B ) West-east slice through a 3D seismic volume (vertical exaggeration, 4.25; see Fig. 3 for location), showing the location of the M w 4.1 earthquake (red star). TVD, total vertical depth from surface. Several north-south–oriented faults can be distinguished (black lines), including one coincident with the location of the earthquake. C, possible channel visible as a “bowtie” structure; SR, edge of Swan Hills Formation reef structure. ( C ) Histogram of seismic event locations versus depth for 20-m bins (red asterisk, M w 4.1 earthquake; red dashed line, treatment well depth). ( D ) Estimated TOC + clay weight percentage derived from well log data.

To test this hypothesis, we analyze an unusually comprehensive dataset from central Alberta, Canada (“Duvernay,” Fig. 1A ). Since December 2013, several notable earthquakes in this region have been linked in space and time with hydraulic fracturing programs ( 29 ). The focus of our study is a M w 4.1 event on 12 January 2016, currently one of the largest magnitude seismic events induced by hydraulic fracturing ( 7 , 17 , 27 , 28 ). Hydraulic fracturing in this region is focused on the Late Devonian Duvernay Formation ( Fig. 2 ), which is composed of multicyclic units of organic-rich shales and fine-grained carbonates ( 30 ). Precambrian basement is inferred to be located ~400 m below the reservoir ( 27 ), while immediately above the Duvernay is the Ireton Formation, which is primarily composed of clays alternated with fine-grained carbonates ( 31 ). These units are conformably overlain by carbonate-rich rocks of the Winterburn and Wabamun groups ( Fig. 2A ).

Here, we consider a third hypothesis that distal dynamic rupture is induced by pore pressure–driven aseismic slip, which propagates from the reservoir depth out to the seismogenic region of the fault ( Fig. 1B ). The idea of an aseismic slip front outpacing a pore pressure diffusion front has been previously proposed ( 10 , 21 , 22 ), with recent laboratory (centimetric scale) and in situ observations at the decametric scale confirming that aseismic slip in pressurized regions of faults can activate slip along nonpressurized fault patches ( 23 – 26 ). Here, we apply the same model but at a much larger scale, i.e., to a kilometer-scale fault rupture that has been linked to hydraulic fracturing ( 27 , 28 ).

Several lines of evidence indicate that earthquake nucleation due to hydraulic fracturing is vertically offset from the injection zone. High-resolution monitoring observations indicate that hypocenters of hydraulic fracturing–induced earthquakes do not coincide with the injection depth ( 14 – 17 ). Moreover, laboratory measurements and in situ experiments indicate that slip on a shale-hosted fault is expected to be stable ( 18 – 20 ). This stability arises, in part, from the rate-state frictional parameters of rocks with high TOC and clay ( 19 ). Together, these observations imply that if the first hypothesis is correct, then dynamic rupture arises from along-fault diffusion of pore pressure, whereas the second hypothesis favors rupture nucleation relatively close to the injection zone.

Several models have been proposed to explain the mechanisms of fault activation by fluid injection. The first hypothesis is an increase in pore pressure within the fault zone, which leads to a reduction in effective normal stress acting on the fault ( Fig. 1B ) ( 7 , 9 , 10 ). Increased pore pressure by hydraulic fracturing is expected because (by design) the fluid pressure within a hydraulic fracture must exceed the in situ minimum stress; however, pressure diffusion away from an induced fracture is strongly inhibited by the low permeability of shales ( 2 , 3 ). This fault activation mechanism therefore requires the existence of a hydrologic connection between a hydraulic fracture system and a fault, either by direct intersection ( 11 ) or, more tortuously, via naturally occurring fractures. An alternative hypothesis posits that fault activation is caused by stress changes arising from poroelastic coupling between hydraulic fractures and the rock matrix ( 12 ). This mechanism is capable of altering fault-loading conditions without any hydraulic connection ( Fig. 1B ), but the magnitude of the stress change is relatively small and diminishes with distance ( 13 ).

Hydraulic fracturing is a process of injecting pressurized fluids into a rock formation to create open fracture systems and is used extensively in the development of low-permeability hydrocarbon reservoirs ( 1 ). Such fine-grained reservoir rocks (often referred to as “shales” for simplicity because of their fine-grained nature) are typically characterized by exceedingly low matrix permeability, typically on the order of 10 −18 to 10 −21 m 2 ( 2 , 3 ), in addition to elevated clay and total organic carbon (TOC). Hydraulic fracture growth is generally accompanied by low moment magnitude (M w < 0) microseismic events ( 1 ), and the occurrence of potentially damaging earthquakes was, until recently, considered unlikely ( 4 ). Hydraulic fracturing–induced earthquakes up to M w 4.8 have now, however, been documented in China, North America, and the United Kingdom ( Fig. 1A ) ( 5 – 8 ).

( A ) Snapshot of overpressure on a vertical fault intersected by hydraulic fractures at the shut-in of stage 23 (numbered black rectangles, stages 20 to 23). Model includes pre-existing overpressure field that developed during hydrocarbon maturation. ( B ) Evolution of effective normal stress, averaged over the along-strike fault section centered on x = 0. Dashed contours show long-term stress evolution (adjusted to the present-day depth), while solid contours show the short-term increase due to hydraulic fracturing. Stage 27, end of stage 26. ( C ) Fault slip shown by gray dashed lines every 10 million years (Ma; long-term geologic creep), by continuous gray lines every 5.5 hours starting from the beginning of stage 23 (accelerated creep), and by red continuous lines every 0.01 s (coseismic slip). Total hydraulic fracturing–induced aseismic slip before the earthquake nucleation shown by the shaded area. ( D ) Moment release as a function of time; onset of accelerated aseismic slip at the beginning of stage 23 leads to a M w 4.1 earthquake early into stage 24.

The spatial and temporal evolution of seismicity leading up to the mainshock, coupled with the interpreted fault geometry and log signatures, are consistent with upward propagation of an aseismic deformation front, followed by initiation of dynamic rupture where the deformation front impinges on unstable regions of the fault ( Fig. 4 ). To simulate slip activation and development on the near-vertical north-south fault observed in the 3D seismic and microseismic data, we performed 1D rate-state friction modeling ( Fig. 5 ) ( 33 , 34 ) constrained using input parameters from a range of available data (fig. S3 and Materials and Methods) and the hydraulic fracturing stage schedule (fig. S4 and table S1). The model was initialized by considering long-term fault creep due to reservoir overpressure generation on a time scale of tens of millions of years ( 35 ). Observed overpressure is mainly confined to the frictionally stable shales (figs. S3 and S5), resulting in preloading of fault stress in the adjacent strata (fig. S5). During hydraulic fracturing, the model shows that accelerated aseismic fault creep results from pressure increases where the hydraulic fractures intersect the fault. The stress transfer arising from aseismic slip substantially outpaces the pore pressure diffusion front, even considering pore pressure–induced permeability increase on the fault ( Fig. 5B and Materials and Methods). The nucleation of dynamic rupture occurs once aseismic slip has propagated ~50 m into the locked carbonate units (~250 m above the fractured well), close to the location of maximum negative (a − b) [friction rate parameter ( 19 ) estimated from the TOC + clay results]. In the numerical simulation, a M w 4.1 earthquake is generated at the beginning of stage 24, exhibiting magnitude and timing that are similar to the observed event (see Materials and Methods for a detailed description of the dynamic rupture propagation). Essential to this coseismic slip model is the propensity of carbonate fault gouge to undergo a near-complete dynamic weakening, owing to a process such as flash heating–induced thermal decomposition of calcite at sliding microasperities (see Materials and Methods) ( 36 – 38 ) or superplastic deformation (fig. S6) ( 39 ). Numerical simulation of slip without dynamic weakening leads to a much smaller event ( Fig. 5D and fig. S7).

Timing of the M w 4.1 earthquake is highlighted in red. Treatment well (red line), stage locations (black triangles), stratigraphy (black dotted lines), and estimated positions of the steeply dipping fault plane and region of increased pore pressure due to hydraulic fracturing (HF) are shown. Creep outpaces the pore pressure diffusion (PPD) front on the fault plane. The mainshock is believed to be located on the same fault plane as the other events but is slightly mislocated horizontally because of the lower frequency content.

Figure 3 shows a map of seismicity recorded during the hydraulic fracturing treatment associated with the M w 4.1 event (see fig. S1 for more details and for cross sections); seismic events are overlain onto a 3D seismic time slice through the mid-Devonian Beaverhill Lake Group, which forms a substrate to the reservoir and contains a reef margin that marks a lithologic boundary between shale and carbonate (see Materials and Methods). The event distribution outlines several rectilinear trends that correlate with linear features in the 3D seismic image ( Fig. 3, A and B ). These trends are very similar to those observed using a dense microseismic monitoring array ( 17 ). We interpret these linear features as anastomosing strands within a regional strike-slip fault system ( 17 , 27 , 28 ). Apparent truncation of the linear features by the reef edge implies that only limited (<50 m) net slip has accumulated on these faults since reef growth in the mid-Devonian. The seismic event distribution is uncharacteristically asymmetric and is concentrated on the east side of the well, suggesting the presence of a lateral stress gradient due to the influence of nearby faults ( 32 ). Figures 3C and 4 show updip progression of seismicity leading up to the M w 4.1 earthquake, culminating in a mainshock and aftershock sequence at a depth of ~3000 m ( Fig. 4 shows propagation upward along the fault plane). The mainshock occurred on a north-south–oriented fault, which is parallel to the treatment well ( 17 ). They are both oriented ~45° from the maximum horizontal stress ( Fig. 3 ) ( 27 ); thus, there is substantial resolved shear stress on the fault [treatment wells are often drilled perpendicular to the maximum horizontal stress ( 1 ); however, in this case, the well orientation was determined by the shape of the land lease boundary].

( A ) Depth slice through the 3D seismic volume showing normalized seismic reflection amplitude at the level of the Beaverhill Lake Group (at 3386-m depth; Fig. 2 ), which directly underlies the Duvernay Formation (gray line, treatment well location; black solid line, profile in Fig. 2 ; red star, M w 4.1 earthquake). Various structural features can be interpreted, including the edge of the Swan Hills reef structure and a series of roughly north-south lineaments. Seismic event locations are overlain, colored by depth. ( B ) Enlargement of (A) around the treatment well [triangles, injection stages (numbered at intervals)]. Focal mechanism of the M w 4.1 event ( 27 ) is shown as a beach ball. Approximate orientation of maximum horizontal stress is labeled (S Hmax ); hydraulic fractures should grow parallel to this direction ( 1 ). ( C ) Depth of the seismic events versus time for the events in the immediate vicinity of the mainshock [red dashed lines, stage start times; black dotted lines, stratigraphy according to a well log ~3 km north (BhL, Beaverhill Lake; Dv, Duvernay; Irtn, Ireton); red star, M w 4.1 earthquake].

Figure 2 shows a profile extracted from the three-dimensional (3D) seismic volume through the location of the M w 4.1 earthquake (see Fig. 3A for profile location). Several subvertical faults are identified (see Materials and Methods). One inferred basement-rooted fault strand coincides with the location of the mainshock, cutting through both the injection zone and the seismogenic region. A graph of TOC + clay weight percentages ( Fig. 2C ), derived from nearby well log data (fig. S3), shows a sharp transition from high values (>30%, corresponding to stable sliding regime) in the Duvernay and Ireton to approximately zero (corresponding to potentially unstable sliding) at shallower depths. Similarly, a depth histogram of seismicity ( Fig. 2C ) shows that most of the seismic activity occurs within the unstable sliding Wabamun and Winterburn carbonate units.

In a broader context, the apparent association of hydraulic fracturing–induced seismicity with faulting in carbonates is consistent with observations of some of the largest natural earthquakes occurring within the sedimentary cover worldwide ( 42 ). Similarly, interconnected, vertically stacked creep and dynamic rupture processes occur at separate depths on some faults such as the San Andreas ( 43 ), where aseismic creep has also been hypothesized as the driving process of natural earthquake swarms ( 44 ). Considering the wealth of detailed geological and operational information acquired by the hydrocarbon industry in areas of injection-induced seismicity, our findings could therefore also be beneficial to understand other induced and natural earthquake systems.

Improved understanding of fundamental processes of fault activation during hydraulic fracturing is key to developing effective monitoring and mitigation strategies and could also help to inform our understanding of natural earthquake triggering. Current mitigation strategies for hydraulic fracturing–induced seismicity in this region are based on the use of a traffic light protocol ( 40 ) [traffic light protocols are feedback systems that allow for an operational response to the nearby occurrence of seismic events exceeding a prescribed set of criteria ( 41 )]. The existing system may not be optimal, as there is an implicit assumption that a large-magnitude earthquake is preceded by smaller precursory events and that changes in injection operations would have an immediate effect on the source process of induced events. According to our model, mitigation of induced seismicity may benefit more from characterizing nearby faults, state of stress (including overpressure), and rock properties (TOC, clay, and carbonate content). Furthermore, the maximum magnitude of an event may be limited by factors that include geometrical considerations (thickness of the unstable layers and vertical distance from injection zone), geological overpressure history, and slip deficit accumulated over the epoch of stable creep in the adjacent stable formations. Deployment of new technologies for direct detection of aseismic slip (such as fiber-optic distributed strain sensors) could provide a useful approach to complement traditional monitoring methods.

MATERIALS AND METHODS

Computing TOC + clay The TOC (total organic carbon) + V clay (volume of clay) weight percentages were estimated from spectral gamma ray (SGR) logs. SGR uranium content is tightly related to the TOC (45, 46), and the SGR thorium reading is a function of the amount of the clay volume present in the reservoir (46, 47). These well log–estimated empirical correlations were calibrated using the closest well with available Duvernay core for analysis, ~17 km southeast from the treatment wellhead. The cored well has available SGR logs and 29 core plugs with measured TOC (from Rock-Eval) and V clay [from x-ray diffraction (XRD)] within the Duvernay Formation. For the correlation of SGR uranium content to TOC (in weight %), a linear function using 0.0033 as conversion factor was applied [i.e., TOC uranium = uranium reading parts per million (ppm) × 0.33; fig. S2A]. This ensured the best correlation between the modeled TOC using well logs and the TOC directly measured on 26 core samples (R2 = 0.73). For V clay estimation, a cross-plot of SGR thorium reading versus measured V clay (from XRD) was built for 20 core plugs taken within the Duvernay Formation (fig. S2B). In the analyzed core, V clay ranged from 3 to 45%, and thorium reading ranged from 3.5 to 13 ppm. The trend line fitting the 20 data points intersects 0% clay at 3 ppm (Th min ) and 100% clay at 21 ppm (Th max ). On the basis of this, we estimated V clay weight % = (Th reading − Th min )/(Th max − Th min ) with R2 = 0.58.

Reflection seismic data Migrated 3D seismic reflection data from the same acquisition program used in this study were discussed in (48). As part of the processing of the reflection seismic data, amplitude normalization was performed, such that the highest observed reflection amplitude had a value of unity. Geologic horizons shown in Figs. 2 and 3 were based on a calibration well that was drilled to a total vertical depth of 3607 m from a surface elevation of 890.9 m (48). The top of Precambrian basement is deeper than local well control and so was correlated using a nearby Lithoprobe deep reflection survey, which has synthetic ties to Precambrian well penetrations (49). Images extracted from the 3D seismic volume can be used to investigate structural features within the region at various depths. We first flattened the reflections at the Wabamun Formation top to highlight geological features that were present before and during deposition of this layer. Figure 2B shows a cross section through the 3D seismic volume displaying seismic amplitudes as a function of depth, which can also be used to visualize structural features, as well as stratigraphy. Figure 3 (A and B) shows a horizontal slice through the 3D seismic volume at a depth of 3386 m, with the positive and negative seismic amplitudes highlighting various structural features discussed below. Although cumulative fault displacement is inferred to be dominantly perpendicular to the horizons, subvertical faults can nevertheless be identified on the basis of localized curvature anomalies (50) or changes in apparent dip of reflections that mark the locations of hinge lines (49). Channel and reef edge features can also be identified in Fig. 2 in a similar manner; however, the interpreted fault strands span multiple stratigraphic horizons and correlate with the linear features in Fig. 3. Our estimate for maximum post–mid-Devonian fault displacement was constrained by the lack of a discernible lateral offset on the reef edge. Figure 3A shows that there is little discernible offset of the mid-Devonian Swan Hills reef margin where it is intersected by the observed faults. Some seismicity forms a cluster aligning parallel to the reef edge (17), suggesting that higher mechanical strength of the reef facies relative to the off-reef facies may have influenced faulting. The theoretical resolution limit for migrated seismic data is equal to the spatial sampling interval; thus, based on the common midpoint interval for the data of 30 m, we estimated that the offset on the reef margin can be no greater than 50 m (allowing for some uncertainty).

Broadband seismic data The hydraulic fracturing operations associated with the M w 4.1 event were monitored using a small network (WSK) of four industry-provided broadband stations, located between 2.4 and 7.6 km from the treatment well. The continuous ground motion recordings from these sensors are publicly available for the period between 00:00 UTC on 1 January 2016 and 1:56 UTC on 19 January 2016. A catalog containing 108 events was detected during this period, with a minimum event magnitude of M L 0.3. Pick times were manually extracted for P- and S-wave onsets, and a grid search was performed to minimize the difference in differential theoretical travel times between each grid point and each pair of receivers with observed differential pick times (51). Theoretical travel times were obtained using a simple ray-tracing tool and the same 1D velocity model that was used to locate the microseismic catalog (17), which was calibrated using perforation shots [explosions used to create holes in the well casing at each stage to allow egress of fracturing fluid into the surrounding medium (1)] with known locations. A 1D velocity model is a reasonable assumption given the relatively horizontal stratigraphy observed in Fig. 2. The grid search was confined to a 8.2 km by 8.5 km by 5 km volume with a 25-m resolution (in east-west, north-south, and depth, respectively), centered around the treatment well. After obtaining the initial locations, the events were subsequently relocated using HypoDD, a double-difference relocation algorithm (52). To obtain high-precision differential pick times between events at each station, the waveforms were band-pass–filtered between 1 and 20 Hz. A 2.5-s window was taken around the catalog pick time (starting 0.05 s before the catalog pick time). For every station, differential pick times were calculated as the lag times associated with the maximum correlation coefficients for the windowed P and S phases for all events cross-correlated with the corresponding windowed phases for all other events. The pick times of each phase were adjusted accordingly, and the square of the maximum correlation coefficient for each phase was used as a weight in the HypoDD algorithm (52). The results (fig. S1) show locations that are consistent with the microseismic locations in (17), and therefore, the same structural features are apparent. The hypocenter depths are also very similar; differences in the apparent dip of the southern imaged structures are likely a manifestation of the station coverage in the broadband network. The one inconsistency is that the M w 4.1 event is slightly offset to the west of a linear cluster of hypocenters in the northern region. Given that the microseismic data show the hypocenter of this event aligning with this cluster and that the frequency content of this event differs from that of the other events that are substantially smaller, which may adversely affect the phase correlations used in the HypoDD algorithm, we believe that this offset is caused by an error in the location of the mainshock. The true hypocenter likely aligns on the same fault strand as the observed cluster of seismicity.

Model fault geometry, ambient stress, and pore pressure A north-south striking subvertical strike-slip fault intersecting the Duvernay Formation and other neighboring stratigraphic units is indicated by reflection seismic data and relocated microseismicity (Fig. 2 and fig. S1). The fault is assumed to terminate in the carbonate (Wabamun) formation updip at the shallowest recorded aftershocks of the M w 4.1 event. Downdip assumed that locking depth is at the basement contact (Fig. 2) (note that this precludes simulation of possible earthquakes in the basement extension of the fault). The total fault height capable of slip is 1017 m. The regional maximum horizontal in situ stress direction is northeast-southwest (27), at ~45° to the north-south striking fault. The minimum stress gradient of 20.78 kPa/m in the proximity of the well, and, by extension, in the proximity of the fault striking parallel to the well ~100 m to the east, is inferred from the hydraulic fracture stage 23 data (see fracture closure stress; table S1) and is within the range of 19 to 22.4 kPa/m suggested previously for the Kaybob Duvernay play (53). The same study estimates maximum stress in the Duvernay to be in the range of 27.1 to 31.7 kPa/m with the median value of 29.5 kPa/m. Corresponding fault normal and shear stresses (neglecting stress rotations near the fault) are σ n = (σ max + σ min )/2 = 25.14 kPa/m and τ = (σ max – σ min )/2 = 4.36 kPa/m. The resulting nominal value of the fault stress ratio τ/(σ n − p hydrostat ) = 0.29 is substantially lower than the nominal friction values of ~0.5 to 0.6 of shale and carbonate-rich fault gouges, and, if taken for the far-field fault loading, would preclude any notable fault slip outside the overpressured reservoir fault section. This deviates from a simple model of a critically stressed crust (54); however, we remark that (i) the faults appear to terminate at the top of the Late Devonian Wabamun Formation (Fig. 2) and therefore represent relic structures that formed at a time when the stress field was likely different from the present day and (ii) sedimentary basins exhibit geomechanical layering with variable stress between stronger and weaker layers, such that if the crust is critically stressed at some depths, it cannot be critically stressed throughout the sedimentary section (55). The low estimated value of the fault stress ratio in the Duvernay is consistent with the shear stress drop there because of the long-term slip along the creeping sections of the fault through overpressured shales or sandstones. Naturally occurring fault creep is the result of the long-term reservoir overpressure associated with the hydrocarbon generation coupled with the stable frictional rheology of the shale fault sections (see discussion of the fault friction and pore pressure regime below). Furthermore, this long-term slip is expected to scale the potentially seismogenic slip of the adjacent locked sections of the fault (through frictionally unstable carbonates and/or basement rocks) when triggered by anthropogenic fluid injection. For the M w 4.1 event, rupture dimensions of ~0.4 km (dip) and ~1 km (strike) are estimated, as outlined by the spatial distribution of the early aftershocks (Fig. 3 and fig. S1), and the shear modulus of the rock of μ = 25 GPa, resulting in an estimated average seismic slip 〈δ〉 of ~ 0.15 m. Using this seismic slip as the proxy for the long-term accumulation of aseismic slip over the adjacent creeping fault section of estimated height h ~ 400 m (twice the depthwise distance from the center of the overpressured Duvernay to the contact with the seismogenic and interseismically locked carbonate), we can estimate the average long-term stress drop along the creeping section Δτ = (2/π) μ 〈δ〉/h ~ 11.9 MPa (or 3.5 kPa/m). This estimate translates to the background shear stress gradient τ b = τ + Δτ = 7.9 kPa/m (fig. S3A) and the nominal value of the background fault stress ratio ~ 0.52. For the ambient pore pressure distribution, we considered an overpressured Duvernay Formation sealed from the formations above and below, where a normal (hydrostatic) pore pressure distribution was modeled. This is consistent with reported drilling experiences above the Duvernay (i.e., in the Ireton) (53). The Duvernay overpressure was modeled as a pore pressure gradient of 17.7 kPa/m [within the window of 15.8 to 18.1 kPa/m (53) and the median value reported for seismogenic Duvernay plays (35)]. Expected transition zones (seals) between the overpressured Duvernay and normally pressurized formations above and below have been suggested to account for potential overpressure-related drilling problems in the Ireton (53) and potential fault-valve action and episodic release of Duvernay overpressure in the bounding formations along vertical strike-slip faults (35), such as the one studied here. We modeled these as the tapered zones between the overpressured Duvernay and normally pressured surrounding formations (fig. S2A).

Fault friction Fault slip takes place in response to changes in the fault frictional strength τ f = f (σ n − p). Evolution of the friction coefficient f was modeled using the slip law formulation of the rate-state theory (56, 57) in the form of (58), which specifies friction f = f o + a ln ( V / V o ) + Θ (1)as a function of the slip rate V = dδ/dt and frictional state variable Θ. Here, f o and V o are reference values, and a is the direct effect coefficient. The evolution of the fault state is governed by the “distance” of the friction f from its steady-state value f ss (V) as follows d Θ / dt = − ( V / L ) ( f − f ss ( V ) ) (2)where L is the state evolution distance. At the steady state, a conventional logarithmic dependence on the slip rate, corroborated experimentally at low slip velocities, was used, i.e., f ss (V) = f LV (V) with f LV ( V ) = f o + ( a − b ) ln ( V / V o ) (3) The parameters of the rate-state frictional model were constrained using laboratory friction data for reservoir shale rocks (19) and carbonates (42, 59), as well as the TOC + clay data from well logs (Fig. 2D and Materials and Methods). Duvernay shale frictional properties (also used for clastic formations downdip of Duvernay) were inferred on the basis of the friction correlation (19) to the clay + TOC gouge content. The TOC + clay content distribution interpreted from the well log analysis was highly irregular through the shales (Duvernay and Ireton), yet we observed to the first order a decreasing TOC + clay content in the updip direction toward the contact with carbonates (above which the TOC + clay was null). As a first-order model, we took the average TOC + clay content values of ~33 to 34% over the Duvernay and clay-rich downdip part of Ireton (<3318 m) and used the corresponding correlation value (19) of a − b = 0.0025 to model the steady-state friction rate dependence uniformly over the Duvernay. For the carbonate-rich gouge of the Wabamun and Winterburn formations, we used a − b = −0.0025 updip of the contact with Ireton shales, as inferred from the hydrothermal data (59) at a representative temperature of ~110°C. The final adopted depth distribution of the a − b parameter is a linear taper between the above two end-members, such that neutral velocity dependence of friction (a − b = 0) takes place at the contact between shales (Ireton) and carbonates (fig. S3B). Reference value of friction coefficient at reference slip rate V o = 1 μm/s was set at f o = 0.5 in the Duvernay (19) and at f o = 0.6 in the carbonates (59) above and in the sandstones below, while a linear taper was used in between these end-members (fig. S3B). The direct effect coefficient was modeled similarly on the basis of the adopted end-member values of a = 0.01 (shales and sandstones) (60) and a = 0.02 (carbonates) (42). For the state evolution distance, we used a uniformly distributed value of L = 100 μm (19, 59, 60).

Strong weakening of carbonate fault gouge at high slip velocity When modeling seismic slip, we recognized the propensity of calcite-rich gouges for flash heating–induced thermal decomposition (36, 37, 61), specifically decarbonation (CO 2 release) and deposition of amorphous carbon (nano-lubrication), resulting in a dramatic drop of the steady-state friction (<0.1), with the slip velocity approaching seismic values. The rate-state frictional framework can accommodate the entire spectrum of frictional responses, from low to high slip velocity, when the steady-state friction expression is amended, e.g., as f ss ( V ) = f w + ( f LV ( V ) − f w ) ( 1 + ( V / V w ) m ) − 1 / m (4)where f w is the flash-heated (minimum) steady-state value of the friction coefficient and V w is a characteristic weakening velocity value, such that the steady-state friction dependence on the velocity reduces to the end-members f LV (V) and f w when V ≪ V w and V ≫ V w , respectively. Exponent m ≥ 1 in Eq. 4 is a fitting parameter, with larger values corresponding to sharper transition between the low-velocity lnV and high-velocity ~1/V dependence of the friction on V [the limit of m = ∞ corresponds to the abrupt transition between the two asymptotic behaviors (62)]. We fit the friction versus slip rate data from the laboratory slip on a marble saw cut (37) with Eq. 4 to infer f w = 0.07, m = 6, and V w,lab = 0.28 m/s (fig. S6, dashed line), while using the low-velocity friction f LV (V) (Eq. 3) characterized by f o = 0.6, V o = 1 μm/s, a = 0.02, and a positive value for the velocity dependence coefficient (a − b) lab = 0.0025 for calcite at room temperature (59) (as opposed to a − b = −0.0025 in situ). The flash heating theory (38) predicts that V w scales with the square of the difference between the weakening temperature T w for the onset of thermal decomposition and ambient temperature T. This allows the in situ value of the weakening velocity to be inferred from the laboratory value, V w = V w,lab (T w − T)2/(T w − T lab )2 = 0.2 m/s, when T lab = 20°C, T = 110°C (in situ temperature), and T w = 600°C [temperature at the onset of decarbonation (37)] are used. The corresponding dependence of the steady-state friction of calcite on the slip velocity under the in situ conditions is shown in fig. S6.

Fault slip modeling To model the fault slip, we solved equations of quasi-dynamic elasticity, which relate the change of fault normal and shear tractions to slip (34). σ ( z , t ) – σ n ( z ) = 0 (5) τ ( z , t ) – τ b ( z ) = ( μ / 2 π ) ∫ ( ∂ δ ( z ′ , t ) / ∂ z ′ ) d z ′ / ( z − z ′ ) – ( μ / 2 c ) ∂ δ ( z , t ) / ∂ t (6)where μ is the elastic shear modulus and c is the shear wave velocity. The shear traction on a slipping fault is given by the fault frictional strength τ ( z , t ) = f ( z , t ) ( σ n ( z ) – p ( z , t ) ) (7)and the evolution of the rate-state–dependent friction is described by Eqs. 1 to 4. We used the spectral method (63) applied to a finite fault of height H without replication (64) and fast Fourier transform to solve Eq. 6 on a grid of N nodes z j = z downdip_end + (j −1) Δz (j = 1,…,N) with uniform spacing Δz = H/N and reduce the problem to a system of dynamic ordinary differential equations (ODEs) for the values of slip at the nodes. Earthquake slip instabilities result in slip rates vastly variable in space and time, rendering the ODE system stiff (33). We overcame the numerical difficulty of solving a stiff system by changing the dynamic variable (time) to the fault-average slip. This allows the transformed ODEs to be solved using a routine solver with adaptive step in Wolfram Mathematica software.

Space discretization, convergence of the numerical solution, and creep versus pore pressure diffusion front dynamics We used the numerical stability considerations (33, 34) for slip on the frictionally unstable (a − b < 0) part of the fault to constrain the spatial discretization Δz ≪ h* = μL/max[(b − a)(σ n – p)] ~ 20 m on the basis of frictional properties of the carbonate units of the fault (fig. S3B). We note that the above criterion does not place any constrain on the grid spacing from the frictionally stable (a − b > 0) part of the fault (through shales and sandstones). We therefore investigated the stability and grid dependence of the numerical method when applied to a simplified setting of aseismic slip on a frictionally stable fault (characterized by the parameters for the Duvernay section of the full fault model) subjected to a constant ambient effective normal stress of 50 MPa. In this numerical example, the fault, initially in the steady state of uniform creep at a rate of 10−10 m/s, undergoes an episode of accelerated aseismic creep driven by 1D pore pressure diffusion p(z,t) − p o (z) = Δp max Erfc(|z|/√4αt) with hydraulic diffusivity α = 0.05 m2/s from a constant overpressure Δp max = 25 MPa source at z = 0. Figure S8 shows the development of the effective normal stress (due to the diffusive pore pressure change), shear stress change, slip rate, and slip along a finite fault of height H = 200 m (−H/2 < z < H/2) in the numerical solution with a large number of grid points N = 29 (fine spatial resolution Δz = 0.39 m). The symmetric crack-like expansion of the accelerated aseismic creep is seen to substantially outpace the pore pressure diffusion, which is quantified in a more direct fashion by tracking the creep and diffusion fronts in time (fig. S8E). Here, the two fronts have been defined as the location of the peak shear stress and the location |z| = 1.3859 √4αt, where pore pressure change is 5% of Δp max , respectively. The runout distances of the two fronts in this example are different by a factor of ~2.5 after 5 hours of fluid injection. Carrying out the numerical solution with a reduced spatial resolution, we observed the grid independence of the results for N ≳ 28 or Δz ≲ 0.8 m (fig. S8), while artificial slowdown of the creep front was evident for a coarser grid. On the basis of this example, we can conclude that, although the numerical solution for slip on the frictionally stable part of the fault remains numerically stable (obtainable) when carried out on a coarse grid, there exists a threshold spatial resolution below which the solution is grid dependent. The latter has to be considered in conjunction with the criteria in (33, 34), when solving for slip on faults with heterogeneous friction (stable and unstable parts). The grid independence of the numerical solution for slip in our full fault model with heterogeneous friction (fig. S3B), as detailed further, has been established by direct computations of the model for N ≳ 211 (Δz ≲ 0.4 m). This spatial resolution threshold is required to fully resolve the time progression of the aseismic slip on the frictionally stable (shales) part of the fault induced by the hydraulic fracture stimulation, while the condition for grid independence of the solution for the slip on the frictionally unstable (carbonates) part of the fault is less restrictive (i.e., the slip there could, in principle, be adequately computed on a coarser grid).

Long-term creep due to the reservoir overpressure generation Given the stable frictional rheology of shales or sandstone fault units, we anticipated that those may have undergone long-term geological creep driven by the fault stress and pore pressure loading. Specifically, we expected that the development of the current-day reservoir pore overpressure, attributed mainly to hydrocarbon generation over tens of millions of years (35), may have amplified the amount of creep in the shale section and, consequently, the stress redistribution along the fault leading to the stress concentration in the locked (frictionally unstable) fault units. We modeled the geologic creep on the fault by applying the background shear and normal stresses, as corroborated in the above, and initially hydrostatic pore pressure at t = −50 million years (Ma; with t = 0 corresponding to the present day) and by ramping up the pore overpressure linearly in time from zero to the current-day profile (fig. S5A). The resulting history of fault slip accumulation (fig. S5D) shows a maximum cumulative current-day slip of ~0.15 m coincident with maximum overpressure (Duvernay), tapering to zero updip into the “locked” carbonate (Wabamun/Winterburn) formations and downdip at the basement contact. Corresponding shear stress redistribution (fig. S5B) leads to unloading in the creeping shales and development of stress concentration in the locked carbonates. By performing additional calculations with earlier initialization time (t < −50 Ma) while maintaining the same overpressure generation time window (−50 Ma < t < 0), we verified that the current-day stress and slip distributions are practically independent of the fault slip history before the overpressure generation.

Fault pore pressure transient due to the hydraulic fracturing perturbation Hydraulic fracturing was completed in multiple stages along the north-south horizontal well with an approximately uniform spacing of ~85 m (fig. S1). In our model, each stage was hypothesized to be a single vertical planar fracture aligned with the direction of maximum horizontal stress (fig. S1B) (27) and contained in height to the reservoir (Duvernay) layer. The hydraulic fractures were assumed to intersect the nearly vertical strike-slip fault inferred some 100 m east of the well. Given the exceedingly small reservoir shale permeability, the fault damage zone (a tabular zone that is approximately meters to tens of meters thick of fractured and/or brecciated rock about the fault core) (65) is the most likely candidate to transmit the high-pressure hydraulic fracturing fluids (66) once intersected by a hydraulic fracture. However, the ambient permeability of the damage zone of a dormant fault is also likely to be sealed, in agreement with the maintenance of high ambient overpressure in the Duvernay (35, 53). Thus, some slippage of the complementary fractures and faults in the damage zone would be required to refresh the fault damage zone permeability (66). This permeability refreshment can be borne by stress changes in the damage zone induced by the slip on the principal fault plane (within the fault core) or due to reduction of the frictional strength in response to the pore pressure perturbation by the hydraulic fracture. Here, we adopted the latter concept and modeled the fault damage zone diffusivity [permeability/(fluid storativity) × viscosity] as a function of the pore pressure perturbation from its ambient value, increasing from a very low, ambient value (5 × 10−5 m2/s) to the “refreshed” value of 0.05 m2/s when pore pressure perturbation reaches or exceeds the 2-MPa threshold (fig. S4A). The choices of the ambient and “refreshed” diffusivity values are motivated by the in situ measurements in deep boreholes on two major faults, Chelungpu (67) (Taiwan) and Longmen Shan (68) (China), respectively. These two faults may provide examples of healed (low permeability) and refreshed or partially healed (high permeability) fault damage zones at the respective times of the measurements, which took place ~6 and ~2 years into the healing process after a major earthquake, respectively. Given the proximity of the intersections of hypothesized hydraulic fractures in stages 23 and 24 with the fault to the epicenter of the M w 4.1 event (fig. S1B), we modeled the 2D pore pressure diffusion along the fault damage zone with pressure-dependent diffusivity from the temporal progression of these and few neighboring stages’ (stages 20 to 26 separated by ~5.5 hours in time and ~85 m in space) intersections with the fault (fig. S4). Each hydraulic fracture was assumed to intersect the fault instantaneously upon the start of fracturing fluid injection into a given stage (as hydraulic fracture propagation time scale is generally small compared to the stage pumping and shut-in duration). The instantaneous shut-in pressure value of ~50 MPa (table S1) was used to approximate the fluid pressure at the hydraulic fracture fault intersection during the stage pumping time of ~2.5 hours, while the condition of no fluid exchange between the hydraulic fracture and the fault was imposed once the stage was shut-in (fig. S4B). The former assumption is likely to somewhat overestimate the actual fluid pressure at the hydraulic fracture fault intersection (due to the viscous pressure drop in the fracture with the distance from the well), while the latter is likely to temporarily underestimate the pressure due to the diminishing but nonzero fluid leak-off from the fracture into the fault after the shut-in. Together, these assumptions provide a first-order approximation of the boundary conditions at the hydraulic fracture fault intersections. The numerical solution of the pore pressure diffusion equation with outlined transient boundary conditions at the hydraulic fracture fault intersections and no-flow conditions imposed at the model far-field boundaries was obtained using the method-of-lines solver in Wolfram Mathematica software. The solution shows a coupled dynamics of diffusion of the pore pressure perturbation from the hydraulic fracture fault intersections, in turn, setting off the diffusion of the ambient (long-term) overpressure by “refreshing” the fault permeability.