坡面降雨入渗产流规律的数值模拟研究

摘要：本文采用运动波理论和两次改进后的Green-Ampt入渗模型建立了坡面降雨入渗产流的动力学模型，并得到了实验资料的良好验证。运用该模型研究了简单坡面上降雨入渗产流的动力学规律，分析了雨强、土壤初始含水量、渗透系数、坡面阻力，以及坡长、坡度等因素对坡面产流过程的影响规律，得出了一些有益的结论。

关键词：入渗 产流 坡面 动力学

1 概述

　　雨水降落在坡面上将产生雨水的聚集并形成坡面水流。坡面水流是土壤水蚀过程的主要动力，搞清产流的动力学特点是进一步研究侵蚀过程规律的基础。坡面水流不同于一般明渠流动，其水深极浅(一般只有几毫米)，沿程不断有质量源和动量源加入，使其随时间和空间有较大的变化。且坡面流的坡度较一般河渠陡得多，边界条件也更为复杂。这些特点使得对坡面水流的研究有相当的难度。

　　坡面产流研究已有很长历史，但对它的数学求解还只有三十多年。60年代后期Woolhiser和Ligget(1967)将运动波模型引入坡面水流研究，大大简化了计算工作，促进了研究的发展。运动波模型是从一维圣维南方程简化而来，其基本假设是水流的能坡和底坡相等，并借助Chezy阻力公式得到流量和水深的关系。Woolhiser和Ligget的研究结果表明在运动波波数k&>10时，运动波模型可以很好地描述坡面水流运动。而实际坡面流的运动波波数一般远大于10(沈冰等，1996)。因此，运动波近似是一种较好的数学描述方式。其后，又有对运动波理论的修正(Ponce，1978，Govindaraju，1988)，保留了水深的沿程变化项，相当于压力梯度，被称为扩散波模型。该模型扩展了适用的参数范围，但并无实质性改进，因此实际应用仍以运动波为主。也有使用完整圣维南方程求解实际问题的(戚隆溪，1997)。土壤入渗过程的研究也有很长历史，从1911年提出概念明确形式简单的Green-Ampt积水入渗模型开始，相继有Horton(1940)，Philip(1957)等模型出现，但G-A模型仍以其简单的形式，明晰的物理概念，良好的扩展性和可信的应用效果受到广泛重视，特别是经过Mein && Larson(1973)和Chu(1978)的两次改进，使其可应用于不均匀降雨的入渗计算，更使它成为最有效和应用最广泛的模型。在国内，G-A模型尚未受到重视，Horton模型曾得到相当广泛的运用，但其参数的物理意义明显不如G-A模型明晰。也有研究者使用更基本的土壤水分运动微分方程，但所需的参数更加难于获取，计算也更为复杂。

　　本文工作旨在建立物理概念明晰的降雨入渗产流综合计算模式，并用以研究简单坡面的产流过程，分析各主要因素的影响和各主要因素的影响和各主要参量的变化规律。以期对坡面产流的动力学规律有清楚的认识。

2 计算模式

　　坡面流运动十分复杂，目前主要采用运动波理论、扩散波或完整圣维南方程进行描述。正如前文所述，运动波近似理论在大多数情况下可以很好地描述坡面流运动过程，且计算简单。因此本文仍采用一维运动波理论，即坡面流基本方程为

　(1)

此处第二式直接使用了水力学中熟知的Chezy公式和Manning公式。其中，x为沿坡面向下的坐标，t为时间(s)，h为水深(m)，q为单宽流量(m²/s)，p为降雨强度(m/s)，此处假设降雨方向垂直向下，i为入渗率(m/s)，S₀为坡面坡度，S₀=sinθ，θ为坡面倾角，n为Manning糙率系数。

　　土壤的入渗过程对坡面流的形成和流动过程影响很大，本文采用形式简单、物理概念明晰的G-A入渗模型，其计算方程为

i=dI/dt=K[1+(θ_S-θ_i)S/I]

