Local stress field inverted for a shale gas play based on focal mechanisms determined from the joint source scanning algorithm

1.   Introduction

Microseismic events and their focal mechanisms have been studied extensively in the oil-gas industry to evaluate hydraulic-fracturing treatment or to characterize induced seismicity due to oil/gas production (Duncan and Eisner, 2010; Maxwell, 2011; Close et al., 2012; Li JL et al., 2012, 2013; Anikiev et al., 2014; Grechka et al., 2016; Kuang WH et al., 2017; Shao J et al., 2019).

The technologies have evolved from inverting for seismic locations and focal mechanisms independently to jointly invert both. For example, Li JL et al. (2011) modeled waveforms to determine locations (X,Y,Z) and focal mechanisms (Strike, Dip, Rake). Willacy et al. (2018, 2019) modeled waveforms to determine locations and focal moment tensors simultaneously. In additions, Rodriguez et al. (2012) also include a fourth term in focal mechanism to include the tensile component. Chambers et al. (2014) developed a method to combine moment tensor reverse time imaging with diffraction stacking and simultaneously solve for source locations and moment tensor components. Zhang X and Zhang J (2016) used a search engine scheme to compare real event waveforms with a pre-computed waveform database. The best match gives the location and focal mechanism of a seismic event. All of these methods work only for waveforms with high signal to noise ratios (S/N). Kushnir et al. (2014) estimated both sets of source parameters based on statistical theory. Anikiev et al. (2014) and Zhebel and Eisner (2015) present a method to invert for moment tensor at each origin time and imaging point and then use the moment tensors to correct polarity before stacking to locate event. Forghani-Arani et al. (2013) used the scanning algorithm to find a best focal mechanism that can suppress the noise.

Liang CT et al. (2016) introduced a method that extends the traditional source scanning algorithms (SSA, Kao H et al., 2006) method to jointly scan source locations and focal mechanisms (jSSA). It differs from Anikiev et al. (2014) and Zhebel and Eisner (2015) in that the focal mechanism parameters are included in the search process while the other two studies solving the moment tensors before scanning locations. Yu YY et al. (2018) conducted a systematical testing using both real field data and synthetic data to evaluate the accuracy of the jSSA and its dependence on noise level, source depths, focal mechanisms and other factors. They found that the surface-based arrays have better constraints on horizontal location errors and angular errors of P-axes (within 10 degrees, for S/N≥0.5). For sources with varying rakes, dips, strikes and depths, the errors are mostly controlled by the partition of positive and negative polarities in different quadrants. More evenly partitioned polarities in different quadrants yield better results in both locations and focal mechanisms. Nevertheless, even with bad resolutions for some focal mechanisms, the optimized jSSA method can still improve location accuracies significantly. The application to field data also yielded not only more events, but also much better delineation of microseismic events.

Both source locations and focal mechanisms can be integrated to study the faulting geometry and stress distribution regionally or locally (Herrmann et al., 2011; Rodriguez et al., 2012; Anikiev et al., 2014; Forghani-Arani et al., 2013; Yang YH et al., 2018). Michael (1984, 1987) developed a linear inversion method to find the stress tensor in a sub-region inside which multiple events are available. Yang YH et al. (2018) and Luo Y et al. (2015) used this method to obtain a gridded stress field for the eastern Tibetan Plateau as well as the region near the Lushan earthquake. Kuang WH et al. (2017) used this method to find the principal stress axes for each fracturing stage. They integrated shear stress, pore pressure and normal stress to compute Coulomb failure stress to find the actual faulting planes. These methods are preferable when the whole region can be divided into sub-regions in which simple stress field exists in individual subregion.

In this study, we applied the jSSA method to a hydraulic fracturing dataset to obtain locations and focal mechanisms of microseismic events simultaneously, and then analyze the maximum compressional stresses of all these events. Finally, a uniformly gridded stress field (50 m by 50 m) was inverted. Thanks to the doubled number of events, plus the stress field extracted from the focal mechanisms, we are able to identify a complex faulting system and complex regional stress distribution.

2.   Method

In this section, we summarize the jSSA method introduced by Liang CT et al. (2016) for the convenience of discussion, and then summarize the method to extract stress attributes based on focal mechanisms.

