河北平原地热流体可采量计算方法及岩溶热储分布规律研究

张德忠， 马云青， 苏永强

(河北省地矿局第三水文工程地质大队，衡水 053000)

0 引言

河北省地热资源利用主要位于平原区，主要开采热储为砂岩热储，且呈现地热井分布集中和利用形式单一的特点。随着能源供需矛盾不断加剧，地热资源开发利用呈现快速增长态势，主要利用形式有地热供暖、疗养洗浴、种植养殖以及工业加工等。截至2015年底，河北省平原区现有地热井1 120眼，年开采量9 085万m3。总体来说，河北省大部分地区不同程度地开发了地热资源，取得了巨大的经济效益和社会效益，为河北省大气污染治理做出了重要贡献[1]。同时，在开发利用过程中，也存在着地热资源勘查程度低、地热水位下降迅速和地热尾水回灌率低等多方面的问题。《河北省地热能开发利用“十三五”规划》《京津冀协同发展规划纲要》和《河北省关于加快推进生态文明建设的实施意见》等政策的发布和实施增大了河北省对地热资源的需求[2]。本文通过热储法、解析法、统计分析法和数值模拟法4种可采量计算方法的对比分析以及实用性评价，开展了河北省平原区地热流体可开采量评价和对比，进一步探索岩溶热储分布规律，寻找有价值的地热开发目标，为下一步勘查开发提供建议和有力的支撑。

1 研究区概况

河北省地处华北平原北部，兼跨内蒙古高原，环围北京和天津，北与辽宁、内蒙古为邻，西靠山西，南与河南、山东接壤，东临渤海。河北省所处大地构造位置，大致以北纬42°线为界，以北属内蒙—大兴安岭褶皱系南缘，以南为中朝准地台。在中朝准地台上又划分为4个二级构造单元： 内蒙地轴、燕山台褶带、山西断隆和华北断坳，河北省平原区主要位于华北断坳中。河北省地层在太古宇—古元古界变质岩系结晶基底之上，不整合沉积中新元古界长城系、蓟县系、青白口系，古生界寒武系、奥陶系、石炭系和二叠系，中生界三叠系、侏罗系、白垩系，新生界第三系和第四系。

河北省地热资源主要是以地下水为载体的水热型地热资源，多以中低温地下热水或温泉形式分布于平原区及山区(图1)。

图1 河北省地热资源分布图

山区地热资源受断裂构造控制，分布在燕山和太行山地幔隆起带，现发现温泉及地热田45处，如隆化县七家温泉(水温97 ℃)、怀来县后郝窑地热田(水温86 ℃)、平山县温塘地热田(水温73 ℃)； 平原区主要热储层为新近系孔隙热储和古生界、中新元古界岩溶裂隙热储，热源主要来自上地幔和深部花岗岩层中的放射性元素，深部热量以热传导和热对流形式传递到地表。地热资源分布面积约为6万km2，占平原区面积的82%。新近系馆陶组孔隙热储分布范围广，埋藏深度1 000～2 000 m，厚度500～1 600 m，地下热水温度40～85 ℃。岩溶裂隙热储主要分布于构造凸起部位，埋藏深度变化较大，地下热水温度50～104 ℃，具有水量大、水温高、易回灌等特点，代表性地热田有牛驼镇、献县、辛集—宁晋等地热田。在现有经济技术条件下，平原区地热水可开采总量约821.9亿m3，按开采期50 a计算，每年可开采量为16.4亿m3，其热量相当于1 460万t标准煤。

2 河北省地热流体可开采量计算方法

进行地热资源勘查的主要目的之一是进行地热资源储量计算，而储量计算应在建立地热系统概念模型的基础上，根据地热地质条件和研究程度不同选择相应的方法进行。概念模型应能反映热储的空间分布、边界条件和渗透特征，研究地热流体的补给和运移规律。地热资源储量计算的重点是地热流体的可开采量，依据地热地质条件和勘查研究程度，目前平原区层状热储地热流体可开采量计算方法主要有热储法、解析法、统计分析法和数值模拟法。每种方法都有各自的适用条件，应根据勘查阶段、开采时间、监测资料情况选用适宜的计算方法[3]。

2.1 地热流体可开采量计算方法

2.1.1 热储法

