栅格新安江-地表地下双人工调蓄分布式水文模型

张 珂，张企诺，陈新宇，李致家，黄鹏年，姚 成，王 晟，周佳奇

(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏 南京 210098；2.河海大学长江保护与绿色发展研究院，江苏 南京 210098; 3.河海大学水文水资源学院，江苏 南京 210098；4.中国气象局-河海大学水文气象研究联合实验室，江苏 南京 210098;5.南京信息工程大学水文与水资源工程学院，江苏 南京 210044)

自20世纪50年代以来，我国兴建了大量水利工程[1]。水库、塘坝等蓄水工程阻拦了部分降雨和地表径流，对场次洪水尤其是中小河流洪水具有明显的蓄洪削峰作用，从而改变了下游径流的时空分布[2]。地下水是水资源的重要组成部分，由于地表水资源匮乏，地下水在我国北方尤其华北地区是主要的供水水源[1]。1980年以后，随着经济社会高速发展，人类活动对下垫面的影响日益显著[3]。地下水的持续超采造成了地下水位的持续下降，包气带随之不断增厚，同时大范围地下水漏斗的形成增大了降雨径流过程的入渗损失量，使流域径流量大幅减少。水利工程和地下水超采等人类活动极大地改变了流域产汇流规律，因此提出相应的模拟方法来科学刻画这种特殊降雨径流关系变得尤为重要，相关研究已经取得了一些进展。王中根等[4]构建了海河流域地表水SWAT与地下水MODFLOW松散耦合模型，模拟结果符合流域管理的应用需求；张建中等[5]建立了模拟地下水埋深变化和灌溉耕地产流特点的蓄满-超渗产流水文模型，在牤牛河等7个华北平原区典型流域的应用取得了较高的模拟精度；李致家等[6]在新安江模型基础上增加了对地表和地下拦蓄的模拟，构建了半分布式新安江-海河模型，有效地提高了海河典型流域的洪水模拟精度。以上研究为构建适用于地表水利工程和地下水超采影响流域的水文模型打下了基础，但是对上述地表和地下两种人为干扰的刻画还不够精细。

与湿润流域相比，半湿润流域降雨时空分布不均，下垫面空间异质性对径流的影响更大，洪水破坏性也更大[7-13]。目前，北方流域的水文模型多为经验模型或半分布式水文模型[14-15]，难以精细刻画流域降雨空间分布不均的特点和土壤植被、地形地貌等下垫面因子的空间异质性[16-18]以及地表水利工程控制区域和地下水埋深的空间分布。针对以上问题，本文在栅格新安江模型[19-20]基础上，引入地表地下双人工调蓄模块，提出了考虑水利工程和地下水超采影响下的降雨径流模拟方法，在此基础上构建栅格新安江-地表地下双人工调蓄分布式水文模型(Gridded Xin’anjiang-dual anthropogenic aboveground and underground regulation distributed hydrological model，GXAJ-DAR模型)，并检验了该模型在受人类活动影响严重的华北典型流域的应用效果。

1 模型构建

1.1 人类活动影响下的降雨径流模拟

本文研究的人类活动影响主要包括地表水利工程运行与地下水超采，将提出相应的方法模拟库塘闸坝对降雨和地表径流的直接拦蓄作用以及地下水漏斗对产汇流的间接拦蓄效应。

1.1.1地表人工调蓄模拟

基于水利普查资料，获取流域库塘闸坝蓄水能力及其空间分布信息，采用“地表调蓄水库”方法刻画库塘闸坝对地表径流的阻拦作用[6,21]。具体计算方法如下：

a.利用ArcGIS提取流域内各库塘闸坝的控制区域，控制区域内栅格产生的部分地表径流将被地表调蓄水库所阻拦，本研究暂不考虑水利工程的调节性能。控制区内的每一个栅格点都会受到水库的影响，但具体影响机制较为复杂，因此对其做概化处理，即以水利工程总库容的一半作为其蓄水容量，将该蓄水容量平均到其控制区内的所有栅格，进而得到控制区域每一栅格的“水库蓄水容量”：

(1)

式中：Rc为单个水利工程控制区域各栅格的水库蓄水容量，mm；V为水利工程防洪库容，万m3；N为控制区域的栅格数量；A为栅格面积，km2。

b.以模拟时段初的土壤含水量为指标计算地表调蓄水库的初始蓄水量：

(2)

式中：Rc0为栅格初始水库蓄水量，mm；W为模拟时段初栅格的包气带土壤含水量，mm；Wm为栅格包气带张力水蓄水容量，mm。

c.网格模拟单元在一次洪水过程中的地表径流拦蓄量估算公式为[1]

