基于Mod flow和时间序列法的汝州市温泉镇地热可开采量预测

左伟，刘运涛，李建林，王树威

（1.河南省地质矿产勘查开发局 第五地质勘查院，河南 郑州 450001；2.河南理工大学 资源环境学院，河南 焦作 454000）

0 引言

地热资源是一种清洁、低碳的可再生能源。地热资源的开发利用可有效优化能源供给结构，对解决全球化的生态环境问题具有深远影响。按照埋藏深度，可将地热资源划分为浅层地热、水热型地热和干热岩［1］。浅层地热一般指地表以下200 m内所赋存的地热资源。我国地热资源分布广泛、储量丰富，浅层地热开发前景十分广阔［2］。就河南省而言，18个省辖市内皆有可采潜力巨大的中低温地热资源［3］，全省地热可采流体量为9.222 393×108m3／a，可利用热能量为196.46×1012kJ／a，约合标准煤6.703 5×106t／a［4］。目前，河南省开采使用的地热资源主要来自被埋藏于新生界与前新生界热储层中的孔隙水、岩溶水和裂隙水［5］。在开发利用中，还存在勘查研究程度低与开发深度浅、地热循环演化机理还没有理清、地热开发利用模式相对简单等问题。

近年来，随着计算机行业不断发展，数值模拟在地下水领域被广泛利用。Visual Mod flow是目前被国内外学者普遍认可并使用的专业化地下水流系统模拟软件［6］，我国最早在20世纪90年代开始应用［7］，在地下水资源量评价、矿井水涌水量预测、地下水位动态研究等方面得到广泛应用［8］。刁维杰等［9］通过建立地下水数值模型，解决了潍坊市北部地区的水资源配置问题；高学平等［10］为定量分析黎河下游的水量平衡情况，构建了耦合HEC-RAS与Mod flow的河流-地下水系统模型，得到了令人满意的结果；刘基等［11］利用Mod flow对矿井涌水量的模拟与预测进行了深入探讨；李文雅等［12］为评价海勃湾水库的地质环境问题，基于Visual Mod flow对蓄水后的海勃湾水库地下渗流场进行了模拟，取得了良好效果；魏传云等［13］、李巍 等［14］、饶磊等［15］均将Visual Mod flow应用于地下水污染防治，为研究地下水环境评价提供可靠依据。但在地热研究中（特别是在地热水资源储量丰富的河南地区），数值模拟应用研究的成果较少［16］，为探索新的有效途径，本文在分析中引入此法。时间序列法［17］是利用井孔水位降深随时间变化的数据定量计算地下水可开采量、预测地下水位动态变化的实用性理论，在地下水动力学领域中广泛使用，也是分析地下浅层热储可开采量的普适性方法。本文利用Visual Mod flow模拟软件建立地下浅层热储数值模型，并与传统地下水动力学法进行对比，深入探究两种分法的优劣和结合的优势，对河南省汝州市温泉镇浅层地热田的可开采量进行科学预测，以期为该区域的浅层地热田合理开发提供依据。

1 研究区概况及地热循环机理

研究区位于汝州市西部，南临汝河，北望箕山，西傍广成湖（涧山口水库），东为汝河冲积平原，占地面积约8.2 km2。区域内地层自老至新出露为太古界、中元古界、寒武系、石炭系、二叠系、三叠系、新近系、第四系。温泉镇地热田盖层由第四系松散层构成，平均厚度为12 m。地下含水系统包括松散岩类孔隙水、碳酸盐岩溶裂隙溶洞水、碎屑岩类裂隙孔隙水和岩浆岩及变质岩裂隙水4种类型，其中碳酸盐岩溶裂隙水是温泉镇地热水的主要来源，主要赋存于灰岩岩溶裂隙含水层中，该含水层岩溶裂隙发育，富水性强，水量丰富。区内热水井水温均大于30℃小于70℃，属于低温地热田，主要应用于采暖、理疗和洗浴。