河北平原热储为圈闭型热储，地热水是在漫长地质历史内形成的，补给缓慢，虽然热能是可再生的，但地热水可视为不可再生的资源。地热流体储存量计算实质是对计算深度内“静储量”的估算，不考虑侧向补给和越流补给的静态储存量。热储法基本原理是将热储概化为平面上无限延伸的、垂向上为一个主要含水层、其上下均为相对隔水层的泥质岩层所组成的3层结构系统。根据地热田的面积和评价基准面的深度，确定热储的几何形状、温度、孔隙度，将评价区划分为若干个子区，为每个子区的各项参数分别赋值，计算每个子区的地热流体储存量，最后把各子区的计算结果累加，得到评价区的总储存量。计算公式为

Q储=Q容+Q弹=Adφ+AμH，

(1)

Q可=Q储·X。

(2)

式中：Q储为地热流体储存量，m3；Q容为地热流体容积储量，m3；Q弹为地热流体弹性储量，m3；Q可为可采地热流体量，m3；A为热储面积，m2；d为热储层厚度，m；H为热储顶板算起的水头高度，m；φ为热储岩石的孔隙度，无量纲；μ为含水热储层弹性释水系数，无量纲；X为可采系数，无量纲。

地热流体可开采量为地热流体储存量乘以可采系数X，可采系数的大小取决于热储层的岩性、厚度、孔隙裂隙发育情况。《华北地热》[4]及《中国地热资源》[5]采用最大允许降深法对地热流体可开采量进行了估算。开采期最大水位降深为150 m时，估算结果为馆陶组热储地热流体可开采量占总储存量的1.0%，明化镇组约占1.2%； 《京津冀油区地热资源评价与利用》[6]采用最大允许降深法对冀中部分地区进行了估算，馆陶组热储地热流体可开采量约占总储存量的1.0%，明化镇组约占1.5%； 《束鹿—宁晋地区地热资源评价研究报告》[7]采用最大允许降深法对1 980 km2进行了估算，馆陶组热储地热流体可开采量约占总储存量的1.8%，明化镇组约占2.5%； 《河北地热》[8]采用最大允许降深法估算，衡水地热田馆陶组热储地热流体可开采量约占总储存量的0.97%[9]。

由此可见，开采系数的取值决定了可开采量计算的精度。根据目前的研究成果，馆陶组热储开采系数约为1%，明化镇组热储开采系数约为1.5%。表1给出了利用热储法对河北平原地热流体储存量和可开采量的计算结果。

表1 河北平原区地热田地热流体可开采量(热储法)

2.1.2 解析法

解析法基本原理是将热储概化为均质、各项同性、等厚、各处初始压力相等的无限承压含水层，采用非稳定流泰斯公式计算单井的开采量、水位随开采时间的变化量，从而计算出在给定压力允许降深下地热流体的可开采量，其计算精度相对较低。计算公式为

(3)

(4)

式中：Q可为地热流体可开采量，m3/d；Q单为单井地热流体可开采量，m3/d；S1为计算区水位降深，m；S2为附加(单井)水位降深，m；R1为开采区半径，m；R2为单井控制半径，m；μ为热储含水层弹性释放系数，无量纲；t为开采时间，d；T为导水系数，m2/d；r为抽水井半径，m。

《河北地热》[8]中地热水可开采量是按照50 a开采期限内，区域水位降深和单井水位降深之和不大于150 m时，求得的地热水最大开采量作为计算地热水的可开采量。表2为利用解析法对河北平原37处地热田地热流体可开采量的计算结果。

表2 河北平原37处地热田地热流体可开采量(解析法)

从表2可以看出，在对馆陶组开展可开采量计算过程中，利用解析法计算得出的结果比热储法计算值偏大2.61倍，明化镇组热储偏大2.04倍。从表3可以看出，解析法计算的地热流体可开采量占总储存量的比例随着计算区面积的增大而减小。由计算公式可以看出，当计算区面积小于100 km2时，全部抽水井影响范围之和占计算区总面积的8%～30%。而当计算区面积大于1 000 km2时，全部抽水井影响范围之和仅占计算区总面积的2%～8%，90%以上的计算区面积地热水位远小于区域水位降深，因此其计算的地热流体可开采量占总储存量的比例较小。

表3 地热流体可开采量占总储存量比例统计表(解析法)

