海底水氣突涌運(yùn)移規(guī)律與麻坑半徑定量評(píng)價(jià)

程章，趙俞成，洪義，王立忠

(浙江大學(xué)建筑工程學(xué)院，浙江杭州，310058)

在海底大陸架斜坡環(huán)境中，天然氣水合物會(huì)分解形成游離狀態(tài)的水氣向上突涌遷移，在海床表面形成與天然氣(自由態(tài)或者溶解態(tài))和天然氣水合物有關(guān)的麻坑[1]。麻坑不僅會(huì)危害到海洋基礎(chǔ)設(shè)施的安全，而且給海洋基礎(chǔ)設(shè)施的選址帶來(lái)困難。通過(guò)海底原位觀測(cè)、聲學(xué)探測(cè)和鉆孔取芯等方法，RIBOULOT等[2]將海底深部流體運(yùn)移形成的麻坑分為I 型麻坑和II 型麻坑這2 種。海底淺埋高壓水氣突涌通常形成II型麻坑[2-3]。

前人對(duì)I型麻坑的形成機(jī)制及模型預(yù)測(cè)的研究相對(duì)完善，而對(duì)于II 型麻坑的形態(tài)演變的理解不足[4-6]。RIBOULOT等[2]根據(jù)沉積物原位樣本和地震反射剖面數(shù)據(jù)，提出V 形區(qū)域假說(shuō)，認(rèn)為在超孔壓作用下聚集區(qū)前端會(huì)形成局部初始裂縫，并在壓力作用下擴(kuò)展，最終突涌至海床表面形成V 形麻坑。SULTAN 等[7]針對(duì)尼日爾三角洲的麻坑及水合物聚集區(qū)的地質(zhì)數(shù)據(jù)，提出V 形區(qū)域內(nèi)存在有大量天然氣水合物和游離天然氣，并認(rèn)為水合物的快速生長(zhǎng)和水合物緩慢溶解是導(dǎo)致II型麻坑發(fā)展的主要機(jī)制。楊志鵬等[8]認(rèn)為海底陸坡上發(fā)育的一部分麻坑是地層小規(guī)模斷層所致，并將其歸結(jié)為斷裂成因[9]。為了模擬II型麻坑的形成過(guò)程，有關(guān)學(xué)者利用等效模擬法，使用有限元軟件分別模擬了麻坑的半徑和水氣遷移路徑[2,10]。

盡管上述研究為理解II型麻坑的形態(tài)發(fā)展提供了一條重要的途徑，但是它們忽略了斜坡海床中原生主應(yīng)力偏轉(zhuǎn)對(duì)于水氣遷移路徑的影響，無(wú)法準(zhǔn)確預(yù)測(cè)現(xiàn)場(chǎng)麻坑的位置及半徑，且采用傳統(tǒng)有限元法等效模擬裂縫擴(kuò)展會(huì)高估突涌路徑的寬度。裂縫在土中沿裂尖的最大主應(yīng)力方向擴(kuò)展，而海床坡角和靜止土壓力系數(shù)等因素會(huì)對(duì)原生主應(yīng)力偏轉(zhuǎn)產(chǎn)生顯著影響，進(jìn)而影響水氣突涌路徑以及麻坑半徑。因此，在預(yù)測(cè)淺埋高壓水氣突涌形成II型麻坑的擴(kuò)展路徑時(shí)，必須要考慮應(yīng)力主軸旋轉(zhuǎn)這個(gè)重要因素。

本文作者基于RIBOULOT等[2]提出的V形區(qū)域假說(shuō)，首先，利用擴(kuò)展有限元和Cohesive 黏聚力單元相結(jié)合的數(shù)值模擬方法，建立二維平面應(yīng)變應(yīng)力-滲流耦合的水氣突涌模型，模擬II型麻坑的形成過(guò)程；然后，根據(jù)數(shù)值模擬結(jié)果，探究海床坡角和靜止土壓力系數(shù)對(duì)裂縫周圍應(yīng)力場(chǎng)、孔壓場(chǎng)以及裂縫擴(kuò)展路徑的影響；最后，綜合考慮海床坡角和靜止土壓力系數(shù)的影響，給出麻坑半徑的評(píng)價(jià)方法，為海洋基礎(chǔ)選址提供一定參考。

