雄安新区容城地热田碳酸盐岩热储采灌数值模拟

马 峰,高 俊,王贵玲,3,刘桂宏,余鸣潇,赵志宏,刘金侠

1.中国地质科学院水文地质环境地质研究所,石家庄 050061

2.自然资源部地热与干热岩勘查开发技术创新中心,石家庄 050061

3.中国地质科学院地球深部探测中心,北京 100037

4.清华大学土木工程系,北京 100084

5.中国石化集团新星石油有限责任公司,北京 100083

0 引言

地热资源是21世纪最具开发潜力的可再生能源之一,可分为水热型和干热岩型[1-3]。雄安新区位于华北平原冀中坳陷中部,水热型地热资源丰富,其也是我国中东部地区地热资源开发利用条件最好的地区之一[4],包括牛驼镇地热田、容城地热田大部分区域以及高阳地热田北部。雄安地热以中、上元古界蓟县系碳酸盐岩为主体热储,具有储层厚度大、温度高、渗透性好、水质优、埋藏浅、易回灌的特点[5-6]。经过多年的地热开发,雄县地热的利用已成规模化,是我国第一个用地热替代燃煤供暖的“无烟城”[7],形成了政企合作、统一管理、集中回灌的“雄县模式”[8]。其中,雄安新区牛驼镇地热田开发利用条件好,而容城地热田相对于东部牛驼镇地热田而言,勘查程度较低。2017年以来,中国地质调查局加大了对容城地热田的勘查力度,在容城凸起及周边不同的构造位置部署地热井10余眼,最大探测深度达到4 000 m,揭露了蓟县系高于庄组热储层,并对其产能进行了评价[9-10];2018年,在核心区钻获了容城地热田产能最大的地热井,单井供暖能力达到30万m2[11];2019年,在启动区钻获了水量112.8 m3/h、井口温度71 ℃的高产能地热井,并揭露了太古宙变质岩基底[4]。近年来,对容城地热田的研究也逐渐增多,如:郭世炎等[12]对容城凸起地热田储层属性与资源潜力进行了分析;马峰等[4]对容城地热田热储空间结构及资源潜力进行了分析,认为容城地热田地热资源赋存条件好,有利于地热资源的规模化开发利用。蓟县系雾迷山组是雄安新区目前开发的主力热储层,其盖层主要为新生界,其地层结构松散,导热性差,具有良好的保温隔热作用[10]。

科学规范地开发利用地热资源,可发挥地热清洁能源的安全、稳定、高效等特点。采灌均衡是将提取完热量的尾水回灌到热储层,实现取热不耗水,这是保持热储压力、延长地热田寿命、防止环境破坏及实现地热资源可持续开发利用的有效模式[7,13-15]。地热流体数值模拟普遍用于地热流体资源的动态预测、开发设计、未来生产能力的评价及地下热储中热水系统的运动规律研究等,已成为地热流体开发和评价的重要技术方法[13-20]。大多数研究地热回灌系统性质的解析方法都基于理想化的假设,比如含水层是水平的、均质的、各向同性的和等厚的等条件,在实际应用中具有很大的局限性。在现实的对井开采系统中,回灌和开采层可能位于热储层不同的深度,且具有不同的垂向和水平渗透性能。因此,采用数值方法解决这些问题具有很大的优势[21-22]。如:彭展翔[22]利用TOUGH2.0软件建立天津雾迷山组热储的温度场数值模型,研究不同采灌条件下热储层温度场的变化趋势,并对采灌方案进行了优化;田光辉[23]利用GMS和TOUGH2.0软件平台分别对天津市温泉度假村砂岩孔隙型热储和基岩岩溶裂隙型热储的可采量进行了评价;刘东林等[24]采用TOUGH2.0软件平台评价天津东丽湖地区雾迷山组地热资源可开采潜力,计算得出了东丽湖地区雾迷山组可新增开采量;阮传侠[25]建立了天津地区雾迷山组地质模型和三维数值模型,对其压力场和温度场进行了模拟计算和预测分析;陈建兵等[26]建立了地热回灌井的温度场和渗流场的三维数值模型,预测西安市三桥地区回灌井成井后井内温度场和渗流场的变化;宋美钰等[27]利用TOUGH2.0 数值模拟软件建立了地热流体的数值模型,对天津地区雾迷山组热储压力场变化趋势进行了预测,分析得出地热资源开发利用可采热流体资源量与实际开采量的关系不大,而与回灌量有更为直接的联系;王宽[28]采用TOUGH2.0 数值模拟软件,根据研究区对井在预设情景下运转30 a的模拟结果,对河南省延津县新近系下部热储压力场、温度场进行模拟,数值模拟在浅层地热能和干热岩研究中得到了很好运用;冯波等[29]对U型地热井供暖潜力进行了数值模拟,认为U型井式闭循环地热系统水平井段长度为400～500 m,循环流速设置为80 m3/h,温度为20 ℃左右时,可实现地热能可持续开采;马子涵等[30]建立了澳大利亚库珀盆地三维分区均质渗透率模型和非均质渗透率模型,分别进行了热储温度场、流场及采热性能变化的研究。