2.1.3 统计分析法

统计分析法是建立开采量和水位降低值之间的相关统计模型，用来预测地热田在定量开采条件下的水位变化趋势。可采用的统计分析法包括相关分析法、回归分析法和时间序列分析等方法。用于预测的模型应具有较高的相关系数，并且预测的时限不应超过实际监测资料的时段长度。

根据计算区地热井水位监测和地热水开采量统计资料，建立累计开采量与区域水位降深的相关方程，即

Q=f(s)。

(5)

式中：Q为累计开采量，万m3；s为区域水位降深，m。

当地热水位下降150 m(区域水位下降值与单井水位降深之和)，单井水位降深按20 m计算，区域水位下降值S为130 m时，计算累计开采量减去已开采量作为可开采量。

选取5处具有10 a左右监测数据的地热田，采用统计分析法对地热流体可开采量进行了计算。从表4可以看出，馆陶组热储地热流体可开采量占总储存量的1.04%～5.00%，与热储法用开采系数计算值相差不大。深州地热田热储层特征较为特殊，其周围断裂较为发育，且多为阻水断裂，地热田以外地热水补给量相对较少，且由于近年来地热水开采量较大，超采严重，造成地热水位迅速下降，2016年其静水位接近100 m，从而计算的地热流体可开采量也相对偏小。

表4 河北平原区地热流体可开采量(统计分析法)

2.1.4 数值模拟法

采用数值模拟法应首先对地热田的地质条件进行概化，建立概念模型。模型应反映热储和盖层的分布规律、地热水补给来源、地热系统热源、重要的热流体和热传递通道、流体动力特征等，为满足地热资源计算需要，宜概化为三维模型。数值模型的求解方法包括有限差分法、有限单元法和边界元法等，均应掌握地热水的补给条件、水动力特征、地热田内温度场特征、热储渗透率、孔隙度、贮存系数等参数，还要收集地热田的监测资料，包括地热田历史开采量、开采井水位和水温变化数据。

经过校正的模型可以用来预测地热田未来一定开采量情况下的压力变化趋势，从而计算出地热流体的可开采量。其理论基础为，在达西定律和质量守恒定律基础上，得到地下水运动基本微分方程为

(6)

式中：Kxx为x方向的渗透系数，m/s；Kyy为y方向的渗透系数，m/s；Kzz为z方向的渗透系数，m/s；h为水头，m；ω为源汇项，表示单位时间从单位体积含水层流入或流出的水量，m3/s，其中流入为正值(表示源)，流出为负值(表示汇)；Ss为贮水率，指含水层地下水水头降低(或升高)一个单位时，由于含水层垂向压缩和地下水的弹性膨胀从单位体积含水层释放(或储存)的水的体积。

本文对辛集市区馆陶组热储进行数值模型计算，将辛集市馆陶组开采井连接，向外扩1 km后的区域作为研究区，面积46.1 km2，垂向上将研究区内馆陶组含水层作为目的层，根据馆陶组热储厚度，将含水层概化为一层。由于研究区不是一个独立完整的水文地质单元，边界范围为人为指定，由于区内各热储之间发育较厚的隔水层，且开采时间较短，使各热储之间水力联系较小，本次不考虑上覆明化镇组和下伏基岩热储对研究区的影响作用，因此模型顶部和底部边界为零流量边界。辛集市馆陶组热储常年开采，已经形成十分明显的水位降落漏斗，整个边界均为补给边界。经过对水文地质条件的概化处理，计算区水文地质概念模型为： 由均质各向同性的承压水、孔隙含水层组成的具有二类边界的二维地下水渗流模型。

本次以2010年11月供暖前的水位高程为模型的初始水位，边界条件公式为

(7)

式中：T为导水系数，m/s；H为煤一已知函数；q为边界上沿法线方向地单位面积流入量，m3/s;n为边界外法线的单位矢量。

计算区平面上共剖分有效矩形网格单元4 592个，网格大小0.1 km×0.1 km，垂向上为一层，模拟期为2010年11月到2015年11月。根据地热水开采时间及地下水动态特征，将一年划分为3个时间段，共15个时间段。

