基于暴雨壤中流机制的山区水文模型研究与应用

张艳军，宋圆馨，董文逊，邱安妮，罗 兰，黄 尔

(1.武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉 430072；2.四川大学 水力学与山区河流开发保护国家重点实验室，四川 成都 610065)

中国山丘区面积广大，山洪灾害频发，给人民的生命财产安全造成了巨大威胁。山区水文模型能有效提高山洪预警效果和山洪灾害防御能力，对推动山洪灾害防御体系从“有”到“好”转变[1]具有十分重要的意义，但是建立精准的山区水文模型十分困难。国家科技部于2017年发布了“山洪灾害监测预警关键技术与集成示范”重点专项申报指南，指南中仅要求“山洪洪峰流量预报精度要由40%提高到50%”[2]，远远低于《水文情报预报规范（GB/T 22482—2008）》[3]对常规水文预报“以实测洪峰流量的20%作为许可误差”的要求，由此可见，对山区小流域进行洪水模拟和预报的难度之大。

山区小流域水文预报的主要难点在于水文资料缺乏，产流机理复杂，降雨分布空间异质性强等。其中，复杂的产流机理是主要因素之一。目前常用于山区水文预报的模型主要基于蓄满、超渗或者两者混合的产流机制。例如，美国国家天气局萨克拉门托预报中心研制的SAC模型[4]的产流机制为蓄满–超渗产流兼容，但由于模型参数复杂，难以调试，限制了其广泛应用[5]。在中国，李琼等[6]将新安江模型、SAC模型及TOPMODEL模型应用于湖北省宜昌市的3个山区小流域，发现3个模型的预报精度都不高，认为简单的蓄满产流和超渗产流模式不足以描述降雨与下垫面的空间不均匀对产汇流的影响。吴金津等[7]将新安江模型、TOPMODEL模型和TVGM模型应用于白沙河山区小流域，发现TVGM模型的效果最好，认为可能是白沙河流域产汇流机制较为复杂，并非简单的蓄满产流，而TVGM模型无明确的产流机制，为系统性的黑箱子模型，相对更贴合山区小流域的实际情况，但同时也掣肘了山洪预报精度的进一步提高。张相芝等[8]在垂向混合产流模型的基础上，在4个山区小流域通过GIS和遥感技术识别小流域基础地貌信息，针对不同的地貌水文单元选择不同产流计算方法，较好地模拟了降水径流的非线性过程，但受制于遥感观测的精度，不能准确归类流域产流特性，从而难以推广。针对山区小流域坡面比降大，汇流历时短、水文响应快等产汇流特点，中国水科院于2014年开发了基于混合产流机制的中国山洪水文模型（CNFF–HM）[9]；之后，翟晓燕等[10]将其应用于安徽省中小流域，Wang等[11]将其应用于福建省山区流域，取得了较高的模拟精度；刘昌军等[12]针对山丘区中小流域地形地貌多样、产流机制时空变化复杂等问题，建立了在平面、垂向、时段上时空变源混合产流模型，模拟精度也显著提升；但上述两个模型对山区小流域暴雨洪水过程的成因机制缺少探究。综上所述，探明小流域下垫面特征及产汇流特性是研制山洪水文模型急需解决的关键技术问题之一[9]。

