Fluid-induced aseismic fault slip outpaces pore-fluid migration

See allHide authors and affiliations

Science  03 May 2019:
Vol. 364, Issue 6439, pp. 464-468
DOI: 10.1126/science.aaw7354

Slip outraces the fluids

Earthquakes generated from fluid injections sometimes appear to occur much faster than would be expected from the speed of fluid migration. Bhattacharya and Viesca used observations from a controlled fluid injection to develop a model to understand the behavior of fluids and aseismic slip. They found that the aseismic slip outpaced the fluid migration, allowing for earthquake triggering to happen much faster and farther away from the initial fluid injection.

Science, this issue p. 464


Earthquake swarms attributed to subsurface fluid injection are usually assumed to occur on faults destabilized by increased pore-fluid pressures. However, fluid injection could also activate aseismic slip, which might outpace pore-fluid migration and transmit earthquake-triggering stress changes beyond the fluid-pressurized region. We tested this theoretical prediction against data derived from fluid-injection experiments that activated and measured slow, aseismic slip on preexisting, shallow faults. We found that the pore pressure and slip history imply a fault whose strength is the product of a slip-weakening friction coefficient and the local effective normal stress. Using a coupled shear-rupture model, we derived constraints on the hydromechanical parameters of the actively deforming fault. The inferred aseismic rupture front propagates faster and to larger distances than the diffusion of pressurized pore fluid.

Fluid injections related to geothermal or oil and gas operations are associated with enhanced local seismicity rates (18). These earthquakes often occur as an enlarging cloud of seismicity migrating away from the point of injection (9, 10). The cloud’s diffusive growth is usually interpreted as representing the diffusion of increased pore-fluid pressure reaching prestressed faults and triggering instabilities (9, 11). However, a growing body of literature suggests that a primary response of faults to fluid pressurization could also be slow, aseismic slip, which, in turn, could transmit the solid stresses that trigger the observed induced seismicity (1216). Large, spatially extended regions of fluid-activated aseismic slip have been inferred in response to deep fluid injection (at depths greater than 1 km) (13, 16). Such large regions of aseismic slip, in some cases, are hypothesized to have triggered moderate-sized earthquakes on faults unlikely to be pressurized by means of fluid migration (16, 17). Yet, whether realistic subsurface hydromechanical conditions permit fluid-activated aseismic slip to transfer potentially earthquake-inducing stresses to regions remote from injection activity has never been tested with data-constrained models of fault slip (18).

Models predict that aseismic slip on fluid-pressurized faults can outpace pore-fluid migration under certain conditions (19). Specifically, theoretical studies of injection-induced earthquake nucleation indicate that the preinjection stress state may influence the aseismic propagation of a rupture front relative to fluid migration (19). However, in the absence of continuous, in situ measurements of active fault deformation and hydrological processes, the existing observations of the spatiotemporal evolution of fluid-activated aseismic slip are only indirect and rely primarily on seismicity migration or surface geodetic observations (16, 17). This leaves the details of both pore-fluid migration and aseismic slip propagation, and their interrelationship, ill-constrained. As a result, these model-predicted scenarios and their underlying physics have had limited opportunities to be satisfactorily tested against direct observations of fault hydromechanical behavior.

Fluid-injection experiments on shallow crustal faults bridge this data gap by activating and directly measuring fluid-induced aseismic slip and pore-pressure evolution at the fault (15). However, we do not know whether observations from a single borehole can constrain models of spatially extended shear ruptures stimulated by pore-fluid migration. We use the experimental data from Guglielimi et al. (15) within a probabilistic inversion framework to estimate the hydromechanical parameters of an activated fault and spatiotemporal evolution of slip. The joint measurements of both pore pressure and slip impose unambiguous constraints on the relationship between pore-fluid migration and aseismic-slip propagation, revealing that the inferred rupture ultimately outpaces pore-fluid migration.

The fluid-injection experiments of Guglielmi et al. (15) activated carbonate faults at depths around 280 m in the Low-Noise Underground Laboratory in southeastern France. They injected water with an increasing rate over a period of 1400 s into the target fault zone from a cross-cutting vertical borehole. They recorded pore pressure, flow rate, shear, and normal displacement of the fault at the borehole during injection. In response to fluid injection, the fault accommodated around 1.2 mm of total slip at the borehole. The first ~0.3 mm of slip was seismically quiescent, but further slip accumulation was accompanied by seismicity recorded at down-hole seismic sensors in nearby monitoring wells. Gugliemi et al. inferred the slip episode was primarily aseismic because (i) the total accumulated seismic moment was inferred to be a small fraction of the moment estimated from the observed slip and approximated rupture area, and (ii) both the slip and pressure time histories appeared insensitive to the onset of seismicity (15).

