Multi-field coupling numerical simulation on delayed reactivation of hydraulic fracturing induced faults：A case study of induced earthquakes in the Fox Creek area of Canada

The Western Canadian Sedimentary Basin is one of the most active regions in the world for hydraulic fracturing-induced earthquakes （Atkinson et al, 2016; Schultz et al, 2016）. Bao and Eaton （2016） elaborated on the spatiotemporal correlation between hydraulic fracturing operations and seismic activity in the Fox Creek area of the Western Canadian Sedimentary Basin, and found that the largest event （M_W3.9） occurred on a fault that appeared to extend from the injection zone to the crystalline basement is a typically delayed triggering earthquake. Gao et al （2022） believed that the triggering of the forementioned M_W3.9 earthquake was due to the existence of complex fluid migration pathways, along which injected fluids spread accompanied by the Duvernay formation, the eastern fault, and the horizontal pathway of the crystalline basement.

Based on the observed seismic catalog and fault information in the Fox Creek area of Canada, we conduct a numerical simulation study on delayed fault activation. This simulation will combine the tectonic background, earthquake distribution, relevant engineering, and lithological parameters of the Fox Creek area, and based on the hydraulic fracturing principle, fracture seepage theory, fault instability criterion, and fluid-solid coupling theory, it will analyze in detail the mechanism and dynamic process of hydraulic fracturing delayed activation fault.

Firstly, the PKN fracture extension model is used to calculate the stress perturbation input term of the injected fluid, identify the specific location of the fault based on the seismic data, and establish a 2D geological model by combining the stratigraphic and tectonic information.

Then, a numerical simulation model of delayed fault activation in porous elastic media is constructed by coupling solid mechanics, fluid seepage law, and fault activation theory.

Finally, the full process of hydraulic fracturing-induced fault activation was numerically simulated using the finite element method, and the evolution characteristics of the fluid-solid coupling field and stress-strain field before and after fault activation were observed by calculating the value of Coulomb stress change （∆CFS）. The simulation was carried out according to the actual construction plan of the fracturing well section and fault activation, and was divided into the following three main stages: the first stage was the water injection stage, where the well section closest to the target fault was fractured and was injected with water for 5 days; the second stage was the stop injection stage, during which the fluid continued to spread and lasted for 15 days; the third stage was the fault activation and subsequent stage, and the activation time was designed to be at the end of the second stage. The fault continued to be calculated for 20 days after activation, and the total simulation time was 40 days.

The results showed that after 5 days of water injection and 15 days of diffusion, the ∆CFS near the western fault continued to increase, verifying that the delayed activation of the western fault was mainly due to the continuous diffusion and accumulation of fluid to the crystalline basement, resulting in changes in the stress state. Based on the calculation results of the PKN model and the actual fracturing parameters （Bao, Eaton, 2016）, the average fluid injection volume of the Duvernay shale layer was 2 015 m³/d, and after conversion, the injection rate was about 23.32 kg/s. We sliced the pore pressure, Y-direction displacement, and ∆CFS at six key time points during the entire simulation process, as shown in Fig.1. Once fluid injection began （Day 1, Fig. 1a, 1g, 1m）, high-intensity pore pressure diffusion occurred near the reservoir injection points, as well as upward and downward Y-direction displacement and corresponding ∆CFS. After 5 days of injection, the pore pressure had spread to deeper areas, and there was obvious fluid accumulation below the surrounding rock layer and above the crystalline basement boundary （Day 5, light-colored area in Fig. 1b）. Corresponding Y-direction displacement and ∆CFS also further increased （Fig. 1h, 1n）. Moreover, due to the downward diffusion of the fluid, there was an obvious high-value area of pore pressure near the lower part and endpoint of the fault inserted into the crystalline basement （Fig.1n）. At this time, the fault in the crystalline basement already showed had an activation trend, but the shear stress on the fault plane had not reached the critical value for sliding and was in a state of creeping or slow sliding. Although the injection source was lost, the fluid still diffused mainly downward under the action of gravity. Because the permeability of the crystalline basement was lower, the fluid accumulated at the upper boundary of the crystalline basement （Fig.1c, 1d）. However, the fault extending to the basement had a relatively large permeability, so the fluid continued to inject into the lower end of the fault, and the corresponding ∆CFS gradually increased （Fig.1o, 1p）. When ∆CFS reached the critical value, the fault was activated. The continuous diffusion of fluid and accumulation along the lower end of the fault described above was the key reason for the delayed activation of the western fault.

In the case when only the presence of the western fault is considered, after 5 days of water injection and 15 days of fluid diffusion, a high ∆CFS value area with 2.56 MPa appeared at the endpoint of the western fault, and a high ∆CFS value area with 1.62 MPa appeared near the M_W3.9 earthquake rupture point. In the case of two faults coexisting, the above ∆CFS high-value areas increased to 2.98 MPa and 2.11 MPa, respectively. This indicates that the ∆CFS concentration area of the eastern fault in the crystalline basement has made the Coulomb stress of the western fault increase to some extent, leading to the western fault more prone to activation. It should be noted that a ∆CFS high-value connection zone has formed between the eastern and western faults in the crystalline basement, which is the result of mutual stress disturbance between adjacent faults （Fig.1p）. When the stress of a certain fault reaches its critical value, it will be activated first.

Numerical simulation calculation results showed that the Coulomb stress increased area generated by the activation of the western fault controlled the M_W3.9 earthquake and its aftershocks, indicating that the actual spatial distribution of earthquakes is consistent with the fault setting and stress evolution results in the model.

In summary, it is very important to simulate the physical mechanism of water injection-induced earthquakes, and if the forward analysis of the relevant mechanism can be carried out in advance, it will provide a scientific basis for predicting the seismic hazard.