1 數(shù)值原理

1.1 應(yīng)力-滲流耦合基本理論

飽和軟黏土中的高壓劈裂過(guò)程會(huì)影響土體的應(yīng)力場(chǎng)，因此，數(shù)值模擬方法采用基于有效應(yīng)力原理的應(yīng)力-滲流耦合分析方法[11-12]。完全飽和的土層應(yīng)力平衡方程可用虛功原理表達(dá)：一定時(shí)間內(nèi)土體上的虛功等于作用在土體上的體力和面力所產(chǎn)生的虛功，土體的平衡方程為

式中：符號(hào)“:”為張量運(yùn)算中的雙點(diǎn)積運(yùn)算；V為土體體積；σ'為土體的有效應(yīng)力張量；pw為孔隙流體壓力；I為二階單位張量；δε˙虛應(yīng)變率張量；S為土體表面積；t為單位面積上的外力矢量；δv為虛速度矢量；f為單位體積上的體力矢量。

對(duì)于不可壓縮流體滲透通過(guò)飽和土體單元，根據(jù)連續(xù)性條件和質(zhì)量守恒定律可知，單位時(shí)間內(nèi)某一土體單元內(nèi)流體質(zhì)量的變化量與該單元流進(jìn)與流出的流體質(zhì)量之差相等，連續(xù)性方程可以寫(xiě)成

式中：t為時(shí)間；ρw為孔隙流體的密度；nw為土體的孔隙率；vw為孔隙流體的滲透速度矢量；n為面S上外法線方向的單位矢量。

1.2 擴(kuò)展有限元方法

擴(kuò)展有限元法(XFEM)由MO?S 等[13]提出，擴(kuò)展有限元法將裂縫視為完全分離的實(shí)體，該實(shí)體獨(dú)立于網(wǎng)格，且無(wú)需在裂縫擴(kuò)展過(guò)程中進(jìn)行裂尖范圍的網(wǎng)格重劃分，大幅減少了計(jì)算量。擴(kuò)展有限元方法通過(guò)引入富集函數(shù)來(lái)表示不連續(xù)位移場(chǎng)[14]：位移不連續(xù)函數(shù)表征裂縫面所引起的位移不連續(xù)；裂尖漸進(jìn)函數(shù)描述裂尖附近應(yīng)力奇異場(chǎng)。

1.3 黏聚力單元區(qū)本構(gòu)行為：裂紋萌生與擴(kuò)展

模擬飽和軟黏土的劈裂破壞時(shí)，Cohesive單元往往會(huì)受到3個(gè)方向的應(yīng)力，采用B-K斷裂準(zhǔn)則來(lái)計(jì)算復(fù)合型裂縫的總臨界應(yīng)變能釋放率Gc，表達(dá)式如下：

式中：Gc為復(fù)合型裂縫總臨界斷裂能釋放率；GIc，GIIc和GIIIc分別為法向、剪切方向和扭轉(zhuǎn)方向的臨界斷裂能釋放率；η為與材料本身特性相關(guān)的參數(shù)。

