砂岩热储温度场对回灌参数的响应机理与规律

Response Process and Mechanism of Sandstone Geothermal Reservoir Temperature to Reinjection Parameters

摘要

研究地热尾水回灌引起的热储层内温度场演化对地热资源的可持续利用具有重要意义.本文通过大型砂槽仿真试验模型的渗透试验、示踪试验和回灌试验结合数值模拟方法，研究了回灌参数与开采井热突破时间的定量关系，并通过非线性拟合和参数敏感性分析讨论了流体粘度与密度对高温流体向低温流体回灌结果的影响，以及回灌参数对开采井热突破时间的影响程度和内在机理与规律.结果显示回灌水在不同渗透率的砂岩层内运移速率不同，开采井热突破时间t分别与Q ^‒0.85、ΔT ^‒0.21和R ^1.4呈线性关系.相关方程和分析结果表明，采灌温差ΔT在大于30 ℃时，其变化对开采井热突破时间t的影响已变得微弱，这是由于ΔT通过影响18.5 ℃等温线在温度过渡区内的相对位置来对开采井热突破时间t产生作用，而试验中采取的高温流体向低温流体回灌产生的误差可以引入粘度修正系数α_μ 修正.

关键词

回灌参数 / 仿真试验 / 数值模拟 / 温度场 / 演化机理 / 地热能 / 地热井.

Key words

reinjection parameters / simulation test / numerical simulation / geo-temperature / mechanism of evolution / geothermal energy / geothermal wells

0 引言

地热资源是能够经济地被人类所利用的地球内部的热能，是一种储量大、分布广、绿色低碳的可再生清洁能源（王伟等，2021；徐梓矿等，2021；康凤新等，2023a）.近年来常规能源日趋紧张，环境问题日益严重，地热逐渐受到人们的重视，掀起一轮地热开发高潮.但是深层地热水的大规模开采会造成热储压力下降，导致地热资源利用率降低，同时尾水的直接排放也会造成水资源的严重浪费和周边环境的热污染以及水化学污染（赵季初，2013；吴立进等，2016；康凤新，2018）.地热尾水回灌对于提高地热资源利用率、减少尾水排放和维持热储压力，实现地热资源的可持续利用具有重要的作用（何满潮等，2003；朱家玲等，2012；康凤新等，2023b）.

由于回灌的地热尾水温度远低于储层温度，回灌必然会造成储层温度的降低.当回灌尾水的冷锋面运移到开采井时，会造成开采井的热突破.随着开采井水温的不断降低，地热资源的利用率也随之迅速降低，直至失去开采价值（Seibt and Kellner， 2003； Obembe et al.， 2016； Kang et al.， 2021）.因此，研究地热尾水回灌引起的储层内温度场变化具有重要意义.为此，雷海燕和朱家玲（2010）运用数学模型结合数值模拟方法，研究了回灌量对热储层温度和压力场的影响；刘志涛等（2019）通过对全井段温度、温度梯度、增温率的监测，定性分析并定量计算了砂岩热储地热回灌井热储温度恢复的来源，发现主要热量来源为外围同层相对高温地层的热传导以及地热水流动；Bodvarsson （1972）通过理论分析的方法计算了均质各向同性多孔介质中冷水回灌造成的热污染范围，并给出了热污染范围与含水层厚度、回灌量和储层热容的关系式；Ganguly et al.（2014）提出了一种描述注入冷水时地热储层瞬态温度分布的解析解，并研究了不同体积注入速率、热导率和多孔介质孔隙率等参数对瞬态热传递现象的影响，结果表明注入速率和热导率对瞬态温度分布有重要影响；Saeid et al.（2014， 2015）通过试验结合数值方法研究了饱和多孔介质中流体密度和粘度对热流的影响，结果表明考虑地层水密度和粘度随温度的变化对预测热突破时间具有重要影响；Mottaghy et al.（2011）和Sippel et al.（2013）运用数值模拟方法建立了海牙西南区域的热储模型，并研究了储层结构和孔隙度、渗透率对热储开采井热突破的影响；Sauty et al.（1980）、Bedre and Anderson（2012）、Franco and Vaccaro （2014）定性地研究了流量、井距、回灌温度和储层几何形状对地热系统寿命的重要性，成果仅适用于进行初步设计，供地热工程师在项目的早期阶段使用.前人对影响开采井热突破时间的各项因素进行了广泛的研究，但是未通过大型模型试验和针对地热采灌系统的可人为调节参数进行全面深入的研究分析.

本文通过大型砂槽仿真模型试验，研究不同回灌温度引起的热储内部温度场演化规律.根据砂槽仿真模型渗透试验、回灌试验和示踪试验结果，建立了砂槽仿真模型的数值模拟模型并标定了数值模型的水动力学和热物性参数，并运用建立的等效砂槽数值模型，研究回灌温度、回灌量和采灌井间距3个回灌参数与开采井热突破时间之间的定量关系.

1 研究方法

1.1　砂槽仿真模型试验研究

1.1.1　砂槽仿真模型设计

