陈 军，刘 意，李 婷

(成都信息工程大学资源与环境学院,成都 610225)

0 引 言

近年来，因人类活动导致气候异常，暴雨洪涝灾害频发。为减少暴雨洪涝导致的财产和人员损失，亟须利用径流汇流模型对流域进行暴雨洪涝灾害预测。基于物理机制的径流汇流模型一般以网格为单元进行汇流模拟，可获取全流域的水流场空间和时间变化，逐渐成为暴雨洪涝灾害预测的重要手段[1-6]。

然而，流域汇流是一个非常复杂的过程，除降水之外，还存在上游河道来水、冰川融水、地下水补给等。由于上述因素难以精确地在流域上按网格单元量化，现有径流汇流模型采用近似的方法进行模拟。以SWAT模型为例，它将多年的降雨数据代入到无水流域中反复迭代，将模型的最终输出作为流域前期水情。这种方式不仅运算量大，还存在模拟误差的迭代问题，因而所得到的前期水情与实际存在较大差距[7-11]；作为当前著名的水动力模型之一，FloodArea由于能提供流域任一栅格单元的水深模拟结果，被广泛应用于暴雨洪涝模拟、山洪风险分析[12-14]，但该模型在径流模拟时，前期水情仅作为输入参数由用户提供。

本文以纳维-斯托克斯方程(N-S方程)为指导，提出一种基于稳定流场的流域径流汇流模拟方法，该方法以稳定流场作为前期水情生成和径流汇流模拟的基本依据。最后，以岷江上域镇江关-校场坝段为例，分析方法的有效性。

1 研究区概况

岷江上游位于四川省北部，流域主体主要位于松潘县和茂县，地形以山地地貌为主，地势北高南低。在西风南支急流、暖湿东南季风以及地形条件的共同影响下，其气候特征主要表现为冬干夏湿，干湿季节分明，属于非地带性的干旱河谷气候。流域年降水量800～1 200 mm，夏季降水高度集中，降水曲线呈单峰变化，镇江关以上春季有极少量的雪水补入，年内径流季节变化与降雨季节相应。岷江上游多年平均年径流量松潘县镇江水文站约为18 亿m3，汇口下不远的五通桥水文站多年平均年径流量约780 亿m3。区域大面积暴雨发生或随洪峰传播，干流峰高量大，严重威胁中下游安全[15]。为降低人类活动对径流汇流模型精度验证的影响，选取岷江上游的镇江关和校场坝两个水文站，将其所在河道位置分别作为流域出口，计算流域重叠范围作为研究流域，如图1所示。流域上游入水口为镇江关水文站，下游出口为校场坝水文站，干流河道总长约35 km，汇流面积约为580 km2。

图1 岷江上游研究区域Fig.1 The study area of the upper reaches of Minjiang River

2 基于纳维-斯托克斯方程的径流汇流模型

纳维-斯托克斯方程(Navier-Stokes equations)是描述黏性不可压缩流体动量守恒的经典运动方程，反映了黏性流体(又称真实流体)流动的基本力学规律，在流体力学中有十分重要的意义[16]。但它是一个非线性偏微分方程，只有在极少数简单流动的情况下才有解析解[17-20]。本文以N-S为指导，探索栅格模式下的径流汇流模拟。

2.1 水流运动的N-S方程

N-S方程是一个用来描述流体运动特征的最基本力学方程组。由于水流运动是同性不可压的流体运动，可以视为是无黏性的流体，N-S公式表达为方程组：

(1)

▽μ=0

(2)

式中：u为速度场；t为时间；ρ为流体密度；p为压强；F是作用于流体上的质量力(外力)。

在径流汇流模拟过程中，流体为实际的水体，采用水深和水速描述水体的属性。N-S方程中通过外力项和压力项求解以模拟水体在重力、摩擦力、压力梯度等作用下的水速大小及变化；通过平流项求解以模拟水流因流动在流域各位置上水深大小及变化。

2.2 栅格模式下的N-S方程分解

栅格模式的径流汇流模拟将流域从空间上和时间上均划分为相等的片段。在空间上，按一定分辨率将流域划分为大小相等的矩形网格，将坡面汇流和河道汇流统一描述为网格的水流运动；在时间上，径流汇流过程分解为相等间隔的时间片。在每一个时间片上，流域径流汇流过程分解为栅格单元为基本单位的模拟计算，上一个时间片栅格单元模拟的水深和水速，作为下一个时间片初始水深和水速。通过时间片的迭代，实现流域的径流汇流过程模拟。

设流域中任一栅格单元(中心栅格)为c，时间片初始时间为tn，终止时间为tn+1，时间片的时间间隔为Δt，则：

Δt=tn+1-tn

(3)

