This paper presents a timely and detailed study of significant injection-induced seismicity recently observed in the Sichuan Basin, China, where shale-gas hydraulic fracturing has been initiated and the aggressive production of shale gas is planned for the coming years. Multiple lines of evidence, including an epidemic-type aftershock sequence model, relocated hypocenters, the mechanisms of 13 large events (M W > 3.5), and numerically calculated Coulomb failure stress results, convincingly suggest that a series of earthquakes with moment magnitudes up to M W 4.7 has been induced by “short-term” (several months at a single well pad) injections for hydraulic fracturing at depths of 2.3 to 3 km. This, in turn, supports the hypothesis that they represent examples of injection-induced fault reactivation.
Xinglin Lei1, Dongjian Huang2, Jinrong Su3, Guomao Jiang2, Xiaolong Wang4, Hui Wang2, Xin Guo4 & Hong Fu5
The geologic reasons why earthquake magnitudes associated with hydraulic fracturing operations are so high in this area are discussed. Because hydraulic fracturing operations are on the rise in the Sichuan Basin, it would be beneficial for the geoscience, gas operator, regulator, and academic communities to work collectively to elucidate the local factors governing the high level of injection-induced seismicity, with the ultimate goal of ensuring that shale gas fracking can be carried out effectively and safely.
When any type of fluid is pressure-injected into an underground reservoir, as is done during fluid waste disposal and shale gas hydraulic fracturing (also referred to as hydro-fracturing or fracking), the pressure of the fluids (pore pressure) underground increases, and the underground stress distribution may change. Stress redistribution and pore pressure changes within and surrounding the reservoir may lead to geomechanical changes, fault reactivation, micro-seismicity, and even damaging earthquakes1,2,3,4,5,6,7. The abnormal increase in seismicity in areas such as the central United States8, 9 and western Canada10, 11 is considered to be a direct result of the rapid increase in injections for several applications including wastewater disposal, enhanced oil recovery (EOR), and shale gas hydraulic fracturing. Such injection-induced seismicity appears to be strongly localized, and moment magnitude (M W ) ≥ 3 events have only been directly tied with wells at a very limited number of sites3, 12, 13.
In recent years, shale gas production has grown rapidly because of the use of horizontal drilling and advancements in hydraulic fracturing technology. Simultaneously, earthquakes caused by the injections used for hydro-fracturing have become a point of significant worldwide public interest14, and a recent study has shown that “short-term” injections for hydraulic fracturing have probably induced earthquakes up to M W ≥ 3 at some sites13.
Shallow earthquakes of a moderate size (M W 3–5) can be destructive3. In addition, the reactivation of pre-existing faults may hinder hydraulic fracturing and lead to the migration of injected fluid into shallower or deeper formations. Therefore, because the technological developments that address such public concerns are essential to the safe conduct of hydraulic fracturing, well operators must develop or be provided with tools that will allow them to manage or forecast the occurrence of destructive earthquakes, and to avoid the reactivation of pre-existing faults, especially those that connect shallower formations.
Herein, we focus on shale gas site blocks in the town of Shangluo and its environs in the Sichuan Basin, China (indicated as ‘A’ in Fig. 1 and Supplementary Fig. 1), referred to hereafter as “the study area”. So far, five M W > 4.0 earthquakes have occurred concurrently with fracturing operations. The largest of these (M W 4.7), which occurred on 28 Jan. 2017, caused significant damages to nearby farmhouses and other constructions (Supplementary Fig. 21). This event was probably the largest earthquake induced by fracking injections reported so far. The remainder of the present paper is organized as follows.
Figure 1 Seismicity of the Sichuan Basin and surrounding regions. The gray and colored dots denote seismic events with magnitudes of at least M L 2.5 and M L 1.5 that have occurred since 1980 and 2014, respectively. Fault data are imported from the digital “Map of active tectonics in China”52. The basement faults are hidden deep faults in the crystal basement and show minor activity53. The oil and gas field outlines are digitized and modified from an open report54. The Shangluo (‘A’, study area of this paper), Weiyuan (‘B’), and other major shale gas blocks in which hydraulic fracturing is in operation are also shown. Area ‘C’ covers the region of wastewater injection induced seismicity in gas-depleted reservoirs3, 20, while seismicity in area ‘D’ was induced by water injection in deep salt wells. The background topography is based on the Shuttle Radar Topography Mission (SRTM3) digital elevation model (DEM) data (http://srtm.csi.cgiar.org/SRTM3). The map was created using the free software GeoTaos_map (developed by Xinglin Lei; https://staff.aist.go.jp/xinglin-lei/) and finished with the software CorelDRAW × 8. (Copyright (c) 2016 Corel Corporation. All rights reserved.)
First, we briefly describe the seismicity and hydraulic fracturing operations in the study area. Then, we present the results of this study, including statistical analyses, moment tensor inversion of the largest earthquakes (MW > 3.5), hypocenter relocation, numerically calculated Coulomb Failure Stress results, and a seismogenic index. These are followed by a discussion and conclusion section. Finally, the Methods section is presented to explain the summarizing methods used in this study.
Seismicity and hydraulic fracturing operation in Shangluo shale gas site
In the Shangluo shale gas site (indicated as ‘A’ in Fig. 1 and Supplementary Fig. 1), the target reservoir for hydraulic fracturing operations is a Silurian mudstone/shale formation that has a burial depth of approximately 2.3 to 3 km. This site corresponds to a wide and flat syncline. In the period from 1970 to Oct. 2008, only 60 earthquakes with M L greater than 2.5 were observed within the study area (Supplementary Fig. 4), indicating that the area has a low level of background seismicity. However, between the start of shale gas prospecting in 2008 and 2012, an increasing event rate was observed (Supplementary Fig. 4a). Based on our survey, we believed that the limited well injections conducted for evaluation purposes might be responsible for the increase (Supplementary Fig. 4).
In the study area, well pads intended for hydraulic fracturing generally have four to eight wellbores with horizontal lengths up to approximately 2500 m (Supplementary Fig. 3). The interval between the lateral portions of two neighbor wells is generally 300 to 400 m16. The lateral portions are commonly drilled in the direction parallel to hydraulic fracture opening (perpendicular to the planes of the fractures), which is perpendicular to the maximum horizontal principal stress, in order to maximize the volume stimulated by the induced fractures. In our study area, borehole breakout data consistently indicate that the maximum horizontal stress axis trends approximately from N100°E to N115°E17, 18 (Supplementary Fig. 3). A multistage hydraulic fracturing technique was applied for treatment. Both single- and multi-well completions are used16.
Figure 2 Epidemic-type aftershock sequence (ETAS) model results for earthquakes occurred in the Shangluo shale gas block (Fig. 1A) from 2014 to 20 Jan. 2017. (a) Mechanism solutions of the largest earthquakes (indicated by circled number and arrows) estimated by the gCAP method using two velocity models. (b) M-t and N-t plots. (c) ETAS parameters and their standard errors (number in brackets) (see Method section and Supplementary Fig. 5 for details). The time-varying forcing rate is plotted with respect to time. The red solid and dashed lines show major hydraulic fracturing time windows, and N# and Y# indicate the associated well pad numbers in the Ning-201 and Ys108 shale blocks, respectively.
Coinciding with the start of systematic horizontal well shale gas hydraulic fracturing, the observed earthquake rate within the study area increased dramatically (Fig. 2b). During the period from Dec. 2014 to 20, Feb. 2017, more than 2,400 M L ≥ 1.0 events were observed there, including four M W ≥ 4.0 events. The largest M W 4.7 event heavily damaged structures in the nearest villages (23 houses collapsed and 548 houses were heavily damaged) (Supplementary Fig. 21). Three factors might be responsible for making this event particularly destructive: (1) a very shallow focal depth of only 1.8 km (Supplementary Fig. 14); (2) a reverse mechanism (Supplementary Fig. 14); and (3) the fact that two villages were located within 2 km of the epicenter (Supplementary Fig. 21.
It should be noted that another large event occurred on 4 May 2017 having a M W of 4.6, focal depth of 2.5 km, and strike-slip mechanism solution (Supplementary Fig. 17), but was not so destructive.
Comparisons of magnitude versus time (M-t) plots for earthquakes and major hydraulic fracturing operations clearly show that significant fluctuations in the seismic event rate correspond with the timing of nearby hydraulic fracturing operations (Supplementary Fig. 4a,b). The frequency-magnitude distribution shows a magnitude of completeness of 1.25 and a b-value of approximately 0.9 (Supplementary Fig. 4c). Note that a b-value close to 1.0 is a sign of fault activation (injection induced or not), but much higher microseismicity values directly caused by fractures opening during hydraulic fracturing operations are commonly observed19.
Characteristics of fracking-induced earthquakes
To examine the statistical features of the observed seismicity, we applied our seismic dataset to the epidemic-type aftershock sequence (ETAS) model, which is useful for statistically separating the total seismicity into external forced background activity and Omori-type aftershocks3, 20,21,22. When used with a time-varying forcing rate, the ETAS model results show that the seismicity is governed by an external force (here, injection-related stress) with a total forcing rate of up to 87% and that Omori-type aftershocks are rare, encompassing only 13% of all earthquakes (Fig. 2). These results show similarities with recorded seismicity induced by injection of wastewater into depleted conventional gas reservoirs at other sites3, 20. In fact, the largest earthquakes in the investigated series did not produce significant aftershocks, as could be expected for tectonic earthquakes. In general, the seismicity rate decreased to a much lower level quickly after hydraulic fracturing ceased. However, in some cases, persistent seismicity continued for several tens of days after completion (Fig. 2). Mechanism solutions of the 13 largest earthquakes (M W > 3.5), based on moment tensor inversion using the generalized cut and paste method (gCAP)23, showed very consistent results (Table 1; Figs 2a and 3), with the best focal depth, corresponding to the minimum misfit error, falling in the range from 1.8 to 4 km (Table 1). Most events showed nearly pure double-couple (DC) mechanisms with very limited non-DC components (less than 4%), and only two events had a negative isotropic (ISO) component greater than 10%, probably due to uncertainty (Table 1).
Table 1. Mechanism solution of the largest earthquakes (MW> 3.5). Columns M1 and FM2 show strike/dip/rake of the two nodal planes. Giso, Gvlcd are respectively the squared ratios of the scalar potency of the ISO and CLVD component to the total scalar potency, representing their relative strengths23. H and MW are central moment depth and moment magnitude of the best solution, respectively. For comparison, ML from the catalog and focal depth (H*) and its standard error determined by the double differential relocation method are also shown. The focal depths are related to a mean surface level. A preferred velocity model (Supplementary Fig. 6M-2) based on results of seismic ambient noise tomography of the Sichuan Basin was used for moment tensor inversion. From this table, we obtained an empirical relation between ML of the catalog and MW of the inversion: MW=ML− 0.22 (R2=0.85).
Two events showed strike-slip motion dominated mechanisms, three events showed strike-slip motion with a significant reverse component, and five events showed pure reverse faulting motion. Other earthquakes produced a reverse faulting mechanism with a lesser strike-slip component. The azimuth of the P-axes of all solutions, which were concentrated in a narrow azimuth range from 90° to 130° (Fig. 2a), agreed very well with the regional maximum horizontal stress (σH ) direction, which trends approximately from N100°E to N115°E as was previously mentioned. The fact that both strike-slip and reverse faulting earthquake were observed suggests that the regional minimum horizontal stress (σh) is close to the vertical stress σV .
The routinely determined earthquake hypocenters (reported by the Yibin Earthquake Mitigation Administration (YEB) and compiled in the China Earthquake Data Center (CEDC) catalog) are scattered around the well pads associated with the hydraulic fracturing operations (Supplementary Fig. 1b). Using phase data manually picked and compiled by the YEB and the Hypocenter Double-Difference (HypoDD) program24, we relocated a small number of earthquake hypocenters (208 of 1,045 events) for the period from 2014 through 2015. The resulting hypocenter locations showed definite improvements but remained imprecise, with estimated errors of up to a few kilometers (Fig. 3a).
Figure 3 Map and section views of relocated earthquakes for two periods: (a) Jan. 1, 2014 through 2015, and (b) since Jan. 1, 2016, with the focal mechanisms of the largest earthquakes (M w ≥ 3.5). The earthquake symbols are colored by date and scaled for magnitude. The gray cross lines show vertical and the horizontal errors given by the relocation program. The red arrows and numbers on the faults indicate the dip direction and dip angle of the fault. The big inward-pointed arrow shows the orientation of the maximum horizontal principal stress (SHmax). This map was created using the free software GeoTaos_map (developed by Xinglin Lei; https://staff.aist.go.jp/xinglin-lei/) and finished with the software CorelDRAW × 8. (Copyright (c) 2016 Corel Corporation. All rights reserved.)