水驱气藏动态储量和水侵量计算新方法*

闫正和 石军太 秦 峰 洪舒娜 白美丽

(1. 中海石油深海开发有限公司 广东珠海 518000； 2. 中国石油大学(北京)油气资源与探测国家重点实验室 北京 102249；3. 石油工程教育部重点实验室(中国石油大学) 北京 102249 )

国内大多数常规气藏都属于不同程度的水驱气藏，边底水活跃的气藏占40%～50%[1]。在水驱气藏的开发过程中，随着气体的产出和地层压力降低，水体逐渐侵入气藏。水体的侵入，一方面可以有效补充地层能量，另一方面也会提高气藏废弃压力、降低气藏采收率。因此，水驱气藏水侵量及储量计算是气藏合理高效开发的基础，对水驱气藏后期生产调整具有重要意义。

常规气驱气藏的动态储量方法[2-5]已相对成熟，但在水驱气藏储量方面适应性较差，如应用常规物质平衡方法计算水驱气藏动态储量往往会导致计算结果偏大[6]。国内外学者针对水驱气藏水侵量及动态储量计算开展了大量研究，推导了若干水侵量计算方法，主要可分为稳态水侵模型[7]和非稳态水侵模型[8-11]两类。但是上述模型假设条件较多，计算过程相对繁琐，且计算所需的水体参数如水侵角、水体内外边界半径等往往难以准确确定，不利于实际应用。因此，一些学者基于物质平衡原理和最优化思想，提出了水侵量计算方法与数学模型[12-19]。如唐圣来 等[12]应用水驱气藏物质平衡方程和变流压下产量响应函数，通过拟合产量、流压、静压等生产动态数据，经过多次循环迭代得出强水驱气藏的动态储量和水侵量；王怒涛 等[13]结合数值反演法，提出了一种自适应遗传算法用于求解地质储量、水侵量、水侵系数以及水体大小；杜凌云 等[14]针对常规物质平衡法中的差值法、图版法、视地质储量法进行分析，建立了综合考虑这3种计算方法的水侵量最优化模型。也有一些学者提出了水驱气藏储量和水侵量简易计算方法[20-22]，基于气井二项式产能方程、试井渗透率、气水相渗曲线的水驱气藏水侵量和储量计算方法[23]，基于水驱特征曲线的储量计算方法[24]。但是，以上最优化计算方法通常需要利用计算机编程来实现，现场应用较为困难，而简易方法大多需要数据拟合得出经验相关式，理论基础薄弱。

因此，笔者基于气藏压力下降引起水体压力下降、水体孔隙体积收缩和水体中水的弹性膨胀导致水侵的机理，应用水体物质平衡原理，推导出了水侵量计算模型，并结合气藏物质平衡方程，建立了水驱气藏动态储量和水侵量计算新方法，并进行了验证和实际应用。该方法可同时计算出气藏的动态储量和水侵量，也可判断产气驱动能量占比的变化，为气藏防水治水提供指导。

1 新方法的提出

水侵量的计算是水驱气藏动态储量计算的先决条件，即在计算水驱气藏动态储量之前，须提前用其他方法计算水驱气藏的水侵量。关于水侵量的计算，一直是水驱气藏动态储量计算的难点。本文将应用物质平衡原理，以水驱气藏的水体为研究对象，将有限水体的弹性膨胀量与水侵量关联起来(图1)，结合水驱气藏的物质平衡方程，形成水体-气藏耦合物质平衡方程。

图1 水驱气藏示意图Fig .1 Schematic diagram of gas reservoir with water influx

水驱气藏物质平衡方程的通式为

GBgi=(G-Gp)Bg+We-WpBw+ΔVp+ΔVw

(1)

式(1)中：G为水驱气藏天然气原始地质储量，108m3；Bgi为原始地层压力下天然气体积系数，无因次；Gp为累计产气量，108m3；Bg为当前平均地层压力下天然气体积系数，无因次；We为水侵量，108m3；Wp为累计产水量，108m3；Bw为水的体积系数，无因次；ΔVp和ΔVw表示孔隙体积的缩小量和气藏中孔隙水的弹性膨胀量，108m3。此处体积单位均采用亿方，便于数值计算。

