河道硬化对傍河地下水源补给结构及范围的影响

刘 昶， 郑明霞， 孙源媛， 苏 婧， 傅雪梅， 丁鸿羽， 席北斗*

1.南昌大学资源环境与化工学院, 鄱阳湖环境与资源利用教育部重点实验室， 江西 南昌 330031

2.中国环境科学研究院, 国家环境保护地下水污染模拟与控制重点实验室， 北京 100012

傍河地下水源因水量充足、水质较好等优点成为我国常见的水源类型. 然而，国内河道整治经常采用护岸与河床固化的方式，形成浆砌三面光河道来满足防洪需求或达到减少河道渗透损失量与河水污染物下渗的目的. 但这一治理方式往往缺乏对河流生态环境的综合考虑[1-2]，河道硬化阻碍了地表水与地下水的相互联系，削减了河水对地下水的补给，降低了傍河地下水的自净能力[3]，使傍河地下水的水量水质受到较大影响，危及地下水安全与可持续发展[4-6].

近年来，河流与地下水之间的交互作用机理和机制逐渐成为研究热点. 冯创业等[7-8]在河流流域内开展了基础性大型入渗试验，采用水量平衡原理计算出河流对地下水的入渗量与河道的入渗能力. 任长江等[9-10]通过流域内不同植被覆盖条件下的地表入渗试验，拟合不同入渗模型，分析总结了不同植被覆盖条件下的土壤入渗规律及地表水对地下水的补给. Marie-Jose等[11]将热红外遥感与直接超孔隙度测量相结合，精确捕捉了一条40 km河段上河水与傍河地下水的交互作用，进一步揭示了地表水与地下水交互作用在水环境系统上的重要性. 研究地表入渗的方法越来越丰富且先进，数值模拟法成为目前最为精确且可靠的研究方法之一. Singh等[12-17]通过建立地下水模型可直观而精确地揭示地表水入渗地下水的状况. Frei等[18]通用耦合河流-含水层的模型研究其交互动力学和交互模式，模型中分析并量化了地表河流对地下水的补给关系，揭示了河水与地下水的交互规律. XIE等[19]在空间上将河水与地下水的渗流作用分为垂向渗流与横向渗流，利用热量作为定量交换示踪剂，结合模型研究河流-含水层的渗流作用. 虽然地表水与地下水之间的交互对于自然界水循环具有非常重要的作用，但由于缺乏对河道硬化造成的地表水入渗量及其对地下水影响的系统、量化的研究，河道硬化在我国河道整治工作中仍然是常见的工程措施，该措施不仅会减少河水对地下水的入渗量，还会使得地下水位下降，改变地下水的补给范围与补给结构，对地下水安全产生威胁[20-24].

以河北省张家口市Y傍河地下水源地为研究对象，建立地下水数值模型，模拟研究区地下水流场，对研究区地下水进行水均衡分析，并运用质点追踪技术对水源地井群的补给区进行刻画. 通过对比河道硬化前后模型计算结果，分析并量化了河道硬化对Y傍河水源地的水均衡状况、补给结构及其补给范围的影响，结合历史与近期水质数据对比分析河道硬化对地下水质条件的影响，以期为相关地下水保护提供参考和依据.

1 材料与方法

1.1 研究区概况

1.1.1自然地理

Y傍河地下水源地位于张家口盆地，其地理位置与水源井分布如图1所示. 张家口盆地为典型的断裂型山间河谷盆地，东西长约90 km，南北宽约21 km，洋河自西向东贯穿其中，研究区是盆地中心的洋河带状平原，周围低山环境，整个盆地的地势由西向东倾斜，坡降约4%，海拔在550～1 000 m范围内. 研究区属北寒温带大陆性半干旱气候区，2004—2018年平均降水量为394 mm，最大年降水量为496 mm(2010年)，最小年降水量为323 mm (2006年)，多年平均气温为5～8 ℃，最低气温为-26 ℃，最高气温为42 ℃.

图1 Y傍河地下水源地地理位置与水源井分布

1.1.2水文地质

