基于蓄满超渗时空动态组合的网格新安江模型

杜若愚，姚 成，刘玉环，李致家，张 珂，朱跃龙

(1.河海大学水文水资源学院，江苏 南京 210098； 2.河海大学计算机与信息学院，江苏 南京 210098)

产流模式主要分为蓄满产流和超渗产流[1]，我国南方湿润地区以蓄满产流为主，北方干旱地区以超渗产流为主，半干旱半湿润地区的产流模式两者兼而有之[2-3]。半干旱半湿润地区降雨时空分布不均，下垫面条件空间异质性大，产流模式混合多变，即蓄满产流和超渗产流随时空变化明显，导致多数水文模型洪水预报结果欠佳[4-5]。

半干旱半湿润地区，是湿润地区与干旱地区的过渡地带，如果降雨落在具有很高下渗能力的土壤，就会像湿润地区那样发生壤中流和地下径流，土壤蓄满后发生饱和地面径流，即蓄满产流；而当降雨落在下渗能力小，再加上土壤保水能力差，就会发生超渗产流，这就在空间尺度上形成了蓄满产流与超渗产流同时发生的现象。而半干旱半湿润地区洪水预报的难点在于确定产流模式，李致家等[6-8]认为产流模式是随时空动态变化的，并在新安江模型的基础上添加了超渗产流结构，完善了产流理论，使其能够适应半干旱半湿润地区复杂条件的混合产流洪水预报；李致家等[9]对基于子流域的蓄满超渗空间组合进行研究，分别采用新安江模型和河北雨洪模型计算不同产流区的产流[10]，但这种产流计算为空间静态组合，因此，Liu等[11-12]研究了基于网格的蓄满超渗时空组合，考虑了产流模式随时空变化的特点，基于TOKASIDE模型构建TOKASIDE-D模型，并成功应用于我国华北地区洪水预报。

TOKASIDE-D模型是基于物理基础的分布式水文模型，其考虑了地形与运动波的蓄满超渗产流、地下水运动、水库入流、调蓄计算和渗漏等模块，核心思路是产流模式随时空的动态变化。由于网格新安江(Grid-XAJ)模型[13-18]自创建以来在全国多数流域得到了广泛应用，因此本文借鉴TOKASIDE-D模型中产流模式时空动态组合的思路，构建基于蓄满超渗时空动态组合的网格新安江(Grid-XAJ-SIDE)模型，并将Grid-XAJ、Grid-XAJ-SIDE和Grid-GA[19](网格格林-安普特)模型应用于半干旱地区绥德流域进行验证和分析，探讨蓄满、超渗产流模式时空动态组合方法对流域洪水模拟结果的影响。

1 模 型 构 建

1.1 Grid-XAJ模型

Grid-XAJ模型以DEM网格为计算单元，主要分为蒸散发、蓄满产流、分水源、汇流4个计算模块[13-18]。模型假设每个网格的降雨和下垫面条件空间分布均匀，不存在张力水蓄水容量分布曲线和自由水蓄水容量分布曲线，只需考虑各要素在不同网格间的变异性。Grid-XAJ模型是蓄满产流型分布式水文模型，广泛应用于我国湿润地区洪水预报。

1.2 Grid-GA模型

Grid-GA模型[19]主要分为蒸散发模块、产流模块和汇流模块3个计算模块，各个模块均在以DEM网格为基础的正交网格内进行。模型认为每个网格的降雨和下垫面条件空间分布均匀，不存在下渗能力分布曲线，Grid-GA模型是超渗产流型分布式水文模型，主要用于干旱半干旱地区洪水预报。

1.3 Grid-XAJ-SIDE模型构建

Grid-XAJ-SIDE模型是在Grid-XAJ模型的基础上，使用径流曲线数-地形指数法划分流域初始蓄超产流区分布，并添加超渗产流模块，同时在产流计算过程中动态识别蓄满网格和超渗网格，分别采用Grid-XAJ模型的产流与分水源和Green-Ampt下渗公式计算产流，最后根据网格间的汇流演算次序，依次将不同径流成分演算至流域出口。

