河南周口龙湖湿地地下水-地表水交互作用模拟及修复决策研究

杨海瑞，徐永新，2，马 念，陈 杰，陶 亮，张志敏，刘 恋

(1.林同棪国际工程咨询(中国)有限公司，重庆 401121；

2.西开普大学 自然科学学院地球科学系，南非 开普敦 8000)

湖泊湿地是位于陆生生态系统和水生生态系统之间的过渡性地带，湿地通过水文循环实现与周边环境系统的物质和能量交换，水文过程制约着湿地环境的生物、物理和化学特征，从而控制着湿地景观格局、生物生长及其他生态功能的形成与演化[1]。

湖泊与地下水的交互作用对湖泊水文循环起重要作用，这种交互作用主要通过地表水与地下水之间的水量交换和水质演变来体现[2-3]。由于这种交互作用，使得地下水成为很多内陆湖泊湿地重要的补给来源。因此，探索地下水与湖泊水量、水质转化机制是理解湖泊湿地水循环机制及其物质、能量与信息传递的关键基础[4-5]，也是水环境治理和生态修复工程决策的前提。

随着全球气候变化和人类活动双重影响的加剧，湖泊湿地地表水-地下水系统的相互作用响应敏感，湿地水文过程更加复杂，湿地生态系统表现出高度脆弱性。在此背景下，地下水对于维持生态系统健康、结构与功能稳定的重要性日益得到重视，关于湿地-地下水交互作用的多学科交叉探讨成为当前国际湿地研究的热点与前沿课题[6]。城市湖泊受人类活动影响更大，由于城市化对湖泊产生负面影响的案例并不罕见，湖泊治理需要对水文关系有非常充分的认识，但湖泊和地下水的相互作用却并不像地表水那样容易被认知，因此，地下水的作用常常被忽视，或者把地下水和地表水分开管理[7-8]。特别是在城市湖泊水环境治理与生态修复工程中，往往较少关注地下水与湖泊水体水量平衡的关系，而相关地下水模型在城市水工程中的应用更是较为少见。

1.1 研究区

龙湖是位于河南省周口市淮阳区的城市湖泊，为典型的中原地区内陆湖泊型湿地，也是中原地区少有的重要城市景观湖。龙湖湿地地处黄淮平原腹地，湖泊东西宽4.4 km，南北长2.5 km，围堤长14 km，面积11 km2，水域面积5.34 km2，通常水深1.5～2.3 m，最大水深超过3 m。

龙湖湿地属淮河流域，主要流经水系为沙颍河水系。龙湖周边河流(渠)共有9条，分别是：淮郑河、七里河、白楼渠、古蔡河、北关沟、新蔡河、南关沟、徐老庄沟等(图1)。其中，淮郑河在龙湖西北侧与龙湖直接相连，曾是龙湖的补给水源来源；
南关沟在龙湖东南侧，与龙湖相通，作为退水通道。近年来，为防止外源污染负荷输入，淮郑河、南关沟和龙湖之间的进出水闸常年关闭，龙湖与周边地表水系现无直接补给关系。

图1 龙湖与周边水系示意

研究区埋深小于50 m的地层划为浅层含水层，主要由全新统、上更新统地层组成，上更新统底板作为浅层含水岩组底板。岩性主要由中细砂、细砂、粉砂、粉细砂、粉土及粉质黏土等组成。降水入渗补给系数一般为0.15～0.25。渗透系数为1～30 m/d，富水性较好，单井涌水量为1 000～3 000 m3/d。

天然条件下，浅层地下水主要接受降水入渗补给，其次为河流等地表水侧渗、渠道侧渗、灌溉回渗。地下水水位埋深2～6 m，总体流向由西北流向东南，水力坡度0.2‰～0.7‰。排泄方式主要为蒸发、人工开采、径流排泄及河流湖泊排泄。

1.2 研究方法

目前，针对湿地地表水-地下水交互作用的研究大多依赖地表水和地下水模型[9-10]。MODFLOW是由美国地质调查局开发出的一套专门用于孔隙介质中三维有限差分地下水流数值模拟的软件[11]。加拿大Waterloo水文地质公司在此基础上，应用现代可视化技术开发研制出Visual MODFLOW[12]，该模型具有合理的菜单结构、友好的界面、功能强大的可视化特征[13]，可作为本研究定量分析的软件支撑。

地表水和地下水的交互作用，不但形成水动力梯度，也形成化学梯度。在龙湖及周边采集水样进行常规离子检测，在水量平衡分析结果上，应用水文地球化学交互模型分析龙湖湿地与地下水之间的水化学联系及潜在风险。

本研究通过建立数值模型对水量平衡做定量计算，为水量、水质风险分析讨论奠定基础，重点介绍数值模拟的方法和过程。

1.3 数值模拟

1.3.1 模型概化