研究区属山间盆地水文地质区，主要分布山地基岩裂隙水和平原区第四系孔隙地下水. 由图2可见，盆地腹部的河流冲积平原第四系厚度可达100～300 m，地下水埋深为3～60 m，含水层岩性以中粗砂、砂砾石为主，单井涌水量10～120 m3h，为HCO3-CaMg型水，历次水文地质勘察得出的含水层渗透系数在1～150 md，近期通过抽水试验得到渗透系数为7.36～163.87 md，孔隙度为0.25～0.32. Y傍河地下水源地处于盆地第四系厚度最厚区域，由3个含水层组和3个不连续的隔水层组成，下部略具承压性，其开采层位主要接受来自上游冲洪积扇地下水的侧向流入补给，其排泄方式以开采和向下游径流为主，地下水流方向大体为自西向东.

图2 研究区水文地质剖面

1.1.3水源地基本情况

Y傍河地下水源地为傍洋河冲积平原水源地，位于张家口盆地洋河带状冲积平原东段的西端，地势开阔平坦. 水源井群共有生产井18眼，沿洋河北岸呈东西向分布两列，距离河岸100 m左右(见图1). 该水源地生产井深96～120 m，含水层由第四系全新统冲积砂砾石孔隙含水层和第四系全新统冲洪积砂砾卵石孔隙含水层构成，补给条件良好，主要接受洋河冲积平原、洋河河水的地表水系补给与万全区城东河和城西河冲洪积扇的地下水补给. 根据建设综合勘察研究设计院有限公司编写的《河北省张家口盆地地下水详查报告》，其最大供水能力为150×103m3d，根据水厂供水量记录计算得到实际平均取水量为54×103m3d，用水高峰期取水量为90×103m3d. 2010年张家口市政府开展洋河综合治理工程，为增强洋河防洪能力，降低河水对地下水的污染，对河道浇筑混凝土进行防渗处理，河道形成了“三面光”状态.

1.2 研究方法

该研究中主要研究方法为数值模拟法，选用美国盐湖城大学环境研究所开发的地下水模拟软件GMS (Groundwater Modeling System)[25-26]来建立研究区地下水数值模型，通过对模型参数和源汇项进行识别校正，使得数学模型能如实反映研究区水文地质条件，模拟的地下水流场能有效拟合实际水位情况，从而揭示研究区地下水的水均衡状况. 在此基础上，使用软件模块MODPATH，对水源井群进行反向质点追踪，通过示踪粒子在地下水中向前 1 000 d的运移轨迹[27]，得到水源井群抽水的捕获区域，即地下水源的补给范围.

地下水流动的数学模型是基于质量守恒定律而建立的，水流控制偏微分方程：

(1)

式中：h为地下水水头，m；Kxx、Kyy和Kzz分别为含水层在x、y、z方向上的渗透系数，md；W为单位体积流量，代表流进汇或来自源的水量，d-1；Ss为含水层孔隙介质的贮水率，m-1；t为时间，d.

1.3 概念模型

1.3.1边界概化

确定模型范围和划分边界时尽可能以自然边界或地貌单元分界线作为模型边界. 该模型南部和北部边界为山区和平原区的交界地带，接受山前地下水侧向径流补给，定义为定流量边界；东部和西部边界分别为盆地地下水流入和流出边界，均定义为定水头边界；地表河流洋河定义为河流边界. 模型范围及边界条件如图3所示.

图3 模型边界条件示意

根据研究区实地勘探结果并结合张家口盆地水文地质资料对模型边界进行确定. 该模型定义西部边界水头值为750 m，东部边界水头值为590 m；根据地下水位数据，应用达西定律计算得到，南北边界侧向流量数值，北部定流量边界划分为三部分，取值分别为5.00×104、10.00×104、15.00×104m3d，南部边界划分为两部分，取值分别为0.45×104、0.30×104m3d.

1.3.2模型结构概化

研究区地下水以水平运动为主，流速缓慢，渗流符合达西定律，其介质是非均质的. 在模拟过程中，岩性相似的相邻含水层通常合并为一个含水层[28]. 据勘探资料及水文地质剖面图，模型最终概化为2个含水层组的非均质各向异性三维稳定流模型. 根据场地高程DEM图和钻孔资料获取各层顶底板标高. 对研究区渗透系数进行分区，如图4所示，其中a区、b区、c区、d区、e区、f区、g区的渗透系数分别赋值为60、6、55、40、120、140、8 md. 根据岩层属性，垂向渗透系数设为水平渗透系数的14，给水度、储水率、有效孔隙度等参数值取对应地层属性的经验值[29-30].

