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
-
摘要:
加拿大西部盆地的福克斯溪（Fox Creek）页岩气开采区自压裂开采以来，地震频度急剧增加，引发工业界和科学界的广泛关注。目前一些典型诱发地震案例的断层活化动力学机制尚未完全厘清。本文以2014年初福克斯溪地区Duvernay地层附近发生的地震群及构造为研究对象，开展地下断层受流体扰动而活化的多场耦合数值模拟研究。文中重点讨论的1号平台附近发生的诱发地震主要沿着两个较为清晰的发震断层发生，最大震级MW3.9事件发生于压裂井下方约1 km处的结晶基底内，并具有典型的延迟触发特点。本文就上述断层的滞后活化现象开展数值模拟研究。首先，基于PKN裂缝扩展模型计算并验证注入流体的应力扰动输入项，根据地震数据识别出断层的具体位置，结合地层和构造信息建立二维地质模型；然后，耦合固体力学、流体渗流定律和断层活化理论搭建多孔弹性介质内的断层活化数值仿真模型；最后，采用有限元方法数值模拟水力压裂活化断层的全过程，通过计算库仑应力改变量（ΔCFS）的值来观测断层活化前后的流固耦合场和应力应变场的演化特征。结果表明，经过5天注水和15天扩散的流体作用，西断层附近的ΔCFS持续增加，验证西断层延迟活化的主要成因是流体逐渐扩散累积到结晶基底并改变了其应力状态。此外，模拟结果表明东断层的存在使得西断层更容易活化。断层位错产生的正ΔCFS区域与诱发地震发生位置高度吻合。本文的数值模拟研究重现了水力压裂活化断层的物理过程，相关机制的正演分析若能事先开展，将可能为地震危害性预测提供科学依据。
Abstract: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 （MW3.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 MW3.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 m3/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 MW3.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 MW3.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.
-
-
![]()
图 11 西断层活化后的∆CFS和地震投影
（a） 只有西断层；（b） 东西断层并存
Figure 11. ∆CFS and aftershock projection after activation of the west fault
（a） Case with only the west fault；（b） Case with co-existence of the east and west faults
![]()
Figure 1. Displacement and stress distribution of the numerical simulation process for activation of the west fault while the two faults are co-existed
Figs. （a）−（f） show snapshots of the pore pressure field at the 6 corresponding time points；Figs. （g）−（l） present snapshots of the displacement field in the Y direction；Figs. （m）−（r） display snapshots of the ∆CFS （variation of Coulomb failure stress） field
![]()
图 1 研究区位置（a）和1号平台两口水平井及附近地震事件的分布（Bao，Eaton，2016）（b）
Figure 1. Study area （a） and distribution of earthquakes near the two horizontal wells of Platform 1 （Bao，Eaton，2016）（b）
![]()
图 2 地震丛集时空特征分布图
两个相对独立的丛集地震分别用浅绿色和浅红色椭圆区域表示，最大震级事件的震源机制解绘制在图中（Wang et al，2016）。（a） 丛集地震的水平投影。红色线段代表1号平台的两口水平井，蓝色箭头指示出逐段压裂方向，蓝色倒三角形为压裂施工阶段P2的最后5天对应的压裂段位（Gao et al，2022）；（b） 垂向剖面。红色虚线表示两条估计断层轨迹，断层位置和倾角参考Bao和Eaton （2016）；（c） M-t图。图中标注了压裂施工阶段P1和P2的时间跨度，其中虚线框表示压裂施工阶段P2的最后5天，对应于图（a）中的蓝色倒三角形
Figure 2. The spatiotemporal distribution of earthquake clusters
Two relatively independent clusters are represented by light green and light red elliptical regions，respectively. The focal mechanism solution of the largest seismic event is depicted in the figure （Wang et al，2016）. （a） The map view of the cluster earthquakes，with red lines representing two horizontal wells of Platform 1，blue arrows indicating the step-by-step fracturing direction，and a blue inverted triangle representing the fracturing segment corresponding to the last 5 days of fracturing construction stage P2 （Gao et al，2022）. （b） A vertical profile，with red dashed lines representing two estimated fault trajectories.The fault location and dip angle are referenced from Bao and Eaton （2016）. （c） The magnitude-time distribution chart，with fracturing construction stages P1 and P2 marked in the figure. The dashed box represents the last5 days of fracturing construction stage P2，corresponding to the blue inverted triangle in Fig.（a）
![]()
图 3 裂缝扩展PKN模型的模拟计算结果
（a） 裂缝长度和裂缝宽度的扩展演化；（b） 裂缝长度与裂缝高度的关系，红色表示裂缝高度的延伸方向，蓝色表示沿着裂缝长度方向的裂缝面；（c） 对应的裂缝流体压力演化图
Figure 3. Simulation results of fracture propagation in PKN model
（a） The evolution of fracture length and width；（b） The relationship between fracture length and height；（c） The corresponding evolution of fluid pressure within the fractures
![]()
图 4 本文数值模拟的模型示意图及主要模型参数
（a） 福克斯溪地区1号页岩气开发平台的水平井、水力裂缝、页岩储层及主要断层的构造位置示意图；（b） 二维数值模拟的模型构建、尺寸及边界条件；（c） 图（b）中关于压裂注入井与断层的放大表示图
Figure 4. Illustration and key model parameters of the numerical simulation in this study
（a） Schematic diagram of the structural positions of the horizontal wells，hydraulic fractures，shale reservoir，and major faults in the Fox Creek area，specifically highlighting Platform 1； （b） The construction，dimensions，and boundary conditions of the two-dimensional numerical model；（c） An enlarged representation of the hydraulic fracturing injection well and the fault within the diagram shown in Fig.（b）
![]()
图 5 西断层活化数值模拟全过程的位移场及应力场分布
图（a）−（f），（g）−（l），（m）−（r）分别为6个时间点对应的孔隙压力场分布、Y方向位移场及∆CFS应力场的切片
Figure 5. Distribution of displacement field and stress field during the numerical simulation process of activation for the west fault
Figs.（a）−（f），（g）−（l） and （m）−（r） show snapshots of the pore pressure field，the displacement field in the Y direction and the ∆CFS （Coulomb failure stress） at the corresponding six time points
![]()
图 6 西断层活化前后的位移场及断层库仑应力分布
图（a）−（d）为断层活化前后的位移场，图（e）−（h）为断层外部区域的∆CFS应力场切片，图（i）−（l）为断层内部的∆CFS等值线切片。图（g）和（k）标注了断层活化前∆CFS高值区的等值线值
Figure 6. Distribution of displacement field and Coulomb stress of the western fault before and after activation
Figs. （a）−（d） show the displacement fields before and after fault activation；Figs. （e）−（h） show the ∆CFS stress field snapshots in the external area of the fault；Figs. （i）−（l） show the ∆CFS contour snapshots in the internal area of the fault. Figs. （g） and （k） mark the contour value of high ∆CFS value area before fault reactive
![]()
图 7 沿断层上盘的孔隙压力变化曲线
曲线颜色对应于不同采样时刻的断层上盘孔隙压力值，蓝色为断层活化前0—20天，红色为模拟的20—40 天，断层活化后
Figure 7. The curve of pore pressure variation along the upper plate of the fault is shown
The curve colors correspond to the pore pressure values of the hanging wall of the fault at different sampling times，with blue representing 0−20 days before reactivation and red representing 20−40 days after reactivation
![]()
图 8 断层活化前（a）、后（b）沿断层上下盘的∆CFS变化曲线
上盘∆CFS对应浅红色曲线，下盘∆CFS为深灰色曲线，每条曲线代表一个模拟采样间隔单位（d）蓝色条带对应断层两侧端部，绿色条带对应代表断层与结晶基底交界区域
Figure 8. The ∆CFS variation curves before （a） and after （b） fault activation along the fault plane
The shallow red curve corresponds to the ∆CFS of hanging wall，while the deep gray curve corresponds to the∆CFS of footwall. Each curve represents a simulated sampling interval （in days）. The blue bands correspond to the ends of both sides of the fault, and the green bands correspond to the area where the fault meets the crystalline basement
![]()
图 9 东西两条断层共存时西断层活化数值模拟全过程的位移场及应力场分布
图（a）−（f），（g）−（l），（m）−（r）分别为6个时间点对应的孔隙压力场分布、Y方向位移场及∆CFS应力场的切片
Figure 9. Distribution of displacement field and stress field of the numerical simulation process for activation of the west fault while the east and west faults are co-existed
Figs.（a）−（f），（g）−（l） and （m）−（r） show snapshots of the pore pressure field，the displacement field in the Y direction and the ∆CFS （Coulomb failure stress） field at the corresponding six time points
![]()
图 10 断层活化前的∆CFS分布
（a） 只有西断层；（b） 东西断层并存
Figure 10. Distribution of ∆CFS before fault activation
（a） Case with only the west fault；（b） Case with co-existence of the east and west faults
表 1 数值模拟的主要参数
弹性
模量
E/GPa泊松比
v密度ρd
/（kg·m−3）孔隙度
ϕ渗透率
κBiot系数
α结晶基底 66 0.20 2 900 0.2 10−17 0.75 页岩层 25 0.25 2 500 0.2 10−15 0.75 上下围岩层 40 0.25 2 500 0.2 10−14 0.75 断层 14.4 0.20 2 500 0.02 10−12 0.80 
下载: 导出CSV
-
惠钢,陈胜男,顾斐. 2021. 流体-地质力学耦合建模表征水力压裂诱发地震:以加拿大Fox Creek地区为例[J]. 地球物理学报,64(3):864–875. doi: 10.6038/cjg2021O0267 Hui G,Chen S N,Gu F. 2021. Coupled fluid-geomechanics modeling to characterize hydraulic fracturing-induced earthquakes:Case study in Fox Creek,Canada[J]. Chinese Journal of Geophysics,64(3):864–875 (in Chinese).
Atkinson G M,Eaton D W,Ghofrani H,Walker D,Cheadle B,Schultz R,Shcherbakov R,Tiampo K,Gu J,Harrington R M,Liu Y J,Van Der Baan M,Kao H. 2016. Hydraulic fracturing and seismicity in the Western Canada Sedimentary Basin[J]. Seismol Res Lett,87(3):631–647. doi: 10.1785/0220150263
Bao X W,Eaton D W. 2016. Fault activation by hydraulic fracturing in western Canada[J]. Science,354(6318):1406–1409. doi: 10.1126/science.aag2583
Biot M A. 1941. General theory of three‐dimensional consolidation[J]. J Appl Phys,12(2):155–164. doi: 10.1063/1.1712886
Catalli F,Meier M A,Wiemer S. 2013. The role of Coulomb stress changes for injection-induced seismicity:The Basel enhanced geothermal system[J]. Geophys Res Lett,40(1):72–77. doi: 10.1029/2012GL054147
Chang K W,Segall P. 2016. Injection-induced seismicity on basement faults including poroelastic stressing[J]. J Geophys Res:Solid Earth,121(4):2708–2726. doi: 10.1002/2015JB012561
Clarke H,Eisner L,Styles P,Turner P. 2014. Felt seismicity associated with shale gas hydraulic fracturing:The first documented example in Europe[J]. Geophys Res Lett,41(23):8308–8314. doi: 10.1002/2014GL062047
Deng K,Liu Y J,Harrington R M. 2016. Poroelastic stress triggering of the December 2013 Crooked Lake,Alberta,induced seismicity sequence[J]. Geophys Res Lett,43(16):8482–8491. doi: 10.1002/2016GL070421
Ellsworth W L. 2013. Injection-induced earthquakes[J]. Science,341(6142):142.
Galis M,Ampuero J P,Mai P M,Gappa F. 2017. Induced seismicity provides insight into why earthquake ruptures stop[J]. Sci Adv,3(12):eaap7528.
Gao D W,Kao H,Wang B,Visser R,Schultz R,Harrington R M. 2022. Complex 3D migration and delayed triggering of hydraulic fracturing-induced seismicity:A case study near Fox Creek,Alberta[J]. Geophys Res Lett,49(2):e2021GL093979. doi: 10.1029/2021GL093979
Gillian R F,Miles P W,Jon G G,Bruce R J,Richard J D. 2018. Global review of human-induced earthquakes[J]. Earth-Sci Rev,178:438–514
Goebel T H W,Brodsky E E. 2018. The spatial footprint of injection wells in a global compilation of induced earthquake sequences[J]. Science,361(6405):899–904. doi: 10.1126/science.aat5449
Holland A A. 2013. Earthquakes triggered by hydraulic fracturing in south-central Oklahoma[J]. Bull Seismol Soc Am,103(3):1784–1792. doi: 10.1785/0120120109
Hulls C K. 2022. Geomechanical Modeling of a Fault During Fluid Injection[D]. London:The University of Western Ontario.
Lei X L,Huang D J,Su J R,Jiang G M,Wang X L,Wang H,Guo X,Fu H. 2017. Fault reactivation and earthquakes with magnitudes of up to MW4.7 induced by shale-gas hydraulic fracturing in Sichuan Basin,China[J]. Sci Rep,7(1):7971. doi: 10.1038/s41598-017-08557-y
Leonard M. 2010. Earthquake Fault Scaling:Self-Consistent Relating of Rupture Length,Width,Average Displacement,and Moment Release[J]. Bull Seism Soc Am,100(5A):1971–1988
Li L,Tan J Q,Wood D A,Zhao Z G,Becker D,Lyu Q,Shu B,Chen H C. 2019. A review of the current status of induced seismicity monitoring for hydraulic fracturing in unconventional tight oil and gas reservoirs[J]. Fuel,242:195–210. doi: 10.1016/j.fuel.2019.01.026
Meng L Y,McGarr A,Zhou L Q,Zang Y. 2019. An investigation of seismicity induced by hydraulic fracturing in the Sichuan Basin of China based on data from a temporary seismic network[J]. Bull Seismol Soc Am,109(1):348–357. doi: 10.1785/0120180310
Nordgren R P. 1972. Propagation of a vertical hydraulic fracture[J]. Soc Petrol Eng J,12(4):306–314. doi: 10.2118/3009-PA
Okada Y. 1992. Internal deformation due to shear and tensile faults in a half‐space[J]. Bull Seismol Soc Am,82(2):1018–1040. doi: 10.1785/BSSA0820021018
Peña Castro A F,Roth M P,Verdecchia A,Onwuemeka J,Liu Y,Harrington R M,Zhang Y,Kao H. 2020. Stress chatter via fluid flow and fault slip in a hydraulic fracturing-induced earthquake sequence in the Montney Formation,British Columbia[J]. Geophys Res Lett,47(14):e2020GL087254. doi: 10.1029/2020GL087254
Perkins T K,Kein L R. 1961. Widths of hydraulic fractures[J]. J Pet Technol,13(9):937–949. doi: 10.2118/89-PA
Rice J R,Cleary M P. 1976. Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents[J]. Rev Geophys,14(2):227–241. doi: 10.1029/RG014i002p00227
Schultz R,Corlett H,Haug K,Kocon K,MacCormack K,Stern V,Shipman T. 2016. Linking fossil reefs with earthquakes:Geologic insight to where induced seismicity occurs in Alberta[J]. Geophys Res Lett,43(6):2534–2542. doi: 10.1002/2015GL067514
Schultz R,Atkinson G,Eaton D W,Gu Y J,Kao H. 2018. Hydraulic fracturing volume is associated with induced earthquake productivity in the Duvernay play[J]. Science,359(6373):304–308. doi: 10.1126/science.aao0159
Segall P,Lu S. 2015. Injection-induced seismicity:Poroelastic and earthquake nucleation effects[J]. J Geophys Res:Solid Earth,120(7):5082–5103. doi: 10.1002/2015JB012060
Skoumal R J,Brudzinski M R,Currie B S. 2015. Earthquakes induced by hydraulic fracturing in Poland township,Ohio[J]. Bull Seismol Soc Am,105(1):189–197. doi: 10.1785/0120140168
Skoumal R J,Ries R,Brudzinski M R,Barbour A J,Currie B S. 2018. Earthquakes induced by hydraulic fracturing are pervasive in Oklahoma[J]. J Geophys Res:Solid Earth,123(12):10918–10935.
Tan Y Y,Hu J,Zhang H J,Chen Y K,Qian J W,Wang Q F,Zha H S,Tang P,Nie Z. 2020. Hydraulic fracturing induced seismicity in the Southern Sichuan Basin due to fluid diffusion inferred from seismic and injection data analysis[J]. Geophys Res Lett,47(4):e2019GL084885. doi: 10.1029/2019GL084885
Wang B,Harrington R M,Liu Y J,Kao H,Yu H Y. 2020. A study on the largest hydraulic-fracturing-induced earthquake in Canada:Observations and static stress-drop estimation[J]. Bull Seismol Soc Am,110(5):2283–2294. doi: 10.1785/0120190261
Wang B,Verdecchia A,Kao H,Harrington R M,Liu Y J,Yu H Y. 2021. A study on the largest hydraulic fracturing induced earthquake in Canada:Numerical modeling and triggering mechanism[J]. Bull Seismol Soc Am,111(3):1392–1404. doi: 10.1785/0120200251
Wang R J,Gu Y J,Schultz R,Kim A,Atkinson G. 2016. Source analysis of a potential hydraulic-fracturing-induced earthquake near Fox Creek,Alberta[J]. Geophys Res Lett,43(2):564–573. doi: 10.1002/2015GL066917
Wang R J,Gu Y J,Schultz R,Zhang M,Kim A. 2017. Source characteristics and geological implications of the January 2016 induced earthquake swarm near Crooked Lake,Alberta[J]. Geophys J Int,210(2):979–988. doi: 10.1093/gji/ggx204
Wells D L,Coppersmith K J. 1994. New empirical relationships among magnitude,rupture length,rupture width,rupture area,and surface displacement[J]. Bull Seismol Soc Am,84(4):974–1002. doi: 10.1785/BSSA0840040974
Zhang H L,Eaton D W,Li G,Liu Y J,Harrington R M. 2016. Discriminating induced seismicity from natural earthquakes using moment tensors and source spectra[J]. J Geophys Res:Solid Earth,121(2):972–993. doi: 10.1002/2015JB012603
Zhu W Q,Allison K L,Dunham E M,Yang Y Y. 2020. Fault valving and pore pressure evolution in simulations of earthquake sequences and aseismic slip[J]. Nat Commun,11(1):4833. doi: 10.1038/s41467-020-18598-z
计量
- 文章访问数: 268
- HTML全文浏览量: 146
- PDF下载量: 237
