裂隙充填型天然氣水合物的各向異性分析及飽和度估算——以印度東海岸NGHP01－10D井為例

王吉亮，王秀娟，錢 進，吳時國

1 中國科學院海洋研究所 海洋地質與環(huán)境重點實驗室，青島 266071

2 中國科學院大學，北京 100049

3 國土資源部海洋油氣資源與環(huán)境地質重點實驗室，青島 266071

1 引 言

天然氣水合物是由水分子與氣體分子組成的一種似冰狀固態(tài)化合物，廣泛分布在世界深水盆地和凍土帶［1］．最近水合物鉆探在美國的墨西哥灣盆地［2］、韓國的郁龍盆地［3］、印度的克里希納－戈達瓦里盆地［4］的細粒沉積物中發(fā)現(xiàn)了裂隙充填型水合物．壓力取芯表明細粒沉積物中天然氣水合物具有兩種明顯不同的產出形態(tài)，一種是孔隙充填，另一種是顆粒驅替．充填在孔隙空間的水合物替代了沉積物顆?？紫堕g的流體，但生成在裂隙內的水合物并不占據孔隙體積，而是迫使沉積層張開形成裂隙產出層狀、脈狀和結核狀的水合物［4－7］，即天然氣水合物占據原來的顆?？臻g，發(fā)生顆粒驅替．裂隙充填型的水合物與孔隙充填型水合物不同，含水合物層的裂隙傾角一般較陡，由于受構造應力作用，裂隙分布與構造主應力有關，一般定向排列，類似于裂隙介質而具有各向異性．裂隙充填型天然氣水合物是僅次于極地和海洋砂巖儲層的分布類型［8］，印度被動大陸邊緣是裂隙充填型水合物較為典型的一個區(qū)域．2006年印度國家天然氣水合物計劃（NGHP01）在印度克里希納－戈達瓦里（Krishna－Godavari，簡稱K－G）盆地NGHP01－10 站位的泥質沉積物中鉆探到高飽和度水合物，壓力取芯和X 射線成像顯示水合物呈結核狀或脈狀充填在高角度的裂隙中［9］．

不同學者提出了多種理論和半經驗模型，利用聲波速度、電阻率測井資料來估算均勻各向同性的孔隙充填型天然氣水合物的飽和度．但是，利用隨鉆電阻率測井基于阿爾奇公式計算的NGHP01－10A井水合物飽和度占孔隙空間80%，而利用壓力取芯計算的水合物飽和度占孔隙空間20%左右［10］，兩種不同方法計算結果相差較大．細粒沉積物中水合物呈脈狀充填在裂隙中，裂隙傾角與區(qū)域構造應力有關，這種定向裂隙導致含水合物層出現(xiàn)各向異性［9］．利用孔隙充填型假設水合物呈均勻分布的各向同性速度模型計算水合物飽和度時產生較大誤差，與壓力取芯計算結果差異較大［5］．關于定向排列裂隙導致的各向異性，前人進行了大量研究，存在多種模型，如層狀介質模型［11］、裂隙嵌于孔隙介質中模型［7，12－15］、周期性薄互層與擴容模型［16］等．不同模型具有不同假設條件，基于這些模型都可以利用速度來研究裂隙充填型水合物的飽和度．但在實際應用中，裂隙方向是影響估算水合物飽和度精度的一個關鍵因素．假設裂隙中完全充填水合物，裂隙充填型水合物儲層可以利用兩端元的層狀介質模型來研究，一為裂隙端元，充填水合物；另一端元是各向同性的飽和水沉積物．Lee等［5，17］利用該模型基于縱波速度分別計算了墨西哥灣和印度NGHP01－10D井裂隙充填型水合物的飽和度．Lee等［5，17］的方法雖然考慮了裂隙引起的各向異性，但在應用中需要假定裂隙沿固定方向，水平裂隙或者垂直裂隙，沒有考慮裂隙傾角變化，計算過程僅用縱波速度或者橫波速度中的一個，導致計算結果與實測壓力取芯結果相比存在一定誤差．本文中我們利用層狀介質模型，假設裂隙中完全充填水合物，在垂直井孔中，波入射角與裂隙傾角有關，考慮裂隙傾角變化，利用縱波和橫波速度同時對水合物飽和度和裂隙傾角進行反演．

