南海神狐海域天然氣水合物降壓開采過程中儲層的穩定性

萬義釗 吳能友 胡高偉 辛 欣 金光榮 劉昌嶺 陳 強

1.青島海洋地質研究所國土資源部天然氣水合物重點實驗室2.青島海洋科學與技術國家實驗室海洋礦產資源評價與探測技術功能實驗室3.吉林大學環境與資源學院 4.中國科學院廣州能源研究所

0 引言

天然氣水合物（以下簡稱水合物）廣泛存在于海底的沉積物和陸地多年凍土帶中，分布范圍廣、資源量巨大、能量密度高，有望成為滿足未來人類能源需求的高效清潔能源[1]。目前，水合物開采的方法主要有降壓法、注熱法、注化學試劑法和二氧化碳置換法等。國際上已經有少數國家[2-6]開展了水合物試采，使用降壓法開采海洋水合物的國家如日本、中國[6]，其中產量最高的是我國于2017年5月10實施的海域水合物試采，累積產氣量為30.9×104m3，平均日產氣量為5 151 m3[6]，本次試采穩定產氣60 d。

總結上述國際天然氣水合物歷次試采經歷可知，總體上水合物的開采效率仍然較低；降壓法是最行之有效的方法。水合物降壓開采通過降低井底壓力在儲層中形成壓降漏斗，當儲層壓力降低到水合物相平衡壓力以下時，水合物開始分解。理論上，井底壓力降低幅度越大，壓降漏斗影響范圍越大，產氣速率就越大[7]，但儲層壓力降低會導致儲層有效應力增大和垂向變形，同時由于海洋天然氣水合物儲層膠結差，強度低，儲層應力增大可能引起儲層失穩破壞，而水合物又在儲層沉積物顆粒間起膠結作用，降壓引起的水合物分解會降低儲層的強度，進一步加大了儲層失穩的風險。因此，儲層穩定性是水合物開采面臨的關鍵問題，是確保水合物開采安全高效的前提。

據最新勘探結果顯示，我國南海北部神狐海域水合物儲層以黏土質粉砂和粉砂質黏土為主，儲層沉積物粒徑總體低于20 μm，是典型的孔隙充填型水合物儲層[8-9]，儲層原位測試滲透率低（＜10 mD）。許多學者利用數值模擬的手段研究神狐海域水合物的開采方法和開采潛力。Zhang等[10-11]建立了神狐海域水合物降壓+注熱開采方法的單一水平井和雙水平井模型，對其產氣能力進行了評價。李剛等[12-14]利用注熱結合降壓以及熱吞吐等方法對神狐海域水合物開采潛力進行了研究，討論了不同井型的開采效果。蘇正等[15]也對神狐海域直井熱激發的開發方法進行了評價。這些數值模擬研究對神狐海域水合物的直井、水平井、雙水平井等開采井型進行了詳細分析，對降壓、注熱、熱吞吐等開采方式進行了模擬。

國內外的學者對水合物開采的儲層變形和破壞也開展了初步的研究。沈海超等[16]將水合物分解效應融合到滲流場與巖土變形場中，建立水合物開采的流固耦合數學模型，對降壓開采過程中近井地帶的儲層穩定性進行了數值模擬。程家望等[17]將儲層沉降和井壁穩定性結合到降壓開采過程中，建立了一維的水合物降壓開采穩定性數學模型，模型不考慮傳熱過程的影響。孫可明等[18]建立了反映水合物分解引起儲層力學性質劣化的熱流固耦合模型，利用ABAQUS二次開發了模型的求解程序，研究了加熱法開采水合物儲層變形破壞規律。

然而，目前針對神狐海域黏土質粉砂和粉砂質黏土水合物儲層開采過程中的力學穩定性問題研究很少[19]。本文基于南海神狐海域水合物的鉆探資料，建立三維的水合物降壓開采的地質模型，綜合考慮儲層中水合物的分解、傳熱、滲流和骨架固體變形的多場耦合過程，建立水合物開采儲層穩定性分析的數學模型及有限元求解方法，獲取水合物降壓開采過程中儲層壓力、溫度、飽和度和應力的時空演化特征，分析神狐海域水合物降壓開采過程中儲層沉降、應力分布和穩定性。