（1）砂槽仿真模型主体结构.砂槽仿真模型长12 m、宽6 m、高6 m，热储介质采用分选及磨圆较好的石英砂岩，清洗过后采用振动平板分层回填，每30 cm为一层，均匀振动，回填期间采集砂岩样品进行密实度测试.底部4~6 m为粗砂岩，2.8~4.0 m为中砂岩，上部回填厚度为2.8 m的黏土作为隔水层.砂槽的砂岩层位于地表以下，黏土层位于地表以上，即试验层段位于地表以下，上部隔水层位于地表以上，具体如图1所示.砂槽仿真模型底面和侧面边界为钢筋混凝土结构，内壁设置防水和保温层，起到砂槽整体的隔水和隔热作用.

（2）砂槽井孔设计.砂槽仿真模型布设有1眼开采井、1眼回灌井、2眼备用回灌井及13个水位监测孔.开采井与回灌井分别位于砂槽横向两端，距离砂槽横向侧边均为1 m，开采井和回灌井均位于砂槽横向中轴，距横向两侧边均为3 m，2眼备用回灌井沿横向中轴回灌井对称分布，距离横向长边均为1.5 m.水位观测孔沿砂槽横向和竖向中轴布置，间隔均为1 m，具体如图2所示，具体井管尺寸参数如表1所示.

（3）砂槽监测点设计.在开采井、回灌井和水位监测孔内设置水位监测探头，采集实时水位数据.开采井和回灌井内安装电磁流量计，采集实时开采量与回灌量数据.砂槽仿真模型内部温度与压力监测点布置如图3和图4所示.垂直于X轴方向布设共5个监测面，间隔3 m，编号a→e；垂直于Y轴方向布设共4个监测面，间隔2 m，编号A→D；垂直于Z轴方向布设共4个监测面，间隔1 m，编号1→4，压力和温度监测点位置为所有面的交点，进行实时温度与压力数据采集.水样取样点位置与砂槽仿真模型内部温度与压力监测点位置相同，数量减半、间隔布置.

1.1.2　试验过程

采用小流量进行砂槽注水，注入速率缓慢且稳定，注水的同时密切监测观测孔水位、开采井与回灌井流量计、水位传感器以及砂槽仿真模型内部压力、温度传感器的数据变化，进行监测设备的调试.直至地热水充满整个热储介质，各测压管内水面保持齐平，在各层压力、温度稳定后开展模型试验.此时砂槽仿真模型的初始水头为5.2 m（以砂槽底部为基准），初始温度15 ℃.

（1）渗透试验.在试验前期准备工作完成后，首先进行砂槽仿真模型的渗透试验.试验分为大降深、中降深和小降深3个阶段，开采井的抽水量与回灌井的回灌量相同，采灌过程同步进行.三个降深采用的流量分别为1.10 m³/h、0.72 m³/h、0.52 m³/h，水位达到稳定需持续8 h，三组试验共历时6 d.

（2）回灌试验.在抽水试验结束，各测压管内水面再次齐平后，开展砂槽仿真模型不同回灌温度的回灌试验.由于砂槽模型的外界环境温度较低，将砂槽仿真模型整体加温后模型边界会有大量热量散失，因此保持砂槽内多孔介质和流体的温度接近外界环境温度，采取高温水进行回灌，减少热量散失.此时砂槽仿真模型初始温度15 ℃.首先进行回灌温度30 ℃的回灌试验，开采井的抽水量与回灌井的回灌量相同，固定为 0.3 m³/h，回灌温度固定为30 ℃，采灌过程同步进行，共历时14 d.30 ℃回灌试验结束后，待砂槽整体温度降为15 ℃，各测压管内水面再次齐平，开始回灌温度50 ℃的回灌试验.将回灌温度调整为50 ℃，采灌流量调整为0.8 m³/h，具体试验过程与30 ℃回灌试验相同，50 ℃回灌试验共历时10 d.

（3）示踪试验.待回灌试验结束，砂槽仿真试验模型水位和水温恢复到初始值后，开展示踪试验，示踪试验的示踪剂选择钼酸铵.开采井的抽水量与回灌井的回灌量固定为0.4 m³/h，采灌过程同步进行，待砂槽仿真试验模型的水位动态稳定后，将 500 g示踪剂钼酸铵，溶入3 L水中，于回灌井中瞬时投放.投放后取水样检测的频率为4 h/次，示踪剂到达监测点后加密为2 h/次，示踪试验共历时6 d.

1.2　数值模拟方法

按照砂槽仿真模型的实际结构建立砂槽的数值模拟模型，根据实际砂槽模型的回灌、示踪等试验数据拟合数值模型的水动力学和热物性参数.运用建立的等效砂槽数值模型，研究砂岩热储温度场对回灌参数的响应过程与机制.

1.2.1　等效砂槽数值模型建立

运用COMSOL Multiphysics数值模拟软件建立砂槽水‒热‒化三场耦合的等效数值模拟模型，模拟砂槽仿真模型在对井采灌条件下储层温度、地下水位和示踪剂浓度的时空变化.COMSOL Multiphysics能够实现有限元分析中的多物理场耦合，被广泛应用于各类工业生产和科学研究中.

（1）多孔介质渗流.对于热储含水层，假设含水层完全饱和，饱和压力流可以用Darcy定律来描述：