但是，山区小流域的产流机制可能既非典型的蓄满产流，也非典型的超渗产流[13]，而是暴雨壤中流理论。吴金津[14]和王俊勃[15]曾对官山河流域的水文资料进行分析，发现官山河流域内会出现降雨量大而径流量小、降雨量小而径流量大的异常情况，进而对官山河黄沟的径流成分进行分析，发现该区域的暴雨壤中流比重远大于地面径流，符合暴雨壤中流机制。暴雨壤中流是由Mirus等[16]在比较上百个不同流域的研究成果后，总结的湿润山区小流域的产流模式。在这种产流模式中，径流主要来源于暴雨壤中流，传统的饱和地表径流和超渗地表径流占河道径流的比重都很小。基岩地形直接控制水流的运动方向、存储和泄流过程[17]，当雨强较小或者土壤前期含水量较小时，土壤深度也会影响到暴雨壤中流的出流过程[18]。近年来，国内外许多学者对这种机制进行了探究。Abdel–Fattah等[19]发现山洪暴发的水文过程与流域的尺度和基岩地形息息相关。Camporese等[20]提出土壤–基岩界面的形状和底部基岩的厚度对于土壤的贮水泄流能力至关重要。Wang等[21]认为地形是影响山坡降雨径流过程的主要因素之一，且山坡纵剖面的曲率相较于平面形状因子对径流过程的影响更大。但由于暴雨壤中流的机理复杂，数据难得，集成该产流机制的水文模型并不丰富。Zhang等[22]基于Dupuit–Forchheimer模型、扩散波方程提出了FSF（fully subordinated flow）模型来模拟暴雨壤中流，但只能在整个山坡运用，不能用于局部山坡。Paniconi等[23]运用1维的HSB（hillslope-storage Boussinesq）模型对山坡壤中流进行模拟，指出了HSB模型可以用简单的山坡蓄水概念模拟复杂的壤中流出流过程,但由于其非线性形式，求解过程往往很困难。Babaali等[24]通过拓展传统的地表Nash单位线模型，开发了能够模拟壤中流和地下流的SNIUH模型，但该模型对于下渗模式的选择很敏感，若选择不当则难以准确预测。目前，对暴雨壤中流的产流过程仍有待进一步研究，以建立相应的流域水文模型。

本文选取官山河流域为典型流域，根据暴雨壤中流的贮水泄流机制提出了3阶段的贮水泄流模型；基于TOPMODEL模型的框架，构建了基于暴雨壤中流机制的山区水文模型（subsurface storm flow-based mountain hydrological model，SSFM）；然后，将这个模型应用于官山河流域场次洪水模拟，并将模拟结果与经典水文模型的模拟结果进行对比分析，得到了满意的模拟结果，提高了山洪预报的精度。

1 研究区概况

官山河流域位于湖北省丹江口市西南部，地处东经110°42′30″～111°34′59″，北纬32°13′16″～32°58′20″。流域总面积465 km2，孤山水文站控制流域面积为322 km2，河道平均坡降0.57%，多年平均流量7.78 m3/s。流域地形以中小起伏山地、丘陵为主，海拔为240～1 606 m，平均高程690 m。官山河流域地势中间低、边缘高，最低点为流域出水口，而且出口处的弯曲、窄河段、上下游卡口区条件，不利于快速泄洪[25]，因此该流域也是山洪灾害、泥石流灾害和滑坡灾害多发的地区，历史上曾暴发著名的1935年的“35.7”暴雨和1975年的“75.8”暴雨。

流域内有大马站、袁家河站和西河站3个气象站，孤山站1个水文站，如图1所示。各站控制面积及权重见表1。

图1 官山河流域地理位置及测站位置分布Fig. 1 Geographical location of Guanshan River basin and distribution of stations

表1 官山河流域各站控制面积及权重Tab. 1 Control area and weight of each station in Guanshan River basin

2 模型研究与改进

2.1 暴雨壤中流

暴雨壤中流（subsurface storm flow, SSF）也称之为壤中水径流[26]。当山区有降水发生时，由于山区土壤下渗率较大，孔隙较多，降水落地后会快速下渗，并积蓄在基岩界面或不透水界面上，使得附近的土壤湿润饱和，进而主要以壤中流的形式产流，形成以暴雨壤中流为主的产流机制。

