Induced Earthquakes Before and After Cessation of Long‐Term Injections in Rongchang Gas Field

For more than three decades, the depleted Rongchang gas reservoir in China's Sichuan Basin was used for the disposal of unwanted water, which resulted in induced earthquakes, with magnitudes as high as 5.2. After all wells were closed, the frequency of seismic activity was observed to decay following a modified Omori law, and since April 2015, seismic activity again began to increase, and a M4.9 earthquake occurred on 27 December. The results of an epidemic‐type aftershock sequence (ETAS) model analysis show that forced seismicity accounted for more than 70% of the total events. For most M ≥ 3.5 earthquakes, including two M ≥ 4 events, the estimated overpressure was lower than the maximum injection pressure. These results, coupled with the fact that postinjection seismic activity has similar characteristics to seismicity during injection, indicate that the injected overpressure fluid was still the driving factor for postinjection seismic activity.


Introduction
Many studies have shown that changes in the subsurface stress field and the fluid pressure generated by injection have the ability to induce seismicity (e.g., Doglioni, 2018;Ellsworth et al., 2019;Foulger et al., 2018;Healy et al., 1968;Keranen & Weingarten, 2018;Lei et al., 2008;Zoback & Harjes, 1997). Research on seismicity induced by different types of injection in different areas has produced some interesting results. If an injection well is connected to a critically stressed fault (e.g., through a fracture), seismicity can be induced almost immediately after injection, as was the case at the Rangely oil field control test site in Colorado, where seismicity occurred less than 1 day after injection (Raleigh et al., 1976). If the raised pressure does not spread to the background value, the seismicity may last for several years after the injection ends. For example, at the Rocky Mountain Arsenal in Colorado, earthquakes continued for many years after injections stopped, due to the low-permeability rock layer boundary (Healy et al., 1968). The critical pressure threshold for initiating fault ruptures varies greatly from site to site. For example, in the Rangely area, the pressure threshold was shown to be approximately 26 MPa (Raleigh et al., 1976), while in the south Sichuan Basin, China, Mw ≥ 3.5 earthquakes have been induced by fluid overpressure ranging from 0.3 to 5.8 MPa (Lei et al., 2019a). In some cases, the maximum magnitude of the induced earthquakes associated with injection has shown progressive increases (Lei et al., 2008). In other cases, however, the largest event occurred in the very early stages of injection (Lei et al., 2013). In regions near injection wells and within the depleted reservoirs, induced seismicity has proven responsive to changes in injection parameters such as the pumping pressure and injection rate (Lei et al., 2008(Lei et al., , 2013. However, for seismicity far from an injection well, poroelastic stress change, in which a temporal change in the injection rate or pressure is damped, has been shown to be more responsive (Segall, 1989).
In recent years, an increase in the magnitude and number of induced events associated with fluid injection at different sites has led to a growing awareness of potential risks. Coupled with its socioeconomic significance, the expansion of fluid injection has led to a surge in scientific interest in induced seismic activity. The Sichuan Basin provided typical cases, in which both long-term and short-term injections have caused noticeable ("felt") induced seismicity, including some destructive earthquakes. Some large earthquakes, such as the M5.7 (in this article M indicates M L for events of M L < 5 and M S for events of Ms ≥ 5) earthquake on 18 December 2018 induced by shale gas hydraulic fracturing (HF) (Lei et al., 2019a), the M6.0 earthquake on 17 July 2019 (possibly induced by long-term injection for salt production) (Lei et al., 2019b), and several earthquakes of M4.3 to M5.4 induced by shale gas HF at the Weiyuan-Rongxian shale gas site (Lei et al., 2020), caused more than 200 casualties including more than 10 deaths and major economic losses. These cases challenge safety implementation measures of related industries and provide opportunities for research on induced earthquakes. Lessons learned from the Sichuan Basin are helpful for both seismology and industrial applications.
The Rongchang (RC) gas field in China's south Sichuan Basin is a well-known site of injection-induced earthquakes ( Figure 1). To date, more than 10,000 surface-recorded earthquakes, ranging up to M5.2, which were clearly injection induced, have occurred in the area (Lei et al., 2008;X. Wang et al., 2011X. Wang et al., , 2012Z. Wang et al., 2018). Curiously, seismicity has remained high even after all of the nearby injection wells were closed in 2013. Both M4.9 and M4.0 earthquakes occurred in 2016. These events injured two people and caused the collapse of 18 houses and extensive damage to more than 622 houses. These events were abnormally destructive, with a maximum hazard level of VI. The VI area extends for more than 10 km (Chongqing Earthquake Agency, 2016). Therefore, seismicity after closing of injection wells is also an important issue that must be addressed.
Following previous studies, the present study mainly considered seismicity in the area after 2008, especially during the period after all injections stopped. During the period from October 2008 to July 2011, a portable seismic network of six stations was installed, which afforded a detailed view. First, we briefly present the history of seismicity, injection, and the results of previous studies. Next, we introduce the data and data processing techniques used in the present study in section 3. Temporal variations in seismicity, statistics obtained using an epidemic-type aftershock sequence (ETAS) model, and stress patterns and overpressures estimated from focal mechanism solutions are presented in section 4. Finally, in section 5, in addition to providing the main conclusion of this study, we analyze the mechanism involved in postinjection seismicity. The results obtained in the present study are important for seismic hazard assessment for sites where long-term fluid injection has taken place.

History of Seismicity and Injection
In the RC gas field, the injection of unwanted water coproduced with natural gas from the RC and nearby gas fields on a routine basis began on 1 July 1978 and lasted for more than three decades. Unwanted water, collected from nearby production wells, was pumped into Permian formations at a depth of approximately 3 km through several deep wells and induced a large number of earthquakes of up to M5.2 (Lei et al., 2008) ( Figure 1b). Multiple lines of evidence, including spatiotemporal correlation between seismicity and injection, and parameters in statistical models, indicate that the seismic activity in the area was caused by industrial water injection. Lei et al. (2008) analyzed the detailed relationship between seismicity and injection rate during the period from 1970 to 2006, and Wei and Liu (2014) investigated the injection/seismicity connection through 2013 (Lei et al., 2008;H. Wei & Liu, 2014). After 2013, due to the inability to maintain a sufficient injection rate at the maximum allowable pressure, all injection wells were shut down, as determined by the operator based on safety and injection efficiency considerations. The maximum pressure was approximately 2.9 MPa at the earliest injection well Luo-4 and 6-9 MPa at other wells. Unfortunately, there are no detailed publicly available pressure data. However, we collected limited data from the pressure gauge installed at the well head, including some wellhead pressure data after injection was stopped (see section 5).

Data
Earthquake catalog data and phase data for relocation from 2008 to 2019 were downloaded from the China Earthquake Data Center (CEDC, http://data.earthquake.cn/). The magnitude of completeness in the earthquake catalog was 1.5 ( Figure 2b). Earthquake waveform data recorded by permanent regional stations were provided by the Data Management Centre of the China National Seismic Network at the Institute of Geophysics, China Earthquake Administration (registration required). In addition, continuous seismograms recorded by six temporary broadband stations operating from October 2008 to July 2011 (X. Wang et al., 2011) were used.

Moment Tensor Inversion
A moment tensor solution was calculated by applying the generalized cut and paste (gCAP) method (Zhu & Ben-Zion, 2013) using the body wave and surface wave information recorded by the permanent broadband stations within an epicenter distance of less than 300 km. The body wave was filtered using a 0.02-0.15 Hz  (Burchfiel et al., 1995;G. Wei et al., 2008;Xu et al., 2012). (b) Magnitude and accumulated number ofM ≥ 2.5 earthquakes that occurred from 1975 through 2019, overlaid with injection windows. The number ofM ≥ 1.5 earthquakes since 2008 is also plotted. Note that the number ofM ≥ 1.5 events declines greatly after injection was terminated but largerM ≥ 2.5 events continue to occur at slightly decreasing rate. band-pass filter, and, for surface waves, the range was 0.02-0.1 Hz. For a set of assumed focal depths, we used the grid search method to scan moment magnitude (step 0.01), strike, dip, and rake angle (step 5°) in order to find the result with the minimum fitting error. The obtained full moment tensor solution was then decomposed into a double-couple component (DC), an isotropic component (ISO), and a compensated linear vector dipole (CLVD) component. Finally, we produced reliable results for 12 M ≥ 3.5 earthquakes occurring between 2008 and 2019 ( Table 1). The mechanism solutions were further used for the inversion of the stress pattern and fluid overpressure required for causing fault reactivation of these earthquakes (see section 3.6 for details).

Hypocenter Relocation
The double-difference method (HypoDD) (Waldhauser & Ellsworth, 2000) was applied in order to relocate the hypocenters of earthquakes observed from 2008 through 2019. The velocity model used herein was that used by Wang et al., (2012) Wang et al., 2012). The maximum hypocenter separation was 10 km, and the number of minimum links per pair was eight. We relocated 696 hypocenters for more than 1,100 earthquakes for the period from 2008 to August 2019. With the aid of phase data obtained from the six temporary stations, the hypocenters of events from October 2008 to July 2011 were relocated with much better precision than for other periods. In total, the focal depth fell within a range of 1-15 km, with a peak at approximately 3 km, coincident with the injection depth.

Detection of Microseismicity
In order to refine our previous results (Z. Wang et al., 2018), we selected the relocated events with M ≥ 1 to serve as template events and applied the match and locate approach, which is based on waveform correlation methods (Peng & Zhao, 2009;Yao et al., 2017;Zhang & Wen, 2015), to detect and locate smaller events. In applying the match and locate approach, a 3-D spatial grid around a template event was searched, and the travel time differences between the template event and potential events were calculated. The differential travel times were used to slide the template waveform window along a continuous waveform in order to calculate the correlation coefficient (CC) for each trace. The average CC and signal-to-noise ratio (SNR) were calculated and used to detect the corrected events. The threshold vales of the CC and SNR for detecting an event were 0.3 and 10, respectively. Since the CC for events with different source time functions was low and the SNR was high, with the combination of a lower CC threshold and a higher SNR threshold, the detection capability was improved and false detection was reduced (Zhang & Wen, 2015). As a result, the magnitude of completeness from October 2008 to July 2011 was reduced from 1.5 to 0.75, which was very helpful for statistical analysis (Figure 2a).

Geophysical Research Letters
ETAS model, the total seismic occurrence rate λ(t) is divided into the background seismicity (also termed forced seismicity) λ 0 (t) and Omori-type aftershocks triggered by previous events v(t): where M min is the minimum magnitude and is greater than or equal to the lower cutoff magnitude of the catalog for completeness (Mc), α is a constant indicating the ability to generate aftershocks (a small α value indicates that small earthquakes dominate triggering), p is a constant representing the decay rate of the aftershock, and c is a constant in the modified Omori's law. Reflecting the irregular injection operations, the background seismic rate also demonstrates irregular patterns and thus could not be represented by a certain functional or empirical form. Following previous studies (Lei et al., 2013(Lei et al., , 2017, we applied a sophisticated algorithm (Marsan et al., 2013;Zhuang et al., 2002) to estimate the irregular forced seismic rate λ 0 (t) and other parameters (which are assumed to be constant over time). The method consists of the following: 1. Initially assuming a constant forcing rate λ 0 (t) = λ 0 , 2. Estimating ETAS parameters (K, α, c, p) that minimize AIC knowing λ 0 , 3. Updating the forcing rate based on these parameters, and 4. Repeating Steps 2 and 3 until convergence of all parameters.
During the third step, the probability ω i that the ith earthquake is a background earthquake was calculated for all events (Zhuang et al., 2002). Then, the forcing rate at the occurrence time t i of the ith event is calculated using a smoothing window for 2 * n e − 1 events centered on the ith event. The appropriate smoothing parameter n e is determined by minimizing AIC, which is defined as follows: where N is the number of earthquakes in the data set.
Finally, the fraction of forced seismicity (f.s.) for a given time window [S, T] is calculated as follows:

Stress Field Inversion and Estimation of Required Fluid Overpressure
Using focal mechanism solutions, it is possible to estimate the mean stress pattern that best matches the fault slips associated with the mechanism solution (Hardebeck & Michael, 2006;Lei et al., 2019aLei et al., , 2019bTerakawa et al., 2012) and to estimate the required fluid overpressure that causes an earthquake by assuming a critical stress state; that is, the most favorably oriented fault is critically stressed. Following a previous study (Lei et al., 2019a), the standard friction coefficient of 0.6 was assumed, and the vertical stress was calculated from the gravity overburdens.

Temporal Variability in Seismicity
We quantitatively investigated the statistical features of the seismicity in the RC gas field from 2008 to 2019. During the detailed observation period from October 2008 to July 2011, a catalog of microseismic events detected by the template matching method described in section 3.4 was also used. The seismic b value, which was estimated by applying the maximum likelihood method to the magnitude frequency distribution, was 0.75 for the period from October 2008 to July 2011 (Figure 2a), which was lower than the value of 0.83 for the total period (Figure 2b). During our study period, the largest earthquake, which had a magnitude of 5.1, occurred on 10 September 2010 (Figure 2c, d). Since the largest event magnitude recorded in the RC field was that of the 1997 M5.2 event, the M5.1 earthquake qualifies as a major event.
The results obtained using the ETAS model, including forced seismicity and other model parameters, are shown in Figure 2f. It was found that the background seismicity accounted for more than 70% of the total seismicity in the area, while Omori-type aftershocks were less important. It is important to note that the 10.1029/2020GL089569

Geophysical Research Letters
basin moves coherently with South China and has a very small present-day overall strain rate (G. Zheng et al., 2017). As a result, the background seismicity driven by tectonic loading is very low and can be ignored in comparison with the rate of induced earthquakes. We therefore refer the background seismicity as "forced seismicity" in this study. The ETAS results are in agreement with the results for the earlier injection stages   (Lei et al., 2008), indicating that the seismicity in the RC gas field was dominated by external force resulting from injections. During the study period, injection was still being performed at three wells: Luo-2, Bao-18, and Bao-46 (Figure 2f). Bao-46 is located more than 10 km from the seismic cluster and had less impact on seismicity in the study area. Following the shutdown of all nearby injection wells in 2013, as seen from the accumulated event number (Figure 2d) and the forced seismic rate, particularly the mean values indicated by blue lines (Figure 2f), the event rate clearly decreased clearly until the middle of 2015. At that point, the event rate rose again, and, on 27 December 2016, a M4.9 earthquake occurred, followed by a M4.0 event on the next day. Since then, seismicity remained at a relatively high level, comparable to that during injections. However, if we focus on earthquakes with M ≥ 2.5, the frequency of events shows no significant changes (Figure 1b). Interestingly, small earthquakes have a clear and rapid response to the termination of water injection, but larger earthquakes show significant sustained activity. As seeing from Figure 2f, the forced seismicity rate after termination of injection can be roughly represented by a modified Omori law for the decay rate of induced seismicity (Langenbruch & Shapiro, 2010).
where t is the elapsed time after the start of injection, t 0 is the time of shutting down injection, R a is the constant seismicity rate during injection, and p is constant. In our case, the injection history was very complex. By focusing on the last stage of injection with the highest seismic rate, the fitting results in Figure 6f are fairly well. In conjunction, the seismic b value is very small and exhibits a slight decreasing trend (Figure 2e). The remainder of this article will further analyze the postinjection seismicity based on other results.

Hypocenter Distribution and Mechanism Solution
The distribution of the injection wells, faults, and relocated earthquakes, including events detected and located by the template matching method, as described in sections 3.3 and 3.4, are shown in Figure 3a.
Overall, the relocated hypocenters were scattered along the anticline bounded bay mapped and unmapped faults that were NE stroked. As mentioned earlier, during the study period, Luo-2, Bao-18, and Bao-46 were active for injection. Of these, only the Luo-2 injection well was surrounded by hypocenters. Therefore, after 2008, water injection in the Luo-2 well has the greatest impact on the seismic activity in the area. Prior to this, water injection in the Luo-4 well played the dominant role in inducing earthquakes (Lei et al., 2008). The Luo-4 well was shut down in 2001 after injecting 700,000 m 3 of wastewater. From the distribution of relocated hypocenters, we found that during water injection mainly in the Luo-2 well, the earthquakes were primarily clustered in three areas, which were centered southwest of Luo-4, southeast of Luo-4, and the north of Luo-2 ( Figure 3a). Interestingly, the numbers of earthquakes near Luo-4 than those near Luo-2. Due to the sparse network, the hypocenter location error before 2008 is relatively large but is generally consistent with the results of the present study (Lei et al., 2008). Regarding this important result, we will give a comprehensive analysis in section 5.
The results of moment tensor inversion show that, except for one event, which shows a normal fault mechanism, other M ≥ 3.5 earthquakes indicate a reverse faulting mechanism, which basically agrees with preexisting faults (detected by seismic exploration) along the wings of the anticline (Figure 3b).
Due to the limited number of stations, the accuracy of the hypocenter depth of earthquakes after 2011 is rather poor, while the CMT depths, especially for the shallower earthquakes, were well constrained by the gCAP method, in which surface waves were used. As shown in the simplified geological section (Figure 3b), only one M ≥ 3.5 earthquake showed a focal depth of approximately 8 km, and all of the other earthquakes fell in the range of 1.6-4.6 km, centered at the depth of the completed gas reservoir in the Permian limestone formations, approximately 2.7 km below the surface. The Mw4.58 and Mw4.02 earthquakes that occurred in December 2016 have centroid depths of 3.0 and 4.5 km, respectively.

Stress Field and Overpressure Required for M ≥ 3.5 Events
Reliable focal mechanism solutions for 12 M ≥ 3.5 earthquakes (Table 1) enabled us to invert the mean stress pattern, including the stress shape ratio (ϕ) for the principal stresses and the directions of the principal stress axes (Figure 4). The resulting ϕ is 0.94, and the maximum horizontal principal stress axis is nearly horizontal (plunge = 2.6°, standard error = 1.7°) with an azimuth angle of 134°. Such a stress regime indicates that reverse faults with a dip angle of approximately 30°, independent of the strike, in this area are favorable to reactivation.  Table 1), faults, injection wells, and seismic stations overlapped on a shaded topographic image. The large arrow indicates the direction of the maximum horizontal principal stress (SHmax) obtained in the present study. (b) A simple geological section overlapped with the earthquake hypocenters, focal mechanism solutions, and faults. The plotted 1-D velocity model obtained by receiver function method (X. Wang et al., 2012) was used for relocation.
By assuming critical stress on favorably oriented faults, the overpressure required to activate the source fault for each M ≥ 3.5 earthquake was estimated. In calculating the overpressure, the stress change caused by solid deformation was ignored. The focal mechanism solutions were projected on a normalized Mohr diagram (Figure 4d), and the estimated overpressure values are also listed in Table 1. Only two events required an overpressure greater than 7 MPa. For the other events, the required overpressure was lower than the maximum injection pressure. With consideration of uncertainty and possible inhomogeneity of the stress field, we believe that fluid pressure, which was raised by long-term injection over more than three decades, was the main driving factor for fault reactivation, and a high level of seismicity occurred during injection and after injection was stopped.

Discussion and Conclusions
Similar to the previous injection-induced seismicity in the study area and other nearby injection-induced seismic activities, the ETAS results showed a smaller value of α and a larger fraction of forced seismicity. As described above, the background seismicity driven by tectonic loading could be ignored, and the forced seismicity was mainly driven by injections and thus was expected to better respond to injection conditions. In agreement with other induced earthquake sequences in the Sichuan base, the fraction of forced seismic activity was very high, over 70%, even after the injection wells were closed. The α value (<1) was significantly smaller than the typical value (>2) for tectonic earthquakes, which is dominated by cascade triggering. Thus, Omori-type aftershock productivity in induced seismicity was very weak, and forced seismicity was dominant.
Following the shutdown of all injection wells by 2013, there was a clear decreasing trend in earthquakes with a seismicity of M < 2.5, but the occurrence rate of M > 2.5 earthquakes remained unchanged, at least for 4 years. The seismic b value decreased from approximately 1.3 to 0.7 from the very beginning of injection to 2006 (Lei et al., 2008). Since 2008, the low b value was maintained, and a further decrease was observed after all wells were shut down. This was close to the b value of 0.72 for Changning Ms6.0 earthquake sequence related to long-term salt mining (Lei et al., 2019b). As we know, the low b value was a high seismic hazard indicator. Together with other results for b values, the Sichuan Basin demonstrated a significantly lower b value than the global mean value of 1.0, which should be considered in seismic hazard models.
Hypocenters were clustered in zones clearly related to preexisting fault zones. When major injections changed from one well to another well, and even after all injection wells were shut down, the spatial distribution did not change much.
Beginning in April 2015, the seismicity, which had been decreasing, once again increased, with two Mw > 4 events occurring in December 2016, followed by multiple aftershocks, which delayed the decline of seismicity. A similar phenomenon was also observed in other areas, such as north-central Oklahoma in the United States, where seismicity has decreased markedly in response to reduced saltwater injection. However, a short-term increase of seismicity was observed following larger magnitude (M > 5) earthquakes after the injection rate had already been significantly reduced (Langenbruch & Zoback, 2016).
The centroid depths obtained by applying the gCAP method were centered at the depth of the depleted reservoir. The estimated fluid overpressure values required to cause the two M4 events were 1.2 and 2.4 MPa, respectively. The overlying shale layer and underlying mudstone layer in the gas reservoir act as a fluid diffusion barrier. The depleted gas reservoir has good sealing properties, and the overpressure fluid resulting from injection can only diffuse along preexisting faults and thus could conceivably persist for a relatively long time. However, the permeability within the reservoir is relatively high. Therefore, especially during later stages, injection at wells within the reservoir could raise the fluid pressure by a large amount. In fact, based on data from the pressure gauge installed at the well head, the well head pressure for Luo-4 (closed in 2001) increased to 3.2 MPa on 27 September 2009 from 0.4 MPa on 1 August 2008. After the Mw7.9 Wenchuan earthquake occurred on 12 May 2008, injection was suspended for a few months and was then carried out extensively at Luo-2 and at other wells 10 km away.
Considering all the results, we suggest that the postinjection seismicity resulted from rupture along preexisting faults driven by pressure diffusion from the overpressure reservoir, in which the fluid pressure was raised by the long-term injection of wastewater. This mechanism is basically the same as that during injection. Until the overpressure decreases to a safe level, there is the potential for a high seismic hazard at sites that experienced destructive earthquakes.
Although the speed of fluid moving out the well-sealed gas reservoir may be lower, as the fluid pressure spreads, delayed seismic activity should generally decline. However, we had no chance to verify this assumption because shale gas development in the synclines at the north and south sides of the RC gas field has begun. Thus far, there are several well pads of shale gas within a distance less than 20 km from the center of the present study area, and HF has been conducted in horizontal wells for evaluation in the past 2-3 years. These activities might have contributed to the increased seismic activity since 2017, and it is important to monitor future seismic activity in this area.