水平井降壓法和熱激法水合物開采對海底邊坡穩(wěn)定性的影響1)

譚 琳 劉 芳

(同濟大學土木工程防災(zāi)國家重點實驗室,上海 200092)

(同濟大學巖土及地下工程教育部重點實驗室,上海 200092)

引言

天然氣水合物(簡稱水合物)是由籠形水分子晶格內(nèi)嵌天然氣分子(主要是甲烷)構(gòu)成的類冰狀晶體化合物,廣泛分布在陸地永久凍土區(qū)和深海海床.據(jù)估計,全球水合物含碳總量是目前傳統(tǒng)化石能源(石油、煤和天然氣)總含碳量的2 倍[1-4],被認為是最有應(yīng)用前景的未來能源之一.在開采海底水合物的過程中,甲烷氣體將大量釋放,可能弓起海床孔壓急劇增長,水合物儲層力學性能顯著劣化,從而誘發(fā)海床失穩(wěn)或海底滑坡[5-7].

目前,已有學者研究了海床溫度變化和海平面升降對水合物富集區(qū)海底邊坡穩(wěn)定性的影響.比如,Gidley 和Grozic[8]通過室內(nèi)模型試驗研究了水合物因升溫分解產(chǎn)生的氣體在坡體中的運移路徑及其誘發(fā)邊坡失穩(wěn)的過程.張旭輝等[9-10]通過離心機試驗研究水合物分解過程中水合物層與上覆層間裂隙的演變,認為水合物與上覆層間的裂隙是導致邊坡失穩(wěn)的重要因素.張建紅等[11]采用離心模型試驗觀測了水合物分解致使孔壓增大所誘發(fā)的海底坡體變形、裂縫發(fā)展與滑裂面形態(tài)的動態(tài)演變.上述研究均表明,因水合物分解而弓起的孔壓上升是導致海底邊坡失穩(wěn)的重要因素,準確預(yù)測孔壓變化是合理評價水合物富集區(qū)海床穩(wěn)定性的關(guān)鍵.為此,不同學者提出了水合物分解過程的孔壓計算模型.如Nixon和Grozic[12]基于不排水假設(shè)提出了水合物分解的孔壓計算模型,并將之弓入極限平衡法分析水合物分解所弓發(fā)的海底滑坡;Sultan 等[13]考慮溫度、壓強、孔隙水化學成分以及孔隙尺寸的影響,基于水合物的熱動力學-化學平衡模型計算水合物分解過程的孔壓,并利用極限平衡法分析了挪威Storegga 大滑坡;Kwon 等[14-15]建立了描述水合物分解過程的熱-流-化學耦合模型求解孔隙水壓力,并基于極限平衡法研究了海床升溫對海底水合物儲層邊坡穩(wěn)定性的影響.

水合物開采誘發(fā)海床失穩(wěn)的研究還相對較少.目前常見的水合物開采方法包括降壓法、熱激法、化學抑制法、二氧化碳置換法、固態(tài)流化法以及不同方法的組合.其中,熱激法最為直接,但能量損失較大,氣體收集困難;降壓法最簡單有效,但開采后期效率降低;目前降壓和熱激結(jié)合的方法被認為是最有前景的長期開采方法.在水合物開采的海床失穩(wěn)致災(zāi)研究方面,蔣明鏡等[16]采用流體動力學和離散元耦合法研究了水合物熱開采誘發(fā)的海底滑坡全過程,分析了開采位置和水合物儲層空間分布的影響,但忽略了水合物開采的實際過程,假設(shè)水合物瞬間分解.Kimoto 等[17-18]建立熱-力-流-化學全耦合數(shù)值模型,采用有限元法研究了降壓加熱結(jié)合法開采過程的海床變形特性.劉鋒[19]采用有限差分法研究了水合物分解過程中坡體的位移演變與滑動過程.Zander等[5]分析了豎井降壓開采對黑海Danube 深海邊坡穩(wěn)定性的影響,認為降壓開采會造成局部海床下沉,但不會降低邊坡的整體穩(wěn)定性.Moridis 等[6-7]研究了水平井降壓和熱開采過程中的海床穩(wěn)定性,認為熱激法在滲透性較差的水合物儲層中會產(chǎn)生足夠高的孔壓致使海床失穩(wěn),降壓法在開采過程中對邊坡穩(wěn)定性無不利影響,但停采后海床穩(wěn)定性由于孔壓增加而有所降低.綜上,水合物開采有可能誘發(fā)海底邊坡失穩(wěn),不同水合物開采方法對海床穩(wěn)定性的影響研究還有待系統(tǒng)開展.本文基于邊坡穩(wěn)定極限平衡分析方法的框架,弓入考慮水合物開采過程的熱-流-化學耦合數(shù)值分析模型,模擬不同開采方法中水合物分解鋒面擴展和瞬態(tài)孔壓演化過程,得到不同開采方法的邊坡穩(wěn)定安全系數(shù)和失穩(wěn)模式的變化規(guī)律.

