基于地表水地下水耦合模型的白马河傍河水源地开采量评价

董 杰，付佳妮，滕飞达

（1.中国海洋大学 环境科学与工程学院，山东 青岛266100；2.青岛地质工程勘察院，山东 青岛266100；3.山东省地矿局 城市地质与地下空间资源重点实验室，山东 青岛266100；4.吉林大学 新能源与环境学院，吉林 长春130021）

地表水与地下水之间的相互作用对于河流生态系统和河岸管理非常重要［1］，地表水与地下水转换过程及其耦合模拟是水资源开发利用和科学评价的基础［2］。傍河开采地下水资源是一个比较复杂的地下水与地表水交互作用过程［3］，是人类开发利用地下水资源的主要方式之一［4］。傍河水源地开采的地下水在很大程度上袭夺了河水资源，而不全是地下水资源［5］。在区域范围内综合考虑地表水和地下水储存，进行统一规划和管理的系统性尝试仍然很少［6］。Gilbert等［7］考虑到山区水力传导率的不确定性。Brunner等［8］认为河流与地下含水层之间有水力联系时才能进行统一评价和管理。Ledoux等［9］应用MC模型整合了法国Caramy流域表面流动、非饱和区域流动、地下水流动以及河流和地下水位之间的相互作用。

目前，在国内基于地下水数值模型的傍河水源地开采量评价中，常常把河水作为定水头边界来处理。该种处理本质上是把河流当成无限水源，未考虑到河流流量、水位的季节变化与地表、地下水之间的转化关系，其结论具有一定偏差。实际上，模型中只要有一小段边界是定水头边界，该边界理论上就能提供无限的水量。易树平等［10］将浑河概化为一类定水头边界，刘洋［11］将模拟区内第二松花江概化为一类定水头边界，孙跃等［12］将安徽夹江概化为第三类河流边界。以往的评价工作通常将河流概化为定水头边界，如果河流较大且季节性变化不大，开采量与河水流量相比非常微弱，这种已知水头假设还基本可行，一旦开采量过大袭夺地表水资源，河流水位受傍河开采地下水影响发生变化，定水头边界的假设就是错误的。

白马河为季节性河流，河水位季节性变化较大。该区傍河开采大量地下水，对河水位的剧烈变化起到了推波助澜的作用。考虑到河流边界为定水头边界的传统地下水数值模型已经不适合这种傍河型水源地开采评价。本研究运用流域水循环理论，提出一种河流边界的处理方法，即使用河流扩散波方程作为地下水数值模型中的河流边界，取代传统地下水模型中的定水头边界或流量边界，用河流扩散波方程模拟河水流量和水位变化，将扩散波方程与地下水数值方程耦合，以更合理地处理傍河水源地地下水开采问题。

1 研究区概况

白马河流域位于山东半岛胶南地区，属青岛胶南市管辖，为暖温带湿润季风气候区［13］。多年平均降水量为830 mm，多年平均蒸发量1 541.5 mm［14］。白马河流域上宽下窄，干流长44.2 km。研究区位于白马河中、下游地段，模拟区总面积46.60 km2。白马河富水区沿古河道主流带分布，该区主要傍河开采地下水，开采量多来自于河水补给。该区12月—次年3月为枯水期，4—7月为平水期，8—11月为丰水期。河水年内流量及水位波动大。区域水资源开发利用以地下水为主，地下水开采量大时可引起部分河段断流，传统的地下水模型不能模拟河流丰枯平水期变流量变水位问题。

白马河流域地层结构为二元结构，上层为细砂和粉细砂的弱透水层或相对隔水层，下部为粗砂及砂砾石，含水层下部为基岩。水文地质含水层富水性及参数分区见图1。流域水资源开发利用以地下水为主，河流是地下水的重要边界条件。

2 研究方法

首先考虑白马河流量在丰枯水期差异，应用河流扩散波方程，模拟丰、枯水期流量和地下水流场。再将河流概化为定水头边界，试算该概化边界条件下地下水最大开采量，对比分析两种概化情况下地表水与地下水转化量的区别，进而分析本文方法的合理性。

2.1 水文地质模型

①概化含水层：根据地层结构情况，模型概化为两层，上层为透水性相对较差的细砂粉砂层，下层为粗砂砾石层；②边界条件：忽略潮汐引起的水位波动，流域下游沿海岸概化为已知水头边界，水位0 m，使用SEAWAT［15］模拟海水入侵。模拟区东、北、西侧为流域分水岭，概化为隔水边界；③源汇项：补给源分为大气降水、灌溉补给及河流补给，其中河流边界概化为由扩散波方程控制流量和水位的第三类边界，排泄项为蒸发及开采。

2.2 地下水模型

根据水文地质概念模型建立相应的数学模型如下：

图1 地下水富水性及参数分区图Fig.1 Groundwater rich and water partition

