黑河地表-地下水交互作用对水量调度的影响

王 玥，徐 磊，赵建世

（1.中国灌溉排水发展中心，北京 100054；2.清华大学，北京 100084）

1 研究背景

黑河发源于祁连山南部向北流至下游的绿洲居延海，潜在蒸发量为1325 mm/a，而降水量仅为190 mm/a，表明该地区气候极为干旱。含水层的地下水位从上游的约200 m到下游的约5 m不等，地表水和地下水交换频繁［1-3］。黑河流域人口约95%聚集在农业发达的中游地区，中游灌溉行为与下游生态系统的关系对整个流域的可持续发展影响重大［4-6］。

为合理利用黑河流域水资源，1992年国务院批准了黑河干流（含梨园河）分水方案。为提高分水方案在实际管理中的可操作性，水利部黄河水利委员会编制了《黑河流域水量分配方案》（简称“97方案”），规定了不同的莺落峡径流量条件下的正义峡下泄水量目标［7］。在丰、平、枯和特枯年份（对于25%、50%、75%和99%频率），上游莺落峡对应径流量为17.1、15.8、14.2和12.9亿m3，下游正义峡的分水控制目标分别为10.9亿、9.5亿、7.6亿和6.3亿m3。其它来水条件下，分配正义峡下泄水量按以上保证率水量直线内插求得。在黑河分水方案实施的过程中，由于分水曲线自身技术特点，出现了丰水年份的分水目标实现压力大、“越丰水越难完成任务”的现象：在枯水年一般比较容易完成水量下泄目标，但在丰水年下泄水量经常达不到目标，尤其当出现来水前枯后丰的年份调度困难更大，为调度管理工作带来了技术挑战。分析这一现象的原因及应对措施是黑河流域管理中亟需回答和解决的问题。

地表-地下水模拟是分析这一问题的重要技术手段。地表-地下水系统模拟的数值方法从1980年代开始建立并迅速发展［8］，目前已经成为流域水循环的主要手段［9-11］。在众多模型方法中，美国地质调查局开发的MODFLOW为地表-地下水模拟提供了重要的平台［12-13］，将不同类型的地表水模拟模型和MODFLOW进行耦合是应用最广泛的方法［14-15］。另外，MODFLOW也开发了地表水模块Stream子程序，为地表-地下水系统模拟提供了便捷的工具。

本文基于MODFLOW平台构建了黑河中游地表-地下水模拟模型，通过情景模拟分析，对黑河流域“越丰水越难完成任务”这一问题背后的物理机制进行了模拟和解释，相关研究结论可以为分水方案的实施提供科学依据。

2 地表-地下水数值模拟模型

本研究采用Visual MODFLOW软件平台构建了黑河中游地表-地下水联合模拟模型，其中地表水模块采用平台自带的Stream子程序。模型可以用于模拟黑河干流莺落峡至正义峡区间的地表水部分以及河道水与地下水的转化关系，该模块还可模拟沿程各灌区从河道取水，并计算正义峡河道流量。

2.1 Visual MODFLOW地下水模拟原理 Visual MODFLOW由加拿大Waterloo水文地质公司基于原MODFLOW软件开发，是目前国际上主流的三维地下水流和溶质运移模拟的标准可视化专业软件。这个完整的集成软件平台将水流模型MODFLOW、平面和剖面流线示踪模型MODPATH和溶质运移模型MT3D同最直观强大的图形用户界面结合在一起［16］，分为输入、运行和输出3大模块。

MODFLOW的地下水流运动规律描述采用有限差分的方法，首先在空间和时间上对研究区进行离散，然后建立研究区域网格点的水平衡关系方程，将相关方程联立形成线性方程组，最后采用迭代的方式进行求解，得到每个网格地下水状态量。在不考虑水的密度变化的前提下，多孔介质中地下水在三维空间中的流动可用偏微分方程表示为：

式中：W为单位体积流量，代表来自源或流进汇的水量，d-1；h为水头，m；Kxx、Kyy、Kzz分别为沿三个坐标轴方向的水力传导系数，m/d；Ss为多孔介质贮水系数，m-1；t为时间，d。

