海底水氣突涌運移規律與麻坑半徑定量評價

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

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

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

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

盡管上述研究為理解II型麻坑的形態發展提供了一條重要的途徑，但是它們忽略了斜坡海床中原生主應力偏轉對于水氣遷移路徑的影響，無法準確預測現場麻坑的位置及半徑，且采用傳統有限元法等效模擬裂縫擴展會高估突涌路徑的寬度。裂縫在土中沿裂尖的最大主應力方向擴展，而海床坡角和靜止土壓力系數等因素會對原生主應力偏轉產生顯著影響，進而影響水氣突涌路徑以及麻坑半徑。因此，在預測淺埋高壓水氣突涌形成II型麻坑的擴展路徑時，必須要考慮應力主軸旋轉這個重要因素。

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

1 數值原理

1.1 應力-滲流耦合基本理論

飽和軟黏土中的高壓劈裂過程會影響土體的應力場，因此，數值模擬方法采用基于有效應力原理的應力-滲流耦合分析方法[11-12]。完全飽和的土層應力平衡方程可用虛功原理表達：一定時間內土體上的虛功等于作用在土體上的體力和面力所產生的虛功，土體的平衡方程為

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

對于不可壓縮流體滲透通過飽和土體單元，根據連續性條件和質量守恒定律可知，單位時間內某一土體單元內流體質量的變化量與該單元流進與流出的流體質量之差相等，連續性方程可以寫成

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

1.2 擴展有限元方法

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

1.3 黏聚力單元區本構行為：裂紋萌生與擴展

模擬飽和軟黏土的劈裂破壞時，Cohesive單元往往會受到3個方向的應力，采用B-K斷裂準則來計算復合型裂縫的總臨界應變能釋放率Gc，表達式如下：

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

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

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

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

JOHNSON 等[18]基于加拿大新斯科舍省3 個場地的現場探頭測量結果，得到飽和細粒沉積物的拉伸斷裂韌度KIc與強度su近似呈線性關系：

式中：su為土體不排水抗剪強度。

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

根據線彈性斷裂力學(LEFM)，裂縫起裂受最大主應力控制，當最大主應力超過預定義的最大閾時，就會形成裂縫。新形成的裂縫方向垂直于裂紋尖端附近的局部最大拉應力方向。

1.4 裂縫面流體流動模型

假設Cohesive 單元內的流體是連續不可壓縮的牛頓流體，流體流動可以當成順裂縫面方向的切向流和垂直于裂縫面的法向流[19-20]。

WANG 等[21]針對不同壓縮系數的裂隙流體，開展了在巖土體中的平面應變裂縫擴展模擬試驗，發現流體的壓縮性對巖土體中裂縫的擴展影響很小。故本研究中不考慮高壓滲流破壞模擬中裂縫內流體的壓縮系數，統一按不可壓縮的水進行模擬。

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

2.1 XFEM網格與邊界條件

圖2所示為平面應變數值模型及網格。取單位厚度的海床建立二維軸對稱平面應變模型。由于麻坑是環形的，因此用軸對稱的幾何模型來進行數值模擬。根據現場地質數據，點劃線為對稱軸，模型長×寬為100 m×150 m，整體水頭為500 m(以D點計)，并沿深度線性增加。其中邊界AB與CD為豎向滑動邊界，AD為自由邊界，均分布線性孔壓；邊界BC豎向位移約束，邊界上的孔壓均勻分布。土體單元采用四節點孔壓單元CPE4P 模擬。單元網格寬度為5 m，下部為規則的長方形網格。裂縫不需要劃分網格。

圖2 平面應變數值模型及網格Fig.2 Plane strain numerical model and grid

2.2 土體本構模型與參數

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

表1所示為數值模擬基本參數[17,24]。模型采用理想線彈性本構，根據彈性力學可知，靜止土壓力系數K0與泊松比v的關系式如下：

表1 數值模擬基本參數表Table 1 Basic parameters of numerical simulation

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

2.3 數值模擬步驟

在給定邊界條件情況下，通過應力迭代法完成斜坡海床模型的初始地應力平衡，進而設置高壓流體注入，開展擴展有限元的應力-滲流耦合分析，模擬深海淺埋高壓水氣突涌運移。

3 淺埋高壓水氣突涌過程中的裂縫擴展規律

3.1 起裂和擴展規律

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

圖3 不同工況對應的突涌起裂壓力Fig.3 Initation pressure of water and gas surge under different conditions

以坡角θ為0°、靜止土壓力系數K0為0.55的平坡海床中淺埋高壓水氣突涌為例，重點分析海床中淺埋高壓水氣遷移過程中的應力場和孔壓場分布規律。圖4所示為高壓水氣突涌過程中土體的主應力和超孔壓場分布圖。由圖4可見：對于平坡海床，初始地應力場沿深度遞增，且最大主應力方向均為豎直方向。由于裂縫內存在流體驅動，因此，裂縫的起裂和擴展會引起周圍土體的應力重新分布，最大主應力方向發生旋轉。裂縫擴展沿著裂尖單元的最大主應力方向，裂縫兩側的土體受到流體擠壓的作用，產生正超孔壓，裂尖附近的土體由于受拉作用會產生負超孔壓。同時，最大主應力方向逐漸變為水平方向，但離裂縫較遠處的土體應力并未發生改變，最大主應力方向仍為豎直方向。由于裂縫擴展過程中土體應力主軸旋轉，最終淺埋高壓水氣突涌路徑和初始地應力場存在夾角。

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

3.2 靜止土壓力系數和海床坡角對麻坑半徑的影響

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

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

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

無論是平坡海床還是斜坡海床，裂縫趨向于沿裂尖的最大主應力方向擴展。由于裂縫內存在流體驅動，裂縫的起裂和擴展會引起周圍土體的應力重新分布，最大主應力方向發生旋轉，并隨著裂縫擴展逐漸趨于水平，但離裂縫較遠處的土體應力并未發生改變。由于裂縫擴展過程中土體應力主軸旋轉，最終導致淺埋高壓水氣突涌路徑和初始地應力場之間存在夾角。

圖6所示為麻坑半徑與靜止土壓力系數、海床坡角的關系圖。由圖6可見：靜止土壓力系數越大，海床坡角越大，對應海床的初始最大主應力方向的傾斜程度更大，導致裂縫擴展更加趨近于水平。對比坡角對于麻坑半徑的影響，可以發現靜止土壓力系數對麻坑半徑的影響程度更大。

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

3.3 麻坑半徑綜合評價

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

圖7 麻坑半徑綜合評價圖Fig.7 Comprehensive evaluation of pockmark radius

羅敏等[25]調研發現海底絕大多數的麻坑直徑落在75～200 m 的區間內。PILCHER 等[26]統計世界各地已觀測到的57 個麻坑區，發現大部分麻坑直徑處于50～270 m 的區間。數值模擬結果表明，麻坑的直徑主要在140～230 m之間，與前人的統計研究結論相一致。

4 結論

1)靜止土壓力系數與海床坡角影響應力主軸旋轉，進而影響突涌路徑。海床坡角越大，靜止土壓力系數越大，土體應力主軸越易發生旋轉，則淺埋高壓水氣突涌形成的麻坑半徑也越大，影響范圍越大。

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

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

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