温泉镇地热田位于北西向断裂带（F1，F2，F3）与北东向构造（F6，F7）的交汇部位（图1）。地热田水热活动严格受断裂构造（断裂）控制，具有断裂深循环型（对流型）地热系统特征。图2主要反映了该区域深循环地热系统热量和质量传递过程。从地下较远处补给深部热水（Tr）沿通道（深度H）上升，上升热水在一定深度（Hi）和来自透水层中冷水（T0）相遇混合，混合水再沿通道上升（温度并不断散失），溢出地表便形成温泉（Ts）。从区域地质构造条件分析，北西向断裂带为九皋山-温泉街断裂和新安-平顶山断裂在此交汇复合形成，九皋山-温泉街断裂和新安-平顶山均为区域性的控热断裂。东西与北西向（F1F3）主干断裂以导水为主、导热为辅的复合断裂带，主要将西北部碳酸盐岩水输送到温泉镇热储层内。北东向（F6）以导热为主、导水为辅的次级断裂带，控制着温泉的形成、储存与出露。

图1 研究区地热田异常区分布范围Fig.1 Distribution range of geothermal field anomalies in the study area

图2 研究区深循环水热系统示意图Fig.2 Schematic diagram of deep-circulation hydrothermal system in the study area

2 地热可开采量预测

2.1 时间序列法预测

2.1.1 时间序列法

依据《地热资源地质勘查规范》（GB11615-2010）规定，单个地热井，确定开采量使用的压力降低值一般不大于0.3 MPa（约合30.6 m）。而对于多井开采的地热田，主要以地热田内代表性监测井保持一定水位年降速条件下的地热田开采量作为一定时限内的可开采量。首先，通过群孔抽水试验资料，对短期开采状态下的主要井孔水位进行模拟，得出基本适应于本地热田的水文地质参数；其次，在一定时期内推算出在适当降深情况下的抽水量，以此抽水量作为可采地热资源。

2.1.2 预测步骤

（1）数学模型及边界条件概化。根据研究区物探解译成果，温泉镇地热赋存区是沿北西向断裂带（F1—F3）为南部北部边界，F6断层以西区域，F6断层为隔水边界，地下水位埋深4.7～7.7 m。其余地段地下水储水能力差。为便于计算，以F6断层为界，数学模型概化为半无限承压含水层。据此，采用干扰井群非稳定流公式计算允许开采量及降深。

（2）计算公式的选择。根据概化的边界条件，按映射原理以F6直线隔水边界为对称轴进行映射，利用承压完整井干扰非稳定流计算，其计算公式为

式中：s为任意点处的降深，m；Qi为第i个井的出水量，m3／d，包括虚拟井（水量为正）；T为导水系数，m2／d；W（ui）为井函数；u′i为虚拟井自变量；ui为井函数自变量，ri为计算点到第i个井（含虚拟井）的距离，m；S为弹性释水系数；t为抽水延续时间，d。

（3）布井方案的确定。根据温泉镇地热富水带的位置，结合研究区水文地质条件，选择F2以南、F3以北、F6以西、F7以东区域布井（图3），单井出水量分别为70，43，30m3／h，总开采量3432m3／d。虚拟井3口，水量3 432 m3／d。

（4）参数验证。将前期推算出的研究区水文地质参数代入干扰井群公式反演，以2018年9—10月群孔抽水试验资料为依据，抽水量分别为70，43，30 m3／h，虚拟井水量与抽水井相对应，以拟合14 d ZK1孔抽水时降深为主。经模拟反演，当T＝253.07 m2／d，弹性释水系数S＝0.006 92时与抽水资料吻合（计算降深15.81 m，抽水试验降深16.508 m，扣除井损0.707 m，实际降深为15.801 m）。

（5）最大降深与可允许开采量的计算。考虑温泉镇地热属于带状热储，地热资源有限，将反演参数代入干扰井群公式计算，初步按100 d计算，累计降深30 m左右。经计算，开采100 d后，ZK1井为水位下降中心区，水位降深33.395 m（表1）。

由表1可知，为满足长期开采的需要，当群井方案中的抽水量为3 432 m3／d且以100 d为期进行计算时，ZK1井的降深基本稳定不变。故而可得温泉镇地热田的可允许开采量为3 432 m3／d。

2.2 基于Mod flow的地热流体数值模拟

2.2.1 温泉镇地热田水文地质概念模型

