2. 江苏省地质矿产局第一地质大队, 南京 210041
2. The First Geological Brigade of the Bureau of Geology and Mineral Resources of Jiangsu, Nanjing 210041, China
0 引言
南京汤山是我国知名的温泉之乡。近年来，随着国民经济的发展和人民生活水平的提高，汤山地区地热水的开采量日益增加，水位降深不断加大，早期的自流井早已不复存在。因此，准确评价汤山地区的地热水资源，合理规划开采地热水已成为当务之急。
目前常见的地热水资源评价方法主要是热储法或解析法，我国《地热资源地质勘查规范》(GB/T11615-2010)[1]中也给出了几种常用的储量计算方法。这些评价方法通常只是计算了地热水的静态储量和热量，缺少对地热水资源的可开采量及开采过程中地热水水温变化的评价[2-7]，难以满足为地热水的规划开发提供必要科学依据的需求。笔者通过对南京汤山地区地热水系统形成的地质、水文地质背景特征的研究[8-9]，建立了南京汤山地区地热水非稳定渗流与热量运移三维耦合数学模型，在对模型进行校正、识别的基础上，预测了汤山地区地热水系统20170401—20220331水位与水温的变化趋势，并依据文献[1]规范，按地热水水位降深不超过50 m，水温变化不超过已有实际观测值最大变幅的10%为控制要求，规划评价了南京汤山地区地热水资源的可开采量。
1 研究区域概况研究区位于南京市东部、江宁区东北端、宁镇山脉西段南缘，属低山丘陵地貌。北部以沪宁高速为界，南部以宁杭公路为界，东部为汤山镇推测断层(F10)，西部为栖霞—汤山头推测断层(F7)，面积约为1.6×107 m2，见图 1。
|
| 图 1 研究区地理位置图(a)及地质简图(b) Fig. 1 Geological position(a) and geological conditions(b) of the study area |
|
|
汤山地区的地下水可分为碳酸盐岩岩溶裂隙水、碎屑岩类及火成岩类裂隙水、松散岩类孔隙水3种类型。其中碳酸盐岩多裸露在山体表层，碎屑岩类及火成岩类裂隙水赋存于碎屑岩类和火成岩类构造裂隙及风化裂隙中，松散岩类孔隙水主要赋存于第四系残积-坡积层灰黄色、黄色亚砂土或亚黏土夹碎石中。
研究区内地质构造极为发育，主要的褶皱构造有汤山—仑山复背斜，断裂构造以东西向区域性压扭性断裂与近南北向张扭性断裂为主。区内北东向的汤山—东昌街断裂破碎带在浅部与北西向的庙山—狼山断裂及一系列的南北向断裂交会，具有良好的控热、控水性，为深部热流向上运移和储存提供了通道和空间。地热水储层为一套寒武系、奥陶系碳酸盐岩，大部分直接裸露于地表，仅南部局部地段存在第四系盖层。热传递以对流和传导两种形式同时存在，且以对流型为主。研究区浅表地温变化范围为16.7~25.5 ℃，平均21.2 ℃，埋深5~10 m的浅水水温为13.5~19.1 ℃；埋深60~ 200 m的碳酸盐岩溶裂隙水水温为18.9~26.2 ℃；埋深200 m以下岩溶裂隙水水温53.0~80.3 ℃。区内地热水的开发主要分布在汤山山体东侧的汤山镇区和西侧的侯家塘附近，共有14口开采井(图 1)。目前实际开采地热水的水井共10口，各井井深为100.0~ 2 000.0 m，温度为30.0~80.0 ℃，各地热水井2016年的水位标高为36~38 m，现仍处于下降趋势，全区总的地热水开采量近年来基本稳定在9.46×105 m3/a。
2 地热水系统概念模型研究区地热水系统四周均为透水边界，接受区域侧向径流的补给或排泄，概化为流量边界。地热水储层寒武系、奥陶系碳酸盐岩在汤山中部出露于地表，向四周逐渐深埋。顶板埋深为0~1 300 m，底板埋深为1 300~3 500m。裸露区直接接受大气降水的入渗补给。受构造作用及岩溶发育深度的控制，系统内构造和岩溶裂隙的分布极为不均，为非均质各向异性。地热水水位除随季节性波动外，还随现状开采呈现下降趋势，为三维非稳定流。
研究区地热水系统顶部接受大气降水的补给，混入冷水，是一热量排泄边界; 底部通过断裂、裂隙，将深部热流向上传递，是一热量补给边界。故顶底部概化为第一类定温度边界，四周均概化为第二类热量边界。
各地热水开采井均按点井处理，大气降水入渗补给按面状入渗补给强度处理。
3 地下水非稳定渗流与热量运移三维耦合数学模型 3.1 地下水三维非稳定渗流数学模型取坐标轴方向与含水层各向异性主渗透方向一致，建立地下水三维非稳定渗流数学模型：
(1) 式中：Kxx、Kyy、Kzz为各向异性主方向渗透系数(m/d)；h为点(x, y, z)在t时刻的水头值(m)；μs为含水层的贮水率；W为源汇项(1/d)；t为时间(d)；Ω为计算域；t0为初始时刻；h0(x, y, z, t0)为点(x, y, z)处初始水头值(m)；h1(x, y, z, t)为t时刻第一类边界(Γ1)上各点的已知水头值(m)；q(x, y, z, t)为第二类边界(Γ1)上单位面积的补给量(m/d)；cos(n, x)、cos(n, y)、cos(n, z)为流量边界外法线方向与坐标轴方向夹角的余弦，其中n为流量边界外法线方向；K为总渗透系数(m/d)；μ为饱和差(自由面上升)或给水度(自由面下降)；θ为浸润曲线外法线与铅锤线间的夹角；Γ3为自由面边界。
3.2 地下水热量运移数学模型根据模型的假定，含水介质骨架与地下水的热动平衡瞬时完成，忽略由温度差引发的上下自然对流的影响，建立地下水三维热量运移数学模型：
(2) 式中：λx、λy、λz为各方向水的热动力弥散系数，由地下水与含水介质骨架的导热系数、地下水的横向与纵向热弥散度、地下水的渗流速度(W/(m·K))计算得出；cw为水的热容量(J/(m3·K))；c为含水介质的热容量(J/(m3·K))；νx、νy、νz为地下水渗流速度分量(m/d)；T为点(x, y, z)处的温度函数(K)；T0(x, y, z)为点(x, y, z)处初始温度值(K)；T1(x, y, z, t)为第一类边界的温度函数(K)；Qc为热源汇项，Qc=cwW(TQc-T)；TQc为源汇项的温度(K)；Q(x, y, z, t)为第二类边界上已知的热量或热流函数(J/(m·d))。
3.3 地下水非稳定渗流与热量运移三维耦合数学模型地下水流运动方程为
(3) 其中：
通过式(3)将式(1)与式(2)耦合，构建研究区地下水非稳定渗流与热量运移三维耦合数学模型(图 2)。
|
| 图 2 模型空间剖分图 Fig. 2 Model space subdivision diagram |
|
|
根据研究区寒武系、奥陶系碳酸盐岩地热水储层、构造和岩溶裂隙、地热水开采井的空间分布特征，将研究区平面上剖分成17 292个三角形单元。垂向上从上向下剖分成6层，上部5层各层厚度均为200.00 m，其余剖分为下部第6层，共7个计算层面，底层埋深的变化根据物探资料解释确定，其余各层面的埋深均随顶板标高的变化而确定，共计121 044个计算单元和74 921个计算节点，见图 2。采用FEFLOW地下水模拟商用软件对模型进行校正、识别，选取20160317—20170317作为模型校正、识别的时段，共分36个应力期，每个应力期为10 d。各层地热水的初始流场根据实测水位给出，初始温度场根据实测资料结合地温梯度计算结果给出，图 3举例说明了模型第一层地热水初始流场和温度场的分布特征。大气降雨量和气温值由实测获得，参数初值通过地热水井单孔抽水试验和热响应试验获得。各开采井的开采层根据井身结构及抽水泵埋深位置确定，大多位于模型的第二层位。将模型表层作为地热水水位和水温的观测层位，各开采井均作为观测井同时进行水位和水温拟合。
|
| x，y为国家2000坐标系的坐标。 图 3 研究区校正识别阶段模型第一层初始流场(a)及初始温度场(b)分布图 Fig. 3 Initial flow field(a) and the initial temperature field(b) of the first layer of the model in the calibration and identification phase in the study area |
|
|
通过对模型校正识别，全区共分27个参数分区(图 4)，各参数分区的参数值见表 1。图 5举例说明了地热水井水位和水温的拟合精度。
|
| Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ为含水层号。 图 4 研究区各含水层(组)参数分区 Fig. 4 Parametric partitioning of each aquifer (group) in the study area |
|
|
| 分区号 | 主轴方向渗透系数/(m/d) | 给水度 | 储水率/(10-6/m) | 热容量/(106J/(m3·K)) | 含水层热传导系数/(J/(m·s·K)) | ||
| Kxx | Kyy | Kzz | |||||
| 1 | 1.30 | 1.30 | 0.90 | - | 1.50 | 2.661 | 2.473 |
| 2 | 4.50 | 4.50 | 3.79 | - | 2.30 | 2.661 | 2.473 |
| 3 | 0.70 | 0.70 | 0.41 | 0.015 | - | 2.343 | 2.025 |
| 4 | 2.60 | 2.60 | 1.98 | - | 2.40 | 2.054 | 2.357 |
| 5 | 3.20 | 3.20 | 2.45 | - | 2.50 | 1.873 | 2.652 |
| 6 | 4.10 | 4.10 | 2.96 | - | 2.12 | 1.413 | 2.921 |
| 7 | 2.30 | 2.30 | 1.80 | - | 1.80 | 1.413 | 2.921 |
| 8 | 2.10 | 2.10 | 1.90 | - | 1.90 | 1.563 | 2.496 |
| 9 | 1.03 | 1.03 | 0.63 | - | 1.80 | 2.661 | 2.473 |
| 10 | 0.60 | 0.60 | 0.35 | - | 1.96 | 2.661 | 2.473 |
| 11 | 1.86 | 1.86 | 1.26 | - | 2.00 | 2.343 | 2.025 |
| 12 | 1.63 | 1.63 | 0.96 | - | 2.10 | 2.054 | 2.357 |
| 13 | 1.26 | 1.26 | 1.04 | - | 2.18 | 1.873 | 2.652 |
| 14 | 0.80 | 0.80 | 0.64 | - | 1.50 | 1.563 | 2.431 |
| 15 | 0.42 | 0.42 | 0.26 | - | 1.10 | 1.413 | 2.921 |
| 16 | 1.37 | 1.37 | 1.03 | - | 1.12 | 2.661 | 2.473 |
| 17 | 1.10 | 1.10 | 0.73 | - | 1.20 | 2.343 | 2.025 |
| 18 | 1.00 | 1.00 | 0.59 | - | 1.38 | 2.054 | 2.357 |
| 19 | 0.62 | 0.62 | 0.37 | - | 0.58 | 1.873 | 2.652 |
| 20 | 0.31 | 0.31 | 0.20 | - | 0.62 | 1.413 | 2.921 |
| 21 | 0.83 | 0.83 | 0.56 | - | 0.57 | 1.532 | 2.314 |
| 22 | 0.52 | 0.52 | 0.34 | - | 0.69 | 1.463 | 2.365 |
| 23 | 0.30 | 0.30 | 0.11 | - | 0.76 | 1.964 | 2.320 |
| 24 | 0.56 | 0.56 | 0.23 | - | 0.52 | 2.146 | 2.460 |
| 25 | 0.52 | 0.52 | 0.34 | - | 0.69 | 1.463 | 2.365 |
| 26 | 0.30 | 0.30 | 0.11 | - | 0.76 | 1.964 | 2.320 |
| 27 | 0.56 | 0.56 | 0.23 | - | 0.52 | 2.146 | 2.460 |
| 注：-为无此参数。 | |||||||
|
| 图 5 研究区水位及水温拟合曲线图 Fig. 5 Water level and temperature fitting curve in the study area |
|
|
从拟合结果来看，计算曲线与实测曲线拟合精度较好，总体变化趋势一致。研究区地热水井水位年最大变幅＞30.00 m，水温最大变幅＞5.0 ℃。各地热水井水位和水温的模拟计算值与实测值的拟合误差均小于变幅的10%；模型计算误差满足精度要求，可用于研究区地热水渗流场与温度场的模拟预测。
5 地热水可开采量规划评价利用上述校正、识别后的地热水非稳定渗流与热量运移三维耦合数值模型，对研究区地热水20170331—20220331逐年的渗流场与温度场变化进行预测。每个应力期为1个月，共计60个应力期。各层的初始流场根据实际观测资料给出，各层初始温度场根据实测资料结合地温梯度计算给出，大气降雨量和气温值采用2006—2016逐月的平均值。
根据汤山温泉管委会要求，区内已有14口井未来全部进行地热水开采，在不改变开采井格局的条件下，根据相关规范要求，以20170331—20220331各地热水井水位降深不能超过50 m，并趋于稳定，水温变化不超过已有实际观测值最大变幅的10%为约束条件，以201603—201703的各地热水井的月平均开采量为基础，结合实际的开采需求，规划各地热水井的开采量。通过多方案对比分析，最终规划评价出了各地热水井的可开采量与水位、水温变化情况，见表 2。图 6举例说明了模型第二层20170401及20220331地热水流场和温度场的分布特征。全区14口地热水井总的可开采量为3.08×106 m3/a。
| 编号 | 井名 | 可开采量/(m3/d) | 第五年末变幅 | |
| 水位降低/m | 水温上升/℃ | |||
| R01 | 寺庄 | 510.0 | 49.90 | 2.5 |
| R02 | 温泉公司2# | 285.0 | 49.80 | 6.9 |
| R03 | 工人疗养院 | 785.0 | 50.10 | 2.8 |
| R04 | 地热公司2# | 1 280.0 | 49.90 | 3.4 |
| R05 | 地热公司1# | 975.0 | 49.98 | 3.6 |
| R06 | 巴厘原墅 | 850.0 | 50.10 | 3.2 |
| R07 | 洋兰 | 560.0 | 50.10 | 3.2 |
| R08 | 温泉公司1# | 1 450.0 | 50.20 | 3.0 |
| R09 | 古泉井 | 195.0 | 50.60 | 2.0 |
| R10 | 中闻集团1# | 165.0 | 50.18 | 5.5 |
| R11 | 中闻集团2# | 125.0 | 50.08 | 6.0 |
| R12 | 蒋介石别墅 | 165.0 | 50.30 | 3.3 |
| R13 | 八三医院1# | 350.0 | 50.20 | 2.8 |
| R14 | 八三医院2# | 755.0 | 50.60 | 3.6 |
|
| a. 20170401，初始流场；b. 20170401，初始温度场；c. 20220331，流场；d. 20220331，温度场。 图 6 研究区预测阶段模型第二层流场及温度场分布图 Fig. 6 Flow field and temperature field of the second layer of the model in the prediction stage in the study area |
|
|
从计算结果(表 2)可以看出，各地热水井的可开采量相差较大。随着热水井的开采，各地热水井的水位逐渐下降，水温均逐渐上升，但上升的幅度不同，这是由地热水储层的非均质性造成的。各地热水井水温的上升，一方面是由于随着地热水水位的下降，大气降雨对地热水的入渗补给相对减弱，另一方面，随着地热水储层深度的增加，地热水储层的水温逐渐上升。说明计算结果和实际情况相吻合。
6 结论1) 建立了南京汤山地区地热水渗流与热量运移三维耦合数学模型，通过对地热水水位和水温的模拟预测，规划评价了地热水资源的可开采量。
2) 南京汤山地区现有14口地热水井的可开采量相差较大，14口地热水井的总可开采量为3.08×106 m3/a。
3) 随着各地热水井的开采，水位的逐渐下降，各地热水井的水温均会逐渐上升，各地热水井水温上升的幅度不同，大部分年内水温上升2.0 ~ 3.0 ℃，个别地热水井水温上升达5.0~7.0 ℃。
| [1] |
地热资源地质勘查规范: GB/T 11615-2010[S].北京: 中国标准出版社, 2010. Geologic Exploration Standard of Geothermal Resources: GB/T 11615-2010[S]. Beijing: Standards Press of China, 2010. |
| [2] |
魏永霞, 唐仲华, 左霖. 数值模拟在地热资源评价及优化开采方案的应用研究[J]. 安徽地质, 2019, 29(1): 74-80. Wei Yongxia, Tang Zhonghua, Zuo Lin. Study on Application of Numerical Simulation in Geothermal Resource Evaluation and Optimization of Exploitation Scheme[J]. Geology of Anhui, 2019, 29(1): 74-80. |
| [3] |
Iqbal K S, Heru B P. A Natural State Model and Resource Assessment of Ulumbu Geothermal Field[J]. IOP Conference Series:Earth and Environmental Science. DOI:10.1088/1755-1315/254/1/012017 |
| [4] |
刘杰, 宋美钰, 田光辉. 天津地热资源开发利用现状及可持续开发利用建议[J]. 地质调查与研究, 2012, 35(1): 67-73. Liu Jie, Song Meiyu, Tian Guanghui. Development Situation of the Geothermal Resources and Suggestion on Sustainable Development Utilization in Tianjin[J]. Geological Survey and Research, 2012, 35(1): 67-73. |
| [5] |
朱喜, 张庆莲, 刘彦广. 基于热储法的鲁西平原地热资源评价[J]. 地质科技情报, 2016, 35(4): 172-177. Zhu Xi, Zhang Qinglian, Liu Yanguang. Evaluation of the Geothermal Resources in the Plain of West Shandong Province[J]. Geological Science and Technology Information, 2016, 35(4): 172-177. |
| [6] |
姜智超.黑龙江省绥化市地热田地热资源评价及合理开发利用[D].长春: 吉林大学, 2015. Jiang Zhichao. The Evaluation and Rational Exploitation of Geothermal Resources in the Geothermal Field in Suihua City, Heilongjiang Province[D]. Changchun: Jilin University, 2015. |
| [7] |
张思桃. 中山虎池围地热田地热资源储量计算与评价研究[J]. 地下水, 2018, 40(6): 33-36. Zhang Sitao. Study on Calculation and Evaluation of Geothermal Resource Reserves in Huchiwei (Weichi Lake) Geothermal Field, Zhongshan[J]. Ground Water, 2018, 40(6): 33-36. |
| [8] |
赵剑畏, 朱士鹏. 南京汤山地下热水的控制因素与资源前景[J]. 地质学刊, 1998(4): 52-58. Zhao Jianwei, Zhu Shipeng. Control Factors and Resources Prospect of Geothermal Water in Tangshan Hill of Nanjing[J]. Journal of Geology, 1998(4): 52-58. |
| [9] |
李爱勇, 朱春生, 杨生. 南京汤山温泉形成条件探讨[J]. 矿产勘查, 2010, 1(6): 546-549. Li Aiyong, Zhu Chunsheng, Yang Sheng. The Research on Formation Condition of Tangshan Warm Spring in Nanjing[J]. Mineral Exploration, 2010, 1(6): 546-549. |
| [10] |
朱春生.南京汤山温泉形成条件与热源探讨[C]//江西省地质学会.2015地学新进展: 第十三届华东六省一市地学科技论坛文集.南昌: 江西省地质学会, 2015: 316-319. Zhu Chunsheng. Forming Condirions and Heat Source of the Tangshan Hot Spring in Nanjing[C]//Jiangxi Geological Society. New Progress of Geosciences in 2015: Collection of the 13th Geosciences Science and Technology Forum of Six Provinces and One City in East China. Nanchang: Jiangxi Geological Society, 2015: 316-319. |
| [11] |
邹鹏飞, 邱杨, 王彩会. 南京汤山温泉区地热水成因模式分析[J]. 高校地质学报, 2015, 21(1): 155-162. Zou Pengfei, Qiu Yang, Wang Caihui. Analyses of the Genesis of Tangshan Hot Spring Area in Nanjing[J]. Geological Journal of China Universities, 2015, 21(1): 155-162. |
| [12] |
栾光忠, 邱汉学. 中低温对流型地热系统的典型成因:南京汤山地热系统的分析[J]. 中国海洋大学学报, 1998, 28(1): 160-164. Luan Guangzhong, Qiu Hanxue. The Type of Low-Medium Temperature Geothermal System of Convection Type:The Genesis Analysis of Tangshan Geothermal System in Nanjing[J]. Periodical of Ocean University of China, 1998, 28(1): 160-164. |
| [13] |
张金华, 魏伟, 杜东, 等. 地热资源的开发利用及可持续发展[J]. 中外能源, 2013, 18(1): 30-35. Zhang Jinhua, Wei Wei, Du Dong, et al. The Development, Utilization and Sustainable Development of Geothermal Resources[J]. Sino-Global Energy, 2013, 18(1): 30-35. |
| [14] |
骆祖江, 杜菁菁. 基于热平衡分析的地埋管地源热泵换热方案模拟优化[J]. 农业工程学报, 2018, 34(13): 246-254, 320. Luo Zujiang, Du Jingjing. Heat Exchange Scheme Simulation Optimization for Ground Source Heat Pump System with Buried Pipes by Thermal Equilibrium Analysis[J]. Transactions of the Chinese Society of Agricultural Engineering, 2018, 34(13): 246-254, 320. |
| [15] |
赵敬波.地下热水流动与热量运移的三维非稳定流数值模拟研究[D].北京: 中国地质大学(北京), 2015. Zhao Jingbo. 3-Dimensional Numerical Modeling of Unsteady Thermal Groundwater Flow and Heat Transport[D]. Beijing: China University of Geosciences (Beijing), 2015. |
| [16] |
周世玲, 骆祖江, 于丹丹, 等. 地下水源热泵系统热平衡预测三维数值模拟[J]. 中国煤炭地质, 2014, 26(9): 34-39. Zhou Shiling, Luo Zujiang, Yu Dandan, et al. 3D Numerical Simulation of Groundwater Source Heat Pimp System Thermal Equilibrium Prediction[J]. Coal Geology of China, 2014, 26(9): 34-39. |
| [17] |
Donaldson I G, Grant M A. An Estimate of the Resource Potential of New Zealand Geothermal Fields for Power Generation[J]. Geothermics. DOI:10.1016/0375-6505(78)90014-7 |
| [18] |
Luo Zujiang, Wang Yan, Zhou Shiling, et al. Simulation and Prediction of Conditions for Effective Development of Shallow Geothermal Energy[J]. Applied Thermal Engineering, 2015, 91: 370-376. DOI:10.1016/j.applthermaleng.2015.08.028 |