1 物理模型

2015年9月，中國地質調查局在我國南海北部陸坡神狐海域完成第3次海域天然氣水合物鉆探航次（GMGS3）。本次鉆探的GMGS3-W19站位水深為1 273.9 m，確定海底以下135～170 m范圍內存在厚度約為35 m的水合物層。根據鉆探資料，建立了如圖1所示的物理模型。水平方向上，模型以井為中心向x方向和y方向延伸400 m。模型的頂面為海底面；水合物賦存于135 m以下，厚度35 m；水合物層上部是135 m的上覆層，底部為厚94 m的下伏層，均不含水合物。打開井段為135～162 m，即打開水合物儲層頂部向下的28 m[20]。

圖1 GMGS3-W19站位模型示意圖

2 數學模型及求解

2.1 基本假設

水合物開采是一個復雜的傳熱傳質過程，包括多相流體在多孔介質中的滲流、熱對流和熱傳導、水合物分解的化學反應以及含水合物沉積物的固體變形。為建立描述水合物開采物理過程的數學模型，需做如下假設：①假設水合物為純I型的水合物，水合物中氣體為100%甲烷，且忽略冰的生成；②在選取的控制體內，保持局部熱力平衡，即沉積物固體顆粒和流體（氣體、液體）的溫度相同；③氣和水在多孔介質中的流動符合Darcy定律；④水合物屬于孔隙充填型，附著在沉積物顆粒表面，且水合物和沉積物顆粒組成連續的復合固體材料，共同受力，其土力學特性為線彈性。

2.2 數學模型

基于以上假設建立綜合考慮熱場、氣水兩相滲流場、沉積物固體變形場和水合物分解的化學場的熱—流—固—化（THMC）的四場耦合模型。

2.2.1 水合物分解動力學的化學場

水合物分解的質量守恒方程：

基于Kim-Bishoni的動力學模型[21]，水合物分解時的產氣速率為：

相應地，水合物分解的產水速率和水合物消耗速率為：

2.2.2 氣水兩相滲流場

水合物沉積物中氣和水在多孔介質中的流動可用連續方程和Darcy定律表示，最終得到氣水兩相滲流的模型方程。

氣相：

水相：

方程（5）、（6）右端項中的mg、mw由水合物分解動力學模型計算，這兩項是耦合滲流場與化學場的關鍵。

2.2.3 熱場

水合物分解過程中的熱場可以由能量守恒方程來表示。考慮儲層的熱傳導和熱對流的水合物分解時能量守恒方程為：

其中

2.2.4 含水合物沉積物固體變形力場

水合物—沉積物復合固體變形的平衡微分方程為：

根據Terzaghi有效應力原理，土體的總應力等于水合物和沉積物顆粒骨架組成的復合固體的有效應力與孔隙壓力之和[22]，則有：

其中

根據線彈性假設，復合固體的應力應變關系為：

其中

2.2.5 輔助方程

除各物理場的控制方程外，還需各物理參數的輔助方程，主要有：

2.2.5.1 絕對滲透率方程[23]

2.2.5.2 毛細管力和氣水兩相相對滲透率方程[24]

其中

2.2.5.3 水合物相平衡壓力方程[25]

2.2.5.4 水合物相變潛熱方程[26]

2.2.5.5 水合物—沉積物顆粒復合固體力學性質方程

三軸實驗結果表明：水合物的存在會增大含水合物沉積物的強度，但水合物對泊松比的影響很小[27]，因此，假設泊松比為常數，彈性模量Em由Santamarina和Ruppel[28]提出的公式來表示：

以上方程構成了水合物開采過程的四場耦合模型。

2.3 有限元求解

經過分析，THMC的四場耦合模型可以分為兩個子系統：包含傳熱傳質的流動系統和固體變形系統。用有限元方法分別對兩個子系統求解。

對流動系統，選取pg、Sw、Sh、T作為獨立的求解變量。其中Sh可以利用方程（1）和方程（4）直接求得，現推導pg、Sw、T有限元方程。

2.3.1 pg的有限元方程

將方程（5）的第一項展開有：

利用氣體狀態方程：

將方程（5）中的密度替換為壓力，可得氣相壓力方程為：

對上式乘以δpg并積分，并用高斯公式降階可得：