具体来说，暴雨壤中流的产流分为3个阶段：1）基岩贮水阶段。基岩或暂态不透水层上的凹洼区尚未填满时，降水需先填洼，不产流，此时可能会出现“雨大水小”现象，即雨量站监测到较大的降雨，但是河道水文站却只监测到较小的洪水。2）快速泄流阶段。基岩或暂态不透水层上的凹洼区已填满时，暂态饱和区扩大，流域产流，将土壤–基岩区的事件前水排出，此时可能会出现“雨小水大”现象，即雨量站监测到较小的降雨，但是河道水文站却监测到了较大的洪水。3）地表泄流阶段。暂态饱和区进一步扩展到地面，此时出现地表径流。对于壤中流的产流机制的模拟，本文采用的是贮水泄流模型，如图2所示。其中，图2（a）描述了3阶段的贮水泄流过程，图2（b）描述了坡面蓄水过程。

图2（a）中：横轴H为土壤–基岩界面以上的饱和地下水位，m；纵轴q为暴雨壤中流的单宽排水率，m2/s。

图2（b）中，假设基岩的斜坡倾角为α，坡长为L，土壤厚度为D（基岩深），有[27]：

图2 贮水泄流机制Fig. 2 Storage and discharge mechanism

式中：V为单宽饱和区体积，m2；q为单宽排水率，m2/s；i为非饱和区向饱和区的输水速度，m/s。

1）基岩贮水阶段。假设在土壤单元上，饱和土壤水水面存在恒定的坡度，且该水力坡度等于基岩的坡度。基岩或暂态不透水层上的凹洼区尚未填满时，降水需先填洼，不产流，假设土壤–基岩界面上的凹洼需要水位至H0才可填满，继而产流，即贮水亏缺水位为H0。考虑到官山河流域是山区小流域，假设降雨空间分布基本一致，初期降水多用于填洼、截留等，因此计算洪水过程线起涨点前的累计雨量，可视为补充贮水亏缺。当H≤H0时，不产流，有：

式中， ωd为土壤有效孔隙度。

2）快速泄流阶段。当土壤–基岩界面上的贮水亏缺补偿完毕，即凹洼填满后，暂态饱和区开始产流。此时当D/cos α≥ H>H0时，有：

式中，Ks为排水系数，m/s。

当饱和区到达地表，即水位到达D/cos α时，暴雨壤中流的单宽排水率为qD，如图2（a）所示。

3）地表泄流阶段。当饱和地下水面抬升到土壤表面后，产生地表径流，当H>D时，有[27]：

式中，Ls为饱和坡长，m。

当饱和地下水面继续抬升至整个倾斜坡面全部饱和时，地下水位为Hs，暴雨壤中流的单宽排水率为qs，如图2（a）所示。

2.2 基于暴雨壤中流的山区水文模型

本文提出的基于暴雨壤中流机制的山区水文模型（subsurface storm flow-based mountain hydrological model，SSFM）借用TOPMODEL（topography-based hydrological model）的基本框架，将流域内任一点的土壤划分为3个相异的含水层：第1个是植被根系层（Srz）；第2个是土壤非饱和层（Suz）；第3个是饱和地下水区，用饱和地下水水面到流域土壤表面的距离Zi表示。基于暴雨壤中流的产流模型框架如图3所示，其中，图3（a）为TOPMODEL模型的基本框架，图3（b）为借用TOPMODEL模型基本框架构建的基于暴雨壤中流机制的产流模型结构。

图3 基于暴雨壤中流的产流模型框架Fig. 3 Schematic diagram of a runoff generation model framework based on subsurface storm flow

如果将研究流域划分成若干单元网格，其水分运动规律如下：降雨先按一定速率下渗进入植被根系层，补给该层土壤缺水量，该层水分会有部分蒸发，满足其贮藏容量和蒸发消耗后的水量能够进一步下渗进入非饱和层土壤。在非饱和层中，水分按重力排水速率Qv下渗进入饱和地下水层，接着经由侧向流动形成基流Qb。随着饱和地下水位逐渐抬高，在流域的下游河道等处便有水分冒出，形成饱和坡面流Qs。所以在该模型中，流域总径流Q是基流与饱和坡面流之和，即：

与TOPMODEL模型相比，本模型引入暴雨壤中流概念，改进之处如下：