u = - \frac{K}{{\rho }_{\mathrm{f}}\mathrm{g}} \left(\nabla p+{\rho }_{\mathrm{f}}\mathrm{g}\nabla z\right)，（1）

其中，u表示Darcy流速（m/s）；K表示热储层的渗透系数（m/s）；P表示流体压力（Pa）；z表示纵坐标（m），ρ _f表示液体的密度.

流体连续性方程可以表示为：

\frac{\partial }{\partial t} \left({\rho }_{\mathrm{f}}\varphi \right) + \nabla \cdot \left({\rho }_{\mathrm{f}}u\right) = Q_{\mathrm{m}}，（2）

其中，ϕ表示多孔介质的孔隙率；Q _m表示质量源项（kg/（m³∙s））.

（2）多孔介质传热.热储层中的渗流传热过程，可以用多孔介质传热方程来描述：

( \rho C_{\mathrm{p}} {)}_{\mathrm{e}\mathrm{f}\mathrm{f}} \frac{\partial T}{\partial t} + {\rho }_{\mathrm{f}} C_{\mathrm{p},\mathrm{f}} u \cdot \nabla T = \nabla \cdot ( k_{\mathrm{e}\mathrm{f}\mathrm{f}} \nabla T )，（3）

其中，（ρC _p）_eff表示常压有效体积热容，对于多孔介质而言，可以表示为：

{\left(\rho C_{\mathrm{p}}\right)}_{\mathrm{e}\mathrm{f}\mathrm{f}} = {\theta }_{\mathrm{s}} {\rho }_{\mathrm{s}} C_{\mathrm{p},\mathrm{s}} + {\theta }_{\mathrm{f}} {\rho }_{\mathrm{f}} C_{\mathrm{p},\mathrm{f}}，（4）

其中，θ _s表示固体材料的体积分数，θ _f表示液体体积分数.

多孔介质的有效导热系数k _eff与固体材料的导热系数k _s以及流体材料的导热系数k _f有关，可以表示为k _s和k _f的加权算术平均值：

k_{\mathrm{e}\mathrm{f}\mathrm{f}} = {\theta }_{\mathrm{s}} k_{\mathrm{s}} + {\theta }_{\mathrm{f}} k_{\mathrm{f}} .（5）

数值模型参数根据实际砂槽仿真模型的几何特性、含水层介质的垂向分布以及采灌井位置、参数设置.水动力场、温度场和化学场边界条件分别为隔水、绝热与无通量边界，初始值分别为初始水头5.2 m（以砂槽仿真模型底板为基准面），初始温度15 ℃，初始浓度0.根据回灌试验的采灌井水位、流量关系和示踪试验钼酸铵的浓度变化规律，对砂槽仿真模型的渗透率进行反演.根据回灌试验中监测点的温度变化趋势，反演砂槽模型的多孔介质导热系数.

1.2.2　回灌参数对热储温度场的影响

以开采井的热突破时间表征热储温度场的变化，将开采温度变化2 ℃设定为发生热突破，研究回灌温度、回灌量和采灌井间距3个回灌参数对开采井热突破时间的影响，并进行参数敏感性分析.

（1）回灌参数与开采井水温变化关系.回灌参数的取值如表2所示，工况一、二、三分别对应回灌温度、回灌量和采灌井间距3个参数对开采井水温变化的影响，分别计算不同回灌参数下开采井的温度随时间变化情况.

（2）回灌温度、回灌量与开采井热突破时间关系.将采灌井流量依次设置为0.1 m³/h、0.2 m³/h、0.3 m³/h、0.4 m³/h、0.5 m³/h、0.75 m³/h、1.0 m³/h、1.25 m³/h、1.5 m³/h、2.0 m³/h、

2.5 m³/h、3.0 m³/h，计算回灌温度分别为25 ℃、30 ℃、…、60 ℃时开采井发生热突破的时间，拟合回灌温度、回灌量与热突破时间关系.

（3）采灌井间距、回灌量与开采井热突破时间关系.按照实际采灌系统运行时的回灌温度，采灌温差通常在30 ℃左右，据此将回灌温度设定为45 ℃，进行采灌井间距、回灌量与热突破时间的关系研究.将采灌井间距分别设置为3 m、4 m、…、 10 m，计算不同回灌量下开采井发生热突破所需时间，拟合采灌井间距、回灌量与热突破时间关系.

（4）参数敏感性分析.回灌温度、回灌量和采灌井间距3个回灌参数的不同取值对开采井热突破时间的影响各异.因此对3个回灌参数进行敏感性分析，研究3个回灌参数对开采井热突破时间的影响程度，回灌参数敏感性研究取值如表2所示.

2 结果与讨论

2.1　试验结果

2.1.1　渗透试验结果

经过3组不同采灌井水位差的砂槽仿真模型回灌试验，得到不同采灌井流量对应的采灌井稳定水位差，如表3所示.

2.1.2　回灌试验结果

回灌流体温度为30 ℃时，选取砂槽仿真模型Z方向1和3截面（分别为图3绿色与红色监测面）的所有监测点，绘制回灌第5天和第10天的温度云图（图5）.回灌井中注入的高温流体逐渐向开采井运移，沿X轴各监测点的温度陆续升高.第10天时，Z3截面显示回灌高温流体的影响范围已到达开采井.由于砂岩层渗透率的差异，Z1截面的高温流体运移速度明显小于Z3截面.