將壓力在單元上用插值函數表示：

代入可得最終氣相壓力的有限元求解方程：

其中

2.3.2 Sw的有限元方程

利用毛細管力方程可得：

將方程（6）中水相壓力用氣相壓力和水相飽和度表示，則有：

上式乘以δpw并積分，降階并利用高斯公式有：

將飽和度在單元上用插值函數表示：

代入可得最終飽和度的有限元求解方程：

其中

2.3.3 T的有限元方程

方程（7）乘以δT并積分可得：

其中

第二個子系統是固體變形模型，其控制方程是式（14）。以x方向為例，對x方向的平衡微分方程乘以δux并積分，降階可得：

由插值函數可得：

可得x方向的有限元方程為：

其中

同理可得y方向和z方向的有限元單元方程。

2.4 模型網格劃分

采用四面體非結構網格對圖1所示的物理模型進行網格劃分，其網格圖如圖2所示。為提高計算精度，在水合物層中加密網格，網格總數為38 068。

2.5 初始條件和邊界條件

圖2 神狐海域水合物開采模型網格圖

根據GMGS3-W19站位的調查結果，海底面的溫度為3.75 ℃，地溫梯度為0.045 ℃/m，按儲層深度折算，儲層初始溫度在縱向上線性分布。地層初始時刻的孔隙壓力隨深度逐漸增加，符合靜水壓力平衡。水合物層的水合物初始飽和度為0.45，上覆層和下伏層全部為水，水飽和度為1。井底的壓力保持定壓力生產，儲層的外邊界為保持原始地層壓力的定壓邊界條件。

初始的地應力分布可由飽和土土體自重折算。模型頂面的水深為1 273.9 m，折算頂部壓力為12.86 MPa，沉積物密度為2 650 kg/m3，則地層應力以25.97 kPa/m的梯度縱向遞增。固體變形場的邊界條件為：上頂面為自由面邊界；儲層下底面為縱向固定條件；側面為水平位移固定條件，即垂直于x=0 m，x=800 m的側面，沿x方向的位移為0；垂直于y=0 m和y=800 m側面，沿y方向的位移為0。

2.6 物性參數

模型水深為1 273.9 m，上覆層厚度為135 m，下伏層厚度為94 m，井底生產壓力為5.0 MPa，水合物層底界初始壓力14.3 MPa，水合物層底界初始溫度為14.0 ℃，其熱物性、滲透率等參數如表1所示。

含水合物沉積物力學性質參數如表2所示。

3 模型驗證

為驗證數值模型及程序的正確性，將模型計算結果與Masuda實驗結果進行對比。Masuda采用Berea砂巖合成水合物進行降壓開采實驗的模型如圖3所示。

Berea巖心長30 cm，直徑5.1 cm，放置于溫度為2.3 ℃的恒溫空氣浴中。巖心的初始溫度（Ti）為2.3 ℃，初始孔隙壓力（pgi）為3.75 MPa，初始水合物飽和度（Shi）為0.443，水飽和度（Swi）為0.351，孔隙度（φi）為0.182，滲透率（Ki）為98 mD。實驗過程中，A端作為出口保持2.84 MPa的恒定壓力。實驗記錄出口的累計產氣量。

表1 南海神狐海域W19站位水合物儲層物性參數表

表2 南海神狐海域含水合物沉積物力學參數表

圖3 Berea巖心實驗模型示意圖

利用本文建立的有限元計算程序計算上述算例，獲得出口端的累計產氣量隨時間的變化。與實驗結果進行對比。累計產氣量對比結果如圖4所示。從對比結果看，筆者建立的多場耦合模型可以與實驗結果較好地吻合，證明了模型和算法的有效性和可靠性。

圖4 本文計算的累計產氣量與實驗結果對比圖

4 結果分析及討論

4.1 產水產氣特征

圖5是井底壓力定為5 MPa條件下產水和產氣速率隨時間變化曲線。由圖5可知，井底壓力降低后，井附近的地層壓力隨之降低，這導致井附近水合物分解，試采井產水產氣速率開始保持一個較高值，之后迅速降低。待產水速率上升時，產氣速率逐漸增加；由于壓力傳導的速度較快，產水速率能很快到達相對穩定的狀態；產氣速率也到達相對穩定的波動狀態。