Rv=min(Rc-Rc0,Rs)

(3)

式中：Rv为栅格地表径流拦蓄量，mm；Rs为洪水过程中栅格地表径流产流量，mm。

非水利工程控制区域的栅格产流量计算不受地表调蓄水库控制；控制区域内进入地表调蓄水库的地表径流将消耗于蒸散发或用于灌溉等。当地表调蓄水库蓄满后，产生的地表径流将不再受到拦蓄影响。

1.1.2地下人工调蓄模拟

基于流域的数字高程地图、地表水位和地下水位资料，采用“地下调蓄水库”方法刻画产流对地下水超采引起包气带变化的非线性响应以及地下水漏斗对产汇流的间接调蓄效应[6,21]。具体计算过程分为两个步骤：

a.当下渗雨量满足包气带张力水蓄水容量时，多余的雨量将进入自由水蓄水库。由于地下水超采的影响，一部分自由水蓄水库中的水将向下渗漏以补充地下水，抬升地下水位。该渗漏过程采用格林-安普特下渗公式[22-23]进行刻画，即:

(4)

式中：f(t)为栅格自由水蓄水库渗漏能力，mm/h；K为栅格饱和水力传导度，mm/h；ψ为栅格湿润锋处土壤吸力，mm；Δθ为栅格土壤饱和含水率与初始时刻含水率之差；F(t)为栅格累计渗漏量，mm；t为时间。

b.自由水蓄水库的渗漏水量将进入地下调蓄水库，累计下渗水量超过地下调蓄水库阈值S0的部分以地下水出流形式(出流系数Kg)形成地下径流。其中，地下调蓄水库阈值S0的计算方法如下：

(5)

式中：Zr为流域出口河道水位，m；Zi为流域出口所在栅格的相对不透水层高程，m；Zg为栅格多年平均地下水位高程，m；ΔH为栅格与流域出口位置的相对高程差，m；ΔH′为模拟时段初地下水位测站几何中心位置的地下水位高程与其多年平均地下水位高程之差，m；n为栅格土壤孔隙度。

1.2 GXAJ-DAR模型

基于上述产流模拟方法，在栅格新安江模型(Gridded Xin’anjiang model，GXAJ模型)[24]基础上构建适用于受人类活动影响流域降雨径流模拟的GXAJ-DAR模型。不透水面上的直接径流与表层自由水蓄水库蓄满后产生的饱和地面径流汇合，水利工程控制区内栅格产生的地面径流经地表调蓄水库拦蓄，超出蓄水容量部分以地面径流形式汇流至流域出口；表层自由水蓄水库中的部分水量将补充地下调蓄水库，达到阈值后，地下调蓄水库将以某一出流系数出流形成地下径流；模型其余部分与GXAJ模型一致，如图1所示。与基于泰森多边形的半分布式新安江-海河模型(XAJ-HH模型)相比，GXAJ-DAR模型考虑了流域降雨空间分布不均的特点和土壤植被、地形地貌等下垫面因子的空间异质性以及地表水利工程控制区域和地下水埋深的空间分布，可以更精细地刻画人类活动影响下的降雨径流过程。

图1 GXAJ-DAR模型结构示意图

GXAJ-DAR模型在GXAJ模型基础上新增以下参数：稳定下渗率fc、土壤孔隙度n、湿润锋处土壤吸力ψ、水利工程防洪库容Vf、流域出口位置相对不透水层高程Zi。fc、n、ψ可根据土壤性质得到，Vf可根据实际水利工程信息获得，Zi采用人工优选法率定得到。姚成[25]研究发现，在分辨率为1 km的栅格内认为蓄水容量分布均匀可以满足洪水模拟精度要求，因此本研究中的GXAJ-DAR模型将不考虑参数和变量在栅格内部的空间变异性问题。

以下将GXAJ-DAR模型与XAJ-HH模型在研究流域进行应用检验，比较分析两个模型的模拟效果。

1.3 模型率定与评估方法

GXAJ-DAR模型的集总式参数和XAJ-HH模型的参数均采用人工优选法在经验取值范围内进行率定。依据GB/T 22482—2008《水文情报预报规范》，选择径流深合格率、洪峰合格率、确定性系数3个指标来综合评判各模型模拟结果(其中确定性系数大于0.5时假定为合格)，并对本文提出的人类活动影响下的降雨径流模拟方法进行适用性评价。

2 模型验证

2.1 研究流域概况

海河流域为受地下水超采严重影响的典型流域，主要分布于我国渤海以西，黄河以北，位于 35°N～43°N、112°E～120°E之间[26]，流域面积达31.78万km2，长度为450 km，其中60%为山区，40%为平原。多年观测数据显示，海河流域年均降水量仅548 mm，由于气候、地形等多种因素的影响，夏季暴雨集中，冬春雨雪稀少，年降水量分布具有很强的地域性及季节差异性，同时还伴随着连续枯水及连续丰水的周期性变化。