2.2 地表水模拟模块及其与地下水的耦合 MODFLOW的Stream模块可以为地表-地下水循环模拟、人类活动影响和地下水资源评价等提供技术支撑［17-19］。本研究采用Stream子程序模拟黑河中游段莺落峡至正义峡的河道水流过程以及地表-地下水转化关系，同时还可以模拟中游各灌区从河道取水带来的河道流量和地下水改变，并最终计算下游正义峡的流量。需要注意的是，Stream模块对地表水模拟的精度不如专业的地表水模拟模型，未来可以将更加精确的地表水模型与MODFLOW结合，改进模拟的精度。

Stream子程序将河道水流与地下水的交换概化为一维流动，其中的重要参数河床水力传导系数可以表示为：

式中：K为河床淤积层渗透系数；L为计算单元河段长；N为河床宽度；M为河床淤积层的厚度；CRIV为河床的水力传导系数［20］。

假设河道与地下含水层间的水头损失仅受到河床淤积层的影响，且河床淤积层处于完全饱和状态，那么河道与地下含水层交换水量可以表述为：

式中：HRIV为河道水位；CRIV为地下含水层河道间的水力传导系数；hi,j,k为河流所在单元的水头值；QRIV为地下含水层与河道之间交换量（河水补给地下水时为正）。当地下含水层水头下降至河床淤积层之下并产生非饱和带的情况下，设河床底部高程为RBOT，则河道与地下含水层的水量交换可以表达为：

从式（5）可以看出，当河道水头HRIV低于RBOT时，河道和地下含水层的水量交换为负值（地下水补给河水）。

同时，Stream子程序可模拟灌区从河道引水的量，灌区引水量可以通过设定河段起点与终点流量来代表。另外，由于Stream是耦合在MODFLOW软件包中的一个子模块，因此可以很好地处理地表地下水之间数据交换和相互作用，计算方便、效率高。

3 黑河中游地表-地下水模型构建

3.1 研究区域概况 黑河流域共发育有30多条河流，可分为东部、中部、西部三个相互独立的子水系（见图1）。本次研究区为东部的黑河干流中游（含梨园河及沿山20余条小支流，如图1中黑框部分所示），研究区域总面积为8717km2。中部子水系包括马营河等小河流，都是浅山短流，流量小且出山后很快被消耗尽，对研究区水循环无影响；西部子水系为讨赖河流域，包括讨赖河和洪水河等，汇合后流入下游金塔盆地，对研究区无影响。

研究区地质构造复杂，边缘被断裂带控制，覆盖第四纪松散地层可存储大量地下水，基底岩性是第三系泥岩；山前堆积物主要是河流冲洪积物，含水层以砂砾石为主，向北岩性逐渐过渡到细土平原区含亚砂土和亚黏土的含水层，含水层逐渐变浅、渗透性减弱。山前南半部的洪积扇最上层为大厚度单层潜水区，下游为细土平原区，逐步由单一潜水含水层区变为多层的潜水-承压水层混合区［21-22］。黑河流域地表地下水交换频繁，出山口流出的地表径流在南部山前大量渗漏补给地下，地下水到了张掖盆地中北部的细土平原带，因含水层渗透性减弱并且地表河床下切，导致地下水大量溢出地表形成地表径流和泉水，其中泉水溢出带主要分布在洪积扇缘和细土平原地，而中游用水灌溉后又有水量补给到地下含水层，形成多次地表-地下交换［21-22］。山前扇群带地下水埋深变化剧烈，扇顶为200～500 m，近山侧为100～300 m。扇中张掖盆地地下水埋深约在50～100 m之间，酒泉东盆地地下水埋深在100～250 m之间。扇缘张掖盆地地下水埋深约为10～20 m，而酒泉东盆地则在80～200 m之间。细土平原区地下水埋深变化也较大，其中张掖盆地大多在5 m以内，而酒泉东盆地南半部则高达10～50 m，而北半部只有1～5 m。

黑河干流水系多年平均出山口径流量为25.11亿m3，其中莺落峡站多年平均径流量为16.19亿m3。黑河流域矿化度小于2 g/L的多年平均地下水资源量为21.76亿m3，其中山丘区为9.36亿m3，平原区为20.32亿m3，二者重复量为7.92亿m3。全流域共有中小型水库57座，总库容为2.76亿m3；干流灌区引水工程96处，设计引水能力约269.8 m3/s，其中直接从干流引水的口门有36处，设计引水能力为223.8 m3/s；机电井11 076眼，其中配套机井9771眼，年提水量5.81亿m3。流域灌溉面积为465.67万亩，其中林草灌溉面积为95.91万亩，农田灌溉面积为369.75万亩，30万亩以上的大灌区有8处共计301.44万亩。