4.2 物理場分布特征

在本模型均質假設條件下，預測得到開采60 d時的壓力變化、水合物飽和度、溫度變化、氣體飽和度等在空間上的變化特征（圖6）。

由圖6-a可知，井壓力降低區域主要集中在試采井附近，生產井中心的壓力最低，壓力比初始地層壓力降低約9 MPa；其水平方向上的影響范圍大致為井附近35 m內，即壓降從井中心的9 MPa到0 MPa的范圍為35 m；在5 MPa生產壓力下，水合物飽和度分解區被限制在生產井附近，在水平方向上，水合物的分解區大約離井2 m范圍內。因儲層滲透率較低，儲層底部的水合物還未完全分解，起到了一定的阻水作用。

由圖6-c可知，開采60 d時，水合物分解的吸熱效應并不能造成溫度在空間上明顯的變化，井中心的溫度最大降低約－4 ℃；沿水平方向，溫度的降低范圍較小，這印證了水合物分解范圍較小的事實。水合物分解產生的甲烷氣體一部分運移到生產井產出，一部分積聚在孔隙空間中，由于氣體飽和度的增加，在兩相滲流的毛細管力作用下，氣體不能完全產出，形成了如圖6-d的氣體飽和度分布。

4.3 力學響應特征分析

在儲層選取兩個點，監測其孔隙壓力、溫度和應力隨時間的變化情況。近井處和遠井處的坐標分別為（x=0.3 m，z=149 m）和（x=8.1 m，z=149 m）。

圖5 GMGS3-W19站位產氣產水隨時間變化曲線圖

圖6 試采60 d后井周物理場變化圖（負值表示比初始值低）

圖7 z為149 m處的近井處和遠井處的力學響應變化曲線圖

從圖7-a可知，由于井底降壓，近井處的孔隙壓力迅速降低后達到穩定值。遠井處孔隙壓力表現為平緩下降趨勢，近井處孔隙壓力低于遠井處。近井處的孔隙壓力降低導致水合物分解，水合物分解吸熱導致其溫度降低（圖7-b）；當水合物分解完畢后，由于周圍傳熱原因，溫度逐漸回升。對于遠井處，其孔隙壓力降低不足以使得水合物大量分解，故其溫度基本保持不變。

由圖7-c可知，孔隙壓力降低引起有效主應力增加。近井處的有效應力因該處孔隙壓力的快速降低而較快地升高，最后保持不變（由于是定井底壓力開采）。遠井處的孔隙壓力降低幅度較小，其有效主應力增加較為緩慢。降壓使得近井處的最大與最小主應力之差比遠井處的大，故近井處剪應力增加較遠井處的明顯。

圖7-d是兩個點的最大和最小有效主應力隨開采時間的變化情況，即有效應力路徑，由圖7-d可知，在t=0時刻，兩點的應力狀態處于沉積物的摩爾—庫倫抗剪強度線[30]的破壞線之外，即表示該處在沉積歷史中已經處于壓縮穩定的狀態。當開采后，其應力路徑均表現為近井處0～1 d內的變化較快，1 d后變化緩慢；遠井處的應力變化滯后于近井處。因兩點的應力路徑均表現為偏離摩爾—庫倫強度線，其沒有靠近或達到破壞線上，表示沒有發生破壞。故基于本模型的初步預測結果顯示，開采60 d內儲層不會發生破壞。

4.4 儲層沉降分析

圖8顯示了開采過程中垂向位移的情況，即生產降壓導致海底沉積物發生的沉降。由圖8-a、c知，生產井降壓造成空間上的沉降形如漏斗；從俯視角度（xOy平面）看，沉降以井位置處為中心呈圓形分布，井位置處的沉降最大。從切面圖（圖8-b、d）可知，在生產井段處的沉降較小，而在生產井段上下的垂向位移變化最大。故生產井段上方附近的沉降量最大，而生產井段下方在滲透壓的作用下發生局部隆起。生產井段上方的土體在上覆應力作用下，其沉降量大于生產井段下方的隆起量。由于生產井段上方的土體在應力作用下整體發生沉降，生產井降壓造成的沉降影響范圍大于孔隙壓力影響范圍。

