Seismic Slip-Pulse Experiments Simulate Induced Earthquake Rupture in the Groningen Gas Field.

Rock materials show dramatic dynamic weakening in large-displacement (m), high-velocity (∼1 m/s) friction experiments, providing a mechanism for the generation of large, natural earthquakes. However, whether such weakening occurs during induced M3-4 earthquakes (dm displacements) is unknown. We performed rotary-shear experiments on simulated fault gouges prepared from the source-, reservoir- and caprock formations present in the seismogenic Groningen gas field (Netherlands). Water-saturated gouges were subjected to a slip pulse reaching a peak circumferential velocity of 1.2-1.7 m/s and total displacements of 13-20 cm, at 2.5-20 MPa normal stress. The results show 22%-81% dynamic weakening within 5-12 cm of slip, depending on normal stress and gouge composition. At 20 MPa normal stress, dynamic weakening from peak friction coefficients of 0.4-0.9 to 0.19-0.27 was observed, probably through thermal pressurization. We infer that similar effects play a key role during induced seismic slip on faults in the Groningen and other reservoir systems.

Introduction

Induced seismicity related to hydrocarbon production or fluid injection in subsurface reservoirs is increasingly widespread (Elsworth et al., 2016; Ellsworth et al., 2019; McGarr et al., 2002; Nicol et al., 2011; Segall, 1989). In the Netherlands' vast Groningen gas field, induced events have occurred since the 1990s (van Thienen‐Visser & Breunese, 2015), reaching magnitude M  = 3.6 (Figure 1). They are the result of reactivation of pre‐existing normal faults that crosscut the reservoir (Dost & Haak, 2007; Van Eijs et al., 2006), and occur within or close to the reservoir interval (Smith, 2019; Spetzler & Dost, 2017)—Figure 1c. A thorough understanding of the mechanical properties of these faults is needed to assess the potential for induced fault reactivation and associated seismic hazard. Specifically, the data on the dynamic frictional behavior of the relevant fault rocks, obtained at realistic imposed displacement and velocity trajectories, are needed to constrain geomechanical modeling of fault rupture and associated seismicity (e.g., Buijze et al., 2017, 2019; Rutqvist et al., 2016; Zbinden et al., 2017).

Figure 1

(a) Map showing the Groningen gas field and hypocenter locations of M w > 1 earthquakes (up to June 2020, based on the public Groningen earthquake catalog from www.rdsa.knmi.com), major faults (gray lines) and the locations of the Stedum‐1 (SDM‐1) and Zeerijp‐2 (ZRP‐2) wells. (b) The principal lithologic units of the Groningen reservoir (depths from SDM‐1). The composition of the simulated fault gouges used here is indicated in wt%—see Hunfeld et al. (2017). (c) Depth distribution of earthquake hypocenters of 87 events (from February 2014 to July 2016), after Spetzler and Dost (2017).

