华北地区分布式蓄超空间动态组合TOKASIDE-D模型研究

刘玉环，刘志雨,2，李致家，张 珂，刘开磊

(1. 河海大学水文水资源学院，江苏 南京 210098；2. 水利部信息中心(水利部水文水资源监测预报中心)，北京 100053；3. 淮河水利委员会水文局(信息中心)，安徽 蚌埠 233000)

洪水是中国最常见的自然灾害之一，特别是在中国华北地区(半湿润)，陡涨陡落的洪水极易造成严重的生命和财产损失[1-3]。在实际生产中，洪水防治重要的非工程措施之一就是利用水文模型进行洪水预报。研究发现，大部分水文模型在湿润地区应用效果较好，但是在半湿润半干旱地区模拟效果不佳[4-6]。这是由于这类地区的下垫面和降水时空分布不均匀，超渗产流和蓄满产流随时空变化的现象尤为明显[7]。另外，华北地区秋冬季降雨较少，用水量较多，地下水开采量大于补给量，地下水位下降。汛期刚开始时，地下水埋深较深，径流在汇流过程中不断渗漏，用于抬升地下水位[8]。因此，华北地区的洪水预报比湿润地区更具挑战性[9]。

目前，经过国内外学者多年的研究，在通过大量的模型改进后，混合产流问题取得了阶段性成果，从单一产流模型如CASC2D模型[10]、SCS模型[11]、TANK模型[12]等，到混合产流模型包括垂向混合径流模型[13-14]、增加超渗的新安江模型[15-16]等，专门针对海河流域的新安江-海河模型以及河北雨洪模型[17]等，再到灵活框架模型如FUSE[18]、FARM[19]及SCCM[20]等。但是，这些改进后模型的产流模式在计算时都是固定不变的，没有系统地动态识别洪水过程中产流模式随时间变化的情况。

李致家等[21]认为混合产流模拟的主要困难是识别随时空变化的产流机制。对于一个特定流域的某场洪水来讲，流域地形地貌、土壤属性、植被等下垫面特征是基本不变的，但降雨和土壤含水量等因子会随时间而动态变化，并且这类因子的变化与产流模式的空间分布是相互关联的。因此，针对这一现象，本文在TOKASIDE模型(topgraphic kinematic approximation and saturation-infiltration double excess grid-based distributed model)的基础上，通过降雨和土壤含水量之间的判定关系，动态识别蓄满产流和超渗产流网格，并添加渗漏模块，构建基于物理基础的分布式蓄超空间动态组合模型(TOKASIDE-D)，并在我国华北地区的北辛店流域进行验证和分析，探究流域产流区蓄满、超渗产流模式空间动态组合方法的应用效果。

图1 TOKASIDE-D模型结构Fig.1 Structure diagram of TOKASIDE-D model

1 模 型 构 建

1.1 TOKASIDE 模型

TOKASIDE模型是刘志雨等[22]提出的基于地形与运动波的蓄满超渗产流机制的网格化分布式模型。该模型是在TOPKAPI模型的基础上发展起来的，进一步考虑了地下水运动、水库入流与调蓄计算等模块以及超渗产流机制。TOKASIDE模型在空间上采用非线性运动波方程法，将每一个计算单元中的水文过程概化为3个“结构上相似”的非线性水库方程，用以描述水流在土壤层、饱和坡面以及河道中运动波方程的融合。3个方程分别代表上层土壤中的排水、饱和土壤(或不透水土壤层)中的地面径流和河道径流。

1.2 TOKASIDE-D模型

TOKASIDE模型在产流计算模块中虽然考虑了流域超渗产流，但是这一考虑是针对流域内所有网格，即扣除蒸发后的净雨先完成超渗过程，后继续下渗进行蓄满过程，这与实际情况中蓄满产流与超渗产流交替主导的现象不符。在华北地区开展径流模拟探究时，需要对蓄满产流和超渗动态交替以及地下水超采的情形进行考虑。针对这一问题，对TOKASIDE模型进行改进，弥补假设条件中缺失的情况，动态识别蓄满产流和超渗产流网格，以提高模型在华北地区的模拟精度。

