南通市地下水压缩开采对水质咸化的控制效应

陆 玮，李 兆，骆祖江

(河海大学地球科学与工程学院，江苏 南京 211100)

南通市地下水资源较为丰富，一直以来以开采优质的第Ⅲ承压水来满足日常生产生活用水需求，但随着地下水的开采和地下水位的下降，从20世纪80年代起出现了水质咸化现象[1-4]，进入90年代后水质咸化现象日趋严重[5]。王琦等[6-7]通过对比分析地下水动力场和水化学场的变化特性，认为过量开采第Ⅲ承压水，导致上覆地层中的咸水入侵进入第Ⅲ承压水是咸化的主要原因。目前地下水开采对第Ⅲ承压水咸化影响的定量分析研究成果较少，难以为第Ⅲ承压水咸化的防控提供科学依据。本文以地下水渗流理论和溶质运移理论[8]为基础，建立了南通市地下水非稳定渗流与溶质运移三维耦合数值模型，通过模拟预测地下水现状开采和压缩开采条件下第Ⅲ承压水的咸化趋势，定量评价了南通市地下水压缩开采对第Ⅲ承压水咸化的控制效应，可为第Ⅲ承压水咸化的科学防控提供参考。

1 地质概况

南通市东临黄海，南望长江，与上海隔江相望，地理位置优越。研究区范围包括崇川区、港闸区和通州区，如图1所示。区域内除有少量基岩出露外，其余地区均为第四纪松散沉积物所覆盖，厚度一般为200～360 m，垂向上相互交错分布着多层砂层，构成了一个错综复杂的含水层系统[9-11]，从上至下依次分为：全新统的潜水含水层(组)、上更新统的第Ⅰ承压含水层(组)、中更新统的第Ⅱ承压含水层(组)、下更新统的第Ⅲ承压含水层(组)。图2为研究区水文地质剖面，各含水层之间以弱含水的黏性土层分隔或直接接触，存在着较为强烈的水力联系，贮存着丰富的地下水资源，其中潜水、第Ⅰ承压水和第Ⅱ承压水均为咸水，第Ⅲ承压水为淡水，水质好，水量大。

图2 研究区A—A′断面水文地质剖面示意图Fig.2 A-A′ hydrogeological profile of the study area

2 概念模型

本次模拟计算平面上以行政区域为界，垂向上包括整个第四纪地层，四周侧向边界均概化为通用水头边界和溶质运移通量边界，底部以第Ⅲ承压含水层底板为界，是一层厚度30～50 m的致密亚黏土，隔水性能良好，概化为隔水边界和溶质运移零通量边界。顶部在接受大气降水补给的同时，地下水又通过其蒸发排泄，在水力学上既是补给边界又是排泄边界，在溶质运移方面概化为通量边界；将大气降水入渗补给和地下水蒸发排泄两者合并考虑，概化为综合有效入渗补给强度。整个计算区含水层无论是平面上还是剖面上，富水性变化较大，概化为非均质层，垂直方向和水平方向的渗透性相差较大，概化为各向异性。各含水层中的地下水除水平运动外还通过层间的黏性土弱含水层或“天窗”发生水力联系，地下水流态为三维非稳定流，地下水开采按点井给出。

3 数学模型

3.1 地下水渗流模型

根据南通市的水文地质概念模型，将主渗透方向设为与坐标轴方向一致，建立了南通市地下水系统三维非稳定渗流数学模型[12-13]，模型的微分方程、初始条件和流量边界分别为

(1)

H(x,y,z,t)|t=0H0(x,y,z,t0) (x,y,z∈Ω)

(2)

(3)

(4)

式中：Kxx、Kyy、Kzz分别为含水层各向异性主渗透方向的渗透系数，m/d；H为点(x,y,z)在t时刻的水头值，m；W为源汇项，L/d；μs为储水率，m-1；t为时间，d；Ω为计算域；H(x,y,z,t)为渗流场水头，m；H0(x,y,z,t0)为初始水头值，m；q(x,y,z,t)为已知流量；cos(n,x)、cos(n,y)、cos(n,z)分别为流量边界外法线方向与坐标轴方向夹角的余弦；μ为饱和差或给水度；qw为自由面单位面积上的大气降雨入渗补给量。式(3)和式(4)分别为第二和第三类流量边界。

