江淮地区地热水资源评价

董 超，杜宁宇，骆祖江*

(1.河海大学地球科学与工程学院，江苏南京 210098；2.河海大学河海里尔学院，江苏南京 210098)

0 引言

地热资源蕴藏丰富且无污染，是一种可再生能源[1-2]。随着石油、煤炭等传统能源逐渐枯竭，合理开发地热资源[3]，对缓解能源束缚和环境压力，提高经济增长的质量具有重要的意义，也是实现大力开发新能源，实现能源清洁发展、安全发展、环境友好发展和可持续发展的一项重要的战略举措[4]。江苏处于中国东部沿海高热流地热异常带，大地热流强，地温梯度高，热储层分布广、厚度大、埋藏条件好，具备形成中温地热资源的地质条件[5-7]。以往的地热水资源评价，往往简化为只评价水量，将水温看作恒定不变，或将水量和水温完全分开来评价，和实际不相吻合。本文以江苏省江淮生态经济区东南部地区地热水系统为例，根据研究区地质、水文地质及地热地质特征，在概化出地热水系统概念模型的基础上，建立了研究区地热水非稳定渗流与热量运移三维耦合数值模型[8]，同时评价了现有地热水井在50m降深条件下的可开采资源量及温度的变化[9]，为该地区地热水的可持续开发利用提供了科学依据，具有较强的工程实际意义[10]。

1 地质及水文地质条件

研究区位于江淮生态经济区东南部的宝应、建湖、高邮、兴化4个县(市)。地形以平原为主，总体地势西高东低。东、南侧以已存在的江淮生态经济区的范围为界，西侧以高邮湖、白马湖等河流湖泊为界，北侧为根据行政区域及自然边界人为划定的边界，总面积约为6 900km2。

研究区地处苏北盆地中部，地层隶属下扬子地层分区，地表均为第四系覆盖，地层分布不全，3 000m以浅分布的地层呈现一老一新特点。即下部为震旦系—下古生界，构成苏北盆地基底，上部为中新生代碎屑岩，分布稳定。

研究区所处的苏北盆地处位于中国东部沿海高热流地热异常区，大地热流强，地温梯度高，具备形成中高温地热资源的地质条件，地热资源成因主要为传导型和中低温对流型[11]。区内地层深部分布有厚度巨大的碳酸盐岩，这些地层在地质历史时期经历了强烈的构造变动和岩溶化作用，岩溶裂隙发育，为热水型地热资源的赋存提供了良好的储存空间。区内断裂构造十分发育，但规模深浅不一，在凸起的边部、凹陷外斜坡带、凹陷边缘断阶带、隆起区的边缘地带与凹陷边缘断阶带均受区域性深断裂控制，具有极好的水热循环条件，也是地热资源富集区，主要发育NE、NNE、EW及NW向四组断裂，NNE和NE向断裂为大、中型新构造断裂构造带，与后期NW向张性、张扭性断裂交汇，构成了深部热流和地下水循环的良好通道[12-13]。

区内地下水根据含水介质类型划分为松散岩类孔隙水、碳酸盐岩类岩溶裂隙水及碎屑岩类裂隙水。松散岩类孔隙水在区内广泛发育，主要赋存在盐城组和第四纪松散堆积地层中的砂层及砂砾层中；碳酸盐岩类岩溶裂隙水赋存于震旦系和古生界碳酸盐岩之中；碎屑岩类裂隙水主要赋存在古近系基岩中。本次对区内6口地热井进行了调查(表1)。

表1 研究区地热水井概况

2 地热水系统概念模型

研究区热水系统顶部一方面接受大气降雨的补给，是一补给边界，另一方面地下水又通过其蒸发，是一排泄边界。地热水通过断裂、裂隙的连通进行循环，水流呈非稳定流状态，区内断裂构造十分发育，规模深浅不一，各含水层均为非均质各向异性，各层之间均发生水力联系，地下水位受大气降水和蒸发季节性变化的影响，为三维非稳定流。由于研究区范围较大，边界距离地热水井较远，因此模型四周均概化为第一类定水头边界，模型底部为隔水边界。