3.2 黑河中游地表-地下水模型构建

3.2.1 模拟范围及含水层条件 本文研究的空间范围东起山丹县祁家店水库、西至酒泉东盆地清水车站，南北以山前为界。由于承压区没有稳定的隔水层，研究区各含水层有不可分割的水力联系，是一个较为完整的地下水盆地（如图1所示）。近几十年来，由于地下水的开采量不断增加，使得承压含水层的连通性增加，不同层之间的水力联系较历史上更为紧密，导致上下层水位相差很小甚至基本相同。基于此，本文将研究区含水层概化为非均质的各项同性单层潜水层，同时由于区内潜水含水层分布广、厚度大，将研究区地下水运移视为二维非稳定流。

图1 黑河流域水系及中游研究区

研究区模型采用正交网格剖分，剖分格距为ΔX=ΔY=1 km，利用GIS中大地坐标的底图导入MODFLOW圈定模拟区域范围。东西向剖分为172（km）格，南北向剖分为155（km）格，垂向上为单层。河流宽度和深度由模型根据地形和流量计算。

3.2.2 模型边界条件 模型边界条件包括侧向边界和垂向边界，具体设置如下。

侧向边界条件：研究区域周边二类流量边界，其中南、东、北边界均为弱透水边界，有较弱的侧向补给［22］：南边界山前断裂带为基岩裂隙水和沟谷潜流补给，而东部民乐和山丹断面以及北部山前断面均为区外侧向流入边界。西边界由于是酒泉西盆地与酒泉东盆地分水岭，可视为零通量边界。

垂向边界条件：上边界是天然潜水面，下边界是侏罗系和第三系泥岩以及砂质泥岩等隔水边界，因此可忽略顶托越流补给。

3.2.3 模拟时段与参数 模型参数率定和验证的时段为1995年1月—2000年12月，计算步长为月，共计72个时段。本文采用ΔX=ΔY=1 km格距的正交网格进行研究区的网格剖分，河床宽度N数据采用代表性断面的实测平均宽度。模拟区域总面积为8717 km2，东西方向可剖分为172格，南北方向可剖分为155格，垂直方向为单层。跟据《黑河流域水资源评价利用及分配方案研究》中提供的不同岩性土的给水度，以及“九五”和“十五”期间黑河中游地区的研究资料提供的水文地质参数，确定研究区导水系数大约在500～14000 m3/d之间，而给水度的范围大约在0.08～0.2之间。本文利用Stream子程序包模拟黑河干流中游段河水与地下水的频繁转化关系，河床底板高程和河床厚度及水力传导系数等相关参数参考文献［23-25］，并通过模型率定与验证进行了调整，莺落峡的实测流量和各灌区实测引水流量来自黑河流域水资源管理局统计资料。1995年初的地下水初始水位根据地下水观测井数据空间插值得到，经过模型率定过程将模型模拟的地下水位与观测井实测值对比，2000年末研究区地下水埋深与观测值拟合良好，以此作为模型计算分析时段2001年初的初始水位值和初始流场。

3.2.4 模型率定与验证 根据上述资料给定研究区的初始参数，经过多次调试和优选，以各观测孔处的模拟水头值与实测值最接近为原则来率定模型参数。模型采用1995—2000年间的实测资料（包括气象水文、地下水位观测资料和水资源利用数据等）对地表-地下水模型进行参数率定和模型验证。注意的对比数据包括：（1）研究区17个观测孔水位动态过程；（2）下游正义峡断面的流量过程。

统计结果表明，模型率定期（1995—1997年），观测井水位值模拟误差在1.0m以内的数据占98%（其中误差小于0.5m的占61%，误差在0.5～1.0m之间的占37%），Nash-Sutcliffe效率系数为0.87；在模型验证期（1998-2000年），观测井水位值模拟误差在1.0m以内的数据占96%（其中误差小于0.5m的占51%，误差在0.5～1.0m之间的占45%），Nash-Sutcliffe效率系数为0.83，表明模型参数率定和验证效果良好。选取其中两个具有代表性的观测井率定和验证期内的拟合过程如图2所示，率定和验证期末（2000年12月末）的地下水位等高线以及所有观测井的对比情况如图3所示。

