傍河水源地地表水地下水转化SFR模型的应用

滕飞达，张 博，梁秀娟，李明乾，杨伟飞，张芷豪

（1．吉林大学地下水资源与环境教育部重点实验室，吉林 长春130021；2．油页岩地下原位转化与钻采技术国家地方联合工程实验室，吉林长春130021；3．吉林大学新能源与环境学院，吉林长春130021）

傍河水源地含水层通常具有富水性强、地下水位埋藏浅、水量稳定、水质优良、易集中开采和便于管理等优点，是地下水开发常见类型［1-2］。河流和地下水之间存在着密切的水力联系和频繁的转化［3］，地下水资源实际上有相当大的一部分是袭夺了河水资源，地下水可开采量的大小和河流的最枯流量、含水层厚度和导水性有关［4］。傍河取水能够激发并增大河水对含水层的入渗补给［5］，以此实现地表水与地下水资源的联合利用。国内外很多学者采用地表水与地下水耦合模型，研究地下水与地表水的交互关系。Jolly等［6］指出湿地地下水和地表水交互作用受地表水水位和地下水水位之差、局部地形和地下水流场等因素的影响。丁元芳等［7］应用MODFLOW对浑河李官堡傍河地下水源地流场进行数值模拟，将浑河概化为一类定水头边界，试算不同的开采方案模拟出相应的流场。王丽雪［8］对第二松花江流域傍河水源地进行地下水数值模拟，将第二松花江及其支流饮马河与伊通河均概化为河流边界。张奇莹等［9］将神木县马镇合河村境内的黄河概化为河流边界、上下游边界概化为一类定水头边界，模拟地下水资源傍河取水开采方案。国内研究通常将河流概化为一类定水头边界或河流边界处理，通过模拟出的流场设计开采方案，没有充分考虑枯、丰水期河流水位与流量的差异及水均衡中地表水与地下水转化量对开采方案的影响。为弥补以往研究的不足，笔者应用SFR模型概化白马河-吉利河，建立地表水与地下水耦合模型，结合地表水与地下水转化量及流场设计开采量方案。

1 研究区概况

白马河-吉利河流域位于青岛市，属暖温带沿海湿润季风区［10］，多年平均气温为12．5℃，年平均相对湿度为74%，年平均降水量为 830 mm［11］。7—9月为丰水期，4—6月、10月为平水期，11月—次年3月为枯水期，丰水期降水量占全年降水量的60%以上，河水丰、枯水期流量差别巨大，在枯水年枯水期会出现断流现象。吉利河干流长39．8 km，流域面积285．1 km2；白马河干流长 44．2 km，流域面积 262．1 km2。 研究区位于白马河-吉利河流域中下游，总面积为148．5 km2。河流北、东、西面为山地，两河河谷成Y形。

白马河第四系松散岩类含水层组按成因可划分为冲洪积层、砂砾石层，坡积-坡洪积黏质砂土夹碎石层、冲积-海积砂层。地下水富水性分区见图1，地下水补给途径为大气降水、灌溉入渗和河流入渗补给。基岩区接受大气降水入渗补给后，地下水由两侧向河谷流动，进入平原。地下水运动方向大体与地表径流方向一致，两河成为主要排泄带。地下水排泄方式主要为人工开采、潜水蒸发和径流排泄。地下水动态季节性变化明显，6月初地下水位下降至最低，7—9月丰水期水位回升速度较快。丰水年水位回升，枯水年水位下降。

图1 地下水富水性分区

2 研究方法

SFR（Streamflow-Routing）模型是美国地质勘探局（USGS）开发的河床汇流模型［12］，原理为省略加速项的扩散波方程。河流边界模型也是由USGS开发的，通过设置上下游河流水位、河床高程及水力传导系数等参数概化河流，与SFR模型的共同点都是通过河床水力传导实现地下水与河水之间的转换，但是它不具备汇流功能，其水位均由观测获得，需要大量长时间的观测数据，不如SFR模型灵活。

白马河-吉利河中下游地区取水来源以地下水为主、地表水为辅，本文综合考虑河水与含水层饱和带水的相互联系，根据河床与地下水之间的河床水力传导系数、河床厚度、水深及断面形状等参数，将SFR模型与MODFLOW模型进行耦合，计算丰、平、枯水期河流流量与上下游水位，使河流概化更符合实际情况、河水与地下水转化量计算更加精确，提高水源地地下水资源量的计算精度；再分别将白马河-吉利河概化为第一类定水头边界模型、河流边界模型，计算3种河流概化情形下的地下水可开采量，并与实际勘查的地下水最大可开采量进行对比分析。

2．1 地下水数学模型

根据水文地质概念模型，建立的地下水数学模型为

式中：h为地下水头，m；kx、ky、kz分别为x、y、z方向的渗透系数，m／d；h1为含水层第一类边界水头，m；h0为含水层初始水头，m；W为开采强度，m3／d；Σ1为含水层第一类边界；D为研究区；μ为给水系数；t为时间，d。

2．2 SFR地表水与地下水耦合模型