1 分析方法

有別于陸地滑坡,水合物開采誘發(fā)海底滑坡伴隨著水合物分解相變所產(chǎn)生的土體強度劣化和孔壓變化使力學特性更加復雜[20-21].鑒于邊坡穩(wěn)定極限平衡分析方法的物理概念清晰,分析框架簡單實用,本文嘗試在該方法框架內(nèi)考慮上述影響,弓入水合物開采過程的瞬態(tài)孔壓場,計算有效應(yīng)力變化,并考慮水合物分解過程水合物儲層抗剪強度的弱化,采用極限平衡分析方法確定水合物開采過程的邊坡穩(wěn)定安全系數(shù),從而建立判定水合物開采誘發(fā)海底滑坡與否的實用方法.圖1 為本文分析流程,具體分析步驟包括兩部分.

圖1 計算流程示意圖Fig.1 Flow chart

首先,采用TOUGH+HYDRATE 程序[22]對水合物開采過程進行熱-流-化學(THC)耦合分析.TOUGH+HYDRATE 是美國勞倫斯-伯克利國家實驗室針對水合物合成/分解過程開發(fā)的非等溫多相流多場耦合數(shù)值模擬器.該模擬器將土骨架簡化為彈性多孔介質(zhì),將水合物合成/分解環(huán)境界定為介質(zhì)孔隙,其中涉及氣、水、水合物和抑制劑4 種物質(zhì),這4種物質(zhì)存在于氣態(tài)、液態(tài)、冰和水合物態(tài)4 種相態(tài).基于Darcy 定律、Fourier 定律和水合物熱動力學穩(wěn)定性質(zhì)[23]建立質(zhì)量守恒與能量守恒方程來描述水合物合成/分解中的滲流、熱傳導和相變過程.該過程中考慮了熱-流-化學多場耦合,即溫度場、孔壓場和化學場(系統(tǒng)中各相態(tài)的物質(zhì)組分)之間的相互作用.孔壓場和溫度場通過水合物的熱動力學穩(wěn)定特性控制水合物的分解/合成,并通過滲流過程中的氣液運移影響化學場;反過來,水合物合成/分解涉及到固、液、氣間的相變會改變孔隙壓力,并通過滲流過程來影響孔壓場,同時水合物合成/分解過程中的放熱/吸熱效應(yīng)會弓起溫度的變化,并通過熱傳導來影響溫度場;孔壓場的變化驅(qū)動滲流過程中氣液運移通過熱對流來影響溫度場,溫度場通過影響各物質(zhì)的流動特性(比如黏滯系數(shù)、密度等)來影響孔壓場.TOUGH+HYDRATE 通過求解質(zhì)量守恒與能量守恒方程來獲取孔壓場、溫度場和化學場等的演變,可以模擬不同布井方案的降壓法、熱激法、抑制劑法以及組合法開采過程,是目前模擬水合物開采最流行的模擬器之一[6,24],并在水合物開采模擬方面得到了大量應(yīng)用[25-26].為了驗證TOUGH+HYDRATE 對模擬水合物分解過程的有效性,本文對文獻[27]的一維水合物分解室內(nèi)試驗進行了模擬,計算模型及參數(shù)與文獻[27]保持一致.該實驗在裝有砂樣的不銹鋼反應(yīng)釜中合成水合物,然后通過降壓法使水合物分解.反應(yīng)釜長50 cm,直徑3.8 cm,砂樣粒徑300～450 μm,滲透率3.0×10-13m2,水與水合物飽和度分別為21.83%和29.61%,初始溫度為1.54°C,壓強為3.535 MPa.實驗中將出口壓強降至0.93 MPa 使水合物分解,同時在出口處收集氣體.一維數(shù)值模擬中將出口單元設(shè)置為恒溫恒壓邊界條件,另一端邊界單元設(shè)置為恒溫、不透水、不透氣條件.如圖2 所示,本文數(shù)值模擬結(jié)果與文獻[27]的試驗及數(shù)值模擬結(jié)果吻合較好.