For SSA or jSSA methods, the area where possible events occur (for example, an orthorhombic region around a perforation point) is divided into 3D grids. Every grid is treated as a virtual source. In addition to the SSA method, which assumes all virtual sources are of explosive types, the jSSA also scans all possible combinations of strikes, dips and rakes for every virtual source. The brightness function (a function to show how “bright” an event is compared to the noisy background, and it can be treated as signal to noise ratio of the stacked trace) ${{b}}_{{i}}$ is then computed as

where ${j}=\mathrm{1,2},\dots ,{N}$ are the indices of N receivers; ${{u}}_{{j}}$ represents the waveform recorded on the jth receiver; ${{t}}_{{i}{j}}$ is the traveltime from the ith virtual source to the jth receiver, while ${{t}}_{{i}}^{0}$ is the origin time of the virtual source; ${{s}}_{{i}{j}}$ represents the waveform polarity from the ith virtual source at (${{x}}_{{i}}^{\mathrm{s}},{{y}}_{{i}}^{\mathrm{s}},{{z}}_{{i}}^{\mathrm{s}}$) to the jth receiver at (${{x}}_{{j}}^{\mathrm{r}},{{y}}_{{j}}^{\mathrm{r}},{{z}}_{{j}}^{\mathrm{r}}$), and it is a function of the locations of sources and receivers, as well as the focal mechanisms of (${{\phi }}_{{j}},{{\delta }}_{{j}},{{\lambda }}_{{j}}$). ${{\phi }}_{{j}},{{\delta }}_{{j}},{{\lambda }}_{{j}}$ are strike, dip and rake for the ith virtual source, respectively. In the application, the explosive source is also scanned by setting ${s}_{ij}=1$ for all receivers. With this modification, Equation (1) converts the waveform data in observed time-space domain (OTSD) to the “source time-space-mechanism-domain (STSMD)”. As pointed out in Liang CT et al. (2016), the times in OTSD and STSMD are the arrival time at receivers and the origin time of seismic sources, respectively. After travel time corrections, polarity corrections and stacking, the R value at each time point is computed:

where ${{R}}_{{k}}$ is the signal to noise ratio (SNR) of the kth time point on the stacked trace. Assuming the ${{M}}_{{s}}$ points after the kth point are signals, while ${{M}}_{{n}}$ points before the kth point are background noises, then the numerator and the denominator of the right hand side of the Equation (2) represent the root mean square amplitudes of the signal and background noise, respectively.

The virtual source that gives the largest R in a time window and R>3.0 are assumed as a solution. The choice of threshold for R is based on synthetic tests and real data application. For example, if the threshold is less than 3.0, the spatial and temporal distribution of obtained events appear to be random.

The velocity model used to compute the ${{t}}_{{i}{j}}$ can be any 3D model with or without anisotropy as long as the ray tracing was conducted in advance and they can be read in for normal moveout (NMO) correction and stacking following the Equation (1).

Follow Stein and Wysession (2003), with focal mechanism parameters (strikes, dips, rakes) = $(\mathrm{\phi },\mathrm{\delta },\mathrm{\lambda })$ obtained from jSSA, the normal vector of the fault ($\overrightarrow{\boldsymbol{N}}$) and slip vector ($\overrightarrow{\boldsymbol{S}}$) can be found by Equations (3) and (4), respectively.

With the assumption that the fault slips in the direction of the maximum shear stress on the fault plane, Michael (1984) established an algorithm to invert the stress tensors for each subarea by fitting source mechanism of all earthquakes in that subarea. The theory is summarized below for the convenience of discussion.

Assuming the stress tensor around the fault is ${\boldsymbol{\sigma }}$ and the slip vector on the fault is $\overrightarrow{\boldsymbol{S}}$ (Equation 3) with the unit or directional vector of $\overrightarrow{{\boldsymbol{s}}}$. The unit or directional normal vector of the fault is $\overrightarrow{{\boldsymbol{n}}}$. The shear stress vector $\overrightarrow{{\boldsymbol{\tau}} }$ can be found by