In the context of natural earthquakes, high‐velocity (HV) friction experiments on both intact and granular rock materials commonly show a large reduction in friction coefficient at slip velocities above ∼0.1 m/s, for which numerous mechanisms have been proposed (see Di Toro et al., 2011; Niemeijer et al., 2012; Wibberley et al., 2008 for reviews). However, the vast majority of HV friction experiments performed to date have been conducted at low normal stress (<10 MPa) and dry or at room‐humidity. Moreover, with the aim of addressing major natural earthquake slip, they mostly feature large displacements (m's) imposed at constant slip velocities of 0.1–2 m/s, generating large amounts of frictional heat. Therefore, many observations reported in the HV friction literature may not apply to magnitude 3–4 induced earthquakes, where fault slip occurs in a short pulse featuring rapid acceleration and immediate deceleration, within displacements of a few decimeters (Heaton, 1990; Tomic et al., 2009; Zoback & Gorelick, 2012).

Only a few experimental studies have employed such a slip‐pulse (Fukuyama & Mizoguchi, 2010; Liao et al., 2014; Spagnuolo et al., 2015). These have demonstrated dynamic weakening in room‐dry, rock‐on‐rock friction tests over just a few cm displacement, at normal stresses ≤10 MPa. However, to what extent this occurs in gouge‐filled faults at upper crustal stress and pore fluid conditions, during induced earthquake slip such as occurring in the Groningen field, has not been investigated.

Here, we address this knowledge gap via a series of rotary‐shear experiments imposing a short‐displacement (13–20 cm), HV slip‐pulse simulating induced magnitude 3–4 seismic slip on confined faults containing water‐saturated gouges derived from the key lithologic units present in the Groningen gas field. We show that even at short displacements, rapid weakening occurs to varying degrees in all materials tested, probably via thermal pressurization of the pore fluid.

Experimental Methods

Sample Materials

The samples investigated here were derived from the Carboniferous shale/siltstone source‐rock, the Slochteren reservoir sandstone, the overlying Ten Boer claystone and the Basal Zechstein anhydrite‐carbonate caprock (Figure 1). Simulated fault gouges with a grain size <50 μm were prepared by crushing and sieving drill core and cuttings derived from these formations, sampled from the Stedum‐1 (SDM‐1) and Zeerijp‐2 (ZRP‐2) boreholes in the seismogenic center of the Groningen field (courtesy of the field operator, NAM). The composition of these gouges was analyzed using conventional X‐Ray Diffraction (XRD) methods (see Hunfeld et al., 2017)—Figure 1b.

High‐Velocity Slip‐Pulse Experiments

A total of 28 experiments were performed on water‐saturated gouges at normal stresses () of 2.5–20 MPa, using the HV rotary shear apparatus installed at the State Key Laboratory of Earthquake Dynamics in Beijing (see Ma et al., 2014). Two types of low thermal conductivity, Ti‐alloy piston assemblies were used, namely (a) a conventional solid‐cylinder assembly (40 mm diameter), and (b) an annular ring‐shear assembly (outer diameter 40 mm, inner diameter 28 mm). In both assemblies, a ∼1 mm thick layer of water saturated gouge paste was sandwiched between the roughened faces (#80 SiC) of the stationary lower piston and the rotating upper piston. An outer Teflon sleeve, and, in the case of the ring assembly an additional inner Teflon cylinder, were used to minimize gouge extrusion during testing. To improve gouge confinement, a small, water‐filled pressure buffer was used to apply a confining pressure of 2 MPa to the outer Teflon sleeve. Temperature and pore fluid pressure were measured using thermocouples and pore pressure transducers embedded in the lower (stationary) piston at the gouge‐piston interface. Further details and a table listing all experiments can be found in the Supporting Information S1

Each sample was pre‐sheared (see Supporting Information S2) through a displacement of 0.45 m at a velocity of 1.5 mm/s (measured at the outer edge of the setup) under a normal stress of 2.5–20 MPa. After pre‐shear, a HV slip pulse was imposed at the same normal stress under undrained conditions, that is, with the closed‐in pore fluid initially at 1 atm. To achieve this, the drive system was programmed to apply a symmetrical, triangular velocity profile characterized by a peak circumferential velocity of 1.5 m/s and 0.15 m total displacement, occurring over ∼0.25 s. The peak velocity was reached within 0.14 ± 0.02 s. Because the drive motor supplied a fixed power input per slip pulse, the peak velocity varied between 1.26 and 1.71 m/s in the solid‐cylinder tests and 1.21–1.70 m/s in the ring experiments, decreasing systematically with increasing normal stress (Table S1). Total displacements attained were 0.131–0.192 m, also decreasing systematically with increasing normal stress.

The data on axial load, torque, and upper piston rotation (vs. machine frame), recorded at 2 kHz, were processed to obtain shear stress (MPa), normal stress (MPa), and apparent friction coefficient (−), versus circumferential displacement and time. No corrections for shear stiffness of the machine were applied in calculating displacement. Pore pressure measurements showed strong spatial and temporal variations during each experiment, precluding determination of effective normal stress and true friction coefficient. To obtain apparent friction coefficient versus displacement curves, the raw displacement (angular rotation) data were smoothed using a moving average with a window of 101 data points (corresponding to a time interval of 50.5 ms). Water‐charged calibration runs with no sample present showed that Teflon seal friction amounted to only ∼2–5% of the total torque measured in real experiments at 5 MPa normal stress. However, since seal friction likely varied from run to run, no seal friction correction was applied.

Results

Friction Coefficient Versus Displacement and Normal Stress

Representative examples of apparent friction coefficient versus displacement curves obtained in experiments at 10 MPa normal stress are shown in Figure 2 for both set‐ups. Similar results were obtained in experiments at different normal stresses. Apparent friction increases abruptly upon initial acceleration. In most cases, a distinct peak occurs within the first 0.01–0.03 m of slip, followed by significant weakening toward a minimum dynamic strength attained at displacements of 0.07–0.14 m. During deceleration, friction generally increases again until the motor stops and displacement ceases, at which point the shear stress supported falls rapidly to zero. The only exception was observed in Ten Boer claystone samples tested in the solid‐cylinder setup (Figure 2a). In these runs, after the initial increase upon slip initiation, the friction coefficient hardly changed during the entire slip pulse—at all normal stresses investigated (2.5, 5, and 10 MPa).

Figure 2

Apparent friction coefficient versus displacement curves obtained at 10 MPa normal stress in (a) the solid‐cylinder and (b) ring setups. Insets show velocity versus displacement profiles for the Slochteren sandstone experiments. Figure S6 shows the individual friction, velocity, and axial displacement data for the same runs.

Peak and minimum dynamic friction coefficients ( and ) obtained in all experiments are plotted versus normal stress in Figure 3. The Basal Zechstein samples showed the highest peak friction (0.76 <  < 0.97), but the lowest dynamic friction (0.17 <  < 0.31). Slochteren sandstone and Carboniferous shale/siltstone samples showed values in the range 0.50–0.76, and values of 0.27–0.46. All experiments performed using the solid‐cylinder and ring setups showed broadly similar behavior, with the exception of the Ten Boer claystone runs. These showed constant friction ( 0.38) in the case of solid‐cylinder experiments, whereas ring experiments showed a peak friction of 0.45 <  < 0.58 falling to a low dynamic minimum of 0.19 <  < 0.39. In the case of the Slochteren sandstone gouges, was insensitive to normal stress, whereas all other materials showed a slight increase with increasing . Minimum dynamic friction coefficients generally decreased with normal stress in all materials. Note that for the phyllosilicate‐rich Ten Boer and Carboniferous gouges, was consistently higher in the ring setup than in the solid‐cylinders. This is likely caused by gouge extrusion (jamming) effects observed at the inner seal in the ring experiments on these gouge lithologies (see Supporting Information S2–S4).

Figure 3

Manually picked peak and minimum dynamic apparent friction coefficients ( and , see legend for definition) versus normal stress for all experiments performed.

The weakening distance d (displacement from to ) lies in the ranges 0.08–0.14 m for the Slochteren, Ten Boer, and Carboniferous gouges, and 0.04–0.10 m in the Basal Zechstein gouges (Figure 2—see also Table S1). No discernible normal stress‐dependence of d was observed in any of the gouge materials tested.

Temperature and Pore Fluid Pressure Measurements

In all experiments, temperature increased rapidly during the HV slip pulse and continued to increase briefly after slip stopped, reaching a peak value ( up to 108°C, Table S1) a few seconds after termination of slip, before slowly decaying back to ambient levels. Pore fluid pressure showed a similarly sharp peak ( up to 1.41 MPa, Table S1) in the case of experiments performed on Slochteren sandstone (Figure 4a) and Basal Zechstein gouges. In the case of the Ten Boer (Figure 4b) and Carboniferous gouges, the pore fluid pressure increase was much more gradual, reaching a maximum (up to 0.28 MPa, Table S1) only tens of seconds to several minutes after the slip pulse. Both temperature and pore pressure showed strong variations with the radial position (solid‐cylinder setup) and azimuthal position (ring setup) for all materials.

Figure 4

Examples of temperature and pore fluid pressure evolution measured during and after HV slip in (a) ring‐setup run LHV1426 and (b) solid‐cylinder run LHV1208. The positions of the thermocouples (T) and pore fluid pressure transducers (P) in the lower (stationary) pistons are indicated, as well as the pore fluid ports (F). The colors of the curves correspond to the respective sensor color, with pore fluid pressure change shown in dashed lines and the temperature change in solid lines.

Assessment of the Role of Thermal Pressurization

The increase in temperature and pore fluid pressure measured in all experiments (Figure 4 and Table S1) indicates that thermal pressurization (TP, Wibberley & Shimamoto, 2005) was active to some degree in all cases, with the largest recorded (1.4 MPa—Slochteren sandstone gouge at 20 MPa normal stress) accounting for a ∼7% reduction in apparent friction. However, since slip localized near the upper piston‐sample interface, ∼0.8 mm above the sensors, in most experiments (see Supporting Information S3), the magnitude of TP is likely much larger than measured. To estimate how much larger, we modeled the temperature and pore pressure evolution within the gouges assuming a slip zone of 100 μm thickness at the top of the samples (as observed post‐test for Slochteren sandstone run LHV1426, see Figure S4). We used a finite element method with COMSOL Multiphysics software, incorporating the theoretical relations for coupled temperature and pore fluid pressure evolution derived for TP due to shear heating by Andrews (2002) and Rice (2006). The heat input was calculated from the measured friction and velocity profiles in experiments performed at 20 MPa normal stress, and the temperature and pore fluid pressure measurements were used to constrain the “thermal pressurization coefficient," (Rice, 2006) (see Supporting Information S5 for details of the model and input parameters). Dilation/compaction during the slip‐pulse is not incorporated in the modeling.

The simulations for the Basal Zechstein gouge (Figure 5, based on run LHV1442) showed a peak temperature of 286°C and (volume‐averaged) pore pressure within the slip zone reaching 8.4–18.4 MPa. Although pore fluid phase transitions were not incorporated in the model, a comparison with the phase diagram of water (cf. Chen et al., 2017) shows that minor pore fluid vaporization would occur within the slip zone where P is highest (see Figure 5c), at least in simulations using a high fluid diffusivity ( = 0.1–0.5 mm2/s), corresponding to an equivalent permeability of the order of 10−19–10−18 m2, in combination with low or intermediate (see Table S3).

Figure 5

Thermal pressurization model results for Basal Zechstein run LHV1442 (ring setup, 20 MPa normal stress), using  = 0.06 MPa/°C and  = 0.1 mm2/s (see Supporting Information S5). (a) Schematic representation of the FEM configuration (shown in more detail in Figure S7). Slip is assumed to localize in a 100 μm thick zone at the top of the 0.8 mm gouge layer. (b) and (c) Temperature and pore fluid pressure distribution at t = 0.21 s after initiation of slip. Vaporization is expected to occur within ∼5% of the slip zone, as indicated. (d and e) Comparison between the modeled temperature/pore fluid pressure and the measurements (at sensors T1, T2, P1, P2—refer Figure 4a), showing the volume‐averaged results for the localized zone as well as for the point corresponding to the sensor positions (green dot in Figures 5b and 5c).

In the case of the Slochteren sandstone (LHV1426, Supporting Information S5), the simulations show similar patterns in temperature and pore pressure to those shown in Figure 5. However, simulations using the full heat input based on the measured friction and velocity profile show dramatic TP with T reaching 364°C and P reaching 22.1 MPa—for which extensive vaporization would occur throughout the slip zone. When using only 30% of the total heat input, which is qualitatively consistent with vaporization buffering the temperature rise (Acosta et al., 2018; Chen et al., 2017), the simulations match the measured T and P much better, resulting in a peak temperature of ∼124°C and pore pressure of 5.2–11.6 MPa within the slip zone.

The reduction in effective normal stress and hence apparent friction coefficient corresponding to the above peak P values would be 42%–92% in the Basal Zechstein and 26%–58% in the Slochteren sandstone gouges. Since ∼72% and ∼53% weakening was observed in runs LHV1442 and LHV1426 (refer Table S1), these results suggest that TP can account for 49%–100% of the dynamic weakening observed in the Basal Zechstein and sandstone gouges.

For the Ten Boer and Carboniferous gouges, the above modeling approach is prone to larger uncertainties due to (a) gouge extrusion in the ring setup, which affects the measured friction and thus the modeled heat generation (Supporting Information S3), and (b) limited accuracy in the determination of the thermal pressurization coefficient (Supporting Information S5). Therefore, no firm conclusions regarding the role of TP in the Ten Boer and Carboniferous gouges can be drawn.

Discussion and Implications

The heat generated during large, rapid earthquake slip (Kanamori & Brodsky, 2004) is widely believed to trigger thermally activated weakening mechanisms such as flash heating (FH) or thermal pressurization (TP), causing low dynamic friction. Theory and experiments indicate that FH is caused by sliding of highly stressed asperity contacts within rock‐on‐rock interfaces (Beeler et al., 2008; Goldsby & Tullis, 2011; Rice, 2006) leading to thermal weakening of the asperities, hence reduction in supported shear stress via dislocation mechanisms (Spagnuolo et al., 2015), dehydration (Brantut et al., 2011), decarbonation (De Paola et al., 2011), or melting (Hirose & Shimamoto, 2005). However, the effectiveness of flash heating in gouges, particularly in the presence of water or at small displacements, remains speculative, because (a) the shear velocity is distributed over many contacts in the localized zone of the gouge layer, and (b) the “quenching” effect of the pore fluid will buffer the temperature rise (Acosta et al., 2018; Chen et al., 2017; Yao, Ma, Niemeijer, et al., 2016; Yao, Ma, Platt, et al., 2016).

In fluid‐saturated faults, model calculations show that TP can strongly reduce the effective normal stress and thus the frictional resistance during large natural earthquakes (Lachenbruch, 1980; Rice, 2006; Sibson, 1973; Tanikawa & Shimamoto, 2009; Viesca & Garagash, 2015; Wibberley & Shimamoto, 2005). However, it is rarely observed directly in experiments (Badt et al., 2020; Faulkner et al., 2011; Rempe et al., 2020), due to the difficulty of reproducing dynamic rupture under realistic stress and pore fluid (pressure) conditions in the laboratory, let alone making measurements of local temperature and pore fluid pressure during slip.

The present experiments, supported by the modeling work, demonstrate that TP, and potentially even vaporization, likely played a significant role in our slip pulse experiments, at least for the quartz‐rich and anhydrite/carbonate‐rich gouges derived from the Slochteren reservoir and Basal Zechstein caprock formations of the Groningen gas field. However, whether other dynamic weakening mechanisms, such as FH, or silica gel and nano‐grain lubrication effects (Goldsby & Tullis, 2002; Han et al., 2010) were also active, cannot be assessed without extensive micro/nanostructural study and mechanistic modeling.

Whether TP was also a dominant weakening mechanism in the experiments on the clay‐rich Ten Boer and Carboniferous gouges remains speculative. In a qualitative sense, the lower hydraulic diffusivity obtained in measurements of bulk fluid transport properties of clay‐rich gouges, versus clay‐poor gouges (Bakker et al., 2016; Chen et al., 2013; Crawford et al., 2008), implies that the present clay‐rich gouges should be more susceptible to TP than the Slochteren and Basal Zechstein gouges. On the other hand, frictional heating in clay‐rich gouges is less efficient compared to clay‐poor gouges, due to a) the lower macroscopic shear stress supported by clay‐rich gouges as a result of a lower intrinsic friction coefficient, and b) the possibility of more distributed shear within a foliated matrix of aligned clay minerals (e.g., Haines et al., 2013; Rutter et al., 1986). The trade‐offs between these effects mean that it is difficult to assess whether TP would be more or less effective in the Ten Boer and Carboniferous gouges, versus the Slochteren and Basal Zechstein gouges. TP cannot be eliminated, however, and could still be the dominant factor controlling the modest dynamic weakening observed in the clay‐rich gouges.