选择河北省保定市清水河流域(38.6°N～39.2°N、114.8°E～115.5°E)作为研究流域，该流域位于海河流域大清河水系，属于半湿润流域。流域出口北辛店水文站多年平均径流量为0.29亿m3(1971—2010年)，多年平均降水量为457.1 mm(1971—2010年)，最大年降水量为 833.7 mm。实测最高洪水位为 19.70 m，最大洪峰流量为1966年的710 m3/s，其次是1989年的 686 m3/s(由于洪水漫滩，实际洪峰流量估计值为 1 500 m3/s)。

清水河流域面积1 650 km2，流域内设有雨量站11个，其中有完整水雨情资料的水文站1个，雨量站6个；设有地下水位测站17个，其中有完整地下水埋深资料的3个(图2)。北辛店水文站上游建有库塘闸坝9座，其中中型水库1座(龙潭水库)、小型水库7座(西显口、北固城水库等)、塘坝1座(水头水库)(图2)。

图2 清水河流域概况

2.2 数据资料

a.水文气象资料。清水河流域的降水量、蒸发量、流量和地下水埋深等水文气象资料由河北省水文局提供，研究所用资料信息包括日尺度和小时尺度的站点降水量和出口流量、小时尺度的地表水位、月尺度的地下水位以及日尺度的蒸发资料，资料年限为1980—2002年。依据2000年前洪峰流量大于25 m3/s、2000年后洪峰流量大于10 m3/s的原则，选取该流域1980—2002年间的12场洪水进行模拟研究(2000年后河道基本断流)，其中8场用于率定，4场用于验证。

b.下垫面资料。研究所需流域下垫面信息包括水利工程控制区域分布和多年平均地下水位相对高程分布(图3)，后者是基于河北省水文局提供的多年地下水位等值线图插值得到的相对于3个地下水位测站几何中心的地下水位高程分布。其余所需下垫面信息与GXAJ模型相同，网格尺度为1 km×1 km。GXAJ-DAR模型的新增土壤参数如表1所示。

(a) 水利工程控制区域分布

表1 GXAJ-DAR模型新增土壤参数值

2.3 模拟精度比较

根据河北省水文局提供的信息，采用清水河流域1980年后10场洪水资料对现有基于降雨径流的经验模型进行参数率定和验证，结果显示洪峰合格率为40%，本研究的模拟统计结果如表2所示。根据GB/T 22482—2008《水文情报预报规范》，当径流深实测值的20%小于3 mm时，即径流深实测值小于15 mm时，许可误差取3 mm。对于清水河流域的12场洪水，除1989072101号洪水外，其余场次洪水的实测径流深均小于15 mm，另外有7场洪水实测径流深小于3 mm，因此两个水文模型对于径流深的模拟误差可基本控制在3 mm以内，径流深模拟合格率均较高，在率定期与验证期均达到100%(表2)；尽管如此，由图4可知GXAJ-DAR模型的径流深模拟误差总体而言小于XAJ-HH模型，经计算，两者的平均相对误差分别为45%和47%。由此可见，GXAJ-DAR模型的径流深模拟效果稍优于XAJ-HH模型。

表2 清水河流域模拟结果统计

(a) 径流深模拟误差

在洪峰模拟方面，XAJ-HH模型在率定期和验证期的模拟合格率分别为37.5%和25.0%，与现有的经验模型相比有所不足，而GXAJ-DAR模型相较于XAJ-HH模型和经验模型在洪峰模拟效果方面有显著提升，率定期和验证期的合格率分别达到62.5%和50.0%(表2)。进一步分析各个场次的洪峰流量模拟误差发现，GXAJ-DAR模型的洪峰流量模拟相对误差明显小于XAJ-HH模型(图4)，经计算，两者的平均相对误差分别为31%和50%。由此可见GXAJ-DAR模型的洪峰流量模拟效果优于XAJ-HH模型和经验模型。

在流量过程模拟方面，验证期GXAJ-DAR模型确定性系数大于0.5的洪水模拟场次多于XAJ-HH模型，两者的确定性系数合格率分别为50.0%和25.0%(表2)，经计算，两者的平均确定性系数分别为0.42和-0.08。图4(c)为两个模型模拟结果确定性系数的分布，其中XAJ-HH模型对第2002062412号洪水的模拟结果与实测过程相差巨大，故该确定性系数异常值点(-5.4)未在图中。总体而言，GXAJ-DAR模型的确定性系数高于XAJ-HH模型，说明GXAJ-DAR模型对流量过程的模拟精度明显高于XAJ-HH模型。