The equality of the slip direction $\overrightarrow{\boldsymbol{s}}$ and the shear stress direction $\overrightarrow{\boldsymbol{\tau} }$ establishes the relationship between the stress tensor $\mathrm{\boldsymbol{\sigma} }$, the normal vector of the fault $\overrightarrow{\boldsymbol{n}}$ and the unit slip vector $\overrightarrow{\boldsymbol{s}}$.

Given the focal mechanisms inverted from jSSA, $\overrightarrow{\boldsymbol{s}}$ and $\overrightarrow{\boldsymbol{n}}$ can be found once $\overrightarrow{\boldsymbol{S}}$ and $\overrightarrow{\boldsymbol{N}}$ are computed using Equations (3) and (4). The stress tensor is linearly related to $\overrightarrow{\boldsymbol{s}}$ and $\overrightarrow{\boldsymbol{n}}$ in the form of Equation (6):

where ${\boldsymbol{G}}_{3l\times 5m}$ is the sensitivity matrix given by Michael (1984) and the components of G are functions of the unit normal vector ($\overrightarrow{\boldsymbol{n}}$) of faults. Where ${l}$ and m are the number of earthquakes and number of blocks, respectively. For shear sources, there are only five independent stress tensor components in Equation (6) with an implicit assumption that the isotropic stress vanishes.

Hardebeck and Michael (2006) introduced the spatial and temporal damping parameters to damp the uncertainties of data and smooth out the abrupt stress variation between neighboring subareas. This technique has been used in many areas to study regional stress fields (e.g. Zhao GZ et al., 2012; Yang YH et al., 2016, Yang YH et al., 2018) and local stress fields in hydraulic fracturing (e.g. Kuang WH et al., 2017). This is the first study to get a uniformly gridded stress field based on the focal mechanism yielded by the jSSA method.

3.   Application

3.1   The locations and focal mechanisms of microseismic events

We applied the jSSA method to find locations and focal mechanisms of microseismic events in one shale gas production field near the Anyue in the central part of the Sichuan basin, China (Figure 1). The basin can be divided into three regions (Lei XL et al., 2017): a depression region to the west of the Longquan fault (LQF), an uplifted region (Central Sichan uplift) between the LQF and the Huayingshan fault (HYSF) and a fold belt to the east of the HYSF. Under the NW-SE oriented compressional stresses from the Tibetan Plateau, a series of NE-SW trending mountain belts, often accompanied by the NE-SW striking faults and folds [e.g. Longmensan fault (LMSF), LQF, HYSF], are observed on the borders of the basin as well as inside the basin. The depth extension of the faults inside the basin is generally less than a few kilometers, but the boundary faults such as the LMSF and the HYSF extend as deep as 20 km. Some faults may cut through the whole crust. The open circles are earthquakes with M>3.0. A linear seismicity (marked by green dashed lines) is located to the SW of the Sichuan basin between LQF and HYSF, and it is associated with a buried fault (Lei XL et al., 2013, 2017). The production field near the Anyue appears to align along with the northeastward extension of the buried fault. One M_S5.1 earthquake occurred to the northeast of the production field. According to Lei XL et al. (2017), the focal mechanism of this earthquake is 223/48/122 for strike, dip and rake, respectively. The strike of the slip is consistent with the buried fault (green dashed line). The buried fault may extend northeastward to form the axis of the Central Sichuan uplift (Liu SG et al., 2021). The shale gas is reserved in Silurian Longmaxi formation (Zeng QC et al., 2018). For the Longmaxi formation, the total organic carbon is about 1.1%–6.3% while the porosity is between 2.8%–5.5% (Dong DZ et al., 2014, 2018; Borjigin et al., 2017).

Figure  1.  Geological backgroup of the study region. The production field (black plus) is located near the Anyue in the central part of the Sichuan basin. The red represents Suining earthquake. The black lines are faults. The major faults are based on Deng et al (2003). LMSFS=Longmenshan fault system, XDF=Xinjing-Deyang fault, LQF=Longquan fault, HYSF=Huayingshan fault. The dominant horizontal stresses are oriented in the NW-SE directions (blue arrows). The thick grey line around the Anyue marks the possible existing fault found by this study. The green dashed line marks a buried fault that may extend northeastward to the Anyue. Blue arrows show the regional maximum compressional stress.

