地下储气库非均质盖层渗漏通道形成过程数值模拟分析

滕振超,张品金,贾善坡,温曹轩

(1.东北石油大学土木建筑工程学院,黑龙江 大庆 163318;2.东北石油大学非常规油气研究院,黑龙江 大庆 163318; 3.长江大学城市建设学院,湖北 荆州 434023)

在天然气地下储气库中,盖层的封闭性直接影响着地下储气库的封存能力[1]。经过学者们近几十年的研究,盖层的封闭能力研究取得巨大进步。截至目前,在常规盖层的封闭机理研究[2-4]和评价方法上都取得了十分丰富的研究成果,盖层封闭性的本质就是渗透性大小,渗透性是影响流体在盖层中流动的重要物性参数。但盖层的封闭性不只与毛细管压力、渗透率等参数有关,盖层岩石的非均质性对于盖层的封闭性也至关重要[5-8]。在复杂的地质演化过程中,盖层受沉积、构造、成岩、地貌等多种因素的影响,盖层中往往夹带发育砂体,使盖层呈现出复杂的非均质性,从而影响到盖层的封闭能力(见图1)。在实际情况下,因为难以准确表述出盖层的非均质特性,往往将盖层视为均质的。由于这些砂体渗透性好且分布具有随机性,所以这些高渗透性砂岩会对盖层的稳定性造成极大的影响,甚至导致盖层突破而产生不可挽回的损失和灾难。

图1 地下储气库示意图

目前,许多学者通过大量的试验模拟研究,从盖层的岩性、厚度、连续性、分布等方面对盖层的稳定性进行了评价,并取得了一定的研究成果[9-10],但针对盖层非均质的研究较少。高帅等[11]通过对层状非均质性盖层进行建模,发现非层状非均质模型能够加速CO2在盖层中运移,层状非均质模型能够减缓CO2在盖层中的运移;郭艳琴等[12]在前人研究的基础上,在垂向叠置和横向叠置上将苏东南地区盒8段砂体进行划分,垂向叠置主要为块状厚层、多层叠置、分段互层和薄互层4个主要类别,横向叠置主要为垂向叠置、孤立砂体、水平搭接和切割叠置4个主要类别,并对各类叠置关系的渗流特征进行了分析;King[13]通过研究重叠砂体之间的连通性,发现存在一个砂地比阈值(即逾渗阈值),当砂地比高于该阈值时,砂体之间会发生连通;Jackson等[14]通过一个砂-泥岩互层的三维模型,得出垂直方向阈值为0.5、水平方向阈值为0.28,即当盖地比低于0.5时,盖层在垂向不连通。

上述研究主要针对储层中砂体之间的连通,储层中砂体的占比相对较高,容易形成渗漏通道。但盖层中砂体的占比较低,所以上述研究对于探究盖层中砂体之间的连通性关系时具有一定的局限性。盖层中流体渗漏通道如图2所示。图2中显示了一种情况,盖层中砂体的占比较低,但其随机分布的位置仍可能会形成流体的渗漏通道。针对盖层中砂体占比较低的情况,本次研究考虑盖层的盖地比小于0.5时,从盖层中砂体的空间分布、渗透系数以及砂体之间距离3个方面分析随机分布的砂体对非均质盖层中流体的渗漏通道的影响。

图2 盖层中流体渗漏通道

1 渗流-应力耦合理论方程

盖层岩石作为一种多孔介质,主要是由固体物质组成的骨架和骨架间的大量微小孔隙组成。流体在岩土中渗流的过程中存在着复杂的受力环境,需要考虑渗流场和应力场的相互耦合作用。一方面,多孔介质的骨架受到外力变形导致其孔隙度改变,从而影响流体的渗透性。另一方面,流体的孔隙压力的变化又会改变多孔介质的应力状态,从而导致多孔介质的物理力学性能发生变化[15]。

1.1 达西定律质量守恒方程

通过渗流力学理论,可以得出单相流体在岩土中渗流的连续性方程:

(1)

