浅层地热能不同开发方案对深层地下淡水咸化效应研究

李 兆，骆祖江，杨 璐，胡 颖

（1．河海大学 地球科学与工程学院，江苏 南京 211100；2．江苏省地质工程勘察院，江苏 南京 211102）

1 研究背景

淡水是人类生产、生活不可缺少的资源，我国655个城市中超过400个以地下水作为饮用水源，然而大量研究显示滨海地区地下淡水资源有咸化的趋势［1-2］，地下淡水资源正在受到威胁，为了进一步保护地下淡水资源，咸水入侵防控成为了当今地下水科学研究的热点。

研究认为地下水开采导致地下水位下降，咸、淡水分界面向内陆侵袭是导致地下水咸化的主要原因［3-4］，在滨海地区开展诸多关于海水入侵的研究［5-7］。然而很多滨海地区地下水受沉积环境影响，天然状态下水质较差，特别是在孔隙地下水常呈现 “上咸下淡”的格局，即浅部含水层贮存咸水，深部含水层贮存淡水［8］，周慧芳等［9］就南通地区地下水咸化进行了系统研究，提出淡水咸化原因有海水入侵、水—岩作用和浅部咸水越流入侵。Ma等［10］针对浅部咸水越流入侵深部淡水这一问题，提出适量开采浅层地下咸水，减小与深层淡水之间的水力梯度，以缓解深部淡水咸化，为咸水入侵防控提供新的思路，但是开采浅部咸水如何利用成为了新的问题。

近年来我国大力发展浅层地热能以实现清洁取暖、节能减排，浅层地热能的开发利用可以有效解决浅层咸水的利用问题。地下水源热泵系统因其造价低、换热功率高而被广泛使用［11］，过去针对地下水源热泵系统的研究大多是针对热量运移［12］，而关于地下水源热泵系统运行对地下水水质影响的研究较少，De Keuleneer等［13］曾就滨海地区利用地下水源热泵系统控制海水入侵，并就不同抽、灌井布井方式进行研究，但是没探讨地下水源热泵系统应用对浅部咸水入侵深部淡水的影响。

本次研究基于地下水渗流、热量运移和溶质运移理论，以江苏省地质工程勘察院（南通分院）科研实验楼地下水源热泵系统为研究对象，建立地下水渗流—热量运移—溶质运移多场耦合数值模型，在考虑合理科学开发浅层地热能的基础上，进一步评价浅层地热能开发对深层地下淡水咸化控制效应。

2 工程概况

科研实验楼位于南通市崇川区，分为主楼（6F）和附楼（2F），采用地下水源热泵系统进行制冷和供暖，根据设计指标，科研实验楼冬季最大负荷345 kW，夏季最大负荷465 k W。地下水源热泵系统工程布置了1口抽水井和1回灌井，受场地限制，抽水井和回灌井之间的距离为58 m。场地布置示意图见图1。

图1 江苏省地质工程勘查院（南通分院）地下水源热泵场地布置示意

2．1 沉积地层 南通市第四系厚度一般可达300 m，其水文地质体结构是砂土与黏性土互层组成的复合含水层系统。区内下更新统（Q1）埋深在190 m以下，厚度60 m左右。由中粗砂、粉细砂构成，为河湖相沉积，含水砂层构成区第Ⅲ承压含水层组；江苏省地质工程勘察院（南通分院）一带缺失中更新统（Q2）砂层，即研究区缺失第Ⅱ承压含水层；区内上更新统（Q3）埋深在40～120 m之间，厚度80～100 m。由中粗砂、粉细砂构成，受两次海浸影响，形成海陆交互相沉积，含水砂层构成区内第Ⅰ承压含水层组；区内全新统（Q4）厚度40～50 m，由粉细砂、亚砂土构成，含淤泥质，为三角洲海陆交互相沉积。上部多为长江冲积土层，属近现代沉积物［14］。