研究区热水系统顶部接受大气降雨的补给，混入冷水，是一热量排泄边界；底部通过断裂、裂隙的联通、循环，将深层的较高温度带到浅层热储含水层，是一热量补给边界，根据增温梯度推算设定为一恒温边界，故顶底部概化为第一类热量边界。计算范围内设定地下水和含水介质骨架的热动平衡是瞬时完成的，将热储含水层四周概化为第二类热量边界。

3 地热水非稳定渗流与热量运移三维耦合数值模型

根据上述概化出的地热水系统概念模型，建立相应的地热水非稳定渗流与热量运移三维耦合数值模型[14]。地热水渗流场与温度场通过地热水流动方程进行耦合计算。

3.1 地热水非稳定渗流数值模型

取坐标轴方向与热储含水层各向异性主渗透系数方向一致，建立如下数值模型：

(1)

式中：Kxx、Kyy、Kzz为各向异性主方向渗透系数，m/d；h为点(x,y,z)在t时刻的水头值，m；W为源汇项(1/d)；t为时间，d；Ω为计算域；h0(x,y,z,t0)为点(x,y,z)处初始水头值，m；q(x,y,z,t)为第二类边界上单位面积的补给量，m/d；cos(n,x)、cos(n,y)、cos(n,z)为流量边界外法线方向与坐标轴方向夹角的余弦；μ为饱和差(自由面上升)或给水度(自由面下降)；Ss为储水率，L/m；Γ1为第一类边界；Γ2为第二类边界；Γ3为自由面边界。

3.2 地热水热量运移数值模型

假设地热水和含水介质骨架的热动平衡是瞬时完成的，并忽略由于温度差引起水的密度不一样而引起的上下自然对流的影响，建立地热水三维热量运移数值模型如下：

(2)

式中：λx、λy、λz为各方向水的热动力弥散系数；cw为水的热容量，J/(m3·K)；c为含水介质的热容量，J/(m3·K)；vx、vy、vz为地热水渗流速度分量，m/d；T0(x,y,z)为点(x,y,z)处初始温度值；T1(x,y,z,t)为第一类边界的温度函数；Γ1为第一类边界；Qc为热源汇项，Qc=cwW(TQ-T)；TQ为源汇项的温度；n为边界外法线向量；Q(x,y,z,t)为第二类边界上已知的热量或热流函数，J/(m·d)；Γ2为第二类边界。

3.3 地热水非稳定渗流与热量运移三维耦合数值模型

地热水流运动方程：

(3)

将(1)式与(2)式通过地热水流运动方程式(3)耦合在一起，构成研究区地热水非稳定渗流与热量运移三维耦合数值模型。

4 模型求解

上述模型采用伽辽金有限元法求解，首先将地下水渗流与热量运移的性代数方程组代入求解方程，其插值称为误差函数或剩余，从某种平均意义上讲希望该误差为零。让微分方程的近似解和精确解之差的总剩余量最小，用近似解来代替精确解。

求解渗流问题的代数方程形式：

(4)

求解热量运移问题的代数方程形式为：

(5)

式中：[K]为总体渗透矩阵；{h}为未知节点水头列阵；[P]为总体储水矩阵；[U]为总体自由面流量补给矩阵；{F}为已知的右端项；[G]为总体弥散矩阵；{T}为未知节点温度列阵；[R]为总体对流矩阵；[M]为总体热容量矩阵；{Z}为已知的右端项。

再对式中的时间采用隐式差分格式进行离散，整理后可以得到：

(6)

(7)

利用上述方程，根据初始条件下的渗流场水头和温度分布，即可计算任意时间的渗流场水头和温度分布。

5 模型识别验证

根据研究区地下水系统结构特征及几何形状，对研究区进行三维剖分。平面上对地热水井周围和主要的导水控水断裂进行网格加密，最终形成每层36 082个三角形单元和18 482个节点。垂向上考虑到地热水井出水层位及层厚的影响，由上往下共剖分8个层位，9个计算层面，共计288 656个单元，166 338个节点。计算域空间单元剖分见图1。

