天然氣水合物水平井降壓開采多相滲流—傳熱—力學(xué)耦合數(shù)值模擬：方法和南海場地應(yīng)用

趙小龍 王增林 趙益忠 左家強 李 鵬 梁 偉 王 冰 陳 雪 雷宏武 金光榮

1.中國石化勝利油田分公司石油工程技術(shù)研究院 2. 中國石化勝利油田分公司3. 巖土力學(xué)與工程國家重點實驗室·中國科學(xué)院武漢巖土力學(xué)研究所 4. 中國科學(xué)院廣州能源研究所

0 引言

天然氣水合物（以下簡稱水合物）廣泛賦存于陸域凍土帶和深海沉積物中，被認為是21世紀最理想的替代能源，吸引世界各國投入相關(guān)研究[1]。2020年，我國首次采用水平井技術(shù)實現(xiàn)水合物日均產(chǎn)氣的重大突破[2]，連續(xù)產(chǎn)氣30天[3]，進一步證實了泥質(zhì)粉砂水合物的開采可行性，也充分說明了水平井等新技術(shù)的應(yīng)用對于提高天然氣水合物的產(chǎn)能至關(guān)重要[4-5]。

在開采水合物時，地層中的固體水合物分解為氣體和水，地層強度降低、吸收熱量，將可能導(dǎo)致地層變形[6]，甚至引發(fā)井壁失穩(wěn)[7]、海底面沉降、海底滑坡等地質(zhì)災(zāi)害[8-9]，其開采過程實際上是含相變行為的多相滲流（Hydrodynamic）—傳熱（Thermal）—力學(xué)（Mechanics）耦合過程（THM耦合）。相比于室內(nèi)實驗和礦場試驗，數(shù)值模擬是較為經(jīng)濟靈活的研究方法，可充分研究各種開采工況下THM耦合響應(yīng)。以往水合物開采數(shù)值模擬研究多關(guān)注不同工況下的產(chǎn)能，對儲層穩(wěn)定性關(guān)注較少[10-12]。儲層穩(wěn)定性是水合物開采所面臨的關(guān)鍵問題之一，也是確保水合物安全高效開采的前提。

目前，國內(nèi)外水合物開采THM耦合數(shù)值模擬研究依然處于起步階段。傳統(tǒng)水合物開采模擬軟件（TOUGH+HYDRATE、STOMP-HYD、MH21-HYDRES和HydResSim等代表性TH兩場耦合軟件）均不能處理水合物開采過程中的力學(xué)變形問題[10]。目前已有水合物開采的THM三場耦合程序較少[13]，且多基于成熟商業(yè)力學(xué)軟件研發(fā)，其中以基于FLAC進行二次開發(fā)的最多。此外，少部分耦合程序基于自行開發(fā)力學(xué)模塊進行[14]。這些程序中刻畫力學(xué)特征的本構(gòu)模型涵蓋了非線性彈性、摩爾—庫倫模型、臨界狀態(tài)模型等彈塑性模型，其中THM耦合方式多采用部分耦合方法[15]，降低了程序的復(fù)雜程度。

Rutqvist等[16]基于TOUGH+HYDRATE程序與力學(xué)程序FLAC，開發(fā)了TOUGH+HYDRATE+FLAC程序，采用摩爾—庫倫模型表征沉積物的力學(xué)性質(zhì)，將水合物強度表達為水合物飽和度的線性函數(shù)，其隨水合物的分解而降低。在墨西哥灣場地的應(yīng)用結(jié)果表明，井周圍儲層的力學(xué)響應(yīng)受儲層內(nèi)的壓力衰減驅(qū)動；在較大的生產(chǎn)壓差下，沉降程度可達數(shù)米[17]；在相同開采時間下，水平井開采將導(dǎo)致更大的海底沉降[18]。Uchida等[19]基于FLAC中的摩爾—庫倫模型，結(jié)合試采產(chǎn)氣資料研究了加拿大Mallik場地水合物開采過程中應(yīng)力變化和儲層變形特征，評估了試采引發(fā)的地層變形。此后，Uchida等[20]基于臨界狀態(tài)力學(xué)本構(gòu)模型，研究2013年日本水合物試采地層變形的力學(xué)響應(yīng)特征[21]，并通過模型再現(xiàn)水合物開采出砂規(guī)律。Zhou等[22]同樣研究了日本水合物試采過程中的力學(xué)響應(yīng)特征。各國學(xué)者還對韓國郁陵盆地[23]、印度K-G盆地[24]等場地的水合物開采中的力學(xué)響應(yīng)特征進行研究。Sun等[25]在COMSOL上開發(fā)了考慮水合物分解的沉積物力學(xué)模擬程序，研究降壓過程中的多物理場響應(yīng)。Sun等[26]基于TOUGH+HYDRATE+FLAC耦合程序，研究了南海神狐海域水合物開采地層和井壁穩(wěn)定性，鉆井過程鉆井液入侵導(dǎo)致儲層屈服破壞，井徑增大[27]。