2 理論模型

2．1 層狀介質的速度模型

超壓流體和氣體使泥質沉積物形成裂隙［18］．圖1給出了裂隙充填型天然氣水合物模型，假設裂隙內完全充填水合物，利用兩個端元的層狀介質模型來模擬裂隙內的水合物和飽和水孔隙介質中的沉積物．模型由I和II兩端元組成，I為裂隙，100%充填水合物，II為孔隙中飽和水的各向同性介質，η1 為裂隙所占的體積分數，假設裂隙中完全充填水合物，φ1＝100%，η2 為各向同性介質所占的體積分數，φ2為飽和水的孔隙度．

圖1 含天然氣水合物層X 射線成像和兩個端元的層狀介質模型［5］Fig．1 X－ray image of gas hydrate bearing sediments and two－component laminated media model［5］

端元Ⅱ中，各彈性參數利用簡化的三相介質模型來計算［5］，利用Voigt－Reuss－Hill［19］平均來計算礦物顆粒間彈性參數，然后用Pride［20］公式來計算干燥沉積物骨架的彈性模量：

式中α稱為膠結因子，是與沉積物固結程度有關的參數，膠結因子與有效壓力和固結程度有關．Ks和μs分別為礦物顆粒的體積模量和剪切模量，Kd和μd分別為干巖石骨架的體積模量和剪切模量，φ為干巖石骨架的孔隙度．Mindlin［21］（1949）認為剪切模量與有效壓力1/3次冪呈正比．Lee［5］假設α與埋藏深度的1/3次方成正比來計算膠結因子．飽和水各向同性介質沉積物的彈性參量可以利用Gassmann［22］流體替換．

含裂隙時橫向各向同性介質的相速度表示為：

式中G是圖1中I和II中任意彈性參數或者彈性參數的組合．橫向各向同性介質中P 波和S波速度用拉梅常數λ和μ計算［23］：

式中φ為入射波與層狀介質對稱軸之間的夾角，λ和μ為拉梅常數．Vp為準縱波速度（簡稱縱波），為水平極化橫波，為垂直極化橫波，也是一種準橫波（簡稱橫波）．由公式（3）計算得到的速度是各向異性介質中的相速度，Thomson［15］給出相速度與群速度之間的關系．利用速度反演水合物飽和度時，通常利用群速度．對于垂向井孔，入射角為0°表示水平裂隙，入射角90°表示垂向裂隙［5］，即入射角與裂隙傾角相關，入射角的大小代表了裂隙傾角的大?。?/p>

Thomson［15］利用γ，δ和ε三個參數來描述弱各向異性介質，這些參數與White模型中各參數間的關系為：

同時，Thomson［15］給出了橫向各向同性介質群速度與相速度之間的關系式為：

式中φg 表示射線方向與層狀介質對稱軸之間的夾角，GVp、GVhs和GVvs分別是縱波、水平極化橫波和垂向極化橫波的群速度．公式（5）中群速度與相速度并不相等，在φg 對應的波前法線角為φ時，可以利用相速度來計算群速度．對于弱各向異性介質，φg和φ滿足以下關系：

式中，α0和β0 分別為φ＝0°時的縱、橫波速度．

2．2 模型分析

在裂隙充填型天然氣水合物層，測井測量的速度與裂隙傾角和水合物飽和度有關，在相同水合物含量時，不同的裂隙傾角對應不同的縱橫波速度．圖2中給出了當裂隙水平和垂直時，縱波和水平極化橫波（Vhs）速度隨水合物體積分數的變化．假設端元II中孔隙度為φ2＝0．6，沉積物礦物組分及百分含量為NGHP01－10D 井的巖芯分析結果，見表1．當水合物體積分數Vh＝0時，表明模型中不存在裂隙充填型水合物，計算的速度為飽和水的孔隙介質的速度；當Vh＝1時，表明模型中全部為裂隙充填型水合物，計算的速度為純水合物的速度．圖2表明縱波和水平極化橫波（Vhs）速度都隨著水合物體積分數增加而增大．在一定的水合物體積分數下，裂隙傾角越大，各向異性介質模型中的縱波橫波速度越大．