利用现有兴热1井、兴热2井、RGM1井、RGS1井、RBQ1井、宝热1井的抽水试验资料对模型进行识别、验证，并根据研究区各地热水井初始水位、温度、1 500m地温分布图及实测地温增温规律，确定模型初始流场和初始温度场。图2举例说明了模型第一层地热水初始流场和初始温度场。各开采井作为观测井同时进行水位和水温拟合。

图1 计算域空间剖分Figure 1 Model space subdivision of computational domain

图2 模型第一层初始流场及初始温度场分布Figure 2 Distribution of initial flow field and initial temperature field in the first layer of modela.初始流场图(m)；b.初始温度场图(℃)

全区共分108个参数分区，上部三层根据地热水井分布和地层岩性特征进行分区，下部地层受断裂带影响，根据断裂带划分为断裂影响带和无影响区。通过反演计算，获得了各层各参数分区的参数值，图3举例说明了第八层的水文地质参数分区，表2举例说明了各参数分区的参数值。各地热水井抽水试验实测数据与计算数据拟合情况见图4。从拟合结果来看，计算曲线与实测曲线拟合精度较好，总体变化趋势一致。通对比对可以发现，模型计算满足有关精度要求[15]，可用于研究区地下水渗流场与温度场的模拟预测。

图3 第Ⅷ层参数分区Figure 3 Model parameters partition of the eighth group

表2 第Ⅷ层各参数分区参数

续表

6 地热水可开采资源量规划评价

利用识别、验证后的地热水非稳定渗流与热量运移三维耦合数值模型，对研究区现有6口地热水井所在的热储含水层2022-03-01至2032-03-01的地下水流场与温度场变化进行预测，每一年分为12个应力期，共计120个应力期。由于研究区地下水开采量少，水位、温度变化幅度小，初始流场和初始温度场均采用模型识别、验证阶段的初始数据，大气降雨量数据采用2000—2020年逐月平均值，气温采用年平均气温15℃。经模型模拟预测，各开采井在2022年3月1日至2032年3月1日未来10a内的水位、温度变化历时曲线见图5。

图4 研究区水位观测值与计算值拟合Figure 4 Fitting results of the observed and calculated water level values in the study area

图5 运行10年水位及温度变化历时曲线Figure 5 Diachronic curve of water level and temperature changes during 10 years of operation

从模型运行得出的结果可以看出，按照50m水位降深确定的开采量进行开采，随着开采的持续，整个研究区水位逐渐下降10a后达到平衡，地下水位处于相对稳定状态；各地热水井温度变化均不超过1℃。其中兴热1井、兴热2井、RGM1井、RGS1井和宝热1井整体温度场相对于初始状态有所上升，RBQ1井整体温度场相对于初始状态有所下降。由于地热水资源的持续开采，温度在局部范围内将发生一定的变化，将会产生冷热堆积，对周围水环境产生一定的影响。

根据模拟计算结果，整个研究区的6口地热水井按照不超过50m降深开采地热水，各地热水井可开采资源量及水位变化统计见表3，合计6口地热水井可开采资源量为5 1290 m3/d。

表3 研究区各地热水井预测水位变化及取水量

7 结论

1)通过建立江淮地区地热水非稳定渗流与热量运移三维耦合数值模型，可以同时预测未来开采条件下水位和水温的变化趋势，并可根据水位和水温两个指标的变化，规划评价地热水的可开采资源量，实现地热水的可持续开采利用。

2)江淮地区现有各地热水井10a内水位降深稳定在50m的条件下，各地热水井地热水的可开采资源量分别为兴热1井3 870 m3/d、兴热2井7 630 m3/d、RGM1井3 720 m3/d、RGS1井890 m3/d、RBQ1井1 480 m3/d、宝热1井33 700 m3/d。全区可开采资源量为51 290 m3/d，各井的温度变化均不超过1℃。