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

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

0 引言

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

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

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

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

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

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

1 多場(chǎng)耦合模擬原理和方法

1.1 非等溫多相滲流

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

表1 TOUGH＋HYDRATE的控制方程表

1.2 力學(xué)

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

表2 FLAC3D的控制方程表

1.3 耦合方法

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

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

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

2 可靠性驗(yàn)證

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

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

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

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

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

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

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

3.1 地質(zhì)概況

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

神狐海域海底溫度介于3.3～3.7 ℃，水合物儲(chǔ)層孔隙度介于32%～55%，有效滲透率介于0.2～20 mD，水合物飽和度高達(dá)64%以上。神狐海域發(fā)育2個(gè)大型水合物礦體[35]，其在空間上厚度不均一，據(jù)地震資料推測(cè)，水合物礦體中心較厚，可達(dá)90 m，水合物富集層厚度向礦體邊緣逐漸減小，游離氣在水合物儲(chǔ)層下部富集[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 概念模型

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

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

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

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

3.3 模型參數(shù)

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

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

相對(duì)滲透率模型公式為：

毛細(xì)管壓力（pcap）模型公式為：

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

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

3.4 模擬結(jié)果分析

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

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

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

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

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

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

圖6 水平井開(kāi)采下孔隙壓力的空間演化特征圖

圖7 開(kāi)采1年時(shí)的壓差和溫度差分布圖

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

圖8 開(kāi)采1年時(shí)的水合物飽和度和氣體飽和度分布圖

3.4.3 地層變形

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

圖9 水平井開(kāi)采過(guò)程中海底（水平井中部正上方）沉降變化圖

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

圖10 水平井開(kāi)采沉降的空間示意圖

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

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

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

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

4 結(jié)論與建議

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

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

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

符 號(hào) 說(shuō) 明

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