圖2 不同裂隙傾角下縱波和水平極化橫波速隨天然氣水合物體積分數的變化Fig．2 Variations between P－wave and horizontally polarized S－wave velocities versus volume fraction of gas hydrate in different fracture dip angles

表1 NGHP01－10D井彈性參數及礦物組分含量［24］Table 1 The elastic parameters and mineralogy content［24］of NGHP01－10D

圖3為縱橫波速度隨入射角變化．模型中孔隙度為φ2＝0．6，水合物體積分數為0．1．在入射角較小時，縱波速度變化不大；隨入射角增大，縱波速度逐漸增加．同樣，水平極化橫波（Vhs）隨著入射角增加而增加．

圖3 水合物體積分數為0．1時，縱波和水平極化橫波速隨入射角變化Fig．3 Variation between P－wave and horizontally polarized S－wave velocities versus incidence angle at volume fraction of gas hydrate equals to 0．1

3 裂隙充填型天然氣水合物飽和度估算

3．1 含天然氣水合物層的特征分析

2006年印度在克里希納－戈達瓦里（K－G）盆地（圖4a和4b）的NGHP01－10、12、13、21等多個站位發(fā)現(xiàn)裂隙充填型天然氣水合物，其中NGHP01－10站位水深1038 m．NGHP01－10A 井是鉆探到的水合物飽和度最高井之一，裂隙充填型水合物厚度達135m，裂隙傾角為41°～88°，縱波速度明顯增高，局部達2000m/s，電阻率達幾十至上百歐米（圖4c和5）．NGHP01－10D 井位于相鄰25m 處，該井位進行了縱橫波速度測量．NGHP01－13站位位于相鄰500m處，電阻率測井出現(xiàn)異常高值，沒有進行聲波測井．NGHP01－12站位位于250 m 處，壓力取芯顯示裂隙充填型水合物在深度70m 處，飽和度達30%．過該井的地震剖面顯示裂隙充填型水合物的BSR并不明顯，含水合物層呈現(xiàn)振幅空白現(xiàn)象，不同站位裂隙充填型水合物厚度和位置略微不同（圖4c）．NGHP01－10站位附近出現(xiàn)的裂隙充填型水合物恰好位于斷層交匯處，并且下部聚集了大量的游離氣．

圖4 （a）印度K－G 盆地地理位置；（b）K－G 盆地水深分布以及站位NGHP01－10，13，12和21位置；（c）三維測線中過NGHP01－10A，13A，12A 和21B井的測線及鉆井位置，綠線為10A 和13A 井電阻率測井；裂隙充填型水合物位于不同深度和不同厚度；地震剖面上BSR 不明顯．Fig．4 （a）Location of K－G Basin；（b）K－G Basin water depth and Site NGHP01－10，13，12，21distribution；（c）The seismic profile passed sites NGHP01－10A，13A，12Aand 21Band resistivity at hole 10Aand 13A．Gas hydrate filled in fractures showing different thickness and different depth．BSR was absent．

NGHP01－10站位有A、B、C 和D 四個孔，其中NGHP01－10A 進行隨鉆測井，但沒有測量橫波速度；在B、C井孔進行了壓力取芯；D 孔先做壓力取芯，然后進行了縱波、橫波、密度等電纜測井（圖5）．由于NGHP01－10D井僅測量了深部含水合物層數據，無法確定飽和水的背景速度，而A 孔測井資料從深度20m 到180m．在深度30～150m 出現(xiàn)縱波速度增加，大于利用簡化的三相介質理論計算的飽和水沉積層速度，表明地層含有水合物．

