摘要
地热资源粗放式开发易引起水位下降，回灌是缓解水位下降的有效措施之一。天津市东丽湖地区碳酸盐岩热储回灌效果好，利用地表水回灌对热储层进行额外补充是可行的方案。本文聚焦地表水长期回灌对储层的影响，利用 COMSOL建立了潘庄凸起构造区渗流数值模型，并对地表水回灌条件下东丽湖地区雾迷山组碳酸盐岩热储未来50年的温度场、压力场进行了预测。结果表明：DLB-01、DLB-02和DLB-03井持续地表水回灌50 a后，不会对周边开采井产生热突破，但DL-71井温度接近设定的热突破界限，因此需留意对附近开采井的水温监测；地表水回灌对周边热储水位提升明显，DL-71井在地表水回灌50 a后，水位上升并稳定在153.5 m左右。
Abstract
The extensive development of geothermal resources is prone to cause the water level to drop. Reinjection is one of the effective measures to alleviate the drop in water level. The reinjection effect of carbonate rock thermal reservoirs in the Dongli Lake area of Tianjin is good. It is a feasible solution to supplement the thermal reservoirs by reinjection of surface water. This paper focuses on the influence of long-term surface water reinjection on the reservoir. A seepage numerical model of the Panzhuang Uplift tectonic area is established using COMSOL, and the temperature field and pressure field of the carbonate rock thermal reservoir of the Wumishan Formation in the Dongli Lake area under surface water reinjection conditions in the next 50 years are predicted. The results show that after 50 years of continuous surface water recharge in Wells DLB-01, DLB-02 and DLB-03, there will be no thermal breakthrough to the surrounding mining Wells. However, the temperature of well DL-71 is close to the set thermal breakthrough limit. Therefore, attention should be paid to the water temperature monitoring of the nearby mining Wells. Surface water reinjection significantly raised the surrounding thermal storage water level. After 50 years of surface water reinjection in Well DL-71, the water level rose and stabilized at about 153.5 meters.
0 引言
地热能是一种清洁、环保的可再生能源，是能够被人类经济开发和利用的地球内部热能资源 （Lin et al.，2021；王文中等，2022）。地热尾水回灌作为地热开发利用中的关键环节，能减缓水位下降、保持热储压力（刘久荣，2003；Kamila et al.， 2021）。有学者（Rivera et al.，2016；Kamila et al.， 2021）提出，实现地热资源可持续开发利用的有效措施之一就是制定科学合理的回灌方案。杨吉龙等（2022）建议积极攻关地热回灌关键技术，加大尾水梯级利用力度，确保应灌尽灌。董寒杰等（2022） 总结天津开发地热的经验，认为“以灌定采”的开发模式值得借鉴。丁正云等（2023）提出，应建立地热井长期动态监测系统等措施保护地热资源。
天津开发利用地热资源始于 20 世纪 70 年代，是中国地热规模化开发利用最早的城市之一。 2021 年，全市地热开采量达到了 5.29×107 m3（殷肖肖等，2024）。由于早期存在大部分只采不灌的单井系统，粗放式的开发导致地下水位快速下降。天津市东丽湖地区拥有丰富的地表水资源，且碳酸盐岩热储易于回灌，可以尝试地表水回灌以对热储进行额外的流体补充。近年来，东丽湖地区在探索地表水回灌方面取得了一定成果，建立了全国规模最大的生产性地表水回灌补给工程，试验最大回灌量达到了 152 m3 /h（沈健，2015；阮传侠等，2017）。刘东林等（2019）、李义曼等（2020）对该地区地表水回灌开展研究，认为外源水回灌至碳酸盐岩热储层可减缓热储压力下降趋势，且不会对储层结构产生破坏。
随着科学技术的发展，数值模拟方法作为一种有效且低成本的手段逐渐被用来研究复杂的地下构造和物理过程（彭展翔，2016；Liu et al.，2017）。刘东林等（2015）使用 TOUGH2 计算出东丽湖地区雾迷山组热储地热资源量，并对该区域地热开发利用提出了合理建议。阮传侠（2018）利用 TOUGH2 对天津地区雾迷山组热储压力场和温度场进行了模拟计算，认为地热回灌能缓解水位下降趋势，保证地热资源的可持续开发利用。宋丹（2020）针对东丽湖地区地表水回灌，通过 TOUGHREACT 拟合了不同回灌水体对储层的影响，认为地表水与地热水按一定比例混合更利于回灌。马峰等（2023）利用 COMSOL 预测了雄安新区容城地热田温度场和压力场的变化趋势，并提出了合理的回灌方案。
以往的研究大多将关注重点放在回灌的近期表现上，地表水长期回灌是否会对储层的温度场和压力场产生影响还需进一步研究。本文利用有限元计算软件 COMSOL 建立了潘庄凸起区渗流数值模型，在地表水回灌条件下，对东丽湖地区雾迷山组热储未来50 a内温度场和压力场变化趋势进行了预测，为该地区地表水回灌工作的开展提供了依据。
1 地热地质特征
本文以天津市东丽湖地区作为研究区，总占地面积约50 km2。为了更好的模拟深层地热流体的运移和开采动态，以地热田的地质断裂边界形成较完整的计算区域，即潘庄凸起为研究区（图1）。
1.1 构造特征
研究区位于潘庄凸起中东部，处于山岭子地热田的东北部。潘庄凸起位于沧县隆起北部，北以大口屯—汉沽断裂为界，南以海河断裂为界，西以天津断裂为界，东以沧东断裂为界。凸起区长轴走向北东，缺失古近系和中生界，新生界下伏中生界— 中新元古界。
1 —中生界；2—寒武系；3—元古宇；4—断裂及推测断裂；5—基岩埋深等值线/m；6—研究区
研究区附近发育的主要断裂有沧东断裂，其总体走向北东—北北东，倾向南东，为北西盘上升的正断层，倾角 30°～60°，顶板埋深 1000～2000 m，是一条隐伏的、控制沧县隆起与黄骅坳陷地质建造的重要断裂（赵苏民等，2007）。
1 —地温梯度等值线/（℃/100 m）；2—研究区
1.2 地温场特征
地质构造格局、水文地质条件等影响着地温场平面分布特征，华北地区地温梯度整体呈高低相间的带状分布（沈健，2015；阮传侠等，2017；王贵玲等，2017）。研究区地温场主要受沧东断裂控制，断裂带附近地温梯度值高（图2），超过5.0℃/100 m，向两侧降低至3.0℃/100 m左右。
2 研究区热储数值模型
2.1 理论方法
（1）一维井筒模型
将地热井简化为一维线单元，考虑地热流体沿井筒轴向的渗流传热过程，井筒内流体与管壁及周围岩体的换热过程均由等效换热系数考虑（赵志宏等，2020；Ma et al.，2022）。一维井筒中不可压缩流体的能量方程为：
(1)
式（1）中：ρf—流体密度（kg/m3）；A —井筒的横截面积（m2 ）；Cp，f —流体的常压热容（J/kg∙K）；T —流体温度（K）；u —轴向流速（m/s）；λ —导热系数（W/ m·K）；fD—达西摩擦因子；di —筒内直径（m）；Qwall— 通过管壁发生的流体-围岩的热交换（W/m）。
达西摩擦因子 fD 与雷诺数 Re、表面粗糙度 e 与管内直径di的比值有关：
(2)
(3)
式（2）、（3）中：μ—流体动力黏度（Pa∙s）。
对于层流（Re<2000），fD 与井筒的表面粗糙度 e 无关，可由Stokes公式计算：
(4)
对于湍流（Re>2000），fD可用Haaland方程计算：
(5)
定义等效传热系数来近似描述地热井内流体与周围岩石之间的传热过程，则通过管壁发生的流体-围岩的热交换表示为：
(6)
式（6）中：Text —井筒-岩体接触面温度（K）。
考虑单位长度内从流体通过井筒的热流相等，可推导出等效传热系数为：
(7)
式（7）中：r0—套管内径（m）；r1—套管外径（m）； r2—水泥砂浆外径（m）；λ1—套管导热系数（W/m∙ K）；λ2—水泥砂浆导热系数（W/m∙K）。
流体在井筒内流动产生的对流换热系数 hint可以通过下式计算：
(8)
式（8）中：Nu—努塞尔数
对于层流，Nu为常数3.66。当井筒内部的流体是湍流状态时（2000<Re<6×106，0.5<Pr<2000），努塞尔数可由下式确定：
(9)
式（9）中：Pr—普朗特数
(10)
（2）盖层传热过程
可用热传导方程来描述盖层传热过程，可表示为（不考虑外部热源）：
(11)
式（11）中：ρr —岩体密度（kg/m3）；Cp，r —岩体比热容（J/kg∙K）；λr —岩体导热系数（W/m∙K）。
（3）储层渗流传热过程
对于深部热储，假设含水层完全饱和，用 Darcy 定律描述饱和压力流：
(12)
式（12）中：u—达西流速（m/s）；k—渗透率（m2）； Pf —流体压力（Pa）；z—垂向坐标（m）。
流体连续性方程可以表示为：
(13)
式（13）中：ϕ—多孔介质的孔隙率；Qm—质量源项（kg/m3 ∙s）。
储层渗流传热过程，可用多孔介质传热方程来描述：
(14)
式（14）中：（ ρCp）eff 表示常压有效体积热容，对于多孔介质而言，可以表示为：
(15)
式（15）中：θs表示固体材料的体积分数，与液体体积分数θf（或孔隙率）的关系为：
(16)
多孔介质的有效导热系数 λeff与固体材料的导热系数λs和流体材料的导热系数λf有关，可以表示为λs和λf的加权算术平均值：
(17)
2.2 回灌渗流模型建立
研究区域位于潘庄凸起构造区，东西约为 41 km，南北约为 43 km，模型垂向上构建至深约 4 km 的雾迷山组地层。根据岩层物理性质的差异，地层自上而下依次为：第四系（Q）、新近系明化镇组 （Nm）、新近系馆陶组（Ng）、奥陶系（O）、寒武系（∈）、青白口系（Qb）和蓟县系雾迷山组（Jxw），几何模型如图3所示。
图3研究区几何模型图
2.3 渗流参数及热工参数选择
根据相关回灌试验数据及收集到的针对本研究区域及邻近区域的文献资料等所记录的地层参数数据，结合抽水试验结果及热突破预期，得出研究区域地层的相关参数（表1）。结合文献资料中的记录和模拟的具体情况，在合理的范围内对部分参数进行了选择及简化，将盖层设置为低渗地层。
表1各地层参数取值
2.4 雾迷山组储层地热井
区内共有雾迷山组地热井 27 眼。DLB-01、 DLB-02、DLB-03三眼井为地表水回灌井，仅用于在丰水期6—9月进行回灌；其余24眼地热井（13眼开采井和11眼回灌井）用于冬季供暖，自11月15日至次年 3 月 15 日共计 121 d，开采和回灌同时进行，其余时间为水位恢复期，不进行采灌。
2.5 网格剖分
在 COMSOL 中，综合考虑计算精度与计算效率之间的平衡，地热井及周围区域的网格精细化，最大单元尺寸 50 m，最小单元尺寸 5 m。雾迷山组储层网格尺寸为细化，最大单元尺寸 3140 m，最小单元尺寸 393 m。其余网格尺寸为常规，最大单元尺寸6390 m，最小单元尺寸798 m；模型共包含约30万个四面体单元，地热井由一维线单元代表（图4）。
图4网格剖分图
2.6 时间离散
区内地表水回灌井，仅在丰水期回灌；其余地热井只在供暖季运行（11月15日—3月15日），其时间周期函数（按累计开采月数算）如图5所示（0表示停采，1表示开采），时间步长为30 d。
图5时间周期函数
2.7 边界条件和初始条件
由于当前缺乏深部地下水补给和排泄量数据，根据当地地下水埋深变化规律，模型的侧向边界设定为水头边界（其中天津断裂处为隔水边界）。模型顶部为地表，与热储层之间有盖层隔开，故将顶部设为隔水边界。模型底部为雾迷山组地层，认为其与下伏地层之间无水力联系，故将底部设为隔水边界。为更精确描述模型边界处温度场的连续变化，将侧向热量边界设置为温度开边界（其中天津断裂处为热绝缘边界）；顶部为地表，取当地年平均气温（14.1℃）设为温度边界；模型底部以根据温度梯度计算得到的地温设为温度边界。
初始水位条件由天津市蓟县系雾迷山组的水位埋深数据插值建立；在分析当地大地热流资料的基础上，根据提供的地温梯度计算模型初始温度场，利用研究区地温梯度及恒温层的温度，计算各个节点的初始温度。
2.8 水位拟合
2021年9—10月，为验证地热井回灌能力，对工作区地热井（DLB-01、DLB-02、DLB-03）开展了回灌试验，回灌流体为系统处理后的地表水，调整模型相关水文地质参数和边界条件，拟合水位与试验动水位埋深数据基本一致（图6）。
图6回灌试验水位拟合与实测水位对比图
a—DBL-01；b—DBL-02；c—DBL-03
3 结果与讨论
利用渗流模型，通过设定地热井回灌方式、采灌水量、回灌温度等参数，以开采井温度下降2℃作为热突破界限，预测采灌条件下的热突破时间。未来 50 年预计采灌情况见表2，DLB-01、DLB-02 和 DLB-03 井仅在丰水期回灌，其余地热井只在供暖季运行。
表2模拟采灌量
注：三眼井仅在丰水期回灌。
3.1 热突破预测
考虑地热梯级利用，回灌温度设定为 25℃，计算时长为50 a，时步为30 d。经模型模拟计算，开采井在50 a内均不会产生热突破。
以地表水回灌井附近热储情况为例，图7a为地表水回灌运行 50 a时，地表水回灌井周边热储温度分布情况。图7b为未进行地表水回灌，周边采灌系统运行 50 a时，周边热储温度分布情况。从图中可以看出，DLB-01、DLB-02 和 DLB-03 井进行地表水回灌不会对周边开采井产生热突破现象。
以距离三眼回灌井较近的开采井DL-71和DL-76 为例（图8）：DL-71 位于 DLB-01、DLB-02 和 DLB-03 三眼回灌井周围，主要受 DLB-03 影响，水温有所下降，t=50 a时温度约为95.5℃，接近设定的热突破界限，建议地表水回灌站投入运行后加强对附近开采井水温的监测；DL-76 距离 DLB-03 井距离为660 m，受影响较小，温度基本稳定。
图750 a后周边热储温度场
a—回灌条件下；b—未回灌条件下
图8周边井未来50 a温度变化
a—DL-71井；b—DL-76井
3.2 水位预测
以开采 50 a范围内，开采井静水位埋深不超过 250 m 作为临界值；回灌井水位上升至距离井口 10 m 作为临界值。经统计，所有开采井在运行 50 a 后的水位均符合要求。
以地表水回灌井附近热储情况为例（图9），从图中可以看出，地表水回灌条件下，DLB-01、DLB-02和DLB-03井周边热储水位上升明显。
距离三眼回灌井较近的开采井DL-71和DL-76 情况如下（图10）：相比无地表水回灌时水位的持续下降，DL-71 在地表水回灌时，水位上升并稳定在 153.5 m 左右；相比无地表水回灌时水位的迅速下降，DL-76在地表水回灌时，水位虽下降，但下降幅度不超过10 m。
图950 a后周边热储压力场
a—回灌条件下；b—未回灌条件下
图10周边井未来50 a水位变化
a—DL-71井；b—DL-76井
4 结论及建议
（1）基于 COMSOL 软件，建立了研究区渗流模型，并对附近热储层温度场、压力场进行了模拟，结果表明：DLB-01、DLB-02 和 DLB-03 井地表水回灌 50 a内，未对周边开采井产生热突破，但DL-71井温度有所下降，因此需加强对周边地热井水温的监测，及时调整回灌方案；附近热储水位上升明显，周边地热井水位的下降趋势得到了一定的缓解。
（2）地热回灌能减缓水位下降、维持热储压力，是地热资源可持续利用的有效措施之一。天津市拥有较为丰富的地表水资源，可持续推进地表水回灌工作的开展。
致谢 感谢清华大学赵志宏老师对本文提出的宝贵意见。