为建立模型，将地下水研究范围扩展到龙湖周边以新运河、周商运河、老黑河、狼牙沟、新蔡河和范丹寺沟为边界包围的区域。这是一个由河流切割的天然地下水系统。模型研究区坐标范围：大地坐标X=20 286 000～20 318 000 m，Y=3 720 000～3 752 000 m，平面面积7.92 km2，底部边界取在地表以下50 m处，顶部边界取在地表。地下水含水层主要为第四系全新统冲积砂层，含水层厚度较大，含水介质为非均质各向同性。研究区水文地质模型概化示意见图2。

图2中，龙湖西北侧新运河、周商运河和老黑河均为季节性河流，概化为第三类河流边界(A2)；
而东南侧狼牙沟、新蔡河和范丹寺沟所处地势较低，对地下水起排泄作用，故概化为第三类排水沟边界(A3)；
龙湖概化为第三类边界(A3)；
研究区南侧至新运河段概化为侧向径流边界(A1)；
研究区顶面为潜水面边界(A4)；
底面为上更新世黏土、亚黏土层，渗透性弱，构成区域隔水底板(A5)。

图2 研究区水文地质概念模型

1.3.2 数值模型建立

在建立水文地质概念模型的基础上，根据有限差分原理，利用Visual MODFLOW建立地下水流数值模型。有限差分原理的基础是时间和空间的离散化。其中，在水平面上采用间距为200 m等间距正交网格将研究区剖分为160行、160列，共剖分出了25 600个单元。从垂向上将研究区划分为1个含水层，为第四系冲洪积砂砾卵石孔隙含水层。

计算采用稳定流建模。根据计算流场形态、水位，并结合地下水监测井所测得的地下水位，对模型进行校正，使模型计算流场与区域地下水运动特征基本一致。

1.3.3 数据来源

本研究收集到的数据包括淮阳县降水序列(1999～2016年)、蒸发序列(1971～2013年)、地下水水位序列(《河南省地下水资料<周口市>》<1995～2016年>，周口市水文水资源局)、淮阳区浅层地下水供水量(《周口市水资源公报》<2008>，周口市水务局)。

实测数据包括：① 龙湖周边约8 km2范围地下水水位统测(2018年5～6月)；
② 龙湖西北侧和东南侧分别选取1眼水井进行了抽水试验，用来校正含水层渗透系数；
③ 2018年5月在龙湖周边7个井位采集了水样，进行了地下水常规离子检测分析。各测井位置见图3，水质采样点见图4。

图3 地下水监测井位分布

图4 水质分析取样点及矿化度值(单位：mg/L)

1.3.4 模型率定

根据水文地质条件，将潜水含水层划分为3个参数分区。地下水水流模型中水文地质参数渗透系数K和给水度μ值主要根据区域环境水文地质调查过程中试验给出的结果，并结合岩性特征和经验值给定初始值，通过模型率定，获得最终水文地质参数。

通过稳定流模拟计算，模拟区地下水稳定流场，见图5。区内地下水流向由西北流向东南，与区域地下水流向一致，与地下水流动特征基本相符。对比实测井孔地下水位与计算地下水位，两者基本一致，见图6，平均误差小于1 m，基本满足模型精度要求，证明所建模型基本可靠。经模型率定，研究区渗透系数取值8～20 m/d，分区见图7。

2.1 湿地地表水-地下水交互作用模式

水头差的存在是水交互作用产生的前提，水力梯度决定了源汇项分布及其作用路径，进而表现为不同的交互作用模式。Jolly等[14]将接触带尺度下湿地划分为4种模式：① 非饱和流-补给型湿地，湿地下垫面与地下水面之间存在不相连的非饱和区间，湿地地表水垂向渗流补给地下水，多见于季节性湿地系统；
② 饱和流-补给型湿地，湿地下垫面与含水层直接连通且湿地水位高于周边地下水，湿地水体因而成为周边地下水的补给来源；
③ 饱和流-排泄型湿地，与②水力梯度相反，四周地下水补给湿地；
④ 饱和流-贯穿型湿地，地下水流场的水力梯度方向连续一致，导致湿地在上游接受地下水补给、在下游排泄至地下水，地下水流“贯穿”整个湿地。

图5 模拟区多年平均稳定流流场示意

图6 计算地下水位与实测地下水位对比

图7 渗透系数分区示意

基于龙湖周边地下水水位统测结果绘制了地下水位剖面图(图8)。由图8可以看出，地下水由西北向东南水位逐渐降低，即由西北流向东南；
龙湖西北侧接受地下水补给，东南侧向地下水排泄，属于典型的“饱和流-贯穿型”湿地。

注：D2,D5,D7,D20,D10,D8表示钻孔编号。

2.2 龙湖湿地地表水-地下水交换量

根据模型地下水均衡计算结果，地下水对龙湖多年平均补给量为617×104m3/a，排泄量为253×104m3/a，净补给量为364×104m3/a。

对龙湖水文过程进行分析。龙湖主要来水包括降水、降水径流、地下水补给、污水排放等，主要出水包括蒸发、向地下水排泄等，年际水量平衡关系如下式所示：

