增强型地热系统热-水动力-力学(THM)耦合模拟——以河北马头营凸起区为例

冯 波,曹云龙,齐晓飞 ,崔振鹏,张兰新

1.吉林大学新能源与环境学院,长春 130021

2.河北省煤田地质局第二地质队,河北 邢台 054001

3.山东省地矿工程勘察院,济南 250014

0 引言

近年来,随着化石燃料大规模使用导致的环境污染等问题的加剧以及“双碳”战略目标的提出,新型替代能源的开发和利用正受到越来越多的关注[1-2]。相较于传统化石能源,地热资源具有较大的环保优势,且具备储量巨大、分布广泛、开采稳定等特点,是一种理想的新型可再生清洁能源[3-6]。

地热能是地球内部以热能形式存在的能量,近年来地热能被广泛地开发利用[7]。干热岩是地热能中的一种类型,是指埋藏在深部不含或微含流体、高温(≥180 °C)岩体内蕴藏的可被当前技术利用的巨大能源(hot dry rock ,HDR)[8-9]。热岩型地热储层岩石基质孔隙度较低,天然裂缝连通性较差,需要对其进行人工改进,以提高热能采收率。经过改造的干热岩储层被称作增强型地热系统(enhanced geothermal systems, EGS)[10]。EGS是干热岩型地热能高效利用的一项重要技术,受到了全球多个国家的高度重视。

在EGS中,通常会通过剪切增产来提高储层渗透率,从而提高地热产量,其中水力效应和热效应对EGS开发的影响不可忽视[11]。可见,热能的开采涉及到裂隙岩体中的传热、渗流和力学变形等多场耦合作用,称为热-水动力-力学(THM)耦合[12-13]。这些多物理过程的空间演化复杂且受大量参数的影响。数值方法被认为是模拟和分析EGS耦合过程的一种可行方法[14]。利用数值方法可以有效模拟裂隙岩体中的流体运移、热交换和应力演化,对于控制采热性能、保证井筒稳定性和安全性、提高地热储层采收率、延长地热储层持续使用寿命等EGS开发的关键技术问题分析具有重要的作用[15]。

在EGS的THM耦合方面, Rutqvist等[16]基于TOUGH-FLAC建立了美国Gysers EGS示范工程的 THM模型,用来分析注水策略和注水效果。在水力剪切的基础上,Xie等[17]提出了一种在 EGS水力压裂时,估计剪切起始点、所需压力以及总剪切增长方向的模型。Dempsey等[18]开发了一个全新的、能反映剪切力作用的数值模式,将其嵌入到 FEHM (finite element和 mass transfer)中,并将其应用于Desert Peak地热田,取得了良好效果。前人也基于THM耦合理论,开展了裂隙渗透性随注水温度、热膨胀系数变化规律的半解析数值模拟研究[19]。

本次研究以马头营EGS开发场地为研究对象,基于该场地2022年注水试验得到的数据,采用TOUGH2Biot模拟软件建立了场地微裂缝流固换热-水动力-力学(THM)耦合模型,模拟了冷水注入下地热储层的响应情况。本文通过将水力响应与现场观测结果(包括井口压力、注入速率)进行比较,对模型进行了校准,分析场地前期注水后EGS响应情况,预测了冷水注入后储层渗透率的变化规律和增产带空间分布范围,以期对指导该地区未来EGS储层增产作业和可持续开发提供理论支持。

1 数值方法

1.1 模拟工具及控制方程

本研究使用TOUGH2Biot模拟器进行THM数值模拟,将Biot固结模型耦合到非等温、多组分、多相传输的TOUGH2模拟器上[20-21]。

TOUGH2代码为模拟地下流系统的热(T)和水动力(H)过程提供了可靠的开源基础。TOUGH2Biot中的TH耦合进程完全继承于TOUGH2。位移、应力和应变可以通过求解扩展力学方程得到[22]。在一次迭代过程中,通过TOUGH2计算流体压力(p)、温度(T)和相饱和度(S),再将计算出来的变量传入机械模型中作为已知变量(例如,位移、应力和应变)用来求解机械模型中的未知变量,机械过程通过应力相关的裂缝渗透率向流体提供反馈。