根据白马河-吉利河水文资料构建SFR耦合模型。河水一方面接受地表汇流的补给，一方面与地下水之间发生交换。河床与地下含水层之间的汇流表示为q地表汇流和q地下汇流两部分。 为简化计算，省略河水的加速运动，选用曼宁公式描述河水运动，于是圣维南方程组可改写为

式中：B为水面宽，m；q地表汇流、q地下汇流分别为地表汇流与地下汇流的流量，m3／d；S0为河床纵比降；Sf为摩阻比降；Q为上游来水流量，m3／d。

方程的源汇项q地表汇流和q地下汇流分别由地表降水径流和地下汇流计算得出。地下水数学模型的边界条件是河流扩散波方程提供的水位和流量，而河流扩散波方程的源汇项又是由地下水流微分方程计算得到的。

2．3 河水与地下水交换模式

河水与地下水交换分为3种模式：一种是地下水位高于河床高程，另一种是地下水位略低于河床高程，这两种均属于饱和渗流机制，河床垂向和侧向渗流量是地下水与河水水头差的函数；第三种是地下水位远低于河床高程，二者之间形成非饱和渗流。无论是哪种模式，河床与含水层的水位差决定着流入和流出量。河道与地下水的非稳定流耦合模型是基于对每一时间步长的河床矩阵与地下水矩阵进行耦合求解，从而提高了交互计算的性能。把地下水流数值方程与扩散波方程联立求解，即得到q地下汇流。河床水力传导系数计算公式为

式中：C为河床水力传导系数；k为河床渗透系数，m／d；w为河床宽度，m；T为河床厚度，m；L为上下游距离，m。

2．4 所用资料

模型所用数据资料包括：①2017年（现状年）流量监测及降水资料（见图2）；②现状年实测地下水位资料，地下水位监测点（25处）分布于研究区全境，流量监测点位于中游（见图3）；③上下游河床形状概化资料（见图 4）。

图2 现状年月降水量、流量过程线

图3 地下水位及流量监测点分布

图4 河道形状概化

3 模型构建与处理

3．1 水文地质概念模型构建

根据研究区的水文地质条件与水动力场，确定概念模型。①含水层概化：把含水层概化为3层，第一层为透水性稍弱的细砂层，第二层为透水性很强的粗砂砾石层，第三层为基岩（见图5）；②水动力条件概化：将研究区概化为非均质各向同性的含水系统，地下水运移符合达西定律，水位随时间变化，地下水流可概化为三维达西流；③边界条件概化：白马河-吉利河下游连通海洋，设置为定水头边界，入海口水位常年稳定，高程为0，概化为一类水头边界，河谷东部、西部和北部均为山地，将地下水分水岭概化为零流量边界；④源汇项概化：将河流边界分别概化为SFR耦合模型、河流边界模型和第一类定水头边界模型3种边界。大气降水和灌溉回渗作面状入渗处理；通过ET模块设定蒸发极限深度和蒸发速率，计算蒸发量；水源地开采排泄作点状和面状排泄处理。

图5 研究区模型概化

3．2 河流动力学概念模型构建

根据现场踏勘，将上、中、下游河床概化为变宽度的梯形河床，河床形状和坡度按照DEM确定。将白马河-吉利河概化为SFR耦合模型，设置现状年流量、上下游水位、河床渗透系数、河床厚度、水深及断面形状等，使河流概化更符合实际情况，再将河流概化为第一类定水头边界和河流边界，对比分析3种不同河流概化条件下模拟出的地表水与地下水转化量、允许开采的地下水量的差异，以判定各河流概化方法的优劣。

3．3 模型识别与验证

利用SFR耦合模型模拟研究区枯水期地下水位，并与实际枯水期水位进行对比。选用现状年2017年1—3月、7—9月两个时段分别进行模型的识别与验证（见图6、图7）。模型识别期计算水位与实测水位拟合误差小于0．5 m的数据点占已知水位数据点数量的75%以上，模型验证期丰水期流场与实测流场非常接近，表明模型具有较高的精度和可信度。

图6 模型识别期实测水位与计算流场

图7 模型验证期实测水位与计算流场

3．4 傍河井点开采量设置

沿河两岸布置130口开采井，当模型初始地下水位下降含水层2／3厚度时即为该研究区地下水可开采量，对3种模型设置相同位置的开采井，计算SFR耦合模型、河流边界模型及第一类定水头边界模型条件下的地下水开采量，对以上3种情况模拟地下水位下降含水层2／3厚度时的地下水最大允许开采量。通过计算，SFR耦合模型、河流边界模型、第一类定水头边界模型模拟的丰、平、枯水期每口开采井开采量分别为700～900、700～800、900～1000 m3／d。

4 讨 论

采用传统的地下水数值模型进行傍河地下水开采量评价，可能会出现一定的偏差。与传统的地下水数值模型相比，SFR耦合模型考虑了河流丰、枯水期流量变化对地下水开采的影响。为了验证SFR耦合模型的优越性，采用河流边界模型和第一类定水头边界模型建立研究区傍河开采地下水数值模型，与SFR耦合模型的水均衡计算结果进行对比（见表1）。