2．2 水文地质条件 工程场地的地下水以松散孔隙水为主，主要赋存在砂层中，根据沉积环境可划分为潜水、第Ⅰ、第Ⅲ承压含水层，各含水层的富水性好，但水质差异较大。含水层之间夹黏土层，分别第Ⅰ、第Ⅲ黏性土弱透水层，水文地质剖面图见图2。潜水主要赋存于全新统地层（Q4）中，含水层下段具有一定微承压特性，水位埋深1．5～2．0 m，水质微咸，TDS在800～1600 mg/L之间，地下水化学类型为HCO3·Cl-Ca·Mg和HCO3·Cl·SO4-Ca·Na，潜水与地表水有着密切的水力联系。第Ⅰ承压水主要赋存于上更新统（Q3）中，水位埋深在3 m左右，水质差，TDS在1100～1800 mg/L之间，水化学类型为Cl·HCO3-Ca·Na。第Ⅲ承压水主要赋存于下更新统（Q1）中，水质好，TDS在700 mg/L左右，水化学类型为HCO3·Cl-Na·Ca。

3 数值模型

3．1 概念模型 以工程场地为中心向四周各扩展1000 m作为计算边界，模型底部以第Ⅲ承压含水层底部作为底部边界。根据土体的性质，模型垂向上共剖分5层，分别为潜水含水层、第Ⅰ弱透水层、第Ⅰ承压含水层、第Ⅲ弱透水层和第Ⅲ承压含水层。由于计算边界与抽、灌井距离远，结合南通地区水文地质条件、渗透系数参考值和影响半径经验公式，地下水源热泵系统的运行过程中抽、灌地下水影响范围小于1000 m，地下水源热泵系统运行对数值模型边界影响很小。因此，四周边界处理为Dirichlet边界，即水力学上处理为定水头边界、热力学上处理为恒温边界、水化学上处理为恒定浓度边界；根据计算域水力学特征，模型顶面处理为自由面边界，底部处理为Neumann边界，即水力学上处理为隔水边界、热力学上处理为热通量边界、水化学上处理为零通量边界。

本项目利用第Ⅰ承压含水层的浅层地热能，需要抽取、回灌第Ⅰ承压含水层的地下水，由于各层存在水位差，各含水层之间存在水力联系，地下水流特征概化为三维非稳定流。地下水源热泵系统运行过程中，假定每个节点的温度同时代表地下水和土颗粒温度，忽略地下水密度随温度发生的变化。3．2 数学模型 各含水层之间的黏性土层作为弱透水层和含水层一起直接参与计算，地下水流为三维非稳定流。取坐标轴方向与各向异性介质主方向一致，地下水系统渗流控制方程如下［15］：

式中：Kxx、Kyy、Kzz为各方向水的渗透系数，m/d；h为地下水位，m；W为源汇项，1/d；Ss为含水层贮水率，1/m。

假设地下水和含水介质骨架的热动平衡是瞬时完成的，即含水介质骨架与周围地下水具有相同的温度，并忽略由于温度差引起水的密度不一样而引起的上下自然对流的影响。地下水系统热量运移控制方程如下［16-18］：

式中：λx、λy、λz为各方向水的热动力弥散系数，W/（m·℃）；T为岩土体温度，℃；cw、cs分别为水和土体的体积热容量，J/（m3·℃）；u为多孔介质中地下水平均流速，m/d；Qc为热的源汇项，J/（m3·d）；βL、βT分别为纵向和横向热弥散度，m；δij为Kronecker记号；ui、uj为地下水平均流速分量，m/d；λ′ij、λp分别为热机械弥散系数和热传导系数，W/（m·℃）。

根据多孔介质流体动力学，溶质迁移的过程中满足质量守恒定律，本次研究针对稀物质在水中迁移过程，可假定流体的密度为定值。取坐标轴方向与主渗透方向一致，则地下水溶质运移控制方程为［19-21］：

式中：Dxx、Dyy、Dzz为各方向的水动力弥散系数分量，m2/d；c为地下水中溶质的质量浓度，mg/L；n为孔隙度；I为溶质源汇项，mg/（L·d）；αL、αT分别为纵向和横向溶质弥散度，m；D′ij为机械弥散系数，m2/d；Dd为分子有效扩散系数，m2/d。

研究表明：当地下水流速较快时，机械弥散系数远大于分子有效扩散系数，因此忽略分子扩散系数的影响［22］。

地下水平均流速分量分别表示为：

根据概念模型，建立与之相应的初始条件和边界条件，所建立的数值模型如图3—4所示。