B-K斷裂準(zhǔn)則和黏聚力牽引分離損傷模型描述了復(fù)合型材料的斷裂行為[15-16]。圖1所示為黏聚力單元損傷模型。圖中δ0為單元初始損傷時(shí)對(duì)應(yīng)的張開(kāi)位移，δf為裂縫起裂時(shí)對(duì)應(yīng)的張開(kāi)位移。由圖1可見(jiàn)：黏聚力牽引分離損傷模型主要由2個(gè)參數(shù)組成，即最大拉應(yīng)力Tmax和總斷裂能釋放率Gc。針對(duì)各向同性的巖土材料，一般可假設(shè)3個(gè)方向的斷裂能釋放率相等[17]，由式(3)可得總臨界斷裂能釋放率Gc等于I 型臨界斷裂能釋放率GIc。基于線彈性斷裂力學(xué)可得到I型臨界斷裂能釋放率GIc與斷裂韌度KIc之間的關(guān)系為

圖1 黏聚力單元損傷模型[16]Fig.1 Cohesive element damage model[16]

式中：GIc為I 型臨界斷裂能釋放率；v為泊松比；KIc為I型斷裂韌度；E為彈性模量。

JOHNSON 等[18]基于加拿大新斯科舍省3 個(gè)場(chǎng)地的現(xiàn)場(chǎng)探頭測(cè)量結(jié)果，得到飽和細(xì)粒沉積物的拉伸斷裂韌度KIc與強(qiáng)度su近似呈線性關(guān)系：

式中：su為土體不排水抗剪強(qiáng)度。

本研究中使用式(5)計(jì)算飽和軟黏土的I型斷裂韌度，并代入式(4)得到相應(yīng)的斷裂能釋放率。

根據(jù)線彈性斷裂力學(xué)(LEFM)，裂縫起裂受最大主應(yīng)力控制，當(dāng)最大主應(yīng)力超過(guò)預(yù)定義的最大閾時(shí)，就會(huì)形成裂縫。新形成的裂縫方向垂直于裂紋尖端附近的局部最大拉應(yīng)力方向。

1.4 裂縫面流體流動(dòng)模型

假設(shè)Cohesive 單元內(nèi)的流體是連續(xù)不可壓縮的牛頓流體，流體流動(dòng)可以當(dāng)成順裂縫面方向的切向流和垂直于裂縫面的法向流[19-20]。

WANG 等[21]針對(duì)不同壓縮系數(shù)的裂隙流體，開(kāi)展了在巖土體中的平面應(yīng)變裂縫擴(kuò)展模擬試驗(yàn)，發(fā)現(xiàn)流體的壓縮性對(duì)巖土體中裂縫的擴(kuò)展影響很小。故本研究中不考慮高壓滲流破壞模擬中裂縫內(nèi)流體的壓縮系數(shù)，統(tǒng)一按不可壓縮的水進(jìn)行模擬。

2 海底淺埋高壓水氣突涌形成麻坑的擴(kuò)展有限元模型

2.1 XFEM網(wǎng)格與邊界條件

圖2所示為平面應(yīng)變數(shù)值模型及網(wǎng)格。取單位厚度的海床建立二維軸對(duì)稱平面應(yīng)變模型。由于麻坑是環(huán)形的，因此用軸對(duì)稱的幾何模型來(lái)進(jìn)行數(shù)值模擬。根據(jù)現(xiàn)場(chǎng)地質(zhì)數(shù)據(jù)，點(diǎn)劃線為對(duì)稱軸，模型長(zhǎng)×寬為100 m×150 m，整體水頭為500 m(以D點(diǎn)計(jì))，并沿深度線性增加。其中邊界AB與CD為豎向滑動(dòng)邊界，AD為自由邊界，均分布線性孔壓；邊界BC豎向位移約束，邊界上的孔壓均勻分布。土體單元采用四節(jié)點(diǎn)孔壓?jiǎn)卧狢PE4P 模擬。單元網(wǎng)格寬度為5 m，下部為規(guī)則的長(zhǎng)方形網(wǎng)格。裂縫不需要?jiǎng)澐志W(wǎng)格。

圖2 平面應(yīng)變數(shù)值模型及網(wǎng)格Fig.2 Plane strain numerical model and grid

2.2 土體本構(gòu)模型與參數(shù)