圖2 水合物一維分解實驗與數(shù)值結(jié)果的對比Fig.2 Comparison between the one-dimensional hydrate dissociation experimental and numerical results

其次,根據(jù)指定時刻的孔壓場計算邊坡的有效應(yīng)力場,并根據(jù)水合物分解鋒面的狀態(tài),確定土體強度參數(shù)(即有效黏聚強度和有效內(nèi)摩擦角)的空間分布,將其輸入商業(yè)軟件SLOPE/W,基于極限平衡分析框架分析邊坡穩(wěn)定性.通過SLOPE/W 提供的最危險滑裂面位置搜尋算法,確定最危險滑裂面位置及相應(yīng)的邊坡穩(wěn)定安全系數(shù),最終得到水合物開采中邊坡穩(wěn)定安全系數(shù)隨時間的變化規(guī)律和最危險滑裂面的演變.

當水合物分解時,海床土體將由初始飽和狀態(tài)過渡到非飽和狀態(tài),不同狀態(tài)的有效應(yīng)力均采用簡化的Bishop 模型[28]計算,即

其中,σ 和σ′分別為土中總應(yīng)力和有效應(yīng)力,u為平均孔壓(如無特殊說明,后文中孔壓均指平均孔壓),采用下式計算

其中,SG和SA分別為孔隙氣相和孔隙液相的飽和度,PG和PA分別為孔隙氣壓和孔隙液壓.

2 數(shù)值建模

2.1 問題描述

海底邊坡以緩坡為主,但峽谷海山等地形區(qū)域存在陡坡,如南海北陸坡白云凹陷區(qū)和珠江口盆地的峽谷區(qū)少數(shù)邊坡達到40°左右[29-30].陡坡穩(wěn)定性差,對外界擾動更加敏感,容易發(fā)生滑坡失穩(wěn).本文參照文獻[16],模擬我國南海北陸坡水合物富集區(qū)的典型45°陡坡在水合物開采過程中的穩(wěn)定性.該海底邊坡的計算網(wǎng)格如圖3 所示,計算域?qū)?000 m,其中邊坡段寬600 m,坡頂水深1000 m,坡底水深1600 m.海床自上而下分別為上覆層、水合物儲層和下臥層,其中水合物儲層埋深100 m,厚150 m,初始水合物飽和度為50%,假設(shè)水合物中的天然氣成分均為甲烷.

假設(shè)采用水平井開采,井軸線垂直于邊坡橫截面,考慮兩種常用開采方法:

(1)單井降壓法.假設(shè)開采井的壓強保持6 MPa持續(xù)開采10 年,隨后停采.該工況重點分析了水平井位置的影響,分別考慮坡肩(A點)、坡中(B點)和坡趾(C點)3 種情況.

圖3 海底邊坡的計算模型Fig.3 Illustration of the submarine slope model

(2)雙井熱激法.注熱井和生產(chǎn)井分別位于B點和D點,兩井的水平距離為40 m.注熱井以0.05 kg/s 的速率注入熱鹽水(鹽度0.035),注水溫度考慮了50°C,70°C,90°C 三種情況;生產(chǎn)井設(shè)置為內(nèi)部固定邊界,開采過程中溫度保持為10.3°C,壓強保持為15.6MPa.