模型建立后，把初始水位、渗透系数、贮水率、地热井开采量及边界条件等初始值代入模型，计算出各观测孔在各时段的水位值，并将其计算水位与实测水位进行比较，反复调整渗透系数、贮水率和边界条件； 当计算值与实测值之间的误差不超过5%时，即可认为调整后的参数符合真实的地热地质环境。本次水位资料选取了研究区14眼馆陶组地热井的观测数据，在95%置信区间内，对模型进行修正拟合。这14眼既是开采井也是监测井，是2010—2015年的连续水位情况，为本次数值模拟反演、拟合的基础数据。水位拟合的时间点选取为每个应力期末，即每年的3月15日、7月15日和11月15日。

本次地热水可开采量的确定原则为在50 a开采期内，馆陶组静水位埋深不超过200 m。辛集市区2015年7月地热水平均水位埋深为95.51 m，今后开采期内静水位最大允许下降速率为3.07 m/a，求得年可开采量为220.97万m3。

对计算区以开采期为50 a，静水位埋深不超过200 m为限，采用解析法和统计分析法对地热流体可开采量进行了计算(表5),结果显示采用解析法计算结果最大，热储法计算结果最小。热储法的开采系数是根据全省地热资源特点给出的平均值，在进行小范围地热资源评价时，其结果偏差较大； 辛集市区计算面积较小,仅为46.1 km2，采用解析法计算时评价结果偏大。采用统计分析法和数值模拟法计算结果适中，可信度相对较高。

表5 地热流体可开采量多种方法计算结果对比

2.2 各计算方法的适用条件

2.2.1 热储法

热储法计算实质是对计算深度内“静储量”的估算。不考虑侧向补给和越流补给的静态储存量，适合于勘查程度低，对热储分布范围、温度分布仅大致掌握，无地热井抽水试验数据的地热田。其可采系数的大小取决于热储层的岩性、厚度、孔隙裂隙发育情况，根据经验确定，人为因素较大。由于可采系数是根据全省范围内地热资源发育特点给出的平均值，因此当进行单个地热田或小范围地热资源评价时，常常得出与实际情况不符的结论。而对于区域性地热资源评价时，却显示出简单准确的特点。

2.2.2 解析法

该方法适用于地下水系统为无限边界局部开采的资源评价，考虑了地下水侧向补给因素，但未考虑水力坡度对侧向补给变化的影响。该方法适合勘查程度较低，有少量地热开采井，有产能试验数据的地热田。计算地热可采资源量时，受到研究区面积、热储类型、热储埋深、热储压力、热储岩性特征等多种因素的影响控制，受控因素不同，地热可采资源量会有很大不同。由计算原理可知，该方法更适用周边无开采井的集中开采区的地热资源评价。

2.2.3 统计分析法

该方法是一种以开采量与水位监测实测数据为基础的评价方法。该方法对资料的准确性及详实性要求较高，适用于勘查程度较高、已开发利用多年、具有多年动态监测资料的地热田。计算初期要考虑资料的合理性、水位监测数值的准确性、开采量调查与实际是否相符等问题，这些都会直接影响计算结果。另外，选择不同的回归方程，计算结果会有所不同，因此该方法更适用于开采时间较长、监测资料齐全的地热开采区，且预测时限不超过监测时长。

2.2.4 数值模拟法

该方法借助计算机程序对研究区域进行剖分，采用有限单元法或有限差分法对剖分单元进行离散，可将动态的水位、水量、补给情况等因素输入模型中，也可对渗透系数、贮水率等进行分区及反演，从而调整不合理的参数值，使其更符合实际的地热地质环境。因此该方法是一种高精度评价方法，适用于勘查程度较高，并且具有一定时期的开采历史，具有较齐全监测资料的地热田。由于热储层分布、地温场分布、边界条件和水文地质参数具有各向异性和非均质性等，拟合后的数值模型很难对热储层进行准确模拟，只能对地热田的地质条件进行概化。经校正的模型可认为是在当前勘查程度下的可靠的模型。辛集市的模拟结果显示，拟合后的部分模型参数和试验取得的参数差距较大，水动力流场图也有一定的差别。因此，根据模型进行预测的结果也同模型的可靠程度有密切的关系[10-11]。

3 河北省深部岩溶热储分布与资源现状