图3 模型剖分平面示意图

图4 三维数值模型示意图

3．3 模型识别、验证 以2017年6月作为模型计算的初始时刻，由于模拟区相对较小，同层初始水位、地下水化学浓度、地温梯度变化较小，模型潜水、第Ⅰ、第Ⅲ承压含水层的初始条件根据实测水位、地下水化学浓度、岩土体温度进行赋值，相邻黏性土弱透水层初始条件通过插值获取，模型各层地下水位初始值从上往下依次为0．87，0．8，0．73，-9．14及-19．01 m；Cl-浓度初始值从上往下依次为2570，2100，1630，879及128 mg/L；平均温度初始值从上往下依次为17．74，18．03，19．74，22．19及24．40℃。

在上述模型的基础上，利用抽水试验反演模型的水力学参数，利用2018年11月16日至2019年3月3日地下水源热泵系统试运行监测数据反演岩土体热物性参数，计算域土力学参数参考土工试验结果，水动力弥散系数参考以往经验系数［23-24］。

采用试错法来反演模型的水力学参数和热物性参数，将不同参数下计算的地下水位与岩土体温度与实测值对比，获取拟合精度最好的参数。对比结果显示：地下水位和岩土体温度的实测值与趋势相同且计算误差小，所建立的模型可信。各层土体参数见表1，地下水位、岩土体温度拟合结果见图5—8。

图5 抽水试验地下水位拟合对比图

表1 各层土体参数一览表

3．4 参数敏感性分析 系统特性X主要由n个因素x＝（x1，x2，…，xn）所决定，分别令各因素在其各自可能范围内变动，分析系统特性X偏离基准状态X*的趋势和程度，以确定这些因素的变化对系统特征的影响程度。若分析xi对特性X的影响，可令其余各参数基准值固定不变，而令xi在其可能的范围内变动，这时系统特性X表现为：

图6 G2观测井埋深55 m处温度拟合对比图

图7 G2观测井埋深90 m处温度拟合对比图

图8 G2观测井埋深125 m处温度拟合对比图

若xi的微小变化引起X的较大变化，表明X对xi敏感，此时xi是高敏感性参数；若xi的变化范围较大而X变化微小时，表明X对xi不敏感，此时xi是低敏感性参数。本次研究的参数灵敏度应为抽水井温度和第Ⅲ承压含水层Cl-浓度，由于本次研究建立的数值模型涉及到的参数较多，且具有不同的量纲。因此，在敏感性分析时定义参数的相对灵敏度来表征参数的敏感性［25］，公式表示如下：

表2 参数灵敏度计算结果一览表

由计算结果可知：对于抽水井热贯通影响纵向热弥散度最大，热容量次之。对于第Ⅲ承压含水层Cl-浓度影响纵向弥散度最大，孔隙度次之。本次研究纵向热弥散度、热容量参数通过现场试验反演获得，孔隙度根据室内土工试验获得，确保了模型的可靠性，而由于模型识别验证阶段是时间较短，第Ⅲ承压含水层Cl-浓度反应不明显，因此各层溶质弥散度参考前期研究成果，以确保模型能客观反映含水层Cl-浓度的变化特征。

4 模拟结果

4．1 原设计方案预测结果 根据地下水源热泵系统的设计运行方案，设计拟采用1组地下水源热泵换热机组，设计负荷见表3，按照应用需求，提供整个科研实验楼所需的夏天制冷、冬天制热的能量，确定工程夏季6—9月开机，共运行122 d；冬季12、1、2月开机，共运行90 d。

表3 科研实验楼工程冷热负荷一览表

科研实验楼设计利用温差为5℃，每天运行10 h，经式（15）和式（16）计算得到夏季循环水量为943．88 m3/d，冬季循环水量为403．88 m3/d，抽汲的地下水以100%进行回灌，利用识别、验证后的数值模型来模拟地下水源热泵系统运行10年第Ⅰ承压含水层温度和第Ⅲ承压含水层浓度变化趋势。计算结果见图9—12。

图9 夏季第Ⅰ承压含水层水位等值线图（单位：m）

图11 第Ⅰ含水层平均温度等值线图（单位：℃）

图12 第Ⅲ承压含水层Cl-浓度等值线图（单位：mg/L）