所有算例均假設(shè)海床表面溫度為4°C,地溫梯度為36.0°C/km,海水鹽度為0.035.坡體的上下邊界恒溫恒壓,兩側(cè)邊界假設(shè)為絕緣邊界(即在邊界處不發(fā)生物質(zhì)與熱量的傳遞).圖4 為邊坡的孔壓和溫度的初始空間分布情況.

圖4 邊坡初始條件Fig.4 The initial conditions in the slope

2.2 計算模型及參數(shù)

水合物開采過程THC 耦合分析的關(guān)鍵計算模型匯總為表1.開采井用透水單元來近似模擬,孔隙率取1.0,滲透率取5.0×10-9m2.表2 為地層的物理力學參數(shù),上覆層主要為黏性土,水合物儲層及下臥層以砂性沉積物為主,參數(shù)取值參考了中國南海北陸坡典型地層條件和已有相關(guān)數(shù)值模擬[31-33],水合物儲層的強度指標參考文獻[34]的試驗結(jié)果,水合物未分解時(對應(yīng)水合物飽和度0.5)有效內(nèi)摩擦角和有效黏聚強度分別取20°和0.5 MPa,水合物分解過程中內(nèi)摩擦角保持不變,黏聚強度隨著水合物飽和度降低而線性減小,當水合物完全分解時黏聚強度降至0.05 MPa.

表1 計算模型Table 1 Models used in this study

表2 地層的物理力學參數(shù)Table 2 Physics and mechanical parameters of strata

3 結(jié)果分析

3.1 降壓開采的計算結(jié)果

3.1.1 孔壓場變化和分解鋒面擴展

圖5 為開采井位于坡體中部時單位長度水平井的產(chǎn)氣曲線.開采初期水合物分解鋒面自井口向外擴散,產(chǎn)氣率迅速增加,在開采后第4 年達到峰值,隨后水合物分解鋒面向井周擴展,分解鋒面處壓強有所增加,產(chǎn)氣率逐漸下降,開采第7 年后產(chǎn)氣率大致趨于穩(wěn)定.在第10 年停采時,單位長度水平井的累計產(chǎn)氣量約為4×105m3(標準狀態(tài)下).

圖5 坡體中部降壓開采的產(chǎn)氣速率QP和產(chǎn)氣量VPFig.5 Profiles of production rate(QP)and cumulative volume(VP)of produced gas per unit length of a depressurized well in the mid-height of the slope

圖6 為開采過程中孔壓場的變化和水合物分解鋒面的擴展情況,其中等值線表示孔壓場,云圖表示水合物飽和度的空間分布.開采初期,由于井口降壓導致井周附近形成局部低壓,促使水合物迅速分解;隨著水合物開采的進行,低壓區(qū)域由井口向四周擴散,弓起水合物在井周持續(xù)分解,分解鋒面由井口逐漸向周圍擴散.

圖6 孔壓變化(單位MPa)和水合物分解鋒面擴展Fig.6 The variation of pore pressure in MPa and propagation of dissociation front during production

3.1.2 海底邊坡穩(wěn)定性

當降壓開采井布置在坡體中部時,圖7 為邊坡安全系數(shù)在開采過程以及停采后的變化情況.在水合物開采之前,邊坡安全系數(shù)約為1.14;隨著水合物的開采,安全系數(shù)逐漸增大,開采后第7 年安全系數(shù)大致穩(wěn)定在1.78 左右;當?shù)?0 年停采時,安全系數(shù)迅速跌落,在停采3 年后穩(wěn)定在1.08 左右,比開采前初始安全系數(shù)降低5%.

圖7 坡體中部降壓開采時邊坡安全系數(shù)變化Fig.7 Evolution of safety factor with a production well installed at the mid-height of the slope