图2 47号观测井和65号观测井的拟合结果

黑河下游正义峡断面的流量过程模拟值与实测值对比结果如图3所示。由于实测数据只有观测井的水位值，图3（a）只给出了模拟流场，无法给出模拟流场和实测流场的对比；图3（b）给出了率定和验证期末观测井模拟和实测水位的对比，表明观测井模拟效果良好。图4给出了下游正义峡断面的模拟流量和实测流量对比，其中模拟流域为Stream模块输出值，模型率定期（1995—1997年）的Nash-Sutcliffe效率系数为0.82，相关系数为0.88；模型验证期（1998—2000）的Nash-Sutcliffe效率系数为0.84，相关系数为0.89，表明模型对地表和地下水交换作用的模型效果整体上较好。

图3 率定和验证期末地下水位等高线与观测井水位对比

图4 正义峡断面流量过程拟合结果对比

4 地表-地下水交互作用及其影响分析

为了分析分水方案调度的情况，利用上述模型模拟了黑河统一调度的实际情景，时段为2001年1月至2006年12月。由于上述模型率定验证的结束时段为2000年12月，此时黑河中游的地下水埋深与观测值拟合良好，可以作为模型模拟的初始水位值。

2001—2006年的模拟结果如表1和图5所示。从结果可以看出，“莺落峡来水越丰时越难完成分水指标”这一现象确实存在。2001—2006年的模拟结果和1992年国务院批准了黑河干流（含梨园河）分水方案中列出的泄水指标分别相差-0.22亿、-0.90亿、-1.76亿、-0.91亿、-1.59亿、-1.17亿m3，在来水较丰的2003、2005和2006年，差异都超过了1.0亿m3。而产生这一现象的主要原因，是由于丰水年河道水位相对较高、地下水向的补给量显著减少，造成“越丰水越难完成分水目标”现象。对比2003年和2001年，地下水向河道的净补给量分别为1.05亿和2.98亿m3，丰水的2003年相对于枯水的2001年减少了1.93亿m3，在这种情况下即使压缩中游用水也无法完成正义峡分水指标。这一现象体现了地下水作为一个天然调蓄池对水资源的调丰补枯作用，对于地表水的调度和管理具有不可忽视的影响。

表1 分水情景模拟 （单位：108m3）

图5 分水情景对比分析

基于上述分析，建议对现行的1992年分水方案进行修正，因为这一方案基本上按照线性关系来考虑不同丰枯年份水量分配，没有考虑的地下水调蓄作用对实际调度过程的影响。如果将这一因素考虑在内，建议将丰水年向正义峡的下泄比例降低5%～10%、平水年的下泄比例保持不变、枯水年的下泄比例增加5% ～10%，这样即有利于方案的实际执行，也有利于中游的用水保障和下游的生态保护。具体的调整方案，建议综合考虑各方面因素进行详细的论证。另外，如果可以在上游修建控制性水库，进行多年丰枯调节，将更有利于整个分水方案的执行和流域的可持续发展。

5 结语

本研究构建了黑河中游地表-地下水模拟模型，利用上述模型模拟了2001年1月—2006年12月黑河统一调度的实际情景，分析了地表-地下水交互作用对水量调度的影响。结果表明，向正义峡下泄的实际水量与1992年分水方案指标相比，在来水较丰的2003、2005和2006年，差异均超过了1.0亿m3，产生了“越丰水越难完成分水目标”的现象。丰水的2003年和枯水的2001年对比表明，由于丰水年份河道水位相对较高，地下水向河道补给量会有所减少了1.93亿m3，进而导致上游向正义峡的下泄比例减少，难以完成分水指标。因此，建议对分水方案进行调整，将丰水年向正义峡的下泄比例降低5%～10%、平水年的下泄比例保持不变、枯水年的下泄比例增加5%～10%，同时在有可能的情况下建设上游多年调节水库，促进分水方案的科学实施和流域水资源的高效利用。

鉴于调度问题的复杂性，后续研究可以在本文地表-地下水模型的基础上，对不同水文条件和经济社会发展情景下的分水方案进行综合对比分析，进一步提出合理分水指标调整方案和调度管理实施方案。