Jin等[14]基于TOUGH2Biot程序[28-29]，根據(jù)線彈性材料假設(shè)，開發(fā)了THM耦合程序TOUGH+hydrate+Biot（也稱hydrateBiot），該程序參與了美國國家能源技術(shù)實驗室（NETL）資助的第二次國際水合物代碼對比研究[13]，程序可靠性得到驗證。Jin等[30]利用該程序?qū)ξ覈虾Ｉ窈Ｓ蛩衔镩_采過程中沉降情況進行研究，結(jié)果表明，水平井開采導(dǎo)致的沉降大于豎直井開采，最大沉降出現(xiàn)在開采射孔段上方，井筒附近側(cè)向變形最大，并研究了多井聯(lián)合開采下的產(chǎn)氣動態(tài)和力學(xué)變形耦合特征。

THM耦合數(shù)值模擬研究加深了對水合物開采過程力學(xué)響應(yīng)規(guī)律的認識[31]，但以往對水平井開采的模型均被簡化為垂直于生產(chǎn)井的剖面來進行模擬，對水平井開發(fā)下真實地層變形空間的演化特征刻畫不足。因此，筆者基于TOUGH+HYDRATE和FLAC程序，開發(fā)更適于水合物礦藏水平井開采的THM耦合程序，采用解析解驗證程序正確性。以我國南海第二次試采產(chǎn)氣量約束模型，研究水平井中長期開采情景下的產(chǎn)氣特征和溫壓演化規(guī)律，分析地層變形和海底沉降演化規(guī)律。

1 多場耦合模擬原理和方法

1.1 非等溫多相滲流

水合物開采涉及復(fù)雜的非等溫多相滲流、傳熱、相變等過程，TOUGH+HYDRATE考慮了四相（水、氣、冰和水合物）四組分（水合物、甲烷、水及鹽等水溶性抑制劑）的水合物生成和分解過程，能模擬復(fù)雜地質(zhì)介質(zhì)中甲烷水合物的非等溫水合反應(yīng)、相態(tài)變化、流體運移和熱量傳遞等過程，可用于降壓法、熱激法、抑制劑注入及其聯(lián)合方法的水合物開采數(shù)值模擬。TOUGH+HYDRATE中通用的質(zhì)量和能量守恒方程見表1[32]。

表1 TOUGH＋HYDRATE的控制方程表

1.2 力學(xué)

力學(xué)計算基于FLAC3D程序二次開發(fā)，F(xiàn)LAC3D是Itasca公司研發(fā)的能進行巖土材料熱流固耦合分析的軟件，廣泛應(yīng)用于巖土、采礦工程分析中。FLAC3D的力學(xué)和熱力學(xué)模塊的控制方程見表2。FLAC3D以節(jié)點計算為對象，力和質(zhì)量集中在節(jié)點上，通過運動方程在時域內(nèi)利用有限差分法和混合離散技術(shù)求解。為保證計算精度和穩(wěn)定性，使用足夠小的時間步來完成運動方程和應(yīng)力—應(yīng)變本構(gòu)方程之間的顯式計算迭代。在一個時間步中，首先調(diào)用運動方程，基于上一時間步的節(jié)點速度和力更新節(jié)點速度；然后將節(jié)點速度用于計算新的應(yīng)變速度和應(yīng)力。對于應(yīng)力—應(yīng)變本構(gòu)關(guān)系，筆者基于FLAC3D內(nèi)置的摩爾—庫倫本構(gòu)模型二次開發(fā)，將水合物沉積強度參數(shù)（體積模量、剪切模量、內(nèi)聚力）表達為水合物飽和度的線性關(guān)系[14,16]。