ΔVp和ΔVw通常采用平均地层压力降幅的简化表达式，这种简化在孔隙压缩系数和水的等温压缩系数数值较小的情况下，可以满足计算要求，但是当其数值较大，简化表达式将造成较大误差，本文将应用孔隙压缩系数和水的等温压缩系数的定义式，严格推导ΔVp和ΔVw的精确表达式。

1.1 气藏中孔隙体积收缩和孔隙水膨胀量表达式推导

应用孔隙压缩系数的定义式

(2)

式(2)中：Cp为气藏孔隙压缩系数，MPa-1；p为平均地层压力，MPa；Vp为孔隙体积，108m3。

对式(2)从原始地层压力到目前平均地层压力进行积分得

Vp=VpieCp(p-pi)

(3)

压力从原始地层压力pi降为p，孔隙体积缩小量ΔVp可表示为

ΔVp=Vpi[1-eCp(p-pi)]

(4)

同理，应用水的等温压缩系数定义式

(5)

式(5)中：Cw为水的等温压缩系数，MPa-1；p为平均地层压力，MPa；Vw为孔隙水体积，108m3。

对式(5)从原始地层压力到目前平均地层压力进行积分得

Vw=VwieCw(pi-p)

(6)

压力从pi降为p，孔隙中原始含水的体积膨胀量ΔVw可表示为

ΔVw=Vwi[eCw(pi-p)-1]

(7)

1.2 水侵量表达式推导

假设水体孔隙体积为Ve，水体孔隙中充满水，因此水体中水的体积也为Ve，水侵量We为水体压力从pi降为p时挤出的水量。在气藏流体产出过程中，气藏和水体间为一非平衡系统，水体压力高于气藏压力，即水体压降小于气藏压降。假设水体压降与气藏压降的比值为β，其数值与气井的平均产气量大小有关。平均产气量越小，平衡条件越容易满足，β越接近1；平均产气量越大，非平衡程度越严重，则β越小(图2)。β的取值大于0而小于等于1。β=1，表明水体和气藏达到平衡，即水体平均地层压力和气藏平均地层压力同步下降。

图2 气藏压降、水体压降与平均产气量之间关系的示意图Fig .2 Schematic diagram for the relationship among reservoir pressure drop,aquifer pressure drop,and average gas production rate

水体中孔隙体积收缩和水体膨胀量可参照式(4)和式(7)的推导过程。假设水体孔隙压缩系数与气藏孔隙压缩系数相同，因为气藏和水体一般处于同一储层，因此该条件不难满足。根据水体孔隙压缩系数和水体中水的等温压缩系数公式，可得水体中孔隙体积收缩量为

ΔVpa=Ve[1-eCpβ(p-pi)]

(8)

水体中水的膨胀量为

ΔVwa=Ve[eCwβ(pi-p)-1]

(9)

那么水侵量为水体中孔隙体积收缩量加上水体中水的膨胀量，表达式为

We=Ve[eCwβ(pi-p)-eCpβ(p-pi)]

(10)

假设水体倍数为n，则Ve=nVpi，则

We=nVpi[eCwβ(pi-p)-eCpβ(p-pi)]

(11)

鉴于水侵初期水体压降只波及到了部分水体，在水侵的过程中，水体压降波及是个动态过程，这种动态变化是气藏中天然气逐渐被采出所导致的，气藏原始状态累计产气量为0，水体压降波及体积为0，水侵量为0；随着累计产气量的增大，水体波及体积增大，直到全部水体被波及，因此选用Gpb进行修正。选用幂函数的原因在于既要满足水体波及体积随着累计产气量的增大而增大的特征，又要满足累计产气量为零时水体波及体积为0。因此在式(11)的基础上加入对累计产气量的修正，表达式为