注： a、b、c、d、e、f、g为对应的渗透系数分区编号.

1.3.3源汇项

降雨入渗量和灌溉入渗量结合蒸发量通过面状补排程序包RCH处理计算. 研究区极限蒸发深度为4 m，根据2019年地下水水位统测结果，计算出地下水位埋深小于4 m的面积为9.49 km2. 模型将降雨补给与入渗补给系数整合，依据张家口市气象资料与研究区水文地质条件，对研究区蒸散发及补给参数进行分区，如图5所示，其中A区、B区、C区、D区、E区、F区的蒸散发及补给参数分别赋值为-0.005 0、0.000 2、0.001 0、0.002 0、0.003 0、0.000 2 md.

注： A、B、C、D、E、F为对应的蒸散发及补给系数分区编号.

河流入渗量通过河流程序包RIV处理计算. 考虑到洋河河床河堤进行过硬化处理，而河流参数主要为河床的水力传导系数，由底积层的岩性、厚度及河床宽度、计算单元的河段长度决定[31]，计算公式：

C=KLWM

(2)

式中：K为河床土壤水力传导系数，md；L为河流映射长度，代表河流在模型中运算的有效长度，m；W为河流宽度，m；M为河床厚度，m.

根据现场地质试验，将洋河分为4个河段赋值，各河段分别赋值为0.5、0.2、0.3、0.1 m2(d·m)，该赋值范围远小于天然条件下的河流水力传导系数，以期刻画河道硬化的实际情况.

水源井开采量通过井流程序包WEL处理计算. 根据张家口市市政水务2014—2018年Y水厂供水量记录，计算得到18口水源井实际的平均单井取水量为3.0×103td，用水高峰期的平均单井取水量为4.7×103td. 考虑后续城市人口增长造成用水量增加，模型内每口井抽水量取值为5.0×103td.

1.4 数学模型求解与参数灵敏度分析

该研究利用GMS里的MODFLOW程序建立了地下水流动模型，采用有限差分的离散方法剖分网格，模拟范围东西方向最长为49.4 km，南北方向最长为25.1 km. 研究区剖分为132行、260列、2层，共 32 014 个活动单元，网格大小为190 m×190 m. 模拟期为2019年2月. 解算器选择稳定性高、收敛性好的共轭梯度解算器PCG2. 完成地下水数值模拟的模型构建和数据输入后，将概念模型映射到MODFLOW数值模型中. 检查无误后保存模型并运行.

通过模型参数的灵敏度分析明确哪些参数对模型输出结果和动态影响较大，进而减小因模型参数的不确定性产生的误差. 对模型中含水层渗透系数、降雨入渗补给系数、河床水力传导系数、定流量边界水力传导系数、河水水位和含水层孔隙度6个主要参数进行灵敏度分析，模拟水源地捕获区最远距离的变幅作为灵敏度分析的依据. 选取各参数分别调整±5%和±10% 对数学模型进行灵敏度分析. 由图6可见，含水层渗透系数与河床水力传导系数是对捕获区最远距离影响的主要敏感因子，其他4个参数影响则较弱，表明Y傍河地下水源地的补给范围主要受到含水层属性与河流入渗的影响. 由于河床水力传导系数是模型研究河道硬化的调整因素，因此含水层渗透系数的获取与赋值大幅影响了模型结果的精度. 通过增加研究区野外试验，获取更密的含水层渗透系数分布，再通过克里金插值法对含水层渗透系数进行更细致合理的赋值，以减小因参数不确定性对模型运行结果可靠性产生的影响.

图6 灵敏度分析结果

2 结果与讨论

2.1 模型识别与验证

水源地周边共有25个水位观测井，其水位数据可较为清晰地描绘水源地的地下水水位状况. 先将水位数据绘制等水位线，将模拟流场等水位线与实测等水位线比对，结果如图7所示.

图7 模拟流场与实测等水位线对比

