坡面水流在流域水文学研究与应用领域中占有十分重要的地位。目前基于物理过程分布式参数流域水文模型开发基本上是以坡面水流的数学物理模拟为基础的(Abbott et al., 1986; Ascough Ⅱ et al., 1997)。同样, 揭示森林植被影响地表径流过程的物理学机制, 也只有从坡面水流运动的基本方程出发, 针对森林植被生态系统的特点, 定量研究林地坡面水流运动的基本规律。

采用水动力学方法研究和模拟森林植被地表径流运动, 最为重要的参数是地表径流的阻力系数。同时, 地表径流阻力系数的估算对于评价水源保护林对地表径流土壤侵蚀、林地地表贮水等方面的影响, 也是十分必要的(沈冰等, 1993; 张洪江等, 1995)。由于坡面水流的水层很薄, 水流的阻力与河渠水力学中的糙率不同, 因此, 一般称坡面薄层水流的阻力系数为"有效糙率"(Engman, 1986; 沈冰等, 1994)。但在野外自然坡面状态下, 要想通过实测流速与水深获取坡面地表径流的有效糙率是极为困难的, 因此, 如何根据实测流量过程线通过寻优计算获取地表径流的有效糙率, 具有重要的意义。

1 林地坡面流基本方程 1.1 林地坡面流基本方程组的形式

对形如图 1的坡面薄层水流, 不考虑雨点冲击力对地表径流的影响, 并近似为定床水流运动(即不考虑水流侵蚀地表土壤引起几何形态变化对阻力的影响), 则坡面单宽水流一维运动连续方程可以表述为:

(1)

若忽略旁侧入流在地表流方向的流速, 动量方程可表达为:

(2)

以上两个方程即为经典的圣维南方程组, 其中:h(x, t)为水深, u(x, t)为水层流速, x为沿水流方向的距离, t为时间, w(x, t)为单位面积单位时间净雨率(或称旁侧入流), g为重力加速度, θ为坡度, S_f为摩阻比降。

Morris(1979)的理论研究证明, 不管坡面的坡度是多大, 也不管采用那一种边界条件, 特别是在天然流域较陡的坡面上, 运动波数k不可能太小, 小坡度近似求得的坡面径流过程线的误差均较小。因此, 令sinθ≈tan θ, cosθ≈1, 重写坡面流圣维南方程组为:

(3)

(4)

式中:S₀为坡度, S₀=tan θ。

1.2 林地坡面净雨率

森林植被生态系统从水文功能作用上可分3层:(1)林冠层; (2)地被物层(3)根系层, 因此, 可以从这3个方面探讨林地坡面净雨率。

降雨进入森林植被生态系统首先经冠层进行第1次分配, 冠层对净雨率的影响可由下述指数衰竭方程来刻画:

(5)

式中:ic(x, t)为林冠截流强度, C_s为林冠饱和截流容量, C₀为林冠初始持水量, k₁为林冠截流强度衰减指数, 当考虑的对象为裸露的坡面时, ic(x, t)=0。

由于地被物层对地表径流的影响可以归结为对水流阻力的影响, 其对林地坡面净雨率的影响可采用与上式相似的方程加以描述。森林植被的根系层对净雨率影响是通过土壤水分入渗强度的影响来体现的。地表处的入渗强度f(x, t)可以通过求解土壤水分运动Richards方程或通过一定的土壤入渗方程来获得。据此, 林地坡面水流的净雨率可表示为:

(6)

式中, r(x, t)为降雨强度; if(x, t)为地被物层截流强度。

2 林地坡面水流的求解

坡面流基本方程是非线性偏微分方程, 其解析解尚难求出, 随着计算机和数值计算技术的迅速发展, 目前可以对坡面流动力波模型求得数值解。但是, 由于方程中的某些参数难以由野外实验求出, 因而在实际水文模型开发中动力波模型数值解法的应用相对较少, 但数值模拟结果对于运动波近似和扩散波近似的精度检验具有十分重要的作用(Woolhiser et al., 1967; Morris, 1980; Govindaraju et al., 1988; 1990;1991)。目前, 运动波模型和扩散波模型近似被广泛应用于基于物理过程分布式水文模型的开发中(文康等, 1991; Beven, 1982; 1985;1989;1991)。

