流體-地質(zhì)力學(xué)耦合建模表征水力壓裂誘發(fā)地震：以加拿大Fox Creek地區(qū)為例

惠鋼， 陳勝男*， 顧斐

1 加拿大卡爾加里大學(xué)化學(xué)與石油工程系， 卡爾加里 T2N1N4 2 中國(guó)石油勘探開(kāi)發(fā)研究院， 北京 100083

0 引言

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

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

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

圖1 加拿大西部盆地震級(jí)大于2.5的地震事件分布(截止至2019-04-30)圓圈代表誘發(fā)地震事件.圓球顯示工業(yè)活動(dòng)誘發(fā)地震的震源機(jī)制(Wang等, 2018).圓球上方數(shù)字代表事件發(fā)生事件.虛線代表落基山脈變形邊緣帶.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.

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

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

1 誘發(fā)地震觸發(fā)機(jī)理及斷層失穩(wěn)準(zhǔn)則

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

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

具體而言，當(dāng)斷層內(nèi)部孔隙壓力或者彈性應(yīng)力足夠高、超過(guò)斷層的破壞強(qiáng)度時(shí)，斷層將發(fā)生滑動(dòng)進(jìn)而誘發(fā)地震活動(dòng).為表征斷層的失穩(wěn)特征，國(guó)內(nèi)外研究者們提出了一系列判別標(biāo)準(zhǔn)，其中最常用的是摩爾-庫(kù)侖破壞準(zhǔn)則.該準(zhǔn)則是確定斷層失穩(wěn)和預(yù)測(cè)強(qiáng)震活動(dòng)后余震分布的有效工具(Catalliet al.，2013).摩爾-庫(kù)倫破壞準(zhǔn)則定義如下：

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

(1)

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

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

(2)

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

+(σ3-σ1)2n2l2.

(3)

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

(4)

(5)

(6)

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

在斷層滑動(dòng)過(guò)程中，斷層周圍的巖石會(huì)發(fā)生彈性變形，導(dǎo)致斷層附近局部應(yīng)力場(chǎng)發(fā)生變化，產(chǎn)生庫(kù)侖破壞應(yīng)力(CFS)正變化和負(fù)變化區(qū)域.研究表明，在CFS正變化區(qū)域，斷層將處于不穩(wěn)定狀態(tài)，斷層發(fā)生失穩(wěn)、滑移的幾率顯著增加(Stiros and Kontogianni，2009；King and Devès，2015).其判別方程基于(1)式得出，Δ為各參數(shù)變化量.

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

(7)

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

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

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

綜上所述，我們可以計(jì)算出壓裂過(guò)程中斷層附近的孔隙壓力及局部應(yīng)力變化，利用摩爾-庫(kù)侖破裂準(zhǔn)則來(lái)判定斷層激活的時(shí)間及空間位置，從而揭示誘發(fā)地震的觸發(fā)機(jī)制.

2 區(qū)域地質(zhì)背景及M3.0地震事件

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

系地層組地層巖性泥盆系Wabamum灰?guī)rWinterburnGraminia灰?guī)rBlueridge灰?guī)rCalmar泥巖和頁(yè)巖Nisku灰?guī)rWoodbendIreton泥巖和頁(yè)巖Duvernay泥巖和頁(yè)巖Beaverill LakeSwan Hills灰?guī)rSlave Point蒸發(fā)巖Elk PointGillwood砂巖Prairie-蒸發(fā)巖Muskeg灰?guī)r寒武系Upper灰?guī)rLower砂巖前寒武系基底變質(zhì)花崗巖

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

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

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

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

圖5 (a) 研究區(qū)位置及大斷層分布,底圖為Duvernay頂部構(gòu)造等值圖； (b) 地震事件及水平井空間分布特征； (c) 壓裂施工與誘發(fā)地震事件時(shí)間順序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

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

3 螞蟻體追蹤技術(shù)識(shí)別斷層

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

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

4 三維地質(zhì)力學(xué)模型與壓裂縫擴(kuò)展

4.1 地層異常高壓與局部地應(yīng)力場(chǎng)

研究表明，F(xiàn)ox Creek地區(qū)的高地層壓力及局部應(yīng)力場(chǎng)對(duì)誘發(fā)地震有重要影響(Eaton et al.，2018；Pawley et al.，2018).本文綜合測(cè)井和壓裂施工資料，對(duì)本區(qū)的地層孔隙壓力及地應(yīng)力參數(shù)進(jìn)行了定量求取.

首先，利用壓裂過(guò)程中的瞬時(shí)關(guān)井壓力來(lái)估算孔隙壓力.本例水平井20級(jí)壓裂段的平均瞬時(shí)關(guān)井壓力及測(cè)試深度分別為56.8 MPa和3387 m.因此孔隙壓力Pp梯度估算為16.77 kPa·m-1，接近前人研究結(jié)果的16.8 kPa·m-1(Eaton et al.，2018).相對(duì)于其他地區(qū)的地層壓力，F(xiàn)ox Creek地區(qū)的地層異常高壓成為誘發(fā)地震事件在該區(qū)頻繁發(fā)生的重要地質(zhì)因素(Pawley et al.，2018).

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

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

圖7 三條斷層的初始應(yīng)力狀態(tài)斜線為庫(kù)侖破壞線.圓點(diǎn)代表各斷層相關(guān)地震事件的震源機(jī)制.摩爾圓內(nèi)的等值線表示激活斷層所需要增加的孔隙壓力(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 巖石力學(xué)及物性參數(shù)的確定

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

表1 各地層彈性參數(shù)及物性參數(shù)統(tǒng)計(jì)Table 1 Statistics of elastic and physical parameters of each formations

4.3 人工壓裂縫擴(kuò)展模擬

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

(8)

(9)

(10)

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

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

5 耦合流體-地質(zhì)力學(xué)有限元模擬

5.1 流體-地質(zhì)力學(xué)耦合方程

本文采用線性孔隙彈性理論來(lái)表征壓裂過(guò)程中的孔隙壓力擴(kuò)散和應(yīng)力擾動(dòng).基于多孔介質(zhì)連續(xù)性方程及固體彈性變形方程，可推導(dǎo)出流體-地質(zhì)力學(xué)耦合方程(Wang，2000)：

(11)

(12)

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

上述方程可以利用有限元方法進(jìn)行求解.本文使用COMSOL Multiphysics軟件進(jìn)行流體-地質(zhì)力學(xué)耦合建模,定量表征壓裂過(guò)程中的孔隙壓力及應(yīng)力變化，再利用摩爾-庫(kù)侖失穩(wěn)準(zhǔn)則來(lái)確定斷層被激活的位置及時(shí)間.

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

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

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

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

5.3 模擬結(jié)果及M3.0地震事件的誘發(fā)機(jī)理

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

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

圖9 (a—b) 壓裂80 h后壓裂層位及基底層位ΔPp分布圖; (c) 地震事件疊合于(b)圖上; (d) 原邊界及大邊界模型震源處(圖9b十字)的ΔPp及ΔCFS隨時(shí)間變化情況.壓裂段及地震事件也在圖中繪出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 討論

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

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

圖10 (a) Duvernay地層頂部的螞蟻體屬性分布.下部圓球顯示南部水平井壓裂過(guò)程中的最大震級(jí)M1.25事件的震源機(jī)制.黑色虛線為D—D′剖面位置; (b) D—D′剖面的地震事件空間分布特征.圓球大小和顏色分別代表地震事件的震級(jí)大小和發(fā)生時(shí)間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

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

7 結(jié)論

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

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