模型构建范围为汝州市温泉镇浅低温地热田。地热田严格受北西向断裂（F1—F5）和北东向断裂（F6—F7）控制，其中F2和F3穿越了寒武系灰岩，构成地热水的活动场所，因此沿北西向断裂带（F2—F3）和北东向构造（F6—F7）之间形成了地热异常区（图1）。为准确模拟温泉镇地热田的地下水状态，模拟过程将南北部边界分别向外延伸至F1和F3。因此，模型的南部边界为F1，北部边界为F3，东部边界为F6，西部边界为F7。

模型主要以赋存在灰岩岩溶裂隙含水层中的碳酸盐岩溶裂隙水进行构建。该含水层受构造控制明显，处在断层带上，多个层段岩石破碎或裂隙较发育，水文地质条件较好，水文地质参数随空间变化明显，垂直方向与水平方向渗透系数不同，故概化为非均质、各向异性。含水层上覆第四系冲洪积物，并有砂质黏土岩和黏土层，厚度较薄，为7～13 m，构成相对隔水层，地下水处于承压状态。综上所述，将其概化为非均质各向异性的承压非稳定流模型。

2.2.2 数学模型

根据上述的水文地质概念模型，可用以下数学模型进行描述：

式中：kxx，kyy，kzz为沿x，y，z坐标轴方向的渗透系数，m·d－1；h为点（x，y，z）在t时刻的水位标高，m；W为源汇项，L·d－1；Ss为单位储水系数，L·m－1；t为时间，d；Ω为承压区域；S1为第一类边界；S2为第二类边界；φ（x，y，z，t）为第一类边界上水位标高；q（x，y，z，t）为流量在时间和空间上的变化函数；n为二类边界内法线向量。

2.2.3 研究区空间离散

应用地下水数值模拟软件Visual Mod flow对研究区进行等距或不等距长方体网格剖分，得到空间离散单元，为提高模拟的精度，对观测井、抽水井、边界等区域进行局部加密，最终在平面上剖分成92×92的矩形网格单元，垂向上为一层，即碳酸盐岩溶裂隙含水层，共计8 464个网格单元，其中有效单元格3 182个，无效单元格5 282个。

2.2.4 边界、源汇项处理

南部边界：F1断裂带位于温泉镇南200 m，该断层将寒武系灰岩阻断，缺少热活动必备的空间条件，起隔水作用，故设为阻水边界；北部边界：F3断裂带同样将寒武系灰岩阻断，起隔水作用，设为阻水边界；东部边界：F6断层切断了所有北西向断裂，使寒武系在断裂以东消失，阻断了热水沿F3断裂向东南的通道，设为阻水边界；西部边界：F7断裂带接受地下水补给，设为一般水位边界（GHB）［18］。

本区地下水流的补给来源为断裂带降水入渗补给和侧向径流补给，以大气降水的垂向补给为主。大气降水经灰岩岩溶裂隙、断破碎带渗入地下深部，汇入F2，F3断裂带，入渗系数取0.5。通过Visual Mod flow中的RECHARGE程序包，可计算降水入渗补给量。侧向径流补给量主要在西部边界，将其概化为补给边界。地下水排泄以人工钻井排泄为主要途径，排泄项采用抽水井的方式表示。

2.2.5 水文地质参数分区

结合温泉镇的水文地质资料及前人研究成果，将研究区划分为5个参数分区（图4），并赋予储水系数Ss和渗透系数k初值，作为模型计算的基础。

图4 参数分区图Fig.4 Parameter partition map

2.2.6 模型的识别与验证

根据群孔抽水试验资料，选取2018年9月2日至2018年9月16日作为模型的校正、识别阶段，整个模拟期共分3个制度期，每个制度期又分10个计算时间步长。此次群孔抽水试验，共设有抽水井5个，分别为ZK1号、3号（中井）、4号（东井）、疗2和疗3。

通过运行该模型，根据观测井预测水位值与实测值的拟合情况，不断地试算参数，直至两者拟合较好，得到各分区的参数值。经校正、识别，各分区的参数识别结果参见表2。

表2 各参数分区参数识别结果Tab.2 Parameter identification results of each parameter partition

