中国海洋湖沼学会主办。
文章信息
- 吴景鑫, 郭秀军, 孙翔, 李宁. 2018.
- WU Jing-Xin, GUO Xiu-Jun, SUN Xiang, LI Ning. 2018.
- 神狐海域天然气水合物层分解井中电学监测效果模拟与分析
- NUMERICAL EVALUATION OF MULTI-ELECTRODE RESISTIVITY LOGGING MONITORING OF GAS-HYDRATE DECOMPOSITION IN SUBMARINE RESERVOIR IN SHENHU AREA, SOUTH CHINA SEA
- 海洋与湖沼, 49(6): 1211-1219
- Oceanologia et Limnologia Sinica, 49(6): 1211-1219.
- http://dx.doi.org/10.11693/hyhz20171100282
-
文章历史
- 收稿日期:2017-11-06
- 收修改稿日期:2018-03-16
2. 山东省海洋环境地质工程重点实验室 青岛 266100;
3. 海洋环境与生态教育部重点实验室 青岛 266100
2. Shandong Provincial Key Laboratory of Marine Environment and Geological Engineering, Ocean University China, Qingdao 266100, China;
3. Key Laboratory of Marine Environment and Ecology, Ministry of Education, Qingdao 266100, China
天然气水合物分解可能导致的海床沉降、麻坑、泥火山、滑坡等地质灾害是天然气水合物开采所面临的主要环境风险之一(Paull et al, 2001)。为保证开采工程安全, 开采过程中需要利用水压测量(Stenvold et al, 2006; Yokoyama et al, 2015)、时瞬地震测量(Hall et al, 2002)、三分量加速度测量(Higley et al, 2005)等技术监测海床沉降; 利用集成式孔压、温度探针测量技术监测表层沉积物物理力学性质(Higley et al, 2005), 利用生产井或监测井中不同深度位置温度、压力测量评价水合物层分解状态(Chee et al, 2014)。除此以外, 对开采过程水合物分解区空间变化实时监测同样重要。在日本南海海槽天然气水合物试采工程中, 采用了深海地震系统(Deep-sea Seismic System, DSS)来实现这一目标。这是一种以生产井为中心布设三维地震接收系统, 通过监测不同时期沉积层地震剖面信息反演水合物层分解区的方法(Takahashi et al, 2011)。
近年来, 关于沉积层电阻率与水合物饱和度关联关系的研究已有诸多进展(王秀娟等, 2010), 利用电测井数据反演水合物饱和度进而推断水合物层分布位置及储量已成为水合物探测的重要手段(Riedel et al, 2006; Wang et al, 2011)。水合物开采过程沉积体系电阻率变化的室内实验也揭示, 水合物层分解会导致沉积体系电阻率整体下降, 同时受局部位置温度变化或流体流动的干扰, 可能产生电阻率异常波动(李小森等, 2013; 李淑霞等, 2010; Priegnitz et al, 2015; Spangenberg et al, 2015)。这些相关研究为利用电阻率法监测水合物层分解区变化提供了物性前提。
基于此, 本文利用监测井设计一种井中电学监测系统来实现水合物层分解状态监测。为检验这种方法的探测能力, 以我国南海神狐海域天然气水合物试采区为对象, 通过归纳总结水合物开采过程水合物层分解和电阻率变化特征, 建立相应地质和电阻率模型, 利用电阻率法正、反演软件计算得到水合物层不同分解状态的电阻率监测剖面图像, 分析剖面图像异常特征, 界定异常图像和水合物层分解区的对应关系。
1 地质模型构建与开采工艺设定 1.1 研究区选取2007年, 我国在南海神狐海域开展了第一次天然气水合物钻探(GMGS-1), 并在SH2、SH3、SH7站位取到水合物样品(Zhang et al, 2007; 吴时国等, 2009)。2015年, 在该海域实施了第三次天然气水合物钻探(GMGS-3), 19个站位的测井资料均显示了天然气水合物的分布, 其中在W11、W17、W18及W19站位实施了取芯钻探(杨胜雄等, 2016), 站位分布见图 1。其中W19站位水合物层厚度大, 饱和度高, 是南海水合物开采的最佳远景区之一(苏丕波等, 2016), 本文将该站位位置作为研究区。
GMGS-3-W19站位水深为1273m。原位水合物层可能存在下伏游离气层(靳佳澎等, 2017), 本文水合物层地质模型构建时不考虑下伏游离气层干扰, 假设电阻率测井所显示的高阻异常均为水合物层。水合物储层分布在海床下134—202m之间。由于该站位与GMGS-1-SH7站位沉积物岩性相近, 都以粉砂质黏土和黏土的互层沉积为主(尚久靖等, 2016), 根据孔隙水淡化估计值得到的水合物饱和度, 将水合物储层由上到下分为三层(苏丕波等, 2016; 郭依群等, 2017):水合物层1(gas hydrate bearing-sediments, 简称GHBS 1);水合物层2(GHBS 2)和水合物层3(GHBS 3)。本文建模所选水合物开采工艺、地质参数均参考Sun等(2015)对SH7站位水合物开采过程模拟选用的参数。各层具体参数如表 1所示。
编号 | 岩性 | 厚度(m) | 水合物饱和度 | 孔隙度 | 渗透率(mD) | 孔隙水盐度(‰) | 温度(℃) |
GHBS 1 | 砂质 黏土 |
31 | 0.60 | 0.41 | 75 | 30.5 | 13.00 |
GHBS 2 | 粉砂质黏土 | 10 | 0.07 | 0.38 | 20 | 30.5 | 13.50 |
GHBS 3 | 黏土 | 27 | 0.12 | 0.37 | 10 | 30.5 | 14.00 |
模型构建中还假设水合物层顶底界面为可发生流体运移和热量交换的定压定温边界, 上覆层和下伏层为透水岩层, 允许发生流体运移与热量交换, 并分别与相邻水合物层具有一致的渗透率、孔隙度等物理属性(金光荣等, 2015)。
1.3 开采工艺设定设定开采工艺为单一垂直井恒压3.0MPa降压法, 生产井过滤器(降压区)位于海床下140—160m, 长度为20m。水合物层分解过程假设为三个阶段:
开采初期, 水合物的分解主要集中于生产井过滤器附近并且分解界面逐渐向前移动, 流体流动的驱动力为压力梯度(Huang et al, 2015; Han et al, 2016), 几乎所有生产井附近的游离气体以及水合物分解释放的气体都直接流动到了生产井, 并没有向上覆层流动, 相应储层分解及游离气聚集模型见图 2a。
![]() |
图 2 水合物开采过程水合物层分解及游离气聚集模型 Fig. 2 Geological model of reservoir decomposition and free gas accumulation during gas production 注: ![]() |
开采中期, 压力梯度和上下盖层的热流交换使得水合物层存在上下两个解离界面, 由于压力梯度的差异, 水合物层上部区域解离界面运动较快, 而气体主要聚集于生产井井口下部区域, 饱和度随压力场的扩散逐渐降低(Sun et al, 2015), 相应储层分解及游离气聚集模型见图 2b。
开采后期, 随着开采的进行, 分解界面逐渐向远离生产井位置运动, 下解离界面运动的主导因素为下伏层的热流交换, 由于远界面压差的减弱以及下伏层的热流交换强于上覆层, 上下解离界面的运动速度发生逆转, 下解离界面的运动速度连续且快于上解离界面, 水合物层上部解离速度减慢, 水合物层下部发生大规模的解离, 上下界面发生融合。界面融合后, 水合物层下部区域的水合物解离范围远大于水合物层上部区域。开采后期, 所生产甲烷主要为溶解态而非气态, 既包括分解产生的含甲烷孔隙水, 同时也包括一部分开采前中期滞留于水合物层的含甲烷孔隙水(Sun et al, 2015; Jin et al, 2016), 甲烷气仅分布于解离界面前沿, 生产井下部的水合物层一旦分解, 其产生的气体会被下伏层的流体溶解或置换, 气体仅沿解离界面呈锲形分布, 水合物层下部无气体饱和度, 两个界面融合后, 融合位置气体饱和度增加, 形成次生水合物, 而次生水合物又导致了气体流动阻力增加发生气体滞留(Huang et al, 2015; Sun et al, 2015)。该阶段水合物储层分解及游离气聚集模型见图 2c。
2 电阻率模型构建 2.1 初始电阻率模型构建根据GMGS-3-W19电测井数据(图 3a)对地层进行初始电阻率赋值, 构建初始电阻率模型(如图 3b)。电测井数据中, 蓝色曲线代表实测电阻率值, 红色曲线代表地层中无水合物赋存时电测井理论上应该测到的电阻率值。其中, 海床面以下0—134m为上覆土层, 电阻率ρ赋值为1.0Ω·m; 134—165m为水合物层1, 电阻率ρ赋值为7.0Ω·m; 165—175m为水合物层2, 电阻率ρ赋值为1.7Ω·m; 175—202m为水合物层3, 电阻率ρ赋值为2.0Ω·m; 202—230m以下为下伏层, 电阻率ρ赋值为1.5Ω·m。
![]() |
图 3 W19电测井数据及沉积层电学模型 Fig. 3 W19 electric logging data and hydrate-bearing sediments resistivity model 注: Z:距海床深度; ρ:电阻率; X:距离 |
(1) 假设沉积层中无水合物存在
依据阿尔奇公式(Archie, 1942), 饱水沉积层电阻率ρ0可表示为:
![](PIC/hyyhz-49-6-1211-E1.jpg)
式中, ρw为孔液电阻率; φ为孔隙度; a和m分别是岩性系数和胶结指数, 参考GMGS-1-SH3资料交会分析结果(王秀娟等, 2013), 分别取0.8和2。
考虑温度的影响, 不同深度沉积层中孔液电阻率ρw为变化值, 利用Arp’s方程进行孔液电阻率计算(Arps, 1953):
![](PIC/hyyhz-49-6-1211-E2.jpg)
式(2)中, ρw1、ρw2分别为温度T1、T2时孔液电阻率。
依据神狐海域GMGS-1期调查资料, 取T1=18℃, ρw1=0.24Ω·m的经验参数作为基准值(Guerin et al, 1999)。根据地温梯度变化规律, 水合物层1、2、3的孔液温度分别取13℃、13.5℃和14℃, 将上述参数分别代入公式(2)可求取不同深度沉积层中孔液电阻率。并进一步利用公式(1)求出不含水合物条件下, 不同深度沉积层电阻率。
(2) 假设沉积层中有水合物存在且均匀分布, 孔隙中只存在孔隙水与水合物, 不存在游离气
此时水合物饱和度sH与式(1)计算的沉积层电阻率的关系可由公式(3)表示:
![](PIC/hyyhz-49-6-1211-E3.jpg)
其中, ρ为水合物储层电阻率, n为饱和度指数。
根据式(1)计算的原始地层电阻率ρ0、测井资料确定的水合物储层初始电阻率与初始水合物饱和度, 可求出水合物层1、2、3的饱和度指数n1、n2、n3分别为1.83、1.82、1.94。
根据式(3)可得到水合物储层1、2、3电阻率与含水合物饱和度的关系曲线, 见图 4。
![]() |
图 4 天然气水合物储层电阻率与水合物饱和度关系 Fig. 4 Relationship between reservoir resistivity and saturation |
对于不同开采时期的GHBS1、GHBS2、GHBS3水合物层, 其饱和度相对初始饱和度降低, 参照水合物开采数值模拟确定不同区域饱和度后(Sun et al, 2015), 电阻率取值根据式(3)代入各水合物层相应饱和度指数所得, 具体参数取值见表 2。
开采阶段 | 模型区域 | 水合物饱和度 | 气饱和度 | 电阻率(Ω·m) |
开采初期 | GHBS 1中部水合物分解区 | 0 | 0.1 | 1.6 |
GHBS 1解离界面区 | 0.4 | 0 | 3.0 | |
开采中期 | GHBS 1分解区上部 | 0 | 0 | 1.3 |
GHBS 1分解区下部 | 0 | 0.1 | 1.6 | |
GHBS 1未分解含气区 | 0.6 | 0.1 | 11.8 | |
GHBS 2未分解含气区 | 0.07 | 0.05 | 1.9 | |
GHBS 3分解区 | 0 | 0 | 1.5 | |
开采后期 | GHBS 1分解区 | 0 | 0 | 1.5 |
GHBS 2分解区 | 0 | 0 | 1.3 | |
GHBS 2未分解含气区 | 0.07 | 0.05 | 1.9 | |
GHBS 3分解区 | 0 | 0 | 1.5 | |
GHBS 3未分解含气区 | 0.12 | 0.05 | 2.2 |
当沉积物孔隙中存在游离气时, 因游离气与水合物均为非导电介质, 故沉积物电阻率ρ1的计算公式为
![](PIC/hyyhz-49-6-1211-E4.jpg)
其中, SG为含气饱和度。
参照水合物开采数值模拟确定不同阶段不同区域含气饱和度后(Sun et al, 2015), 根据式(4)计算相应电阻率, 见表 2。根据水合物层分解模型和不同区域电阻率取值构建分解过程水合物层电阻率模型, 见图 5。
![]() |
图 5 不同开采阶段水合物储层电阻率模型 Fig. 5 The electrical model of reservoir decomposition at different stages 注: a:开采初期; b:开采中期; c:开采后期; Z:距海床深度; X:距离 |
电阻率监测系统布设于生产井右侧20m远的监测井中, 监测系统主要由井筒电极阵列、采集站、电源模块、传输通讯线及上位机组成, 采集站和电源模块设置在水密舱内, 见图 6。采集站包括串行电阻率采集单元和主控计算机, 主控计算机通过通讯电缆与平台或者生产船上的上位机连接。串行电阻率采集单元包括与电极一一对应的状态开关, 状态开关由四个继电器组成, 分别用于将电极设置成A、B、M、N四个状态, 状态开关接受主控计算机的控制进行工作。井筒电极阵列需垂直铺设在监测井中, 且井筒电极阵列的电极与监测井的井壁紧密贴合, 电极阵列朝向生产井方向, 电极阵列由多个单节点电极导杆连接组成, 单节点电极导杆包括绝缘的导杆本体和套设在导杆本体外表面上的电极。
![]() |
图 6 水合物层分解井中电阻率监测系统示意图 Fig. 6 The electrical well for monitoring hydrate decomposition |
本文设定监测井中布设128个电极, 极距1m, 采集层数40进行模拟。全空间电场正演计算采用有限元法, 将电阻率监测井左侧沉积层划分为50×512矩形网格块, 单位电极距四剖分, ∆z=0.25m, 横向网格宽度∆x=0.5m, 见图 7。不同网格块依据前文所建电学模型进行电阻率赋值, 选取采集装置后经正演计算得到不同网格节点处电位值, 进而计算得到水合物层分解不同阶段沉积层视电阻率剖面。反演过程采用最小二乘法, 所有模型均选择5次反演迭代结果进行讨论, 反演结果中RMS误差小于1%。本文正、反演计算采用GeoTomo公司的RES2Dmod ver. 3.02软件系统及RES2Dinv ver. 4.05软件系统实现。
![]() |
图 7 水合物层模型及正演网格剖分 Fig. 7 Reservoir model and forward mesh generation |
图 8为利用井中电阻率监测系统, 采用不同装置对天然气水合物储层进行探测, 所得理论剖面图像。图 8显示虽然不同装置探测剖面图像特征不同, 但储层分界面和电性特征都有良好的反映, 储层1表现为明显的高阻条带。偶极装置(Dp-Dp)、二极装置(P-P)、温纳-施隆贝格装置(W-S)、温纳装置(W)探测图像中均存在一些因方法引起的伴随异常, 这些异常的存在会影响对图像的分析判断, 相对而言, 偶极装置对水合物层分界面的探测更清晰, 各水合物层电性特征也最为稳定; 同时, 偶极装置与二极装置数据盲区也相对较少。故在实际探测中采用偶极装置。
![]() |
图 8 不同采集装置探测水合物储层理论剖面图像 Fig. 8 Monitoring results of data collection with different electrode arrays 注:![]() |
利用偶极装置对图 5所示不同分解阶段水合物储层电阻率模型进行探测, 理论探测图像见图 9。
![]() |
图 9 不同开采阶段水合物储层偶极装置理论电剖面图像 Fig. 9 Resistivity sections of different stages of hydrate decomposition in reservoir 注: ![]() |
图 9显示在开采初期, 随着分解区的形成反映水合物储层1的高阻条带中出现了低阻异常区, 异常区尺寸和分解区范围相当; 开采中期, 水合物层1分布区中的低阻异常区范围扩大, 同时在水合物层3分布区出现了低阻异常区; 开采后期, 分解区表现为贯穿三个水合物层的低阻“ < ”型异常条带, 条带分布范围和分解区范围相当。不同阶段的气体汇聚区和次生水合物分布区在图中未表现出明显异常特征。
4 监测能力评价为压制探测电阻率剖面中装置效应和水合物储层背景的影响, 更清晰显示不同分解阶段水合物储层空间分布特征, 将不同分解阶段水合物层探测电阻率值ρd与相应位置的原始水合物层探测电阻率值ρb做比值, 得到相对电阻率。绘制不同分解阶段水合物层相对电阻率图像, 见图 10。
![]() |
图 10 不同开采阶段水合物储层相对电阻率图像 Fig. 10 Relative resistivity sections of different reservoir decomposition stages 注: ![]() |
对比图 10和图 5可以看出, 不同分解阶段水合物层分解区都表现为明显的相对电阻率低值带, 分解区边界对应相对电阻率值急剧增大位置。选取相对电阻率1的等值线作为分解区边界, 和实际分解区边界误差不超过0.5m。不同阶段的气体汇聚区和次生水合物分布区在图中同样未表现出明显异常特征。
5 结论(1) 天然气水合物储层电阻率随水合物饱和度增大呈指数增大。开采过程水合物储层电阻率变化复杂, 分解区电阻率减小, 气体汇聚区和周边围岩电阻率差异不大, 次生水合物生成区电阻率变大。
(2) 利用井中电阻率系统进行水合物层分解状态探测时, 偶极装置是首选探测装置。在偶极装置探测剖面上不同饱和度水合物层表现为不同的层状电阻率异常带, 异常带界面和水合物层界面对应; 分解区表现为低阻异常区, 异常区范围和分解范围具有对应性。气体汇聚区和次生水合物分布区在图中未表现出明显异常特征。
(3) 利用探测结果处理形成的相对电阻率剖面可较准确的划定分解区边界, 误差不超过0.5m。
王秀娟, 吴时国, 王吉亮, 等, 2013. 南海北部神狐海域天然气水合物分解的测井异常. 地球物理学报, 56(8): 2799–2807 |
王秀娟, 吴时国, 刘学伟, 等, 2010. 基于电阻率测井的天然气水合物饱和度估算及估算精度分析. 现代地质, 24(5): 993–999 DOI:10.3969/j.issn.1000-8527.2010.05.022 |
苏丕波, 尉建功, 2016. 神狐海域首次发现Ⅱ型天然气水合物. 中国地质调查成果快讯, 1(21): 23–26 |
李小森, 冯景春, 李刚, 等, 2013. 电阻率在天然气水合物三维生成及开采过程中的变化特性模拟实验. 天然气工业, 33(7): 18–23 |
李淑霞, 夏唏冉, 郝永卯, 等, 2010. 电阻率测试技术在沉积物-盐水-甲烷水合物体系中的应用. 实验力学, 25(1): 95–99 |
吴时国, 董冬冬, 杨胜雄, 等, 2009. 南海北部陆坡细粒沉积物天然气水合物系统的形成模式初探. 地球物理学报, 52(7): 1849–1857 DOI:10.3969/j.issn.0001-5733.2009.07.019 |
杨胜雄, 张明, 梁金强, 等, 2016. 神狐海域钻探发现大型天然气水合物矿体. 中国地质调查成果快讯, 3(21): 6–8 |
尚久靖, 郭依群, 龚跃华, 等, 2016. 神狐海域天然气水合物钻探区沉积体系研究. 中国地质调查成果快讯, 1(21): 19–22 |
金光荣, 许天福, 刘肖, 等, 2015. 天然气水合物降压热激法模拟开采方案优化研究. 中南大学学报(自然科学版), 46(4): 1534–1543 |
郭依群, 杨胜雄, 梁金强, 等, 2017. 南海北部神狐海域高饱和度天然气水合物分布特征. 地学前缘, 24(4): 24–31 |
靳佳澎, 王秀娟, 陈端新, 等, 2017. 基于测井与地震多属性分析神狐海域天然气水合物分布特征. 海洋地质与第四纪地质, 37(5): 122–130 |
Archie G E., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Transactions of the AIME, 146(1): 54–62 DOI:10.2118/942054-G |
Arps J J., 1953. The effect of temperature on the density and electrical resistivity of sodium chloride solutions. Journal of Petroleum Technology, 5(10): 17–20 DOI:10.2118/953327-G |
Chee S, Leokprasirtkul T, Kanno T et al, 2014. A Deepwater sandface monitoring system for offshore gas hydrate. In: Proceedings of the Offshore Technology Conference. Houston, Texas: Offshore Technology Conference |
Guerin G, Goldberg D, Meltser A., 1999. Characterization of in situ elastic properties of gas hydrate-bearing sediments on the Blake Ridge. Journal of Geophysical Research, 104(B8): 17781–17795 DOI:10.1029/1999JB900127 |
Hall S A, MacBeth C, Barkved O I et al, 2002. Time-lapse seismic monitoring of compaction and subsidence at Valhall through crossmatching and interpreted warping of 3D streamer and OBC data. In: Proceedings of 2002 SEG Annual Meeting. Salt Lake City, Utah: Society of Exploration Geophysicists, 1696-1699 |
Han H, Wang Y, Li X S, et al, 2016. Experimental study on sediment deformation during methane hydrate decomposition in sandy and silty clay sediments with a novel experimental apparatus. Fuel, 182: 446–453 DOI:10.1016/j.fuel.2016.05.112 |
Higley P, Woolsey J R, Goodman R, et al, 2005. Support of Gulf of Mexico Hydrate Research Consortium:Activities to Support Establishment of a Sea Floor Monitoring Station Project. United States: The University of Mississippi, 9-14 |
Huang L, Su Z, Wu N Y., 2015. Evaluation on the gas production potential of different lithological hydrate accumulations in marine environment. Energy, 91: 782–798 DOI:10.1016/j.energy.2015.08.092 |
Jin G R, Xu T F, Xin X, et al, 2016. Numerical evaluation of the methane production from unconfined gas hydrate-bearing sediment by thermal stimulation and depressurization in Shenhu area, South China Sea. Journal of Natural Gas Science and Engineering, 33: 497–508 DOI:10.1016/j.jngse.2016.05.047 |
Paull C K, Dillon W P., 2001. Natural Gas Hydrates:Occurrence, Distribution, and Detection. Washington D C, USA:American Geophysical Union Geophysical Monograph Series,: 211–235 |
Priegnitz M, Thaler J, Spangenberg E, et al, 2015. Characterizing electrical properties and permeability changes of hydrate bearing sediments using ERT data. Geophysical Journal International, 202(3): 1599–1612 DOI:10.1093/gji/ggv245 |
Riedel M, Long P E, Collett T S., 2006. Estimates of in situ gas hydrate concentration from resistivity monitoring of gas hydrate bearing sediments during temperature equilibration. Marine Geology, 227(3-4): 215–225 DOI:10.1016/j.margeo.2005.10.007 |
Spangenberg E, Priegnitz M, Heeschen K, et al, 2015. Are laboratory-formed hydrate-bearing systems analogous to those in nature?. Journal of Chemical & Engineering Data, 60(2): 258–268 |
Stenvold T, Eiken O, Zumberge M, et al, 2006. High-precision relative depth and subsidence mapping from seafloor water pressure measurements. SPE Journal, 11(3): 380–389 DOI:10.2118/97752-PA |
Sun J X, Ning F L, Li S, et al, 2015. Numerical simulation of gas production from hydrate-bearing sediments in the Shenhu area by depressurising:the effect of burden permeability. Journal of Unconventional Oil and Gas Resources, 12: 23–33 DOI:10.1016/j.juogr.2015.08.003 |
Takahashi H, Asakawa E, Hayashi Tet al, 2011. Development of ocean bottom multi-component seismic system for methane hydrate dissociation monitoring. In: Proceedings of the Offshore Technology Conference. Houston, Texas, USA: Offshore Technology Conference |
Wang X J, Wu S G, Lee M, et al, 2011. Gas hydrate saturation from acoustic impedance and resistivity logs in the Shenhu area, South China Sea. Marine and Petroleum Geology, 28(9): 1625–1633 DOI:10.1016/j.marpetgeo.2011.07.002 |
Yokoyama T, Nakatsuka Y., 2015. Monitoring system of seafloor deformation for methane hydrate production test. The Japanese Geotechnical Society, 63(2): 26–29 |
Zhang H Q, Yang S X, Wu N Y, et al, 2007. Successful and surprising results for China's first gas hydrate drilling expedition. Fire in the Ice, 7(3): 6–9 |