中心栅格c的径流汇流模拟需要依据tn时刻中心栅格及8邻域栅格单元b的水深和水速，模拟得到tn+1时刻中心栅格的水深和水速。

(4)

(5)

式中：Δx、Δy为邻域栅格坐标相对中心栅格坐标的偏移，单位为像素；vmax表示时间片上允许的最大速度标量。

中心栅格剩余水量比例Sc为：

(6)

表1给出了栅格单元不同的水速方向上，流出水量和剩余水量与初始水量之间的关系。从表1可见，以式(5)和式(6)为基础的多流向法平流项模拟算法满足质量守恒定律。

表1 不同水速方向上流出水量和剩余水量关系表Tab.1 The relation table of outflow and remaining water amount in different water velocity directions

依据式(5)和式(6)，平流计算后中心栅格水深为：

(7)

(8)

依据式(5)和式(6)，当模拟网格水速大小等于vmax时，时间片上的网格水流全部流出。因此，vmax决定了水质模拟的时间片间隔，如下式：

(9)

式中：C为栅格分辨率。

从式(9)可见，vmax越大，时间片时间间隔越小，时间精度越高，但1 h需要迭代的次数越多；另外，为保证水量平衡，当栅格单元c的模拟速度大于vmax时，设置为vmax。

(3)中心栅格c的水动力模拟。水动力模拟对应于N-S方程的压力项和外力项求解。设中心栅格高程为Hc，邻域栅格高程为Hb，则水体高度差ΔHb为：

(10)

(11)

式中：α为与水体密度、重力加速度、摩擦力等相关的正数常量，通过流域实际水文数据进行率定。

地表水流由于受到地表摩擦力的影响，使水流速度发生垂直变化。在径流汇流模拟过程中，当水深低于给定阈值时，将水速乘以衰减系数以减小水流速度。衰减系数ε近似定义为：

(12)

式中：dmin和dmax分别为水流摩擦力的下界深度和上界深度，m；σ为摩擦比例系数，通过流域实际水文数据进行率定。

3 基于稳定流场的径流汇流模拟

3.1 稳定流场的基本概念

假定流域上游持续不间断地提供稳定水源，通过径流汇流模型不断迭代，流域出口处的水深和水速最终趋于稳定。此时的流域水流场称为稳定水流场。

显然，现实中不存在稳定水源，稳定流场也不可能存在。但在流域上游设定稳定供应水源，通过模型迭代，如果能在流域出口处最终模拟出稳定的水深和水速，则径流汇流模型具有一定的合理性。

设上游来水点稳定供水量为ξ(单位为m3/h)，则考虑上游来水的栅格单元水量收支由式(4)修改为式(13)：

(13)

式中：ξc为时间片上栅格单元分摊的上游来水量。

根据上游来水的定义，在上游来水点附近，ξc>0，其余位置上ξc=0，并且ξc满足：

ξ=∑ξcN

(14)

式中：N为小时需迭代的时间片数，由式(9)换算得到。

合理的径流汇流模型是稳定水流场建立的关键。稳定水流场创建过程中，采用式(15)定量描述稳定流场是否创建完成。

(15)

式中：W为流域下游控制点水深；n为当前已模拟的小时数；h为前推小时数；A为累积误差阈值。

3.2 径流汇流模型完善与稳定流场生成

直接按N-S方程分解得到的径流汇流模型存在“汇流水流方向异常”问题。在流水不畅的河道区域，如果模拟水速过大，上游栅格水流快速向中心栅格聚集，中心栅格水流又不能及时向下游流动，导致短时间内产生较大的水体高差，进而对周围8邻域栅格产生更大速度增量。随着迭代次数的增加，这些位置上很容易产生水流方向异常，甚至出现相邻上下游栅格流速方向相反并频繁交替的现象，导致水深和水速变化剧烈。表2展示了模型完善前某栅格单元典型时间片上的栅格水速和水深变化。从表2可见，现有模型违背了现实径流汇流规律，即使在上游稳定供水情况下也难以生成相对稳定的水流场(如图2所示的模型完善前流域出口水深和水速曲线)，从而不能精确模拟流域的径流汇流过程。

为避免上述问题，在径流汇流模型中添加水流混合模拟。水流混合模拟由栅格单元水速增量修正和邻域栅格单元水流混合模拟两步组成。

(1)栅格单元水速增量修正。假定水深大于1 m的栅格水流与周围水团混合后水速增量减小，并与水深呈反比关系。其修正公式为式(16)：

(16)

(2)邻域栅格单元水流混合模拟。假定每一个时间片上，各栅格单元均与8邻域栅格水流混合，混合过程中遵循动量守恒定律。首先计算中心栅格和8邻域栅格的总动量和总水量，两者相除得到中心栅格水流混合后的水速。