圖8 開采井段（135～162 m）引起沉降空間分布狀態圖

圖8 中在30 d時，井位置處的最大沉降量為0.032 m（即32 mm）；而海底面沉降約為0.014 m，海底面沉降范圍（＞5 mm）的半徑約為166 m；隨開采進行，60 d時生產井位置的最大沉降量為0.035 m，海底面沉降約為0.018 m，海底面沉降范圍（＞5 mm）的半徑約為232 m。沉降量和沉降范圍均隨生產進行逐漸增大。

圖9-a是生產井位置海底面的沉降量（垂向位移）隨時間的變化情況。在前15 d，沉降約0.009 m（即9 mm）。隨后，因孔隙壓力在空間上逐漸趨于動態平衡狀態，海底面處的沉降速率降低。60 d后沉降逐漸變得緩慢，最終達到的0.018 m。可見前15 d的沉降占到總沉降量的1/2，故定壓開采條件下，其沉降主要發生在開采前期。

圖9-b是生產井位置處不同時刻的沉降量(垂向位移)情況，由圖9-b可知，由于降壓導致的應力變化主要集中在生產井周圍，故降壓點附近發生較大位移；又因降壓點下方底部是固定不動的，不發生垂向位移，降壓點底部的隆起部分將因邊界效應使得其隆起量為0。模型頂部是自由面，可以自由移動，海底面以下一定深度范圍內的沉積物是整體下沉。

4.5 滲透率對儲層沉降的影響

滲透率是影響氣水運移和壓力影響范圍的關鍵因素。圖10是不同滲透率下，生產井位置處海底面的沉降隨時間變化關系圖。

由圖10可以看出，滲透率較低時，壓降范圍較小，海底面沉降的下降速率基本保持不變；而當滲透率增加時，壓降范圍增加，沉降先以較高的速率發生，之后沉降速率降低，最后海底面以較低的沉降速率發生沉降。對于滲透率分別為1.0 mD、2.5 mD和5.0 mD，以60 d時的沉降量為基準，其沉降一半所需要的時間分別為24 d、15 d和9.5 d。隨滲透率增加，沉降速率加快，達到相同沉降量的時間提前。

4.6 井底壓力對儲層沉降的影響

井底壓力大小直接地影響地層中孔隙壓力分布范圍，引起骨架有效應力的變化，進而影響儲層的沉降。圖11是不同生產壓力下，生產井位置處海底面的沉降在生產過程隨時間的變化規律。從圖11可以看出，在生產前期，不同生產壓力下的沉降基本一致，差異較小；待進入穩定產氣速率階段后，海底面的沉降逐漸產生差異。在60 d時，其沉降量分別為0.016 m、0.018 m和0.020 m；其沉降一半所需時間約為15 d。生產壓力降低，沉降速率增加，達到相同沉降量所需時間提前，但其影響程度比滲透率的影響程度小。

圖9 生產井位置的海底面沉降量隨時間和儲層深度的變化圖

圖10 生產井位置的海底面沉降隨滲透率變化圖

圖11 生產井位置的海底面沉降隨生產壓力變化圖

5 結論

1）以南海北部神狐海域的天然氣水合物鉆探資料為基礎，建立了水合物單一垂直井降壓開采的物理模型，利用非結構網格對模型進行劃分；考慮水合物開采過程中的傳熱傳質和沉積物變形過程，建立了熱—流—固—化的四場耦合模型，基于非結構網格，采用有限元方法對模型進行求解，獲得儲層的孔隙壓力、溫度、飽和度和應力的時空演化特征。

2）神狐海域水合物儲層滲透率較低，降壓開采時壓力降低的影響范圍有限，主要局限在井筒附近范圍內，水合物的分解范圍也較小。

3）水合物開采過程中，儲層中孔隙壓力的減小導致了有效應力的增加，有效應力的增加主要在井筒附近，且井筒附近的最大和最小主應力之差比遠井處的大，故近井地帶的剪應力較大，易發生剪切破壞。開采60 d井筒附近應力路徑線遠離摩爾—庫倫的強度包線，這表明在本模型的假設條件下，儲層不會發生破壞。

4）有效應力增大導致儲層的沉降，開采60 d，儲層最大沉降為32 mm，海底面最大沉降為14 mm，且儲層沉降主要發生在開采的早期。