表1 水均衡计算结果万 m3／d

4．1 SFR耦合模型与定水头边界模型对比分析

4．1．1 地表水与地下水转化量

从表1可知：丰水期SFR耦合模型、第一类定水头边界模型计算的河水补给地下水量分别为3．08万、7．55万 m3／d，地下水转化为河水的量分别为 2．66 万、5．67万m3／d；平水期SFR耦合模型、第一类定水头边界模型计算的河水补给地下水量分别为 2．39万、9．46万 m3／d，地下水转化为河水的量分别为 2．30 万、5．78万m3／d；枯水期SFR耦合模型、第一类定水头边界模型计算的河水补给地下水量分别为2．27万、9．99万m3／d，地下水转化为河水的量分别为 2．40 万、4．72 万m3／d。

由表1可知，第一类定水头边界模型计算的河流补给量是最大的，且枯水期河水对地下水的补给量最大，丰、平、枯水期都是河流补给量大于地下水转化为河水的量，即地表水补给地下水。SFR耦合模型计算的地表水与地下水转化较为均衡，丰水期河流补给量略大于地下水转化为河水的量，为地表水补给地下水；平水期河流补给量与地下水转化为河水的量基本持平；枯水期河流补给量略小于地下水转化为河水的量，为地下水补给地表水。定水头模型计算的枯水期河水补给地下水的量大于丰水期河水补给地下水的量，是“无限供水的边界”，不符合白马河-吉利河季节性河流特征。综上，SFR耦合模型相对于定水头边界模型更能准确计算出地表水与地下水的转化量，计算结果更为合理。

4．1．2 模拟开采量

由表1可知，丰水期SFR耦合模型、第一类定水头边界模型计算的可开采量分别为8．39万、10．82万m3／d，平水期SFR耦合模型、第一类定水头边界模型计算的可开采量分别为 7．4 万、9．24 万 m3／d，枯水期SFR耦合模型、第一类定水头边界模型计算的可开采量分别为 7．05 万、8．8 万 m3／d。 由表 1 可知，第一类定水头边界模型模拟的傍河开采量比SFR耦合模型的大，第一类定水头边界模型计算的最大允许开采量为3 511万m3／a，SFR耦合模型计算的最大开采量为2 779万m3／a，参照《青岛市胶南县吉利河-白马河地区供水水文地质勘察报告》及《白马河-吉利河地下水模型及脆弱性安全性评价》等以往研究资料，研究区平水年地下水可开采量为2 870万m3／a，SFR耦合模型比第一类定水头边界模型计算结果更符合实际。

4．2 SFR耦合模型与河流边界模型对比分析

4．2．1 地表水与地下水转化量

由表1可知，丰水期SFR耦合模型、河流边界模型计算的河水补给地下水量分别为3．08万、2．83万m3／d；平水期SFR耦合模型、河流边界模型计算的河流补给地下水量分别为 2．39 万、1．77 万 m3／d，地下水转化为河水的量分别为 2．3 万、3．57 万 m3／d；枯水期SFR耦合模型、河流边界模型河计算的水补给地下水量分别为 2．27 万、1．42 万 m3／d，河水补给地下水量分别为 2．40 万、2．85万 m3／d。 河流边界模型计算的丰、平、枯水期地下水转化为河水的量与SFR耦合模型的相当，但地下水转化为河水的量较大，即地下水补给河水。与河流边界模型相比，SFR耦合模型具备汇流功能，地表水与地下水转化量更为均衡，与研究区丰、平、枯水期地表水与地下水实际转化关系更为相符。

4．2．2 模拟开采量

由表1可知，丰水期SFR耦合模型、河流边界模型计算的可开采量分别为 8．39 万、7．40 万 m3／d，平水期SFR耦合模型、河流边界模型计算的可开采量分别为 7．40 万、7．00 万 m3／d，枯水期 SFR 耦合模型、河流边界模型计算的可开采量分别为 7．05万、6．77万m3／d。综上，SFR耦合模型计算的开采量在丰、平、枯水期均略大于河流边界模型的，SFR耦合模型计算的开采量为2 779万m3／a，河流边界模型计算的开采量为2 576万m3／a（远小于研究区平水年地下水可开采量2 870万m3／a）。与河流边界模型相比，SFR耦合模型计算结果更符合实际。

5 结 论

本文综合考虑白马河-吉利河在丰、平、枯水期的流量差异，应用SFR地表水与地下水耦合模型，使河流概化更符合实际情况，地表水与地下水的转化量计算更加精确。SFR地表水与地下水耦合模型比河流边界模型更能精细概化河道特征，计算的开采量更接近实际情况。在枯水期部分河段断流时，无河水补给地下水，而第一类定水头边界模型为无限补给，使计算的地下水资源量明显偏大，与实际不符。综上，SFR地表水地下水耦合模型比传统的地下水数值模型更适合于傍河地下水水源地地表水与地下水转化量计算及开采量评价。