ΔW=(P-E)+(Roverland-roverland)+

(Runderground-runderground)+(Qin-Qout)

式中：ΔW为龙湖蓄水水量，m3；
P为湖面降水量，m3；
E为湖面总蒸发散量，m3；
Roverland为湖区地表流入量，m3；
roverland为地表流出量，m3；
Runderground为地下径流补给量，m3；
runderground为地下径流排泄量，m3；
Qin为污水溢流入湖量，m3；
Qout为地表水取水量，m3。

水量平衡计算如表1所示。由表1可知，龙湖总水补给量为1 466×104m3，其中地下水补给量6.7×104m3，占龙湖来水总量的42%，是最大补给来源，而降水量是第二大补给来源，占龙湖补给量39%。

表1 龙湖多年平均水量平衡计算

2.3 龙湖水系现有交互模式下龙湖水量风险分析

从前文可知，龙湖水体水量平衡和水位波动很大程度受周边地下水动态波动和降水量的直接影响。由该地区历年(1971～2017年)降水量的分析可知，当地降水量存在年际变化，但总体相对稳定，没有明显减少或增加趋势(图9)。由此可以推断，地下水水位变化是影响龙湖水位波动的主要原因。

图9 降水量年际变化趋势分析

一定区域的地下水位通常主要随着降水量或地下水开采量变化。在降水量相对稳定的情况下，浅层地下水开采则是地下水位主要影响因素。资料显示，2009～2016年龙湖周边浅层地下水水位总体呈下降趋势(图10)，而从淮阳地区浅层地下水供水量年际变化看，研究区历年浅层地下水开采量均维持在较高水平(图11)，相关区域存在地下水超采问题[15-16]，可见，地下水超采是中原地区水资源面临的一大威胁，也间接为龙湖带来补给水量短缺的风险。

2.4 水文地球化学交互概念模式及风险分析

对7个取样点水样进行检测分析，主要检测结果分布见图4。湖水TDS为937～1 143 mg/L，周边地下水TDS均小于1 000 mg/L。

图12为龙湖及周边地下水水样点Piper三线图。由图4、图12可以看出，龙湖西北侧补给区，沿径流方向TDS呈上升趋势，地下水化学类型逐渐由HCO3-Na·Ca·Mg型转变HCO3-Na型；
龙湖东南侧排泄区，沿径流方向TDS降低，总体呈“低-高-低”的变化，表明湿地地表水与地下水之间存在明显的水文地球化学联系，并形成独特的水文化学循环模式。在龙湖水体中，TDS受蒸发浓缩作用明显。

图11 淮阳区浅层地下水供水量年际变化

注：2009～2016为年份，1～12为月份。

2.5 修复策略

(1)龙湖周边地下水与龙湖水量交互作用密切，地下水位较大幅度波动会影响湖泊的水量平衡。2016年之前的数据表明龙湖区域地下水有持续下降的趋势，这将对龙湖湿地的生态平衡造成威胁，迫切需要加强对该地区及周边市县更大范围地下水开采的管控。

注：2009～2016为年份，1～12为月份。23～25号，27号为地下水长期监测井编号。

(2)对龙湖及周边水系环境进行全面治理，在达到一定条件后重新恢复龙湖与周边地表水系的联系，一方面扩大龙湖地表水补给范围，另一方面通过水系连通改善龙湖水质，消除蒸发浓缩作用带来的盐碱化风险。在极端条件下，可通过湖外水系通道调沙颍河水补给龙湖。

(3)在全球气候变化大背景下，区域降水量变化会加大龙湖湿地水量平衡的不确定性，有必要提前制定相应对策。建立龙湖水系地表水、地下水监控体系，通过对水位、水量、水质、矿化度等参数及其相互响应机制和对生态系统影响的分析，制定具有针对性的管理策略并能适时调整。

本文以河南省周口市龙湖湿地作为研究区，利用收集的水位监测、统测资料、水质检测资料，应用Visual MODFLOW构建了地下水流数值模型。对龙湖地表水-地下水交互作用模式、交换水量及水文地球化学交互模式进行了分析研究。得出如下结论：

(1) 龙湖周边地下水由西北流向东南属于典型的饱和流-贯穿型湿地。

(2) 龙湖湿地与浅层地下水关系十分密切，地下水补给量占龙湖来水总量的42%，地下水是龙湖湿地最主要的补给来源。地下水位较大幅度波动会影响湖泊的水量平衡，应加强对该地区及周边市县地下水开采的管控。

(3) 沿龙湖“西北侧-龙湖-东南侧”方向，TDS含量呈现“低-高-低”的变化趋势，龙湖湿地地表水与地下水之间存在明显的水化学联系。同时由于蒸发浓缩作用，湖泊也存在盐碱化风险。

(4) 通过模型计算量化地下水-湖泊水量交换，明确湖泊湿地水量平衡机理，为制定风险管控对策奠定基础，这种方法可以为同类型湖泊湿地水平衡研究和修复决策提供参考。