目前雄安地区已开展了大量的数值模拟研究工作。例如:王树芳等[31]开展了牛驼镇地热田回灌与示踪试验研究,建立了水平裂隙型介质模型,模拟结果显示,400 m 以内的井间距可能会导致开采井温度快速下降,而 1 500 m 的井间距比较安全;刘久荣等[32]通过回灌试验与数值模型,提出了雄县牛驼镇地热田最优开发利用模式为集中开采和集中回灌模式;庞菊梅[33]利用FEFLOW(finite element subsurface FLOW system)软件对牛驼镇地热田南部地区进行水-热耦合数值模拟,结果显示,在现有的开采和回灌条件下,地热田雾迷山组储层开采井的地下水位呈下降趋势,10 a内地下水位下降约为 15 m,且在地热水开采的中心区域出现地下水位降落漏斗;胡秋韵等[34]通过分析容城凸起地热资源成藏模式,基于COMSOL软件模拟了不同采灌量对可采地热资源量的影响,通过在各开采区内合理布置开采和回灌井,评价了容东地区采灌均衡条件下的地热资源量。

雄安地区地热资源的规模化、可持续开发需要制定科学合理的采灌方案。本文基于COMSOL Multiphysics软件,建立了容城地热田地区整个蓟县系碳酸盐岩热储地质模型和数值模型,对单井开采模式下蓟县系热储未来50 a内温度场和压力场变化趋势进行了预测;并针对下降趋势制定了回灌方案,分析不同井间距对热储温度场和压力场的影响程度,以期为雄安地区地热资源合理开发、可持续开发提供依据。

1 研究区地热地质概况

雄安新区位于渤海湾盆地冀中坳陷北部,冀中坳陷内部凹凸相间的构造格局和断裂分布为热量的传导提供了有利条件[35-38]。研究区的凹陷区有巨厚的低热导率沉积层,而凸起区碳酸盐岩的高热导率地层为热流的侧向运移提供了有利条件[39],从而导致了热量向凸起区聚集,并在区内形成了容城地热田、牛驼镇地热田和高阳地热田,本次研究的重点区域为容城地热田(图1)。

F1. 太行山山前断裂;F2. 大兴断裂;F3. 容城断裂;F4. 牛东断裂;F5. 南苑--通州断裂;F6. 大兴断裂;F7. 燕山山前断裂。据文献[38]修改。

容城地热田主要地层包括新生界第四系、新近系,中元古界蓟县系、长城系和太古宇。蓟县系由铁岭组、洪水庄组、雾迷山组、杨庄组和高于庄组组成,其中雾迷山组为蓟县系目前开发的主力热储,在整个容城区域均有分布(图2)。新生界是碳酸盐岩热储的主要盖层,其中第四系厚度为340～420 m,新近系明化镇组厚度为400～500 m,新近系明化镇组底界埋深在755～1 000 m之间。蓟县系厚度变化较大,为400～2 500 m,最大埋藏深度超过3 500 m。

图2 雄安新区容城地热田雾迷山组热储埋深等值线

2 容城地热田数值模型

2.1 数值模拟方法

一般的数值模拟方法应用地热井简化模型。对于细长的地热井结构,其长度与半径之比可达 104以上,这是热储建模中一个极具挑战性的数值问题。为了克服这一数值难题,Al-Khoury等[40-41]和Saeid等[42]提出了一种用一维线单元模拟地热井中传热过程的伪三维井筒模型。本文进一步扩展了此模型,引入一维线单元对该三维井筒结构进行简化,考虑沿井筒轴向的渗流传热过程,井筒内流体与盖层的热交换通过等效换热系数来近似考虑[34,43](图3),其他主要假设条件[44]包括:1)储层处于完全饱和状态,且不考虑储层内气液相变过程;2)盖层处于完全干燥状态,且不考虑不同储层之间的水力联系;3)井筒内流体沿轴向流动,且同一深度流速沿径向处处相等;4)考虑井筒内流体性质如密度、黏度、热传导系数、比热容等与温度的相关性,但不考虑井筒内气液相变过程。

r1. 井径;r0. 套管后井径;r2. 灌浆后井径;Δr. 井径差;T0. 井壁温度;T1. 套管井壁温度;T2. 灌浆井壁温度;Text.外部温度;h. 深度。