在TOKASIDE模型基础上，通过降雨和土壤含水量之间的判定关系，动态识别蓄满产流和超渗产流网格，添加渗漏模块，构建基于物理基础的TOKASIDE-D模型，模型结构见图1(图中p为降雨量，pe为降雨强度，θt为网格土壤含水率，θS为网格田间持水率，f为稳定下渗率，Q为模拟流量)。将流域网格初步划分为“蓄满网格”和“超渗网格”两类。蓄满网格和超渗网格分别采用TOKASIDE的土壤水非线性水库方程和Green-Ampt下渗公式[23]计算，另外针对地下水超采出现的渗漏问题，在模型中添加渗漏模块。

1.2.1 基本假设

基于TOKASIDE模型的假设，根据华北地区下垫面情况，现添加以下假设条件：(a)网格内的降雨、土壤性质、土地利用等特征均匀分布，不存在下渗能力分布曲线。(b)土壤分成2层(饱和层、下渗层)计算。我国北方半干旱地区的特色是土层厚，年降雨量仅在400～800 mm，并且多集中于7—9月，所以其土壤基本不可能达到饱和状态。但其表层土壤还是会达到饱和状态，所以此处将土壤分成2层进行模拟计算[24]。蓄满和超渗发生在饱和层(土壤水模块)；下渗层则为无限下渗土层(渗漏模块)，主要针对地下水超采形成的渗漏问题。(c)蓄超网格划分：河道附近地势平坦，土壤水含量相对较高，河道网格(三级以上)全部为蓄满网格；非河道网格可为蓄满网格或超渗网格。(d)蓄超网格转化：非河道网格在计算过程中，通过判断网格土壤含水量是否达到田间持水量、降雨是否超过下渗率，动态调整网格的蓄满或超渗属性。

1.2.2 土壤水模块修改

在原始TOKASIDE模型中，网格计算分为非河道网格和河道网格。从图1可以看出，在2种网格内，均需要完成土壤水、地表水计算，河道网格另需进行河道汇流演算。因此，产流机制的转换主要发生在非河道网格的土壤水模块。

降雨未发生时，在计算过程中当非河道网格土壤含水率未达到田间持水率时，为超渗网格，达到田间持水率后转为蓄满网格。降雨发生时，降雨强度超过网格下渗率为超渗网格，若网格为蓄满则修改为超渗；否则原网格属性不变。蓄满产流主要包含饱和地面径流和壤中流，而超渗产流计算产生超渗地表径流。

a. 当土壤含水率未达到田间持水率时(θ<θs)，计算网格为超渗网格：

(1)

式中：Vs——网格单元的土壤含水量，m3；X——网格单元尺寸，m；Quo、Qus——从上游区域进入当前网格单元i的地表、地下径流量，m3；hout——时段末土壤水深度，m；Ks——饱和水力导水率，m/s；θ0——初始土壤含水率，m3/m3；Sf——湿润锋面处土吸力，m；I——累积下渗量，m。

如果降雨强度pe≤f，超渗地表径流为hexf=0(hexf为超渗土壤水出流深度，m)，下渗量全部补充土壤含水量：

hout=hin+pet-hexf

(2)

如果pe>f，超渗地表径流为hexf=pe-f，下渗量全部补充土壤含水量至田间持水量，产生壤中流hsoil：

hsoil=hout-hmax

(3)

式中：t——时间，s；hin——时段初土壤水深度，m；hmax——土壤水最大深度，m。

b. 当土壤含水率达到田间持水率时(θ≥θs)，计算网格转换为蓄满网格，其计算流程如下：

(4)

式中：Cs——当地的传导率，m/s；αs——土壤特性决定的参数。

时段初土壤水深度hin=hmax；壤中流hsoil=hout-hmax；饱和地表径流hexf=hout-hsoil-hmax；更新时段末土壤水深度hout=hmax。