回灌流体温度为50 ℃时，同样绘制Z1、Z3截面所有监测点第5天和第9天的温度云图（图6）.回灌井中注入的高温流体逐渐向开采井运移，但由于回灌温度和回灌量的增加，流体向开采井的运移速度和监测点的温度增加幅度显著上升.在第5天时，Z1和Z3截面的高温流体影响范围都已到达开采井；第10天时，回灌高温流体基本已填充整个Z3截面所在的粗砂岩层.由于砂岩层渗透率的差异，Z1截面的高温区域范围仍然小于Z3截面.

2.1.3　示踪试验结果

绘制砂槽仿真模型内部Cc1、Cc3两个水样取样点（具体位置见图4）和开采井出水口的钼酸铵浓度变化图（图7）.Cc1点第0.7天开始检测到钼酸铵浓度升高，第2.5天钼酸铵浓度到达峰值；Cc3点第0.5天开始检测到钼酸铵浓度升高，第2天钼酸铵浓度到达峰值，到达峰值后浓度逐渐降低；开采井第2.5天开始检测到钼酸铵浓度升高，钼酸铵浓度分别在第4天和第4.6天达到峰值.由于砂岩层渗透率的差异，Cc3点示踪剂到达和达到峰值的时间要快于Cc1点，峰值浓度也明显高于Cc1点.

2.2　数值模拟结果

2.2.1　等效砂槽数值模型参数

建立等效数值模型时，水文地质参数选取非常重要，渗透率就是其中之一.渗透率的大小取决于含水层介质.由于开采井与回灌井外侧1 m处为隔水边界，承压水井的Dupuit公式与Theis公式均不适用，无法据此求得仿真试验场的各项水文地质参数（如渗透系数、影响半径、单位涌水量等），因此利用渗透试验、回灌试验与示踪试验的监测数据，对砂槽模型的水文地质参数和热物性参数进行拟合，拟合效果如图8所示，拟合结果如表4所示.

2.2.2　回灌参数与开采井水温变化关系

为研究不同回灌参数对开采井水温变化的影响，运用建立的等效砂槽数值模型，分别计算了不同回灌温度、回灌量和采灌井间距条件下开采井的出水温度变化，并绘制了相应的关系曲线分别如图9、图10和图11所示，其中采灌温差ΔT为回灌温度减储层初始温度.

不同采灌温差ΔT时开采井的出水温度随时间的变化情况如图9所示.不同采灌温差ΔT时开采井出水温度的变化趋势基本相同，在回灌的高温地热水到达开采井后开采水温迅速升高，而后升温速度逐渐放缓并逐渐趋近于回灌水温度，开采井的热突破时间t随着采灌温差ΔT的升高而减小.不同回灌量Q时开采井的出水温度随时间的变化情况如图10所示.开采井出水温度的变化趋势与图9相同，开采井的热突破时间t也随着回灌量Q的升高而减小.不同采灌井间距R时开采井的出水温度随时间的变化情况如图11所示，开采井出水温度的变化趋势依旧与图9相似，开采井的热突破时间t随着采灌井间距R的增大而增大.

2.2.3　回灌量、采灌温差与热突破时间关系

为研究回灌温度、回灌量对热突破时间的影响，运用建立的等效砂槽数值模型，分别计算了不同回灌温度、回灌量条件下开采井的热突破时间，绘制回灌量与热突破时间的关系（图12）和采灌温差ΔT与热突破时间的关系曲线（图13）.随着回灌量Q和采灌温差ΔT的增大，开采井发生热突破所需时间逐渐减小，且随回灌量Q增大而减小的速度更快.

2.2.4　采灌井间距、回灌量与热突破时间关系

根据实际工程的平均回灌温度，将采灌温差ΔT设置为30 ℃，运用等效砂槽数值模型，分别计算了不同回灌量和采灌井间距条件下开采井的热突破时间，并绘制采灌井间距、回灌量与热突破时间的关系曲线（图14和图15）.随着采灌井间距R的增大，开采井发生热突破所需时间迅速延长.

3 参数敏感性分析

为确定砂岩热储温度场对回灌参数的响应过程与机制，估算地热对井采灌系统的使用寿命，进行了回灌温度、回灌量和采灌井间距3个回灌参数的敏感性分析，分别确定回灌参数对开采井热突破时间的影响程度.

3.1　采灌温差ΔT

回灌的流体温度代表了返回储层的热量.当采灌井间距R=10 m时，计算了12种不同回灌量条件下，采灌温差ΔT分别为10 ℃、15 ℃、…、45 ℃时开采井的温度变化模型，如图13所示.随着回灌温度的升高，开采井的热突破时间提前，但随着回灌温度的继续升高，其对开采井热突破时间的影响迅速减小.当回灌温度以30 ℃为中心上下变动5 ℃时，开采井热突破之前的温度变化如图16所示，开采井水温开始出现降低的时间基本不随回灌温度的改变而变化，该回灌温度对开采井热突破时间的影响程度在6.9%之内.选取采灌井间距R=10 m，Q=1.0 m³/h和2 m³/h时的采灌温差与热突破时间关系进行非线性拟合，结果如图17和式（6）、（7）所示，开采井热突破时间t与ΔT ^-0.21呈线性关系，t与ΔT ^-0.21的线性关系系数随Q的增大而减小.

