摘要：为了提高地表地下双重人类活动干扰下的流域洪水模拟精度，基于调蓄水库概念建立了考虑地表水利工程和地下水超采影响下的网格化降雨径流模拟方法，通过在栅格新安江模型基础上引入地表地下双调蓄模块，构建了栅格新安江-地表地下双人工调蓄分布式水文模型，将该模型同半分布式新安江-海河模型在海河典型流域进行应用比较。结果表明：栅格新安江-地表地下双人工调蓄分布式水文模型在清水河流域的模拟精度高于新安江-海河模型;基于栅格的双调蓄结构能根据水利工程控制区域和地下水埋深的实际空间分布对地面和地下径流人工调蓄过程进行准确模拟;栅格新安江-地表地下双人工调蓄分布式水文模型对自由水库蓄水容量参数的取值较新安江-海河模型更加合理;栅格新安江-地表地下双人工调蓄分布式水文模型能有效提高清水河流域的洪水模拟精度。
关键词：水利工程;地下水超采;降雨径流;栅格新安江-地表地下双人工调蓄分布式水文模型;洪水预报;海河流域
自20世纪50年代以来，我国兴建了大量水利工程[1]。水库、塘坝等蓄水工程阻拦了部分降雨和地表径流，对场次洪水尤其是中小河流洪水具有明显的蓄洪削峰作用，从而改变了下游径流的时空分布[2]。地下水是水资源的重要构成，由于地表水资源匮乏，地下水在我国北方尤其是华北地区是主要的供水水源，地下水的开采和利用成为支撑生产生活的重要保障[1]。1980年以后，随着经济社会高速发展，人类活动对下垫面的影响日益显著[3]。地下水的持续超采造成地下水位的持续下降，包气带随之不断增厚，同时大范围地下水漏斗的形成增大了降雨径流过程的入渗损失量，从而使流域径流量大幅减少。水利工程和地下水超采等人类活动极大地改变了流域产汇流规律，因此提出相应的模拟方法来科学刻画这种特殊降雨径流关系变得尤为重要。初步研究已经取得了一些进展：王中根等[4]构建了海河流域地表水SWAT与地下水MODFLOW松散耦合模型，模拟结果符合流域管理的应用需求;张建中等[5]建立了模拟地下水埋深变化和灌溉耕地产流特点的蓄满-超渗产流水文模型，在牤牛河等7个华北平原区典型流域进行应用并得到较高的模拟精度;李致家等[6]在新安江模型基础上增加了对地表和地下拦蓄的模拟，构建了半分布式新安江-海河模型，有效地提高了海河典型流域的洪水模拟精度。以上研究为构建适合地表地下人为干扰严重的流域的水文模型打下基础，但是对地表地下人为干扰的刻画还不够精细。
与湿润流域相比，半湿润流域降雨时空分布不均，下垫面空间异质性对径流的影响更大，洪水破坏性也更大[7-13]。目前，北方流域的水文模型多为经验模型或半分布式水文模型[14,15]，难以精细刻画流域降雨空间分布不均的特点、土壤植被、地形地貌等下垫面因子的空间异质性[16-18]以及地表水利工程控制区域和地下水埋深的空间分布。针对以上问题，本文在栅格新安江模型[19,20]基础上，引入地表地下双人工调蓄模块，提出了考虑水利工程和地下水超采影响下的降雨径流模拟方法，在此基础上构建栅格新安江-地表地下双人工调蓄分布式水文模型(GriddedXin'anjiang-DualAnthropogenicAbovegroundandUndergroundRegulationDistributedHydrologicalModel，GXAJ-DAR模型)，并检验了该模型在受人类活动影响严重的华北典型流域的应用效果。
研究方法
人类活动影响下的降雨径流模拟
本文研究的人类活动影响主要包括地表水利工程运行与地下水超采，将建立相应的方法模拟库塘闸坝对降雨和地表径流的直接拦蓄作用以及地下水漏斗对产汇流的间接拦蓄效应。
地表人工调蓄模拟
基于水利普查资料，获取流域库塘闸坝蓄水能力及其空间分布信息，采用“地表调蓄水库”方法刻画库塘闸坝对地表径流的阻拦作用[6,21]。具体计算方法如下：
利用ArcGIS提取流域内各库塘闸坝的控制区域，控制区域内栅格产生的部分地表径流将被地表调蓄水库所阻拦，本研究暂不考虑水利工程的调节性能。控制区内的每一个栅格点都会受到水库的影响，但具体影响机制较为复杂，因此对其做概化处理，即以水利工程总库容的一半作为其蓄水容量，将该蓄水容量平均到其控制区内的所有栅格，进而得到控制区域每一栅格的“水库蓄水容量”，即R=10VcNA(1)式中：R为单个水利工程控制区域各栅格的水库蓄水容量，mm;V为水利工程防洪库容，万m3;N为控制区域的栅格数量;A为栅格面积，km2。
以模拟时段初的土壤含水量为指标计算地表调蓄水库的初始蓄水量，即æWöRc0=Rcç1-W÷(2)èmø式中：Rc0为栅格初始水库蓄水量，mm;W为模拟时段初栅格包气带土壤含水量，mm;Wm为栅格包气带最大蓄水容量，mm。
网格模拟单元在一次洪水过程中的地表径流拦蓄量估算公式为[1]R=min(RRc0,R)(3)式中：Rv为栅格地表径流拦蓄量，mm;Rs为洪水过程中栅格地表径流产流量，mm。
非水利工程控制区域的计算栅格产流量不受地表调蓄水库控制;控制区域内进入地表调蓄水库的地表径流将消耗于蒸散发或用于灌溉等。当地表调蓄水库蓄满之后，产生的地表径流将不再受到拦蓄影响。
地下人工调蓄模拟
基于流域的数字高程地图、地表水位和地下水位资料，采用“地下调蓄水库“方法刻画产流对地下水超采引起包气带变化的非线性响应以及地下水漏斗对产汇流的间接调蓄效应[6,21]。具体计算过程分为两个步骤：
当下渗雨量满足包气带张力水蓄水容量，多余的雨量将进入地表自由水蓄水库。由于地下水超采的影响，一部分自由水库中的水将向下渗漏以补充地下水，抬升地下水面。该渗漏过程采用格林-安普特下渗公式[22,23]进行刻画，即æyDqöf(t)=Kç1+F(t)÷(4)èø式中：f(t)为栅格自由水库渗漏能力，mm/h;K为栅格饱和水力传导度，mm/h;y为栅格湿润锋处土壤吸力，mm;Dq为栅格土壤饱和含水率与初始时刻含水率之差;F(t)为栅格累计渗漏量，mm;t为时间。
自由水库的渗漏水量将进入地下调蓄水库，累计下渗水量超过地下调蓄水库阈值S0的部分以地下水出流系数Kg出流形成地下径流。其中，地下调蓄水库阈值S0的计算方法如下：ì(Zï(-Z)´n,Zr+DH>Zi')DH>Zg'+DH'S0=íZr+DH-Zg-DHn,Zr+DH>Zg+DHZi+DH(5)ï0,Z+DH'>Z+DH式中：Zr为流域出口河道水位，m;Zi为流域出口所在栅格的相对不透水层高程，m;Zg为栅格多年平均地下水位高程，m;DH为栅格与流域出口点的相对高程差，m;DH'为模拟时段初地下水位站几何中心位置的地下水位高程与其多年平均地下水位高程之差，m;n为栅格土壤孔隙度。
GXAJ-DAR模型
基于上述产流模拟方法，在栅格新安江模型(GXAJ)[24]基础上构建适用于受人类活动影响流域降雨径流模拟的GXAJ-DAR模型。不透水面上的直接径流与表层自由水库蓄满后产生的饱和地面径流汇合，水利工程控制区内栅格产生的地面径流经地表调蓄水库拦蓄，超出蓄水容量部分以地面径流形式汇流至流域出口;表层自由水库中的部分水量将补充地下调蓄水库，达到阈值后，地下调蓄水库将以某一系数出流形成地下径流;模型其余部分与GXAJ模型一致，如图1所示。与基于泰森多边形的半分布式新安江-海河模型(XAJ-HH)相比，GXAJ-DAR模型考虑了流域降雨空间分布不均的特点、土壤植被、地形地貌等下垫面因子的空间异质性以及地表水利工程控制区域和地下水埋深的空间分布，可以更精细地刻画人类活动影响下的降雨径流过程。
GXAJ-DAR模型在GXAJ模型基础上新增参数如下：稳定下渗率fc、土壤孔隙度n、湿润锋处土壤吸力y、水利工程防洪库容Vf、流域出口位置相对不透水层高程Zi。fc、n、y可根据土壤性质得到，V可根据实际水利工程信息获得，Z采用人工优选法率定得到。姚成[25]研究发现，在DEM分辨率为1km的栅格内认为蓄水容量分布均匀可以满足洪水模拟精度要求，因此本研究中的GXAJ-DAR模型将不考虑参数和变量在栅格内部的空间变异性问题。
以下将GXAJ-DAR模型与XAJ-HH模型在研究流域进行应用检验，比较分析2个模型的模拟效果。
模型率定与评估方法
GXAJ-DAR模型的集总式参数和XAJ-HH模型的参数均采用人工优选法在经验取值范围内进行率定。依据GB/T22482—2008《水文情报预报规范》，选择径流深合格率、洪峰合格率、确定性系数3个指标来综合评判各模型模拟结果(其中确定性系数大于0.5时假定为合格)，并对本文提出的人类活动影响下的降雨径流模拟方法进行适用性评价。
研究流域与数据资料
研究流域
海河流域为受地下水超采严重影响的典型流域，主要分布于我国渤海以西，黄河以北，位于35°N～43°N、112°E120°E之间[26]，流域面积达31.78万km2，平均长度450km，其中60%为山区，40%为平原。多年观测数据显示，海河流域年降水量仅548mm，由于气候、地形等多种因素的影响，海河流域夏季暴雨集中，冬春雨雪稀少，年降水量分布具有很强的地域性及季节差异性，同时还伴随着连续枯水及连续丰水的周期性变化。
选择河北省保定市清水河流域作为研究流域，该流域位于海河流域大清河水系，属于半湿润流域。流域出口北辛店水文站多年平均年径流量0.29亿m3(1971—2010年)，多年平均年降水量457.1mm(1971—2010年)，最大年降水量833.7mm。实测最高洪水位19.70m，最大洪峰流量为1966年710m3/s，其次是1989年686m3/s(由于洪水漫滩，实际洪峰流量估计值为1500m3/s)。
清水河流域面积1650km2，流域内设有雨量站11处，其中有完整水雨情资料的水文站1处，雨站6处，设有地下水位测站17处，其中有完整地下水埋深资料的3处(图2)。北辛店水文站上。