3.2 溶质运移模型

在不考虑溶质吸附、化学反应以及生物降解的情况下，将坐标轴方向设为与溶质运移方向一致，建立了南通市地下水系统三维非稳定流溶质运移数学模型，模型的微分方程[14]、初始条件和边界条件分别为

(5)

ρ(x,y,z,t)|t=0ρ0(x,y,z,t0) (x,y,z∈Ω)

(6)

x,y,z∈Γ2)

(7)

g(x,y,z,t) (t>0;x,y,z∈Γ3)

(8)

式中：ρ为流体中溶质的质量浓度，mg/L；n为孔隙度；I为源汇项；Dxx、Dyy、Dzz为沿坐标轴方向的水动力弥散系数分量，m2/d；De为分子有效扩散系数，m2/d；vx、vy、vz为介质中流体沿坐标轴方向的平均线性速度，m/d；αL、αT分别为纵向和横向弥散度，m；ρ(x,y,z,t)为溶质的质量浓度分布，mg/L；ρ0(x,y,z，t0)为已知溶质的质量浓度分布，mg/L；Ω代表整个模型区域；f(x,y,z,t)为Γ2边界上已知的弥散通量；cos(n,x)、cos(n,y)、cos(n,z)分别为边界Γ2外法线方向与坐标轴方向夹角的余弦；g(x,y,z,t)为已知函数，表示Γ3边界法线方向上的对流弥散总通量。式(7)和式(8)分别为第二和第三类边界条件。

4 模型的识别验证

采用地下水三维模拟软件Visual Modflow[15-17]进行模拟计算，在平面上剖分为159×100的网格单元，如图3(a)所示，其中有效单元为8 332个，无效单元为7 568个。剖面上从上至下分别将潜水、第Ⅰ承压、第Ⅱ承压、第Ⅲ承压含水层以及各含水层之间的黏性土弱含水层剖分成独立的层位，共分7层，图3(b)(c)分别为第35行和第43列网格剖分结果。选取2014年8月31日至2016年8月31日和2016年8月31日至2017年12月31日分别作为模型的识别和验证阶段，把每个月作为一个应力期，识别阶段共分为24个应力期，验证阶段共分为17个应力期，每个应力期为一个时间步长。

(a) 平面剖分

(b) 第35行剖分

图3 研究区网格剖分
Fig.3Gridmeshingofthesimulationfield

研究区分布有85口开采井，10口水位观测井，17口水质监测井，基本覆盖全区。各含水层的初始等水位线图由实测资料得到，各含水层之间黏性土弱含水层的初始等水位线图由其上下含水层插值得到。选取地下水中分布广泛且较为稳定的氯离子作为地下水咸化的模拟因子[18]，各含水层的氯离子初始质量浓度场由实际监测资料给出，各含水层之间黏性土弱含水层的氯离子初始质量浓度场由其上下含水层插值得到。参数的初始值根据前人资料结合经验给出[19]，经过识别、验证，模型共分为59个参数分区，其中潜水含水层4个，第Ⅰ黏性土弱含水层4个，第Ⅰ承压含水层8个，第Ⅱ黏性土弱含水层4个，第Ⅱ承压含水层5个，第Ⅲ黏性土弱含水层14个，第Ⅲ承压含水层20个。综合有效入渗补给系数为2×10-6。图4为第Ⅲ承压含水层水文地质参数分区，表1为第Ⅲ承压含水层各分区水文地质参数。图5和图6分别为识别阶段和验证阶段第Ⅲ承压含水层水位及水质监测井的拟合情况，可见实测值与计算值吻合良好，表明反演计算所得的参数准确、可靠，模型模拟结果符合实际情况。