Despite these uncertainties, our findings imply that TP is likely to be an effective weakening mechanism not only in large, natural earthquakes (Cubas et al., 2015; Noda & Lapusta, 2013), but also in relatively small, induced earthquakes such as in Groningen. However, even for short slip‐pulses, the effectiveness of TP depends on the fluid transport properties (fluid diffusivity, permeability—see Supporting Information), which may be different in the pre‐sheared gouges tested here versus natural faults (Faulkner et al., 2010). In addition, the larger‐scale fault structure (e.g., step‐overs, intersections), and off‐fault processes (damage zone, drainage), as well as shear‐induced dilatancy, are key factors that might limit the effectiveness of TP in natural faults. In the present experiments, the samples mostly showed a transient dilation in the order of 1%–3% of the total layer thickness (0.8 mm) during HV slip (see Figure S6). However, it cannot be assessed whether this dilation occurred in the shear band, bulk gouge, or both, and therefore its influence on the magnitude of the TP effect is uncertain (see Supporting Information S5). Therefore, care must be taken in extrapolating the present findings to either induced or natural M3‐4 earthquakes. Nonetheless, the rapid weakening observed in our slip pulse experiments offers a feasible mechanism for the larger induced events in Groningen, as well as for modest natural earthquakes, including fore‐ and aftershocks in tectonically active terrains.

Dynamic modeling studies of induced earthquakes in reservoir settings, including the Groningen field, typically assume a linear slip‐weakening model for individual faults featuring a uniform, fault‐wide friction drop from 0.6 to values in the range 0.2–0.5, to facilitate earthquake generation (e.g., Buijze et al., 2017, 2019; Rutqvist et al., 2016; Zbinden et al., 2017). However, there has been little physical basis to underpin the choice of dynamic friction values for the various lithologies involved. The dynamic weakening data and slip‐weakening trajectories obtained in the present experiments provide much‐needed experimental constraints for such dynamic rupture simulations and broadly support the assumptions made by Buijze et al. (2019) and Buijze (2020) in work on the Groningen field.

Supporting Information S1

Click here for additional data file.