5）儲層滲透率和降壓開采的井底壓力對儲層沉降影響明顯。儲層滲透率越大、井底壓力降壓幅度越大，儲層沉降量越大，且沉降的速度越快。

致謝：感謝日本東京大學Shigemi Naganawa博士提供的Masuda實驗數據。

符 號 說 明

t表示時間，s；pw、pg分別表示水相、氣相壓力，Pa；pc表示毛細管壓力，Pa；pe表示水合物平衡壓力，Pa；Sh、Sg、Sw分別表示水合物、氣相、水相的飽和度；ρh、ρw、ρs、ρg分別表示水合物、水相、沉積物顆粒、氣相的密度，kg/m3；ρm表示水合物和沉積物顆粒組成復合固體的密度，kreac表示水合物分解速率常數，mol/(m2·Pa·s)；Ars表示單位體積儲層水合物分解的表面積，m－1；mw、mh、mg分別表示水合物分解產水速率、單位體積水合物分解速率、單位體積儲層中水合物分解的產氣速率，kg/(m3·s)；Nh表示水合物數；MCH4表示甲烷氣體的摩爾質量，kg/mol；MH2O表示水相的摩爾質量，kg/mol；Mh表示水合物摩爾質量，kg/mol；Krg、Krw分別表示氣水兩相滲流的氣相、水相的相對滲透率；K、K0分別表示沉積物多孔介質、不含水合物時的絕對滲透率，m2；Krgo、Krwo表示無水合物時的氣相端點、油相端點的相對滲透率；n表示滲透率遞減指數；ng、nw分別表示氣相、水相相對滲透率指數；Sgr、Swr分別表示氣相、水相相對滲透率端點的飽和度； μg、μw分別表示甲烷氣體、水相黏度，Pa·s；φ表示含水合物沉積物孔隙度；g表示重力加速度，m/s2；T表示儲層溫度，K；Cw、Cg、Ch、Cs分別表示水相、氣相、水合物相、沉積物顆粒相的比熱容，分別表示水相、氣相相對于控制體的速度，m/s；λs、λg、λw、λh分別表示沉積物顆粒、甲烷氣相、水相、水合物熱傳導系數，W/(m·K)；Qh表示水合物分解的相變潛熱，W/m3；σ表示含水合物沉積物的總應力張量，MPa；σ'表示復合固體的有效應力，MPa；σc0表示參考圍壓，MPa；α表示Biot系數；ε表示應變張量；I表示單位向量； 表示表示沉積物固體變形位移，m；Gm、lm分別表示水合物和沉積物顆粒組成復合固體的拉梅常數；Em、Es0、Eh分別表示復合固體、不含水合物沉積物在參考圍壓σc0下、純水合物的彈性模量，MPa；vm表示復合固體泊松比；b表示σc圍壓條件下，沉積物彈性模量的敏感系數；c表示水合物彈性模量系數；d表示水合物飽和度非線性效應系數；A0～A5，B1、B2表示常系數。

[ 1 ] Klauda JB & Sandler SI. Global distribution of methane hydrate in ocean sediment[J]. Energy & Fuels, 2016, 19(2): 459-470.

[ 2 ] Dubreuil-Boisclair C, Gloaguen E, Bellef l eur G & Marcotte D.Non-Gaussian gas hydrate grade simulation at the Mallik site,Mackenzie Delta, Canada[J]. Marine and Petroleum Geology,2012, 35(1): 20-27.

[ 3 ] Uddin M, Wright F, Dallimore S & Coombe D. Gas hydrate dissociations in Mallik hydrate bearing zones A, B, and C by depressurization: Eあect of salinity and hydration number in hydrate dissociation[J]. Journal of Natural Gas Science and Engineering,2014, 21: 40-63.

[ 4 ] Schoderbek D, Martin KL, Howard J, Silpngarmlert S & Hester K.North Slope hydrate fi eld trial: CO2/CH4exchange[C]//OTC Arctic Technology Conference, 3-5 December 2012, Houston, Texas,USA. DOI: http://dx.doi.org/10.4043/23725-MS.