分别挑选率定期和验证期洪峰流量最大和最小的两场洪水做典型模拟流量过程展示，见图5。对于两场大洪水(图5(a)(b))，两个模型在洪峰、洪量和洪水过程三方面的模拟精度均较好，差异不显著。对于两场小洪水(图5(c)(d))，在率定期，GXAJ-DAR模型的洪峰模拟精度高于XAJ-HH模型；在验证期，GXAJ-DAR模型的洪峰和洪量模拟精度均明显高于XAJ-HH模型。由此可见，GXAJ-DAR模型在模拟小洪水方面相比于XAJ-HH模型具有更显著的优势。

(a) 1989072101号洪水(率定期)

综上所述，GXAJ-DAR模型的综合模拟效果优于XAJ-HH模型和现有的经验模型，可达到丙级预报精度。GXAJ-DAR模型分布式的参数和结构能更准确地刻画半湿润地区产汇流过程对流域下垫面特征和降雨、蒸发、土壤湿度等水文气象因子的时空响应，有效提高洪水模拟精度。

2.4 模型参数对比

两个模型的集总式参数率定结果如表3所示。清水河流域的土壤含水量较湿润地区低，降雨过程中包气带往往无法蓄满，且降雨空间分布不均，因此，产流过程对张力水蓄水容量Wm和自由水蓄水容量Sm的空间分布较为敏感[6]。XAJ-HH模型在各个流域分块内采用同样的张力水蓄水容量分布曲线和自由水蓄水容量分布曲线来刻画蓄水容量的空间分布，由于清水河流域地形起伏较大，各个分块的下垫面特性差别明显，因此这种方法的计算结果不能很好地反映流域实际情况；GXAJ-DAR模型则根据地形指数和土壤类型来计算每个栅格内的平均Wm，在此基础上根据植被类型计算栅格内的平均Sm，能较好地反映流域下垫面实际物理特性。

表3 清水河流域模型参数率定值

2.5 地下拦蓄量空间分布

图6和图7分别为洪量最大的1989072101号洪水和洪量最小的1997073108号洪水期间地下拦蓄量的空间分布。1989072101号洪水期间，由于累计降水量较大且前期土壤含水量较高(图6(a))，除流域东南角因初始包气带和自由水蓄水库缺水量相对较大而导致地下拦蓄量相对较小外，流域绝大部分区域的地下拦蓄量均达到较高水平(图6(c))，其中流域的西北部的沙壤土区和部分中部粉壤土区的下渗能力相对较高，使得这部分的地下拦蓄量又明显大于其他区域。

(a) 包气带缺水量

(a) 包气带缺水量

1997073108号洪水期间的地下拦蓄量空间分布呈现出很大的差异。该场降雨雨量小，降雨初期流域大部分区域的张力水和自由水蓄水容量均较低(图7(a)(b))，因此入渗雨量无法满足下垫面包气带和自由水蓄水库缺水量。由于以上原因，流域的地下拦蓄量呈现较低水平，许多区域甚至出现无地下拦蓄现象，地下拦蓄作用主要集中在流域的西北区域和水系附近区域(图7(c))，究其原因，这些区域的初始土壤湿度较大，因此入渗水量足以使其达到田间持水量，多余部分便可进入自由水蓄水库开始产流与渗漏过程。

3 结 论

a.在模拟精度方面，XAJ-HH模型在清水河流域的模拟效果与现有的经验模型相比有所不足；GXAJ-DAR模型的径流深模拟效果稍优于XAJ-HH模型，在洪峰流量和流量过程方面的模拟精度较XAJ-HH模型和经验模型有明显提高，尤其在模拟小洪水方面具有更显著的优势。这说明GXAJ-DAR模型能更精细地刻画降雨径流过程对流域下垫面特征和降雨、蒸发、土壤湿度等水文气象因子的时空响应，有效提高洪水模拟精度。

b.模型结构与参数方面，GXAJ-DAR模型中栅格尺度的地表和地下双调蓄结构能基于水利工程控制区域和地下水埋深的实际空间分布对地面径流和地下径流进行拦蓄；GXAJ-DAR模型对下垫面参数(Wm、Sm)的计算方法能更好地反映流域下垫面特性的实际空间分布，Sm平均值的率定结果更为合理。

c.从应用层面来看，本文提出的GXAJ-DAR模型能有效提高清水河流域人类活动影响下的洪水模拟精度，为地表地下人类活动双重干扰下的洪水预报提供理论和技术支撑，尤其在我国北方流域具有较大的推广应用价值。