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引言
近年来，与页岩气开发相关的水力压裂诱发地震案例在美国（Ellsworth，2013；Holland，2013；Skoumal et al，2015，2018）、加拿大（Atkinson et al，2016；Schultz et al，2016，2018）、欧洲（Clarke et al，2014；Galis et al，2017；Li et al，2019）、中国（Lei et al，2017；Meng et al，2019；Tan et al，2020）等多个国家和地区被广泛观测到。诱发地震现象最为突出的区域主要分布在北美的加拿大阿尔伯塔（Alberta）省福克斯溪（Fox Creek）地区（Wang et al，2016，2017）及不列颠哥伦比亚（British Columbia）省圣约翰堡（Fort St. John）和道森克里克（Dawson Creek）地区（Peña Castro et al，2020；Wang et al，2020，2021）、美国得克萨斯州（Texas）和美国俄克拉荷马州（Oklahoma）等地区，这些地区的页岩气开发直接造成附近地域地震频度和震级的显著提升，随之而来的地震危害性攀升不可忽视。
Ellsworth （2013）指出流体注入导致的孔隙压力和孔隙弹性应力变化是诱发地震的两种核心机制。流体注入使得断层的孔隙压力增加和有效正应力减小，摩尔圆左移相交于破裂包络线导致断层破裂失稳，是诱发地震的最主要和最常见的机制（Atkinson et al，2016；Schultz et al，2016）。惠钢等（2021）基于流固耦合理论对福克斯溪地区2015年2月8日发生的ML3.0地震事件展开数值模拟研究，结果表明压裂后断层周围的孔隙压力和局部应力场发生了显著变化，直接导致其活化，而此处孔隙压力的增加是断层活化的首要因素。另外，将作用到储层岩体的孔隙弹性应力耦合到流体计算中也是主要的研究方法。Deng等（2016）基于多孔线弹性理论和流固耦合模型计算了福克斯溪地区水力压裂作业引起的应力场演化过程。Segall和Lu （2015）及Chang和Segall （2016）采用类似方法建立多场耦合模型，研究了废水注入对不同类型断层的孔隙弹性应力扰动。本文拟借鉴上述经典模型及方法开展多孔弹性介质中的多场耦合数值模拟研究.
探索诱发地震的致灾机理是地质灾害预测与防治的重要基础工作，而断层受扰活化的动力学数值模拟则是重要的手段和方法。当受扰断层的剪应力小于接触面的摩擦力时，断层将保持静止状态；当二者相等时，断层将出现准动态的滑动趋势；当断层面某处的剪应力超过静摩擦力后，断层将产生局部破裂和活化，进而可能发生规模性滑动，导致诱发地震（Ellsworth，2013）。
加拿大西部盆地是全世界水力压裂诱发地震研究最为活跃的区域之一（Atkinson et al，2016；Schultz et al，2016），其中Duvernay地层曾发生水力压裂诱发地震的较大震级事件（Wang et al，2016，2017）。加拿大不列颠哥伦比亚地区的几次水力压裂诱发地震事件也引起了广泛关注。Wang等（2020）计算了2015年8月17日发生的MW4.6地震的应力降并对最大可能震级进行了详细分析，认为该次地震发生在一条先存断层上，其最大震级受控于地质构造条件。其后续数值模拟结果表明压裂注液区到目标断层之间很可能存在一条高渗透率的流体通道（Wang et al，2021），该结论与同一地区Montney地层的MW4.2水力压裂诱发地震类似（Peña Castro et al，2020）。Bao和Eaton （2016）详细论述了加西盆地福克斯溪地区6个压裂平台的注水作业与地震活动的时空关联性，发现震级最大事件（MW3.9）发生在一条看似从注水区域延伸到结晶基底的断层上。该断层附近陆续发生的地震事件表明其可能被工业注水持续激活，但是其具体活化过程并未在已有研究中进行严格的正演分析。Gao等（2022）认为上述MW3.9地震的触发是因为存在复杂的流体运移通道，首先是注入的液体沿Duvernay地层内部的水平向通道向东平移，到达东侧断层后使其发生无震滑移并触发较深部的一系列地震，精定位后的地震序列更是证明东断层到西断层之间在结晶基底附近可能存在一条近水平的流体通道，该通道促发了西断层末端MW3.9地震的发生。上述分析为复杂构造下的诱发地震活动成因提供了全新视角，但并未得到地震成像或应力计算的佐证。然而，较深部地震的延迟发生可能是压裂停止后流体缓慢扩散所直接导致。
本文拟基于加拿大福克斯溪地区已经观测到的地震目录和断层信息，拟开展断层受扰延迟活化的流固耦合数值模拟研究。模拟中将发震构造简化为两条主要断层，忽略已有研究中提出的含流体通道的复杂构造环境（Gao et al，2022）。在遵循实际施工进程和参数的前提下，重点分析页岩层内的压裂注水对储层下方主要断层的延迟活化过程。本模拟将结合福克斯溪地区的区域构造背景、地震分布、压裂注水工程和岩性参数，根据水力压裂原理、裂隙渗流理论、断层失稳准则和流固耦合理论详细分析水力压裂延迟活化断层的机制和动力学过程。
1. 研究区域与数据
1.1 研究区域及井位介绍
本文的研究区域位于加拿大阿尔伯塔省的福克斯溪地区（图1），在该地区曾观测到众多与页岩气开发有关的诱发地震（Schultz et al，2016）。重点关注的Duvernay地层在2014年12月至2015年4月水力压裂作业期间于6个不同井位发生了时空关系迥异的诱发地震序列，而2015年1月23日更是在1号平台附近发生了一次走滑-逆冲机制的MW3.9地震。诱发地震的丛集展布方位及震源机制解均为本文断层几何形态的模拟提供了依据（Bao，Eaton，2016；Hulls，2022）。
![]()
Figure 1. Displacement and stress distribution of the numerical simulation process for activation of the west fault while the two faults are co-existedFigs. （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上述MW3.9地震及其前震和余震均发生在1号水平井附近，地震数据剖面图显示出两条明显的、急剧倾斜的地震活跃带，一直从Duvernay地层内的注入带延伸到结晶基底的上部（Bao，Eaton，2016；Wang et al，2016）。我们重新梳理了这些地震和井轨迹，如图2所示。基于时空关系聚合效应，能够判断出两个相对独立的丛集地震。浅绿色区域内的丛集地震主要发生于2014年12月23日至2015年1月15日之间，最大震级为发生在结晶基底之内的MW3.2地震，丛集地震的分布显示了其东侧断层的大致形态，估计倾角为75°—85°。浅红色区域内的丛集地震陆续发生于2015年1月10日至3月30日之间，最大震级为1月23日发生于结晶基底的MW3.9地震，丛集地震的分布显示其东侧断层的倾角约为80°—85°，地震丛集持续时间长、震级大，其时空关系存在显著的前震余震特征。MW3.9地震的破裂位于压裂井正下方、整个压裂段的末端，但是在压裂结束后15天才发生，伴随一定的前震和尾部较长的余震（图2c中粉红色和红色圆圈）。值得注意的是，图2a中红色线段表示的两口水平井均为自北向南实施压裂，北段压裂时触发绿色区域的丛集地震，但同时南部近井也有小震发生。考虑到1号平台施工过程中93%的压裂液未能返排，其附近断裂带可能发生了持续的增压和激活（Bao，Eaton，2016）。最大震级MW3.9地震位于前寒武纪结晶基底上部的深处，随后的序列则向较浅处靠近注入带的位置迁移。
![]()
图 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. 研究方法与原理
本文的主要研究内容是数值模拟Duvernay地层内对应于2015年1月23日Mw3.9地震的西侧断层在压裂流体扰动下的活化过程，主要涉及水力裂缝扩展及其与断层的流固耦合、断层活化等关键步骤。本节将阐述数值模拟与模型构建的方法原理及实现路径。
2.1 人工压裂缝扩展模拟
水力压裂裂缝的尺度、扩展及闭合都对工程效能具有重大影响，本节首先使用经典的PKN模型来描述水力裂缝的扩展过程。PKN模型最初由Perkins和Kein提出（1961），后因Nordgren （1972）添加滤失项而改进，并因此得名。该模型基于垂直平面应变理论建立线弹性裂缝力学方程、缝内液体流动方程及连续方程，并在给定的边界条件下进行求解，是一个椭圆形等高裂缝模型，其表达式为：
$$ \begin{array}{c}L （ t ） ＝0.654{\left[\dfrac{G{q}_{0}^{3}}{ （ 1－\nu ） \mu {h}^{4}}\right]}^{1/5}{t}^{\,{4}/{5}} \text{，} \\ W （ t ） ＝2.5{\left[\dfrac{ （ 1－\nu ） \mu {q}_{0}^{2}}{Gh}\right]}^{1/5}{t}^{1/5} \text{，} \\ P （ t ） ＝2.5{\left[\dfrac{{G}^{6}\mu {q}_{0}^{2}}{ （ 1－\nu { ） }^{4}{h}^{6}}\right]}^{1/5}{t}^{1/5}\text{，} \end{array} $$ (1) 式中，$L （ t ） $为t时刻裂缝扩展长度，$ W （ t ） $为t时刻裂缝扩展宽度，$ P （ t ） $为t时刻流体压力，G为剪切模量，q0为注入速率，μ为注入流体的黏度，h为压裂缝缝高，v为泊松比。采用Python语言对式（1）进行模拟和可视化，得到裂缝扩展的缝长、缝宽和压力等关系（图3）。作为示例，这里实施的模拟计算时间为10 min。模拟计算采用加拿大Duvernay地层实际压裂注水参数作为输入（2 015 m3/d），裂缝的长宽高取值由PKN模型结果表达式计算求得，流体注入稳定后其孔隙压力均值约为（50±15） MPa （图3c）。据此，我们将该值作为后续多场耦合模拟的流体输入关键参数。
![]()
图 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 clustersTwo 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）2.2 多孔弹性介质流固耦合模拟理论
多孔弹性介质常用于研究流体和断层的交互力学作用（Segall，Lu，2015；Bao，Eaton，2016；Deng et al，2016；Goebel，Brodsky，2018）。假设压裂流体在介质中服从达西渗流定律，基于Biot方程描述固体孔隙弹性应力及流固耦合的质量守恒和力学平衡（Biot，1941；Rice，Cleary，1976），有如下多孔弹性介质流固耦合数学模型：
$$ \nabla \sigma ＋ {\rho }_{{{\mathrm{b}}}}g＝0 \text{，} $$ (2) $$ {\rho }_{{\mathrm{f}}}S\frac{ {\text{∂}} P}{ {\text{∂}} t} ＋ {\rho }_{{\mathrm{f}}}\alpha \frac{ {\text{∂}} {\varepsilon }_{v}}{ {\text{∂}} t} ＋ \nabla \cdot \left[-{\rho }_{\mathrm{f}}\frac{k}{{\eta }_{\mathrm{f}}} （ \nabla P-{\rho }_{\mathrm{f}}{{{g}}} ） \right]＝0 \text{，} $$ (3) $$ \varepsilon ＝\frac{1}{2} （ \nabla {\boldsymbol{u}} ＋ \nabla {{\boldsymbol{u}}}^{{\mathrm{T}}} ） \text{，} $$ (4) 式中，P为流体压力，ρf为流体密度，S为流体存储系数，k为渗透率，α为Biot系数，ηf为流体动态黏度，ρb为岩石体积密度，$\nabla {\boldsymbol{u}} $为位移场，ε为体积应变，g为重力加速度。
2.3 库仑应力计算
通过模拟断裂带周围孔隙压力和地应力的交互变化（Catalli et al，2013），基于摩尔-库仑准则计算断层附近的库仑应力改变量（Okada，1992）：
$$ \mathrm{\Delta }\mathrm{C}\mathrm{F}\mathrm{S}＝\Delta \tau －\mu （ \Delta {\sigma }_{{\mathrm{n}}}－\Delta P ） \text{，} $$ (5) 式中，∆P为孔隙压力变化量，∆τ为断层上的剪切应力变化量，∆σn为有效正应力变化量，∆CFS为库仑应力改变量，该值为正则表明断层更趋向于失稳。
2.4 数值计算步骤、模型构建及相关参数
本研究的模拟按照实际压裂井段的施工计划和断层活化进行计算，主要分为以下三个阶段：第一阶段为注水阶段，向距离目标断层物理位置相对接近的井段压裂注水，总注水时间折合为5天；第二阶段为停止注水阶段，此时流体将继续扩散，持续15天；第三阶段为断层活化及后续阶段，活化起始时间设为第二阶段结束时刻，断层活化后继续计算20天，总模拟时间共计40天。
基于Duvernay地层的实际构造、压裂参数和地震活动数据（图2），本文按照以下几何尺寸设置模型：边长为10 km，模型中包含了Duvernay页岩层、上下围岩、结晶基底等三类地层；西断层长600 m，厚10 m，东断层长400 m，厚10 m，东西断层倾角均约为83°；注入井位于西断层上方的储层内部，垂直展布于平面向内方向（图4）。采用PKN模型计算出的孔隙压力作为流体注入扰动。模型左、右、底边界都采用辊支撑（即第二类/自由边界条件），上界为自由表面、边界孔隙压力为0，其它边界无流体流动，在断层上下盘两侧分别设置监控线来观测相关参数在模拟过程中的变化，监测线与断层上下盘的距离接近于0。需要说明的是，本模型用二维垂直截面来清晰表达推测断层的线性特征（图2b），断层倾角的设定则综合考虑了地震丛集的分布特征和MW3.9最大震级事件的震源机制解。Zhang等（2016）计算得出MW3.9地震的震源特征为走滑逆冲机制类型，其滑动面（节面Ⅰ）的反演结果为走向176°，倾角83°，滑动角135°，滑动角135°表示其走滑分量和逆冲分量大约各占50%。二维截面模拟主要考虑YOZ剖面的逆冲机制部分（图2），因此无法体现垂直于模型平面向内的走滑分量，后续将走滑逆冲机制断层简称为“逆断层”。本文假设模拟平面（即YOZ面）内的断层尺度为“长度”，而垂直于YOZ面向内方向的断层尺度为“宽度”。
![]()
图 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下面阐述本文模型及模拟设定的断层物理尺寸和位错距离的理论计算依据。Leonard （2010）在Wells和Coppersmith （1994）提出的“震级-破裂尺度”经验公式基础上，凭借更完备的数据集提出了隐伏逆断层的长度L与矩震级MW的经验关系式：
$$ {M}_{\mathrm{W}}＝1.67\mathrm{l}\mathrm{g}L ＋ 4.24 \text{.} $$ (6) 通过式（6）估算MW3.9地震所对应的断层长度约为625.8 m，证明本模型设置的断层近似长度600 m具有一定合理性。
进一步采用Leonard （2010）提出的中小震级断层位错关系分别计算上下盘的相对滑动距离。对于震级小于MW5.0的地震，地震矩与断层长度符合M0∝GL3，断层位错距离$ \overline {D} $为
$$ \overline {D} ＝C{L}^{2－\beta } \text{，} $$ (7) 式中：L为断层破裂长度，单位为m；C为经验常量，震级较小的逆断层一般取值为3.7×10−5 （Leonard，2010）。根据式（7），对于长度为600 m的断层，发生MW3.9地震所产生的滑动位移约为0.029 6 m。由于震源机制解中逆冲机制分量约占50%，因此在模拟中施加的断层位错距离为总滑移距离的50%，为0.014 8 m。数值模拟其它相关参数详见表1。
表 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 3. 模拟结果及分析
本文的数值模拟工作分为两部分，首先分析Duvernay页岩层下方西断层独立存在时的活化过程，然后进一步研究东西两条断层同时存在时，西断层的活化过程，并讨论东断层对西断层活化的影响。数值计算主要分为三个阶段：一是压裂流体持续注入5天；二是停止注入15天，实际数据表明此时断层附近的∆CFS将达到临界值，因此模拟中断层活化的开启阀门受此时间控制；三是断层活化并进行有限滑动，直至20天后模拟结束。本节详细描述相关模拟结果。
3.1 独立西断层活化过程模拟
由于西断层位于Duvernay页岩层压裂注入点的正下方，并且延迟发生了一次位于结晶基底内的MW3.9最大震级事件，因此其活化过程可能与压裂流体的注入和持续渗透直接相关。
基于PKN模型计算结果和实际压裂参数（Bao，Eaton，2016），Duvernay页岩层的平均流体注入量为2 015 m3/d，经过换算，注入速率约为23.32 kg/s。本文使用该平均速率对两个压裂井进行源头设置，持续输入到附近地层之中。根据实际水平井压注施工日程和断层被激活的流体注入范围，将注水总时长设定为5天，停止注液15天后断层滑移，并继续计算观测20天。我们对整个模拟过程中的6个关键时间点的孔隙压力、Y方向位移和∆CFS进行了切片展示（图5）。流体注入一经开始（第1天，图5a，5g，5m），储层注液点附近就产生了高强度孔隙压力扩散、向上和向下的Y方向位移以及对应的∆CFS。上述变化发生在压裂注水点附近，最初注液阶段对下部断层的影响并不明显。注水持续5天后停止，此时孔隙压力已经扩散至更深部区域，围岩层下方与结晶基底边界上方出现明显的流体聚集（第5天，图5b的浅色区域）。对应的Y方向位移和∆CFS也出现进一步增强（图5h，5n）。并且，由于流体向下扩散，插入结晶基底的断层下部及端点附近出现明显的高值区（图5n）。此时结晶基底内的断层部分已经存在活化趋势，但断层面的剪切应力并未达到滑移临界值，处于蠕滑或者慢滑移状态。
![]()
图 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天后则进入模拟的第二阶段，此时停止注液，持续15天。由于流体注入源头消失，水平井眼附近的孔隙压力因失去了持续增加的动力而逐渐减弱。但是，流体依然在重力作用下发生向下为主的扩散。由于结晶基底的渗透率更低，使得流体聚集到结晶基底上边界处（图5c，5d）。然而，延伸到基底的断层则具有相对较大的渗透率，因此流体持续灌注到断层下端，这一现象能够被清晰地观测到，且在本阶段最后一天（第20天）尤为明显（图5d）。上述流体渗透特性使得处于结晶基底内的断层内的孔隙压力逐渐增大，有效正应力逐渐减小，根据式（5），其对应的∆CFS也逐渐增大（图5o，5p）。当达到临界值后，断层则发生活化。上述流体的持续扩散及沿断层下端聚集，正是西断层延迟活化的关键原因，本模拟过程充分展示了这一点。
当模拟进行到第20天时，断层下端的∆CFS达到临界值，设定的断层阀门打开，此时断层打破之前的蠕滑状态，进入破裂和快速滑移状态（Zhu et al，2020）。模拟中假设断层的总滑动量和滑动方向受限于MW3.9最大震级事件，符合逆断层动力学特征。
由于逆断层活化，断层附近的位移方向发生了显著变化，图5k和5q分别展示了断层活化后Y方向位移变化与∆CFS的瞬时变化。断层活化后，随着断层滑动的停止和愈合，其附近的位移场和应力场基本维持既有分布和大小，趋于稳定（图5l，5r）。
为了更好地呈现断层活化前后其内部及附近区域的∆CFS和上下盘的位移情况，我们截取了相关模拟结果的切片继续进行展示（图6）。显然，断层活化前的位移场主要体现了压裂流体的注入和扩散对附近地层的形变影响（图6a−c），而活化后断层上下盘则体现出清晰的逆断层位移特征，即上盘向上移动、下盘向下移动（图6d）。此时，该逆断层的活化也导致∆CFS正值集中在断层附近（图6h），∆CFS的分布均表现出与断层活化前的注液阶段（图6e，6f）和停注扩散阶段（图6g）截然不同的辐射花样。位于结晶基底的断层内部的∆CFS也表现出对应的变化特征，从注液开始（图6i）到停止（图6j），再到流体扩散聚集（图6k），其正负值持续增加，直至达到临界值，断层活化。
![]()
图 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 faultFigs.（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在流固耦合模拟中，我们分别在断层上下盘设置了监控线，便于对应力应变等关键参数进行过程分析。图7绘制了断层上盘的孔隙压力变化曲线，横坐标对应从上至下沿断层上盘的距离，曲线颜色对应于不同采样时刻。由于断层下盘的孔隙压力变化曲线在形态和数值上都与上盘相似，故只展示上盘结果。图中由深到浅的蓝色曲线对应于0—20 d的断层上盘孔隙压力值，由浅到深的红色曲线则对应于20—40 d的孔隙压力值。结果表明，断层上盘活化前后的孔隙压力分布具有明显的规律性：围岩和结晶基底界面存集了大量向下扩散渗流的流体，使其孔隙压力值始终处于高位；断层活化前，位于围岩上半部分的孔隙压力值呈U形分布，靠近注液点的上端和渗透率变化界面的下端具有更高的孔隙压力值，而位于结晶基底的下半部分的孔隙压力则快速降低；对应于前5天注液阶段的断层上端孔隙压力值（前5条蓝色曲线左端）逐渐上升，停止注液后的孔隙压力值则逐渐降低；断层活化后，其上端到中部的孔隙压力值呈线性增加特征，下半部分趋势与活化前的快速降低吻合；总体上围岩中的断层面孔隙压力高于结晶基底，与渗透率减小具有直接关系。
![]()
图 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 activationFigs. （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图8给出了断层上下盘的∆CFS值分布曲线，图中的浅蓝色条带表示接近断层上下端点的区域，浅绿色条带表示结晶基底界面附近区域。从图8a中可知断层活化前上下盘的∆CFS值随着流体的注入和渗透而逐渐增加，特别是结晶基底部分较为明显。断层上下盘的∆CFS值由红色和蓝色实线绘制，高值区主要集中在断层两端，且正负有分。断层中段的上盘∆CFS为正，下盘为负，值域区间主要在±0.5 MPa之间，断层活化后其上下盘的∆CFS值则几乎保持稳定，体现于不同时间采样点的曲线基本重复（图8b）。断层两端的∆CFS值达到5—6 MPa，这与断层模型两端的几何形态和约束有关，其值大小不能完全反映实际情况。活化前的孕震阶段结晶基底以下的断层∆CFS值更高，体现出更容易活化的特征，事实上后续的MW3.9地震正发生于这一区域。
![]()
图 7 沿断层上盘的孔隙压力变化曲线曲线颜色对应于不同采样时刻的断层上盘孔隙压力值，蓝色为断层活化前0—20天，红色为模拟的20—40 天，断层活化后Figure 7. The curve of pore pressure variation along the upper plate of the fault is shownThe 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 reactivation3.2 东断层存在时西断层的活化过程模拟
本节进一步分析东西两条断层同时存在时，西断层的活化过程以及东断层对西断层活化的影响。数值模拟采用完全相同的步骤计算三个阶段的孔隙压力、Y方向位移和∆CFS值，并进行6个关键时间点的切片展示（图9）。整个模拟过程中流体的注入、扩散、聚集模式与仅有西断层的情况类似，但是部分流体也渗透到了东断层之中（图9c−f）。这一流体渗透模式也使得东断层下部（结晶基底部分）出现了∆CFS的显著变化，但其集中程度和数值显然小于西断层（图9m−p）。实际观测表明东断层对流体的响应更为迅速，率先活化（Bao，Eaton，2016；Gao et al，2022），这与东断层的岩石物理性质及强度有关。本模拟虽然得出流固耦合作用下东断层的活化趋势，但其活化过程本文不作分析，重点关注的仍是西断层在流体作用下的延迟活化过程。西断层活化后东断层下端点的辐射花样也发生了细微变化，体现出断层之间的相互触发作用，尽管后期并未观测到东断层附近的新发地震。
![]()
图 8 断层活化前（a）、后（b）沿断层上下盘的∆CFS变化曲线上盘∆CFS对应浅红色曲线，下盘∆CFS为深灰色曲线，每条曲线代表一个模拟采样间隔单位（d）蓝色条带对应断层两侧端部，绿色条带对应代表断层与结晶基底交界区域Figure 8. The ∆CFS variation curves before （a） and after （b） fault activation along the fault planeThe 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本节的模拟结果表明西断层附近区域流固耦合场的演化模式与存在东断层的情况大体类似。但是东断层的存在是否会使得西断层更易活化呢?已有研究表明可能存在两条断层之间的流体通道，使流体通过东断层传导到西断层并使其成核活化（Gao et al，2022）。实际观测中东断层先被及时激活，而西断层则是延迟活化，其机制可能是东断层已有位错、流体扩散和流体通道的共同作用。基于建模的复杂度和可靠性原因，本研究暂不考虑潜在流体通道，只从流体注入和扩散角度分析两条断层共存时的相互影响，且重点关注西断层的活化过程，因为其对应该时段该区域最大震级事件。因此，我们将西断层活化前一天的库仑应力场进行对比分析，并标注出∆CFS集中点的数值（图10）。结果表明，在仅有西断层的情况下，经过5天注水和15天扩散的流体作用，西断层端点出现2.56 MPa的∆CFS高值区，MW3.9地震破裂点附近出现1.62 MPa的∆CFS高值区。而在两条断层并存的情况下，上述∆CFS高值区分别增加至2.98 MPa和2.11 MPa。这表明东断层在结晶基底中的∆CFS集中区一定程度上增加了西断层的库仑应力，使得西断层更容易活化。注意到东西断层之间在结晶基底内形成了一个∆CFS高值区联通带（图10b），这是邻近断层间应力相互扰动的结果，当某一断层应力达到其临界值后将会率先活化。
![]()
图 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-existedFigs.（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 points4. 讨论
本文模拟结果表明西断层在活化过程中存在三处变化显著的位置：第一处是靠近注液点的断层近端（沿断层走向0—50 m处），此处流体扩散到达的时间更短、孔隙压力更高；第二处是断层远端部分（沿断层走向550—600 m处），其机理是断层尖端的应变累积和应力集中；第三处是沿断层400—450 m处，此厚度非常小的矩形面，端点对应矩形短边，上下盘对应矩形的长边。断层贯穿两个不同岩石物理属性的地层，周围应力场和形变值集中分布。本研究假设断层是一个矩形，当断层上下盘（矩形长边）发生相对滑动时，短边会出现局部变形以达到符合震级的相对滑动距离，但是有限元方法的连续性使其仍然保持连接，而无法实现上下盘的实际分离。这种约束使得图8b中断层左右两端出现了库仑应力变化的异常高值区。第一处和第三处断层位置是断层的上下端点，其变化会受到模型本身设置的限制。此外，活化前位于结晶基底内的断层附近出现更为集中的库仑应力增加，这是因为模型中结晶基底具有更高的弹性模量，即该结晶基底的刚度和脆性更强，更容易引起相对更大震级的地震。这与本文研究区的实际观测相符合，东断层和西断层附近发生的较大震级地震均发生于结晶基底内部。
本文将注水速率设定为23.32 kg/s，这与加拿大福克斯溪地区Duvernay地层的实际平均施工强度相符。由于该段压裂时的返排率低至7%，说明大量流体存留在地下构造缝隙之中，因此并未在模拟中考虑流体的返排影响。经过测试，发现改变流体注入速率也会显著影响断层附近的库仑应力变化值，从而改变活化时间。当增加或减小注入量时，断层被激活的时间也会相应地减小或增加。这一分析具有一定的实际意义。如果能够事先计算出“安全可控诱发震级范围”的流体注入量或速率，就能为压裂施工提供有价值的建议。同理，这一规律可以应用到风险发生后，调整施工参数以调控地震风险当中。
将西断层附近发生的地震事件投影到其活化后库仑应力变化∆CFS的场图中，得到图11中的结果。图11a和11b分别展示了只有西断层以及东西断层并存的两种情况，图中灰点表示2015年1月23日MW3.9事件之前的地震，浅绿点表示该事件及其之后的地震。显然，MW3.9地震及其余震和后续诱发地震几乎完全位于西断层活化产生的库仑应力增加区，且震源深度与震级存在正相关关系（假设深度为正值）。该结果表明本文的多场耦合数值模拟结果具有合理性和可解释性。这里的库仑应力改变量∆CFS主要通过式（5）计算，除了表达断层位错产生的本身应力场外，还包含了孔隙压力的影响。模型中两个注液点附近的库仑应力呈现出集中趋势，并且随流体扩散到断层附近，这使得库仑应力场变得复杂。另外，东侧灰点所指示的诱发地震也位于东侧断层的库仑应力增加区（图11b）。
![]()
图 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最后需要指出本研究存在的局限性。文中模拟计算结果显示断层上下盘位移和库仑应力集中区能够与实测地震相吻合，虽然能够体现出模拟的合理性和主要的断层活化机制，但是仍有较多细节需要进一步探索：第一，本文仅从压裂注水和流体扩散的角度分析西断层活化过程，并未考虑前震与MW3.9主震的成核关系，后续需要结合已有研究（Gao et al，2022），从时空关系和机理上进行更多更细致的分析；第二，本文用于计算的模型为二维模型，且只考虑了断层的逆冲机制部分，这样不免忽略掉了一些三维空间计算应该考虑的走滑分量；第三，两条断层的实际活化顺序是先东断层而后西断层（图2），而东断层具有远场触发的特点，西断层则具有延迟触发的特点，其特殊性还需要在进一步搜集相关数据的基础上继续研究；第四，断层的临界滑动特征还需要进一步分析，其孕震、成核、破裂、滑动、愈合过程仍要结合速率状态方程等手段进行深入讨论。上述问题值得在今后的进一步工作中系统研究。
5. 结论
本文针对加拿大阿尔伯塔省Duvernay地层附近两条明显的、高度倾斜的地震活跃带（断层）展开数值模拟研究。根据丛集地震展布和震源机制信息，在模型中假设两条断层的产状和基本性质。在PKN裂缝扩展模型和多孔弹性理论的基础上，依据当地地质地层构造数据建立模型，并进行多物理场耦合运算。结果表明，在仅有西断层的情况下，经过5天注水和15天扩散的流体作用，西断层附近出现了2.56 MPa的∆CFS高值区，而在两条断层并存的情况下，∆CFS高值区增加至2.98 MPa。这说明东断层在结晶基底中的∆CFS集中区一定程度上增加了西断层的库仑应力，使得西断层更容易活化。因此，我们认为东断层的存在会改变西断层被激活的时间和强度，这与性质相似断层之间的相互触发有关。同时，模拟计算出西断层活化产生的库仑应力增加区覆盖了MW3.9地震及其余震，表明实际地震空间分布与模型中的断层设定及应力演化结果相符。本研究能够一定程度上解答中强型诱发地震的“尾随效应”难题，为注液诱发地震的危害性预测提供新的物理视角。最后，对本研究存在的不足进行了总结。
-
![]()
图 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
计量
- 文章访问数: 282
- HTML全文浏览量: 146
- PDF下载量: 239