表2 FLAC3D的控制方程表

1.3 耦合方法

多相滲流（H）、傳熱（T）和力學(xué)（M）三場耦合主要分為兩種方法：全耦合法和部分耦合法。全耦合法需要同時求解幾類方程，由于這幾類方程的特征差別較大，在地下流動系統(tǒng)多場耦合模擬中并不多見，全耦合法一般僅用于傳熱和多相滲流兩場的耦合。部分耦合法即在一個時間步長里按照一定的順序序列計算，計算完成該時間步長后考慮其對系統(tǒng)的影響，然后進行下一步時間步長計算，如此反復(fù)直到計算停止。本文TOUGH+HYDRATE中傳熱和多相滲流采用全耦合法，力學(xué)計算采用部分耦合法（圖1-a），步驟如下：①采用TOUGH+HYDRATE計算出流體壓力、溫度和水合物飽和度；②傳遞給FLAC3D進行力學(xué)計算（包括更新力學(xué)參數(shù)與水合物飽和度的關(guān)系）；③更新孔隙度和滲透率，傳遞給TOUGH+HYDRATE進行下一步的計算，如此反復(fù)進行。

圖1 TOUGH+HYRATE與FLAC3D耦合及其調(diào)用邏輯圖

TOUGH+HYDRATE+FLAC3D模擬器的主程序為TOUGH+HYDRATE，基于FLAC3D開發(fā)的力學(xué)計算模塊在TOUGH+HYDRATE內(nèi)部中調(diào)用。調(diào)用分為3個階段：初始化階段、計算階段和結(jié)束階段。在初始化階段，讀入與TOUGH+HYDRATE對應(yīng)的力學(xué)計算網(wǎng)格數(shù)據(jù)，設(shè)置力學(xué)模型和參數(shù)，設(shè)置輸出；在計算階段，TOUGH+HYDRATE單步內(nèi)計算完畢后把結(jié)果通過文件形式傳遞給FLAC3D進行力學(xué)計算，并把結(jié)果寫出文件供TOUGH+HYDRATE下一步計算使用；結(jié)束階段輸出時間變化的結(jié)果。具體調(diào)用邏輯見圖1-b。

2 可靠性驗證

2.1 滲流—力學(xué)（HM）耦合

為驗證TOUGH+HYDRATE+FLAC3D耦合程序的正確性，筆者選用Terzaghi一維固結(jié)模型來對比，其解析解參見本文參考文獻[29]。該模型描述了土層頂部堆載后垂向上發(fā)生的排水固結(jié)過程。由圖2可見，數(shù)值解和解析解結(jié)果幾乎一致，驗證了程序的正確性。

圖2 一維固結(jié)模型沉降模擬結(jié)果和理論解對比圖

2.2 熱—力學(xué)（TM）耦合

為驗證程序傳熱—力學(xué)耦合正確性，筆者選取一維熱傳導(dǎo)引起的變形問題[29]進行對比，具體解析解參見本文參考文獻[29]。該模型描述了地表施加低溫后土層在傳熱作用下發(fā)生的垂向變形過程。由圖3可見，數(shù)值解和解析解結(jié)果幾乎一致，表明程序可以刻畫水合物分解吸熱過程引發(fā)的變形過程。

圖3 一維熱傳導(dǎo)數(shù)值模擬結(jié)果和理論解對比圖

3 南海場地應(yīng)用

3.1 地質(zhì)概況