圖8 給出了典型時刻(圖7 的時刻A至E)最危險滑裂面的位置.開采前(時刻A),最危險滑弧經(jīng)過坡腳(見圖8(a));開采1 年后(時刻B),受降壓的影響,井口周圍處的土體孔壓降低,有效應(yīng)力增大,雖然井口處水合物分解導致分解區(qū)黏聚強度下降,但影響范圍尚小,土體抗剪強度整體增強,最危險滑弧向深部發(fā)展(見圖8(b)),相應(yīng)地,邊坡的穩(wěn)定安全系數(shù)有所提高;開采5 年后(時刻C),隨著水合物進一步分解,低壓區(qū)由井口向四周和深部繼續(xù)擴散,最危險滑弧由深部轉(zhuǎn)為淺部(見圖8(c)),位于水合物儲層分解區(qū)上方;開采10 年后,最危險滑弧位于深部水合物儲層之下(見圖8(d)).上述分析結(jié)果表明,在采用降壓法開采水合物的過程中,雖然水合物分解會導致分解區(qū)土體黏聚強度有所喪失,但是影響范圍有限,因孔壓降低、有效應(yīng)力提高所弓起的摩擦強度提高占主導,水合物分解區(qū)的土體整體抗剪強度反而有所增加,最危險滑弧不經(jīng)過水合物分解區(qū)域,邊坡安全系數(shù)較開采前有所提高.

圖8 坡體中部降壓開采時最危險滑裂面變化Fig.8 The evolution of potential slip surface with a depressurized well installed at the mid-height of the slope

停采后,水合物分解區(qū)的孔壓隨之升高并逐漸恢復到靜水壓狀態(tài),由于水合物分解所弓起的黏聚強度降低對邊坡穩(wěn)定性的不利影響開始變得顯著,水合物分解區(qū)成為抗剪強度薄弱的部位,最危險滑弧通過水合物分解區(qū)域(圖8(e)),邊坡安全系數(shù)低于水合物開采前的初始值.

3.1.3 降壓井位置的影響

圖9 為降壓井位置對產(chǎn)氣效率的影響.當井口壓強相同時,開采井位于坡趾時產(chǎn)氣率最高.坡趾的初始孔壓較高,啟動坡趾開采時的壓強變化梯度較大,降壓區(qū)擴散較快,水合物分解速率較大.隨著開采井位置上移至坡肩,產(chǎn)氣效率逐漸降低.

圖9 降壓開采位置對產(chǎn)氣率QP的影響Fig.9 Effect of the production well location on the production rate,QP

降壓井位置對安全系數(shù)的影響見圖10.不同布井工況下,水合物分解均弓起分解區(qū)土體的黏聚強度下降,但降壓所導致的摩擦強度增加更為顯著,因此開采過程中邊坡穩(wěn)定性均有所增強;當開采井位于坡體中部時,安全系數(shù)提高最為顯著.停采后,由于開采位置不同,導致分解區(qū)(抗剪強度弱化的區(qū)域)空間位置不同,邊坡的安全系數(shù)也有所區(qū)別,當開采井位于坡體中部時,停采后的安全系數(shù)最低.圖11對比了不同開采位置情況下最危險滑弧的位置.開采10 年后,當開采位置在坡體中部或坡肩時,最危險滑弧位于水合物分解區(qū)下方的深部土體,出露位置在坡趾附近;當開采位置在坡趾時,最危險滑弧位于水合物分解區(qū)上方的淺部土體,出露位置在坡趾上方.停采后10 年,自坡體中部開采的工況下,最危險滑弧通過水合物分解區(qū)(圖8(e)),安全系數(shù)低于開采前的初始值;坡趾和坡肩開采工況下,安全系數(shù)較坡體中部開采工況大,其中,當坡趾開采時,由于停采后坡體孔壓未完全恢復到初始狀態(tài)(圖12(b)),安全系數(shù)仍略高于開采前的初始值.

圖10 降壓開采位置對邊坡安全系數(shù)的影響Fig.10 Effect of the well location on the safety factor of the slope

圖11 開采10 年后的最危險滑弧位置Fig.11 The potential slip surface after 10-year production with a depressurized well installed

圖12 停采10 年后的最危險滑弧位置Fig.12 The potential slip surface after 10-year termination of production with a depressurized well installed

3.2 開采方法的影響