图4 第Ⅲ承压含水层水文地质参数分区Fig.4 Hydrogeological parameter division for aquifer No.Ⅲ

表1 第Ⅲ承压含水层各分区水文地质参数
Table1HydrogeologicalparametersofeachzonesforaquiferNo.Ⅲ

分区渗透系数/(m·d-1)KxxKyyKzz有效孔隙度储水率/m-1弥散度/m纵向横向150.050.05.00.4400.00095070.7220.020.02.00.4430.00065080.8317.517.51.60.2960.000960212.1422.022.02.20.4320.00040080.8550.050.05.00.4380.00095080.8660.060.06.00.4420.00065080.8725.025.02.50.4440.00065080.888.08.00.80.4280.000400151.597.07.00.70.4420.00009590.9101.01.00.10.2060.000500323.21115.015.01.50.2850.000850252.51218.018.01.80.4330.00160090.91310.010.01.00.4400.00160090.91430.030.03.00.4410.00040090.91512.012.01.20.4410.00045090.91625.025.02.50.4480.000910121.21720.020.02.00.4420.00040090.9181.51.50.20.4410.00160090.9198.08.00.80.2000.000500353.52030.030.03.00.2350.000500292.9

(a) 识别阶段

(b) 验证阶段图5 地下水水位拟合结果Fig.5 Fitting graph of the groundwater level

(a) 识别阶段

(b) 验证阶段图6 地下水氯离子质量浓度拟合结果Fig.6 Fitting graph of the groundwater chloride ion concentration

5 地下水水位与咸化趋势预测

5.1 现状开采条件下地下水水位与咸化趋势

按2017年现状开采布局及开采量对2018—2034年的地下水水位变化和地下水咸化趋势进行模拟预测。现状共有85口开采井，其中第Ⅰ承压含水层分布有59口开采井，总开采量为444.21万m3/a；第Ⅲ承压含水层分布有26口开采井，总开采量为179.51万m3/a，其中港闸、崇川两区开采量为89.51万m3/a，通州区开采量为90.00万m3/a。降水量采用2007—2017年的年平均值，同时考虑区内蒸发作用，确定综合降雨入渗补给系数为2×10-6。图7和图8分别为2018年1月1日第Ⅲ承压含水层等水位线和氯离子质量浓度场。

图7 第Ⅲ承压含水层等水位线(单位：m)Fig.7 Flow field for aquifer No.Ⅲ(unit: m)

图8 第Ⅲ承压含水层氯离子质量浓度场(单位：mg/L)Fig.8 Chloride ion concentration for aquifer No.Ⅲ(unit: mg/L)

a. 地下水水位预测。图9为2024年和2034年12月31日第Ⅲ承压含水层预测等水位线，可以看出第Ⅲ承压含水层地下水水位整体上出现了回升现象。2018—2034年，最靠近降落漏斗的观山音街道附近水位回升最快，水位上升7.86 m，狼山镇街道水位回升4.12 m，其余地区水位也呈缓慢上升趋势。

b. 地下水咸化预测。图10为2024年和2034年12月31日第Ⅲ承压含水层预测氯离子质量浓度场，可以看出第Ⅲ承压含水层中氯离子质量浓度大于 250 mg/L 的面积从2018年初的147.35 km2增加到2034年底的 355.17 km2，其中质量浓度大于270 mg/L的面积达到30.67 km2。2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承压水咸化速率分别为19.31 km2/a、12.76 km2/a、11.49 km2/a和9.59 km2/a。

(a) 2024年12月31日

(b) 2034年12月31日图9 现状开采条件下第Ⅲ承压含水层预测等水位线(单位：m)Fig.9 Predicted flow field for aquifer No.Ⅲ under current mining conditions(unit: m)

(a) 2024年12月31日

(b) 2034年12月31日图10 现状开采条件下第Ⅲ承压含水层预测氯离子质量浓度场(单位：mg/L)Fig.10 Predicted chloride ion concentration for aquifer No.Ⅲ under current mining conditions(unit: mg/L)