将25个水位观测数据导入模型进行校准. 根据校核目标的填充颜色查看模型校正情况，当模型校正目标颜色条变小且为绿色时，可认为模型各项参数调至最优. 由图8可见：在拟合效果图中校核目标的条形有17个点为绿色，8个点为黄色，未出现红色点位，说明观测值与计算值的差值在校核置信范围区间内；计算水位拟合误差平均值为0.48 m，模型拟合效果达到要求，水流模拟结果可靠，能较好地反映地下水情况.

注： 图(a)中绿色或黄色点位为观测井，不同颜色表示水位拟合校核结果，其左侧误差条位于点上方时表示计算水位高于实测水位，反之则表示计算水位低于实测水位. 当观测值与计算值差值在校核置信范围(0～100%)时，用绿色充填；当误差超出置信区间范围(100%)但小于200%时，则用黄色充填；当误差大于200%时，则用红色充填.

2.2 模型运行结果

通过对模型进行合理调参，参数修正并不断迭代[32]，最终使拟合效果达到要求，水流模拟结果可靠，得到能较好地反映地下水情况的模型. 由图9可见，运行结果可以观察到研究区地下水流场情况，其主要受地形与地层构造控制，地下水流向主要为由西向东，流向在盆地东侧转向东南方向，地下水水头在580～780 m之间.

注： X、Y、Z分别表示空间坐标轴正东向、正北向和垂直水平向上方向.

2.3 河道硬化对傍河水源地补给区范围影响分析

为了研究水源井群取水的补给范围，应用MODPATH模块模拟计算出水源井群抽水的补给区. 基于已建好的MODFLOW水流模型进行质点追踪.

对于无硬化状态下的模拟，通过调整河流参数来模拟分析，主要调整河流传导系数(河床土壤饱和导水率)与河水水位. 在无硬化处理状态下，河水水位取值参考硬化处理前的历史资料，河流传导系数取值参考实地测量参数结果并代入式(2)计算获得，使河流传导系数取值更接近未硬化处理的自然河流，将河流传导系数设定为20～30 m2(d·m).

研究区含水层介质主要为砂卵砾石，孔隙度取值为0.3. 该模型对每口井设置示踪粒子16个，分别对硬化和无硬化的情况进行 1 000 d时长捕获模拟，结果如图10所示.

图10 硬化与无硬化状态下数值模拟水源地1 000 d补给区对比

由图10可见：无硬化状态下 1 000 d补给区上游部分偏向河流，为沿河流分布的狭长型区域，最长距离为 5 062 m，最宽距离为 1 002 m，面积为4.74 km2；硬化状态下 1 000 d补给区较无硬化状态下整体向远离河流方向偏离，其宽度更宽，呈上游窄下游宽的“喇叭形”，沿河流方向上最长距离减少了271 m，垂直河流方向的最宽距离增加了210 m，面积增加了0.77 km2. 根据不同状态下等水位线对比发现，硬化状态时，补给区地下水水位下降了2～6 m.

河道经硬化处理后，水源井群补给范围发生了改变，河水与地下水水力联系变弱造成了其补给范围向远离河流的方向偏移，来自北部上游地下水的补给比例增大，其补给范围的面积也进一步扩大，给水源地的用水安全带来新的潜在风险. 观察等水位线变化可知，硬化处理后河水对河流附近地下水补给作用减弱；河流北岸的等水位线方向也受到改变，北岸地下水流向向东南方向偏移；水源井群上游等水位线变密，下游等水位线变疏，说明井群上游地下水流速变快，下游地下水流速变慢，在地下水水位下降的情况下，水源地取水对当地地下水流场的扰动影响变得更为强烈.

2.4 河道硬化对傍河水源地水均衡影响分析

河道硬化条件下，对井群 1 000 d补给区域进行水均衡分析. 由表1可见：井群 1 000 d补给区地下水系统的总补给量为251.00×103m3d，其中，上游边界流入补给量为主要补给来源，占总补给量的82.88%；降雨入渗补给量与河流入渗补给量则分别占总补给量的14.31%和2.81%. 井群 1 000 d补给区地下水系统的总排泄量为251.00×103m3d，下游边界流出排泄量与地下水开采量为主要排泄项，共计占总排泄量的98.61%；蒸散发排泄量占总排泄量的1.39%；河流排泄量为0. 井群 1 000 d补给区地下水系统深层对浅层含水层的越流水量为3.48×103m3d，浅层对深层含水层越流水量为3.90×103m3d.