本次建模采用成熟的地质体多场耦合有限元计算软件COMSOL Multiphysics。COMSOL为一款通用的工程仿真软件平台,通过在一个模块中添加其余模块的因变量,实现各个模块间的任意耦合。其优势在于用户可以耦合多个物理场、化学场,并提供强大的开放接口供用户修改控制方程。

2.2 模型范围

数值模型计算范围水平方向上为容城地热田,面积约为300 km2;垂直方向上,根据岩性差异从上往下主要考虑第四系、新近系和蓟县系,根据容城已有地热井的相关钻井资料,地热井井深在1 700 m左右,模型深度设为2 000 m。第四系和新近系的底板埋深根据已有的钻孔资料,利用surfer软件插值后导入COMSOL软件中,建立的容城地热田地质模型如图4所示。

Q. 第四系;Nm. 明化镇组;Jx. 蓟县系。

2.3 热储层参数

研究区蓟县系热储层参数按照表1取值。

表1 模型参数取值列表

深部热储中的流体密度(ρf)、黏度(μ)、导热系数(kf)和比热容(Cp,f)随温度(T)变化情况按如下经验公式[45-46]计算:

ρf=838.47+1.40T+0.003T2+3.72×10-7T3;

(1)

μ=1.38-0.02T+1.36×10-4T2-4.64×10-7T3+

8.90×10-10T4-9.08×10-13T5+3.85×10-16T6;

(2)

kf=-0.87+0.009T-1.58×10-5T2+7.98×10-9T3;

(3)

CP,f=12010.15+80.41T+0.31T2-5.38×10-4T3+

3.62×10-7T4。

(4)

2.4 网格剖分

在COMSOL Multiphysics中,综合考虑计算精度与计算效率之间的平衡,地热井及周围区域的网格精细化,最小单元尺寸5 m,边界及断裂面网格细化,模型共包含约26万个单元(图5)。

图5 雄安新区容城地热田地质模型网格剖分示意图

2.5 时间离散

由资料①可知,地热井只在每年11月15日至翌年3月15日4个月供暖期内处于开采或回灌运行状态,其余时间地热井均处于停运状态,故地热井流量可表示为

(5)

式中:q为地热井流量,m3/h;n为开采周期数;t为运行时间,月。

2.6 模型验证

1)水位埋深拟合

容城区目前有两眼地热井的水位监测数据,分别是领秀城1井和领秀城2井。领秀城1井的水位监测从2016-07-01—2017-10-30,领秀城2井的水位监测则从2013-11-15—2017-10-30。通过调整模型的水文地质参数和相关边界条件,使模拟水位与实测水位基本吻合(图6),可利用拟合得到的参数进行容城区的采灌优化设计。从水位分布云图(图7)可以看到,区内因抽水而形成了漏斗,回灌则能有效地恢复储层压力。

a. 领秀城1井,2016-07-01—2017-10-30;b. 领秀城2井,2013-11-15—2017-10-30。

图中黑色箭头表示热流运移方向。

2)温度拟合

根据钻井资料(1)河北雄安新区管理委员会综合执法局.雄安新区地热资源保护与开发利用规划(2019--2025年),雄安新区:河北雄安新区管理委员会,2019.,利用surfer软件进行插值后得到的模型地温梯度如图8所示,基岩地温梯度基本在1.5 ℃/hm左右。领秀城1井的监测水温为51 ℃,领秀城2井的监测水温为49 ℃。通过调整相关的热物性参数,使得计算得到的两眼监测井的温度与实测温度一致,从而确定热储层参数,如表2所示。

表2 研究区地层参数取值

图8 研究区新生界地温梯度

3 结果与讨论

3.1 热储开采与回灌数值模拟

为控制雾迷山组热储层地热水水位的持续下降和地热尾水排放污染,地热供暖尾水必须回灌,以保证采灌均衡。研究区采用“一采一灌”对井模式,目前已有11眼开采井和11眼回灌井,拟部署9对对井。按照图9的井位分布设定,部署钻孔的位置处也匹配1眼回灌井。井间距初步设计为600 m。开采量和回灌量均设置为100 m3/h,区内已有地热井的流量根据抽水试验,按最大涌水量考虑,具体取值如表3所示。通过调整模型的水文地质参数和相关边界进行模拟,结果如图10、11所示。图10为模拟不同开采年份的热储温度分布。从图10可以看出,随着开采年限的增加,回灌冷锋面在地下渗流场的作用下逐渐向开采井运移,从而导致某些开采井发生热突破,开采温度逐渐降低。本模型以开采井的平均温度降低1 ℃作为热突破的临界值,不同开采井的温度变化曲线如图11所示。图11中曲线的振荡幅度与时间周期有关。统计开采井的热突破时间可知:谷丰3井的热突破时间最长,为25 a;谷丰2井的热突破时间次之,为13 a,这是由于谷丰2井靠近回灌井,热突破时间相对较短;领秀城1井的热突破时间最短,约10 a左右,这是由于领秀城1井靠近回灌井领秀城4井,为避免快速发生热突破可以采用减小采灌量的方法。