We=GpbnVpi[eCwβ(pi-p)-eCpβ(p-pi)]

(12)

式(12)中：We为水侵量，108m3；Gp为累计产气量，108m3；b为水体波及体积修正因子，无因次。若b=0，表明水体波及速度极快，即不考虑水体波及的动态过程。若b=1，表明水体波及体积和累计产气量是正比关系。

1.3 水驱气藏动态储量和水侵量计算方法

将ΔVp的表达式(式(4))、ΔVw的表达式(式(7))和We的表达式(式(12))代入式(1)可得

GBgi=(G-Gp)Bg+GpbnVpi[eCwβ(pi-p)-eCpβ(p-pi)]-WpBw+Vpi[1-eCp(p-pi)]+Vwi[eCw(pi-p)-1]

由于

那么

(13)

式(13)中：Swi为气藏原始含水饱和度。

整理得

GpBg+WpBw=G{(Bg-Bgi)+

(14)

若令

F=GpBg+WpBw

(15)

Eg=Bg-Bgi

(16)

(17)

(18)

E=Eg+Efw+Eaq

(19)

那么

F=GE

(20)

式(15)～(20)中：F为产出流体的地下体积，108m3；Eg为气体弹性膨胀能，无因次；Efw为气藏中孔隙体积收缩和气藏中水的弹性膨胀能，无因次；Eaq为水体孔隙体积收缩和水体中水的弹性膨胀能，无因次；E为天然气产出综合弹性膨胀能，无因次。

应用多个实测压力点，以F为纵轴，E为横轴做图，先设置b=0，试算水体倍数n和β的值，β值通常取1，主要改变n值，使所有数据都近乎落在同一条直线上，且截距尽可能接近0。若截距大于0，则增大n的取值，若截距小于0，则减小n的取值。最终得出最接近零截距的n值即为最终的n值，对应直线斜率值即为该水驱气藏的动态储量。

将式(20)变形为

(21)

应用多个实测压力点，以F/E为纵轴，Gp为横轴做图，先设置b=0，试算水体倍数n和β的值，使所有数据都近乎落在同一水平线上，对应的纵轴数值即为该水驱气藏的动态储量。若直线向下偏离水平线，则增大n值，若直线向上偏离水平线，则减小n值，直到直线非常接近水平线，得出最终的n值。若在b=0的情况下，无论如何调整n和β，都无法满足水平线，那么再调整b值。

以上2种确定方法可结合使用，在F与E的直角坐标图中，使所有数据都近乎落在同一条直线上，且截距尽可能接近零；在F/E与Gp的直角坐标图中，使所有数据都近乎落在同一水平线上；且F与E的直角坐标图中直线点的斜率近乎等于F/E与Gp的直角坐标图中水平线纵轴数值。那么该数值即为水驱气藏动态储量G。

根据最终确定出的b、n和β的值，可以计算出水侵量，公式如下：

(22)

计算出动态储量和水侵量之后，就可以评价天然气开采过程中各种驱动能量在整个驱动能量中占比的变化，将式(19)变形可得

(23)

1.4 水驱气藏动态储量和水侵量计算方法简化

对于非异常高压的水驱气藏，若Cp和Cw数值很小，约小于5×10-4MPa-1，则式(14)可以简化为

GpBg+WpBw=G[(Bg-Bgi)+

(24)

令nβ=n′，则最终表达式可以简化为

GpBg+WpBw=G[(Bg-Bgi)+

(25)

根据式(15)～(19)，Efw+Eaq将简化为

Efw+Eaq=

(26)

式(26)中：n′为等效水体倍数。通过这样的简化，将简化前的3个调整参数b、n和β减少成了简化后的2个调整的参数b和n′。

根据最终确定出的b和等效水体倍数n′的值，可以计算出水侵量，公式如下：

(27)

2 新方法的验证

首先用理论数值模型对本文方法进行验证，然后利用文献中已有的实例进行验证。

2.1 理论数值模型验证

采用CMG软件构建理论数值模型，水驱气藏埋深3 000 m，原始地层压力为30 MPa，储层温度为100 ℃，气层厚度为32 m，孔隙度为0.2，渗透率为8 mD，原始含水饱和度为0.5，孔隙压缩系数为0.000 32 MPa-1，水的等温压缩系数为0.000 445 9 MPa-1，水的体积系数为1.030 15，天然气相对密度为0.58，气藏外边界半径为400 m，水体模型选用Fetkovitch模型，水体外边界与气藏外边界的比值分别设置为1.7和1.9，那么计算得出对应的水体倍数为1.89和2.61。由于CMG软件中设定的地表温度为15.6 ℃，因此计算得出该气藏天然气原始体积系数为0.0042 98，采用容积法计算得出该气藏的天然气原始地质储量为3.742 2×108m3。应用建立的CMG理论数值模型，定产气量5×104m3/d以一个月的时间步长生产了10年，输出累计产气量、累计产水量、平均地层压力、平均视压力(p/Z)和累计水侵量的动态数据。水体外边界与气藏外边界的比值为1.7和1.9的2个案例，以年为时间步长的动态数据和水侵量计算表如表1和表2所示。

表1 水体外边界与气藏外边界的比值为1.7时动态数据与水侵量计算表Table 1 Calculation table of dynamic data and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.7

表2 水体外边界与气藏外边界的比值为1.9时动态数据与水侵量计算表Table 2 Calculation table of dynamic data and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.9

然后假设天然气原始地质储量、水体倍数和水侵量未知，已知原始含水饱和度、孔隙压缩系数、水的等温压缩系数和水的体积系数，并采用CMG数值模拟模型输出的累计产气量、累计产水量和平均视压力(p/Z)数据，应用本文建立的方法，反求该气藏的天然气原始地质储量和水侵量动态。若数据点满足本文提出的直线关系和水平线关系，且求出的动态储量与模型真实储量(容积法计算得出)相差不大，水侵量计算结果与理论数值模型输出结果接近，则证明本文提出的水驱气藏动态储量和水侵量半解析方法是合理和可靠的。

理论数值模型验证过程中，b=0且β=1就能很好地拟合出直线和水平线，因此对于理论数值模型，可以认为水侵初期就波及到了整个水体，水体和气藏实时达到平衡状态。水体边界与气藏外边界的比值为1.7和1.9的拟合图如图3、4所示。可以看出这2个例子F与E的直线关系非常好，且几乎穿过原点(图3a、4a)；F/E与Gp的水平线完全水平(图3b、4b)。调整确定的水体倍数n分别为1.858 6和2.605 2，计算的天然气动态储量分别为3.725×108m3和3.747×108m3，与理论数值模型真实天然气地质储量3.742 2×108m3相差不大。如表1、表2、图3c和图4c所示，新方法计算的水侵量与理论数值模型输出的真实水侵量几乎完全吻合。理论数值模型验证结果表明，本文建立的水驱气藏动态储量和水侵量计算新方法在理论上是正确的。

图3 水体外边界与气藏外边界比值为1.7时动态储量与水侵量计算结果Fig .3 Calculation results of reserves and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.7

图4 水体外边界与气藏外边界比值为1.9时动态储量与水侵量计算结果Fig .4 Calculation results of reserves and water influx when the ratio of water body outer boundary to gas reservoir outer boundary is 1.9

2.2 实例验证

实例的数据选自文献[20]，气藏原始地层压力为78.789 MPa，储层温度为111.85 ℃，综合压缩系数为0.000 542 69 MPa-1，原始含水饱和度为0.562 2[25]，假设水的体积系数为1.071 36，水的等温压缩系数为0.000 235 MPa-1，那么计算出孔隙压缩系数为0.000 105 5 MPa-1，累计产气量、累计产水量、压力、视压力和实际水侵量等动态数据[20]见表3所示。应用本文方法，计算出Eg、Efw和Eaq、E和F的值(表3)。

表3 某实例水驱气藏动态数据[20]与水侵量计算表Table 3 Dynamic data and water influx calculation table of a water drive gas reservoir

应用本文方法，F与E拟合图、F/E与Gp拟合图和水侵量计算结果如图5所示。可以看出，F与E的直线关系非常好，且几乎穿过原点(图5a)；F/E与Gp的水平线近乎水平(图5b)。调整确定的水体倍数n为5.03，水体压降与气藏压降的比值β为1，b值为1.04，求得该气藏天然气动态储量为19.232×108m3，与文献中该气藏实际天然气地质储量19.321×108m3相差不大。如表3和图5c所示，新方法计算的水侵量与文献中该气藏的实际水侵量吻合程度也很高。该实例验证结果表明，本文建立的水驱气藏动态储量和水侵量计算新方法在实际应用中也是可靠的，可以推广应用。

图5 文献中气藏实例动态储量与水侵量计算结果Fig .5 Calculation results of reserves and water influx of the gas reservoir in the literature

3 新方法的应用

番禺某气田分布着多个中高孔、中特高渗边底水气藏，气藏A被A1井控制，气藏B被B1井控制，气藏C被C1井控制。这些井所控制的气藏的物性参数见表4。单井的累计产气量、累计产水量、平均地层压力和平均视压力数据见表5。

应用本文方法计算出各井的Eg、Efw和Eaq、E和F的值(表5)。F与E拟合图、F/E与Gp拟合图和水侵量计算结果如图6所示。可以看出，3口井F与E的直线关系都非常好，且都几乎穿过原点(图6a)；3口井F/E与Gp的曲线接近水平线(图6b)。3口井调整确定的水体倍数n、水体压降与气藏压降的比值β和b值，以及求得的各井控制的动态储量见表4。本文方法计算的各井水侵量见表5和图6c。A1、B1和C1井产气驱动能量占比变化曲线如图6d、e和f所示。可以看出A1井和B1井天然气膨胀能和水体能量相当，而C1井水体能量是天然气生产的主要驱动能量，一直维持在90%左右。此外，A1井和B1井水体能量在开发过程中是变化的，A1井从最初的45%上升为53%，水体能量越来越强；B1井水体能量从最初的48%降低为42%，水体驱动能量逐渐变弱。

表4 番禺某气田各井物性参数表Table 4 Physical property parameters of three wells in Panyu gas field

表5 番禺某气田3口井动态Table 5 Production performance of three wells in Panyu gas field

图6 番禺某气田3口井动态储量与水侵量计算结果Fig .6 Calculation results of reserves and water influx of three wells in Panyu gas field

4 结论与建议

1) 本文建立的水驱气藏动态储量和水侵量计算新方法，可同时计算出气藏的动态储量和水侵量，适用于同一压力系统的气藏、区块或单井，除原始地层压力外，还需要至少2个静压值才能计算，数据点越多，且采出程度越高，储量和水侵量计算结果越可靠。

2) 综合产出流体地下体积F与综合弹性能E指示曲线和F/E与Gp指示曲线，增加了拟合的精度和解释结果的唯一性，理论数值模型和实例验证均表明：F与E指示曲线确实符合穿过原点的直线，斜率即为气藏的动态储量，F/E与Gp指示曲线的确为一水平线，纵轴数值即为气藏的动态储量，且水侵量计算结果与理论数值模型和实际水侵量吻合程度很高，从理论和实践均证明所建立的水驱气藏动态储量与水侵量计算方法是正确的。

3) 建立的水驱气藏动态储量和水侵量计算方法也可判断天然气弹性能、气藏岩石和孔隙水弹性膨胀能、水体能量在整个驱动能量中的占比的变化，评价水驱气藏水侵能量强弱，为气藏防水治水措施提供依据。

4) 对于异常高压水驱气藏，水体中的溶解气的弹性膨胀能量不可忽略，考虑水体中溶解气的影响将是异常高压水驱气藏动态储量和水侵量计算的下步研究方向。