流體-地質力學耦合建模表征水力壓裂誘發地震：以加拿大Fox Creek地區為例

惠鋼， 陳勝男*， 顧斐

1 加拿大卡爾加里大學化學與石油工程系， 卡爾加里 T2N1N4 2 中國石油勘探開發研究院， 北京 100083

0 引言

近年來，隨著水平井多級壓裂技術在頁巖氣開發中的廣泛應用，水力壓裂引發的人為地震活動顯著增加(張東曉和楊婷云，2015；Atkinson et al.，2016； Lei et al.，2019；Schultz et al.，2017).全球一些較大震級的誘發地震事件，例如在中國四川盆地發生的2017-01-28筠連4.9級地震，2018-12-18興文5.7級地震，2019-09-07威遠5.4級地震和加拿大西部盆地發生的2019-03-04 Red Deer 4.2地震事件，都被證實與附近的水平井壓裂作業有時空相關性(Lei et al.，2017，2019,2020；Schultz and Wang，2020).

水平井多級壓裂被認為是一項有效開發頁巖氣資源的前沿技術.在水力壓裂過程中，壓裂液被加壓并泵入目標地層，造成巖石的拉伸破壞，產生特定方向的壓裂縫(一般平行于最大水平主應力).然后將配伍的支撐劑注入到新生壓裂縫中，使其在后續油氣開采過程中保持張開狀態，從而使頁巖氣被高效快速采出.一般而言，該壓裂過程會產生不連續的、瞬時的微地震事件，而且其震級通常小于0，多數在-1級以下(Chen等，2018).然而，當水力壓裂縫在地下擴展并溝通潛在的斷層或者大裂縫時，就會誘發較大震級的地震活動.一些中等強度(0.5

加拿大西部盆地具有豐富的頁巖氣資源，水平井多級壓裂技術在該盆地取得了巨大成功(Atkinson et al.，2016).近年來，該盆地發生了很多中等強度的地震事件.這些事件都被證實與附近的水力壓裂作業密切相關(Bao and Eaton，2016；Eaton et al.，2018；Schultz and Wang，2020).圖1顯示了加拿大西部盆地記錄的M>2.5地震事件.其中，震級2.5被認為是加拿大西部盆地誘發地震事件的長期監測閾值(Schultz et al.，2017).圖1還顯示出誘發地震事件群具有明顯的聚集特征，即沿落基山脈變形邊緣帶 (圖1虛線)分布.其中的震源球指示地震事件的震源機制(Wang et al.，2018).圖中顯示Fox Creek地區是該盆地誘發地震事件的頻發區域之一.前人研究表明，該地區誘發地震與地層孔隙壓力(Eaton and Schultz，2018)、縱向接近基底的距離(Pawley et al.，2018)、平面接近生物礁邊界(Schultz et al.，2017)及斷層應力狀態(Zhang et al.，2019)等地質因素密切相關.此外，Fox Creek地區單口水平井壓裂液注入總體積達到32944 m3.根據壓裂液注入體積與誘發地震震級的經驗線(McGarr，2014)，該注入體積所對應的地震事件最大震級可以達到M4.0.綜上可知，Fox Creek地區水力壓裂誘發的地震事件受到地質和工程因素的共同制約(Pawley et al.，2018；Schultz et al.，2018).

圖1 加拿大西部盆地震級大于2.5的地震事件分布(截止至2019-04-30)圓圈代表誘發地震事件.圓球顯示工業活動誘發地震的震源機制(Wang等, 2018).圓球上方數字代表事件發生事件.虛線代表落基山脈變形邊緣帶.Fig.1 Map of historical seismicity of ML≥2.5 up to 2019-04-30 in Western CanadaThe cycles denote the induced earthquakes. The magnitude-scaled “beach balls” show the focal mechanisms of induced earthquakes (Wang et al., 2018). Four digits above the beach ball represent the nucleated time for associated earthquakes. The dashed line indicates the deformation edge of the Rocky Mountains.

近年來，國內外學者采用不同的研究方法對誘發地震進行了表征.國外方面，Bao和Eaton(2016)研究了誘發地震活動的觸發模式，提出壓裂過程中的應力變化能夠激活斷層.Lele等 (2017)建立了壓裂誘發地震的概念模型，認為水力壓裂縫與附近斷層溝通、并引起斷層內部壓力上升是導致斷層活化并誘發地震事件的主要原因.國內方面，張致偉等(2012)研究了四川地震活動與注水井壓力與注水量的關系, 認為誘發地震頻次、強度與注水量呈正相關.何登發等(2019)利用高精度三維地震資料探索誘發地震的形成機制，認為四川長寧地區的地震活動與現今青藏高原東南緣的大地構造背景密切相關.本文在Lei等(2017)研究的基礎上，綜合分析地質、巖石力學及流體力學等特征，來表征加拿大西部盆地水力壓裂誘發的地震事件.

本文采用一種綜合地質、巖石力學及流體力學的研究方法，對加拿大西部盆地Fox Creek地區2015年2月8日發生的M3.0誘發地震事件進行精細表征.首先，建立區域地質、完井、測井、三維地震及壓裂施工等資料庫.其次，利用該區高分辨率三維地震勘探資料, 采用螞蟻體追蹤技術識別潛在斷層.然后，根據測井曲線和施工數據等定量求取巖石力學參數及地應力，建立三維地質力學模型，計算出水力壓裂縫的幾何形狀和縫內流體壓力的分布特征.最后，建立流體-地質力學耦合模型，計算水力壓裂過程中斷層附近的孔隙壓力及局部應力變化，利用摩爾-庫侖破裂準則判定斷層激活的時間與空間位置，從而揭示本次誘發地震事件的觸發機制并提出相應的風險控制對策.

1 誘發地震觸發機理及斷層失穩準則

一般認為，誘發地震的觸發機制主要包括兩方面(見圖2).一種是注入流體通過高滲透單元(如壓裂縫)直接溝通斷層，造成斷層內流體壓力增加；另一種是在流體注入或者開采過程中，彈性應力從目標地層傳遞到斷層，引起斷層面應力變化 (Ellsworth，2013).上述兩種情況均導致斷層失穩并發生移動，從而誘發地震事件(Healy et al.，1968； Raleigh et al.,1976).

圖2 流體注入或采出誘發地震示意圖(Ellsworth,2013)Fig.2 Schematic diagram of fluid injection or extraction-induced seismicity(Ellworth,2013)

具體而言，當斷層內部孔隙壓力或者彈性應力足夠高、超過斷層的破壞強度時，斷層將發生滑動進而誘發地震活動.為表征斷層的失穩特征，國內外研究者們提出了一系列判別標準，其中最常用的是摩爾-庫侖破壞準則.該準則是確定斷層失穩和預測強震活動后余震分布的有效工具(Catalliet al.，2013).摩爾-庫倫破壞準則定義如下：

CFS=τ+μ(σn+pp),

(1)

其中，CFS為斷層面上的庫侖破壞應力，MPa；μ為摩擦系數，無因次；τ為剪切應力(滑動方向為正)，MPa；σn為正應力(拉伸方向為正)，MPa；Pp為孔隙壓力，MPa.在三維應力場中，斷層面(l,m,n)上的剪切應力τ及正應力σn計可以通過最大，垂向及最小主應力(σ1，σ2和σ3，MPa)計算，其計算公式為(Lei et al.，2020)

σn=σ1l2+σ2m2+σ3n2,

(2)

τ2=(σ1-σ2)2l2m2+(σ2-σ3)2m2n2

+(σ3-σ1)2n2l2.

(3)

為方便運算，本文對上述公式進行了二維近似.在二維摩爾圓中，斷層面左右側的剪切應力τ及有效正應力σn的計算公式為

(4)

(5)

(6)

其中，τl和τr分別為斷層左、右側的剪切應力，MPa；σ1與σ3分別為斷層面承受的最大、最小主應力，MPa；β為斷層走向與最大水平主應力的夾角.

在斷層滑動過程中，斷層周圍的巖石會發生彈性變形，導致斷層附近局部應力場發生變化，產生庫侖破壞應力(CFS)正變化和負變化區域.研究表明，在CFS正變化區域，斷層將處于不穩定狀態，斷層發生失穩、滑移的幾率顯著增加(Stiros and Kontogianni，2009；King and Devès，2015).其判別方程基于(1)式得出，Δ為各參數變化量.

ΔCFS=Δτ+μ(Δσn+Δpp).

(7)

圖3 (a—b)壓裂過程中孔隙壓力變化和應力變化引起的庫侖應力失穩Fig.3 (a—b)The Coulomb Failure tiggered by changes of pore pressure and in-situ stress

摩爾圓可以直觀顯示上述判別準則.通常根據斷層面的剪切應力及有效應力來繪制摩爾圓及庫侖破壞線.具體來說，在壓裂過程中，壓裂縫擴展并溝通地層的潛在斷層時，注入流體會增加斷層內的孔隙壓力，從而降低斷層面上的正應力，使斷層向摩爾-庫侖失穩向移動，直至到達庫侖破壞線引起斷層滑動(見圖3a).在這種情況下，注入流體通常會以超過1 km·d-1的速率在已溝通的斷層內擴散，因此會很快觸發地震事件(Tadokoro et al.，2000).

另一種情況是壓裂縫沒有直接溝通斷層.但壓裂液的注入引起地層孔隙壓力增大，并導致區域巖石彈性變形，從而改變斷層附近的局部應力場.在這種情形下，摩爾圓也將接近庫侖破壞線并激活斷層(見圖3b).其中，垂直主應力變化是由負載壓力變化和Biot系數決定，而水平主應力變化則與巖石的泊松比有關，并且通常小于垂直應力變化(Haddad et al.，2020).

綜上所述，我們可以計算出壓裂過程中斷層附近的孔隙壓力及局部應力變化，利用摩爾-庫侖破裂準則來判定斷層激活的時間及空間位置，從而揭示誘發地震的觸發機制.

2 區域地質背景及M3.0地震事件

Fox Creek地區位于加拿大西部盆地中部(見圖1).該地區地層發育較全，自上至下依次為泥盆系、寒武系及前寒武系花崗巖基底，上述地層的總厚度超過2000 m(圖4a).其中Duvernay地層的頁巖氣資源為該地區主要的油氣開發對象.該地層廣泛發育富有機質頁巖，上覆泥巖為主的Ireton地層，下覆灰巖為主的Swan Hills地層，地層之間均為整合接觸.Duvernay地層厚度約為40 m，埋藏中深低于海平面約2520m.地層平均有效孔隙度為6.6%，平均滲透率為392 nD(3.92×10-20m2)，平均總有機質(TOC)含量為4.5%(Creaney and Allan，1990；Weir et al.，2019).

系地層組地層巖性泥盆系Wabamum灰巖WinterburnGraminia灰巖Blueridge灰巖Calmar泥巖和頁巖Nisku灰巖WoodbendIreton泥巖和頁巖Duvernay泥巖和頁巖Beaverill LakeSwan Hills灰巖Slave Point蒸發巖Elk PointGillwood砂巖Prairie-蒸發巖Muskeg灰巖寒武系Upper灰巖Lower砂巖前寒武系基底變質花崗巖

Fox Creek地區作為誘發地震事件的高發區域，其地震活動最早可追溯到2013年12月(Schultz et al.，2017).此后，隨著水平井多級壓裂技術的廣泛應用，該區誘發地震活動頻繁發生.其中有超過10次震級大于3.0的地震事件被證實與水力壓裂作業密切相關(Bao and Eaton，2016；Eyre et al.，2019).這些地震事件主要歸因于該地區特定的地質因素和工程因素，如Duvernay地層高孔隙壓力、與Precambrian基底的距離、與Swan Hills生物礁邊緣的距離、被激活斷層的初始應力狀態等地質因素及壓裂液的累計注入體積等工程因素(Eaton et al.，2018；Pawley et al.，2018).截止至2018年12月，本地區有606口水平井進行了水力壓裂作業.統計表明，平均每口水平井壓裂級數為28.7段，水平段長度為2135 m，壓裂液注入速率為9.5 m3·min-1，累計注入體積達到32944 m3.如前文所述，該注入體積可誘發震級超過M4.0的地震事件(Mcjarr，2014).

本文以2015年2月8日發生的M3.0地震事件群為例，對誘發地震的觸發機制進行綜合研究.該地震群發生在Fox Creek地區中部的一口壓裂水平井附近.圖5a表明，研究區與本區發育的大型斷層都有一定的距離(Pawley et al.，2018).其中與最近斷層的距離已超過5 km，表明本次事件與大型斷層無直接關系.因此，需要利用高精度的三維地震勘探資料，識別發育規模相對較小的斷層，從而揭示構造與誘發地震的關系.

2015年2月5日—14日，該水平井從南向北實施了水力壓裂作業.根據壓裂參數的統計結果，該水平井共壓裂級數20段，水平段長度1.5 km.平均施工壓力55.8 MPa，平均壓裂液注入速率為8.7 m3·min-1，平均每個壓裂段注入1904 m3，累計注入量為38078 m3.區域地震臺站監測結果顯示，壓裂過程中共有109起地震事件發生(Bao and Eaton，2016).圖5b顯示了該誘發地震群與水平井壓裂段的空間分布特征.可以看出，大部分地震事件發生在水平井的東側，以及Duvernay地層下方200～1300 m范圍內.

現場資料表明，壓裂開始40 h后，南部地震群在水平井筒東側約200 m處被激活.其中包括本次研究的M3.0事件.該事件的震源分析表明(圖5b震源球)，被激活的斷層呈近N-S走向展布(Schultz et al.，2017).南部地震群的發生持續了40 h，直到2月8日結束(圖5c中第6級壓裂段結束).隨后，中部地震群被激活，結束后，北部地震群在距井筒東側300 m處被激活，其震級范圍為M1.5～M2.9(圖5b—c).

圖5 (a) 研究區位置及大斷層分布,底圖為Duvernay頂部構造等值圖； (b) 地震事件及水平井空間分布特征； (c) 壓裂施工與誘發地震事件時間順序Fig.5 (a) Location of the study area and distribution of large faults.The base map dnotes the elevation of top Duvernay Formation； (b) Spatial distribution of the induced earthquakes and horizontal well； (c) Time sequence of fracturing operation and induced events

利用地震震級與地震事件頻率關系擬合曲線的斜率值(b值)，可識別地震事件群的觸發類型.一般而言，對于天然斷層激活所誘發的地震事件群，其b值一般小于1(Burridge and Knopoff，1967).而水力壓裂活動所誘發的地震事件群，通常b值要大于1.5(Eaton et al.，2014).本次地震事件群的b值為0.84，為典型的天然斷層活化特征.根據水平井壓裂縫的位置推斷，本次地震群可能是由于壓裂縫擴展并溝通了斷層的基底部分，從而激活斷層并誘發一系列地震事件.為進一步明確其觸發機制，需要結合測井、三維地震、巖石力學及應力場、施工數據等資料對本次誘發地震進行綜合表征.

3 螞蟻體追蹤技術識別斷層

該地區進行過三維地震勘探，因此可用來探尋構造特征與誘發地震事件的關系.在三維反射地震解釋中，螞蟻體跟蹤技術通常用于識別小型斷層，其工作流程為：(1)井震結合建立關鍵井合成地震記錄，提取工區相關層位；(2)對三維地震數據進行預處理，包括中值濾波、結構平滑、高斯空間濾波和帶通濾波處理；(3)從地震屬性中提取Chaos、方差和傾角偏差，進行地震反射邊緣檢測；(4)將濾波方法應用于地震解釋后處理，開展螞蟻體自動跟蹤，刻畫出斷層的空間分布(Pedersen et al.，2002; Hui et al., 2021).

圖6 (a—c) 三個地層的螞蟻體屬性切片及地震事件 (黑點) 疊合圖; (d) A—A′地震剖面斷層解釋結果及斷層兩側的地震事件.A—A′剖面位置見圖5bFig.6 (a—c) Horizontal cross-section view of ant tracking property and seismic event (black spot) in three formations; (d) Interpretation of faults bounded by related seismic events in the A—A′ seismic profile.A—A′ section position is shown in Fig.5b

圖6a—6c分別為Duvernay、Cambrian、Precambrian地層的螞蟻體屬性切片及地震事件疊合圖.可以看出，大部分地震事件集中在基底，并與基底復雜的斷層分布特征大致相同.基于螞蟻體屬性特征，分析斷層與地震事件之間的空間關系，共識別出3條近垂直的天然斷層(圖6d).三條斷層均由Precambrian基底向上延伸至Duvernay地層頂部.其中，斷層1的走向和傾角分別為175°和89°，與前人的震源反演結果一致(Schultz et al.，2017).斷層1高度約為1300 m，延伸距離為810 m.其發育規模表明，該斷層被激活后具有誘發最大震級4.8級地震的潛力(Zoback and Gorelick，2012).斷層2和斷層3走向大致相同，為北東—南西向，延伸長度分別約為820 m和590 m，傾角分別為89°和88°.根據前人的研究結果(Yehya et al.，2018；Fan et al.，2019)，本次研究將三條斷層帶的厚度設置為45 m，包括5 m的斷層核及兩側各20 m的斷層破裂帶.圖6d還顯示出壓裂誘發的地震事件在縱向上沿三條斷層兩側分布.結合圖5b—5c,這些地震事件與水平井的某些壓裂段具有時間與空間相關性.這表明，壓裂過程中水平井壓裂縫擴展并溝通了三條基底斷層，引起斷層失穩并誘發了M3.0地震事件.

4 三維地質力學模型與壓裂縫擴展

4.1 地層異常高壓與局部地應力場

研究表明，Fox Creek地區的高地層壓力及局部應力場對誘發地震有重要影響(Eaton et al.，2018；Pawley et al.，2018).本文綜合測井和壓裂施工資料，對本區的地層孔隙壓力及地應力參數進行了定量求取.

首先，利用壓裂過程中的瞬時關井壓力來估算孔隙壓力.本例水平井20級壓裂段的平均瞬時關井壓力及測試深度分別為56.8 MPa和3387 m.因此孔隙壓力Pp梯度估算為16.77 kPa·m-1，接近前人研究結果的16.8 kPa·m-1(Eaton et al.，2018).相對于其他地區的地層壓力，Fox Creek地區的地層異常高壓成為誘發地震事件在該區頻繁發生的重要地質因素(Pawley et al.，2018).

根據井眼垮塌和鉆井誘導縫等資料，結合Shen等(2019)的研究結果，確定本地區自Ireton至基底地層的最大水平主應力SHmax方向均為NE45°.同時，為便于耦合模擬，本文對最大及最小應力軸進行了水平近似處理.垂直主應力Sv梯度則由附近直井的密度測井數據求得，為24.6±0.25 kPa·m-1.其中0.25為分析誤差，設定為計算值的1%(Zoback，2007).根據施工過程中的閉合壓力，估算本區最小水平主應力Shmin梯度為20.92±0.56 kPa·m-1.最后，利用經驗公式SHmax= 3Shmin- 2Pp(Zoback,2007)，估算SHmax梯度為29.16±1.68 kPa·m-1.主應力計算結果與前人研究結果基本一致(Eyre et al.，2019；Zhang et al.，2019；Shen et al.，2019).

通過比較三應力張量的大小(SHmax≥Sv≥Shmin)，可以確定本區發育的斷層類型為走滑斷層.根據公式(4)—(6)，繪制出三條被激活斷層的摩爾圓，以表征其初始應力狀態.如圖7所示，如果忽略壓裂過程中的局部應力變化，則激活斷層1、斷層2和斷層3分別需要增加3.5、2.0、0.5 MPa的斷層內部孔隙壓力.

圖7 三條斷層的初始應力狀態斜線為庫侖破壞線.圓點代表各斷層相關地震事件的震源機制.摩爾圓內的等值線表示激活斷層所需要增加的孔隙壓力(MPa).Fig.7 Original stress state of three faultsThe pink line represents the Coulomb failure line.The dots denote focal mechanisms of related seismic events.The contours within the Mole circle show the increased pore pressure (MPa) required for the fault activation.

4.2 巖石力學及物性參數的確定

利用巖石力學參數可以定量表征壓裂過程中巖石彈性變形對于斷層激活的影響.研究區附近的一口直井測量了從Duvernay地層到Muskeg頂部地層的縱波(VP，m·s-1)及橫波(VS，m·s-1)測井數據.從Muskeg底部到Precambrian地層的速度數據則采用前人的研究成果(Ronald et al.，2019).可據此計算各地層的動態泊松比(υ,無因次)和楊氏模量(E，GPa)，再結合三軸應力實驗得到的靜態參數對其進行校正.此外，根據附近取心井的巖心化驗數據，對各地層的孔隙度和滲透率進行了統計(Weir et al.，2019).表1統計了從Duvernay地層到Precambrian基底各地層的彈性參數及物性參數.其中靜態泊松比范圍為0.19～0.30, 靜態楊氏模量范圍為48～70 GPa, 孔隙度和滲透率分別為0.02～0.07和5.2×10-20～6.46×10-18m2.根據以上結果，建立了整合應力場、巖石力學參數和物性參數的三維地質力學模型.

表1 各地層彈性參數及物性參數統計Table 1 Statistics of elastic and physical parameters of each formations

4.3 人工壓裂縫擴展模擬

當對水平井實施水力壓裂作業時，壓裂縫擴展方向通常與水平最大主應力SHmax平行.對于水力壓裂縫擴展的表征，通常是利用施工數據和地質參數，進行壓裂過程中的凈壓力(施工壓力減去最小主應力及摩擦阻力等)的歷史擬合，從而模擬壓裂縫的空間擴展.本文采用PKN模型模擬壓裂縫的三維實時擴展.在PKN模型中，斷裂面處于平面應變狀態，斷裂界面為橢圓形，裂縫高度為固定值.假定忽略濾失效應，則在某一時刻壓裂縫半長和寬度的計算公式為 (Yew and Weng,1997)

(8)

(9)

(10)

其中，L(t)為t時刻壓裂縫半長，m；W(t)為t時刻壓裂縫寬度，m；P(t)為t時刻壓裂縫內初始壓力, MPa；G為剪切模量, GPa；qo為注入速率，m3·min-1；μ為注入流體黏度，cp；h為壓裂縫縫高，m.

模擬結果表明，壓裂縫半長為104～166 m，平均為131 m(圖5b).裂縫寬度為0.0135～0.0173 m，平均為0. 0166 m.根據裂縫滲透率與裂縫寬度的經驗公式，計算裂縫滲透率約為7.5×10-12m2，遠超過Duvernay地層基質滲透率3.94×10-19m2.上述壓裂縫屬性特征將用于流體-地質力學耦合模型.

5 耦合流體-地質力學有限元模擬

5.1 流體-地質力學耦合方程

本文采用線性孔隙彈性理論來表征壓裂過程中的孔隙壓力擴散和應力擾動.基于多孔介質連續性方程及固體彈性變形方程，可推導出流體-地質力學耦合方程(Wang，2000)：

(11)

(12)

其中，f(x,t)為單位體積應力，N·m-3；q(x,t)為流體注入速率，m3·s-1；u為平移向量，m；ν為泊松比，無因次；ε為體積應變，無因次；α為Biot系數，無因次；M為Biot模量，無因次；k為滲透率，μm-2；η為動態流體黏度，cp.

上述方程可以利用有限元方法進行求解.本文使用COMSOL Multiphysics軟件進行流體-地質力學耦合建模,定量表征壓裂過程中的孔隙壓力及應力變化，再利用摩爾-庫侖失穩準則來確定斷層被激活的位置及時間.

5.2 建立有限元耦合模型及初始化設置

三維模型在x、y、z方向的尺寸分別設為2.5 km×3.5 km×1.8 km.z底部設為-4.1 km.根據實際地層特征，自Ireton至Precambrian地層厚度依次設置為160、40、100、100、200、1200 m.水平井深度設為-2520 m.壓裂縫沿NE45°擴展，其幾何參數及屬性采用前文計算結果.三條解釋斷層嵌入到模型中，其斷層帶厚度設為45 m，包括中間5 m的致密斷層核部分及兩側各20 m的斷層破裂帶(Fan et al.，2019).其中斷層核及破裂帶滲透率分別設為經驗值1×10-20m2、1×10-14m2(Yehya et al.，2018).同時，模型中各地層物性及巖石力學參數設定為表1對應數值.壓裂液密度設定為1200 kg·m-3，壓縮系數設定為4.6×10-10/Pa，動態黏度設定為0.4 mPa·s.圖8a—8b為耦合模型的三維視圖及網格劃分.

圖8 (a) 三維耦合模型初始化示意圖; (b) 模型網格劃分.在壓裂縫及斷層附近加密網格Fig.8 (a) 3D view of poroelastic model initialization; (b) 3D mesh using triangular elements. The mesh surrounding the fractures and faults are refined

在耦合模型中還設置了初始條件和邊界條件.假設模型的初始流體系統處于流體靜力平衡狀態，則可以設置初始孔隙壓力p(t=0)=0和初始應力張量σ(t=0)=0.因此,耦合模擬結果即為孔隙壓力和應力張量的相應變化值.將各壓裂段的壓裂液注入體積轉化為各段凈注入壓力，并作為其耦合模擬過程中的驅動力，其計算公式如(10)所示.壓裂縫延伸長度與寬度如公式(8)和(9)所示.凈注入壓力的作用時間同各壓裂段實際施工時間保持一致.前一段壓裂結束后，壓裂縫末端延伸凈壓力接近于0(接近于地層孔隙壓力)，后一段壓裂開啟，如此循環至最后一級壓裂結束.模型邊界設定為固定邊界，即頂部處于無牽引狀態，側向和底面邊界的位移保持為零.為便于耦合模擬，本文對最大及最小應力軸進行了水平近似處理.即垂直上覆壓力為變量，最大及最小主應力取決于上覆壓力，因此三軸應力簡化為單軸應力狀態.設定耦合模型受控于線性彈性力學和達西定律，并在壓裂縫及斷層附近加密網格(圖8b).通過固體力學與多孔介質流體流動的耦合模擬，計算出壓裂過程中的應力擾動和孔隙壓力變化.

5.3 模擬結果及M3.0地震事件的誘發機理

圖9a—9c顯示壓裂80h后壓裂層位(Duvernay)及震源位置所在基底(Basement)的孔隙壓力平面變化情況.可以看出，斷層1兩側破裂帶區域的孔隙壓力變化明顯.表明致密斷層核對流體流動的阻礙作用不明顯.圖9d為MW3.0事件震源位置處(圖9b十字)ΔPp及ΔCFS隨時間變化的情況.MW3.0地震時間發生時，ΔPp及ΔCFS瞬時值分別達到3.6 MPa及2.2 MPa.基于有效正應力變化值(Δσn-ΔPp=-3.12 MPa)和剪切應力變化值(Δτ=0.14 MPa)，圖7所示的斷層1將達到失穩狀態，可以斷定斷層1被激活.同樣，壓裂過程中孔隙壓力及局部應力的變化使得斷層2和3也被激活.值得注意的是，孔隙壓力增加仍然是影響MW3.0誘發地震事件的主要因素(圖9d).

另外，為檢測耦合模型是否受到設定邊界的影響.將模型尺寸擴大為5 km×7 km×1.8 km(原模型體積的4倍).模擬結果如圖9d所示.可以看出，新模型震源位置的ΔCFS及ΔPp與原模型差別不大(變化值7%左右).因此在原模型尺度大小及邊界條件的設定下，邊界大小對于本次模型的結果輸出影響不大.

圖9 (a—b) 壓裂80 h后壓裂層位及基底層位ΔPp分布圖; (c) 地震事件疊合于(b)圖上; (d) 原邊界及大邊界模型震源處(圖9b十字)的ΔPp及ΔCFS隨時間變化情況.壓裂段及地震事件也在圖中繪出Fig.9 (a—b) ΔPp at the Duvernay Formation and basement, respectively, at t=80 hours after the onset of fracturing operations; (c) Induced events overlapped Fig.9b; (d) Temporal ΔPp and ΔCFS at the nucleation position of MW3.0 earthquake. Fracturing stages and induced events are also depicted

6 討論

目前，針對Fox Creek地區誘發地震的風險控制對策主要體現在施工控制方面，如優化新鉆水平井與已知斷層的空間位置、控制壓裂液注入速率及累計注入體積等(Schultz et al.，2018).其中，增大水平井筒與已知斷層的距離被證實可以有效地降低地震風險.

如圖10a所示，根據螞蟻體反演屬性，本例中的斷層1向南延伸3 km，呈SE-NW走向.2015年11月，即北部水平井(本例)完井9個月之后，在斷層西側1.5 km左右完鉆一口NS向、水平井段長度為3 km的新水平井.自11月11日至29日，對該水平井實施壓裂施工作業，共壓裂級數31段.壓裂液注入速率及累計注入體積分別為9.5 m3·min-1和38303 m3，與北部水平井壓裂施工參數8.7 m3·min-1和38078 m3大致相同.

圖10 (a) Duvernay地層頂部的螞蟻體屬性分布.下部圓球顯示南部水平井壓裂過程中的最大震級M1.25事件的震源機制.黑色虛線為D—D′剖面位置; (b) D—D′剖面的地震事件空間分布特征.圓球大小和顏色分別代表地震事件的震級大小和發生時間Fig.10 (a) The horizontal cross-section of ant tracking attributes showed the simulated pore pressure changes at the Basement of Fault 1 at t=80 hours after the onset of HF operations; (b) Vertical F-F′ section view of earthquakes distributions. The balls denoted induced earthquakes, colored by time and scaled by magnitude

地震監測結果表明，壓裂過程中在該井以東600～1500 m處誘發了一系列震級較小的地震事件，其中最大為M1.25事件(圖10a).圖10b比較了南、北水平井壓裂過程中的兩次主要地震事件的空間分布特征.可以看出，南部水平井誘發的事件震級明顯變小，而且深度主要集中于壓裂地層(Duvernay)上下50 m范圍內，這與北部地震事件主要集中在基底的情況形成鮮明的對比.此外，南部最大地震事件M1.25的震源機制表明(圖10a下部震源球)，被激活斷層的走向為NE 28.2°±8.3°(Zhang et al.，2019).這與北部被激活斷層(斷層1)的近南北走向顯著不同.因此可以推斷，由于南部水平井與斷層1位置較遠，其壓裂施工并未激活斷層1.相反，在南部水平井和斷層1之間或存在一個小規模斷層(圖10a斷層4).南部地震群的產生可能由于該斷層與南部井的壓裂縫連接而被激活.此外，模擬結果表明，南部井壓裂過程中的CFS變化僅為0.15 MPa，遠小于北部井的1.53 MPa.由此可見，優化水平井與已知斷層的距離可以有效降低誘發地震活動的風險.該現場措施為優選水平井井位、降低地震風險提供了可靠的現場依據.

7 結論

本文以Fox Creek地區2015年2月8日發生的M3.0地震事件為例，基于螞蟻體識別的潛在斷層及三維地質力學模型，通過流體-地質力學的耦合模擬，明確了水力壓裂施工過程中孔隙壓力及局部應力的變化特征，揭示出本次地震事件的觸發機制，并提出相應的風險控制對策.具體包括：(1)利用螞蟻體追蹤方法識別出三條由Precambrian基底向上延伸至Duvernay頂部的近垂直斷層；(2)基于三維地質力學模型，模擬出水平井壓裂縫半長約為131 m；(3)部分壓裂縫與解釋斷層相溝通，導致部分注入流體沿斷層高滲透破裂帶向下擴散，在基底地層激活斷層并誘發MW3.0地震事件.其中孔隙壓力的增加是斷層失穩的主要因素；(4)增大水平井與已知斷層之間的距離，被現場證實可以有效地降低地震風險.因此在致密儲層進行水平井鉆井及壓裂施工作業之前，明確地下斷層的分布狀態至關重要.

致謝圖1地震活動目錄來自加拿大Alberta地震活動記錄(www.inducedseismicity.ca/catalogues/).MW3.0誘發地震事件及壓裂施工數據來自Bao和Eaton(2016).本次研究使用的三維地震數據由加拿大TGS公司提供.該研究得到Canada First Research Excellence Fund (CFREF)的資助.感謝卡爾加里大學地質系David Eaton教授的幫助.感謝兩位審稿人的寶貴意見.