栅格型新安江模型的参数估计及应用

姚 成,纪益秋,李致家,刘开磊

(1.河海大学水文水资源学院,江苏南京 210098;2.南通市港闸区幸福街道办事处,江苏南通 226012)

基于栅格的分布式水文模型是目前水文研究的热点之一.1969年,Freeze等[1]发表的“一个具有物理基础数值模拟的水文响应模型的蓝图(FH69蓝图)”被认为是分布式水文模型研究的开始.过去受制于技术的限制,许多流域信息、特征不能及时准确地获得,地理信息无法及时准确地获得和使用,使得分布式水文模型发展缓慢.而数字高程模型(DEM)的出现为分布式水文模型的构建创造了条件,分布式水文模型取得了飞速发展.比较具有代表性的分布式和半分布式流域水文模型有SHE[2],TOPKAPI[3]和TOPMODEL[4]模型等.随着信息、遥感、计算机等技术的发展,基于DEM的分布式水文模型以其充分考虑水文要素和各种参数空间变化的特点正成为水文模型的发展趋势,而新安江模型[5]作为一个半分布式的概念性模型,自研制以来,就得到了广大学者的认同,在中国的洪水预报中也得到了广泛成功的应用,因此研究栅格型新安江模型具有一定的现实意义[6-7].本研究利用DEM,提取了水文模拟与参数估计所需要的流域地貌特征,并根据新安江模型与扩散波模型理论,构建了栅格型新安江(Grid-Xin'anjiang)模型.同时,本文结合土壤类型、植被覆盖及利用等信息,对空间变化模型参数的估计方法进行了一定的研究与验证,并将Grid-Xin'anjiang模型应用于安徽省屯溪流域洪水模拟计算.

1 Grid-Xin'anjiang模型原理

Grid-Xin'anjiang模型以流域内每个DEM栅格作为计算单元,并假设在栅格单元内的降雨和地貌特征、土壤类型以及植被覆盖等下垫面条件空间分布均匀,模型只考虑各个要素在不同栅格之间的变异性.Grid-Xin'anjiang模型先计算出每个栅格单元的植被冠层截留量、河道降水量和蒸散发量,然后再计算出栅格单元的产流量,并采用自由水蓄水库结构对其进行水源划分,即划分为地表径流、壤中流以及地下径流3种水源,最后再根据栅格间的汇流演算次序,依次将各种水源演算至流域出口.在进行次洪模型计算时,采用一维扩散波模型,且假设在原来的坡地栅格内也存在一个“虚拟河道”,每个栅格单元的壤中流和地下径流直接流入河道,采用河道汇流的计算方法,则栅格间的扩散波汇流演算只包括坡面汇流和河道汇流2种方式.考虑到日径流模拟时对汇流演算精度要求不高且时间步长较大,为了增加模型的运行效率,Grid-Xin'anjiang模型采用简便的线性水库法和滞后演算法进行栅格间汇流计算.

1.1 植被冠层截留

植被冠层截留是指降雨在植被冠层表面的吸着力、承托力及水分重力、表面张力等作用下储存于其表面的现象[8].在一次降雨过程中,植被冠层对降雨的累积截留量[9]可表示为

式中:Icum——植被冠层的累积截留量;flc——植被覆盖率;Pcum——累积降雨量;Cvd——植被密度的校正因子;Scmax——植被冠层的截留能力,即植被冠层的最大截留量.

1.2 河道降水

假设在每一个河道栅格内,降雨分布均匀,河道形状不发生改变,则可以根据栅格单元内河道部分所占的面积比例进行河道降水的计算.河道降水Ich可表示为

式中:Lch——河道长度;Wch——河道断面最大过水面积所对应的水面宽;Agr——栅格单元的面积;P——时段降雨量.其中,LchWch表示河道所占面积.对于流域内的坡地栅格,即使有“虚拟河道”存在,也不考虑河道降水的影响,即在坡地栅格上Ich=0.

1.3 蒸散发

将每个栅格单元内的土壤分为上层、下层和深层3层,每一层对应的张力水蓄水容量分别为 W UM,W LM和W DM.在栅格单元实际蒸散发计算时,冠层截留量按蒸散发能力蒸发,当截留水量小于蒸散发能力时,则采用3层蒸散发模型[10].

1.4 单元产流及分水源

土壤蓄满表示的是土壤含水量达到田间持水量,而不是饱和含水量.新安江模型引进张力水蓄水容量分布曲线来考虑土壤含水量面上分布不均的问题,而在Grid-Xin'anjiang模型中,由流域地貌特征以及土壤、植被等下垫面条件来确定任意一个栅格单元的张力水蓄水容量WM,而在栅格单元内暂不考虑张力水含水量分布不均的问题.将计算时段内栅格单元的实测降雨先扣除相应时段的蒸散发、植被冠层截留、河道降水后,再考虑上游入流是否补足当前单元的土壤含水量,即可得到实际用于产流计算的时段雨量Pe.当Pe≤0或Pe+W0≤WM时,R=0;当 Pe+W0＞WM时,R=Pe+W0-WM.其中:R为时段产流量;W0为栅格单元实际的张力水含量.

在Grid-Xin'anjiang模型中,任意栅格单元内的产流量R均被划分为3种水源:地面径流Rs、壤中流Ri和地下径流Rg.所需参数包括表层土自由水容量SM、表层自由水含量对壤中流的出流系数Ki、表层自由水含量对地下水的出流系数Kg.

1.5 汇流演算

Grid-Xin'anjiang模型在进行短时段栅格间扩散波汇流演算时,假设任意栅格单元都由坡地和河道组成,即原来的坡地栅格上也存在一个“虚拟河道”,地下径流与壤中流都直接汇入河道或“虚拟河道”中,因此栅格间的汇流就由坡面汇流及河道汇流组成,均采用扩散波模型.在模型进行产汇流计算时,考虑了栅格间的水量交换问题以及河道排水网络的影响.即:如果当前计算单元的土壤含水量处于未蓄满状态,则上游栅格的入流量首先补充当前栅格的土壤含水量,直至其蓄满为止;如果当前栅格有河道存在,则地表径流的出流量将先按一定的比例汇入河道,然后再汇至下游栅格.

Grid-Xin'anjiang模型的坡面水流以及河道水流运动均利用一维扩散波方程组来描述,在进行栅格间汇流演算时,坡面水流运动方程需在每个栅格单元上进行离散,其中的连续性方程为

式中:hs——坡面水流的水深;t——时间;Agc——栅格单元的面积;Qs——栅格单元的地表径流流量;Qsout——栅格单元的地表径流出流量;Qsup——上游栅格入流量.

本文采用基于两步MacCormack算法[11-12]的二阶显式有限差分格式进行坡面与河道水流扩散波方程组的求解.

2 研究流域概况

选用湿润地区的屯溪流域作为研究流域.该流域位于安徽省境内皖南山区,流域面积为2692.7 km2.屯溪流域地势西高东低,最大、最小以及平均海拔高程分别为1398m,116m,380m,相对高差较大.该流域雨量充沛,多年平均降雨量约为1800mm,降水在年内年际分配极不均匀,汛期内的降雨量一般占年总雨量的60%以上.屯溪流域内植被良好,主要包括常绿针叶林、落叶阔叶林、混合林、森林地、林地草原、牧草地与作物地,土壤类型主要为黏壤土.本研究采用USGS提供的精度为1km(30″)的DEM高程数据提取了屯溪流域的地形地貌特征.屯溪流域水系及站点分布见图1.

图1 屯溪流域水系及站点分布Fig.1 River network and distribution of gauging stations in Tunxi Basin

3 Grid-Xin'anjiang模型参数估计

3.1 蒸散发参数

在湿润流域,平均的上、下层张力水容量一般可以取20mm与60mm,而整个包气带张力水容量常以120mm作为估值[5,10].以此为基础,对与屯溪流域W UM与W LM的估计分别取 W UM≈0.167W M;W LM≈0.5WM.参数C与栅格单元的植被覆盖率有关,在植被密集地区可取0.18,因此可令C=0.18flc,flc可通过相关公式求得.本研究暂未考虑蒸散发折算系数K在不同栅格单元的高程修正问题,而是认为它在流域内空间分布均匀,因此该参数主要与测量水面蒸发所用的蒸发器有关.对于国内普遍采用的E-601蒸发器而言,可取K≈1.

3.2 产流及分水源参数

单元产流及分水源计算中WM与SM[13]可表示为

式中:WM——栅格单元张力水蓄水容量;SM——栅格单元自由水蓄水容量;θs——饱和含水量;θfc——田间持水量;θwp——凋萎含水量;La——包气带厚度;Lh——腐殖质土层厚度.

式(4)与式(5)中的θs,θfc,θwp均可以根据栅格单元的土壤类型通过查土壤参数统计表获取,因此只要知道每个栅格单元的La与Lh即可获得WM与SM在流域的空间分布.在自然界中,影响包气带厚度的因素较多,很难进行直接推求.L a与 L h可通过与地形指数及土壤类型对应的土壤水分常数进行估算,可假定地形指数大的地方包气带较薄而地形指数小的地方包气带较厚.这与实际情况基本相符.一般而言,在湿润地区地形指数大的地方大多位于河道附近,而这些区域的地下水埋深较浅,包气带相对较薄;相反,地形指数小的地方基本位于流域的上游山坡,远离河道,包气带相对较厚.因此,可以假设流域上地形指数最大的栅格单元对应的张力水蓄水容量最小,而地形指数最小的栅格单元对应的张力水蓄水容量最大.图2为估计的屯溪流域包气带厚度以及腐殖质土厚度的空间分布.

图2 屯溪流域土壤厚度分布Fig.2 Distribution of thickness of soil in Tunxi Basin

K i与K g这2个参数属于并联参数,其和K i+K g代表的是自由水出流的快慢,应与单元的土壤类型有关,而自由水指的是饱和含水量与田间持水量之间那部分可以在重力作用下自由流动的水,因此可以将 θs与θfc作为衡量自由水出流快慢的指标.Ki/Kg表示的是壤中流与地下径流的比,此比值可以通过 θwp来反映[14].具体计算公式为

式中:moc——自由水出流综合影响因子;mr——自由水出流校正系数,根据已有研究成果,在估计模型参数时,可以取mr=1.moc可直接通过新安江模型分水源中的结构性约束(Ki+Kg=0.7)[5]选取.具体步骤是:先给moc赋予一初值,然后根据每个栅格单元的土壤类型由土壤参数统计表和式(6)确定栅格单元对应的Ki+Kg,再统计出流域内所有栅格Ki+Kg的均值并与0.7进行比较,由此对moc进行调整,使其尽量满足该结构性约束.

3.3 汇流参数

当Grid-Xin'anjiang模型进行扩散波汇流演算时,模型参数的估计主要是基于流域地貌特征以及河道断面信息.其中,坡面汇流的曼宁糙率系数n h根据栅格单元的植被类型由植被参数统计表确定.而每个栅格单元河道(包括虚拟河道)汇流的曼宁糙率系数nc[15]由下式进行计算:

式中:Ad——栅格单元的上游汇水面积,km2;Soc——河道坡度,可由DEM提取得到;n0——系数,n0可以由流域出口点的糙率及坡度代入式(7)反算出.地表径流汇入河道的比例也是在生成的研究流域水系基础上,采用面积比例法进行计算[7].在判断河道形状时,主要是根据流域内实测站点的断面数据,反演出每一栅格单元的河道形状.以屯溪流域为例,定义某一断面宽度指数α,考虑到随着上游汇水面积的增加,越到流域下游其河道过水断面面积应该越大,且变化比较明显,因此可认为 α在流域内应当是变化的,即河道断面的尺寸是空间变化的.根据流域内已有的实测断面资料分析,即可获得α的空间分布.

4 应用结果分析

屯溪流域具有1982—2001年以及2006年的日降雨资料,有40场洪水的时段降雨资料,屯溪出口站有对应的流量资料.本文将1982—1995年的日资料和对应的次洪资料用于模型参数的率定,剩余资料用于模型参数的检验,对于模型参数估计方法的验证采用率定期资料进行,对于屯溪流域降雨资料的空间插值采用距离平方反比法.

采用率定期内的资料对模型估计参数进行验证.其中,日模型模拟的径流深以及洪峰相对误差的合格率分别为100%和79%,平均相对误差水平分别为3.7%和12.6%,确定性系数的均值为0.91.次洪模型模拟的径流深、洪峰相对误差以及峰现时差的合格率分别为91%,48%和83%,平均相对误差水平分别为9.1%,19.9%和2.2h,确定性系数的均值为0.87.图3为应用估计参数模拟的2次洪水过程线比较.

图3 Grid-Xin'anjiang次洪模型估计参数模拟结果比较Fig.3 Comparison of estimated parameters of Grid-Xin'an jiang model for floods

由估计参数的模拟结果可以看出,根据土壤、植被及地貌等下垫面条件估计出的Grid-Xin'anjiang日模型参数值较为合理,对于屯溪流域率定期内的日径流模拟精度较高,因此直接将日模型的参数估值用于该流域检验期内的日径流模拟,检验期内模拟的径流深与洪峰相对误差合格率均为91%,平均相对误差水平分别为10.1%和7.0%,确定性系数的均值为0.90.对于Grid-Xin'anjiang次洪模型而言,需要率定检验的主要是扩散波汇流方法中的坡面糙率系数nh,率定方法选择人工试错法.对于次洪模型中其他参数则不再进行率定.通过率定扩散波汇流方法的糙率系数,对于率定期内洪水而言,Grid-Xin'anjiang次洪模型模拟的径流深、洪峰和峰现时间的平均相对误差水平分别为7.7%,7.9%和1.8 h;合格率分别为96%,96%和91%,确定性系数的均值为0.91.对于检验期内洪水而言,Grid-Xin'anjiang次洪模型模拟的径流深、洪峰和峰现时间的平均相对误差水平分别为8.0%,9.7%和2.6h;合格率均为94%,确定性系数的均值为0.90,模拟结果较好.图4为应用率定参数模拟的2次洪水过程线比较.

图4 Grid-Xin'anjiang次洪模型率定参数模拟结果比较Fig.4 Com parison of calibrated parameters of Grid-Xin'an jiang model for floods

5 结 语

以新安江蓄满产流模型与一维扩散波汇流模型为基础,构建了适用于湿润流域的分布式Grid-Xin'anjiang模型,并在模型进行产汇流计算时,考虑了栅格间的水量交换问题以及河道排水网络的影响.根据流域的地貌特征、土壤类型及植被覆盖等下垫面特性,对基于物理机制的模型参数估计方法进行了一定的研究与验证.从Grid-Xin'anjiang模型在屯溪流域的应用情况可以看出,模型的应用效果较好,模拟精度较高.为了使Grid-Xin'anjiang模型可以更好地应用于不同类型的流域,还应当加入超渗产流、融雪模型等其他水文过程的计算方法,对于栅格间壤中流与地下径流的汇流演算方法也有待进一步改进.为了使模型更好地用于缺乏资料地区的洪水模拟及预报,还需要进一步加强对模型的应用研究,开展模型在其他流域以及与其他模型应用的对比分析,同时也需要进一步完善模型参数的估计方法.

[1]FREEZE R,HARLAN R.Blueprint for a physically-based,digitally-simulated hydrologic response model[J].Journal of Hydrology,1969,9:237-258.

[2]ABBOTT M,BATHURST J,CUNGE J,et al.An introduction to the European Hydrological System-Système Hydrologique Européen,“SHE”,1:history and philosophy of a physically-based,distributed modelling system[J].Journal of Hydrology,1986,87:45-59.

[3]TODINI E,CIARAPICA L.The TOPKAPI model[M]//SINGH V.Mathematical models of large watershed hydrology.Colorado:Water Resources Publications,2001:471-504.

[4]BEVEN K,KIR KBY M.A physically based variable contributing area model of basin hydrology[J].Hydrology Science Bulletin,1979,24(1):43-69.

[5]赵人俊.流域水文模拟[M].北京:水利电力出版社,1984.

[6]李致家,姚成,汪中华.基于栅格的新安江模型的构建和应用[J].河海大学学报:自然科学版,2007,35(2):131-134.(LI Zhijia,YAO Cheng,WANG Zhong-hua.Development and application of grid-based Xin'anjiang model[J].Journal of Hohai University:Nature Sciences,2007,35(2):131-134.(in Chinese))

[7]YAO Cheng,LI Zhi-jia,BAO Hong-jun,et al.Application of a developed Grid-Xin'anjiang model to Chinese watersheds for flood forecasting purpose[J].Journal of Hydrologic Engineering,2009,14(9):923-934.

[8]芮孝芳.水文学原理[M].北京:中国水利水电出版社,2004.

[9]ASTON A.Rainfall intercep tion by eight small trees[J].Journal of Hydrology,1979,42:383-396.

[10]赵人俊,王佩兰.新安江模型参数的分析[J].水文,1988(6):2-9.(ZHAO Ren-jun,WANG Pei-lan.Analysis of the parameters of the Xin'anjiang model[J].Journal of China Hydrology,1988(6):2-9.(in Chinese))

[11]MACCORMACK R.Numerical solution of interaction of shock wave with laminar boundary layer[J].Lecture Notes in Physics,1971,8:151-163.

[12]WANG M,HJELMFELT A.DEM based overland flow routing model[J].Journal of Hydrologic Engineering,1998,3(1):1-8.

[13]李致家.水文模型的应用与研究[M].南京:河海大学出版社,2008.

[14]KOREN V,SMITH M,WANG D,et al.Use of soil property data in the derivation of conceptual rainfall-runoff model parameters[C]//Proceedings of the 15th Conference on Hydrology.Long Beach:AMS,2000:103-106.

[15]KOREN V,REED S,SMITH M,et al.Hydrology laboratory research modeling system(HL-RMS)of the national weather service[J].Journal of Hydrology,2004,291:297-318.