c. 当降雨强度pe>f时，蓄满网格修改为超渗网格，则计算模块由蓄满模块转为超渗模块，按照式(1)～(3)进行计算。

1.2.3 渗漏模块

华北地区地下水超采严重，地下水位埋深大，建模时考虑垂向土壤深层的渗透。当土壤水的含水量超过蓄水能力时，穿过饱和层以上的厚土壤层进行垂直向下运移。假定上层土壤层的下渗率是土壤水含量的一个函数，则下渗水量为

(5)

式中：Pr——下渗水量，m/s；Ksv——土壤垂向饱和导水率，m/s；θsat——土壤饱和含水量,m3/m3；αp——取决于土壤类型的指数(沙石:αp≅11；泥土:αp≅25)。

1.2.4 地表水模块

土壤表层饱和后的剩余降水构成地表水模块的输入，地表坡面径流计算与土壤部分的描述相似，公式为

(6)

式中：Vo——地表径流量，m3；ro——由水量平衡方程求解得到的余下水量(或土壤水回流量)，m/s；Co——与地表水流曼宁公式相关的系数；αo——使用曼宁公式得出的指数，取3/5。

1.2.5 河道演算模块

河道径流的运动波近似方程采用运动学方法进行描述，在这种方法中动量方程通过曼宁公式近似写成：

(7)

式中：Vc——单元网格的河道水量，m3；rc——旁侧入流，包括地表径流量和土壤水流量，m2/s；Quc——来自上游单元河道网格的流量，m3/s；S0——河底坡降；nc——河道曼宁糙率，s/m2/3；Cx——湿周，m。

1.2.6 模型参数

TOKASIDE-D模型的主要参数包括计算单元土壤类型对应的土壤厚度L(m)、土壤横向饱和水力导水率Ksh(m/s)、垂向饱和水力导水率Ksv(m/s)，土地利用类型对应的地表曼宁糙率ns，河道分级对应的河道曼宁糙率nc等。Green-Ampt下渗公式的参数：饱和水力导水率Ks(m/s)，湿润锋面处的土壤吸力Sf(m)。在构建TOKASIDE-D模型时，将TOKASIDE模型与Green-Ampt下渗公式中物理含义相同的参数合并，比如：TOKASIDE模型的垂向饱和水力导水率Ksv和Green-Ampt下渗公式的饱和水力导水率Ks，从而减少模型参数的数量。TOKASIDE-D是基于物理基础的水文模型，在进行参数优化率定时，大部分参数可根据土壤和土地利用性质得来，而其余参数使用人工优选法和SCE-UA自动优选法[25]结合率定得到。

2 研究流域与数据资料

2.1 研究流域

北辛店水文站为海河流域大清河水系清水河的控制站，控制流域面积1 650 km2(图2)。该站实测最高洪水位为19.70 m，最大洪峰流量为710 m3/s，多年平均年径流量为0.29亿m3，多年平均年降水量为457.1 mm。该流域属太行山山前倾斜冲积平原区，地势西高东低，缓缓倾斜，土壤以沙壤土为主。流域西部以山地丘陵为主，几乎全为次生林或次生温带灌草丛，东部为平原区，地势平坦开阔，土质疏松，土层深厚，已开垦为农田。

图2 北辛店流域地形及站点分布Fig.2 Topography and station distribution of Beixindian Watershed

近年来，受人类活动影响，海河流域下垫面条件发生了较大变化，流域产汇流规律也发生一定程度的改变。受地下水过量开采的影响，海河流域地下水埋深较深，不易蓄满，产流机制呈现出“先超后蓄”的特点，超渗地面径流产流机制发挥了重要作用。

2.2 数据资料

本文所需下垫面数据为数字高程地图(DEM)、土壤类型和土地利用类型数据。DEM从地理空间数据云网站下载，选择GDEMV2 DEM数字高程数据产品，分辨率为30 m。土壤数据来源于联合国粮农组织(FAO)，1∶100万土壤栅格数据资料。土地利用数据可从美国地质勘探局(United States Geological Survey，USGS)下载，数据为分辨率为1 km的矢量文件。