t_{\mathit{Q}\mathrm{0.5}} = \mathrm{9.36} ∆ {\mathit{T}}^{-\mathrm{0.21}}，（6）

t_{Q\mathrm{2.0}} = \mathrm{5.93} ∆ T^{-\mathrm{0.21}}，（7）

其中，t为开采井的热突破时间，单位为d，下标Q0.5和Q2.0表示流量分别为0.5 m³/h和2 m³/h；ΔT为采灌温差，单位为℃.

3.2　采灌井间距R

将对井采灌系统中的采灌井放置在彼此之间的最佳距离处可以大幅提高储层的利用率.采灌井间距增加可以提供更大的可开采储层范围，因此使用寿命更长.然而，随着采灌井间距的增大，系统的采收率会降低.

为了定量分析采灌井间距R对开采井热突破时间t的影响程度，研究了采灌温差ΔT=30 ℃时，12种不同回灌量条件下，采灌井间距R分别为3 m、 4 m、…、10 m时开采井的温度变化模型，如图14所示.随着采灌井间距R的增加，开采井的热突破时间t延后，热突破时间t的增加速率也随采灌井间距R的增大有所增加.当回灌量Q=1 m³/h时，采灌井间距R以9 m为中心上下变动1 m，开采井热突破之前的温度变化如图18所示，开采井水温开始出现降低的时间随采灌井间距R的改变而变化，该采灌井间距对开采井热突破时间的影响程度达到20.1%，远大于采灌温差ΔT的影响.选取Q=0.5 m³/h和 1 m³/h时采灌井间距R与热突破时间t进行拟合，结果如图19和式（8）、（9）所示，开采井的热突破时间t与采灌井间距R ^1.4呈线性关系，线性关系系数随Q的增大而降低，与式（6）、（7）的系数变化规律相同.

t_{\mathit{Q}\mathrm{0.5}} = \mathrm{0.36} R^{\mathrm{1.4}}，（8）

t_{\mathit{Q}\mathrm{1.0}} = \mathrm{0.21} R^{\mathrm{1.4}}，（9）

其中，t为开采井的热突破时间，单位为d，下标Q0.5和Q1.0表示流量分别为0.5 m³/h和1 m³/h；R为采灌井间距，单位为m.

3.3　回灌量Q

回灌量是地热对井采灌系统运行后可以调节的参数，它可以对系统寿命产生直接影响，是影响系统寿命最重要的因素之一.为了定量分析回灌量的参数敏感性，研究了0.1 m³/h、0.2 m³/h、0.3 m³/h、0.4 m³/h、0.5 m³/h、0.75 m³/h、1.0 m³/h、1.25 m³/h、1.5 m³/h、2.0 m³/h、2.5 m³/h、3.0 m³/h共12种情况，如图12和图15所示.随着回灌量Q的升高，开采井的热突破时间迅速减小，对比图13结果，回灌量Q对开采井的热突破时间的影响明显大于采灌温差ΔT.当采灌井间距R=10 m，采灌温差ΔT=30 ℃时，回灌量Q以1.0 m³/h为中心上下变动0.25 m³/h，开采井热突破之前的温度变化如图20所示，开采井水温开始出现降低的时间随回灌量Q的改变而产生较大变化，该回灌量的上下波动对开采井热突破时间的影响程度达到36.6%，远大于采灌温差ΔT的影响.选取采灌井间距R=10 m，ΔT=10 ℃和45 ℃时的回灌量与热突破时间关系进行非线性拟合，结果如图21a和式（10）、（11）所示，开采井热突破时间t与Q ^-0.85呈线性关系，且t与Q ^-0.85的线性关系系数随采灌温差ΔT的增大而降低，但降低的幅度较小，符合对ΔT的敏感性分析结果.同样，对采灌温差ΔT=30 ℃，R=3 m和R=10 m时的回灌量与开采井热突破时间进行拟合，结果如图21b和式（12）、（13）所示，当采灌井间距R改变时，开采井热突破时间t依然与Q ^-0.85呈线性关系，证明了t与Q线性关系的合理性.

t_{\mathrm{10}} = \mathrm{5.87} Q^{-\mathrm{0.85}}，（10）

t_{\mathrm{45}} = \mathrm{4.26} Q^{-\mathrm{0.85}}，（11）

t_{\mathit{R}\mathrm{3}} = \mathrm{1.85} Q^{-\mathrm{0.85}}，（12）

t_{\mathit{R}\mathrm{10}} = \mathrm{4.83} Q^{-\mathrm{0.85}}，（13）

其中，t为开采井的热突破时间，单位为d，下标10和45表示采灌温差分别为10 ℃和45 ℃，下标R3和R10表示采灌井间距分别为3 m和10 m；Q为回灌井回灌量，单位为m³/h.

4 讨论

4.1　流体粘度与密度的影响