神狐海域位于南海北部珠江口盆地白云凹陷，水深 1 000 ～ 1 500 m[33]，是我國重要的海域水合物賦存區(qū)，也是我國兩次水合物試采的試驗區(qū)。神狐海域沉積地層自老到新可劃分為古新統(tǒng)神狐組、始新統(tǒng)文昌組及下漸新統(tǒng)恩平組陸相地層，晚漸新統(tǒng)珠海組，中新統(tǒng)珠江組、韓江組、粵海組，上新統(tǒng)萬山組和第四系瓊海組。文昌組、恩平組發(fā)育成熟—高成熟烴源巖，可為水合物生成提供氣源。區(qū)內(nèi)流體底辟、氣煙囪廣泛發(fā)育，與高角度斷裂溝通深部烴源巖，成為水合物富集成藏的流體運移通道[34]。

神狐海域海底溫度介于3.3～3.7 ℃，水合物儲層孔隙度介于32%～55%，有效滲透率介于0.2～20 mD，水合物飽和度高達64%以上。神狐海域發(fā)育2個大型水合物礦體[35]，其在空間上厚度不均一，據(jù)地震資料推測，水合物礦體中心較厚，可達90 m，水合物富集層厚度向礦體邊緣逐漸減小，游離氣在水合物儲層下部富集[36]，部分區(qū)域可能以三相混合層存在。第二次試采先導(dǎo)孔GMGS6—SH02資料顯示[3,37]，水合物層厚45.6 m，平均有效孔隙度為37.3%，平均水合物飽和度為31%，平均滲透率為2 mD。三相混合層厚24.6 m，平均有效孔隙度為34.6%，平均水合物飽和度為11.7%，氣體飽和度為13.2%，平均滲透率為6.63 mD。游離氣層厚19 m，平均有效孔隙度為34.7%，平均水合物飽和度為7.3%，平均滲透率為6.8 mD。

3.2 概念模型

將南海第二次試采場地概化水平延伸地層，模型頂面對應(yīng)于海底面（圖4）。模型模擬范圍為1 223.6 m×1 071 m×486 m 的區(qū)域。模型垂向上是由頂板（厚213 m）、水合物儲層（厚27 m）、游離氣層（厚33 m）、底板（厚213 m）等地層構(gòu)成，垂向總厚度為486 m。水平井埋深為241.5 m，位于游離氣層頂部；水平井在平面上位于模型中部，并沿y方向延伸，總長度為300 m；y方向上的模型邊界距離水平井兩端各385.5 m。在x方向上，模型以水平井井筒為中心各側(cè)向延伸611.8 m。

圖4 南海水合物水平井開采概念模型圖（僅繪制了x正方向的一半）

模型的初始壓力場服從靜水壓力，其水合物儲層底部的初始壓力為14.7 MPa。模型的初始溫度場根據(jù)水合物儲層底部靜水壓力和相平衡條件計算得相平衡溫度為14.80 ℃，然后根據(jù)地溫梯度4.02 ℃/100 m推測得到。對于水合物儲層和游離氣層，通過將溫壓平衡后的初始水飽和狀態(tài)分別更改為水合物與水兩相和氣水兩相共存條件并利用軟件自平衡得到。模型側(cè)向為非滲透絕熱邊界，模型頂?shù)诪槎囟▔哼吔?。假設(shè)模型受等向的初始應(yīng)力作用，即初始時水平和垂向的初始應(yīng)力相等[16]，初始應(yīng)力以25.48 kPa/m隨深度增加[26]。對于力學(xué)邊界，模型側(cè)面和底面為法向方向為0的位移邊界，模型頂面為自由位移邊界，可在開采過程中自由下沉。

模型在x方向上，水平井處網(wǎng)格大小為0.1 m，此后網(wǎng)格大小隨距水平井筒距離而呈倍數(shù)增加，增長因子為1.2，最大為102 m，共剖分77份。y方向上，300 m水平井被剖分為30份，其余網(wǎng)格的大小隨距水平井的距離而呈倍數(shù)增加，增長因子為1.2，最大為74 m，共剖分52份。z方向上，水合物儲層和游離氣層的網(wǎng)格大小均為3 m，頂?shù)装宓木W(wǎng)格大小隨之增加，最大為頂?shù)椎?8 m，共剖分48份。整個模型被剖分為 77×52×48 = 192 192 個網(wǎng)格。

3.3 模型參數(shù)