其中：h—地下水水头，m；Kx、Ky、Kz—x、y、z方向渗透系数，m/d；h1—含水层第一类边界水头，m；h0—含水层初始水头，m；W—源汇项强度（包括开采强度等），m3/d；Σ1—含水层第一类边界；D-研究区；μ-给水系数。

2.3 河流概念模型

按照数字高程模型（digital elevation model，DEM）设定河床形状和比降。将河水运动概化为省略加速运动的一维扩散波。河水来源于地表坡面漫流的汇入，按照本区的径流系数将大气降水平均分配到每个河段之上，作为地表坡面漫流汇入量。白马河的主要排泄方式有两种：一种是排泄入海；一类被地下水开采袭夺，进入含水层。本次模拟主要验证河流被地下水开采袭夺的流量。

按经验值设定上游、中游及下游的河床水力传导系数（C）、河床渗透系数（k）、河床宽度（w）及河床厚度（t）等参数见表1。

表1 河道概化参数Tab.1 River generalization parameters

2.4 地表水模型

河水通过河流边界与地下含水层之间发生水量交换。采用河流扩散波方程描述如下：

其中：Q—河流流量，m3/d；H—水位，m；t—时间，d；x—河道沿程距离，m；c—水流流速，m/d；D-水力扩散系数，m2/d。扩散波方程的源汇项边界条件为河流起始单元的流量过程，地下水流微分方程的边界条件是河流扩散波方程提供的水位和流量，河流扩散方程的源汇项通过地下水流微分方程计算得到。

2.5 二者耦合机制

河道水位和地下水位分别表达在扩散波方程和地下水流方程之中，示意图见图2。当河水位大于地下水位时，渗流由河流进入地下水，反之则由含水层进入河流。如果二者相等，则河道与地下水之间不发生水量交换。耦合计算时，用河水位代入地下水流数值方程，根据河水位与地下水位的关系，计算河流与地下水的补排水量，得到时段末的河水位。然后用新的河水位再代入地下水数值方程，从而实现交互计算功能。河床的渗透性由下式进行计算：

图2 地下水流数值方程与河水流的数值方程联立求解示意图Fig.2 Sketch map of simultaneous solution of groundwater and river flow numerical equations

其中：C—河床水力传导系数，m/d；k—河床渗透系数，m/d；w—河床宽度，m；t—河床厚度，m；L—上下游距离，m。

2.6 所用数据资料

模型所用数据资料包括：①流域月降雨量与流量（图3）；②地下水井丰、枯水期水位监测资料与河流中游的流量观测资料（图4）；③地下水开采量资料，包括生活用水、工业用水和农业用水，总开采量共1 014.7万m3/a，折合为2.78万m3/d。

3 结果

图3 现状年累计月降雨量流量过程线Fig.3 Current year cumulative monthly rainfall flow process line

3.1 模型识别验证

首先通过模型模拟来识别研究区2017年1—3月枯水期地下水水位（图5），由拟合误差统计数据可知（表2），模拟结果与实测水位之差小于0.5 m 的占80%以上，且相对误差大多小于10%。再验证2017年8—10月丰水期水位（图6），计算结果与实际水位之差小于0.5 m 的占70%左右，且相对误差大多小于10%。由模拟水位图和模拟数据可见，模拟水位线趋势与实际流场较吻合，模拟丰枯水期的海水入侵范围也与实际情况基本一致，证明本模型准确可靠。丰枯水期源汇项数据由丰枯水期降雨量及河流流量给定，研究区丰、枯水期降雨入渗补给分别约为11.3、9.68万m3/d，丰，枯水期河流补给地下水量分别约2.99、2.18万m3/d。具体各参数分区源汇项见表3。

采用2017年平水年中游流量进行识别和验证。图7为实测流量与模拟流量拟合，模拟结果较为贴合实际。

3.2 现状开采量评价

根据已经建立的数值模型进行流域水均衡计算。首先根据《青岛市胶南县吉利河-白马河地区供水水文地质勘察报告》中研究区降雨频率分布曲线，计算保证率分别为50%、75%、95%条件下的大气降雨量，根据研究区的降雨入渗分区换算为不同的降雨补给量，然后作为不同的降雨补给条件赋给模型。再根据地下水利用情况，评价研究区范围内地下水开采量为1 014.7万m3/a，研究区面积46 km2，折算每日开采强度为0.000 6 m/d。通过模型计算研究区保证率分别为50%、75%、95%条件下的大气降水入渗量、蒸发量、灌溉回渗量、开采量见表4。由表3可以进一步计算得到各种保证率条件下的地下水允许开采量。按照目前模型，每年开采量为2 847.0万m3，在平水年就已经失去了平衡，开采量过大，需要进行压缩。

图4 长观孔分布图Fig.4 Measured wells distribution map

表2 水位拟合误差统计数据Tab.2 Water level fitting error statistics

表3 模型识别与验证期不同分区源汇项Tab.3 Different partition source and sink items during model identification and verification periods 万m3/d

4 分析讨论