其中:ρf为流体密度;u为渗流速度;Qm为内部及外部的流体源;φp为岩石孔隙度;t为时间。

达西运动方程为

(2)

其中:K为流体的渗透系数;μ为流体粘滞系数;∇pf为驱替压力梯度;g为重力加速度;z为垂向坐标。

弹性储存方程为

其中:S为弹性储存系数;pf为孔隙流体压力。

考虑岩石变形对渗透性的影响,其质量守恒方程为

(3)

其中:ε为体积应变;α为Biot系数。

Biot系数可以通过实验测量,也可以根据流体模量和排水模量定义为

(4)

其中:排水模量Kd总是小于流体模量Ks,因此Biot系数总是有界的εp≤αB≤1。αB并不取决于流体的性质,而是取决于多孔介质的性质。对于较软的多孔介质,Biot系数接近1 (即Kd≪Ks);而对于刚性岩石,它接近于孔隙率(即Kd≈(1-εp)Ks)。在本模型中Biot系数的取值为0.7[16]。

1.2 渗流-应力耦合变形场方程

多孔介质的控制方程为

-∇·σ=ρavg，

(5)

其中:σ为应力张量;ρav为总密度。

根据弹性力学理论,假定岩土体为弹性体,在平面应变条件下,简化为

(6)

其中:E为杨氏模量;υ为泊松比;σii为应力分量;εii为应变分量。

平面应变分析采用小变形,εxx、εyy、εzz和剪切应变εxy、εyz、εxz与位移u和υ有关。

(7)

将方程(6)、方程(7)代入控制方程(5),得到COMSOL Multiphysics求解的方程。

2 数值建模

2.1 建模及边界处理

研究以国内某储气库为背景,通过建模分析非均质盖层的稳定性。该储气库储层为厚层灰色粉砂岩,以细砂岩为主,圈闭类型为背斜,储层埋深1 125～1 320 m,单层厚度5～28 m,累计厚度71 m,根据测井岩心的实测数据,平均孔隙度34%,渗透率896 mD,直接盖层为灰色、紫红色泥岩,单层厚度3～14.5 m,累计厚度106.5 m。

根据上述储气库地质背景建立了数值模型，模型示意图见图3。模型水平向范围1 000 m;盖层厚度110 m;储层厚度30 m;假定盖层中存在如图3所示的随机砂体,砂体之间均间隔0.1 m。为了更好地观测盖层中垂向的超孔压分布规律,设置3条测线：测线1、测线2、测线3,测线位置见图3,模型地层参数见表1。

图3 模型示意图

表1 地层参数

模型通过Comsol软件的达西定律模块和固体力学模块进行双向耦合。固体力学的边界条件为:左右两侧边界为法向约束边界,底部边界为固定边界,顶部施加上覆岩层荷载,荷载值为12.12 MPa;流体力学边界条件为:两侧为非渗透边界,底部为非渗透边界,顶部为孔隙压力边界。底部中心处设置一座注入井,以3×10-4m3/s的流速持续注入CO2,流体持续注入时间为3.5×105s。该模型主要探讨流体在盖层中渗漏通道的形成过程,为简化计算过程、节约计算时间,所以仅考虑二维平面上的单相渗流。

2.2 初始地应力平衡

在自然状态下,地层中会存在未受工程扰动的天然应力,这种力被称为地应力,或者岩体初始应力。在岩土工程数值模拟中,为了较为精确地还原岩土体的实际情况,需要求解该岩土体的初始应力场状态,这一过程即为地应力平衡。因此,地应力平衡是各类岩土工程数值模拟的第一步,对初始应力场的处理是后续的数值模拟能否成功的重要条件之一[17]。

根据有效上覆压力计算公式确定盖层上覆压力[18]。在该模型中盖层埋深1 210 m,所以施加在盖层上覆岩层荷载为12.12 MPa。在流体注入前,采用超孔压的方式对模型进行地应力平衡,获得的地应力平衡下的孔隙压力见图4，地应力平衡下的垂向位移见图5。初始地应力条件下的垂向位移数量级为10-16m,位移较小;孔隙压力随深度自上而下成线性增加,可以认为初始地应力平衡结果较好。