[ 5 ] Yamamoto K, Terao Y, Fujii T, Terumichi I, Seki M, Matsuzawa M, et al. Operational overview of the fi rst oあshore production test of methane hydrates in the Eastern Nankai Trough[C]. Oあshore Technology Conference, 5-8 May, Texas, USA.

[ 6 ] 陳惠玲, 朱夏. 南海神狐海域天然氣水合物試采60天關井[N].中國國土資源報, 2017-07-10(1).Chen Huiling & Zhu Xia. The trial production well of gas hydrate from Shenhu area, South China Sea, shut down on the 60thday of production[N]. China Land Resource Report, 2017-07-10(1).

[ 7 ] Kurihara M, Sato A, Ouchi H, Narita H, Masuda Y, Saeki T, et al.Prediction of gas productivity from eastern Nankai trough methane-hydrate reservoirs[J]. SPE Reservoir Evaluation & Engineering, 2009, 12(3): 477-499.

[ 8 ] 梁金強, 王宏斌, 蘇新, 付少英, 王力峰, 郭依群, 等. 南海北部陸坡天然氣水合物成藏條件及其控制因素[J]. 天然氣工業,2014, 34(7): 128-135.Liang Jinqiang, Wang Hongbin, Su Xin, Fu Shaoying, Wang Lifeng, Guo Yiqun, et al. Natural gas hydrate formation conditions and the associated controlling factors in the northern slope of the South China Sea[J]. Natural Gas Industry, 2014, 34(7):128-135.

[ 9 ] 于興河, 王建忠, 梁金強, 李順利, 曾小明, 沙志彬, 等. 南海北部陸坡天然氣水合物沉積成藏特征[J]. 石油學報, 2014,35(2): 253-264.Yu Xinghe, Wang Jianzhong, Liang Jinqiang, Li Shunli, Zeng Xiaoming, Sha Zhibin, et al. Depositional accumulation characteristics of gas hydrate in the northern continental slope of South China Sea[J]. Acta Petrolei Sinica, 2014, 35(2): 253-264.

[10] Zhang KN, Moridis GJ, Wu NY, Li XS & Reagan MT. Evaluation of alternative horizontal well designs for gas production from hydrate deposits in the Shenhu area, South China Sea[C]//International Oil and Gas Conference & Exhibition, 8-10 June 2010,Beijing, China. DOI: http://dx.doi.org/10.2118/131151-MS.

[11] 胡立堂, 張可霓, 高童. 南海神狐海域天然氣水合物注熱降壓開采數值模擬研究[J]. 現代地質, 2011, 25(4): 675-681.Hu Litang, Zhang Keni & Gao Tong. Numerical studies of gas production from gas hydrate zone using heat injection and depressurization in Shenhu area, the South China Sea[J]. Geoscience,2011, 25(4): 675-681.

[12] 李剛, 李小森, Zhang Keni, Moridis GJ. 水平井開采南海神狐海域天然氣水合物數值模擬[J]. 地球物理學報, 2011, 54(9):2325-2337.Li Gang, Li Xiaosen, Zhang KN & Moridis GJ. Numerical simulation of gas production from hydrate accumulations using a single horizontal well in Shenhu area, South China Sea[J]. Chinese Journal of Geophysics, 2011, 54(9): 2325-2337.

[13] 李剛, 李小森, 陳琦, 陳朝陽. 南海神狐海域天然氣水合物開采數值模擬[J]. 化學學報, 2010, 68(11): 1083-1092.Li Gang, Li Xiaosen, Chen Qi & Chen Zhaoyang. Numerical simulation of gas production from gas hydrate zone in Shenhu area, South China Sea[J]. Acta Chimica Sinica, 2010, 68(11):1083-1092.

[14] Li G, Moridis GJ, Zhang KN & Li XS. Evaluation of gas production potential from marine gas hydrate deposits in Shenhu area of South China Sea[J]. Energy & Fuels, 2010, 24(11): 6018-6033.

[15] 蘇正, 何勇, 吳能友. 南海北部神狐海域天然氣水合物熱激發開采潛力的數值模擬分析[J]. 熱帶海洋學報, 2012, 31(5): 74-82.Su Zheng, He Yong & Wu Nengyou. Numerical simulation on production potential of hydrate deposits by thermal stimulation[J]. Journal of Tropical Oceanography, 2012, 31(5): 74-82.