Grid-XAJ-SIDE模型主要分为蒸散发、蓄超网格动态识别、产流、分水源、汇流5个计算模块，模型计算流程如图1所示(图中Pt为降雨量，Ipt为降雨强度，Wt为网格土壤含水量，WM为网格张力水蓄水容量，ft为下渗能力，t为计算时段，n为总时段数)。

图1 Grid-XAJ-SIDE模型计算流程Fig.1 Calculation flow-process of Grid-XAJ-SIDE model

Grid-XAJ-SIDE模型基于网格新安江模型假设，并添加以下假设条件：①网格内降雨、土壤性质、土地利用等特征均匀分布，不存在下渗能力分布曲线；②流域内网格在计算过程中，依据土壤含水量是否达到田间持水量，降雨是否超过下渗能力，动态识别蓄满网格或超渗网格；③每个网格的土壤垂向分布均匀。

1.3.1 蒸散发计算

Grid-XAJ-SIDE模型采用3层蒸散发模型计算流域实际蒸散发量。3层蒸散发模型将每个网格的土壤分为上层、下层和深层，对应的土壤蓄水容量分别为WUM、WLM和WDM，具体计算公式详见文献[1]。

1.3.2 初始蓄超网格判定

Grid-XAJ-SIDE模型基于径流曲线数和地形指数，对流域的蓄满网格和超渗网格的初始空间分布进行确定，径流曲线数源自SCS模型，表示某种土壤水分条件下的曲线数，主要用来描述地面超渗产流发生的难易程度，径流曲线数越大，发生超渗径流的可能性越大。初始蓄超网格判定采用的数据主要包含DEM、土壤类型和土地利用空间分布。划分方法如下：计算每个网格内的径流曲线数和地形指数，再根据径流曲线数进行初始分类[20]，当径流曲线数小于60时，将网格划分为蓄满网格，否则划分为超渗网格。在此基础上，再加入地形指数进一步修正，当蓄满网格的地形指数小于7，将该网格修改为超渗网格；当网格的地形指数大于25，将该网格修改为蓄满网格，最终得到整个流域的初始蓄超网格空间分布。

1.3.3 蓄超网格动态识别

与初始蓄超网格判定不同，蓄超网格动态识别是Grid-XAJ-SIDE模型在产流计算过程中，根据当前时段网格单元的土壤含水量、降雨强度与土壤下渗能力之间的关系，动态识别蓄满网格和超渗网格，是一种动态过程；而初始蓄超网格判定是在模型开始计算之前，基于径流曲线数-地形指数法，根据下垫面信息进行判定，是一种静态过程。

蓄超网格动态识别的主要原则：当降雨强度大于下渗能力时，该网格为超渗网格；当网格的土壤含水量达到田间持水量后，该网格则为蓄满网格。其中，流域初始土壤含水量的值由日模型计算得到。

1.3.4 产流及分水源计算

1.3.4.1 蓄满产流及分水源

Grid-XAJ-SIDE模型对于蓄满网格采用蓄满产流模式，即在降雨过程中，直到土壤含水量达到田间持水量才能产流，而在此之前，所有来水均用于补充土壤含水量。将计算时段内网格单元的实测降雨先扣除相应时段的蒸散发量，即可得到实际用于产流计算的时段雨量Pe，则

(1)

式中：R为时段产流量，mm；W0为上一时段土壤含水量，mm。

在Grid-XAJ-SIDE模型中，每个网格单元的R均被划分为3种径流成分即地表径流Rs、壤中流Ri以及地下径流Rg。与产流计算一样，在进行分水源计算时，认为每个网格内自由水蓄水容量分布均匀。分水源计算公式为

Ri=KiS

(2)

Rg=KgS

(3)

(4)

式中：Ki为自由水含量对壤中流的出流系数；Kg为自由水含量对地下水的出流系数；S为自由水含量，mm；SM为自由水蓄水容量，mm。

1.3.4.2 超渗产流

对于超渗网格，采用Green-Ampt下渗公式计算网格下渗能力：

(5)