表2 模型完善前后典型栅格水速和水深变化Tab.2 The water speed and depth change of typical grid before and after the completion of runoff model

从表2可见，模型完善后，典型栅格的水流不断累积，水深稳定增加，水速逐渐减小，水流方向保持稳定，符合径流汇流的实际规律。

以岷江上游2018年5月20日12时镇江关水文站瞬时流量作为岷江干流入流量(ξc=86.4)，在h=10、A=0.05的条件下,迭代15 h后，稳定水流场生成，如图3所示的模型完善后流域出口的水深和水速曲线。相对于模型完善前流域出口在上游稳定供水条件下仍然有约0.08 m幅度的水位变化(图2)，完善后的模型水深和水速均趋于稳定，为基于稳定流场的前期水情生成和径流汇流模拟奠定了基础。

图2 模型完善前流域出口水深与水速变化Fig.2 The change of water depth and velocity at pour point before the model was modified

图3 模型完善后流域出口水深与水速变化Fig.3 The change of water depth and velocity at pour point after the model was modified

3.3 基于稳定流场的径流汇流模拟

为客观模拟流域前期水情，依据流域汇流时间，在上游河道输入模拟时段前数小时的流量序列，通过模型迭代建立水流场作为流域前期水情。在模拟流域前期水情过程中，首先将前期流量序列的第一个小时流量代入模型作为上游稳定供水，创建流域的稳定流场；在此基础上，将前期上游流量序列和面雨量序列逐小时迭代，模拟流域的前期水情。

以流域前期水情为基础，将模拟时段逐小时上游流量序列和面雨量序列代入模型中迭代，模拟流域的径流汇流过程。

为提高径流汇流模拟精度，利用校场坝水文站水位数据对模型参数分别率定。在率定某一个参数时，其余参数设为固定值，取与实际水位误差最小的参数值作为率定结果。率定结果如表3所示。

表3 模型参数率定结果Tab.3 The calibration results of model parameters

4 径流汇流模拟精度评价

为客观评价径流汇流模拟的精度，将径流汇流模拟按无稳定流场无上游来水(A)、无稳定流场有上游来水(B)、稳定流场(C)3种方式，以2018年5月21日12时至2018年6月15日11时的长时间序列(600 h)为例，分别做出3种模拟方法的流域出口校场坝的水位曲线，并与实际水位曲线进行对比，制作小时水位曲线(图4)和误差统计表(表4)。其中，A和B均未在模型中模拟水流混合；B和C均采用24 h前期上游流量和面雨量模拟前期水情。

从图4可见，由于未考虑水流混合模拟，A和B模拟出的流域出口水位存在剧烈波动，与实际水位变化情况并不吻合，这与直接分解N-S方程产生的“汇流水流方向异常”有关；作为本文的径流汇流模拟方法，C模拟的结果无水位剧烈波动，与实际水位曲线最接近。逐小时统计模拟水位与实际水位误差，求取3种方法的平均水位误差及与实际水位曲线的相关系数。从表4可见，基于稳定流场的径流汇流模拟的平均水位误差最小，与实际水位曲线的相关系数达0.86，为高度正相关，其模拟结果具有一定的可信度，能在一定程度上指导暴雨洪涝灾害预测工作。

图4 典型时段校场坝模拟水位与实际水位对比Fig.4 The comparation between simulated water level and actual water level of Jiaochangba

表4 3种径流汇流方法的水位误差统计
Tab.4 The water level error statistics of threerunoff simulation methods

模拟方法平均水位误差/m相关系数A0.310.53B0.230.72C0.180.86

5 结 语

针对前期水情误差影响径流汇流模拟精度的问题，提出了一种基于稳定流场的径流汇流模型。首先探讨了栅格模式下N-S方程的平流项、压力项和外立项的具体实现，建立了初步的径流汇流模型。实验发现，直接分解后得到的模型存在“汇流水流方向异常”问题。通过在径流汇流模拟中添加水流混合环节，不仅解决了该问题，也为稳定流场创建奠定了基础。然后，将径流汇流模拟前的上游流量序列作为上游来水量，通过模型迭代模拟流域前期水情；最后，通过岷江上游镇江关至校场坝流域径流汇流模拟验证，发现本文方法更真实地反映了流域径流汇流过程，模拟精度得到了提升，在流域暴雨洪涝预警预测中具有一定的应用价值。

然而，模型只是原型的模拟，不可能反应原型的一切特征，尤其像流域径流汇流这样复杂的现象。本文将上游入流作为流域前期的唯一水源，并未考虑流域冰川融水、地下水等产汇流因素，仅能得到近似的前期水情。另外，本文选取的实验流域无明显的闸坝调蓄，人类活动影响较小。要进一步提高流域径流汇流模拟精度，需要在模型中进一步考虑冰川融水、人类活动等因素。

□