Our physical model for this experiment has two coupled components: (i) a hydrological model for fluid flow into the rock formation (Fig. 1A) and (ii) a model for quasi-static rupture of a planar fault (Fig. 1C). We presumed that fluid flow was fault-parallel and radially outward from the borehole, on the basis that the slip surface was embedded within a narrow fault damage zone, which provided a compliant conduit with high fluid permeability relative to the surrounding host rock (Fig. 1A) (20). The hydrological model has two parameters: permeability k (m2) and hydraulic storage Ss (m−1). The model response to borehole injection is the spatiotemporal evolution of pore-fluid pressure p within the damage zone, which we solved numerically subject to the variable injection rate of the field experiment. We used an adaptive Markov chain Monte Carlo (MCMC) routine to obtain the full posterior probability distributions for permeability and specific storage (20) given the pore-pressure data of (12) (Fig. 1B). The hydrological model provides a spatiotemporal evolution of pore-fluid pressure away from the borehole and along the fault plane. Our model for fault slip is that of a shear rupture driven by the evolving fluid pressure and concomitant local decrease in frictional strength (Fig. 1C). We account for the decrease by assuming that shear strength is the product of the local effective normal stress, the difference between the fault normal stress σ and the local pore-fluid pressure p, and a coefficient of friction (17). Our rupture model parameters included the magnitude of the initial fault shear stress τ, initial normal stress σ, and the host-rock shear modulus Go. We parametrized the friction coefficient by hypothesizing it to either have a constant value f or to linearly weaken with slip from a peak value of fp to a residual value fr over a slip scale Dc. The slip at the borehole that we predicted with this model was compared to the observed slip (Fig. 1D) to obtain Bayesian constraints on the model parameters using the same MCMC routine applied to the hydrological data.

Fig. 1 Model schematics and in situ observations of pore pressure and slip.

Data are from (15). (A) Fluid is injected at an instantaneous volume rate Q(t) through the open section of a borehole into a permeable fault damage zone of width b, permeability k, elastic modulus G, and hydraulic storage Ss. The host rock is characterized by permeability ko and shear modulus Go. The map shows France, with LSBB indicating the location of the Low-Noise Underground Laboratory. (B) The volume rate injection history (gray) and the observed pressure response at the borehole (black). l/min, liters per minute. (C) Shear slip on the planar fault is driven by the background shear stress τ and facilitated by the decrease in effective normal stress σ¯=σ p(r,t) as pore pressure p(r, t) builds up, where r is the radial distance from the borehole center and t is the time from the start of injection σ is the fault-resolved normal stress. The rupture has an instantaneous slip parallel length L(t). (D) The observed fault slip (black) and fault opening (orange) at the borehole. The cumulative count of accompanying microseismicity is shown as red bars.

We compared our hydrological model with the data and inferred that increased permeability was likely to have accompanied fault slip, and we constrained the magnitude of the increase. We found that in response to the volume injection rate history, the pore pressure at the well (Fig. 2A) shows multiple instances where a nondecreasing injection rate accompanied rapid pressure drops [most prominently at 618 and 940 s; fig. S2, see (20)]. Our hydrological model, with constant permeability and specific storage, cannot accommodate these observations. As a first-order modification to accommodate this behavior, we allowed the model permeability to undergo step changes at the times corresponding to these anomalous pressure drops. The pressure history of our enhanced model captured the observed pressure history with plausible parameter values (Fig. 2A). The permeabilities derived from this four-parameter model (Fig. 2) required a 60% permeability enhancement across these pressure drops. The permeability magnitudes we determined are consistent with estimates from previous pulse injection tests carried out in the test area (~10−12 m2) (21). The magnitude range we determined for specific storage (10−3 to 10−5 m−1) also overlaps with the range of in situ estimates from fault damage zones (10−4 to 10−5 m−1) (22, 23). We chose the minimum root-mean-square misfit model from the posterior probability distributions of the parameters (Fig. 2B) as the appropriate description of the spatiotemporal pore-pressure evolution that stimulates the shear rupture. This particular model of ours has a hydraulic diffusivity of 0.04 m2 s−1.

Fig. 2