为了克服砂槽仿真模型难以整体加温和试验过程中边界热量的散失，试验采用了高温回灌来研究开采井水温变化和热突破时间.由于高温流体向低温流体流动与低温流体向高温流体流动的作用是相反的，而流体温度的变化会引起粘度和密度的变化，这个相反的过程对研究结果可能会产生一定的影响.因此取采灌温差ΔT分别为10 ℃、20 ℃、30 ℃和40 ℃时，讨论高温回灌和低温回灌两个相反过程对研究结果的影响，回灌参数取值如表5所示.

将流体密度ρ固定为1 000 kg/m³，流体动力粘度μ设置为随温度变化（式14），研究高温回灌和低温回灌时流体粘度对开采井温度变化和热突破时间的影响，结果分别如图22和图23所示.高温回灌和低温回灌时开采井的温度变化趋势完全相反，相较于高温回灌的温度升高速率，低温回灌时的温度降低速率相对较低.在同一采灌温差ΔT下，低温回灌时开采井达到热突破的时间要大于高温回灌，且随着采灌温差ΔT的增大，低温回灌与高温回灌时开采井达到热突破的时间差也在增大.这是由于随着采灌温差ΔT的增大，砂槽仿真模型内流体的动力粘度和回灌流体的动力粘度差也在不断增大，导致用高温流体向低温流体流动来反向模拟低温流体向高温流体流动过程的误差不断增大.因此引入粘性修正系数α_μ 来表示低温回灌与高温回灌时开采井热突破时间比值，并与不同采灌温差ΔT对应的流体动力粘度差Δμ进行拟合，结果如图24和式（15）所示.低温回灌与高温回灌时开采井达到热突破的时间比α_μ 与砂槽水初始温度对应的动力粘度和回灌水的动力粘度差Δμ呈对数函数关系.

（14）

其中，μ为水的动力粘度，单位为Pa·s；T为水的温度，单位为℃.

{\alpha }_{\mathit{\mu }} = \mathrm{l} \mathrm{n} \left(\mathrm{1.78}\times {\mathrm{10}}^{\mathrm{3}}∆\mu +e\right),（15）

其中，α_μ 为粘性修正系数，即低温回灌与高温回灌时开采井达到热突破的时间比；Δμ为砂槽水初始温度对应的动力粘度和回灌水的动力粘度差，单位为Pa·s.

同样固定流体动力粘度μ等参数，研究流体密度ρ在温度变化条件下，高温回灌和低温回灌对开采井温度变化和热突破时间的影响差异，结果显示流体密度ρ随温度的变化对开采井热突破时间的影响小于动力粘度μ对开采井热突破时间影响程度的5%，与动力粘度μ相比可忽略不计.以上结果表明，采取高温流体向低温流体流动来替代低温流体向高温流体流动会产生一定的误差，这个误差主要来源于流体的粘度随温度变化，可以引入粘度修正系数α_μ 来修正这一误差.

因此保持采灌温差ΔT=30 ℃不变，同时考虑回灌流量Q和采灌井间距R对开采井热突破时间t的影响，拟合结果如图25和式（16）所示.在对井采灌系统设计时，根据设计使用年限可以确定回灌量Q和采灌井间距R的合理搭配，保证开采井在设计使用年限内不会发生热突破.考虑到模型试验和实际工程储层参数之间的差异，在实际使用时，式（16）中的系数应该按照实际热储的参数进行一定的调整.

t = \mathrm{0.19} {\alpha }_{\mu } R^{\mathrm{1.4}} Q^{-\mathrm{0.85}} .（16）

4.2　回灌参数对热储温度场的影响机理

热储层的孔隙度和初始温度等固有参数会对对井采灌系统的运行有较大影响，但是这些参数是人为不可控的.回灌量、回灌温度和采灌井间距是重要的人为控制参数，可以对其进行优化，以获得最长的使用寿命和最高的能量生产效率.敏感性分析结果表明，采灌温差ΔT和回灌量Q的增加会加快开采井的热突破时间t，采灌井间距R的增加会延缓开采井的热突破时间.ΔT在绝对值较小时，其变化对于开采井的热突破时间t影响显著，随着ΔT的增加，其变化对开采井热突破时间的影响程度迅速降低，当ΔT增大到30 ℃左右时，其上下波动1 ℃，对开采井热突破时间变化的影响已经降低到1.4%左右，基本可忽略不计.相较而言，回灌量Q和采灌井间距R对于开采井热突破时间t的影响一直很大，说明不同回灌参数对热储温度场的影响机理不同.

取开采井发生热突破时井筒周边50 mm厚的薄层，研究其微观温度分布如图26所示，可以发现不同采灌温差ΔT时井筒附近微观温度分布基本无变化，都是在18.5 ℃等温线前缘到达开采井壁时，开采井发生热突破.绘制回灌第2.5天，不同采灌温差ΔT时采灌井所在轴线和平面温度分布，如图27所示.储层内部温度场可以划分为回灌高温区、初始低温区和温度过渡区3部分.随着采灌温差ΔT的增大，过渡区的宽度缓慢增长且增长速率越来越小.当ΔT从30 ℃增大到40 ℃，过渡区宽度仅增加了2%，可以认为采灌温差ΔT大于30 ℃时，温度过渡区的宽度已经趋于稳定，当采灌温差ΔT继续增大，只有过渡区内的温度梯度大小随之线性增加.此时温度过渡区内18.5 ℃等温线会随着温度梯度的增大逐渐向过渡区前缘移动且前移的距离逐渐减小，因此随着采灌温差ΔT的增加，18.5 ℃等温线到达开采井井壁的时间逐渐提前，但提前的时间越来越短，符合前节采灌温差与热突破时间的规律关系.