从实际流域在天然降雨条件下的水文模拟角度分析, 由于天然降雨及土壤入渗具有时空变异性(雷志栋等, 1988; Skaggs, 1982), 净雨率在求解域上是随时间空间变化的, 这就导致了坡面流的非恒定性, 在这种条件下求解坡面流基本方程就更为复杂。因此, 一般分布式水文模型开发均假定在计算网格内降雨和入渗不随空间变化(Nearing et al., 1989; Dawes, 1997; Baffaut et al., 1997; Liu et al., 1997), 则可用某一入渗强度方程描述土壤水分入渗过程, 将其与降雨过程耦合, 并与基本方程联立求解, 就可获得坡面流的径流过程线。对较小尺度的坡面水流而言, 忽略地表径流水深对入渗强度的影响和入渗的空间变异后, 可以用同一入渗模型来模拟, 这样就可以将入渗强度随时间变化对坡面流的影响考虑在内, 从而使得渗透坡面水流的求解就大为简化。

坡面流采用运动波近似, 方程(4)就可采用坡面流的单宽流量与水深的函数关系来表达:

(7)

式中:α、m为与水流形态有关的参数。

将渗透坡面土壤水分运动简化为入渗率只随时间变化的Smith-P入渗模型; 将坡面流圣维南方程组简化为一维运动波模型, 从而大大地简化了求解过程。综上, 渗透坡面林地地表径流求解模式如下:

式中:t_p为地表开始积水的时间, f ∞为稳定入渗率, t₀为临界时间常数, A, a为系数, C_fs、C_f0、k₂分别为地被物层饱和截流容量、初始持水量和截流强度衰减指数。模型求解的初始条件为干床条件, 即:

(9)

上边界条件为:

(10)

对于以上方程组, 可采用一步Lax-Wendroff有限差分格式计算坡面流求解域每一点的水深, 方程(8a)的差分方程可写为(Kibler et al., 1970; Woolhiser, 1974):

(11)

相应地, 土壤水分入渗方程的差分值可由(8c)、(8d)式求得, 计算时, t_p可由下式确定:

(12)

运动波模型的差分方程由二阶近似, 其稳定性可由下式判断:

(13)

模型求解的误差为二阶差商0(x²)。

至此, 可以通过上述方法, 对坡面水流进行数值模拟计算, 求出坡面水流出口(x=L₀)的流量过程线。

3 林地坡面流有效糙率的计算方法

当林冠截流强度和土壤水分入渗强度模型参数已知时, 求解以上差分方程共有2个未知数α、m。由于坡面水流阻力特征与水流流态有关, 若对于特定的水流流态, m为一定值, 则可选一定的时间步长Δt和距离步长Δx, 通过叠代运算, 求得x=L₀(出口处)的流量过程Q_c(t)=q(t)*B(B为降雨小区宽度), 以方程(14)为目标函数:

(14)

式中:Q₀、Q_c分别为实测流量和计算流量值; d为实测流量过程的点数。当目标函数的F为最小值时, 利用式(13)检验计算的稳定性, 若此式成立, 此时的α即为所求的最优值。对于水流流态不定的情形, 亦可通过这一方法优选参数α、m值。然后根据一定的阻力函数关系就可求得不同林地的地表径流阻力系数(有效糙率系数), 如用出口处的水深通过满宁公式就可求得有效糙率值。计算过程如下:

地表径流水流阻力的谢才公式:

(15)

谢才系数C不随水流雷诺数或水深而变化。式中R为水力半径, 对于坡面薄层水流而言, 可用水流深度h来代替, 对于运动波近似, S₀=S_f, 则上式可写为:

(16)

谢才系数用满宁公式来确定, 即:

(17)

坡面林地的有效糙率n就可由正式计算:

(18)

4 人工模拟降雨条件下不同林地地表径流阻力的计算和分析 4.1 实验标准地与实验方法 4.1.1 林分特征

野外人工模拟降雨试验地点在密云水库上游流域内。试验林分包括刺槐、油松人工林和蒙古栎林3种水源保护林, 每种水源保护林林分调查结果见表 1。

4.1.2 林地土壤水分入渗方程