研究流域水文数据包括降雨及流量数据(时段为1 h)，由河北省水文水资源勘测局提供。气温数据来自中国国家气象数据中心。流量数据来自北辛店水文站，降雨资料来自黄龙寺、顺平等6个雨量站(图2)。选用北辛店流域1980—2002年的12场次洪水数据，其中1980—1996年的9场为率定期，1997—2002年的3场为验证期。

3 TOKASIDE-D模型的应用检验

将TOKASIDE模型、Green-Ampt模型、TOKASIDE-D模型用于北辛店流域进行应用检验。依据传统的水文模拟与预报精度评定准则，并参考《水文情报预报规范》[26]的规定，结合半干旱地区洪水特征，选择以下评价指标：径流深误差合格率，该误差以实测径流深的20%作为许可误差判定预报径流深是否合格；洪峰误差合格率，以实测洪峰的20%作为许可误差；峰现时差合格率，以峰现时间小于3 h为许可误差；平均确定性系数，评价洪水实测过程与预报过程之间的拟合程度。北辛店流域3种模型参数率定的结果见表1，模拟结果见表2。

表1 北辛店流域各模型的优化参数

表2 北辛店流域3个模型的模拟结果统计指标

由表2可得：径流深误差合格率方面，率定期TOKASIDE-D模型最高(56%)，而Green-Ampt模型最低(33%)，验证期3个模型均只有1场合格，合格率33%；洪峰误差合格率方面，率定期TOKASIDE-D模型也是最高(78%)，TOKASIDE模型最低(44%)，验证期TOKASIDE与Green-Ampt模型合格率均为33%，而TOKASIDE-D模型高于两者。峰现时差方面，TOKASIDE模型和TOKASIDE-D模型的率定期和验证期合格率一致；但均低于Green-Ampt模型的合格率，说明TOKASIDE-D模型在模拟洪峰时间方面与原TOKASIDE模型差异不大，土壤入渗后补充土壤水的过程延长了饱和地面径流的出现时间，进而使得模拟洪峰时间出现滞后。TOKASIDE-D确定性系数比原模型有提高，说明改进模型对洪水形状的模拟与观测值更接近。

此外，选取2个典型的洪水事件1996080420和1980081700进行详细分析，其洪水实测过程线与模拟过程线如图3所示。首先，从2个洪水事件的实测值可以看出，2个洪水的洪峰值差不多，前者洪水起涨快，退水慢；后者的洪水过程陡涨陡落，呈尖瘦形态。从实测过程线可以大致推断出1996080420洪水的降水量大于1980081700，前者降水更多的是入渗到土壤中，而后者洪水明显由超渗产流主导，降水入渗量少，基本上发生了超渗产流。

图3 北辛店流域部分洪水事件实测过程线与模拟过程线对比Fig.3 Comparison of observed and simulated hydrographs of partial flood events in Beixindian Watershed

对比3个模型的模拟过程线，在1996080420洪水事件中(图3(a))，TOKASIDE-D模型的模拟洪水过程线在起涨阶段以及洪峰值效果最好，模拟洪水的整体形态与实测洪水最接近。TOKASIDE模型和Green-Ampt模型的模拟洪水起涨速度都很快，但模拟洪峰值都比实测洪峰大得多，且退水速度快，与实测退水过程偏差较大，说明这2个模型在模拟时将大部分净雨拦在了地表，转为超渗地表径流，入渗补充土壤含水量的部分偏少；而TOKASIDE-D模型在涨水阶段的土壤蓄水量在退水时被“释放”出来。但相对实测的退水流量还是偏小，这主要是下渗的水量在渗漏模块的作用下继续下渗补充地下水，从而出现部分水量不足的情况。在1980081700洪水事件中(图3(b))，3个模型的模拟洪水过程线都呈现陡涨陡落的形态，与实际洪水过程线拟合较好。这是因为1980081700的洪水过程陡涨陡落，呈尖瘦形态，是比较典型的超渗产流洪水过程，观察降雨过程，暴雨中心在流域中下游，因此得以快速汇流至流域出口，并没有太多的下渗损失。另外，洪水发生在1980年，人类活动和水利工程相对较少，比较接近天然的降雨径流过程，因此3种模型在该场次洪水模拟的效果相对较好。