DownLoad: Full-Size Img PowerPoint

Figure 2 shows the receiver array used in the surface-based monitoring project. 1215 one component geophones with 4.5 HZ natural frequency are deployed along 10 branches. The average spacing is about 20 meters. The elevations of geophones range from 300 to 420 meters. The horizontal well is aligned along the nine perforations (shown as black dots). The altitudes on the southern side of the horizontal well are generally higher (more in red) than the northern side (more in blue) and this is also observable on Figure 1.

Figure  2.  Surface-based monitoring system. Black solid dots are perforations, denoted as P1 to P9 from right to left, respectively. The pluses are geophones with color-coded altitudes ranging from about 300 (blue) to 420 (red) meters.

DownLoad: Full-Size Img PowerPoint

The perforation shots are calibrated to obtain an average velocity model to locate events. The average velocity of v_P=4770 m/s was used for the jSSA locating. For the region with 3D model available, travel times can be computed using 3D raytracer first before being used directly in the jSSA method.

In this study, the scanning ranges in Y and Z directions for all stages are [–150, 350] meters and [–2500, –2150] meters, respectively. The scanning range in X direction is [–300, 300] meters around each perforation point. The scanning intervals are 20 meters for all three directions. Scanning ranges for focal mechanism parameters are [0, 360] degrees for strike, [0, 90] degrees for dip and [0, 180] degrees for rake. The angular interval of 20 degrees is used.

Microseismic events for stages 1–9 are shown in Figure 3 as beachballs. For the stage 1 (Figure 3a), a NW-SE striking linear seismicity passing the perforation 1 is observed. And the seismicity is much denser to the northwest than the southeast of the horizontal well. We interpret this as a fracture generated in the first stage and it is named Frac1. Fractures with similar strike are also observed to pass the perforations 2 and 3 (Frac2 and Frac3 on Figures 3b and 3c, respectively). A close look will find that the Frac1 is still active during the fracturing treatment of the stage 2, while the Frac2 is still active at the fracturing time of the stage 3. The fracture associated with the stages 4 and 5 (Frac4 on Figures 3d, 3e) doesn’t pass through the perforations 4 or 5. Instead it passes through the middle of the perforations 4 and 5. Interestingly, a thin linear seismicity along the Frac4 is also observed in stage 3 (Figure 3c) when this region hasn’t been fractured yet. It also exists in the stage 6 (Figure 3f). This may indicate that the Frac4 is likely to be an existing local fracture that has been activated in the stages 3–6. For the stage 6, a new fracture passing the perforation 6 (Frac5 on Figure 3f) is formed. For stages 7–9, all events are clustered around the horizontal well without any clear linear pattern.

Figure  3.  Locations and focal mechanisms (beachballs) of microseismic events. (a)–(i) are seismic events for the stages 1–9 respectively. Red stars are perforations. The green solid lines show the major fractures created by the hydraulic fracturing treatment while the red dashed lines show the events may be associated with early stages.

DownLoad: Full-Size Img PowerPoint

All 3144 events are replotted in Figure 4a to show the overall pattern of the seismicity. For clarity, Figure 4b shows 1541 large events with R>3.5 (black beachballs) and 189 “explosive” sources (green dots). Based on the stress regime characterization criteria used by Zoback (1992), the faulting types for majority of events are strike-slips with only a small portion of events in normal or imploding/explosive types.

Figure  4.  Locations and focal mechanisms (beachballs) of microseismic events. (a) All events for R>3.0. Open stars are perforations. The open circles are explosive/implosive sources. Black lines show the major fractures created by the hydraulic fracturing treatment. Beachballs in black, magenta, blue, red, green, orange, purple, yellow, and cyan are events occurred during stages 1–9 (count from right), respectively. (b) Events with R≥3.5. The vertical profiles along the dashed lines are displayed in Figure 10.

DownLoad: Full-Size Img PowerPoint

All five fractures are aligned in NW-SE direction. Overall, the fractures are better defined (more events) to the northern side than southern side of the horizontal well. The events to the southern side are much more scattered and sparse. These contrasts are more prominently illustrated on Figure 4b. Events are distributed much denser at several spots near the horizontal well. The last three stages are dominated by clusters of events that show little linear feature. The events appear to be limited within 100 meters on both sides of the perforations for the last three stages while the events of the first six stages extend further northward and southward. The numbers of events increase significantly from stage 6.