水均衡分析(见表1)表明，水源地地下水的主要补给来源于上游地下水的流入，而河流入渗补给量与降雨入渗补给较少；地下水的去向主要为下游边界流出和地下水开采，地下水未对河流进行排泄. 相比于无硬化处理的情况，硬化处理状态下的河水对地下水的河流入渗流量由53.85×103m3d减至7.05×103m3d，减幅为86.91%；地下水对河水的河流排泄量由13.31×103m3d减至0，减幅为100.00%，地下水停止向河水排泄. 此外，上游边界流入补给量减少了2.07×103m3d，减幅为0.99%；下游流出排泄量减少了35.56×103m3d，减幅为18.95%.

表1 不同条件下井群1 000 d补给区地下水均衡表

研究显示，河道硬化不仅直接导致河流入渗量大幅减少，还间接使地下水的上游边界流入量与下游边界流出量随之减少；地下水位的下降使地下水静水位低于河床高程，导致地下水不再对河流进行排泄.

2.5 河道硬化对傍河水源地水质影响分析

研究区地下水水化学类型以HCO3-Ca型为主，河水pH、硬度及氨氮浓度均明显高于地下水. 地下水水质例行监测数据表明，河道硬化前研究区傍河地下水pH为7.8～8.1，总硬度为465～588 mgL(以CaCO3计)，氨氮浓度为0.02～0.30 mgL. 实地采样测试发现，河道硬化后近年来地下水pH为6.73～7.82，属中性水；矿化度为300～790 mgL，属低、中矿化度水；总硬度为295～325 mgL(以CaCO3计)，属软水或稍硬水；氨氮浓度在河道硬化后均未检出，这可能由于研究区地下水在河道硬化后切断了河水向地下水的入渗途径，使得河水中的污染物质，如氨氮等不能进入地下水，说明河道硬化对于防止河水对地下水的污染起到了一定的阻断作用，使地下水水质与理化环境发生部分改变；另一方面，河道硬化后地下水位不断下降，溶解氧与氧化还原电位监测数据均出现下降趋势，表明研究区地下水化学环境向还原性方向演化.

因此，河道硬化虽然能有效减少地表水污染物的入渗，但同时也使地下水位下降，地下水环境发生改变，进而可能导致地下水自净能力降低[33-34]，地下水脆弱性随之增加，导致地下水环境生态系统的抵抗力稳定性变弱[2,24]，对地下水源的用水安全带来诸多不利因素与潜在风险.

3 结论

a) 通过数值模型法模拟分析河道硬化条件下的傍河地下水源，结果表明，水源井群 1 000 d的补给范围呈上游窄下游宽的“喇叭形”，其最长距离 4 791 m，最宽距离 1 212 m，面积为5.51 km2；对补给范围进行水均衡计算，稳态下地下水总补给量与总排泄量均为251.00×103m3d (不包括层间越流量)，主要的补给项为边界流入与降雨入渗的补给，分别为208.04×103与35.91×103m3d，占比分别为82.88%与14.31%；主要的排泄项为边界流出与地下水开采的排泄，分别为152.12×103与95.40×103m3d，占比分别为60.60%与38.01%.

b)通过模型模拟河床无硬化的自然状态，分析河道硬化造成的影响. 对比发现，河道硬化使得其地下水总补给量与总排泄量减少了48.87×103m3d，河水对地下水的入渗流量减少了46.79×103m3d，入渗量减幅达86.91%，并且地下水停止了对河水的排泄；河道硬化同时使得水源井群 1 000 d补给范围在沿河流方向上的最长距离减少了271 m，垂直河流方向上的最宽距离增加了210 m，面积增加了0.77 km2，补给区向远离河流方向上偏移，相比自然河道条件下，补给区地下水水位下降了2～6 m.

c) 河道的硬化处理会改变河水与地下水的交互关系，使得傍河地下水的水力结构与补排条件发生改变，其傍河地下水源井群的补给范围发生偏移，地下水水位下降，井群上游地下水流速变快，下游地下水流速变慢，水源地取水对当地地下水流场的影响变得更为强烈；同时其地下水水质与理化环境均发生改变，傍河地下水自净能力降低，地下水脆弱性随之增加，给傍河地下水带来新的潜在风险.