河北省岩溶热储层主要赋存于古生界和中新元古界地层中。其岩溶裂隙发育，构成深部热水储集层。基岩热储被新生界沉积层掩盖，有利于储集层的聚热和保温，在基岩隆起(古潜山)部位构成有重要开发利用价值的地热田。由于岩溶热储具有温度高、易回灌和可持续性好等特点，岩溶热储将成为今后地热开发的主要热储，但其勘查程度较低，本文仅采用热储法对地热流体储存量和可开采量进行计算评价。

3.1 古生界寒武系—奥陶系热储

寒武系—奥陶系热储(4 000 m以浅)分布面积约1.96万km2。其中奥陶系热储层由中奥陶统峰峰组、磁县组与马家沟组和下奥陶统亮甲山组及冶里组的碳酸盐岩地层构成，厚度为600～800 m。主要分布于沧县台拱部分地区和冀中台陷牛驼镇、临清台陷的宁晋及黄骅台陷的黑龙村孔店等地。热储岩溶发育程度因其顶面盖层不同有较大差异，单井涌水量150～1 500 m3/d，水温可达34～73 ℃，是本区富水性较好的热储之一。

寒武系热储层由上寒武统凤山组及下寒武统府君山组碳酸盐岩地层构成，钻探揭露厚度为30～70 m，岩性为褐灰色细晶白云岩、灰质白云岩，主要分布于冀中台陷的霸州、临清台陷的宁晋和黄骅台陷的马头营台凸，热储温度为70～78 ℃(图2)。

图2 寒武系—奥陶系热储分布图(4 000 m以浅)

3.2 中新元古界蓟县系、长城系热储

河北省平原区蓟县系、长城系热储(4 000 m以浅)分布面积约1.53万km2。其中蓟县系雾迷山组热储是一套浅海相沉积，岩性为白云岩、中部夹泥质白云岩，总厚度达300～1 000 m以上，主要分布在牛驼镇断凸、沧县台拱、宁晋断凸一带，埋深800～2 000 m，高阳台凸埋深3 000 m左右。该热储大部分地区直接被新近系及古近系所覆盖，经历了漫长的地质时期剥蚀风化和淋滤作用，风化壳厚度较大，一般为100 m，最厚为200 m，岩溶裂隙发育，连通性好。单井涌水量 400～1 500 m3/d，井口水温达到 60～80 ℃，最高可达104 ℃。该热储层是本区最好的基岩热储层。

长城系高于庄组热储层主要分布于牛驼镇、河间等地，厚度为420～980 m，岩性为灰色白云岩、角砾岩和泥质白云岩，单井涌水量 400～1 500 m3/d，井口水温达到 50～80 ℃(图3)。

此外，根据最新资料研究，在沧州任丘、河间以及衡水的深州等地，深度6 000 m左右有望获取超过150 ℃的高温地热资源，将为河北省地热发电奠定资源基础，是实现河北省地热发电突破的关键部位[12-15]。

图3 蓟县系、长城系热储分布图(4 000 m以浅)

3.3 岩溶热储资源现状

根据最新调查评价成果显示，河北省岩溶热储(4 000 m以浅)地热资源量为152.5×1016kJ，地热流体储存量为3 035亿m3，地热流体可采量为151亿m3(表6)，按开采期50 a计算，年可开采量约3亿m3。截至2015年底，河北省共有岩溶热储地热井183眼，其中寒武系—奥陶系热储地热井26眼、长城系—蓟县系热储地热井157眼，约占河北省地热井总数的14%，主要分布于廊坊、保定和沧州3市，年开采量约1 300万m3，约占年可开采量的4.3%。

表6 平原区岩溶热储地热资源量统计表(热储法)

3.4 岩溶热储回灌

2012年河北省地热资源开发研究所在牛驼镇地热田开展了岩溶裂隙热储回灌试验，回灌试验为同层热储一采一灌无压回灌模式，抽灌井距离为410 m。回灌井井深2 002 m，利用热储层为蓟县系铁岭组和雾迷山组岩溶裂隙热储，岩性为白云岩、灰岩，岩溶裂隙较为发育。回灌主要设备包括过滤器、排气罐、消气器、温度计、压力表和电磁流量计，回灌水温在40 ℃左右。经过52 d共1 248 h的回灌试验，总回灌量达8.3万m3，最大回灌量可达83.15 m3/h，验证了牛驼镇地热田内蓟县系铁岭组和雾迷山组岩溶裂隙热储的回灌能力。回灌过程中，抽水井的水位、水量和水温未受到回灌影响，也进一步证明了地热回灌的可行性，回灌水温不会在短时间内达到抽水井而形成热突破[16-17]。