与原TOKASIDE模型对所有网格设置超渗“门槛”不同的是，TOKASIDE-D模型通过对每时刻的土壤含水量和降雨等动态因子进行“监控”，动态调整网格的产流机制。以1996080420和1980081700这2个典型的洪水事件为例，展示洪水计算过程中产流机制的动态变化过程，如图4和图5所示。从图4可以看出，1996080420洪水从8月4日20:00开始降水，一直到5日的2:00，降雨8 h，流域基本还是超渗状态，在3:00流域出口附近网格基本蓄满，在随后的3 h内，网格蓄满数量大量增加，直到6日5:00达到最大，此时模拟洪水也达到了峰值，随后，网格土壤含水量开始下降，洪水过程处于退水阶段，蓄满网格数量逐渐减少。在这里也可以解释，TOKASIDE-D模型在1996080420洪水的模拟中，网格蓄水程度高，退水时网格中有大量的壤中流和地下径流产生，从而在退水阶段模拟优于其他2个模型。从1996080420洪水的实测过程可知其由超渗主导，这一点从图5也可以看出，整个洪水过程蓄满网格数量确实远远少于1996080420洪水，超渗网格占主导地位，进而模拟的洪水过程线呈现尖瘦形态。

图4 1996080420洪水的蓄满网格与超渗网格动态分布Fig.4 Dynamic distribution of saturation-excess and infiltration-excess grids for flood event 1996080420

图5 1980081700洪水的蓄满网格与超渗网格动态分布Fig.5 Dynamic distribution of saturation-excess and infiltration-excess grids for flood event 1980081700

综上可知，Green-Ampt模型在北辛店模拟效果相比其他2个模型较差，主要因为产流机制不符合当地的特征。北辛店流域降雨初期土壤水含量低，容易发生超渗产流，而后土壤入渗能力强，土壤蓄水量大，转为蓄满产流主导。由于Green-Ampt模型为超渗产流模型，净雨大多发生产流，转为超渗地表径流，导致模拟洪水过程远大于实测过程。TOKASIDE模型虽然表现尚可，但是其模拟洪峰与实测洪峰相差较大；而TOKASIDE-D模型结合了蓄满产流和超渗产流2种产流机制，并且根据土壤含水量以及降雨的时空分布，动态调整产流机制，能更加准确地把握北辛店流域的产流机制的过程，适应性更强。该地区的渗漏严重，虽然增加了渗漏模块，但部分场次洪水模拟结果依旧较差，在后续研究中，进一步探究该流域的存在的问题以及独有特征，寻找合适的解决方法。

4 结 语

在基于物理基础的分布式模型TOKASIDE模型的基础上，通过降雨和土壤含水量之间的判定关系，动态识别蓄满产流和超渗产流网格，并添加渗漏模块，构建了基于物理基础的分布式蓄超空间动态组合TOKASIDE-D模型，并在北辛店流域进行试验应用。研究结果表明，与原始模型(TOKASIDE模型、Green-Ampt模型)相比，TOKASIDE-D蓄超空间动态组合模型在北辛店流域为代表的华北地区洪水模拟精度有较为明显的提升。

TOKASIDE-D模型能够很好地结合蓄满模型和超渗模型的优点，并且根据土壤含水量与降雨等因子的变化，动态调整产流机制，并针对地下水超采考虑地下水的渗漏。随着模型结构与原理的改进完善，TOKASIDE-D模型能够准确捕捉降雨、土壤含水量等动态因子对产流模式的控制效应，在我国北方地区的中小流域，有着很大潜力的研究价值以及较高的应用价值。