力学模型假设岩石可以作为弹性材料移动,并遵循胡克定律的广义版本。以应力平衡方程、相容方程和应力-应变关系为基础,结合有效应力规律,以位移(如wx、wy和wz)为主要未知变量,得到常用的Biot固结模型。考虑温度影响,建立了扩展的Biot力学模型[22],见表1。

表1 三维扩展Biot力学模型

1.2 裂隙岩体渗透率演化数学模型

在本研究中,储层裂缝模型通过随机建模生成。在模型的每个网格中存在许多假想的裂缝。将不同倾角、不同倾向裂缝的渗透率转换成渗透率张量,用来表示一组相同倾角、相同倾向裂缝的等效渗透率[23]。令A=cosβi,B=sinγi,C=sinβi,D=cosγi,则式(1)表示为

(1)

式中:k为渗透率(μm2);n为网格数目;bi为第i个网格的等效水力裂缝孔径(m);li为第i个网格的缝中流体流动距离(m);βi和γi分别为第i个网格裂缝面趋势和倾角(°)。

由式(1)得到的9个渗透率分量对于每个网格来说都是不同的,这些渗透分量(例如kx,ky和kz)在全局框架中被分解,以表示裂缝介质的等效渗透率。

断裂面的应力状态决定了其产生剪切破坏的能力,所以确定断裂面的应力状态对于研究裂缝的破坏潜力非常重要。根据前人[17]研究,作用在断裂面上的有效正应力和剪应力的关系为:

(2)

(3)

式中:j、m、s分别为裂缝面法线相对于主应力轴σx、σy、σz的方向余弦;σ′x、σ′y、σ′z分别为主应力轴的有效法向应力(Pa);τn为作用在给定网格中断裂上的剪应力(Pa)。

采用Mohr-Coulomb破坏准则定义破碎岩石的抗剪强度,公式为

Fc=|τn|-μsσ′-c。

(4)

式中:Fc为抗剪强度(Pa);μs为静摩擦系数;c为内聚力(Pa)。

由式(4)可知,当Fc> 0时,表示作用在断裂面上的应力状态满足剪切破坏准则。所以前人提出了基于静/动摩擦模型来评估剪切位移的假设[24],如式(5)所示:

(5)

式中:wd为剪切位移(m);τex为超剪切应力(Pa);Kf为剪切断裂刚度(Pa/m);μd为动摩擦系数。

岩石在剪切破坏后物理性质会发生变化,而且在覆盖层荷载的作用下,岩石会产生剪切位移,这会导致岩石的渗透率发生改变。根据前人[25]的研究成果,剪切位移与裂隙渗透率的增加有如下经验关系:

(6)

式中:Δk为渗透率增加幅度百分比;Δkmax为渗透率增加幅度百分比最大值;w5和w95分别为最大渗透率增加5%和95%对应的剪切位移。

利用公式(1)—(6),可以评价裂缝性地热储层与冷水注入相关的渗透率演化。

2 模型建立

2.1 研究区概况

研究区位于中国河北省乐亭县马头营凸起区(图1)。马头营凸起区北部为燕山褶皱带,西南为黄骅坳陷,东南为胶辽隆起区,且被乐亭凹陷、柳南洼陷、石臼坨凹陷所包围。区内存在大而深的断裂,其是热量和物质通过地壳上升的优先通道[27]。

据文献[26]修编。

工作区内地层从老到新依次为:太古宇(Ar)、寒武系(∈)、奥陶系(O)、中生界(Mz)、古近系(E)、新近系(N)和第四系(Q)。新生界为局部盖层,岩性主要为变粒砂岩和黏土岩。马头营凸起区热流值大于75 mW/m2,略高于全球平均热流值[28],地温梯度一般在30 ～ 50 ℃/km之间,中部地区最高可达70 ℃/km左右[29-30]。

研究区目标储层为太古寺组变质花岗岩,主要埋深为4 200～4 500 m,最高温度略高于150 ℃。主要矿物成分为钾长石(30%,体积分数,下同)、斜长石(33%)、石英(20%)、绿泥石(15%)和黑云母(2%)[31]。