Again, it’s even more obvious to find, for the first six stages, events occurred during later stages (defined by time) are often found in the fractures created in earlier stages. For example, events occurred during stage 2 (magenta) are also found along the Frc1; events occurred during stage 3 (blue) are also found along the Frc2–3; events occurred during stage 4 (red) are mostly confined to the Frc2–4; events occurred during stage 5 (green) are mostly confined to Frc4; events occurred during stage six (orange) are widely distributed from Frc3–5. From stage 7 to 9 (purple, yellow and cyan), events are distributed all over between X=300 m to X=900 m.

Another point to make is that the explosive sources are mostly found during the last four stages (green dots on Figure 4b). Caution should be taken though, for a specific focal mechanism, if all or majority of receivers are located in one quadrant of the usually four quadrants of the radiation pattern of P waves, no polarity corrections are involved to yield the strongest stacking energy (Equation 1), therefore it is falsely regarded as an explosive source by the program. Given such a large-aperture array, the possibility for such cases should be small though. Another possibility to see “explosive/implosive” source is that the isotropic components for some events are more dominant than the shear components.

Yu YY et al. (2018) have done extensive tests to show the accuracy of the focal mechanisms found by the jSSA method. Figures 5 and 6 present an example to show the reliability of the polarities predicted by the inverted focal mechanism. Figure 5 shows a focal mechanism with predicted polarities for all receivers. P wave polarities are positive for receivers along the branches B1, B9 and B10, but are negative or mostly negative for receivers along the branches B4–8. The rear parts of the branches B2 and B3 are negative, which are opposite to the other receivers on the same branch. After using the predicted polarities shown in the Figure 5 to correct the waveform for each trace, i.e, multiplying each trace on Figure 6a by 1 (positive) or –1 (negative), Figure 6b shows that the first arrivals are all positive after the corrections.

Figure  5.  P wave polarity of receivers based on the inverted focal mechanism. The red and blue symbols represent positive and negative polarities, respectively. On the background, the grey triangles and reversed black triangles are positive and negative polarities, respectively. The black lines are the projected fault or auxiliary plane. The green triangle marks the epicenter. The 10 branches are labeled as B1–B10, respectively.

DownLoad: Full-Size Img PowerPoint

Figure  6.  Seismograms of one seismic event before and after polarity corrections. (a) Original waveform and (b) waveforms with polarity corrected by the inverted focal mechanism shown in Figure 5.

DownLoad: Full-Size Img PowerPoint

Because focal mechanisms alone cannot distinguish between primary faulting planes and auxiliary planes, it is not straightforward to use strikes, dips and rakes to assess the fracture geometry. However we can compute three principal stresses (P-, T- and B-axes). Yu YY et al. (2018) had done extensive testing using both synthetic and real data to estimate the angular errors of P-axes extracted from focal mechanisms yielded from jSSA. They found that for sources located at a depth of 2000 meters, the angular errors are less than 10° for about 80% of the events.

3.2   Stress field obtained by inversion

The maximum compressional stresses computed directly from focal mechanisms often contains too much small scale (or high frequency) local stress perturbations that are too complex to interpret. Another problem is that every P-axis is extracted independently, and therefore errors in data and errors due to the undesirable receiver coverage may result in unexpected results for individual P-axis. Very often, a stress map with uniform gridding and smooth variation are more useful for engineering purpose. Therefore, a stress inversion method (Hardebeck and Michael, 2006; Martínez-Garzón et al., 2014) is used to extract a uniform 2D stress field.

A damping scheme is adapted to limit the influences of uncertainties in focal mechanisms and to make stress tensor distribution more smoothing. But the damping always tradeoffs with data misfit as shown in Figure 7. To balance the smoothness of the results and the data misfit, the damping number of 0.8 at the corner of the curve is chosen.

Figure  7.  The tradeoff between model uncertainty and data misfit. Each circle is corresponding to a damping factor (the number next to each circle). The plus marks the damping factor chosen in this study.