图4 地应力平衡下的孔隙压力(单位:Pa)

图5 地应力平衡下的垂向位移(单位:m)

3 结果处理及分析

根据上述建立的概念模型,分别建立工况A、工况B、工况C和工况D 4个数值模型对非均质盖层进行分析。在其他参数相同的情况下,工况A为理想盖层,其内部不含砂岩,旨在探讨盖层岩石的微观非均质性对流体在盖层中渗流的影响;工况B探讨盖层中随机砂体的空间位置对盖层中流体渗漏通道形成的影响;工况C探讨非均质盖层中砂体的渗透系数对盖层中流体渗漏通道形成的影响;工况D探讨非均质盖层中砂体间的距离对盖层中流体渗漏通道形成的影响。

3.1 理想盖层封闭

在理想盖层条件下,为探究盖层岩石的微观非均质性对流体在盖层中的渗流影响,将工况A分为A-1、A-2分别进行讨论。工况A-1盖层岩石为非均质岩石,盖层岩石渗透系数为5×10-20～5×10-18m2;工况A-2盖层岩石为均质岩石,盖层岩石渗透系数为2.525×10-18m2。不考虑盖层中砂体的影响,注入流体后不同时刻工况A-1与工况A-2的孔隙压力分布对比云图如图6所示，不同时刻下测线1处的孔隙压力线图如图7所示。

图6 工况A孔隙压力云图对比(单位:Pa)

图7 工况A测线1处孔隙压力对比

根据图6可以看出，随着流体的不断注入,储层中的孔隙压力不断增大;当储层流体压力达到一定值后,流体开始向盖层中渗入,因为盖层的渗透性较差,导致孔隙压力急剧减小;工况A-2孔隙压力在水平方向上分布较为均匀,工况A-1孔隙压力在水平方向分布呈现出“凹凸不平”状,这是盖层的非匀质性导致的盖层渗透能力的差异。通过图7显示的储气库垂向上孔隙压力线可以看出,孔隙压力传递到距离注入井100 m左右的位置时,孔隙压力为0 MPa。此时两种工况测得的孔隙压力有稍许差异,主要是因为盖层的非均质性导致。此时工况A表现出储气库盖层的最佳封闭场景。

3.2 盖层渗漏通道的形成

工况B为假设盖层中含有随机砂体,且砂体距离盖-储交界面处及砂体之间距离均为0.5 m。流体持续注入储气库,经过3.5×105s后盖层中孔隙压力分布云图如图8所示。3.5×105s时刻下,工况A和工况B同一测线测得的孔隙压力对比图如图9所示。

从图8可以看出，随着时间的增加,虽然砂体1与储层不连续,但因为盖层中的砂体具有良好的渗透能力,能够产生压力差,所以流体注入后源源不断的进入盖层中的砂体1中,流体经过砂体1不断流入砂体2,以此类推直到流入砂体5,此时盖层中形成了突破盖层的气体渗流路径。这是因为受盖层的封闭作用影响,流体会向运移动力和阻力差最大的方向移动,即优势渗漏通道方向。通过对比图9中工况A和工况B的测压线发现,当盖层中形成运移时,储层中的孔隙压力明显降低,这是因为流体进入到盖层中,导致储层中孔隙压力变小,此时盖层中的砂体已经成为局部储层空间,此时砂体的存在已经严重影响到盖层的封闭能力。

图8 工况B在不同时刻下的孔隙压力云图(单位:Pa)

图9 工况A、B沿各测线的孔隙压力对比

3.3 渗透率对渗漏通道的影响