5.2 压缩开采条件下地下水水位与咸化趋势

南通市地下水压缩开采方案中共分为两个阶段，第一阶段为2018—2019年，研究区分布有55口开采井，其中第Ⅰ承压含水层分布有31口开采井，总开采量为352.05万m3/a，相比现状开采减少92.16万m3/a；第Ⅲ承压含水层分布有24口开采井，总开采量为154.60万m3/a，相比现状开采减少24.91万m3/a，其中港闸、崇川两区开采量为 74.50万m3/a，通州区开采量为80.10万m3/a。第二阶段为2020—2034年，在第一阶段基础上不改变开采井的分布，继续减少每口井的开采量，其中第Ⅰ承压含水层总开采量减少为301.35万m3/a，相比现状开采减少142.86万m3/a；第Ⅲ承压含水层总开采量减少为108.90万m3/a，相比现状开采减少70.61万m3/a，其中港闸、崇川两区开采量为 55.80万m3/a，通州区开采量为53.10万m3/a。降雨入渗量同现状开采条件。

a. 地下水水位预测。图11为2024和2034年12月31日第Ⅲ承压含水层预测等水位线，可见第Ⅲ承压含水层地下水水位整体上是回升的。从2018—2034年，最靠近降落漏斗的观山音街道附近水位回升最快，水位上升8.00 m，相比现状开采水位上升了0.14 m，狼山镇街道水位回升4.31 m，相比现状开采水位上升0.19 m。

b. 地下水咸化预测。图12为压缩开采条件下2024和2034年12月31日第Ⅲ承压含水层预测氯离子质量浓度场，可以看出第Ⅲ承压含水层中氯离子质量浓度大于250 mg/L的面积从2018年初的147.35 km2增加到2034年底的329.21 km2，其中质量浓度大于 270 mg/L 的面积为1.76 km2。2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承压水咸化速率分别为17.50 km2/a、11.92 km2/a、9.93 km2/a和7.52 km2/a。

由以上分析可知，现状和压缩开采两种条件下，第Ⅲ承压含水层地下水位均逐渐上升，与上覆含水层之间的水头差均逐渐减小，因此上覆含水层的咸水入侵速度均逐渐减缓。相比于现状开采，压缩开采条件下第Ⅲ承压水咸化速率明显减慢，且咸化面积明显减少。由此可以判断，采用压缩开采方案能有效控制南通市第Ⅲ承压水的咸化问题。

(a) 2024年12月31日

(b) 2034年12月31日图11 压缩开采条件下第Ⅲ承压含水层预测等水位线(单位：m)Fig.11 Predicted flow field for aquifer No.Ⅲ under compression mining conditions(unit: m)

(a) 2024年12月31日

(b) 2034年12月31日图12 压缩开采条件下第Ⅲ承压含水层预测氯离子质量浓度场(单位：mg/L)Fig.12 Predicted chloride ion concentration for aquifer No.Ⅲ under compression mining conditions(unit: mg/L)

6 结 论

a. 保持现状开采条件下，2034年底第Ⅲ承压含水层中氯离子质量浓度大于250 mg/L和270 mg/L的面积达到355.17 km2和30.67 km2，2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承压水咸化速率分别为19.31 km2/a、12.76 km2/a、11.49 km2/a和9.59 km2/a。

b. 压缩开采条件下，2034年底第Ⅲ承压含水层中氯离子质量浓度大于250 mg/L和270 mg/L的面积分别达到329.21 km2和1.76 km2，2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承压水咸化速率分别为17.50 km2/a、11.92 km2/a、9.93 km2/a和7.52 km2/a。

c. 压缩开采条件下，相比于现状开采条件2018—2034年第Ⅲ承压含水层中氯离子质量浓度大于250 mg/L和270 mg/L的面积分别减少了25.96 km2和28.91 km2，咸化速率明显减慢，采用压缩开采方案能有效控制南通市第Ⅲ承压水的咸化问题。