模型所用參數(shù)列于表3[3,37]中，所有地層的孔隙度為0.32。儲層的水合物飽和度為0.34，滲透率為2 mD；游離氣層的孔隙度為0.32，游離氣飽和度為0.25，滲透率為6 mD；頂?shù)装宓貙拥臐B透率均為6 mD。模型其他水熱物性、相對滲透率、毛細管壓力參數(shù)見表3。

表3 模型所用水熱物性參數(shù)表

相對滲透率模型公式為：

毛細管壓力（pcap）模型公式為：

我國近年來獲得較為豐富的南海神狐海域沉積物力學(xué)實驗數(shù)據(jù)[12,38-42]，其與日本海洋水合物場地較為相近。本次研究所采用海洋沉積物力學(xué)數(shù)據(jù)見表4[17,26]。

表4 模型所用力學(xué)參數(shù)表

3.4 模擬結(jié)果分析

3.4.1 氣水產(chǎn)出特征

圖5顯示了水平井開采的氣水產(chǎn)出速率和產(chǎn)出量。圖5-a顯示，模型計算獲得的30天累計產(chǎn)氣量為81×104m3，這與第二次試采工程的產(chǎn)氣量（30天為86×104m3）相當，表明了本模型概化的可靠性。開采1年后，累計產(chǎn)氣量可達1 190×104m3，累計產(chǎn)水量可達 3 450×104m3，氣水比（產(chǎn)氣量與產(chǎn)水量之比）為0.35。整個開采過程中的累計產(chǎn)氣量和累計產(chǎn)水量隨時間呈線性增加，氣水比變化不大。

圖5 水平井開采的氣水產(chǎn)出動態(tài)特征圖

從圖5-b中可見，水平井降壓立即導(dǎo)致生產(chǎn)井附近游離氣和地層水迅速流入生產(chǎn)井，初始時刻的產(chǎn)氣速率和產(chǎn)水速率較高，可達4.5×104m3/d和700 m3/d。高生產(chǎn)速率導(dǎo)致水平井周游離氣濃度降低，產(chǎn)氣速率迅速降低并達到最低點；當較遠位置的游離氣或水合物分解氣運移到達生產(chǎn)井后，產(chǎn)氣速率逐漸回升。產(chǎn)水速率呈現(xiàn)快速下降趨勢并逐漸趨于穩(wěn)定，表明地層中逐漸建立了較為穩(wěn)定的動態(tài)壓力平衡。在開采約50天后，產(chǎn)水速率穩(wěn)定，其產(chǎn)氣速率在短暫回升狀態(tài)，也逐漸趨于穩(wěn)定。可見，產(chǎn)氣速率的回升得益于地層中壓降范圍的擴展。但隨著地層中游離氣的不斷產(chǎn)出，地層中的游離氣含量不斷減少，產(chǎn)氣速率在約250天后出現(xiàn)輕微下降趨勢。

3.4.2 流體參數(shù)空間分布演化特征

圖6顯示了沿水平井延伸方向上的孔隙壓力分布特征，水平井降壓導(dǎo)致生產(chǎn)井周孔隙壓力降低，低壓區(qū)域隨開采進行而不斷向外擴展；結(jié)合50天后表現(xiàn)穩(wěn)定的產(chǎn)水速率和井周壓力擴展范圍不明顯，這間接證明孔隙壓力已經(jīng)達到了動態(tài)平衡狀態(tài)。圖7-a、c清晰地顯示，在開采1年后，壓降的影響范圍已經(jīng)穿過游離氣層，到達底板；而由于水合物儲層有效滲透率較低，壓降影響范圍向上穿切作用較弱，未明顯穿過水合物儲層頂部。側(cè)向上，沿著（圖7-a）和垂直于（圖7-c）水平井剖面上的壓降影響范圍并不相同，垂直于水平井的剖面上壓降影響范圍較大。

圖6 水平井開采下孔隙壓力的空間演化特征圖

圖7 開采1年時的壓差和溫度差分布圖

