The 2019 MW 5.7 Changning Earthquake, Sichuan Basin, China: A Shallow Doublet With Different Faulting Styles

The increased seismic activity of the last ~10 years in Changning county of Sichuan Province comprised just small (mostly ML < 5.0) injection‐induced earthquakes. The MW 5.7 earthquake on June 17, 2019, is the largest event ever reported there. Moment tensor of the mainshock was remarkably dominated by a compensated linear vector dipole. We resolve its fine structure showing it was a doublet, allowing approximation by a thrust‐ and strike‐slip subevent. The mainshock nucleated as thrust faulting, which (together with the largest aftershocks) can be linked with previously known reverse faults, favorably oriented to regional stress field. Contrarily, the strike‐slip segment of the mainshock, less favorably oriented, was probably facilitated by elevated pore pressure due to previous injections. Shallow active strike‐slip faulting, not yet mapped in the region, is a new feature, important for future hazard assessment.


Introduction
Earthquake doublets provide clues for understanding complex fault systems (e.g., Danré et al., 2019;Hicks & Rietbrock, 2015;Lay et al., 2013). Particularly challenging are closely spaced doublets whose subevents have unequal focal mechanisms, separated just by a few kilometers and a few seconds, intuitively suggesting a causal relationship between the involved faults. Global low-frequency moment tensor (MT) solutions provide only an overall characterization of such events, often featuring large non-double-couple (non-DC) components. This paper makes use of superior quality of broadband Chinese seismic networks and reveals one such a rare closely spaced doublet in a salt-mining district of Sichuan Basin. It points to the existence of segmented faulting, especially to an unmapped shallow strike-slip (SS) fault, important for future seismic hazard assessment in this densely populated region of the world.
Sichuan Basin belongs to major gas and oil resources in China. The annual natural gas production should be as high as 5 × 1,010 m 3 in 2030 (Ma, 2017). The basin, important also for its agriculture, supports a population of over 100 million. Three shale-gas fields in the southern part of the basin, Shangluo, Zhaotong, and Changning, were reported as sites of induced seismic events (Lei et al., 2017(Lei et al., , 2019aMeng et al., 2019). The Sichuan Basin also belongs to the main salt-producing areas, using the solution mining method, with injection wells drilled to~3-km depths. Increased seismicity has been recognized in the last~10 years, likely related to the long-lasting water injection in the Changning salt mines (Sun et al., 2017). Given that potential seismogenic faults are poorly known, events of magnitude greater than 6.0 cannot be ruled out (Lei et al., 2019a).
Indeed, on June 17, 2019, according to the China Earthquake Networks Center (CENC), an M 6.0 earthquake occurred in Changning county (see Figure 1). Hereafter, according to current practice in China, we use a generic symbol M, denoting either local magnitude M L or an empirical estimate of the surface-wave magnitude M S = 1.13, M L −1.08 (Chen et al., 2014), for events below or above M L~5 , respectively. The earthquake killed 13 people. More than 200 were injured, and a large number of buildings were damaged (Yi et al., 2019). This is another M > 5.0 earthquake occurring in 2019 in the Gongxian-Changning region, Sichuan Province, after the M L 5.3 earthquake on January 3, 2019, studied by Lei et al. (2019a). Only nine M ≥5.0 earthquakes have been previously known to occur in the Changning and adjacent areas: Yibin M 5.5 in 26 BC, Gaoxian M 5.5 in 1610, Nanxi M 5.0 in 1892, Zigong M 5.7 in 1896, Weiyuan M 5.0 in 1905, Jiangan M 5.0 in 1936, Zigong M 5.0 in 1954, Fushun M 5.0 in 1959, and Yibin M 5.4 in 1996(China Earthquake Administration, 1995. The June 17, 2019, M 6.0 earthquake occurred at a distance of a few kilometers from the salt-mining wells (Lei et al., 2019b). The mainshock was followed by an aftershock sequence, including four M 5+ events, forming a clear line (subvertical planar) structure,~20 km long, trending southeast-northwest ( Figure 1b). This trend coincides with the major tectonic element of the area, the Changning anticline. There are many small faults in the anticlinal structure, often associated with a typical caprock/fold system, most of them being reverse faults (Qian & Tang, 1992;Sun et al., 2017).
The 2019 M 6.0 event is noteworthy for a very large compensated-linear-vector-dipole (CLVD) component of its single-point source model; for example, Global Centroid Moment Tensor (GCMT) Project reported the CLVD of −98%. Later, in the paper, we demonstrate that this significant non-DC character of the event may explain great scatter between the strike/dip/rake angles reported by several agencies, for example, Figure 1. The 2019 Changning earthquake sequence. (a) Sichuan Basin, China. Focal mechanism of the mainshock, provided by two agencies (GCMT and IEF) and this paper. (b) The mainshock (M 6.0) and four M 5+ aftershocks are shown by stars and the focal mechanism symbols; the smaller aftershocks are plotted as pink points. The earthquakes were relocated in this paper using the CENC phase data and HypoDD method (Waldhauser & Ellsworth, 2000). For the two largest shocks, most deviating from a double-couple (DC) model, their focal mechanism diagram is split in two: the major and minor DC. Major tectonic elements are plotted as lines: anticlines (blue), synclines (green), and the known active faults (red). Inset in panel (b) demonstrates the previous 10-year activity of the region, M > 0 (CENC). (c) The first-motion polarities of the mainshock (compressions in red and dilatations in green) clearly disagree with the full-MT solution, but they are very well explained by the thrust-fault nucleation subevent; the diagrams in panel (c) are not scaled with the moment. thrust faulting (TF) of GCMT and SS of Institute of Earthquake Forecasting (IEF, China Earthquake Administration), shown in Figure 1a. It also explains why first-motion polarities do not fit with the full MT; see Figure 1c (top).
While the contemporary evidence that injection can cause earthquakes is undoubted (Ellsworth, 2013;Grigoli et al., 2018), very little is known about relatively large earthquakes occurring in regions previously activated by injection activities and thus far experiencing only smaller events. How is their rupture nucleating, propagating, and eventually causing aftershocks? Therefore, we focus on the source process of the M 6.0 event and its M 5+ aftershocks. We resolve spatiotemporal complexity of the mainshock and causative faults. Finally, we demonstrate that the SS segment of the mainshock is suboptimally oriented to the regional stress field, indicating possible pore pressure effects.

Source Models of the Mainshock
Our waveform data of the M 6.0 mainshock consist of 12 broadband records at epicentral distances~120-210 km (see Figure S1 in the supporting information). The station selection was based on data quality (low noise and absence of disturbances; Zahradník & Plešinger, 2010) and azimuthal coverage. Several velocity models were tested, providing qualitatively the same results, for example, Crust 2.0 (Bassin et al., 2000). The waveform modeling is performed with established methods (Sokos et al., , 2016Zahradnik et al., 2017) and recently developed uncertainty assessments, detailed in Text S1.
The single-point source models pointed to a highly non-DC solution. A typical full MT is presented for the range of 0.03-0.06 Hz in Figure 1 and Tables S1 and S2. Without loss of generality, hereafter, we demonstrate results for a regional velocity model (Zhao & Zhang, 1987), shown in Table S3, the same model as used by Sun et al. (2017). The centroid depth is definitely very shallow, 3-5 km. In this depth range, the DC part of the full MT is low (~10%), and it varies from TF to SS faulting (see Figures 2a and S2). The formally best-fit solution is found at the depth of 4 km, consisting of DC = 20%, CLVD = −68%, and ISO = −12%. The GCMT solution was similar in terms of its high negative CLVD percentage (−98%) but placed the event much deeper at 12 km. The waveform fit of our model is very good (variance reduction [VR] = 0.90, Figure S1a). A formal decomposition of MT into two DC parts (Jost & Herrmann, 1989) yields the moment ratio of the minor DC to major DC equal 0.61, with strike/dip/rake angles (hereafter s/d/r, in degrees) = 336/48/86 and 114/86/−1, respectively.
At shallow depths and relatively low frequencies, the MT resolution is limited and must be accompanied by an assessment of uncertainties and parameter trade-offs (Halló et al., 2017(Halló et al., , 2019Vackář et al., 2017). Using these new methods, validated also on nuclear tests (Liu, Li, Zahradník, Sokos, & Plicka, 2018;Liu, Li, Zahradník, Sokos, Liu, & Tian, 2018) and demonstrated in Figure 2b, we find that nodal planes are strongly nonunique. While the P-axis is very well resolved (P_azimuth 69°, P_plunge 3°), the T-axis is arbitrary in the orthogonal plane. Well resolved also are the ISO percentage (almost vanishing), the large negative CLVD%, and the low DC%. The event is basically a CLVD whose major axis, almost horizontal, is the only axis well resolved ( Figure S3).
Extending the frequency range of the waveform inversion, the source appears to be composed of several subevents. Here we present a typical example of the results for the frequency range of 0.03-0.10 Hz; a further increase would deteriorate waveform fit due to the inherent inaccuracy of the velocity model. In this extended range, the single-point full-MT solution decreases VR from 0.90 to 0.80, due to increased complexity of waveforms (Figures S1a and S1b). Assuming that the rupture process might have included several tectonic faults, we focus on seeking subevents in DC-constrained mode. After experimentation with various geometries of the trial source grids, we present multiple-point source models for a horizontal line grid at a depth of 4 km ( Figure 3). We choose azimuth of 120°, close to the aftershock alignment. Similar results were obtained for a line deepening from 2.4 to 6.1 km in the northwestern direction. Using iterative deconvolution , the model suggests two subevents ( Figure 3a): a dominant SS, slightly NW shifted relative to the epicenter, and a weaker thrust fault, shifted to SE. Their moment ratio TF/SS is low, 0.35. Perturbing the parameters (e.g., the line azimuth, grid increment, or frequency), the space-time separation and size of the subevents are varying. However, the two focal mechanisms remain very stable (Table S1), close to the formal major/minor DC components of the single-point source model. Another strong indication of the stability is provided by the modified iterative deconvolution in which the moment release is artificially "slowed down," providing models closer to finite-fault slip inversions (Zahradník & Gallovič, 2010). Although the DC-constrained focal mechanisms are free for all calculated subevents, we again obtain just two source types: the SS and TF (Figure 3b). Here, in Figure 3b, the TF rupture again appears earlier than SS, shifted to SE, but the moment ratio TF/SS is larger,~1.2.
The best insight into space-time separation and moment ratios of subevents is provided by a joint search for the TF-SS source pairs (Zahradnik & Sokos, 2014), using non-negative least squares (Lawson & Hanson, 1974). Here (Figures 3c-3e) we fix the s/d/r angles of the subevents and use frequency 0.03-0.10 Hz, as above. Twelve of the two-point DC-constrained models shown in Figure 3c have VR > 0.80, greater than the single-point full-MT solution. A stable feature is the position of the SS event, preferred in points 7 to and 9, and the significant involvement of the thrust-fault subevent: The moment ratio of the TF/SS subevents varies from 0.3 to 1.2. Their tensor sum yields CLVD between −46% and −90%, in good agreement with low-frequency single-point full-MT model (Table S1). An exact thrust-fault subevent position and its timing are not unambiguously resolved (Figures 3d and S4). Nevertheless, most likely, the SS subevent is delayed and northwestward shifted~3-9 km relative to the TF subevent. The TF subevent can be situated at the hypocenter (point 7). The formally best-fitting two-point model is the pair (SS,TF) with position (8,6), characterized by VR = 0.84. Its TF/SS moment ratio is 0.72, the subevent separation of 6 km, and the momentrate function (Figure 3e) indicates a 2-s time delay of SS relative to TF. If considering just the SS subevent of this pair, VR would drop to 0.72 ( Figure S1d).
To prove that the mainshock source process started with TF, we employ independent information about the initial (nucleation) part of the source process, provided by first-motion polarities. We visually inspected more than 120 stations and kept the 109 most reliable readings (Figures 1c and S5). These polarities have been projected onto focal sphere using the same velocity model as in the waveform inversion, with several trial depths. For an arbitrary depth between 3 and 5 km, the observed polarities are well fitted with the TF mechanism. At the same time, neither the SS mechanism nor the full-MT single-point model can (a) Correlation between observed and synthetic waveforms versus depth. The DC focal mechanism symbols (depth of 1-10 km) are color coded relative to the DC percentage. Shown for the depth of 3-5 km is also the full moment tensor. (b) Uncertainty shown by random sampling of the posterior probability density of model parameters (nodal lines and P-T axes, the source-type plot (Hudson et al., 1989), and histograms of the DC, CLVD, and ISO). Due to shallow depth and the remarkably non-DC character, only the P-axis is resolvable. It demonstrates a limited physical meaning of the single-point low-frequency models, if presented without uncertainties. explain the data. Worth to mention that if directly inverting the polarities (without any waveform inversion), using FOCMEC code (Snoke, 2003), we obtain a mechanism just negligibly deviating from our TF, Kagan angle <20°(see Figure S6). Therefore, we conclude that the mainshock was a closely spaced TF + SS doublet, nucleating as a thrust. The term "doublet" accentuates the two involved focal mechanisms rather than two distinct episodes. Our approximate model does not rule out a possibly continuous source process.

Discussion
The mainshock M 6.0 was followed in 41 minutes by an M 5.1 aftershock, shifted~10 km to NW. Then, 8 hours later, an M 5.3 occurred close to the largest shock. Finally, 5 and 17 days later, the M 5.4 and M 5.6 aftershocks, respectively, appeared northwestward of the M 5.1; see Figure 1b. The earthquakes were "jumping forth and back" in the SE-NW direction along the mapped anticline. Most of the observed smaller aftershocks followed the same trend (also Figure 1b), without any systematic space-time migration, perhaps except the dominance of the latest events at NW ( Figure S7). We have shown that significant perhaps  (a), formally shifted for the presentation, color coded relative to time, and circle sized according to moment. (e) The moment-rate time function for the best-fitting source pair (8,6) with the moment ratio TF/SS = 0.72, CLVD = −83%, VR = 0.84. The blue and red lines refer to SS and TF, respectively; the sum is shown in green. For more details, see Figure S4. The mainshock undoubtedly comprised both the thrust-and strike-slip faulting. even dominant part of the mainshock was an SS subevent, whose one nodal plane (strike~120°) agrees with this NW-SE structure. Therefore, we suggest that the mainshock included an unmapped left-lateral SS, aligned with the mapped anticline structure. The GCMT solution (Figure 1a) does not indicate an SS faulting geometry in the best-DC nodal lines. Nevertheless, a formal decomposition of the GCMT solution into the major and minor DC parts does provide an almost same SS and TF mechanisms as found in this paper, with the TF/SS moment ratio of 1.04. Note however that our physically based retrieval of the two subevents is deeper than the formal decomposition. Besides the SE-NW trend, the aftershock distribution indicated also shorter faults, for example, an N-S trending structure, seen in Figure 1b near the M 6.0 epicenter.
Seismic activity parallel with the anticline has been well documented by the 2008-2018 seismicity (Figure 1b, inset). In Sichuan Basin we face a caprock/fold system in the uppermost 3-6 km, slowly slipping on the underlying basement (a "decollement" process). The anticlines of the caprock/fold are rich in symbiotic faults, whose majority are steeply dipping thrust faults, and earthquakes often occur on such faults inside the caprock (Qian & Tang, 1992). This is consistent with the observed event depths in our study (Table S1); except one M 5+, the events are shallower than 5 km. However, the cited authors also claimed that the caprock/fold system cannot accumulate a large amount of strain, which contrasts with the relatively large magnitude of the 2019 earthquake.
A closer look at the 10-year seismicity reveals some deviations from the general SE-NW seismicity trend. It points to a complex system of active faults, also indicated by deviatoric focal mechanisms of the M 5+ events (Figure 1, Table S1), whose nodal planes are not parallel with the NW-SE anticline. Three TF mechanisms have a very high DC percentage. All these mechanisms share similar azimuth of the P-axis as the mainshock, 47-84°. This finding is in agreement with geological evidence that the shale and sandstone rocks of the anticline include small faults whose majority are thrust faults (Sun et al., 2017). In this sense, we interpret the TF nucleation of the mainshock as an expected rupture process.
According to our MT analysis, one of the M 5+ aftershocks is different (M 5.6 in Figure 1b and Table S1). It is the most delayed and the deepest one, occurring on July 4, near the NW termination of the supposed mainshock rupture. This M 5.6 TF event has again a significant non-DC character, but, opposed to the mainshock, its CLVD is positive (~80%). As such, its formal major/minor DC decomposition shares a common T-axis, and both DC parts are thrust faults (Figure 1b). Due to smaller magnitude, this event is more challenging than the mainshock; thus, we cannot prove whether it actually consisted of two subevents or not. The event should deserve more attention elsewhere, possibly illuminating the region where the major anticline is changing from linear to a bended structure. Not only such a bend could have stopped the mainshock rupture, but also it might be generally less active as suggested by the 2008-2018 activity in Figure 1b (inset).
Finally, in Figure 4, we compare the 2019 sequence with the regional stress field. The stress was calculated from 39 focal mechanisms of 2008-2016 (Liu et al., 2019). Using the given stress, similarly as Hainzl et al. (2016), we resolve shear and normal tractions on both nodal planes and identify the planes with greater shear instability (defined by  (Liu et al., 2019;Ran et al., 2008;Wang et al., 2015;Yi et al., 2019;Zhu & He, 2014). The upper inset shows the pre-2019 regional stress field (Liu et al., 2019)-the principal axes (azimuth/plunge) and the shape ratio R. Faults as in Figure 1. The bottom inset is a ternary diagram (Frohlich, 1992). (b) The events of panel (a) are displayed on Mohr's diagram, using their nodal planes of the greater shear instability. The shear and normal tractions correspond to the vertical and horizontal coordinates, respectively. The fracture criterion for the friction coefficient of 0.6 is shown by the tangent line. (c) The same events as in panel (b) but using the more stable nodal planes. Note that the stress field is not calculated here. Focal mechanisms serve just for resolving tractions in the a priori known stress field. It is inferred that the left-lateral SS subevent of mainshock (M 6.0-SS) belongs to a few earthquakes that are most distant from the tangent line, that is, the least favorably oriented to the stress field. Those six events are depicted by the dotted ovals, and all of them are SSs or very close to SS mechanism (red dots in the ternary diagram). Their fractures could be assisted by the pore pressure diffusion due to water injection activities of the previous years (Sun et al., 2017). Vavryčuk, 2014). Importantly, the SS subevent of the mainshock (M 6.0-SS) belongs to a few six events, all of the SS type, least susceptible to failure (most distant from the Mohr's circle envelope), that is, less favorably oriented to the stress. Thus, we recall the water injection effects of the previous years (Sun et al., 2017) and suggest that the SS component of the mainshock was likely facilitated by the elevated pore pressure. Additional effects, for example, poroelastic coupling and time-dependent nucleation (Segall & Lu, 2015), would require specific modeling that is beyond the scope of the present paper. Static Coulomb stress perturbation caused by TF fault, allowing rupture of SS fault, is inconclusive ( Figure S8). Due to the space-time proximity of the two subevents, a likely scenario appears to be a continuous evolution of a single stressfluid-fault system, involving two segments of different faulting styles. From the longer-term perspective, the 2019 earthquake sequence might have been affected by the M 5.7 Xingwen earthquake (December 18, 2018, 28.24°N, 104.95°E), thought to be induced by hydraulic fracturing in the southerly adjacent Shangluo shale gas area (Lei et al., 2019b).

Conclusions
The Changning county in Sichuan Basin, China, has been known for the frequent occurrence of small (M L < 5) injection-induced seismic events, particularly well documented in the 2008-2018 period. Generally, earthquakes follow a major tectonic element-the SE-NW trending anticline. The anticline is a caprock/fold system, 3-6 km thick, slowly slipping on the underlying basement. The earthquakes occur inside this depth range. The system was considered unable to accumulate large strains, but the recent 2019 M W 5.7 (M 6.0) event indicates that moderate earthquakes are possible.
We analyzed MTs of the mainshock and all major (M 5.1-5.6) aftershocks, finding shallow centroid depths and various focal mechanisms. The aftershocks appeared to be basically thrust faults. This finding is consistent with previous studies that declared small thrusts to be typical for the region. The mainshock was composed of two subevents: TF and SS. This complex faulting geometry caused the notable non-DC character of the single-point low-frequency source models and explains why focal mechanisms varied among several reporting agencies. The significant left-lateral SS subevent, striking parallel with the anticline and aftershocks, is a new finding. Long seismoactive SS faults have not been mapped in the region yet and should be considered in the seismic hazard assessment.
The mainshock nucleated as a thrust fault, well oriented to the regional stress field. Then, within a few seconds and a few kilometers, it progressed toward the northwest as an SS rupture. The SS faulting, less favorably oriented to regional stress than TF, was probably facilitated by the elevated pore pressure due to previous water injections. The adopted two-point model is an approximation of a possibly continuous rupture process, involving a smooth transition between the TF and SS segments. A deeper insight into the observed fault-stress-fluid interactions would require extensive finite-source dynamic modeling (Galis et al., 2017;Gallovič et al., 2019).