1）增加植被冠层截流部分；

2）在饱和地下水区，将从非饱和区下渗来的水量分成暴雨壤中流和基流两部分；

3）使用基岩（或暂态不透水层）埋深代替原模型中的Szm（非饱和区最大蓄水深度），进而提高模型精度；

4）土壤–基岩界面上暂态饱和区的出流由贮水泄流模型的公式计算。

2.2.1 产流

1）冠层截留

假定，降雨先落在冠层上，当截留量达到冠层的最大截留能力时，降雨穿过冠层，剩下的降雨才继续落在地表。同时，冠层蒸发只发生在潮湿的叶面上，干燥的叶面不蒸发。该过程使用的具体公式如下：

式中：St+Δt为t+Δt时刻的植截留量，mm；St为t时刻的植截留量，mm；P0为降雨量，mm；F为叶面积比率；Imax为冠层的截留能力，mm；Ei为每个网格单元的实际蒸发量，mm。

2）植被根系区下渗蒸发模式

在模型计算中，假设植被根系层的缺水量即为Horton超渗产流中土壤表层缺水量。假设截留后净雨以Horton模式下渗，下渗后剩余的净雨形成初步坡面径流。下渗的雨量先对植被根系区进行填充、蒸发，填充蒸发后如果有剩余，则全部流入下一层非饱和区。该过程使用的具体公式如下：

超渗产流的计算公式：

式中：ft为下渗率，m/h；fc为稳定下渗率，m/h；f0为土壤含水量为零时的下渗率，m/h；Srz为植被根系缺水量，m；Srmax为植被根系最大缺水量，m。

降雨在坡面上的产流量为：

式中，qr为下渗量，m。

更新该层土壤状态，用以提供下一时刻计算的初始状态。公式如下：

式中：td为时间参数，h/m；Suz为非饱和区土壤含水量，m。这样，各点实际下渗水量为：

4）暂态饱和区的暴雨壤中流

通过大量的实验表明[28–29]，在上一层土壤中，下渗的水除了对土壤进行填充外，另一部分水通过土壤中的大孔隙，以快速的壤中流的形式侧向排出。基流为：

2.2.2 汇流

本文中，坡面汇流过程采用Nash瞬时单位线法，表达式如下：

式中，k为蓄泄系数，n为线性水库的个数，Γ为伽马函数，e为欧拉数。暴雨壤中流的汇流过程采用基于DEM的等流时法，即求出各单元上产生的净雨及此单元流到出口断面的时间，就可以得到流域出口断面的出流过程[30]，其计算公式如下：

式中，ti为流经网格单元所用的时间，Li为流程长度，Sv为流速。

3 模拟结果分析

本文利用TVGM模型、新安江模型、TOPMODEL模型及SSFM模型对官山河流域1973—1983年及2010—2017年12场洪水进行模拟，其中，前8场为率定期，后4场为验证期，选取遗传算法率定参数，模拟结果见表2。

表2 官山河流域洪水模拟结果Tab. 2 Flood simulation results of the Guanshan River Basin

由表2可知：从洪峰误差角度来看，TVGM模型、新安江模型、TOPMODEL模型、SSFM模型总体的洪峰平均误差分别为–29.61%、–51.74%、–29.08%、–24.82%。因此，从洪峰误差的角度来看，SSFM模型的模拟效果相对较好，TVGM和TOPMODEL模型模拟的效果一般，新安江模型的模拟效果相对较差。

从峰现时间误差来看，新安江模型总体的峰现时间误差合格率达到75%；SSFM模型总体的峰现时间误差合格率达到66.7%；TOPMODEL模型和TVGM模型总体的峰现时间误差合格率都较低，不足50%。新安江模型和SSFM模型在峰现时间误差上的模拟效果都尚可，而TOPMODEL和TVGM模型的峰现时间误差模拟相对较差。