I=Kt+S(θ_S-θ_i)ln(1+I/S(θ_S-θ_i)(2)

　　其中 K为土壤饱和导水率(渗透系数)(m/s)，θ_S为土壤饱和含水率，即有效孔隙率(%)，θ_i为土壤初始含水率(%)，S为土壤吸力(m)，I为累积入渗量(m)。

　　经典的Green-Ampt模型是干土积水入渗模型，其前提是在整个入渗过程中地表始终有积水。Mein && Larson 1973年将其推广应用至降雨入渗的情况。设有稳定的雨强p，只有p大于土壤的入渗能力时，地表才能形成积水。而在降雨的初始阶段，全部降雨都渗入地下。由G-A模型知，入渗率是随累积入渗量的增加而减小的。设想当累积入渗量达到某一值时，i=p，此时开始积水，称此累积入渗量为I_p。因此由G-A模型入渗公式可以导出开始积水时的I_p值

I_p=(θ_S-θ_i)S/(p/K)-1　(3)

　　开始积水时间由t_p=I_p/p给出。因此整个过程的入渗率可表示为

式中的I为积水开始后的累积入渗量(包含未积水时段的入渗量在内)。由于不是由t=0开始积水，I的计算须采用修正后的公式

K[t-(t_p-t_S)]=I-S(θ_S-θ_i)ln[1+I/S(θ_S-θ_i)] t&>t_{p　　　　　　(5)}

　　t_S表示假设由t=0开始积水，到入渗量I=I_p(或i=p)时所需时间，可理解为一个虚拟时间，可计算如下

Kt_S=I_p-S(θ_S-θ_i)ln[1+I_p/S(θ_S-θ_i)]　　　　(6)

　　改进的主要思想是将整个过程假设为从一开始就是积水入渗，这样该曲线在积水后部分相对于实际入渗曲线将向左平移t_p-t_S，将这条曲线向右平移t_p-t_S，再加上积水前的入渗强度等于降雨强度的关系，就得到真实的入渗过程。

　　但稳定降雨在实际应用中远不能满足需要，Chu(1978)将Mein && Larson(1973)改进的G-A模型再作推广，使其可应用于变化的降雨过程。基本作法是，对每个计算时段将地表状态分为四种情况：

　　1.开始无积水，结束无积水

　　2.开始无积水，结束有积水

　　3.开始有积水，结束有积水

　　4.开始有积水，结束无积水

　　在每一时段开始，已知降雨总量与入渗总量，剩余总量。根据两个因子判断时段结束时是否有积水。

　　若时段开始无积水，使用因子c_u，若时段开始有积水，使用因子c_p，其表达式为

其中M代表θ_S-θ_i，P(t_n)代表t_n时刻降雨总量，R(t_n-1)代表t_n-1时刻剩余总量。可以证明，时段结束时积水与否与此两因子的正负等价。当i&<K时，始终无积水，不用此两因子判断。

3 模型求解

　　模型的求解也分为两部分。运动波方程是一个非线性的对流型方程，求解采用了简单的一阶显式迎风格式。经过比较，一阶显式迎风格式对此问题的计算结果在幅值、相位、守恒性几个方面的综合效果较其他一些格式为优。我们曾使用过Preissmann四点偏心格式，它在计算此类坡面薄层水流运动波方程时耗散过大，表现为达到平衡产流的时间偏快。我们也曾试过中心差分格式和高阶迎风格式，结果发现标准的中心差分格式仍然和在普通的线性对流型方程中的表现一样是绝对不稳定的。一些变形的中心差分格式能得到较好的结果，但在守恒性和稳定性上不如一阶迎风格式。高阶迎风格式的结果也与此类似，其中二阶迎风格式实际计算时稳定的Currant数比一阶迎风小一半。除此之外，一阶迎风格式的计算是最为简单的，边界条件也易于处理。

　　将式(1)第一式写为

其中F代表通量，F=F(h)=q=1/nh^5/3S₀^1/2，Sr为源项，Sr=pcosθ-i。

　　计算格式形式按有限体积法写为

　　一阶迎风格式的通量写为

　　[t_n，t_n+1]时段入渗方程的求解首先根据两个因子判断t_n+1时刻有无积水，若无积水可直接由R(t_n+1)=R(t_n)得到Sr=0。若有积水，按照方程(5)给出的关系直接用Newton法求解代数超越方程。若是均匀降雨、坡度不变且土壤的物理参数也不变，对所有的网格可以采用同样的积水时间，可以直接根据(5)式求解。若由于各种参数的变化造成各网格的产流时间不同，则还须考虑产流网格向相邻未产流网格汇入造成该网格的既定入渗曲线的改变。这时c_u，c_p的公式中还须加入相应的汇入量。

4 模型验证

　　Morgali and Linsley(1965)曾在长为72 feet(合22m)，坡度S₀=0.04的坡面上进行了无入渗降雨产流实验，降雨强度为3.66 inches/hour(1.55mm/min)。实验中分别采用了光滑和粗糙两种表面条件。运用本文模型对无渗透坡面上的降雨产流过程进行模拟计算。几种条件下单宽流量随时间变化过程的数值模拟结果与实验结果的对比如图1，图2，图3所示。所取阻力系数值均按原作者进行数值计算时所取值。结果表明，当阻力系数取值较准确时，数值解与实测数据符合较好。需要说明的是，实验中粗糙表面所设粗糙颗粒远大于在一般坡面上的尺度，因此其阻力系数也较一般无植被覆盖的平坦坡面大得多。退水过程中由于已无雨滴打击的影响，其阻力系数有所减小。

　　包括渗透过程的降雨产流过程计算我们采用了Lima(1992)的实验数据来验证。该实验在长1m，宽0.5m，坡度S₀=0.1的土质坡面上进行。降雨强度为0.03741mm/s。本文计算中涉及到的实测土壤参数和根据实测土壤参数与降雨总量率定的参数为：K=1.67×10^-6m/s，θ_S=0.506，θ_i=0.0107，S=0.02m。计算结果和实验结果见图4。整个产流退水过程均符合较好。与原作者用土壤水分运动微分方程求得的结果也很接近。

5 坡面降雨产流规律

上述模拟验证结果表明本文所建立的模型能够较好地模拟坡面降雨入渗产流的水动力过程。因此，我们进一步运用该模型进行数值模拟实验，分析各种因素对产流过程的影响。

5.1 降雨强度、土壤入渗率、初始含水率对产流过程的影响

　　通过改变Lima实验中初始给定的参数值，分别对Lima实验在不同降雨强度，不同初始含水率(不同饱和度)和不同土壤渗透系数下的产流过程进行数值预测，结果分别如图5、图6和图7所示。图例中的系数分别代表与原始实验参数的比值(后文中图例意义相同)。限于篇幅，本文仅给出了单宽流量的变化过程，而未一一给出相应各种情况的流速、水深和切应力(τ=γhS₀)的变化曲线。除特别说明，它们的变化规律与单宽流量的规律相同。　　

　　各模拟结果的一个共同规律是，产流过程都有一个快速增长的初始阶段，然后增长速度迅速降低，逐渐接近平衡产流状态。达到平衡需要很长时间。而降雨停止后的退水过程则又变得很快。

图5 不同降雨强度产流过程
Fig.5 Runoff process with various rain intensities

　　图5的计算结果表明，降雨强度增大，坡面流单宽流量随之增大；坡面流的其它各水动力因子(流速、水深和切应力)均随之增大。而产流开始的时间和产流的初始阶段随雨强的增加逐渐缩短。退水过程也随降雨强度的增加而略有延长。平衡时的流量几乎和雨强成正比。根据模型，水深和流速与流量分别有3/5和2/5次幂的关系，因此它们与雨强也有近似3/5和2/5次幂的关系。切应力则与水深有同样的规律。

　　从图6结果可以看到，土壤初始含水率越高，或者说土壤初始非饱和度越低，坡面上的产流量越大，各动力学参量(流速、水深和切应力)也相应越大，且产流开始时间和达到平衡的时间也有所提前。

　　图7所示的模拟计算结果表明，随着土壤渗透系数的减小，产流过程的变化规律基本类似于其随土壤含水率的变化规律。很显然，土壤初始非饱和度低或渗透系数蝎使土壤在降雨产流的初始阶段吸收的水量少，因此产流量较大。但在接近平衡时，土壤的初始非饱和度影响不大；渗透系数因为是饱和状态的入渗能力值，则会影响到平衡产流量。

　　另外，从图1，2，3的计算结果还可以发现坡面阻力系数对产流过程也会产生一定影响。一般随着阻力系数的增加，产流达到平衡的时间和退水时间均会延长，即起到一种延迟作用。而坡面阻力系数对于产流的开始时间，产流量，以及接近产流平衡时其他各动力学参量则没有太大影响。

5.2 坡面长、坡度对产流的影响

　　图8为不同坡长条件下的产流过程计算结果。表明随着坡长的增加，出口处的产流量随之增大，出口处坡面流的其他各水动力因素也均增大。这与降雨强度增加有类似的作用，但对产流开始时间没有影响。

　　坡度变化对产流的影响情况则比较复杂。图9为不同坡度下包含入渗过程的产流过程模拟计算结果。图10、图11中还给出了出口处水深、流速和切应力随坡度的变化。

　　坡度对产流的影响是多方面的。图9中给出了5°～60°范围内的8条产流过程曲线。表明随着坡度的增大，产流有加快和减少的两种趋势。结合图10可以看到，最大单宽流量(产流停止时刻)仅在坡度很缓时有极小的增加，坡度大于5°后均呈减小趋势。这是因为在雨水垂直于水平面降落的假设下，坡度增大使得实际承雨面积减小，或者可以理解为在相同的坡长上有效降雨强度随坡度增加而减小。由于土壤的渗透系数也会随坡度变化，我们的模型中采用了蒋定生(1998)关于渗透系数随坡度增加而减小的关系，这会使得产流量有随坡度增加的趋势。但坡度大于5°后渗透减少的作用不如承雨面积减少的作用大，其综合的结果主要仍是使产流量减少。其次从图10可知，坡度增大水深是减小的。主要的原因是产流量的减少。由运动波模型可知，在相同流量下坡度增加将使水深减小。因此，坡度增加的总体作用是使得水深减小。

　　图11则表明随着坡度的增加，坡面流的流速和切应力呈先增后减的变化趋势，其间有一个达到最大值的临界坡度。图中给出了从1°～60°范围内的出口处流速和水深的最大值(降雨停止时刻)的计算结果。两者的临界坡度值均约在40°～50°之间。由于流量和水深随坡度增加同时减小，而在不同阶段二者变化的速率不一样，流量减小的速度逐渐增加，是上凸的曲线，水深的变化则正相反，因此流速具有这种变化趋势。切应力的计算直接用了公式τ=γhS₀，因此是水深h和坡度sinθ值的变化规律共同决定了切应力的变化趋势。

　　上述计算是雨强相当大(2.245mm/min，属特大暴雨)，或雨强与渗透系数的比值很大(p/K=22.45)的情况。若保持渗透系数不变，将雨强减小至1mm/min和0.5mm/min，所得结果与上述规律仍然一致，u，τ的临界坡度范围也无明显变化，仅单宽流量略有增加的范围有所扩大，在p=0.5mm/min时为15°以下，增幅仍然很小。

　　坡度在产流中的上述作用将有助于理解土壤侵蚀现象中侵蚀量的临界坡度问题。

6 结论

　　本文建立了坡面上降雨入渗产流的动力学模型，用实验资料验证了模型的有效性。并通过模拟计算，分析了简单坡面上降雨入渗产流的动力学规律。得到以下几点结论。

　　(1)仅降雨强度增大，所有参量(单宽流量，水深，流速，切应力)均会增大。单宽流量与降雨强度基本呈正比关系。水深和流速与雨强分别有近似3/5和2/5次幂的关系，切应力与水深有同样的规律。

　　(2)仅渗透率增大，所有参量减小。

　　(3)仅坡长增加，所有参量均增大。单宽流量与坡长也基本成正比。

　　(4)仅土壤初始非饱和度(饱和含水量和前期含水量之差)增加，所有参量均减小。

　　(5)坡度增加，将会带来几个方面变化。坡长相同时承雨面积将减小，所接受的雨量将减少；另一方面，坡度增大，渗透率将减小，导致产流增加。而在同样的产流量下，坡度增大水流的水深减小流速将增大。综合起来，坡长不变时流量随坡度增大主要是减小的，仅在坡度较缓时有很少增加，水深随坡度增大也减小，但流量变化是凸曲线，水深变化是凹曲线。流速和切应力将先增加，在达到某一临界坡度时变为减小，两者的临界坡度并不一样，但都约在40°～50°之间。

参 考 文 献

[1] J.L.M.P.de Lima，Model KINNIF for overland flow on pervious surface.In Overland Flow A.J.Parsons && A.D.Abrahams(eds)69-88 UCL Press 1992.

[2] Govindaraju，R.S.，On the diffusion wave model for overland flowWater Resour.Res.vol.24 NO.5 734-754 1998.

[3] Horton，R.E.An approach toward a phisical interpretation of infil tration-capacity.Soil.Sci.Soc.Am.Proc.399-417 Vol.5 1940.

[4] Philip，J.R.，The theory of infiltration.Soil.Sci.84，254-264，1957.

[5] Ponce，V.M.，R.－M.Li and D.B.Simons， Applicability of Kinematic and diffusion models，J.Hydraul. p. Amer.SocCivil. Eng. Vol. 104，353－360，1978.

[6] R.G.，Mein and C.L.，Larson，Modeling infiltration during a steady rain.Water Resour.Res.vol.9 NO.2 384-394 April 1973.

[7] SHU TUNG CHU，Infiltration during an unsteady rain，Water Resour。Res.vol.14 NO.3 461-466 June 1978.

[8] Woolhiser，D.A.， and Ligget，J.A.，Unsteady，one dimensional flow over a plane-The rising hydrograph，Water Resour.Res.，vol.3，NO.3，753-771，1967.

[9] 陈国祥，谢树楠，汤立群。1996 黄土高原地区流域侵蚀产沙模型研究。黄土高原水土保持。(孟庆枚主编).黄河水利出版社，1996.

[10] 蒋定生。黄土高原水土流失与治理模式。中国水利水电出版社，1997，9.

[11] 戚隆溪，黄兴法。坡面降雨径流和土壤侵蚀的数值模拟。力学学报，vol.29.NO.3，1997.

[12] 沈冰。地表水文有限元模拟。西北工业大学出版社，1996.