DownLoad: Full-Size Img PowerPoint

Figures 8a and 8b are the inverted maximum compressional stresses and the angular uncertainty estimated from bootstrapping algorithm, respectively. In fact, majority of the angular uncertainty are less than 30°.

Figure  8.  Inverted maximum compressional stresses and uncertainties. (a) The inverted stress field; (b) the angular uncertainties estimated by 1000 times of boot strapping tests. The region is divided into grids with uniform spacing of 50 m in both X and Y directions. Solid dots are perforations for stages 1–9. The dashed black line is the interpreted existing fault that intersects the horizontal well.

DownLoad: Full-Size Img PowerPoint

The statistics of azimuths and inclinations of the inverted stresses are shown on Figure 9. The stress field is much smooth and simple. The maximum compressional stresses to the north side of the horizontal well appear to be striking systematically NW-SE direction. The fractures are better developed to the northern side of the horizontal well, and strikes of the fractures Frac1–5 (Figures 3 and 4) are consistent with the inverted maximum compressional stress. Stress field near and around the horizontal well is much more complex with prominent NE-SW striking maximum compressional stress and this is consistent with the regional dominant maximum compressional stresses striking in 70°. The maximum compressional stress axes to the south side of the horizontal well are dominantly in north-south direction. The inclinations are mostly less than 15°.

Figure  9.  Rose diagram of azimuths and inclinations of inverted maximum compressional stresses. (a) The azimuths of the maximum compressional stresses; (b) the inclinations of the maximum compressional stresses.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Using the jSSA method, we are able to locate 3144 events in the studying area. The event number is almost doubled compared to the traditional SSA method. The increased event number with extra information of focal mechanisms make it possible to conduct a comprehensive and high resolution stress analysis in the fractured area.

Majority of events are strike-slip faulting. But there are 189 events are identified as “explosive/imploding” sources. However if the polarities of the first arrivals recorded on the array are mostly located in one quadrant of the four quadrants of a complete radiation pattern, it is also falsely recognized by the program as an “explosive/imploding” source. However, the jSSA method is not able to distinguish between the true and false “explosive/imploding” source.

Based on seismic locations, focal mechanisms and P-axes (Figures 3-4), five fractures created by the treatment can be identified. The strikes of the fractures are all in NW-SE direction and are consistent with the inverted stress filed shown in Figure 8a. Overall, the fractures are well developed to the northern side of the horizontal well. The events to the southern side are more scattered with less large events with R≥3.5 (Figure 4b).

Figure 10 shows the vertical projection of seismic events along the Frc3 and Frc4 in Figure 4. The most prominent feature is that a group of large events form a dipping plane extending from –2150 m to –2500 m (limited by the depth range of scanning). To its south side, events are generally small and are mainly confined at the depth of the perforations (–2300 m); to its northwest side, the events are more abundant in all depths. Vertical profiles along other lines show similar features.

Figure  10.  The vertical projection of seismic events along the Frc4 (a) and Frc3 (b) in Figure 3. Horizontal axis shows distances along the line starting from the intersection of the Frc3 or Frc4 and the horizontal axis in Figure 4. The solid dots are the projections of perforations P3–5, as denoted by texts beside the dots. Only events with normal distances to the surface traces of the two fractures less than 40 meters are plotted. The thick dashed lines are the projected existing fault on two planes. The sizes of the circles are proportional to the R values of the events.

DownLoad: Full-Size Img PowerPoint

Together with other evidences, a schematic interpretation shown in Figure 11 is proposed. We suggest that an existing fault may intersect the horizontal well with a small angle. These evidences include: (1) large events are clustered around the horizontal well from west to east (Figure 4b); (2) the group of P-axes in SWW-NEE directions are largely distributed in the vicinity of the horizontal well; (3) stress field near and around the horizontal well has a prominent group of NEE-SWW striking maximum compressional stresses (Figure 8).

Figure  11.  The interpretation map of the study region. The Block I (light green) was separated from the Bock II (light blue) by a NW dipping fault. The Block II is cut to the depth of the horizontal well (thick dashed line). The Block I is more fragile than the Block II to be fractured, therefore the microseismic events are more abundant in Block I. The maximum compressional stress (MCS) are slightly different in the two blocks.

