大清河淀西平原地下水数值模型研究

李 巍，郭抒燕，薛丽娟，丁晓梦

（1.水利部海河水利委员会，天津 300170；2.天津市龙网科技发展有限公司，天津 300170）

我国华北地下水水资源短缺、地下水超采问题由来已久[1-4]，出现一系列环境、生态问题[5-7]。数值模型是研究地下水问题的有效手段之一，韩忠等人通过对MODFLOW 的二次开发，解决了面状补给的问题，丰富了源汇项包，提高了模型运行效率，形成优化后的MODFLOW-RAW 主程序、RAW 面状源汇项包和改进后的RCH、EVT 程序包[8，9]。基于大清河淀西平原地下水农业开采数据以面状为主的特点，为满足地区发展和国家工农业发展规划的规定及地下水资源的可持续发展利用，有必要基于MODFLOW-RAW 对大清河淀西平原地下水资源做出新的评价。而建立大清河淀西平原地下水数值模型是研究、管理地下水的重要技术手段和有效方法。

1 研究区概况

1.1 水文地质

大清河淀西平原区主要由山前冲洪积倾斜平原、冲洪积扇和中部冲积平原构成，山前冲洪积平原、冲洪积扇呈扇状交错分布于山前，含水砂层主要由砂砾石、粗砂、中砂、中细砂等各类砂、砂砾石组成，含水层下部无连续隔水层，垂向水力联系良好，常因下部含水层的开采而疏干，因此在山前常将这一含水层组与下部承压含水层组看作是统一含水层组。

平原区地下水系统的上游为大清河冲洪积扇子系统，分为4个亚级，包括拒马河冲洪积扇系统、瀑河-漕河冲洪积扇系统、唐河-界河冲洪积扇系统、大沙河-磁河冲洪积扇系统；下游为大清河古河道带子系统。

1.2 含水组划分

由区域水文地质剖面图以及查阅历史资料可知，大清河淀西平原分4个含水岩组。

（1）第一含水层组。底界面埋深一般小于50 m，在山前平原冲洪积扇和扇间、扇前地带具有不同的水文地质特征。在冲洪积扇地区，含水层粒度大、厚度30～50 m、垂向连续性强，属单层或双层结构，透水性强，导水系数多大于5 000 m2/d，单位涌水量为20～30 m3（/h·m）。含水层直接裸露，或者被薄层砂质黏土覆盖，具有强入渗补给和储存条件，且山前地区与山区河谷含水体相连，具有侧向径流补给条件。

（2）第二含水层组。底界埋深一般120～210 m，在山前平原，与第一含水层组之间缺乏稳定的隔水层，二者之间具有较好的水力联系。自西向东发育2～3 套中细砂～中粗砂～砂砾石韵律层，含水层透水性与导水性均比第一含水层组强。目前，由于2个含水层组混合开采，人工加强了二者之间的水力联系。

（3）第三含水层组。底界埋深250～310 m，几乎全部为淡水。山前平原含水层呈扇状、扇裙状展布，由3～4套中细砂～中粗砂～砾石韵律层组成，下段含水层更受到不同程度的风化。典型平原区南部扇体内单井涌水量20～50 m3（/h·m），扇间地带单井涌水量10～20 m3（/h·m）。在大型扇体内部与上覆第二含水层组之间无连续分布的隔水层，二者水力联系良好，在其他地段一般都有5～10 m 的黏土或砂质黏土分布，水力联系变弱。

（4）第四含水层组。底界为第四系基底，由于少有勘探井揭穿本含水岩组底板，底板埋深及地下水厚度不清。该层地下水水力性质均为承压水，TDS较第三含水层组略有增高。

1.3 地下水开发利用

由于多年的开采，形成了保定市区漏斗、一亩泉漏斗和高蠡清漏斗。保定市区漏斗形成于1974 年低水期，但季节性非常明显，进入20 世纪80 年代以后，就转为常年性漏斗，漏斗中心地下水位埋深及影响面积从1980 年的15.14 m 和136.3 km2发展到2000 年的32.0 m 和241.0 km2，20 a 间最大埋深增加了16.83 m，影响面积增大了104.7 km2。保定市区漏斗的发展主要原因是市区自备井的无序取水。

平原区内包括四大灌区，分别为房涞涿灌区、易水灌区、唐河灌区、沙河灌区。灌区多引水库或地表河流水系水进行农业灌溉，但是近年来，由于降雨减小地表径流被上游水库拦截等因素影响，灌区采用地下水的比例有大幅度增大，这对地下水模型的建立与验证产生一定影响。