圖5 電纜測井測量的NGHP01－10D 井P波和S波（藍線）速度和隨鉆測井測量的NGHP01－10A井P波速度（黑線）；利用簡化的三相介質模型［5］計算的飽和水背景速度（紅線）．Fig．5 P and S wave velocities （blue lines）from wireline logging at Hole NGHP－10D and P wave velocity from logging while drilling at Hole NGHP01－10A（black line）．The water－saturated baseline velocity calculated using simplified three－phase model［5］（red line）．

3．2 天然氣水合物飽和度估算

3．2．1 各向異性速度模型

基于2．1節(jié)的速度模型，當裂隙為垂直和水平時，利用縱波或橫波速度都可以計算水合物飽和度．由于水合物密度略微小于海水的密度，利用密度計算孔隙度近似等于地層總孔隙度［25］，得到的孔隙度為飽和水孔隙度與裂隙充填水合物孔隙度的綜合體現(xiàn)．當裂隙內充填水合物體積分數（Vh）增加時，飽和水孔隙度（φ2）降低，即

φt由密度測井計算得到：

式中，ρg 為顆粒骨架密度，取2．75g/cm3，ρw 為海水密度，取1．03g/cm3，ρb 為沉積物密度，由密度測井獲得．

裂隙充填型水合物飽和度（Sh）為

在垂直井孔中，當同時測量了縱波和橫波速度，利用2．1節(jié)中縱波和橫波速度模型能夠同時反演水合物飽和度和裂隙傾角．圖6 給出了本文利用NGHP01－10D 井縱波和橫波速度估算水合物飽和度和裂隙傾角的反演流程圖．在垂直井孔中，僅能測量到橫向極化橫波速度Vhs．利用公式（3）—（6）和表1中NGHP01－10D 井巖性資料，基于縱波速度、橫波速度和縱橫波速度聯(lián)合計算了水合物飽和度（圖7）．假設裂隙傾角為水平時，利用縱波速度計算的水合物飽和度為25%～40%之間，平均飽和度為34%左右；而利用橫波速度計算的水合物飽和度為25%至75%之間，平均飽和度為60%左右．在水平裂隙時，橫波速度計算的飽和度遠大于縱波速度計算結果．假設為垂直裂隙時，縱波速度計算飽和度為10%至25%之間，平均飽和度為20%左右，橫波計算飽和度為5%～15%之間，平均為8%左右．在垂直裂隙時，橫波速度計算飽和度小于縱波速度計算結果．利用縱波和橫波速度聯(lián)合計算時，計算的水合物飽和度為10%至25%之間，平均飽和度為24%左右，裂隙傾角變化范圍為75°～85°（圖8）．圖7給出NGHP01－10B和10D 不同深度的壓力取芯計算的水合物飽和度．從該圖可以看出，假設為垂直裂隙，利用縱波速度計算飽和度與縱橫波速度聯(lián)合計算的水合物飽和度與壓力取芯計算結果相吻合，在該井位，縱橫波速度聯(lián)合計算的飽和度吻合更好些．在10D 井沒有測量裂隙傾角，10A 井測量裂隙傾角主要分布范圍為60°～88°．

圖6 縱波和橫波速度聯(lián)合計算水合物飽和度流程圖Fig．6 Flow charts of calculating gas hydrate saturation using both P－wave and S－wave velocities

圖8 NGHP01－10A 井測井得到的裂隙傾角和利用縱波和橫波聯(lián)合反演得到的NGHP01－10D井裂隙傾角Fig．8 Fracture dips at Hole NGHP01－10A from logging and fracture dips calculated fromVp and Vs at Hole NGHP01－10D

3．2．2 各向同性有效介質模型（EMT）

有效介質模型（Effective Medium Theory）假設天然氣水合物是固體骨架的一部分，水合物生成降低了地層孔隙度、膠結干燥沉積物骨架，影響骨架的體積模量和剪切模量［23］．不含水合物干燥沉積物骨架的體積模量和剪切模量利用修改的Hashin－Shtrikman－Hertz－Mindlin理論［26－28］計算．骨架的礦物組分及含量見表1，臨界孔隙度為0．36．模型的理論基礎是水合物充填在孔隙空間，具有各向同性．圖7給出了利用縱波速度基于EMT 模型計算的NGHP01－10D水合物飽和度（綠色），計算的水合物飽和度為25%～60%之間，計算結果略微大于假設為水平裂隙時，縱波速度計算的水合物飽和度，遠大于垂直裂隙和壓力取芯計算的水合物飽和度．