各水位观测井拟合精度如图5所示。

由图5可以看到，在抽水试验过程中，各水位观测孔均分布在95%的置信区间内，达到了模型识别和验证的要求。模拟结果表明，水文地质概念模型基本合理，试算结束后的水文地质参数和各源汇项与实际情况基本相符，可以认为本次所建立的数值模型基本反映了模拟区的地热流体运动规律。可将之应用于地热流体流场和可开采量动态预测等相关研究。

图5 地下水位观测井拟合精度Fig.5 Groundwater level fitting accuracy of the observation well

2.2.7 研究区可开采量预测

采用数值模拟方法确定地下水开采量，应在地下水总补给和总排泄量均衡的前提下，以使地下水得到充分的开采利用，水位不出现持续下降为标准。通过调节抽水井抽水强度，使水位只在某一范围变动或趋于稳定，没有持续下降的趋势，将这种条件下的开采量作为最大可开采量。综上所述，在长期开采条件下，将各水位观测井水位逐渐趋于稳定时的模拟开采量作为温泉镇地热田的最大可开采量。以此为约束条件，对研究区未来5年（2018年9月17日至2023年9月17日）的地下水最大可开采量和地下水位动态变化进行模拟预测（图6）。

图6 模拟期各水位观测井水位变化曲线Fig.6 Water level change curves of each water level observation well in the simulation period

模拟结果表明，当ZK1井开采量为1 200 m3／d，1号井开采量为480 m3／d，3号井开采量为720 m3／d，总开采量为2 400 m3／d时，各观测井水位变化逐渐趋于稳定，各水位观测井水位变化情况如图6所示。因此，确定温泉镇地热田地热流体的最大可开采量约为2 400 m3／d。

2.3 讨论

赋存在地下浅层热储中的地热流体受断层边界、源汇项等复杂地质环境的影响，其运移循环机理是多种复杂变量相互叠加作用后的综合显现，具有特定的水文地质物理机制。时间序列法是Thesis公式在地下水问题研究中的经典应用，也可作为研究浅层地热水问题的一种方法。它将抽象复杂的动态物理过程转变成具有概化性的数学模型，通过简洁的数学方法对地热开采后的水位作出评价；Visual Mod flow数值模拟法就是在严格的物理成因基础上，将研究区的各类限制条件输入到所构建的模型中，以此对地下热储层进行研究分析。前者利用计算机模拟地热动态变化特征，更加生动、详实，且以模型精度为保证，具有更高的可信度。但需要大量而准确的野外勘查数据为支撑，适用于后期室内阶段的分析研究；后者将复杂的地下水流问题转变为数学建模理论，仅从井孔水位降深随时间变化关系的角度出发进行探究，减少了对水文地质资料的依赖程度，为预测水位变化提供便捷、有效的途径。但真实的地热流体处于受各类外界变量影响的非线性地质环境系统中，此种方法考虑因素单一，不能完全反映地热水流的运动状态，只能对可开采量进行简单粗略地估计，更适合开采现场的试验观测。本文将室内解析反演法与室内数值模拟法相结合，通过两种方法的互补优势，为地热资源的实际生产开发提供有效的参考依据。

近年来，随着计算机信息技术在地热资源领域中的不断应用，相比于利用传统时间-降深序列的可开采量预测，涵盖多角度地质物理特性的数值模拟法对浅层热储流体的模拟将会得到长足发展。

3 结论

（1）汝州市温泉镇地热田为浅层低温地热田，地热类型属断裂深循环型地热系统。热储层由石炭系、寒武系灰岩破碎带组成，地下热水为裂隙承压水。

（2）采用时间序列法得到研究区的可允许开采量为3 432 m3／d；利用Visual Mod flow数值模拟法预测未来5年水位降深达到稳定状态后，地热流体的最大可开采量为2 400 m3／d。结合研究区水文地质条件，汝州市温泉镇浅层地热田的地热流体的最大可开采量应处于2 400～3 432 m3／d之间。浅层地热中应用Visual Mod flow数值模拟法，实现了现场试验技术与室内研究方法的有机结合，可为地热可开采量的确定提供有效支撑。