Model comparisons with observed pore-pressure history at the well in response to the imposed injection history. (A) The best fit to the pore-pressure time history, with two step changes in permeability (the four-parameter hydrologic model) at 618 and 940 s, is shown in red. The observed pore-pressure history (black curve) and imposed injection history (gray curve) are shown. (B) Posteriors for permeability are shown in blue for k1, orange for k2, and green for k3. (C) Posterior for the storage coefficient. The vertical dashed lines in (B) and (C) show the parameter values for the best-fit model shown in (A). The values are k1= 0.8 × 10−12 m2, k2= 1.1 × 10−12 m2, k3= 1.3 × 10−12 m2, and Ss= 2.15 × 10−4 m−1.

The inferred fault permeabilities show a monotonic increase (k3 > k2 > k1) without any a priori constraints on the parameter interrelation. Fault-zone permeability varies with shear slip on reactivated faults both in the laboratory (2426) and in the field (15, 27). Our estimates of the step changes in permeability indicate, however, that permeability enhancement saturates even as slip continues to accumulate at an accelerated rate. The posteriors for k3 and k2 overlap to a greater extent than those for k1 and k2 in Fig. 2B. This saturation in permeability enhancement, instead, seems to more closely follow the observed leveling of the fault normal displacement after about 1000 s of injection reported in (15) (Fig. 1D).

Examining the history of slip at the borehole provides evidence for a slip-weakening description of fault strength, with constraints on the change in friction coefficient and characteristic slip distance, as well as the host rock shear modulus and in situ stress state. We solve for the evolution of slip at the borehole using our best-fit hydrological model. We found that the slip acceleration cannot be captured with the slip predicted by a constant friction coefficient model [fig. S4, see (20)]. A potential source of slip acceleration is the weakening of the fault that accompanies sliding. Laboratory tests of rocks under triaxial compression suggest that the postfailure strength of initially intact carbonate and granitic rocks reduces approximately linearly with slip (2830). We choose such piecewise-linear slip weakening as a simple description of reduction in fault strength. With a slip-weakening friction coefficient, our modeled quasi-static slip at the borehole captured the prominent features of the slip history. Our best-fit model (pink curve in Fig. 3A) reproduced the initial slow accumulation of slip transitioning into the subsequent rapid acceleration beyond 1180 s of injection. However, this fit does not explain the initial, slower acceleration before 1180 s. An excellent fit to the first 1180 s of slip alone (yellow curve in Fig. 3A) requires a Dc approximately three times as small (0.10 as opposed to 0.37 mm, Fig. 3F) and a Go around twice as large (20.71 as opposed to 11.84 GPa, Fig. 3C).

Fig. 3

Fits to the observed slip at the borehole in response to fluid injection. Data are from (15). (A) The model time history (purple and orange curves) reflects the slip at the center of nearly elliptical ruptures whose strength is determined by the local effective normal stress and a piecewise-linear slip-weakening friction coefficient with peak strength fp, residual-to-peak strength ratio fr/fp and slip-weakening distance Dc. The shaded regions around the colored lines represent 10% of the accepted samples, with shading density representing frequency of occurrence, and the blue squares are the observed slip. The purple curve fits the entire time history of slip and the orange curve the first 1180 s only. The three contour insets show the spatial distributions of slip corresponding to the purple curve at the indicated times. The dashed blue contours mark the locations where increases in pore pressure are 2, 1, and 0.5 MPa in outward order. (B) The linear slip-weakening model derived from the fits in (A). (C to G) The Bayesian posteriors corresponding to the fits in (A). The reported values of Go and τ are derived assuming fp = 0.6. The vertical dashed lines represent the least misfit models shown in (A).

Our data-inferred model has a rupture front that accelerates ahead of a region of substantial pore-fluid pressure increase. The three insets in Fig. 3A show the spatial extent of the slip-weakening rupture at three instances of time. The top-right inset shows that the slip-parallel extent of the aseismic rupture near the end of fluid injection is about three times the size of the region which has been pressurized by 1 MPa or more. The advancement of aseismic slip ahead of fluid migration is due to the large inferred background stress relative to peak strength (τ/fpσ ~ 0.64) and the large overpressures sustained over the last 600 s of injection. The inferred rupture properties (τ > frσ) suggest that rupture propagation might have ultimately turned dynamic if fluid pressurization was continued and no other mechanisms were initiated (18, 19, 31). Such runaway ruptures might greatly exceed any estimate of the maximum induced earthquake magnitude derived assuming that the stimulated rupture is spatially limited by the extent of the fluid-pressurized volume (31, 32). On the other hand, if nucleation is otherwise prohibited, such large aseismic ruptures could continue to transmit stresses to regions outside the fluid-pressurized regions. We note that the conditions that allow such rapid growth of aseismic slip are not dependent on the absolute stress levels and, therefore, could reasonably be expected to apply to depths greater than that of the Guglielmi et al. (15) experiment.