为了研究盖层中砂体的渗透系数对盖层中流体运移的影响,将工况C分为工况C-1～C-6进行讨论。工况C-4为各个砂体的渗透系数不变,工况C-1～C-3为各个砂体的渗透系数分别降低一个到三个数量级,工况C-5、C-6为各个砂体渗透系数分别提高一个和两个数量级,通过观察盖层中孔隙压力的大小,分析了砂体的渗透性对流体运移的影响。不同工况下各个砂体的渗透系数见表2,3.5×105s时刻下,工况C下测线1、2、3测得的孔隙压力线图如图10所示。

表2 不同工况下各个砂体的渗透系数

图10 工况C沿各测线的孔隙压力对比

根据图10可以得出,砂体的渗透系数越大,流体越容易从储层渗入砂体中,导致砂体中的孔隙压力就越大,砂体的渗透系数和砂体中的孔隙压力呈现出非线性关系,当渗透系数增加到某一数值时,孔隙压力将不再增长。当砂体的渗透系数减小时,砂体5中孔隙压力逐渐减小到0 MPa,同时储层的孔隙压力增大,这是因为砂体渗透性降低,流体积聚在储层中导致储层的孔隙压力增大。

3.4 距离对渗漏通道的影响

盖层渗透能力较差,砂体之间距离的大小也会影响着流体在盖层中的运移,在工况D中,通过改变砂体与储层之间的距离,研究各个砂体之间距离对流体运移通道的影响。

工况B中,砂体1距离盖层-储层交界只有0.5 m;在工况D-1中,砂体1距离盖层-储层交界面增加到1.5 m;在工况D-2中,砂体1距离盖层-储层交界面增加到2.5 m。3.5×105s时刻下,工况D的孔隙压力云图见图11。3.5×105s时刻下,工况D不同测线的孔隙压力线图见图12。

图11 工况D的孔隙压力云图(单位:Pa)

根据图12可以分析得出,增加砂体1到储层的距离后,储层中的孔隙压力迅速上升,盖层中的各个砂体中的孔隙压力迅速减小,虽然只增加0.5 m,但是砂体1中的孔隙压力减少一半。这是因为随着砂体1与储层的距离增大,流体难以渗入砂体,流体积聚在储层中,导致储层的孔隙压力升高。通过工况D分析得出:因为盖层岩石的渗透性较差,砂体之间的距离越大,流体渗流的阻碍越大,从而盖层的封闭效果就越好。

图12 工况D沿各测线的孔隙压力对比

4 结论

通过借助多物理场分析软件COMSOL Multiphysics,建立二维渗流场与应力场耦合分析模型并进行计算分析,研究储气库注入流体后,含有砂体的非均质盖层中孔隙压力的分布情况,得到如下结论:

(1) 在砂体体积占盖层体积比例较低的情况下,盖层中流体渗漏通道的形成往往是各种类砂体共同作用下的结果,单一种类的砂体难以形成运移,例如“坝”状砂体的存在会导致流体的垂向运移,层状砂体的存在会导致流体的侧向运移,所以对盖层的封闭性评价时,必须同时考虑盖层微观和宏观上的岩石分布特征。

(2) 盖层中砂体的渗透性的好坏对于流体运移通道的形成至关重要。当砂体渗透性越好,流体越容易从储层渗入砂体,盖层中越易形成流体的运移通道。当砂体渗透性较差时,对流体的运移具有阻碍作用,此时盖层中难以建立起流体的运移通道。

(3) 盖层中随机分布砂体的位置对于优势运移通道的形成起到关键作用。盖层岩石的渗透性较差,当盖层中砂体之间距离过大时,流体在砂体间的渗流阻力越大,不利于运移通道的形成。因此当砂体间距离较近时,需要重视流体优势渗漏通道的形成。

(4) 盖层中砂岩的连通性和连续性是盖层封闭质量的关键,即使砂体所占比例较低,但仍然能够形成运移通道。如果砂体之间形成了流体的运移通道,那么盖层的封闭能力将受到影响。

(5) 研究仅在二维平面上进行单相的流体注入模拟分析,但实际上三维空间的情况下流体的渗漏通道更为复杂,并且在储气库建库过程中会牵扯到多相流体的驱替过程,对于以上的不足下一步还需要更为细致的研究。