表3 研究区地热井流量取值列表

图9 研究区已有井位分布

a. 领秀城1井;b. 领秀城2井;c. 谷丰2井;d. 谷丰3井。

3.2 无回灌情况下

在没有回灌井的情况下,对研究区已有的11眼开采井和拟部署的9眼开采井进行模拟。没有回灌井的情况下开采井温度通常不会降低,但长期的开采势必会造成热储水位的下降,我们以热储层地热流体开采最大允许埋深值150 m为约束条件,用来计算研究区地热资源的最大开采量。模拟期从2017年开始,模拟周期50 a,计算结果如图12所示,可以看到,在无回灌井开采40 a的情况下,在地热井比较集中的区域水位埋深较大,并且因抽水形成的漏斗区域逐渐扩大,相对孤立的地热井则会形成小范围的抽水漏斗。通过监测地热井的水位埋深变化,可以得到地热井的水位变化曲线(图13)。从图13可以看到,在40 a左右水位埋深降低至150 m,每年水位下降约1.13 m。

图13 研究区开采井水位变化曲线

3.3 有回灌情况下

在有回灌的情况下,采用“一采一灌”方案的对井系统,即1眼开采井匹配1眼回灌井,这样热储水位由于回灌的作用,基本保持了压力的平衡,唯一需要考虑的是在回灌条件下开采井的热突破问题。我们以开采井的热突破为限制条件,用来计算研究区地热资源的最大开采量。计算过程为,在已有开采井和回灌井的基础上(图9),井间距分别考虑600、800和1 000 m 等3种情况,先按照图14所示的布井方案在研究区域按最大数量布置对井系统;然后计算对井系统开采井的热突破时间,同时考虑在不同开采量下的开采温度变化,计算时长为50 a。计算结果如图15所示。最大开采量的统计结果如表4所示。

表4 研究区不同井间距对井系统的最大开采量计算

a. 井间距600 m;b. 井间距800 m;c. 井间距1 000 m。

a. 井间距600 m;b. 井间距800 m;c. 井间距1 000 m。

由表4可知:当井间距为600 m时,新区布井数量为138对,随着开采量的增加(100、200到300 m3/h),热突破时间从50 a降低到21 a,其中开采量为300 m3/h时,新区最大开采量为5.166 4×107m3/a,折合标准煤1.765 7×107t;当井间距为800 m时,新区布井数量为80对,随着开采量的增加(300、400到500 m3/h),热突破时间从51 a降低到28 a,其中开采量为500 m3/h时,新区最大开采量为5.257 4×107m3/a,折合标准煤1.796 8×107t;当井间距为1 000 m时,新区布井数量为50对,随着开采量的增加(400、500到600 m3/h),热突破时间从50 a降低到30 a,开采量为600 m3/h时,新区最大开采量为5.348 3×107m3/a,折合标准煤1.827 8×107t。计算结果表明:增大对井间距虽然减少了可布置的井对数量,但是延长了热突破时间,增加了开采量,导致能开采更多的热量。若以开采温度降低1 ℃为热突破标准,则在热突破时间较长的情况下能开采更多的热量。即在50 a 的开采年限里,开采温度最终下降3和5 ℃ 的热突破时间相差不多,最终开采的热量相差不大(图11)。

4 结论

1)通过COMSOL建立了容城地热田群井回灌数值模型,根据已有的地质资料和相关的监测数据对模型进行了识别和验证,温度和水位拟合结果较好,说明模型可以用于后续的计算。根据现有的采灌井布置方案,领秀城1井的热突破时间最短,谷丰2井的热突破时间较短,建议这两眼井可以采用减小采灌量的方法,以延长热突破时间。

2)采灌均衡情况下,采用“一采一灌”方案的对井系统,考虑井间距分别为600、800和1 000 m 等3种情况,增大对井间距。虽然可布置的对井数量减少了,但是其开采量的增加及热突破时间的变慢能开采更多的热量,若以开采温度降低1 ℃为热突破标准,在开采量较小,热突破时间较长的情况下能开采更多的热量。在50 a的开采年限里,开采温度最终下降3 ℃和5 ℃的热突破时间相差不大,最终开采的热量相差不多。建议在实际工程中,通过数值模拟方法选择适宜的对井间距,保证开采量的增加及减慢热突破时间,开采更多的热量。