式中：Qcool、Qheat分别为夏季和冬季循环水量，m3/d；Pcool、Pheat分别为建筑的冷、热负荷，kW；Pinput-cool、Pinput-heat分别为地源热泵机组制冷和热输入功率，kW；ΔT为利用温差，℃；toperation为每天运行时间，h；cw为水的体积热容量，J/（m3·℃）。

若按照设计方案运行，10年后第Ⅰ承压含水层中抽水井处水体温度为20．59℃，上升了0．85℃；第Ⅲ承压含水层Cl-浓度最大为156．11 mg/L，增大了28．11 mg/L。计算结果显示，位于长江三角洲地区的南通市由于夏季炎热，冬季温和，夏季的制冷周期和制冷负荷高于冬季制热周期和制热负荷，按照既有方案运行地下水源热泵系统会出现热失衡问题，即出现轻度热贯通现象［26］；深部淡水水质咸化，即第Ⅲ承压淡水Cl-浓度升高，而由于第Ⅲ承压含水层上部弱透水层厚度较厚，Cl-浓度增大幅度有限。

4．2 改进方案预测结果 为了科学合理开发浅层地热能，控制地下水系统环境恶化，对江苏省地质工程勘察院（南通分院）科研实验楼地下水源热泵系统运行方案进行比选，即在抽、灌井位置和建筑冷、热负荷不变的基础上，通过调整利用温差和夏季灌采比来控制热贯通和深层地下淡水咸化，利用已建立的数值模型，模拟不同运行方案下的地下水位、岩土体温度、第Ⅲ承压含水层咸化发展变化趋势，分析改进方案对地下水环境恶化的控制效应。

以往研究显示，增大利用温差和减小灌采比有利于控制热贯通的发生［27］，本次研究采用增大利用温差和减小夏季灌采比的改进方案进行分析，将夏季灌采比调整为80%和60%，冬季灌采比仍维持100%，以尽量减小夏季和冬季注入含水层的热量平衡；同时将利用温差调整为7和9℃，以减小含水层循环水量，弱化抽、灌井之间的相互影响，评价地下水源热泵系统运行10年地下水环境的变化，具体方案及模拟预测结果见表4、图13和图14。

表4 改进运行方案和模拟预测结果一览表

图13 抽水井温度变化曲线

图14 第Ⅲ承压含水层Cl-浓度变化曲线

计算结果显示采用方案8（利用温差9℃、夏季灌采比60%）的方案抽水井温度和第Ⅲ承压含水层Cl-浓度上升幅度最小，在井间距、冷热负荷不变的条件下，将利用温差由5℃分别提高至7和9℃运行地下水源热泵系统10年，抽水井温度分别降低了0．06和0．13℃，第Ⅲ承压含水层Cl-浓度最大值分别降低了5．58和9．09 mg/L；将夏季灌采比分别从100%降低至80%、60%，运行地下水源热泵系统10年，与原运行方案相比较抽水井温度分别降低了0．29和0．53℃，第Ⅲ承压含水层Cl-浓度最大值分别降低了2．73和5．20 mg/L。

5 结论

以江苏省地质工程勘察院（南通分院）地下水源热泵系统为例，建立孔隙地下水渗流—热量运移—溶质运移多场耦合模型，研究地下水源热泵系统运行过程中地下水位、Cl-浓度和岩土体温度的变化趋势。获得的主要结论如下：（1）按照原设计方案运行10年，由于夏季冷负荷高于冬季热负荷，第Ⅰ承压含水层出现轻度热贯通现象；同时在水动力作用下，第Ⅰ承压含水层的咸水入侵深部第Ⅲ承压含水层淡水，Cl-浓度升高相对较小。（2）增大利用温差和减小夏季灌采比能有效缓解热贯通和深层淡水咸化的发展趋势，增大利用温差对深层地下淡水咸化的控制效应更为显著，而减小夏季灌采比对热贯通的控制效应更明显。（3）通过模拟计算可以发现，调整利用温差和夏季灌采比可实现地下水源热泵系统运行方案改进，缓解热贯通，对第Ⅲ承压含水层咸化也存在影响，改进运行方案的本质是通过调整循环水量改变地下水动力场，进而控制地下水系统环境的恶化。