2 基于优化MODFLOW的淀西平原地下水数值模型

以下通过建立大清河淀西平原地下水数值模型，并结合南水北调受水区调水方案，预测在不同地下水开发利用条件下的地下水位、水量变化，为指导地下水开发利用、地下水管理提供科学支撑。

2.1 MODFLOW 的优化

针对大区域模型中源汇项数据处理的难题，韩忠等人在MODFLOW2005 源程序的基础上进行了优化，开发了MODFLOW-RAW，并利用模型降雨、蒸发和开采数据处理的新方法，简化模型的数据处理过程，缩短模型运行时间[8，9]。其中，包括降雨子程序包（RCH）、蒸发子程序包（EVT）的处理方法改进和面状源汇项程序包（RAW）的开发，MODFLOW-RAW 系列程序包不仅通过改变数据存储方式减少了人工处理数据的工作量、加快了模型数据的读取速度，且压缩后的源汇项文件更有利于数据的修改及模型的调试。

同时，借助RAW 程序包具可同时处理最多4 个面状源汇项的平行处理能力，可将区域农业开采、灌区开采、灌区回灌和白洋淀入渗分别计入，使模型的前、后处理更加清晰、便捷。

2.2 水文地质概念模型

针对淀西平原区第四系松散含水层组建模。该区域地质条件复杂，目前的水文地质工作基础尚不足以详细刻画控制全区的三维物理渗流结构，能够取得的较为详细的地质调查成果为前述的4层结构的含水层系统，对区域地下水系统结构已经做了相当的概化。本次地下水数值模型对区域地下水含水层组进行进一步的概化，将第一含水层组和第二含水层组合并作为浅层地下水系统，主要模拟潜水循环运动规律，同时也包含与潜水含水层紧密相关的微承压含水层的循环在内；将第三层和第四层合并作为统一的深层地下水系统，主要模拟区域第三层承压水的循环运动规律，同时也包含循环量较小的第四层承压水的循环在内。

含水层厚度大、分布广、地下水流以水平方向为主，垂向运动较弱但不可忽视，地下水运动空间概化为三维流；源汇项和水位均随时间发生变化，表现为明显的非稳定流特征；水文地质参数随空间变化，体现非均质性，在同一点的水平方向上，参数无明显方向性，可视为各向同性。

综上，由于第一含水层组和第二含水层组水力联系紧密，且实际开采多以混采井居多，故将第一含水层组和第二含水层组合并概化为第一含水层，将第三含水层组作为第二含水层，第四含水层组作为第三含水层。其中，第一含水层为潜水含水层，第二、三含水层为承压含水层。综合考虑地下水开发利用层位的特点，参照含水层的发育程度、含水层的渗透性、地下水水力性质、地下水动态特征等因素，将模型在垂向上概化为三层的非均质各项同性、空间三维结构的非稳定地下水流动系统。

2.3 地下水数学模型

根据水文地质概念模型，对非均质各向同性、空间三维结构的非稳定地下水流动系统，可用如下方程的定解问题来描述：

式中：Ω为渗流区域；h为含水层的水位标高（m）；Kx、Ky、Kz分别为水平和垂向渗透系数（m/d）；ε为含水层的源汇项（1/d）；h0为含水层的初始水位分布（m）；S为含水介质的储水率（1/m）；Kn为边界面法向方向的渗透系数（m/d）；τ1为渗流区域的侧向和下边界；q为τ1边界的流量（m/d），流入为正，流出为负，隔水边界为0；n为边界面的法线方向；Q为流量（m3/s）；K为渗透系数（m/d）；H为外扩边界水头（m）；H0为模型边界水头（m）；a为单个单元格（网格）长度（m）；b为单个单元格（网格）宽度（m）；L为渗流路径长度（m）；t为时间（d）。

2.4 地下水数值模型

2.4.1 模型离散

（1）时间离散。模拟期为2000—2010 年，为反映年内的周期变化，以月单位作为一个应力期，每个应力期为一个时间步长。

（2）空间离散。根据区域的含水层结构、边界条件和地下水流场特征，将模拟区第一层剖分为220行、160 列的1 km×1 km 规则网格，第一层有效单元共12 554个，下部第二、三层7 983个。

2.4.2 定解条件