人工模拟降雨标准地土壤水分入渗Smith-Parlange方程参数是由针头式模拟降雨机在恒定雨强条件下的实测数据拟合的(于志民等, 1999), 其结果如表 2所示。

4.1.3 人工模拟降雨实验方法

人工模拟降雨采用摆喷式中型人工模拟降雨机, 降雨高度为4 m, 试验面积是由长6 m、宽2 m的钢板围成的径流小区, 降雨机具的雨强率定和流量测定方法见于志民等(1999)。

4.2 坡面流运动波近似与参数计算

在认为坡面水流为层流, 即水流流速与水深呈线性关系的条件下, 于志民和王礼先(1999)计算了本研究区不同林分的巴津系数。以林分结构、土壤类型、地表枯落物等相似为依据选取本次实验标准地的巴津系数计算人工模拟降雨坡面流的运动波数, 其值在103数量级, 远远大于坡面水流运动波近似判别标准, 因此, 可以用运动波模型全项近似坡面流(Lighthill et al., 1965)。

由于人工降雨高度较林冠低, 同时计算选用了标准地湿润后的第2次和第4次降雨, 林地地被物层处于近饱和状态, 因此在根据式(8)给出的方程组计算不同实验标准地地表流的有效糙率时, 可以不考虑林冠和地被物的截流损失, 即方程(8a)中的i_c和i_f均为零。计算的时间步长Δt=60s, 距离步长Δx=100 cm。以刺槐人工林P1R4、P2R2、P3R2 3次降雨试验, 油松林人工林P4R3降雨试验、蒙古栎林P5R2降雨试验的实测流量过程, 分别对这5个试验标准地的阻力进行叠代计算, 计算结果如表 3。

从表 3中可以看出, 刺槐林地的地表阻力最大, 有效糙率n在0.1931~0.2131之间, 其次是蒙古栎林地, 为0.1701, 油松林地为0.1423。从参数寻优的结果来看, m值在2.76~3.090之间, 此值接近层流时的m=3。当m=1时, 坡面流为线性水库, 当m=3/2时为100%紊流谢才阻力, 当m=5/2时为100%紊流满宁阻力, 当m=2时为67%紊流谢才阻力(Ponce et al., 1997)。从率定的有效糙率来看, 坡面人工林地地表径流阻力系数在典型低草草原(0.10~0.20)和兰草草皮(0.17~0.48)地表的满宁系数之间(Woolhiser, 1974), 与陈奇伯等人(1994)通过实验研究的长江花岗岩地区松栎混交林、马尾松纯林的地表径流阻力系数较为接近。造成刺槐林地有效糙率较大的原因有3个方面:1是由于当时营造的刺槐林是为了获得薪柴, 造林密度较大, 达5000株/hm²; 2是由于刺槐具有很强的根孽能力, 形成较为密集的刺槐灌丛分布在刺槐林之间; 3是刺槐林枯枝落叶层较厚。相比之下, 蒙古栎林由于枯枝落叶层较厚, 但林下植被较刺槐林稀疏, 其阻力系数较刺槐林小。油松林为常绿针叶林, 林下枯枝落叶层较薄, 一般在0.5 cm左右(于志民等, 1999), 其林地地表径流的阻力系数最小。

为检验参数率定的可靠性, 选取刺槐人工林Plo t1 P1R9, Plot3, P3R3;油松人工林Plot4 P4R4;蒙古栎林Plot5 P5R3 4次人工模拟降雨试验, 采用表 3中的最优参数值对其进行数值模拟计算, 并与实测流量过程进行比较, 将计算结果与实测结果分别点绘在图 2、图 3、图 4和图 5中。

从图中可以看出, 通过上述模型和参数计算的流量与实测流量十分吻合。由于退水过程实测点数较少, 图中没有给出退水段的模拟结果。采用运动波近似刻画陡坡坡面流的运动过程是可行的。

5 结语

通过以上的分析和计算可以看出, 针对林地地表径流运动形成过程和运动特点建立的渗透坡面林地地表径流运动基本方程, 采用数值模拟和叠代计算可以获取地表径流的有效糙率值.对于已知糙率的情形下则可使用本模型和计算方法模拟地表径流运动过程.模型和计算方法可以用于基于物理过程分布式参数流域水文模型开发中.