为进一步探明深层高温地热资源,前人在马头营地热田打了多口探井。具有该岩性特征的探井温度测井曲线如图2所示。地下4 000 m最高温度可达200 ℃,温度超过150 ℃的区域超过600 km2。该地热储层为单塔子群白庙组(Arb),其主要由灰白色变粒岩及浅粒岩组成,这套地热储层温度较高、岩性致密、渗透率低、无流体,符合增强型地热储层的特征。地热区储层岩性由下至上为太古宇单塔子群白庙组变粒岩及浅粒岩、新近系泥岩、第四系松散沉积[32]。

图2 马头营隆起带钻孔深度剖面岩石岩性和温度分布特征

根据前人使用XRMI(哈里伯顿电成像)技术测得的M-2井天然裂缝遥测图(图3)分析,研究区地下4 200～4 500 m处存在大量天然裂缝,天然裂缝密度为0.02～0.26条/m,裂缝宽度在几十μm量级。地层最大主应力方向为SEE—NWW[33]。

a. 已识别裂缝的蝌蚪图;b. 裂缝群倾角方位;c. 裂缝群倾角走向。图3b、3c中百分数表示裂隙数的百分比。

2.2 研究区循环试验

从2022-06-25开始,在马头营M-1井进行了一系列不同注入压力和温度的注水测试,如图4所示,根据井口压力(WHP)将增产方案分为2个阶段:阶段Ⅰ,2022-06-25—2022-07-23, 井口压力总体维持在24 MPa,井口注入温度为35～40 ℃;阶段Ⅱ,2022-07-24—2022-08-26,井口压力较高,变化范围在3 MPa左右,井口注入温度为40～45 ℃。

图4 研究区Ⅰ、Ⅱ阶段现场增产处理过程中井口压力、注入速率的变化

本次模拟的时变注入过程基于M-1井的流量监测数据进行设置。由图4可知,M-1井的实际注入过程总体上可分为2个阶段,前期注入速率为14.08 kg/s,后期注入速率为19.47 kg/s。在整个注入过程中,由于设备维修、电动机过热等原因,M-1井共计停泵29次,其中4次时间较短(2 h以内),在本次模拟中忽略不计。故本次模拟的注入过程共考虑了25次停泵,并对间隔较短的停泵过程进行了合并。

2.3 概念模型建立

研究场地主要目的储层岩性为花岗岩,储层内含有大量天然裂隙,储层岩石物性参数见表2。储层距地表4 100 m,花岗岩储层内部根据裂隙的分布可划分为3段,其中,上段和下段均为不透水基岩,各厚100 m,中段为人工裂隙储层,厚300 m。花岗岩储层孔隙度小,含水量低,无补给来源;花岗岩储层底部为恒温恒压边界,注水过程中无位移,注水井筒忽略流体沿井筒的摩擦压力损失,其余边界采用开放流动边界。概念模型如图5所示。

图5 研究区目的层概念模型

表2 研究区储层岩性参数与裂隙变形参数

2.4 数值模型建立

2.4.1 模型空间离散

根据前人钻井测井资料和钻孔成像结果提供的信息,建立了马头营干热岩场地系统的概念模型(图5)。模型顶部位于地表以下4 100 m处。模型几何尺寸水平方向为800 m×800 m,垂直方向为500 m。z轴与最大主应力(σV)方向一致。x轴和y轴分别与中间主应力(σH)和最小主应力(σh)对齐。中心注入区(200 m×200 m×200 m)网格尺寸为10 m,周边区域网格尺寸增大至50 m。注入井位于区域中心,注入井半径为0.1 m。与现场测试结果一致,注入段位于花岗岩储层,厚度为50 m(图6)。

图6 研究区计算网格剖面图

模型在纵向上划分为3段。其中,上段和下段均为不透水基岩,各厚100 m,中段为人工裂隙储层,网格尺寸为10 m,厚度300 m。

2.4.2 定解条件

根据储层监测的实际条件,设定储层初始温度为160 ℃。区域中心初始流体压力为45 MPa,与现场得到的初始流体压力相等。压力剖面是根据静水压力指定的。初始地质力学条件来自于前人研究[32]。最大主应力为垂直应力梯度,约为25 MPa/km,对应于上覆地层的体积密度。中间主应力方向约为正东方向(x轴,σH),梯度为19 MPa/km。最小主应力方向约为正北向(y轴,σh),梯度为28 MPa/km[34]。