Aseismic slip on fluid-activated faults can outpace pore-pressure diffusion under fairly general frictional strength prescriptions, provided that the fault is nearly critically stressed where the fault-resolved shear stress approaches the preinjection fault strength. We demonstrate the underlying mechanical principle with a simple, constant-friction model, wherein the fluid-driven rupture front also enlarges diffusively. We considered the in- or antiplane slip of a fault with a uniform preinjection state of stress (Fig. 4A), with strength determined by a friction coefficient f (20). In this scenario, a source of fluid that diffuses along the fault with a hydraulic diffusivity α can drive an aseismic rupture front that also propagates in a diffusive manner. Here, the rupture length grows proportionally to the diffusive length scale for pore-fluid migration,αt. The constant of proportionality, λ, depends on the preinjection state of stress and strength of the fault (Fig. 4B). Under conditions approaching a critical state of stress, the rupture front will quickly advance well beyond a region of substantial pore-pressure increase with λ1. Observations from studies of fractures cross-cutting boreholes indicate that the subset of fractures that are hydraulically conductive are also likely to be critically stressed (3335). This observation, combined with our model above, suggests that large overpressures created through subsurface fluid injection may activate aseismic slip over regions considerably larger than the pressurized region and transmit solid stresses to regions remote from the injection point. Accordingly, the observed diffusion of induced microseismicity away from injection sites may reflect stress perturbations due to aseismic slip of a fluid-conducive fault fracture network, rather than the diffusion of increased pore-fluid pressures (9, 11).

Fig. 4 Elementary model of aseismic slip front outpacing fluid migration.

(A) Model schematic of an enlarging sliding region of a fault (red) in response to a fluid source within a permeable fault zone. a(t), rupture half-length; x, along-fault distance from the injection point. (B) The red curve shows aseismic slip front amplification as a function of the fault stress parameter, assuming a constant friction coefficient f and a fluid source with constant pore-pressure increase Δp. σ¯ is the preinjection effective normal stress. The blue dashed lines show asymptotic behavior under end-member fault stress regimes. The front of the sliding region leads fluid diffusion (λ > 1) if the fault is critically stressed (20).

Our study also provides valuable in situ constraints on the physics underlying rupture models in general. Such models of spatially extended fault slip almost exclusively use either remote seismological or laboratory-derived estimates for fault constitutive parameters, in particular for the determination of fault strength. These parameters critically influence fault stability and its hazard potential. Some of these parameters, like the slip weakening distance Dc, show orders of magnitude of variation between field-scale seismological [Dc ≈ 0.5 to 1 m, (36, 37)] and centimeter-scale laboratory measurements [Dc ≈ 10−6 to 10−3 m, (38, 39)]. Our estimates of Dc, and other fault constitutive parameters, show instead that parameter values derived from in situ measurements at the field scale broadly agree with those derived in the laboratory (40). This provides the first proof of concept that such well-instrumented “natural-fault laboratories” could help constrain the broader issue of scale dependence of fault constitutive behavior and, as a result, increase our confidence in the realism of physics-based models of earthquake hazard.

Supplementary Materials

Materials and Methods

Supplementary Text

Figs. S1 to S4

References (4258)

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: The authors thank F. Cappa for providing the high-resolution pore-pressure time history from the field experiments and details of the experimental setup. Funding: This work was supported by U.S. Geological Survey (USGS) grant G17AP00016, National Science Foundation (NSF) grant EAR-1653382, and the Southern California Earthquake Center (SCEC). SCEC is funded by NSF Cooperative Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. The experimental data were obtained with funding by the Agence Nationale de la Recherche (ANR) through the HPPP-CO2 project under contract ANR-07-PCO2-0001 and the HYDROSEIS project under contract ANR-13-JS-06-0004-01. Author contributions: P.B. contributed to methodology, formal analysis of data, and writing. R.C.V. contributed to conceptualization, methodology, and writing. Competing interests: The authors declare no competing interests. Data and materials availability: All data used in this work has been made available in the supplementary materials of (15). The Fortran library developed for Bayesian inversion is freely available at Zenodo (41).
View Abstract

Stay Connected to Science

Navigate This Article