式中：f(t)为土壤下渗能力，mm/h；Ks为饱和水力传导度，mm/h；Ψ为湿润锋处的土壤吸力，mm；Δθ为饱和含水率与初始含水率之差；F(t)为累计下渗量，mm。计算时无需考虑下渗能力分布曲线，根据超渗产流原理，当降雨强度大于土壤下渗能力时，以下渗能力下渗，超出下渗能力的部分降雨形成地表径流，计算步骤如下：① 根据式(5)计算各网格下渗能力，其中，初始土壤含水率由日模型计算得到，初始土壤F(t=0)赋极小值0.000 01，即认为初始时刻的下渗能力无限大；②计算每个网格的R：

(6)

③计算该网格当前时段的下渗量：i=Pe-R。该网格τ时段的累积下渗量为：F(t=τ)=F(t=τ-1)+i(t=τ)，其中i为当前时段下渗量，mm。将F(t=τ)代入式(5)，计算下一时段的土壤下渗能力。每个网格存在一个最大累积下渗量Fmax，当土壤累积下渗量大于最大累积下渗量时，土壤含水量达到田间持水量，多出的水量形成地下径流。④重复步骤①，依次计算流域内每个网格在各时段下的径流量。

1.3.4.3 土壤含水量更新

Grid-XAJ-SIDE模型中，蓄满网格土壤含水量按照蓄满产流算法进行更新。超渗网格土壤含水量更新方式：根据超渗产流原理，当降雨强度小于下渗能力时，降雨全部下渗到土壤中，地面不产生积水；当降雨强度大于下渗能力时，降雨分为两部分，一部分按照下渗能力下渗到土壤中，补充土壤含水量；另一部分则留在了地面，成为超渗地表径流。设网格当前时段下渗量为i，上一时段上层、下层和深层的土壤含水量分别为WU0、WL0和WD0，三层土壤的蓄水容量分别为WUM、WLM和WDM，按照先上层后下层的次序更新土壤含水量。入渗量优先补充上层土壤含水量，若上层土壤达到蓄满，多余的水量继续下渗补充下层土壤含水量，如果下层土壤同样已经蓄满，多余的水量继续补充深层土壤含水量，直至整个土层土壤含水量达到田间持水量。

1.3.5 汇流计算

Grid-XAJ-SIDE模型汇流分为逐网格坡面汇流和河道汇流2个阶段计算。坡面汇流采用一维扩散波方程计算，使用基于两步法MacCormack算法[21]的二阶显式有限差分格式进行扩散波方程组的求解；河道汇流采用基于网格的马斯京根河道汇流演算法[22]进行计算。

1.3.6 模型参数

Grid-XAJ-SIDE模型以Grid-XAJ模型为基础，并添加了超渗产流计算模块，实现了蓄超产流时空动态组合的定量模拟。由于下垫面的空间分布具有异质性，研究流域的网格参数与下垫面的土壤、土地利用、高程等分布具有紧密相关性，本文涉及的3个模型均为分布式水文模型，其网格的上层张力蓄水容量WUM、下层张力蓄水容量WLM、张力水蓄水容量WM、自由水蓄水容量SM、地下径流出流系数KG、壤中流出流系数KI、土壤总孔隙度θe、土壤有效孔隙度θs、湿润锋处土壤吸力ψ和饱和水力传导度Ks等参数的取值也为分布式。其中WUM、WLM、WM、SM、KG和KI均根据Grid-XAJ参数及其空间分布估计方法进行计算[23-24]；每个网格的θe、θs、ψ和Ks等超渗产流参数，均与土壤类型相关，根据Rawls等[25-27]的研究确定每个网格的参数值；蒸发折算系数K、深层蒸散发系数C、最大累积下渗量IO、地下径流消退系数CG、壤中流消退系数CI、河网蓄水消退系数CS、滞后时间Lag和河道径流逐网格马斯京根法演算参数Kech、Xech等参数则基于流域实测水文资料，采用人工优选法在经验取值范围内进行率定。

2 模 型 验 证

2.1 研究流域