初始裂縫與水平方向的夾角β在0°～45°范圍內(nèi)[3,22]，本研究β取45°。初始裂縫長(zhǎng)度取10 m。天然氣水合物聚集區(qū)(concentrated hydrate zone，CHZ)的埋深范圍在10～800 m[22]。為了減少計(jì)算量，本研究中V 形裂縫的埋深取60 m。根據(jù)CHILLARIGE等[23]統(tǒng)計(jì)結(jié)果，海床坡角在0°～10°之間。故參數(shù)分析中取坡角θ為0°，2.5°，5.0°，7.5°和10.0°。

表1所示為數(shù)值模擬基本參數(shù)[17,24]。模型采用理想線彈性本構(gòu)，根據(jù)彈性力學(xué)可知，靜止土壓力系數(shù)K0與泊松比v的關(guān)系式如下：

表1 數(shù)值模擬基本參數(shù)表Table 1 Basic parameters of numerical simulation

根據(jù)飽和軟黏土的泊松比取值范圍，參數(shù)分析中靜止土壓力系數(shù)K0取0.50，0.55，0.60，0.65，0.70和0.80。數(shù)值模型中土體是均質(zhì)分布的，擴(kuò)展有限元模型使用摩爾-庫(kù)侖破壞準(zhǔn)則和理想線彈性本構(gòu)來(lái)描述每個(gè)土層的應(yīng)力-應(yīng)變行為。假設(shè)注入流體在裂隙內(nèi)的濾失系數(shù)很小，取1×10-9。飽和軟黏土的抗拉強(qiáng)度一般很小，可忽略不計(jì)，但為了保證數(shù)值模型的收斂性，本研究中將抗拉強(qiáng)度取2 kPa。由于缺少對(duì)海底飽和軟黏土斷裂韌度的測(cè)量，海底淺埋高壓水氣突涌數(shù)值模型的土體臨界斷裂能釋放率采用半經(jīng)驗(yàn)公式進(jìn)行推算。

2.3 數(shù)值模擬步驟

在給定邊界條件情況下，通過(guò)應(yīng)力迭代法完成斜坡海床模型的初始地應(yīng)力平衡，進(jìn)而設(shè)置高壓流體注入，開(kāi)展擴(kuò)展有限元的應(yīng)力-滲流耦合分析，模擬深海淺埋高壓水氣突涌運(yùn)移。

3 淺埋高壓水氣突涌過(guò)程中的裂縫擴(kuò)展規(guī)律

3.1 起裂和擴(kuò)展規(guī)律

圖3所示為海床中不同坡角和靜止土壓力系數(shù)K0條件下對(duì)應(yīng)的突涌起裂壓力。由圖3可見(jiàn)：在初始裂縫埋深相同的情況下，海床坡角越大，靜止土壓力系數(shù)越大，突涌發(fā)生所需要的壓力則更大。起裂壓力幾乎隨著靜止土壓力系數(shù)增大而線性增加。當(dāng)坡角等于0°，靜止土壓力系數(shù)K0從0.5增加到0.8 時(shí)，起裂壓力從510 kPa 增加到740 kPa；而當(dāng)靜止土壓力系數(shù)等于0.5 時(shí)，坡角從0°增加到10°，起裂壓力從510 kPa 增加到620 kPa。對(duì)比坡角為0°和10°工況下的初始裂縫處側(cè)向壓力可知，軟黏土中的突涌起裂壓力主要取決于初始裂縫所受到的側(cè)向壓力，約60%的起裂壓力需用于克服側(cè)向壓力，使裂縫張開(kāi)。

圖3 不同工況對(duì)應(yīng)的突涌起裂壓力Fig.3 Initation pressure of water and gas surge under different conditions