圖13 對比了在坡體中部實施單井降壓開采和雙井熱激開采的產(chǎn)氣率以及不同注水溫度熱激法開采的產(chǎn)氣率.熱激法的產(chǎn)氣速率在開采前期高于降壓法,但迅速降低,開采4 年后,產(chǎn)氣率低于降壓法.由于熱激法會導致孔壓上升,抑制水合物分解,因此產(chǎn)氣率曲線在早期出現(xiàn)波動.對比不同注水溫度對應(yīng)的產(chǎn)氣率,可以看出注水溫度越高產(chǎn)氣率越高.

圖13 不同開采方法的井口產(chǎn)氣率QPFig.13 Production rate under different production methods

圖14 為不同開采方法的邊坡安全系數(shù)變化情況.在開采過程中,熱激法的邊坡穩(wěn)定性明顯差于降壓法,并且注水溫度升高會使最危險時刻提前.注熱過程中,由于溫度升高,流體(尤其是氣體)膨脹產(chǎn)生超孔壓,注水溫度越高,超孔壓越高(圖15).注水溫度為90°C 的情況下,開采1 個月后井口周圍最高超孔壓達到3.8 MPa(圖15(c)),有效應(yīng)力顯著減小,同時水合物分解導致分解區(qū)黏聚強度下降,從而造成土體的抗剪強度明顯降低,邊坡安全系數(shù)比開采前初始值降低7%.圖16 給出了注水溫度為90°C 情況下開采后1 年和10 年兩個典型時刻的坡體超孔壓分布情況.隨著開采的進行,因水合物分解使土體滲透性提高,且流體由開采井抽出,孔壓逐漸消散,邊坡安全系數(shù)稍有回升,但仍低于初始值(圖14).如圖17 所示,與降壓開采工況不同,熱激開采工況下邊坡的最危險滑弧始終通過水合物分解區(qū).

圖14 不同開采方法的邊坡穩(wěn)定安全系數(shù)Fig.14 Evolution of safety factor with different production methods and well in the middle of the slope

圖15 熱激法不同注水溫度開采30 天后超孔壓場與最危險滑裂面的位置Fig.15 The field of excess pore pressure and the position of potential slip surface after 30-day thermal stimulation with injected hot water of different temperatures

圖16 熱激法的超孔壓場Fig.16 The evolving field of excess pore pressure under thermal stimulation

圖16 熱激法的超孔壓場(續(xù))Fig.16 The evolving field of excess pore pressure under thermal stimulation(continued)

圖17 熱激法最危險滑裂面的演變Fig.17 The evolution of the potential slip surface under thermal stimulation

4 結(jié)論

本文基于極限平衡法的分析框架,考慮水合物開采過程的邊坡瞬態(tài)孔壓及抗剪強度的變化,分析了水合物開采井位置與開采方法對邊坡穩(wěn)定安全系數(shù)的影響.主要結(jié)論如下:

當采用降壓法開采時,在開采過程中,邊坡穩(wěn)定性主要受孔壓降低所致的摩擦強度升高的影響,邊坡安全系數(shù)有較大提高;停采后,孔壓恢復,水合物分解所致黏聚強度下降的影響凸顯,邊坡安全系數(shù)顯著下降,低于開采前的初始值.

降壓開采井的位置影響產(chǎn)氣率和邊坡穩(wěn)定性.若井口壓強相同,隨著降壓開采井由坡趾上移至坡肩,產(chǎn)氣效率逐漸降低;若在坡體中部布設(shè)降壓井,開采過程的邊坡穩(wěn)定性最好,但停采后的邊坡穩(wěn)定性最差.

當采用雙井熱激法開采時,在升溫過程可能產(chǎn)生較大超孔壓,導致坡體內(nèi)有效應(yīng)力降低,邊坡安全系數(shù)顯著降低,若開采策略不當,存在誘發(fā)性滑坡的風險.

需要指出的是,本研究假設(shè)邊坡土體的滲透率和儲層水合物飽和度的空間分布均勻,當實際地層中土體滲透率和水合物分布不均勻時,即使采用降壓法開采也可能造成局部孔壓上升而顯著影響邊坡的整體穩(wěn)定性,土層性質(zhì)的空間變異性影響還有待進一步研究.