[16] 沈海超, 程遠方, 胡曉慶. 天然氣水合物藏降壓開采近井儲層穩定性數值模擬[J]. 石油鉆探技術, 2012, 40(2): 76-81.Shen Haichao, Cheng Yuanfang & Hu Xiaoqing. Numerical simulation of near wellbore reservoir stability during gas hydrate production by depressurization[J]. Petroleum Drilling Techniques,2012, 40(2): 76-81.

[17] 程家望, 蘇正, 吳能友. 天然氣水合物降壓開采儲層穩定性模型分析[J]. 新能源進展, 2016, 4(1): 33-41.Cheng Jiawang, Su Zheng & Wu Nengyou. A geomechanical stability model analysis of hydrate reservoir for gas hydrate production by depressurization[J]. Advances in New and Renewable Energy, 2016, 4(1): 33-41.

[18] 孫可明, 王婷婷, 翟誠, 羅慧玉. 天然氣水合物加熱分解儲層變形破壞規律研究[J]. 特種油氣藏, 2017, 24(5): 91-96.Sun Keming, Wang Tingting, Zhai Cheng & Luo Huiyu. Patterns of reservoir deformation due to thermal decomposition of natural gas hydrates[J]. Special Oil & Gas Reservoirs, 2017, 24(5): 91-96.

[19] Sun JX, Zhang L, Ning FL, Lei HW, Liu TL, Hu GW, et al. Production potential and stability of hydrate-bearing sediments at the site GMGS3-W19 in the South China Sea: A preliminary feasibility study[J]. Marine and Petroleum Geology, 2017, 86: 447-473.

[20] Yang S, Zhang M, Liang J, Lu J, Zhang Z, Melanie H, et al. Preliminary results of China's third gas hydrate drilling expedition:A critical step from discovery to development in the South China Sea[J]. Fire in the Ice, 2015, 15(2): 1-5.

[21] Kim HC, Bishnoi PR, Heidemann RA & Rizvi SSH. Kinetics of methane hydrate decomposition[J]. Chemical Engineering Science, 1985, 42(7): 1645-1653.

[22] Biot MA & Willis DG. The elastic coeきcients of the theory of consolidation[J]. Journal of Applied Mechanics, 1957, 15(2):594-601.

[23] Masuda Y, Naganawa S & Ando S. Numerical calculation of gas production performance from reservoirs containing natural gas hydrates[J]. SPE Journal, 1997, 29(3): 201-210.

[24] Sun XF & Mohanty KK. Kinetic simulation of methane hydrate formation and dissociation in porous media[J]. Chemical Engineering Science, 2006, 61(11): 3476-3495.

[25] Duan ZH & Sun R. A model to predict phase equilibrium of CH4and CO2clathrate hydrate in aqueous electrolyte solutions[J].American Mineralogist, 2006, 91(8/9): 1346-1354.

[26] Esmaeilzadeh F, Zeighami ME & Fathi J. 1-D modeling of hydrate decomposition in porous media[C]//Proceedings of World Academy of Science: Engineering & Technology, 2013, 43: 648.

[27] Soga K, Lee SL, Ng MYA & Klar A. Characterisation and engineering properties of methane hydrate soils[C]//Proceedings of the 2ndInternational Workshop on Characterisation and Engineering Properties of Natural Soils, 29 November-1 December 2006,Singapore. London: CRC Press, 2006.

[28] Santamarina JC & Ruppel C. The impact of hydrate saturation on the mechanical, electrical, and thermal properties of hydrate-bearing sand, silts, and clay[C]//Proceedings of the 6thInternational Conference on Gas Hydrate, 6-10 July 2008, Vancouver, British Columbia, Canada.

[29] Clarke M & Bishnoi PR. Determination of the intrinsic rate of gas hydrate decomposition using particle size analysis[J]. Annals of the New York Academy of Sciences, 2000, 912: 556-563.

[30] Rutqvist J, Moridis GJ, Grover T, Silpngarmlert S, Collett TS& Holdich SA. Coupled multiphase fluid flow and wellbore stability analysis associated with gas production from oceanic hydrate-bearing sediments[J]. Journal of Petroleum Science &Engineering, 2012(92/93): 65-81.