有效降雨在定义上等于实际降雨扣除降雨损失所剩下的部分,在表现形式上为形成地表径流的那一部分降雨。降雨损失一般包括蒸发蒸腾、植物截留、填洼、土壤入渗等方式,因此有效降雨强度I e 的计算应由流域上的实际降雨强度I以及蒸发蒸腾、植物截留、填洼、土壤入渗等降雨损失所决定。
在黄土高原沟壑区,降雨产流方式以超渗产流为主,即当流域降雨强度超过土壤入渗能力时,便产生径流。一般情况下每次降雨都历时不长,降雨损失也主要为植物截留和土壤入渗两种方式,蒸发蒸腾和填洼损失量一般都很小,可以忽略。
本模型有效降雨计算中主要考虑植物截留和土壤入渗两种降雨损失,故而有效降雨的各组成成分以及相互之间的物理过程见图2-4。
图2-4 有效降雨示意图
植物截留过程与土壤的渗透过程很相似,即当降雨被植物林冠截留后,产生初始截留强度;当植物林冠截留饱和后,仍具有一定的截留能力,称稳定截留强度。它们分别相当于土壤下渗过程中的初始下渗率和稳定下渗率。因此,这里采用Horton的渗透方程来描述植物截留损失的物理过程 [64] 。
式中 J t ——降雨过程中任意时刻的植物截留强度,mm/h;
J c ——植物林冠稳定截留强度,mm/h;
J 0 ——植物林冠初始截留强度,mm/h;
α——植物林冠特性系数,%,取值与降雨强度和植物种类(林分)有关;
t——降雨历时。
研究表明:初始截留强度J 0 与降雨强度I和郁闭度A直接相关,且有J 0 =AI;而稳定截留强度J c 与降雨强度I和林冠特性系数α有关,可由试验观测结果计算得到 [64] 。
本书所基于的产流模式是超渗产流模式,这也是我国黄土高原地区最普遍适用的产流模式。超渗产流意味着只有当有效降雨强度超过土壤入渗能力时,地表才能形成径流。超渗产流模式示意图见图2-5。
图2-5 超渗产流模式示意图
而土壤入渗能力则采用物理概念明晰的Green-Ampt方程 [65] 来计算:
式中 f——土壤下渗能力,m 3 /s;
F——土壤累积入渗量,m 3 ;
K——土壤饱和导水率,m/s;
t——时间,s;
S F ——土壤湿润锋面处土壤水吸力,m;
θ s ——土壤饱和含水率,m 3 /m 3 ;
θ i ——土壤初始含水率,m 3 /m 3 。
土壤入渗过程模拟的关键是确定积水开始时刻t p 。由Green-Ampt公式可知,土壤下渗能力f是随累积入渗量F的增加而减小的,当土壤下渗能力f逐渐下降到等于实际降雨强度I,即f=I时,此时地表开始产生径流,定义这个时刻为积水开始时刻t p ,此时的累积入渗量为F p 。
在降雨强度I恒定的情况下,可直接运用Green-Ampt土壤入渗方程计算。
在积水开始t p 时刻,由于有f=I关系成立,因此由式(2-2)可以求出:
则积水时刻:
因此,下渗过程可表示为:
式中F为积水开始t p 时刻之后的累积入渗量,由于不是从t=0开始积水,根据Mein & Larson的研究 [66] ,F的计算须采用如下修正后的公式,用牛顿迭代法求解:
t s 表示假设由t=0就开始积水到累积入渗量F=F p 时所需的时间,计算如下:
由于Green-Ampt土壤入渗方程只适用于恒定降雨强度的情况,然而在实际情况下,降雨强度一般多为非恒定降雨强度的情况。为了能适用于更为普遍的非恒定降雨强度下土壤入渗能力的计算,本书运用预测校正法对Green-Ampt方程的迭代计算进行改进,以拓展Green-Ampt土壤入渗方程的适用范畴。具体改进方法如下。
将非恒定降雨强度下的降雨过程分解成由若干个降雨强度恒定的Δt时段组成,Δt也就是分布式模型中的时间步长。
(1)在n=1时(即第一个Δt时段内)。由于在Δt时段内,降雨强度I恒定,因此可按恒定降雨强度情况来计算。首先预测在第一个Δt时段内还未开始积水或处于积水开始临界,则应有f=I关系成立,所以:
如果 <t,则说明地表已开始积水,预测不成立,应校正有:
如果t p ≥t,则预测成立,有:
(2)在n时刻(n>1)。同样首先预测在此时段内还未开始积水或处于积水开始临界,则应有f=I关系成立,所以有F(n)=F(n-1)+I(n)Δt。
假如 <t,则说明地表已经积水,预测不成立,应校正有:
式(2-14)中t s 表示假设由t=0就开始地表积水到F=F p 时所需的时间,按式(2-9)计算。
假如 ,则预测成立,有:
式中 F(n)——n时刻的累积入渗量,其他变量含义同上。
需要说明的是,为了保证计算的精度,时间步长Δt不能取得太大(最好不超过100s),不然会因为时间间隔太大而影响计算的精度,并引起较大的误差。
模型中每个网格单元的地表径流计算,采用水量连续平衡方程(Beasley et al.) [67] ,表示为:
式中 W ——单元中所滞留的水量,m 3 ;
t——时间,s;
W i ——进入单元格的水量,m 3 ;
W o ——流出单元格的水量,m 3 。
而在Δt时间内单元格所滞留的水量(以体积计,下同)又可进一步表示为:
式中 A(i,t)和A(i,t-Δt)——分别代表t时刻以及t-Δt时刻垂直于径流方向的过水断面面积,m 2 ;
Δx——网格空间步长(正方形网格),m;
Δt——模型时间步长,s。
进入单元格的水量W i 包括有效降雨量和从相邻单元汇入当前单元的水量之和,可表示为:
式中 I e (i,t)——当前时刻有效降雨强度,m/s;
Q(u,t-Δt)——t-Δt时刻相邻单元汇入当前单元的流量,m 3 /s;
u——相邻8个网格单元中汇入当前单元的那些网格单元。
流出单元格的水量W o 即为进入下一相邻单元格的水量,可表示为:
式中 Q(i,t)——当前时刻流出单元格的流量,m 3 /s。
由于在本书中实行了流域的“沟坡分离”,其中把流域的水流汇集主流路单元定义为沟道,其余单元定义为坡面。所以,地表径流计算过程也相应地分为如下的坡面径流计算和沟道径流计算两部分。
(1)坡面径流计算。在坡面单元的地表径流计算中,由于坡面径流深度较小,并常在整个网格单元上漫流,因此把垂直于径流方向的坡面水流断面面积A概化成以网格空间步长为边长的矩形断面,坡面径流示意图见图2-6。
此时坡面径流深度h可采用下式计算:
式中 A——垂直于径流方向的矩形过水断面面积,m 2 ;
Δx——网格空间步长(正方形网格),m。
坡面单元流速v采用谢才(A.Chezy)公式计算:
式中 n——曼宁糙率系数,根据流域下垫面因子和土地利用类型的不同而选取相应的值,具体参考L.F.Huggins [68][69] 等人的成果;
S 0 ——网格单元地表坡度比降;
h——径流深度,m。
图2-6 坡面径流示意图
(2)沟道径流计算。在沟道单元的地表径流计算中,由于沟道是水流汇集和水沙输移的主要通道,其径流深度较大,并且沟道过水往往只占整个网格单元的一部分,相应的沟道径流示意图见图2-7。
图2-7 沟道径流示意图
由于水流并不在整个网格单元上漫流,因此沟道过水断面面积A并不能像坡面那样概化成分布于整个网格单元边长的矩形断面,而必须根据沟道的水位-流量(即h-Q)关系,运用谢才公式Q= 来确定沟道上各级实测流量Q所对应的过水面积A的值,即A-Q对应曲线,再与水量连续平衡方程式(2-19)~式(2-22)相结合,采用牛顿迭代法求解出A(i,t),并通过谢才公式求得h和v。
模型中整个流域的地表径流计算原理为:程序自动判断当前网格单元是坡面网格单元还是沟道网格单元,并自动调用相应的地表径流计算公式计算出当前网格单元的地表径流,然后通过程序所生成的水流汇集网络图,并运用水量连续平衡方程,将流域中的各个网格单元联系起来进行水流汇集计算,即可确定小流域任意网格单元在时间和空间分布上的产汇流过程。