由圖8-a可知，水平井降壓僅引起水合物儲層底部的小部分水合物分解。水合物分解范圍很小，分解量較少，這是因為水合物存在導(dǎo)致儲層有效滲透率較低，故降壓易導(dǎo)致氣層中流體的流動。水合物分解導(dǎo)致水平井周圍的溫度降低（圖7-b、d），游離氣層溫度受其頂部降溫傳熱影響而降低。水合物開采導(dǎo)致游離氣層中的氣體向水平井聚積（圖8-b），而水平井中部聚體的氣體大于水平井兩端。開采1年后，游離氣層底部的氣體飽和度降低，表明該位置處的氣層出現(xiàn)虧空趨勢，這可能與水平井產(chǎn)氣速率呈現(xiàn)的降低趨勢有關(guān)。水合物儲層底部水合物分解的氣體能很快被抽出，故水合物儲層底部未出現(xiàn)明顯的氣體積聚。

圖8 開采1年時的水合物飽和度和氣體飽和度分布圖

3.4.3 地層變形

圖9顯示了海底沉降的變化情況，水平井降壓開始后，海底呈線性下沉。整個沉降速率可以53天劃分為兩個階段，沉降速率的改變與產(chǎn)水速率變化具有一定的相關(guān)性，這是因為產(chǎn)水速率變化直接反應(yīng)整體系統(tǒng)中孔隙壓力是否會進入壓力動態(tài)平衡狀態(tài)。在本模型采用的力學(xué)參數(shù)下，開采60天預(yù)測得到海底沉降為0.16 m，而開采1年后的海底沉降量達到約 0.52 m。

圖9 水平井開采過程中海底（水平井中部正上方）沉降變化圖

由沉降的空間分布圖（圖10）可知，水平井開采導(dǎo)致了空間上大范圍的地層變形（垂向位移）。水平井上方和側(cè)向的地層發(fā)生沉降，而底部底板發(fā)生一定程度的隆起?？拷骄戏降牡貙映两底畲?，而向上延伸到海底過程中不斷減小，這表明水平井開采過程易于發(fā)生地層變形和井筒破壞的位置處于水平井及其上方附近。

圖10 水平井開采沉降的空間示意圖

開采60天時（圖10-a、c），沿水平井方向剖面上的沉降范圍大于垂直于水平井方向的剖面。故在生產(chǎn)井埋深附近，其z平面上垂向位移分布是以水平井為長軸的橢圓形；但這種疊加作用隨距水平井距離而減弱，海底面上的沉降呈近圓形分布。開采到365天時，沉降的空間范圍加大，沉降量增加，海底面上的沉降變?yōu)闄E圓形分布。

圖11-a清晰地顯示了平面位置上對應(yīng)于水平井中部和端部處不同深度上的沉降分布規(guī)律，圖中Line A 為截面x=611.8 m 和y=535.8 m 的交線，平面位置對應(yīng)水平井中部；Line B為截面x=611.8 m和y=385.8 m的交線，平面位置對應(yīng)水平井端部。開采60天時，水平井開采導(dǎo)致井周地層有效應(yīng)力增加，孔隙壓縮而發(fā)生沉降，最大沉降位置位于生產(chǎn)井上方附近，因水平井以淺產(chǎn)生的壓降受限于水合物儲層中，水合物儲層中變形最大；水平井底部由于滲透力作用而出現(xiàn)一定隆起，這是開采初期，水平井壓降影響范圍不大，但水平井降壓導(dǎo)致其上覆總應(yīng)力降低，從而引起下伏層位發(fā)生膨脹現(xiàn)象。隨開采進行，沉降量不斷加大，而其底部由于壓降擴展影響，有效應(yīng)力增大明顯，從而轉(zhuǎn)變?yōu)槌两第厔荨Ｋ骄胁康奈灰拼笥谒骄瞬?，這是由于水平中部壓降強烈，受到水平井兩端壓降的影響，壓降疊加效應(yīng)使得其垂向位移偏大。

圖11 不同深度的垂向位移圖