圖7 NGHP01－10D井基于水平和垂直裂隙傾角，利用縱波和橫波速度估算的水合物飽和度及利用有效介質模型計算的水合物飽和度與壓力取芯計算的飽和度對比Fig．7 Gas hydrate saturations estimated from P－wave and S－wave velocities in horizontal and vertical fracture dip angles and gas hydrate saturation calculated from effective medium theory at Hole NGHP01－10Dand pressure cores

4 結 論

裂隙充填型天然氣水合物層具有明顯的各向異性，含水合物層縱波速度、橫波速度與裂隙傾角和水合物含量有關．在估算水合物飽和度時，裂隙對裂隙水合物層的橫波速度影響大于對縱波速度的影響．在低水合物含量和低裂隙傾角時，裂隙充填型水合物層的縱波速度變化不大，而孔隙充填型水合物層的縱波速度隨水合物飽和度增加而增加．在一定的裂隙充填型水合物含量時，隨著裂隙傾角的增大，縱波速度和水平極化橫波波速增大．

基于不同速度，利用層狀介質模型計算NGHP01－10D井天然氣水合物飽和度不同．利用水合物呈均勻分布的有效介質模型，計算的水合物飽和度遠大于壓力取芯計算的水合物飽和度．當水合物充填在裂隙中時，假設為水平裂隙時，利用縱波速度估算的水合物飽和度與利用有效介質模型假設水合物充填在孔隙空間估算的飽和度相差不大．假設為垂直裂隙時，利用縱波速度計算的水合物飽和度與壓力取芯相接近．橫波對裂隙傾角變化非常敏感，所以假設傾角水平或垂直時，利用橫波估算得到的水合物飽和度結果偏離壓力取芯結果較大．但是利用縱波和橫波速度聯(lián)合計算的水合物飽和度與壓力取芯結果吻合相對較好，水合物飽和度為10%～25%，平均為24%，裂隙的傾角在60°～90°，裂隙傾角較陡．裂隙傾角與區(qū)域構造環(huán)境有關而呈現(xiàn)定向特性，表明裂隙充填型水合物與孔隙充填型水合物不同，含水合物層具有各向異性．

致 謝 感謝印度國家水合物計劃01（NGHP01）航次的所有工作人員，文章中的測井數據和壓力取芯資料均是由他們采集獲得．感謝印度國家地球物理研究所對地震資料進行處理．

（References）

［1］ Kvenvolden K A．Gas hydrates－geological perspective and global change．ReviewofGeophysics，1993，31（2）：173－187．

［2］ Cook A，Guerin G，Mrozewski S，et al．Gulf of Mexico Gas Hydrate Joint Industry Project Leg II：Walker Ridge 313 LWD Operations and Results，2010：1－24．

［3］ Park K P，Bahk J J，Kwon Y，et al．Korean National Program Expedition Confirm Rich Gas Hydrate Deposits in the Ulleung Basin，East Sea：Fire in the Ice Methane Hydrate Newsletter，2006：6－9．

［4］ Collett T S，Riedel M，Cochran J，et al．Geologic Controls on the Occurrence of Gas Hydrates in the Indian Continental Margin：Results of the Indian National Gas Hydrate Program（NGHP）Expediton 01Initial Reports，2008：1－135．

［5］ Lee M W，Collett T S．Gas hydrate saturations estimated from fractured reservoir at Site NGHP－01－10，Krishna－Godavari Basin，India．J．Geophys．Res．，2009，114：B07102，doi：10．1029/2008JB006237．

［6］ Wojtowitz G，Zevos A，Clayton C R I．Numerical modeling of grain－displacing gas hydrate morphologies．//Proceedings of the 7th International Conference on Gas Hydrates（ICGH 2011）．Edinburgh，Scotland，United Kingdom，2011：8．