从表2中的径流深误差来看，TVGM模型、新安江模型、TOPMODEL模型、SSFM模型总体的径流深平均误差分别为–30.83%、26.87%、–18.43%、–9.67%。SSFM模型在径流深方面的模拟效果要明显优于其他3个模型。

TVGM模型、新安江模型、TOPMODEL模型及SSFM模型在官山河流域进行山洪模拟的典型洪水过程线的结果参见图4～8。由图4～8可见：各场次洪水的雨洪响应明显。大部分洪水，例如图4、5、7和8洪号“19730410”“19730430”“19770718”和“19800802”模拟过程中，实测流量过程线表现为单峰洪水，TVGM模型、新安江模型、TOPMODEL模型及SSFM模型模拟的流量过程线与实测洪水过程线趋势一致，但模拟流量普遍偏小。图6中洪号为“19740521”的洪水为多峰洪水，SSFM和TOPMODEL模型模拟的多峰洪水与实际洪水涨落过程相似，模拟的洪峰流量偏大，但在许可误差范围之内。

图4 洪号为“19730410”的典型场次洪水的模拟结果Fig. 4 Simulation results of typical flood “19730410”

图5 洪号为“19730430”的典型场次洪水的模拟结果Fig. 5 Simulation results of typical flood “19730430”

图6 洪号为“19740521”的典型场次洪水的模拟结果Fig. 6 Simulation results of typical flood “19740521”

综合分析TVGM模型、新安江模型、TOPMODEL模型及SSFM模型在官山河流域的洪水模拟效果，认为新安江模型的模拟效果最差，SSFM模型的模拟效果最好。究其原因，因为新安江模型的产流机制是蓄满产流机制，而SSFM模型的产流机制引入了暴雨壤中流机制，在诸如官山河的山区小流域，地形条件复杂多变，后者更能反映地形地貌的空间异质性。

图7 洪号为“19770718”的典型场次洪水的模拟结果Fig. 7 Simulation results of typical flood “19770718”

图8 洪号为“20170902”的典型场次洪水的模拟结果Fig. 8 Simulation results of typical flood “20170902”

4 结 论

1）本文提出了3阶段的贮水泄流模型，构建了基于暴雨壤中流机制的山区水文模型SSFM。以官山河流域为典型，选取12场实测洪水对其进行模拟，并将SSFM的模拟结果与经典水文模型进行对比，发现SSFM模型取得了更好的模拟效果。此外，SSFM模型中，通过确定湿润山区小流域山坡上的土壤–基岩界面的贮水亏缺补偿过程，推导的3阶段贮水泄流模型可以解释山区小流域“雨大水小，雨小水大”的反常现象，这说明可以在暴雨壤中流的“贮水亏缺”和“快速泄流”方面开展深入研究。

2）从洪水3要素来看，TVGM模型、新安江模型、TOPMODEL模型、SSFM模型4种模型在官山河流域的模拟效果最好的是SSFM模型，其次为TOPMODEL模型。SSFM模型在官山河流域的山洪模拟中适用性更好，平均洪峰误差仅为–24.82%，平均径流深误差仅为–9.67%，可以较好地模拟官山河流域的洪水过程，可用于该流域的山洪预警预报。

3）在以官山河流域为例的山区小流域，由于TOPMODEL模型计算了流域的地形指数，在山洪模拟中可以体现山区小流域的空间异质性，所以模拟效果相对较好。而SSFM的模拟效果比TOPMODEL模型进一步提高，说明模型改进是成功的，也反向证实了官山河流域的产流机制符合暴雨壤中流机制。

4）中国山洪预报预警的精度有待进一步提高。尽管SSFM模型相较于TOPMODEL模型在洪峰、径流深和峰现时间误差方面的模拟精度都有所提升，达到了国家科技部发布的指南中“山洪洪峰流量预报精度要由40%提高到50%”[2]的要求，但是模拟效果仍有很大的提升空间。若要进一步提高山洪预报预警的精度，需要继续深入研究山区小流域复杂的产汇流机制，以进一步改进本文的SSFM模型。