选择半干旱地区的绥德流域作为研究流域，流域高程、水系及站点分布如图2所示。绥德水文站位于陕西省绥德县，东经110°14′，北纬37°30′，建站于1959年6月，流域控制面积3 893 km2。该站多年平均降雨量为486 mm，多年平均气温为9.7 ℃，实测最大流量为2 350 m3/s，相应水位为818.59 m。绥德流域地形东西跨度大，主要以山地为主，总体呈现西高东低的趋势；土壤以壤质砂土为主。

图2 绥德流域地形、水系及站点分布Fig.2 Topography, river and station distribution of Suide Watershed

2.2 数据资料

Grid-XAJ-SIDE模型所需下垫面数据主要包括DEM、土壤类型和土地利用类型。DEM数据选择由美国太空总署(NASA)与国防部国家测绘局(NIMA)联合测量的SRTM数据，分辨率为30 m，通过地理空间数据云(http://www.gscloud.cn)下载。土壤类型数据来源于联合国粮农组织(FAO)和维也纳国际应用系统研究所(IIASA)构建的世界土壤数据库HWSD (harmonized world soil database version)1∶100万土壤栅格数据。土地利用类型数据通过美国地质勘探局(United States Geological Survey，USGS)下载，数据是分辨率为1 km 的矢量文件。

以1 h为计算步长建立模型，选取2010—2018 年站点资料较为齐全且洪峰较大的15 场洪水资料进行计算，将收集到的实测降雨和流量资料插值为1h步长数据，其中前10场洪水资料用于模型参数率定，后5场洪水用于模型验证。

2.3 验证结果与分析

将Grid-XAJ模型、Grid-GA模型和Grid-XAJ-SIDE模型应用于绥德流域进行应用检验和对比分析时，网格大小均采用1 km×1 km，既可以更好地与土壤、植被等数据的网格大小保持一致，也可以在保证应用精度的前提下，提高模型运算效率。Grid-GA模型和Grid-XAJ-SIDE模型对于超渗产流的计算，计算时段会影响模拟结果，但绥德流域原始资料大多为2 h步长，目前受限于流域内观测资料的精度，3个模型计算步长均选择1 h。

依据GB/T 22482—2008《水文情报预报规范》[28]评价模型模拟结果，选择径流深误差、洪峰相对误差和确定性系数等洪水模拟指标进行综合评价。3个模型人工优选参数见表1，率定和验证的精度见表2。

表1 绥德流域各模型参数值

表2 绥德流域各模型模拟结果对比

对于径流深的模拟，率定期Grid-XAJ-SIDE模型、Grid-XAJ模型和Grid-GA 模型模拟合格率分别为90.0%、80.0%和100.0%，验证期合格率分别为80.0%、80.0%和100.0%，模拟结果均较好，主要由于绥德流域实测径流深普遍偏小，15场洪水中，只有20170723号洪水实测径流深大于15 mm，且3个模型对于径流深的模拟值与实测值误差多数在3 mm以内，依据GB/T 22482—2008《水文情报预报规范》[28]，3个模型的径流深合格率都比较高。绥德流域地处半干旱区，洪水历时短、涨落快、总径流深小但洪峰较大，因此准确模拟洪峰更为重要。在洪峰合格率方面，Grid-XAJ-SIDE模型率定期和验证期模拟结果均最好，分别为70.0%和60.0%，相比于Grid-XAJ模型(率定期60.0%，验证期40.0%)和Grid-GA模型(率定期30.0%，验证期60.0%)有明显提高。3个模型对于峰现时间的模拟结果相同。3个模型的确定性系数普遍较差，说明在绥德这样的半干旱流域，确定性系数不能用来描述陡涨陡落的洪水过程。

20130726号洪水事件中(图3(a))，Grid-XAJ-SIDE模型和Grid-XAJ模型模拟结果接近实测过程，此次洪水的前期土壤蓄水程度较高，缺水量小，是以蓄满产流为主导的洪水过程，初始土壤饱和度如图4(a)所示。而Grid-GA模型模拟结果与实测过程相差较大，这是因为Grid-GA模型认为土壤是半无限土柱，即不会达到“蓄满”状态，入渗到土壤中的降雨不会再产生径流。在20180710号洪水事件中(图3(b))，Grid-XAJ-SIDE模型和Grid-GA模型模拟结果明显优于Grid-XAJ模型，这是因为此次洪水的前期土壤蓄水程度低，缺水量大，不易蓄满，是以超渗产流为主导的洪水过程，初始土壤饱和度如图4(b)所示。此次洪水过程中Grid-XAJ-SIDE模型和Grid-GA模型在洪水起涨阶段模拟结果相近，且与实测过程接近，但在退水阶段，Grid-GA模型退水过程过快，与实测过程偏差较大。Grid-XAJ-SIDE模型认为，即使是超渗产流主导的网格，当下渗量补充土壤含水量达到“蓄满”状态后，下渗的水量会产生壤中流和地下径流，因而比Grid-GA模型更接近实测退水过程。

图3 绥德流域部分洪水实测过程线与模拟过程线对比Fig.3 Comparison of observed and simulated hydrographs of partial flood events in Suide Watershed

图4 绥德流域部分洪水场次初始土壤饱和度分布Fig.4 Initial soil saturation distribution of partial flood events in Suide Watershed

Grid-XAJ-SIDE模型通过对逐时刻的土壤含水量和降雨等因子进行实时“监控”，动态调整网格的产流方式。以20130726号洪水和20180710号洪水过程为例，展示在模型计算过程中蓄超网格的动态变化，如图5所示。20130726号洪水过程的整体蓄水程度较高(图4(a))，容易蓄满，因此本场洪水从7月26日16:00开始降雨之后，蓄满网格数量便迅速增多(图5(a))，到7月27日2:00达到最大，此时模拟洪水达到顶峰，之后进入退水阶段。从图4(b)和图5(b)可以看出，20180710号洪水过程整体蓄水程度不高，且下渗较小，由于降雨强度大，历时短，随着降雨的进行，从开始降雨(20180710T10:00)到洪峰时刻(20180711T14:00)蓄满网格减少，整个洪水过程蓄满网格数量少于超渗网格，是以超渗网格为主导的洪水过程，进而模拟的洪水过程也呈尖瘦形态。

图5 绥德流域部分洪水蓄满和超渗网格空间动态分布Fig.5 Spatial dynamic distribution of saturation-excess and infiltration-excess grids in Suide Watershed

Grid-XAJ-SIDE模型以网格为计算单元进行产汇流计算，并基于下垫面特性的分布信息推求模型参数的空间分布，因此模型在输出流域出口断面流量的同时，可以输出任意网格单元的流量过程。以绥德流域20180710号洪水过程为例，由图6可以看出，此次洪水事件前期主要发生在流域上游，随着洪水过程的进行逐渐转移至中下游，以洪水过程流量空间动态分布为依据，不仅可以反映流域各点的流量过程，还可以推测暴雨中心的转移过程，也从侧面反映出半干旱地区降雨空间分布的不均匀程度。

图6 绥德流域20180710号洪水过程流量空间动态分布Fig.6 Spatial dynamic distribution of discharge for flood event 20180710 in Suide Watershed

通过3个模型的对比验证与分析，Grid-XAJ-SIDE模型结构完善，结合了蓄满型产流和超渗产流模型的优点，在下垫面条件比较复杂的半干旱地区提升了洪水模拟精度，对流域内的洪水发生过程实现了精细描述。

3 结 语

在网格新安江模型(Grid-XAJ)的基础上，基于径流曲线数和地形指数划分流域初始蓄超分布，构建基于蓄满超渗时空动态组合的网格新安江(Grid-XAJ-SIDE)模型，并以半干旱地区绥德流域实测数据对模型进行了验证。研究结果表明，与Grid-XAJ模型和Grid-GA模型相比，Grid-XAJ-SIDE模型根据土壤含水量和降雨等动态产流影响因子的实时变化，更加精确地描述了网格内的产流模式，更好地结合了蓄满产流模型和超渗产流模型的优点，使洪水模拟精度有了较为明显的提升，对半干旱半湿润地区洪水预报具有积极意义。