综上所述，由于回灌水温度和含水层初始温度的不同，冷热水交界面处会发生热传导和对流，在回灌高温区和初始低温区之间形成温度过渡区.回灌量Q可以代表回灌水从回灌井向四周扩散的速度，决定了温度过渡区整体运移距离的大小，采灌井间距R决定了温度过渡区向开采井运移的总路程，以上二者决定了温度过渡区到达开采井的时间，而采灌温差ΔT决定的是过渡区的宽度和过渡区内温度梯度大小，ΔT通过影响18.5 ℃等温线在温度过渡区内的相对位置来对开采井热突破时间t产生作用.因此采灌温差ΔT对开采井热突破时间t的影响是在回灌量Q和采灌井间距R的基础上作用的，其影响程度会受到回灌量Q和采灌井间距R的制约，因此采灌温差ΔT的变化对开采井热突破时间t的影响程度随着ΔT绝对值的增大而逐渐减小.在实际对井采灌系统中，采灌温差ΔT普遍都保持在30 ℃以上，因此ΔT的变化对于开采井热突破时间t的影响微弱，与回灌量Q和采灌井间距R相比可忽略不计.因此，可以采取降低回灌温度的措施增大热储的采收率和能量生产效率.

5 结论

本文通过大型砂槽仿真模型试验和数值模拟，研究了回灌引起的热储层内部温度场的演化规律，对回灌温度、回灌量和采灌井间距3个回灌参数与开采井热突破时间的关系进行了定量研究，并讨论了回灌参数对开采井热突破时间的影响程度及其内在机理与规律，得到结论如下.

（1）在回灌井中注入的高温流体逐渐向四周扩散，引起回灌井周边的温度场升高.在热水锋面到达开采井抽水形成的降落漏斗范围后，高温流体向开采井的运移速度加快，热水锋面向开采井延伸.由于仿真模型不同砂岩层渗透率的差异，内部温度场的演化情况有所不同，上部中砂岩层高温流体运移速度明显小于下部粗砂岩层.

（2）回灌的高温地热水到达开采井后开采水温迅速升高，而后升温速度逐渐放缓并逐渐趋近于回灌水温度.开采井热突破时间t随回灌流量Q和采灌温差ΔT的增大而缩短，随采灌井间距R的增大而延长.

（3）采灌温差ΔT在绝对值较小时，其变化对于开采井的热突破时间t影响显著；当其升至30 ℃以上时，其变化对开采井热突破时间t的影响微弱，指示可以通过降低回灌温度的方式提高热储的采收率.与采灌温差ΔT相比，回灌量Q和采灌井间距R的变化对于开采井热突破时间t的影响程度始终较大，指示为防止过早出现热突破，应严格控制回灌量和采灌井间距.开采井热突破时间t分别与Q ^-0.85、ΔT ^-0.21和R ^1.4呈线性关系.

（4）通过高温流体向低温流体流动来反向模拟低温流体向高温流体流动的过程会产生一定的误差，这个误差主要来源于流体的粘度随温度变化，可以引入粘度修正系数α_μ 来修正这一误差，粘度修正系数α_μ 与砂槽水初始温度对应的动力粘度和回灌水的动力粘度差Δμ呈对数关系.

（5）回灌量Q和采灌井间距R决定了温度过渡区到达开采井的时间，采灌温差ΔT决定了温度过渡区的宽度和过渡区内温度梯度大小，ΔT对开采井热突破时间t的影响是在回灌量Q和采灌井间距R的影响基础上作用的.当18.5 ℃等温线到达开采井壁时开采井发生热突破，随着采灌温差ΔT的增加，18.5 ℃等温线逐渐向温度过渡区的前缘逼近，到达开采井井壁的时间因此逐渐提前，但随着ΔT绝对值的增大，ΔT的增加引起的提前时间逐渐缩短.

参考文献

[1] Bedre, M. G., Anderson, B. J., 2012. Sensitivity Analysis of Low-Temperature Geothermal Reservoirs: Effect of Reservoir Parameters on the Direct Use of Geothermal Energy. Transactions-Geothermal Resources Council, 36 2: 1255-1261

[2] Bodvarsson, G., 1972. Thermal Problems in the Siting of Reinjection Wells. Geothermics, 1(2): 63-66. https://doi.org/10.1016/0375-6505(72)90013-2

[3] Franco, A., Vaccaro, M., 2014. Numerical Simulation of Geothermal Reservoirs for the Sustainable Design of Energy Plants: A Review. Renewable and Sustainable Energy Reviews, 30: 987-1002. https://doi.org/10.1016/j.rser.2013.11.041

[4] Ganguly, S., Mohan Kumar, M. S., 2014. Analytical Solutions for Transient Temperature Distribution in a Geothermal Reservoir Due to Cold Water Injection. Hydrogeology Journal, 22(2): 351-369. https://doi.org/10.1007/s10040-013-1048-2

