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

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

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

0 引言

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

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

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

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

Jin等[14]基于TOUGH2Biot程序[28-29]，根據線彈性材料假設，開發了THM耦合程序TOUGH+hydrate+Biot（也稱hydrateBiot），該程序參與了美國國家能源技術實驗室（NETL）資助的第二次國際水合物代碼對比研究[13]，程序可靠性得到驗證。Jin等[30]利用該程序對我國南海神狐海域水合物開采過程中沉降情況進行研究，結果表明，水平井開采導致的沉降大于豎直井開采，最大沉降出現在開采射孔段上方，井筒附近側向變形最大，并研究了多井聯合開采下的產氣動態和力學變形耦合特征。

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

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

1.1 非等溫多相滲流

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

表1 TOUGH＋HYDRATE的控制方程表

1.2 力學

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

表2 FLAC3D的控制方程表

1.3 耦合方法

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

圖1 TOUGH+HYRATE與FLAC3D耦合及其調用邏輯圖

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

2 可靠性驗證

2.1 滲流—力學（HM）耦合

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

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

2.2 熱—力學（TM）耦合

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

圖3 一維熱傳導數值模擬結果和理論解對比圖

3 南海場地應用

3.1 地質概況

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

神狐海域海底溫度介于3.3～3.7 ℃，水合物儲層孔隙度介于32%～55%，有效滲透率介于0.2～20 mD，水合物飽和度高達64%以上。神狐海域發育2個大型水合物礦體[35]，其在空間上厚度不均一，據地震資料推測，水合物礦體中心較厚，可達90 m，水合物富集層厚度向礦體邊緣逐漸減小，游離氣在水合物儲層下部富集[36]，部分區域可能以三相混合層存在。第二次試采先導孔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 概念模型

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

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

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

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

3.3 模型參數

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

表3 模型所用水熱物性參數表

相對滲透率模型公式為：

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

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

表4 模型所用力學參數表

3.4 模擬結果分析

3.4.1 氣水產出特征

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

圖5 水平井開采的氣水產出動態特征圖

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

3.4.2 流體參數空間分布演化特征

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

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

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

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

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

3.4.3 地層變形

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

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

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

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

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

圖11-a清晰地顯示了平面位置上對應于水平井中部和端部處不同深度上的沉降分布規律，圖中Line A 為截面x=611.8 m 和y=535.8 m 的交線，平面位置對應水平井中部；Line B為截面x=611.8 m和y=385.8 m的交線，平面位置對應水平井端部。開采60天時，水平井開采導致井周地層有效應力增加，孔隙壓縮而發生沉降，最大沉降位置位于生產井上方附近，因水平井以淺產生的壓降受限于水合物儲層中，水合物儲層中變形最大；水平井底部由于滲透力作用而出現一定隆起，這是開采初期，水平井壓降影響范圍不大，但水平井降壓導致其上覆總應力降低，從而引起下伏層位發生膨脹現象。隨開采進行，沉降量不斷加大，而其底部由于壓降擴展影響，有效應力增大明顯，從而轉變為沉降趨勢。水平井中部的位移大于水平井端部，這是由于水平中部壓降強烈，受到水平井兩端壓降的影響，壓降疊加效應使得其垂向位移偏大。

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

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

4 結論與建議

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

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

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

符 號 說 明

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