以坡角θ為0°、靜止土壓力系數(shù)K0為0.55的平坡海床中淺埋高壓水氣突涌為例，重點(diǎn)分析海床中淺埋高壓水氣遷移過(guò)程中的應(yīng)力場(chǎng)和孔壓場(chǎng)分布規(guī)律。圖4所示為高壓水氣突涌過(guò)程中土體的主應(yīng)力和超孔壓場(chǎng)分布圖。由圖4可見(jiàn)：對(duì)于平坡海床，初始地應(yīng)力場(chǎng)沿深度遞增，且最大主應(yīng)力方向均為豎直方向。由于裂縫內(nèi)存在流體驅(qū)動(dòng)，因此，裂縫的起裂和擴(kuò)展會(huì)引起周圍土體的應(yīng)力重新分布，最大主應(yīng)力方向發(fā)生旋轉(zhuǎn)。裂縫擴(kuò)展沿著裂尖單元的最大主應(yīng)力方向，裂縫兩側(cè)的土體受到流體擠壓的作用，產(chǎn)生正超孔壓，裂尖附近的土體由于受拉作用會(huì)產(chǎn)生負(fù)超孔壓。同時(shí)，最大主應(yīng)力方向逐漸變?yōu)樗椒较?，但離裂縫較遠(yuǎn)處的土體應(yīng)力并未發(fā)生改變，最大主應(yīng)力方向仍為豎直方向。由于裂縫擴(kuò)展過(guò)程中土體應(yīng)力主軸旋轉(zhuǎn)，最終淺埋高壓水氣突涌路徑和初始地應(yīng)力場(chǎng)存在夾角。

圖4 水氣突涌主應(yīng)力和超孔壓分布規(guī)律Fig.4 Principal stress and pore pressure distribution of water and gas surge

3.2 靜止土壓力系數(shù)和海床坡角對(duì)麻坑半徑的影響

圖5(a)～(d)所示為坡角為0°，靜止土壓力系數(shù)分別為0.8，0.7，0.6 和0.5 的平坡海床中淺埋高壓水氣突涌路徑和最大主應(yīng)力矢量圖。由圖5(a)～(d)可見(jiàn)：靜止土壓力系數(shù)為0.8時(shí)，裂縫的擴(kuò)展方向會(huì)更加接近水平，淺埋高壓水氣突涌路徑更長(zhǎng)，形成的麻坑半徑也更大，半徑幾乎達(dá)到100 m，此時(shí)裂縫擴(kuò)展對(duì)周圍土體應(yīng)力的影響范圍也更大。

圖5(d)～(f)所示為K0為0.5，坡角分別為0°，5°和10°這3 種工況下海床中淺埋高壓水氣突涌路徑和最大主應(yīng)力矢量圖。由圖5(d)～(f)可見(jiàn)：坡角為10°時(shí)，裂縫整體的擴(kuò)展方向更趨近于水平，導(dǎo)致水氣突涌形成的麻坑半徑更大，達(dá)到80 m左右。

圖5 水氣突涌擴(kuò)展路徑及最大主應(yīng)力矢量圖Fig.5 Propagation path and maximum principal stress vector of water gas surge

無(wú)論是平坡海床還是斜坡海床，裂縫趨向于沿裂尖的最大主應(yīng)力方向擴(kuò)展。由于裂縫內(nèi)存在流體驅(qū)動(dòng)，裂縫的起裂和擴(kuò)展會(huì)引起周圍土體的應(yīng)力重新分布，最大主應(yīng)力方向發(fā)生旋轉(zhuǎn)，并隨著裂縫擴(kuò)展逐漸趨于水平，但離裂縫較遠(yuǎn)處的土體應(yīng)力并未發(fā)生改變。由于裂縫擴(kuò)展過(guò)程中土體應(yīng)力主軸旋轉(zhuǎn)，最終導(dǎo)致淺埋高壓水氣突涌路徑和初始地應(yīng)力場(chǎng)之間存在夾角。