[5] He, M. C., Liu, B., Yao, L. H., et al., 2003. Study on the Theory of Seepage Field for Geothermal Single Well Reinjection. Acta Energiae Solaris Sinica, 24(2): 197-201 (in Chinese with English abstract).

[6] Kang, F. X., 2018. Comprehensive Evaluation of Geothermal Clean Energy in Shandong Province. Science Press, Beijing (in Chinese).

[7] Kang, F. X., Shi, Q. P., Ma, Z. M., et al., 2023a. Genetic Mechanism of the Karst Geothermal Reservoir in Buried Uplifts of Basins: A Case Study of Heze. Acta Geologica Sinica, 97(1): 221-237 (in Chinese with English abstract).

[8] Kang, F. X., Zhao, J. C., Huang, X., et al., 2023b. Heat Accumulation Mechanism and Resources Potential of the Karst Geothermal Reservoir in Liangcun Buried Uplift of Linqing Depression. Earth Science, 48(3): 1080-1092 (in Chinese with English abstract).

[9] Kang, F. X., Zhao, J. C., Tan, Z. R., et al., 2021. Geothermal Power Generation Potential in the Eastern Linqing Depression. Acta Geologica Sinica-English Edition, 95(6): 1870-1881. https://doi.org/10.1111/ 1755-6724.14877

[10] Lei, H. Y., Zhu, J. L., 2010. Modeling of Exploitation and Reinjection of Porous Medium Geothermal Reservoir. Acta Energiae Solaris Sinica, 31(12): 1633-1638 (in Chinese with English abstract).

[11] Liu, Z. T., Liu, S., Song, W. H., 2019. Change Characteristics of Geothermal Field for Geothermal Return Water Reinjection of Sandstone Reservoir in the Northern Shangdong. Acta Geologica Sinica, 93(S1): 149-156 (in Chinese with English abstract).

[12] Mottaghy, D., Pechnig, R., Vogt, C., 2011. The Geothermal Project Den Haag: 3D Numerical Models for Temperature Prediction and Reservoir Simulation. Geothermics, 40(3): 199-210. https://doi.org/10.1016/j.geothermics.2011.07.001

[13] Obembe, A. D., Abu-Khamsin, S. A., Hossain, M. E., 2016. A Review of Modeling Thermal Displacement Processes in Porous Media. Arabian Journal for Science and Engineering, 41(12): 4719-4741. https://doi.org/10.1007/s13369-016-2265-5

[14] Saeid, S., Al-Khoury, R., Nick, H. M., et al., 2014. Experimental-Numerical Study of Heat Flow in Deep Low-Enthalpy Geothermal Conditions. Renewable Energy, 62: 716-730. https://doi.org/10.1016/j.renene.2013.08.037

[15] Saeid, S., Al-Khoury, R., Nick, H. M., et al., 2015. A Prototype Design Model for Deep Low-Enthalpy Hydrothermal Systems. Renewable Energy, 77: 408-422. https://doi.org/10.1016/j.renene.2014.12.018

[16] Sauty, J. P., Gringarten, A. C., Landel, P. A., et al., 1980. Lifetime Optimization of Low Enthalpy Geothermal Doublets. Advances in European Geothermal Research. Springer, Dordrecht, 706-719. https://doi.org/10.1007/978-94-009-9059-3_64

[17] Seibt, P., Kellner, T., 2003. Practical Experience in the Reinjection of Cooled Thermal Waters Back into Sandstone Reservoirs. Geothermics, 32(4-6): 733-741. https://doi.org/10.1016/s0375-6505(03)00071-3

[18] Sippel, J., Fuchs, S., Cacace, M., et al., 2013. Deep 3D Thermal Modelling for the City of Berlin (Germany). Environmental Earth Sciences, 70(8): 3545-3566. https://doi.org/10.1007/s12665-013-2679-2

[19] Wang, W., Fu, H., Xing, L. X., et al., 2021. Crack Propagation Behavior of Carbonatite Geothermal Reservoir Rock Mass Based on Extended Finite Element Method. Earth Science, 46(10): 3509-3519 (in Chinese with English abstract).

[20] Wu, L. J., Zhao, J. C., Li, A. Y., et al., 2016. Key Issues of Geothermal Resource Exploitation and Utilization in the Depression Area of Northern Shandong Province. Geology and Exploration, 52(2): 300-306 (in Chinese with English abstract).

[21] Xu, Z. K., Xu, S. G., Zhang, S. T., 2021. Hydro-Geochemistry of Anning Geothermal Field and Flow Channels Inferring of Upper Geothermal Reservoir. Earth Science, 46(11): 4175-4187 (in Chinese with English abstract).

[22] Zhao, J. C., 2013. Lubei Geothermal Tail Water Reinjection Experiments in Sandstone Reservoir. Shandong Land and Resources, 29(9): 23-30 (in Chinese with English abstract).

[23] Zhu, J. L., Zhu, X. M., Lei, H. Y., 2012. Analysis of Impact of Pressure Compensation between Geothermal Wells on Reinjection Effeciency. Acta Energiae Solaris Sinica, 33(1): 56-62 (in Chinese with English abstract).

基金资助

国家自然科学基金项目(42072331;U1906209)

泰山学者工程专项经费(tstp20230626)