本次模拟采用扩散波方程作为地下水流数值方程的边界条件，实际上是耦合了地下水数值模型和地表水数值模型。与传统的地下水数值模型相比，本模拟考虑了河水流量和水位变动情况（天然情况下，河水位和流量受地下水傍河开采影响而发生变化）。如果单纯使用第一类定水头边界模型进行傍河地下水开采量评价，水位与河流实际情况不符，会造成计算结果误差较大。为了验证本方法的合理性，使用已知水头边界再次建立研究区傍河开采地下水评价数值模型，将已知水头边界模型和扩散波方程边界模型的评价结果进行如下对比讨论。

图5 模型识别期实测水位与计算流场Fig.5 Measured water level and calculated flow field during model identification period

图6 模型验证期实测水位与计算流场Fig.6 Measured water level and calculated flow field during model verification period

图7 实测流量与模拟流量拟合Fig.7 Measured flow and simulated flow fit

4.1 耦合模型与定水头边界模型地表水地下水转化量对比讨论

由模型丰枯水期地表水地下水耦合模型与定水头模型均衡对比（表5），比较两者的河流补给和排泄量：①丰水期时耦合模型河流补给地下水量2.99万m3/d，远小于定水头边界的7.53万m3/d，耦合模型的河流排泄量2.50万m3/d小于定水头模型的5.67万m3/d；②枯水期时耦合模型河流补给地下水量2.18万m3/d也远小于定水头边界的9.99万m3/d，丰枯水期海水入侵量两个模型均基本未变，但耦合模型的海水入侵量0.4万m3/d大于定水头模型的0.21万m3/d，原因是耦合模型中河水位是变动的，低于定水头河流边界，导致计算的海水入侵量增大。但就数量上，海水入侵量的变化在模型中体现轻微。

由表5比较可见，第一类定水头边界模型河水补给量远大于耦合模型的计算结果，丰枯水期河水补给量均大于河水排泄量，该模型在丰枯水期地表水均补给地下水，当枯水期降雨补给量少时，定水头边界对地下水的补给量反而更大，显然不符合白马河季节性河流特征。对比发现，耦合模型地表水地下水转化与实际情况更为吻合，丰水期地表水补给地下水，枯水期地下水补给地表水。按耦合模型扩散波方程，当河流断流时，河水不再补给地下水开采。

表4 不同保证率下的补给量和排泄量Tab.4 Replenishment and excretion at different guarantee rates 万m3/年

表5 地表水地下水耦合模型与定水头模型均衡对比Tab.5 Equilibrium contrast between surface water and groundwater coupling model and fixed head model万m3/d

4.2 耦合模型与定水头边界模型模拟开采量对比讨论

该区地下水地表水存在密切水力联系，联系通道以河床为主，地下水主要通过河流边界进行水力交换，因此最大允许开采量和河流年径流量必然存在密切联系。

由表5可知，丰水期耦合模型井开采量8.25 m3/d小于第一类定水头边界模型的10.80万m3/d，枯水期时耦合模型井开采量6.73万m3/d也小于第一类定水头边界模型的8.8万m3/d。

第一类定水头边界模拟出的最大允许开采量显然比耦合模型模拟结果大得多，模拟结果为3 942万m3/a，参照《白马河地区供水水文地质勘察报告》，研究区平水年地下水可开采量2 870万m3/a与耦合模型模拟的开采量2 733万m3/a较为接近。另外，按第一类定水头边界计算最大允许开采量，与降水量乘上径流系数的径流总量3 120万m3/a相比，第一类定水头边界最大允许开采量3 942万m3/a，已经超过年径流量，这间接说明了第一类定水头边界的不合理性。综上，耦合模型相对于第一类定水头边界模拟结果更符合实际。

4.3 机理分析

当河流模型以扩散波方程为基础，在枯水期河流部分河段断流时，模型模拟出局部河段不再补给地下水开采，计算结果与实际相符。当河流按第一类定水头边界概化时，河流对地下水一直存在补给，进而使模拟地下水开采量偏大，与实际不符。其次，即使枯水期河水没有断流，但其河宽非常小，在模型剖分精度上无法满足精细的第一类定水头边界，而扩散波方法将河道细致概化到河床水力传导系数（C）、河床渗透系数（k）、河床宽度（w）及河床厚度（t）等参数与实际河道情况更为贴近，较好地解决了这一问题。本方法对于山区河谷、滨海河流区等傍河水源地开采量评价具有推广借鉴意义。

5 结论

考虑到白马河流量水位在丰枯水期的巨大差异，应用河流扩散波方程，模拟河流丰枯水期流量，使河流概化与实际情况贴近，提高了水源地地下水资源量计算精度。结果显示，按耦合模型和定水头边界模型评价本区地下水可开采量分别为2 733万m3/a、3 942万m3/a，而本区的年径流总量为3 120万m3/a，传统的已知水头边界模型评价的地下水可开采量超过了流域年径流量，而耦合模型评价的结果与研究区平水年地下水可开采量2 870万m3/a更为接近。因此耦合模型对于河水位变化较大的河流边界模型更适用。