（1）初始水位。以2000 年1 月的地下水水位资料为基础，参考长序列数据，绘制初始流场，采用内插法和外推法获取各含水层的初始水位。

（2）边界条件。模型根据水文地质条件，可将区域划分为2 类边界，即流量边界和通用水头边界。区域西部山前侧向流入边界，月平均径流量为2.11×106m3；东部的边界以行政区划分，定为通用水头边界。模拟边界，如图1所示。

图1 淀西平原模拟边界示意

2.4.3 源汇项处理

以区县为单位，蒸发量采用改进的MODFLOWRAW 中的EVT 蒸发程序包处理，根据蒸发极限埋深、蒸发速率及模拟水位确定。降雨入渗采用改进的MODFLOW-RAW 中的RCH 程序包处理。灌区开采和回灌以及其他开采包括生活工业用水均采用面状补排处理方式。开采山前侧向流量采用GMS中的WELL井流程序包处理。具体应用中，在每个网格中心点上设置一个注水井，然后将线状的流量数据都分配到网格中心点上，将侧流量以点井的形式导入模型中。这种源汇处理方式效率高，便于调参过程中源汇数据的修改。补给排泄项对应的程序包，详见表1。

表1 模型中源汇项的数据对照

2.4.4 模型识别与校验

通过模拟水位与实测水位的对比来检验长序列模型的可靠性。对研究区内19 个观测点共计2 267个数据进行了误差分析，其平均绝对误差为2 m。其中，11 个观测孔的平均绝对误差小于1.5 m，4 个观测孔的平均绝对误差介于1.5~2 m，4个观测孔的平均绝对误差大于2 m。综上，模拟区内误差2 m 之内的观测孔约占80%，整体拟合效果较好。拟合误差分布，详见表2。

表2 拟合误差分布

2.5 模型应用

基于拟合验证后的大清河淀西平原地下水数值模型，分别建立现状和南水北调实施后2 种开采条件下的10 a预测模型，通过方案对比，定量分析地下水调控对地下水资源的影响。南水北调压采条件与现状开采条件下各市地下水储变量统计，详见表3。

表3 南水北调压采条件与现状开采条件下各市地下水储变量统计 108 m3/a

从储变量变化上对2 个方案进行比较，可以看出，南水北调压采实施后，浅层地下水储变量亏损减少4.19×108m3/a，约下降32%，浅层地下水亏损的状况得到了较好缓解。浅层地下水-30 m 水位等值线的范围有了明显缩小，面积约从0.62×104km2缩减到0.41×104km2，下降34%，降落漏斗程度减轻。深层地下水的储变量也有所减少，但由于区域内对深层地下水的开采相比浅层要小很多，因而减少幅度不大。总体上而言，在南水北调压采实施后由于减少了人工开采量，有效地控制了现有降落漏斗的扩大和加深，对受水区的地下水恢复起到了很大作用，同时对缓解受水区的地下水亏损状况和平原区的地下水紧张局面、涵养和保护地下水资源起到了很好的积极作用。

3 结论

（1）对MODFLOW 的源程序进行了优化应用。使用改进后的MODFLOW 程序模块，实现了基于面状处理方式的各源汇项的分类处理，在保证水均衡与模拟精度要求合理的前提下，可以较好地提高地下水流数值模型的运算速率。通过模型的外部建模方式更好地实现了各个运行文件参与平台的集成以及后期数据的分项调整功能，弥补了MODFLOW 源汇项处理过于简单的问题。

（2）构建了水文地质概念模型。通过水文地质条件分析确定了淀西平原为非均质各向同性、空间三维结构、非稳定地下水流系统的水文地质概念模型，确定了模型的范围结构、边界条件、水文地质参数及各源汇项。

（3）建立了长时间序列地下水流数值模型。在概念模型的基础上建立了大清河淀西平原的地下水流数值模型，运用模拟期地下水长期观测孔水位、各年各层地下水水位等值线对模型进行识别验证，通过这些实测数据对区域内地下水开采量进行估算和评价。总体看来，浅层含水层拟合效果良好，较好地反映了实际水文地质条件。

（4）建立了不同情境下的预测模型，预测分析了地下水流场变化趋势。从总体的预测结果来看，在南水北调压采方案实施后地下水环境得到了较好的改善，有效地控制了现有降落漏斗的扩大和加深，对缓解受水区地下水的亏损状况和平原区的地下水紧张局面、涵养和保护地下水资源起到了很好的积极作用，对受水区的地下水恢复起到了很大作用。