DownLoad: Full-Size Img PowerPoint

Based on these observations, we conclude that an existing fault zone dips northwestwardly with a high angle (~60°) and intersects the horizontal well with a small angle. This fault is probably the northeastward extension of the buried fault to the southwest of the region (dashed green line on the Figure 1). This fault may also serve as a structural boundary. The structure to NW is more fragile to be fractured to yield more events, contrary to its SE side, where less events are located and events are mostly confined in the depth of horizontal well (Figure 10). The maximum compressional stress field to the NW and SE of the fault strike dominantly in NNW-SSE and N-S directions, respectively. The altitudes of receivers on Figures 1 and 2 also show a topographical contrast between the SE (relatively high) and NW.

Because the five fractures (Figure 4) created by hydro-fracturing treatments intersect with the existing fault with an acute angle, the fluids may easily flow along the fault to the fractures created in the earlier stages. Therefore, the events created in the later stages are often found in area fractured in earlier stages (Figures 3b, 3c, 3d, 3f). For the last three stages, the existing fault zone is probably wider to produce large cluster around the well and the events are much larger (Figures 4a and 4b), but the fractured area is limited.

Langenbruch and Shapiro (2015) suggested that elastic-rock heterogeneity controls the occurrence of fluid-injection-induced seismicity. The seismic events occurred preferentially in rock sections characterized by low Poisson’s ratio and high Young’s modulus. Microseismic activity in this study may indicate different rock properties on different sides of the existing fault. Unfortunately, we do not have detailed information to confirm this argument.

Keep in mind though, the working area in this study is only 1600 m by 600 m in scale, therefore the stress field determined is very local compared to the regional stress field. So it is reasonable to see large discrepancy between the stress field determined in this study and the regional stress as shown on Figure 1. However, availability of stress maps determined from other techniques, such as wellbore break outs etc, would be highly preferred to check the reliability of this technique. Nonetheless, the method used here at least presents a new way to extract stress field as a substitute to other more expensive technologies.

Be aware that the current version of jSSA solves the focal mechanisms as double couples (DC). However, tensile components are also observed in hydraulic fracturing and fluid injection experiments (Zoback, 2007; Rodriguez et al., 2012; Grechka et al., 2016) even though the DC components are often dominant for most of the cases. In addition, to solve the tensile component, a fourth parameter (slope) will need to be scanned in the grid search scheme and it will be even more time-consuming than the current version that deals with only strikes, dips and rakes. Therefore, we would like to leave this to future studies.

Given the array to monitor tectonic faults becoming denser and denser, the method presented here has the potential to extract detailed stress distributions associated with a fault system and further to shed light on the dynamic process of an earthquake in a micro scale.

5.   Conclusions

We present a uniform stress field inverted from the microseismics obtained by the joint scanning of source location and focal mechanisms, i.e, the jSSA method. The accuracy and efficiency of the method has been fully investigated in details in early publication. In this study, we take one step further to invert for a uniform stress field based on the events and focal mechanisms yielded from the jSSA method. Based on the above discussions, we conclude: (1) Due to the availability of more events and focal mechanisms, the jSSA method makes it possible to obtain a more complex stress field and fracture system; the method presents a new way to extract stress field as a substitute to other more expensive technologies. (2) Five hydraulic-fracturing induced fractures are identified; the events associated with these fractures are generally smaller in size. (3) One existing fault or weak zone is also identified; events associated with it generally have higher stacked energy and are more densely populated; this fault is probably the extension of the buried fault to the southwest of the shale gas production field. (4) The existing fault may also serve as a structural boundary where the rocks to the NW side are easier to be fractured while events on the other side are sparse with relatively low stacked energy. (5) The existing fault also divides the inverted stress field into two regimes: the maximum compressional stress field to NW and SE strike dominantly in NW-SE and N-S directions, respectively.

Acknowledgement

We are grateful for the two unanimous reviewers for many constructive comments and suggestions that have helped to improve this paper. This work was partially supported by the National Key R&D Program of China (Nos. 2018YFC1503401, and 2016ZX05023‐004), National Natural Science Foundation of China (Nos. 4217040570, and 41674059).