参考场地研究调查资料[32]可知,4 000 m深度上储层的温度及压力并无明显变化,因此将模型的温度和压力条件分别设置为160 ℃、40 MPa。为了模拟不同的注入过程,将时变Dirichlet边界条件应用于注入井。忽略了沿井筒的摩擦压力损失,可以有效地将井口高压传导到注入段。模型域的大部分边缘(即恒压恒温)采用开放流动边界,底部边界采用零热量、零质量流动边界条件。我们施加了一个限制,零位移法向底部和横向边界面,而顶边界允许自由移动。

2.4.3 初始渗透率

储层的天然裂隙分布基于前人钻孔研究(图 3)。表3列出了不同层段的裂缝密度、倾向和倾角。压裂层段天然裂隙密度为0～5条/5m,隙宽在几十μm数量级。统计显示,天然裂隙与区域构造方向一致,平均走向为76°,倾角平均为60°,与区域张性应力背景产生的裂隙结果非常一致。为了表现垂向方向上裂隙密度的差异,压裂模型每个网格中生成的裂隙数目根据表3中的裂隙密度确定。每个网格裂隙的倾向和倾角随机规律一致,按照表3裂隙密度生成。

表3 储层随机裂隙生成参数设置

由前人[26]研究可知,初始裂缝孔径b= 2 μm。根据式(1)计算初始渗透率kx、ky、kz的分布,结果如图7所示。裂缝渗透率计算模型表明,裂缝性地热储层渗透率各向异性是裂缝各向异性的函数。

图7 研究区初始裂缝渗透率在x、y、z方向上的空间分布三维图

2.4.4 储层性质及参数

储层的主要参数均通过室内实验测量,静摩擦系数和动摩擦系数的选取依据前人[35]的实验结果。根据式(5)和式(6),可通过Δkmax计算剪切断裂刚度,式(6)中的参数参考水力学实验结果[25]。THM模型的主要参数值如表2所示。

2.5 模型验证

为了解系统是否正确响应流体注入,以及边界和初始条件是否正确设置,需要对模型进行验证,以确保模型的正确性[36]。本次模型校准是通过2022-07-28的注入试验进行的,选取了注入试验中的一部分数据,注入速率从14.08 kg/s增加到了19.47 kg/s,在这期间注入经过了1次停泵,停泵时注入速率变为0。把这部分数据的注入流速、温度等数据放入建立好的模型中,利用实测注入流速模拟储层改造过程中井口压力的变化趋势,通过调整储层的物性参数将实测数据与模拟结果拟合。由于本次模拟过程中场地岩石的物性参数已经由实验室测得,数据稳定可靠,所以在模拟过程中主要调整的是由随机裂缝介质模型在储层中产生的随机裂隙;并且由于在裂隙的随机产生过程中设置的裂隙产状与实际场地大致一致,所以在模型校正时应尽量不改变裂缝的倾向和倾角。综上,在本次模型校正的过程中主要是调整随机裂隙生成模型中裂隙的隙宽,通过反复校正,模拟压力曲线与实测压力曲线吻合良好(图8)。

图8 研究区2022年进行的注入试验中观测和模拟井口压力变化对比图

3 结果与讨论

3.1 现场注入试验历史拟合

3.1.1 阶段Ⅰ

储层的压裂主要包括2个阶段。阶段Ⅰ注水速率为14.08 kg/s,持续29 d,如图9所示。模拟曲线与实测数据基本吻合,在注入早期,注水流速极不稳定,由于注水流速波动的时间间隔很短,小于模型所设置的最小时间步长;所以在模拟的过程中将其中注入速率上下波动的部分进行了合并处理,即按照平均注入速率进行设置。在29 d内井口停注多次,实测井口压力和模拟压力都呈现上下波动的趋势,说明在注入过程中裂隙之中已经充满水,且并未出现新的裂隙和裂隙开度增大的情况。

图9 研究区阶段I的模拟和实测井口压力对比图

3.1.2 阶段Ⅱ

在注水阶段Ⅱ,注水速率为19.47 kg/s,持续34 d,如图10所示。模拟曲线与实测曲线基本吻合,由于这个阶段注入速率变大,可以看到这个阶段的井口压力要高于阶段Ⅰ。在注入56 d以前,井口压力的变化规律同阶段Ⅰ一样,随着注入泵的注和停,井口压力分别呈现上升和下降的趋势。56 d以后,井口的压力减小,说明此刻储层内形成了新的裂隙或者是裂隙开度增大。

