下伏气的天然气水合藏开采模拟初探.doc

下伏气的天然气水合物藏开采模拟初探 白玉湖*1李清平1李相方2 1 中海石油研究中心,北京,100027 2 中国石油大学北京石油天然气工程学院,北京,102249 摘要建立了开采下伏气的天然气水合藏的三维数学物理模型,模型中系统地考虑了气-水-水合物-冰相多相渗流过程、水合物分解动力学过程、水合物相变过程、冰-水相变过程、热传导、对流过程、渗透率变化等对于水合物分解以及产气和产水的影响。分析了下伏气的水合物藏在降压开采过程中气层和水合物层中的压力,温度,气、水、冰、水合物等饱和度随时间和空间的变化规律,得到了不同生产压力条件下的产气速度、产水速度,水合物分解产气速度,水合物瞬时产气比率,累计产气比率等参数的变化规律。发现,在开采具有上覆水合物层的气藏时,水合物能够提供相当的产气量,能够明显提高气井产量,延长气田稳产时间。 关键词自由气藏天然气水合物藏数值模拟降压开采 1前言 天然气水合物是由天然气和水在一定的温度和压力条件下,形成的一种非化学计量型的、似冰状的笼型晶体。保守估计,全球天然气水合物中有机碳含量占全球总量的53.3,是传统化石能源中总碳量的两倍[1]。我国南海天然气水合物的储量大约在600-700亿吨[2], 2007年南海天然气水合物的成功取样,初步证实了专家关于南海蕴藏着丰富的天然气水合物的论断。地层中水合物聚集类型分为三种[3];1水合物层上覆在含有自由气和水的气藏之上;2水合物上覆在自由水层之上;3水合物层的底层和盖层均为非渗透层。对第一种类型而言,可在开采气藏的同时开采水合物,是现阶段最有望实施开采的水合物储层。西伯利亚Messoyakha气田就是此种埋藏类型,大约有36的天然气产量来源于上覆的水合物层[4]。加拿大西北部地区的Mallik油田、美国阿拉斯加Prudhoe湾/Kuparuk河区域以及日本海Nankai海槽已经进行了现场的钻井和测试,显示了水合物层下面具有自由气层的特征。我国南海水合物成功取芯表明南海极可能存在着第一种类型的天然气水合物藏。本文开展下伏气的天然气水合物藏开采研究。 Holder等[5]较早开展了下伏气的天然气水合物藏的开采研究,Burshears[6],考克斯[7]拓展了Holder的工作,并给出了一些定性的认识。但限于当时对水合物的认识,模型较为简单,假设较多。随着对水合物研究的不断深入,水合物开采模型取得很大进展,从简单的采用经典Stefan方程描述水合物分解过程[8-11]到复杂的包含动力学过程的多相渗流模型。Yousif等[12]引入了Kim-Bishnoi[13]分解动力学模型,考虑分解后生成水的流动以及变渗透率的情况,但未考虑温度变化对于水合物分解的影响。Moridis[14]等建立了水合物开采的通用处理程序EOSHYDR2,该模型是在通用多相、多组分、热流动多目的模拟器TOUGH2基础上改进的新模式,能够模拟典型的或自然的水合物沉积下的热行为、非等温气体的释放、流 国家863计划2006AA09A209资助 *联系方式Tel. 010-84523729,Email byh_ 体的流动行为等。Tang 等[15]用TOUGH-Fx/Hydrate 对水合物分解过程进行了模拟,并与其开展的实验进行了对比,但模型中未考虑冰的影响。Sun 等[16]提出了考虑冰影响的水合物分解模型,并对一维和二维条件下的水合物分解进行了模拟。总之,目前考虑结冰及冰融化的模型并不多见,而在降压法开采过程中结冰是很重要的问题。本文采用所提出的模型,对下伏气的天然气水合物藏的降压法开采进行分析。 2模型建立 2.1 物理模型 图1给出了降压开采下伏气的天然气水合藏的示意图。在开采气藏时,生产井钻穿水合物层到达自由气藏,通过开采水合物层之下的游离气来降低储层压力,使得与天然气接触的水合物变得不稳定而发生分解。还可以通过控制气井的产量,从而可以控制水合物的分解速度,提高了水合物开采安全性。图2为根据图1建立的三维物理模型。 2.2数学模型 降压法开采水合物藏涉及到复杂的物理机制,从现阶段的实际情况出发,模型中引入了一些基本假设1只考虑甲烷气形成的SI 型水合物,不考虑盐影响;2冰的成分为纯水; 3水相和气相的渗流符合达西定律;4质量传输中,不考虑分子扩散及水动力学扩散。考虑三组分水合物、甲烷气和水四相气相、水相、水合物相、冰相的控制方程为 气相控制方程 g g g g g g g g s k p m q t φρρμ∂∇⋅∇∂ 1 水相控制方程 图1 下伏气的天然气水合物藏开采示意图 图2下伏气的天然气水合物藏开采物理模型 Fig.1 The diagrammatic sketch of production of gas from hydrate reservoir underlain by gas zone Fig.2 The physical model of production of gas from hydrate reservoir underlain by gas zone w w w w I I w w w w []k s s p m q t ρφρρμ∂∇⋅∇∂ 2 水合物相控制方程 h h h s m t φρ∂-∂ 3 能量方程 g g w w w w pg g h h I I e w g t t k k C T TC p TC p K T m H m H q t ρρμμ∂∇⋅∇∇⋅∇∇⋅∇-∆∆∂ 4 其中 w w w g g Vg h h h I I I r r 1t C s C s C s C s C C φρρρρρφ- w w g g h h I I r 1t K s K s K s K s K K φφ- g rg f g g 000g ew w ,,2ln kk p p q x x y y z z r r πρδμ---- w rw f w w 000w ew w ,,2ln kk p p q x x y y z z r r πρδμ---- 饱和度方程w g h I 1s s s s 5 毛管力方程c w g p p p - 6 水相和岩石状态方程 w w0w w w01c p p ρρ-,w g w0g0 0122p p p p c φφφ- 7 边界条件 gp 0,0,0,p t p ,gp 0,0,0,T t T g g g gn g g 0k v p ρρμ-∇,w w w wn w w 0k v p ρρμ-∇,n ,,xyz 初始条件g 0i |t p p ,h 0hi |t s s ,wi 0wi |t s s ,0i |t T T 其中g ρ、w ρ、h ρ、I ρ、r ρ分别为气相、水相、水合物、冰相和岩石骨架的密度;φ为孔隙度;g s 、w s 、h s 、I s 分别为气相,水相、水合物和冰相的饱和度;w p 、g p 、f p 分 别为水相、气相及生产井底的压力;g k 、w k 分别为气相、水相的渗透率;h m 、g m 、w m 分 别为水合物分解速率、水合物分解的产气和产水速率;w C 、Vg C 、pg C 、r C 、h C 、I C 分别为水的比热、气体的等容比热、气体的等压比热、岩石骨架的比热、水合物以及冰的比热;w K 、g K 、r K 、h K 、I K 分别为水、气、岩石骨架、水合物和冰的导热系数; w c 、c φ分别为水相、岩石孔隙的压缩系数;w0p 、g0p 分别为水相和气相的参考压力;w r ,ew r 分别为生产井的半径,供液半径;,,x y z δ为δ-函数;0 x ,0y ,0z 为生产井点的坐标;gn v ,wn v 分别为气相和水相的渗流速度;h H ∆为水合物分解的相变潜热;e q 为从储层的底层和盖层所传递的能量;I H ∆为水结冰时的相变潜热。由于水合物的特殊性需要补充的方程为 1水合物分解动力学模型 采用Kim-Bishnoi 模型描述水合物的合成和分解动力学过程[13] g d s eq m k A f f - 8 其中s A 为比面,d k 为水合物的分解速率常数,f ,eq f 为当地气体逸度和平衡时气体逸度。 2绝对渗透率模型 随着水合物的分解,储层绝对渗透率发生变化,采用幂率模型估计当地绝对渗透率[17], 2e e 0000e 11k k βφφφφφφ⎛⎫- ⎪-⎝⎭ 9 k 为当地的绝对渗透率;0φ为整体孔隙度;0k 为水合物完全分解后的绝对渗透率,是与孔隙度0φ相对应的渗透率,即最大渗透率;e φ为有效孔隙度,e 0h I 1s s φφ--,该变量把水合物和冰的饱和度考虑在内;β为确定绝对渗透率随孔隙度变化而变化的指数。 3相对渗透率模型和毛管力模型 采用修改的Brooks-Corey 模型来描述相对渗透率和毛管力[18] G 0e*rg rg g n k k s ,W 0e*rw rw w n k k s 10 c e*c ce w n p p s - 11 其中,0 rg k 、0rw k 为相对渗透率曲线的末端值;G n 、W n 分别为与气相、水相相对应的指数;ce p 为进口毛管压力,c n 表示与孔隙结构相关的指数。由于水合物的特殊性,对饱和度作了修正,采用有效流动空间上的有效饱和度的定义方法,则 e e g gr e*g e e wr gr 1s s s s s ---,e e e*w wr w e e wr gr 1s s s s s ---,g e g g w s s s s ,e w w g w s s s s gr e gr g w s s s s ,e wr wr g w s s s s 4水合物及水的相平衡模型 水合物相平衡曲线的数学表达式如下所示[14] -32-53-74-105e 53exp-43.89211734346280.776302133739303-7.2729142703050210 3.8541398590072410273.2-1.0366965682883410 1.0988218047530710exp-1.941385044645610 3.3101821339792610-22T T T T K T T p T ⨯⨯⨯⨯⨯⨯⨯⨯23-44-85.55402644938060.0767559117787059273.2-1.30465829788791108.860653166875710 T T T K T T ⎧⎪⎪⎪⎪⎪⎨⎪⎪⎪≥⎪⨯⨯⎪⎩ 12 在实际水合物存在的压力温度范围内几兆帕到几十兆帕之间,压力对水冰点的影响非常微弱[19],因而,水-冰A-I 的相平衡曲线可以采用一条直线的形式表示,即 qp T T 13 因此,在考虑冰水相变存在时,水合物藏中可能存在相态组合系统包括冰-水合物IH 、气-冰GI 、气-冰-水合物GIH 、水-水合物WH ,水-气WG 、水-气-水合物WGI 、水-冰-水合物WIH 、水-冰-气WIG 、水-气-水合物-冰WGHI 九种相态系统组合,可见,水合物藏开采过程中相态系统变化非常复杂。 5水合物相变吸热模型 水合物分解是一个吸热过程,每千克水合物分解需要吸收的热量为[20] H AT B ∆ 14 其中A 、B 为常数,1050.A J kg K -,3527000B J kg 。 3模型求解 当考虑冰的存在对水合物分解的影响时,水合物藏中存在众多的相态组合系统。如果仍然采用传统方法对系统进行求解就会导致很大的计算量。而基本变量转换方法却可以很好地处理此类问题[16]。所谓的基本变量就是从离散的控制方程中直接求解出的基本未知量,对每一个网格而言,基本变量的选择都是从系统所涉及的变量中选取,如压力、各相饱和度、温度等。在降压法开采水合物藏的系统中,压力总被当作基本变量,饱和度水、气、水合物、冰基本变量则是根据各相所处的相态来确定。其中基本变量的选取是依据系统内水相和冰相转换而确定的,表1中列出了含水相系统的三种相态变化情况。 表1 相态及与之对应的基本变量 Table 1 The phase statuses and their corresponding primary variables in water phase Cases Phase status Primary variables Pha1 water W g w h ,,,p s s T Pha2 Ice-water WI g w h I ,,,p s s s Pha3 ice I g h I ,,,p s s T 在每一个时间步内,首先确定每一个网格所处的相态组合情况,然后确定所要求解的基本变量,从而可确定所要求解的控制方程。当网格处于Pha1时,需要求解方程1、2、3、4,方程2中不用考虑冰影响的;当网格处于Pha2时,不用求解方程4;当网格处于Pha3时,需要求解方程1、2、3、4,方程2中不考虑水影响的。采用有限差分方法对控制方程进行离散,对流项采用迎风格式,扩散项采用二阶中心差分格。先显式求解水合物的饱和度,然后隐式求解气相压力,再求解相关相饱和度,最后采用隐式方法求解能量方程得到温度的分布,具体的离散形式及求解步骤见本文作者文献 [21]。模型中基本参数表2所示。此外上覆水合物层的厚度为30米,初始水合物饱和度为0.62,水饱和度为0.1,初始温度280K ,初始压力为6.0MPa ;下伏气层的厚度为30米,初始气相饱和度为0.9,水相饱和度为0.1,初始压力为6.0MPa ,初始温度为280K 。 表2 模型中的基本参数取值 Table 2 The values of main physical variables Parameters Values Parameters Values w c 1/Pa 5.0e-10 G n 1.5 c φ 1/Pa 8.0e-10 W n 4.0 w K W/K.m 0.56 β 3.0 g K W/K.m 0.07 0φ 0.3 r K W/K.m 3 0k m 2 3.97e-14 h K W/K.m 0.49 i p Pa 6.0e6 I K W/K.m 3.4 I ρ kg/m 3 900 w C J/kg.K 4211 r ρ kg/m 3 2.5e3 r C J/kg.K 840 g μ Pa.s 1.0e-5 h C J/kg.K 1800 g r m 0.1 Vg C J/kg.K 2206 h ρ kg/m 3 910 I C J/kg.K 2100 w μ Pa.s 1.0e-3 i T K 280.0 0i k 2 /mol m Pa s ⋅⋅ 3.6e4 Y X Z pg 5E064.8E064.6E064.4E064.2E064E063.8E063.6E063.4E063.2E063E062.8E062.6E062.4E062.2E06 Y X Z sh 0.60.550.50.450.40.350.30.250.20.150.10.05 Y X Z sw 0.460.440.420.40.380.360.340.320.30.280.260.240.220.20.180.160.140.12 Y X Z t 281.5 281 280.5 280 279.5 279 278.5 278 277.5 277 276.5 276 275.5 275 274.5 274 273.5 Y X Z si 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 Y X Z sg 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 图3 在开采第116天时水合物层和气层中各相饱和度、压力及温度的分布。其中图标中pg、sw、sh、sg、si和t分别代表压力、水、水合物、气、冰饱和度及温度。生产井位于图中的最左边的角上。 Fig. 3 The distributions of pressure, temperature, and saturations on the 116th day in overlying hydrate zone and underlain free gas zone in 3-D hydrate reservoir by depressurization in which the production well is located at the corner of the left hand side. Here, symbols pg, sw, sh, sg, si and t denote pressure, water saturation, hydrate saturation, gas saturation, ice saturation, and temperature, respectively. 4结果讨论 图3以井底压力1MPa为例,给出了三维下伏气的水合物藏在开采第166天时,水合物层和气层中压力、水、水合物、气、冰饱和度及温度的分布规律,其中生产井位于最左边的角上。在计算中发现,在降压开采的初期阶段,压力降主要在底部气层中传播,这是因为气层的渗透率要高于水合物层。当气层压力降低至水合物相平衡压力之下时,水合物层底部紧邻气层的水合物开始发生分解,产生气体补充到气层中,水合物层基底气相饱和度增加,压力降落开始向水合物层中传播,在水合物层/气层交界面上,存在着较大的垂向压力梯度和温度梯度,正是这种温度梯度和压力梯度驱动着热量和流体的流动,加速水合物的分解。在降压过程中,水合物层/气层交界面首先开始结冰,然后结冰区域逐渐向顶层扩展。水合物分解总体上是从水合物底层向顶层不断延伸,但也存在着井筒周围的径向分解,因此,若为了增加气井的产量则应在水合物层段进行完井生产,增加水合物分解的面积和速度。如果水合物层作为保护气藏的盖层,从安全角度出发,则不应该对水合物层进行完井生产,以免水合物该层较早分解失去封闭作用。 图4给出了井底压力为1MPa时,总产气速度和水合物分解产气速度随时间的变化规律。总产气速度是指生产井单位时间内的产气量m3/d,而水合物分解产气速度则是指单位时间内由于水合物分解而产生的气体的量m3/d。可见,在生产的初期,总产气速度很高,然后开始下降进入稳产阶段。而水合物分解产气速度在初期呈逐渐上升趋势,这是因为气层气体的不断产出,气体压力逐渐下降,从而水合物分解的驱动压力相平衡压力和气层压力之差逐渐增大,因而水合物分解逐渐增多。然而,由于水合物分解吸热导致了水合物层温度降低,使得水合物的相平衡压力降低,从而水合物分解驱动力并不是一直随着气层压力降低而逐渐增加的,而是在某个阶段达到最大值。图中在40天时即达到了水合物最大的分解速度,随后水合物分解进入了平稳的阶段,水合物产气速度基本维持在某个定值附近。当气 层压力接近井底生产压力时,水合物分解消耗了大量的能量,水合物层温度降低导致了相平衡压力的相应降低,分解驱动力则逐渐减小,分解产气速度降低。图中在170天附近,水合物分解产气速度迅速下降,总产气速度也在此时下降。图5给出了井底压力为1MPa 时,水合物瞬时产气比率和累计产气比率随时间的变化规律。水合物瞬时产气比率是指水合物分解产气速度和总产气速度的比值。累积产气比率是指水合物分解的累积产气量和总产气量的比值。从图中看出,瞬时产气比率在40天时逐渐增加至0.5,然后进入稳定区域,约在170天瞬时产气比率迅速下降至0.2左右。从累积产气比率看,水合物分解提供的气体占总产气量的35,说明上覆的水合物层能够很好地延长气藏的寿命,并增加气藏的产量。 图6和图7分别给出了井底压力为1MPa ,3MPa 和5MPa 条件下总产气速度和产水速度的比较。可以看出,随着井底压力的降低,总产气速度、产水速度等都增加。从图7可见,当井底压力为5MPa 时,在生产前期阶段0-20天产水速度为零,原因有两方面,一方面地层内的压力梯度较小,相比于气体流动能力,水相流动能力要差,而且,气藏中水相饱和度本身就很低,这进一步降低了水相渗透率和水相流动能力;此外,初期阶段水合物也没有发生分解,这两个因素导致了开采初期井产水速度较低。在井底压力为1MPa 时,到生产后期,由于部分水结冰,以固态形式存于地层孔隙之中,所以产水速度在170天时降低幅度较大。 图4井底压力1MPa 时水合物分解产气瞬时速度和总产气速度的对比 图5 井底压力1MPa 时水合物产气比率和累计产气比率 图6不同井底压力条件下产气速度的比较 图7不同井底压力条件下产水速度的比较 Fig. 4 The comparisons of total gas and dissociation gas rate under bottom pressure of 1.0MPa Fig. 5 The comparisons of instantaneous and cumulative ratios of dissociation gas volume under bottom pressure of 1.0MPa Fig. 6 The comparisons of gas rate under different Fig. 7 The comparisons of water rate under different 图8给出了井底压力分别为1MPa ,3MPa 和5Mpa 时水合物分解速度的比较。可见,水合物分解速度曲线基本上可分为三个阶段,第一个阶段为分解速度上升阶段,在这个阶段内,由于下伏气层不断产出,压力不断降低,水合物分解速度逐渐增加,当水合物分解驱动力至最大值时,水合物分解产气速度达到最大值。第二个阶段为平稳产气阶段,在此阶段,水合物分解速度较为稳定,本计算中由于几何模型较小,水合物分解产气速度略有下降。第三个阶段为产气迅速下降阶段,在这个阶段,水合物分解产气速度迅速下降,这是因为前两个阶段水合物的分解消耗了地层大量的能量,造成地层温度较低,因而水合物分解的相平衡压力降低,加之气层压力的不断消耗,共同使得水合物分解的驱动压差减小。图9和给出了井底压力为1MPa ,3MPa 和5Mpa 时水合物分解产气比率的比较。可以看出,井底压力越低,水合物分解产气比率达到最大值的时间越短,在井底压力为1MPa 、3MPa 和5MPa 时,达到最大分解产气比率所需要的时间分别约为40天、60天和200天。井底压力为1MPa 时,水合物分解比率最高值超过0.5。在计算中也发现,在同样的开采时间段内,井底压力越低,水合物累计产气比率越高。可见,在开采具有上覆水合物层的气藏时,水合物能够提供相当的产气量,能够明显提高气井产量,延长气田稳产时间。 5结论 建立了下伏气的天然气水合物藏开采的三维数学物理模型,模型中系统地考虑了气-水-水合物-冰相多相渗流过程、水合物分解动力学过程、水合物相变过程、冰-水相变过程、热传导、对流过程、渗透率变化等对于水合物分解以及产气和产水的影响。分析了下伏气的水合物藏在降压开采过程中气层和水合物层中的压力,温度,气、水、冰、水合物等饱和度随时间和空间的变化规律。发现在降压开采的初期阶段,压力降主要在气层中传播,当气层压力降低至水合物相平衡压力之下时,水合物层底部近邻气层的水合物开始发生分解,产生气体补充到气层中,压力降落开始向水合物层中传播,在水合物层/气层交界面上,存在着较大的垂向压力梯度和温度梯度。降压过程中在水合物层/气层交界面结冰,然后结冰区域逐渐向顶层扩展。水合物分解总体上是从水合物底层向顶层不断延伸,但也存在着井筒周围的 图8不同井底压力水合物分解产气速度比较 图9 不同井底压力水合物分解产气比率比较 Fig. 8 The comparisons of dissociation gas rate under different bottom pressure Fig. 9 The comparisons of instantaneous ratio of diss- ociation gas volume under different bottom pressure 径向分解,因此,若为了增加气井的产量则应在水合物层段进行完井生产,增加水合物分解的面积和速度。如果水合物层作为保护气藏的盖层,从安全角度出发,则不应该对水合物层进行完井生产,以免水合物盖层较早分解失去封闭作用。经分析发现,在开采具有上覆水合物层的气藏时,水合物能够提供相当的产气量,能够明显提高气井产量,延长气田稳产时间。参考文献 1.Kvenvolden K A. A primer on the geological occurrence of gas hydrate. Geol. Soc. London, 1998, Spec. Publ. 137, 930 2.吴时国,姚伯初.天然气水合物赋存的对地质构造分析与资源评价. 北京科学出版社, 2008 3.Moridis G. J. and Timothy S. Collett, Strategies for gas production from hydrate accumulations under various geological and reservoir conditions. Proceedings, TOUGH Symposium 2003, Lawrence Berkeley National Laboratory, Berkeley, California, May 12-14,2003 4.Makogon, Y, F. Hydrate of natural gas, PennWell Publishing Co. Tulsa, Oklahoma, 1981 5.Holder Gerald D., Patrick F. Angert, Simulation of gas production from a reservoir containing both gas hydrate and free natural gas, 1982, SPE11105,1-4 6.Burshears M, Obrien T J, Malone R D. A multi-phase, multi-dimensional, variable composition simulation of gas production from a conventional gas reservoir in contact with hydrates. 1986, SPE 15246, 449453 7.考克斯J.L. Cox编美,曾昭懿,吕德本译,天然气水合物性质、资源与开采. 北京石油工业出版社,1988,11,98-122 8.Makogon Y F. Hydrate of hydrocarbons. PennWell Publishing Co. Tulsa, Oklahoma,1997. 9.Ji C, Ahmadi G., Smith D H. Constant rate natural gas production from a well in a hydrate reservoir. Energy Conversion and Management, 2003, 44, 24032423. 10.Naval Goel, Michael W, Subhash S. Analytical modeling of gas recovery from in situ hydrates dissociation, Journal of Petroleum Science and Engineering, 2001, 29, 115127. 11.Goodarz A., Chuang J, Duane H S. Numerical solution for natural gas production from methane hydrate dissociation. Journal of Petroleum Science and Engineering, 2004, 41, 169385. 12.Yousif M H, Abass H H, Selim M S. Sloan E D. Experimental and theoretical investigation of methane-gas-hydrate dissociation in porous media. SPE 18320, 1991, 6976. 13.Kim H C, Bishnoi P R, Heidemann R A and Rizvi S S H. Kinetics of methane hydrate dissociation. Chem. Eng. Sci., 1987, 427 16451653. 14.Moridis G J. Numerical studies of gas production from methane hydrates. SPE 87330, 2002,111. 15.Liang-guang Tang, Xiao-Sen Li et al. control mechanisms for gas hydrate production by depressurization in different scale hydrate reservoirs. Energy Fuels, 2007, 211227233 16.Sun X F, Kishore K M. Kinetic simulation of methane hydrate ation and dissociation in porous media. Chemical Engineering Science, 2006, 61 34763495 17.Civan F C. Scale effect on porosity and permeability kinetics, model and correlation. AIChE J., 2001, 472 271287 18.Lake L W. Enhanced Oil Recovery. Prentice-Hall Inc. Upper Saddle River, NJ, 1989. 19.Sloan E D. Clathrate hydrates of natural gases, Second edition. Marcel Deckker Inc., New York, 1998. 20.Kamath V. Study of heat transfer characteristics during dissociation of gas hydrate in porous media. PhD thesis, University of pittsburgh, Pittsburgh, 1983 21.Bai Y H, Li Q P, Y u XC, Feng G Z. Numerical study on the dissociation of gas hydrate and its sensitivity to physical parameters, China Ocean Engineering, 2007, 214625636 The numerical simulation on gas production from hydrate reservoir underlain by a free gas zone BAI Yuhu1, LI Qingping 1, LI Xiangfang2 1 Technology Research Dept., China National Offshore Oil Corp., Research Center, Beijing 100027, China 2 Faculty of Petroleum Engineering, China University of Petroleum, Beijing 102249, China Correspondence should be addressed to Bai Yuhu emailbyh_ Abstract A 3-D mathematical model of gas produced from hydrate reservoir underlain by a free gas ation is established, in which the effects of the flow of multiphase fluids, the intrinsic kinetic process of hydrate dissociation, the endothermic process of hydrate dissociation, ice-water phase equilibrium, the variation of permeability, the convection and conduction on the hydrate dissociation and gas and water production. The parameter variations, such as temperature, pressure, and saturations of water, ice, gas and hydrate etc. with space-time are respectively elucidated during the gas production from free gas ation. The relation between such parameters as gas rate, water rate, dissociation gas rate, instantaneous gas ratio of hydrate dissociation, cumulative gas ratio of hydrate dissociation, and time are respectively presented. The results show that the overlying gas hydrat