［7］ Ghosh R，Sain K，Ojha M．Effective medium modeling of gas hydrate－filled fractures using the sonic log in the Krishna－Godavari basin，offshore eastern India．J．Geophys．Res．，2011，115：B06101．

［8］ Boswell R，Collett T S．The gas hydrates resource pyramid：Fire in the ice．MethaneHydrateNewsletter，2006，6（3）：5－7．

［9］ Cook A E，Goldberg D．Extent of gas hydrate filled fracture planes：Implications for in situ methanogenesis and resource potential．Geophys．Res．Lett．，2008，35（15）：L15302．

［10］ Collett T S，Riedel M，Cochran J，et al．Geologic Controls on the Occurrence of Gas Hydrates in the Indian Continental Margin：Results of the Indian National Gas Hydrate Program（NGHP）Expediton 01Initial Reports．San Antonio：AAPG Annual Convention，2008．

［11］ White J E，Angona F A．Elastic wave velocities in laminated media．JournalofAcousticSocietyofAmerica，1955，27（2）：310－317．

［12］ Eshelby J D．The determination of the elastic field of an ellipsoidal inclusion，and related problems．Proc．R．Soc．A，1958，241（1226）：376－396．

［13］ Walsh J B．New analysis of attenuation in partially melted rock．J．Geophys．Res．，1969，74（17）：4333－4337．

［14］ Hudson J A．Wave speeds and attenuation of elastic waves in material containing cracks．Geophys．J．Int．，1981，64（1）：133－150．

［15］ Thomsen L．Weak elastic anisotropy．Geophysics，1986，51（10）：1954－1966．

［16］ Yang D H，Zhang Z J．Poroelastic wave equation including the Biot/squirt mechanism and the solid/fluid coupling anisotropy．WaveMotion，2002，35（3）：223－245．

［17］ Lee M W，Collett T S．Pore－and fracture－filling gas hydrate reservoirs in the Gulf of Mexico Gas Hydrate Joint Industry Leg II Green Canyon 955 H well．Marineand Petroleum Geology，2012，34（1）：62－71．

［18］ Sassen R，Losh L S，Cathles L，et al．Massive vein－filling gas hydrate：Relation to ongoing gas migration from the deep subsurface in the Gulf of Mexico．MarineandPetroleum Geology，2001，18（5）：551－560．

［19］ Hill R．The elastic behaviour of a crystalline aggregate．Proc．Phys．Soc．A，1952，65（5）：349－354．

［20］ Pride S R，Berryman J G，Harris J M．Seismic attenuation due to wave－induced flow．J．Geophys．Res．，2004，109：B01201．

［21］ Mindlin R D．Compliance of elastic bodies in contact．Journalof AppliedMechanicsTransactions，1949，71：259－268．

［22］ Gassmann F．Elasticity of porous media．VierteljahresheftderNaturforschendenGesselschaft，1951，96：1－23．

［23］ White J E．Seismic Waves－Radiation，Transmission，and Attenuation．New York：McGraw－Hill，1965．

［24］ Helgerud M B，Dvorkin J，Nur A，et al．Elastic－wave velocity in marine sediments with gas hydrates：Effective medium modeling．Geophys．Res．Lett．，1999，26（13）：2021－2024．

［25］ Lee M W．A simple method of predicting S－wave velocity．Geophysics，2006，71（6）：F161－F164．

［26］ Dvorkin J，Nur A．Elasticity of high－porosity sandstones：Theory for two North Sea data sets．Geophysics，1996，61（5）：1363－1370．

［27］ Dvorkin J， Nur A． Time－average equation revisited．Geophysics，1998，63（2）：460－464．

［28］ 王秀娟，吳時國，劉學偉．天然氣水合物和游離氣飽和度估算的影響因素．地球物理學報，2006，49（2）：504－511．

Wang X J，Wu S G，Liu X W．Factors affecting the estimation of gas hydrate and free gas saturation．ChineseJ．Geophys．（in Chinese），2006，49（2）：504－511．