图10 研究区阶段Ⅱ的模拟和实测井口压力对比图

为了进一步解释这种现象,选取离井口水平距离10 m的网格(Ele_1)进行单独分析。由于冷水的注入,导致该单元的温度降低,压力升高,有效应力降低,这种现象可以用摩尔图和应力圆的运动来表示(图11):在注水阶段开始时,Ele_1中裂隙的初始应力状态还未达到破坏包络层,表现为井口压力随着注水停注变化;持续注水到达第56天,由于热应力和水力应力的变化,Ele_1中裂隙的应力状态移动到破坏包络层,所以表现为Ele_1的裂隙被破坏,Ele_1的渗透性被改变。

总的来说,所构建模型能够较为准确地再现M-1井注入水之后储层的应力响应。经过63 d的储层注入增产处理,该模型能够证明储层内部已经发生水力剪切破坏,注入井周围渗透率显著提高。

3.2 储层渗透率变化规律

前文已经提到,经过2个阶段的注入,储层井附近的网格已经发生破坏,本模型模拟了现场水力压裂之后储层破坏的范围(图12)。由于该研究区内最大水平主应力与中间水平主应力相差不大,所以模型表现为渗透率沿最大与中间主应力变化差别不大。沿最大水平主应力方向(约正东向)扩展的最大值为22 m,沿垂直方向有中等扩展,沿最小主应力方向(约正北向)有最小扩展。这主要是因为裂缝的最佳定向面取决于最大水平主应力方向。

a. z=235 m 的x-y平面;b. y =415 m 的x-z平面;c. x=405 m的y-z平面;d. 三维增产带区域。

在整个场地循环注水完成之后,注入井附近储层渗透率全面提高。如图12所示,注水井附近储层最大渗透率增加了29倍,x方向上的平均渗透率增加了1.50倍,y方向的平均渗透率增加了1.10倍,z方向的平均渗透率增加了1.05倍,平均渗透率增加倍数不高的原因是注水井注水时间短,增产带扩展的空间距离不大,只占储层的很小一部分。图12d为场地循环注水试验后模拟预测储层的三维增产带范围,可见渗透增强区域的体积约10万m3。

4 储层在不同注入模式下增产带预测

4.1 模拟方案

前文模拟中校正了模型的参数,使之与实际场地参数吻合。从得出的结果可以看出,由于增产时间不长和间歇注入等原因,场地储层的增产范围扩展距离不长,增产带还有进一步提升的可能。根据马头营研究区实际场地条件,注水井M-1与抽水井M-2相距200 m,所以想达到最理想的生产状态,空间增产带最好沿应力最大的方向扩展,尽可能达到200 m。根据前人[4]研究,裂缝性地热储层渗透率的提高主要是由注入冷却和压力升高共同作用引起的,并且在致密储层初期采用冷水增产是提高注水井附近储层渗透率的较好选择。因此有理由认为,合理地增大注入压力和降低初期注水温度将有利于储层渗透率的提高。

考虑温度、压力的影响,并且结合实际场地条件,设计如表4所示的模拟方案。表4的9个方案中,方案a—c探讨了注入压力对增产带扩展的影响,方案d—f探讨了温度对增产带扩展的影响,方案g—i探讨了注入速率对增产带扩展的影响,时间都设置为100 d。

表4 储层增产模拟方案

4.2 模拟结果

为了深入了解流体压力对EGS演变的作用,本次研究构建了3种不同注入压力下的热-水动力-力学耦合模型,具体的模拟参数见表4中方案a—c。在模拟过程中,逐步提升了井口压力并保持模型的热物理特性恒定不变。 模拟结果见图13。

图13 模型预测100 d后不同井口压力下z=230 m x-y平面的增产带范围

根据图13的模拟结果可知,在以20 MPa定压注入的时候,储层的增产带范围扩展到近100 m,且注入井附近裂隙开度增长了约34倍;随着压力增长到40 MPa,储层的增产带范围扩展到近200 m,且注入井附近裂隙开度增长了约47倍。