河北省现开展岩溶裂隙型热储生产性回灌的地区有2处，分别为雄县和容城，回灌井共28眼。目前，雄县城区已有地热井38眼，其中地热开采井22眼，地热回灌井16眼，全部开采蓟县系雾迷山组岩溶裂隙热储，地热水主要用于小区、单位建筑物供暖，供暖期120 d。2014—2015年供暖期地热水开采总量为383.59×104m3，回灌总量约383.59×104m3，实现了蓟县系雾迷山组岩溶裂隙热储地热水100%生产性回灌。

4 结论

(1)在现有经济技术条件下，河北省平原区地热水总可开采量达821.9亿m3，其热量相当于7.3亿t标准煤。2015年开采量为9 085.5万m3，利用的热量相当于82.2万t标准煤，减排二氧化碳196.13万t、二氧化硫1.39万t、氮氧化物0.49万t、粉尘0.66万t、灰渣0.08万t，可节省治理费用约2.3亿元，经济效益与环境效益十分显著。

(2)地热流体可开采量计算方法有多种，在地热田前期勘查阶段可采用热储法进行地热资源量估算； 开采初期可采用解析法初步评价地热流体可开采量，作为开采量控制的依据，防止地热资源超采； 在大规模开采阶段，地热田达到详查或勘探程度，可采用统计分析法和数值模拟法进行地热流体可开采量评价。

(3)河北省岩溶热储具有分布广、温度高、易回灌和可持续性好等特点，岩溶热储将成为今后地热开发的主要热储层。下一步应加强深部地热资源探测和地热尾水回灌技术研究。

参考文献：

[1] 刘忠凯,吕文斌,冯来全,等.河北省地热资源调查评价与规划报告[R].衡水:河北省地矿局第三水文工程地质大队,2015.

[2] 苏永强.河北省地热能开发利用"十三五"规划报告[R].衡水:河北省地矿局第三水文工程地质大队,2016.

[3] 张银妹,赵阳,刘江涛,等.河北省地热资源可开采量评价研究报告[R].衡水:河北省地矿局第三水文工程地质大队,2016.

[4] 陈墨香.华北地热[M].北京:科学出版社,1988.

[5] 陈墨香,汪集旸,邓孝.中国地热资源——形成特点和潜力评估[M].北京:科学出版社,1994.

[6] 阎敦实,于英太.京津冀油区地热资源评价与利用[M].武汉:中国地质大学出版社,2000.

[7] 邓孝.束鹿—宁晋地区地热资源评价研究报告[R].石家庄:中国地质科学院水文地质环境地质研究所,1990.

[8] 张德忠,刘志刚,卢红柳,等.河北地热[M].北京:地质出版社,2013.

[9] 赵静,王旭开,万力，等.深层地热资源评价中的回收率问题[J].地质科技情报,2008(6):89-92.

[10] 张中祥,韩建江,徐健,等.沉积盆地型层状热储可采资源量计算方法探索[J].水文地质工程地质,2005,32(1):73-77.

[11] 孙颖,许辉熙，刘久荣，等.中低温地热田的地热资源计算评价[J].安徽农业科学，2009，37(14)：6535-6537.

[12] 李郡,王全凯,王琰,等.河北省地热资源现状调查评价与区划报告[R].衡水:河北省地矿局第三水文工程地质大队,2015.

[13] 王永波,丁文萍,田月,等.河北牛驼镇地热田高温地热水成因分析[J].城市地质,2016,11(3):59-64.

[14] 李宜程.山东省中深层地热资源开发潜力分区评价[D].济南:山东建筑大学,2016.

[15] 吕东亮,林黎,林建旺,等.岩溶热储地热能可持续开发潜力的模糊综合评价[J].勘察科学技术,2008(4):20-23,52.

[16] 刘福东，黄加力，李子果，等，河北省牛驼镇地热田地热尾水回灌试验[R].衡水市：河北省地矿局第三水文工程地质大队.2012.

[17] 周阳，邓念东，王凤，等.浅层地热能适宜性分区结构的分形原理[J].中国地质调查，2017，4(1)：18-23.