圖11-b顯示了沿水平井（y方向）和垂直于水平井（x方向）兩個方向上的海底沉降情況，圖中y方向上的水平距離為y方向原坐標與76 m之和，便于將最大位移中心移動到同一位置（橫坐標）進行比較。由圖11-b可見，平面上對應(yīng)于水平井降壓中心處的海底沉降最大，且兩個方向上的位移相近。海底沿水平井延伸方向的海底沉降大于垂直于水平井延伸方向，且其差值隨開采進行而不斷增大。這是由于沿水平井和垂直于水平井兩個方向的垂向位移疊加效應(yīng)的結(jié)果。

4 結(jié)論與建議

1）基于多相滲流—傳熱—力學(xué)多場耦合理論，開發(fā)了TOUGH＋HYDRATE和FLAC耦合程序，通過一維排水固結(jié)和熱傳導(dǎo)變形解析解驗證耦合程序的正確性。

2）基于我國南海第二次試采產(chǎn)氣量約束水平井開采模型（開采壓力7 MPa、井長300 m），得到水平井開采的氣水產(chǎn)出動態(tài)特征及其壓力等空間演化特征。經(jīng)預(yù)測，開采1年后的產(chǎn)氣量可達1 190×104m3，壓降在游離氣層中的影響范圍較大，產(chǎn)氣來源于游離氣層和水合物層。

3）通過THM耦合模型，基于海洋沉積物力學(xué)實驗數(shù)據(jù)，預(yù)測試采60天時的海底沉降約為0.16 m，開采1年后可達到約0.52 m，增加幅度超過3倍，水合物中長期開發(fā)過程中的海底沉降問題不可忽視。海底沉降呈線性增加，空間上的沉降呈橢圓形分布。水平井及其上方的變形（垂向位移）最大，是最易發(fā)生地層破壞的位置。后期建議結(jié)合海底地形、峽谷展布、礦體空間厚度不均等特征，研究中長期開采下的海底變形和滑坡風(fēng)險。

符 號 說 明

V表示單元體積；m3；Γ表示單元表面積，m2；Mk表示組分k的質(zhì)量或能量累積項，kg/m3或J/m3；Fk表示組分k的質(zhì)量或能量通量，kg/(m2·s)或J/(m2·s)；n表示單位法向量；qk表示組分k的質(zhì)量或能量源匯項，kg/(m3·s)或J/(m3·s)；t表示時間，s；β＝A、G、H、I分別表示液相、氣相、水合物相和冰相；k＝w、m、h、i分別表示組分水、甲烷、水合物和鹽等抑制劑；φ表示孔隙度，無量綱；Sβ表示β相的飽和度，無量綱；ρβ表示β相的密度，kg/m3；Xβk表示β相中組分k的質(zhì)量分數(shù)，無量綱；θ表示能量組分；ρR表示巖石骨架密度，kg/m3；CR表示巖石骨架比熱容，J/(kg·K)；T表示溫度，℃；Uβ表示β相的內(nèi)能，J/kg；Qdiss表示水合物相變能量，J/m3；K表示地層固有滲透率（簡寫表示滲透率），m2；Krβ表示β相的相對滲透率，無量綱；μβ表示β相黏滯度，Pa·s；pβ表示β相壓力，Pa；g表示重力加速度，m/s2；λ表示含水合物地層復(fù)合熱傳導(dǎo)率，W/(m·K)；hβ表示β相的比焓，J/kg；ρm表示巖土材料或沉積物的平均密度，kg/m3；v表示節(jié)點速度矢量，m/s；H表示指定材料的函數(shù)；σ'表示有效應(yīng)力，Pa；σ表示總應(yīng)力，Pa；表示無窮小應(yīng)變速率張量，無量綱；ε表示無窮小應(yīng)變張量，無量綱；I表示單位張量；α表示Biot系數(shù)，無量綱；p表示孔隙流體壓力，Pa；u表示位移矢量，m；Δεe表示彈性應(yīng)變增量；Δεp表示塑性應(yīng)變增量；ΔεT表示熱應(yīng)變增量；βT表示熱膨脹系數(shù)，1/K。Fil表示l節(jié)點在i方向的不平衡力分量，N；ml表示l節(jié)點的集中質(zhì)量，kg；KrA、KrG分別表示液相、氣相相對滲透率，無量綱；SA*、SG*分別表示歸一化液相、氣相飽和度，無量綱。