储层增产带范围随着注入压力的增大而增大,这主要归因于2个因素:1)井口压力的提升会显著减小储层的有效应力,从而进一步触发已有裂缝的剪切破裂;2)压力梯度的增强能够把冷水更有效地推进到热储层中,进而扩大了增产带的范围。 因此,三维THM耦合模拟的结果表明,在地热储层中,注入冷却和压力增高两者共同作用导致了渗透率的提高。在水力剪切的情况下,现有裂缝再次剪切激活,主要是由热应力和环境应力共同驱动的(即注入压力小于最小主应力)。然而,适中的注入压力可以诱发现有裂缝的剪切破裂,并且扩大增产区域。本文模拟结果与Yuan等[11]在美国Raft River地区进行水力压裂增产数值模拟得到的结果类似,即随着注入压力增大,储层增产带的范围也随之扩大。

同时,为了深入研究热应力对EGS演变的影响,我们对3种不同注入温度下的THM耦合模型进行了模拟,具体参数见表4中方案d—f。在模拟过程中,井口注入温度从20 ℃增加到60 ℃,并保持注入压力不变。模拟结果见图14。

从图14的模拟结果中可以发现,在注入温度为20 ℃时,储层的增产带范围扩展最大,达160 m,最大开度变化为46倍。随着注入温度提高,储层的增产带范围逐渐变小。但是值得注意的是,在方案e和方案f的注入条件下,增产带的范围相差其实并不大(图14b、c),主要差别在于注入井周围裂隙开度变化幅度。如以40 ℃注入时,注入井附近裂隙被压裂的幅度更大,开度变大的范围也比以60 ℃注入时更大。

注入温度越低,储层增产带的扩展越大。随着温度的逐步上升,储层增产带逐渐变小,这主要是由于注入冷水引起冷却收缩,减小了有效应力。另外,井口周围渗透率的提高也为冷水深入储层提供了正反馈;因此,适度降低注入温度可增强储层的水力压裂能力,从而更好地破坏裂隙。然而,过低的温度会导致水的黏度增大,进而降低流动性;故在实际EGS发展过程中,应避免过长时间地使用过冷水。在方案e和方案f模拟结果的启示下可以得出,在实际工程应用中,如果只考虑水力压裂的穿透效果,注入温度的影响考虑优先级可以适当降低。合理的应用手段是在水力压裂致密储层的初期适当地注入冷水,这是提高储层渗透率的有效手段。根据Cheng等[37]的室内试验研究,不同温度储层水力压裂在增产带扩展范围方面无明显差异,推测其结论与本文模拟差异原因可能是由于温度变化本身对增产带范围变化影响不大,模拟是结果实际过程的抽象与理想化的结果,可能会导致与实际结论产生微小误差。

由于实际场地大多数时间为定速率注入,所以为了探究注入速率对场地储层增产带的影响,我们进行了3个不同注入速率下的THM耦合模拟,具体参数见表4中方案g—i。在改变注入速率的过程中同时保持模型的热物理特性恒定不变。模拟结果如图15所示。

图15 模型预测100 d后不同注入流速下z=230 m x-y平面的增产带范围

从图15可以看出,储层增产带的范围随着注入速率的增大而增大,这是由于在注入过程中井附近储层的压力增大,即当注入速率增加时,地热井中的压力也会随之增加。高压力可以推动地层中的岩石产生裂隙,这样地热流体可以在岩石中流动,增产带的范围随之扩大。并且,持续较大的流量注入使储层中压力增大的范围变大,导致更广阔的地区产生裂隙。同时,高压力导致周围岩石中的微小裂隙进一步扩大,形成新的地热流路径。本文的模拟结果与Cheng等[37]的室内试验研究结果一致,即当注入速率增加时,储层积攒的能量也会随之提高,从而增加储层裂隙扩展范围。

5 结论与建议

1)干热岩储层在水力压裂增产过程中,若处于特定的拉应力状态下,渗透率沿最大水平主应力方向增加幅度最大,沿垂直方向增加幅度中等,平行于最小主应力方向增加幅度最小。

2)水力压裂过程中,注入流体的温度对储层改造的影响较大,注入流体的温度降低,可以提高工作流体的穿透性,使储层增产带扩展范围变大;因此在实际生产过程中,在水力压裂的初期合理采用冷水增产是提高储层渗透率的有效手段。

3)储层增产带的范围随着注入速率的增大而增大。