圖6所示為麻坑半徑與靜止土壓力系數(shù)、海床坡角的關(guān)系圖。由圖6可見(jiàn)：靜止土壓力系數(shù)越大，海床坡角越大，對(duì)應(yīng)海床的初始最大主應(yīng)力方向的傾斜程度更大，導(dǎo)致裂縫擴(kuò)展更加趨近于水平。對(duì)比坡角對(duì)于麻坑半徑的影響，可以發(fā)現(xiàn)靜止土壓力系數(shù)對(duì)麻坑半徑的影響程度更大。

圖6 麻坑半徑與靜止土壓力系數(shù)和海床坡角的關(guān)系圖Fig.6 Relationship of pockmarks radius with coefficient of earth pressure and seabed slope angle

3.3 麻坑半徑綜合評(píng)價(jià)

圖7所示為麻坑半徑的綜合評(píng)價(jià)圖。由圖7可見(jiàn)：麻坑半徑與海床坡角、靜止土壓力系數(shù)呈正相關(guān)。海床坡角越大，靜止土壓力系數(shù)越大，則形成的麻坑半徑越大，淺埋高壓水氣突涌影響的范圍越大。取坡角等于2.5°時(shí)，當(dāng)靜止土壓力系數(shù)K0從0.5 增加到0.8，麻坑半徑從72 m 增加到105 m；取K0為0.6，當(dāng)坡角從0°增加到10°，麻坑半徑從72 m增加到96 m。

圖7 麻坑半徑綜合評(píng)價(jià)圖Fig.7 Comprehensive evaluation of pockmark radius

羅敏等[25]調(diào)研發(fā)現(xiàn)海底絕大多數(shù)的麻坑直徑落在75～200 m 的區(qū)間內(nèi)。PILCHER 等[26]統(tǒng)計(jì)世界各地已觀測(cè)到的57 個(gè)麻坑區(qū)，發(fā)現(xiàn)大部分麻坑直徑處于50～270 m 的區(qū)間。數(shù)值模擬結(jié)果表明，麻坑的直徑主要在140～230 m之間，與前人的統(tǒng)計(jì)研究結(jié)論相一致。

4 結(jié)論

1)靜止土壓力系數(shù)與海床坡角影響應(yīng)力主軸旋轉(zhuǎn)，進(jìn)而影響突涌路徑。海床坡角越大，靜止土壓力系數(shù)越大，土體應(yīng)力主軸越易發(fā)生旋轉(zhuǎn)，則淺埋高壓水氣突涌形成的麻坑半徑也越大，影響范圍越大。

2)靜止土壓力系數(shù)對(duì)淺埋高壓水氣突涌路徑的影響比坡角的影響更大。當(dāng)坡角等于0°，靜止土壓力系數(shù)從0.5 增加到0.8 時(shí)，麻坑半徑從70 m增加到100 m；而當(dāng)K0為0.5，坡角從0°增加到10°時(shí)，麻坑半徑從70 m增加到80 m。

3)在初始裂縫埋深相同的情況下，起裂壓力隨著靜止土壓力系數(shù)和坡角增大而增加，其中靜止土壓力系數(shù)的影響更大。當(dāng)坡角等于0°，靜止土壓力系數(shù)K0從0.5 增加到0.8 時(shí)，起裂壓力從510 kPa 增加到740 kPa；而當(dāng)靜止土壓力系數(shù)等于0.5，坡角從0°增加到10°時(shí)，起裂壓力從510 kPa增加到620 kPa。

4)應(yīng)力-滲流耦合分析方法結(jié)合了擴(kuò)展有限元與Cohesive 黏聚力單元，能夠考慮地層原生主應(yīng)力影響，即靜土土壓力系數(shù)和海床坡角的影響，模擬出水氣突涌過(guò)程中的運(yùn)移路徑以及應(yīng)力場(chǎng)、孔壓場(chǎng)變化，實(shí)現(xiàn)水氣突涌運(yùn)移過(guò)程的模擬以及對(duì)形成的麻坑半徑進(jìn)行綜合評(píng)價(jià)。