坡面降雨入滲產(chǎn)流規(guī)律的數(shù)值模擬研究

1、坡面降雨入滲產(chǎn)流規(guī)律的數(shù)值模擬研究摘要：本文采用運(yùn)動(dòng)波理論和兩次改良后的Green-Apt入滲模型建立了坡面降雨入滲產(chǎn)流的動(dòng)力學(xué)模型，并得到了實(shí)驗(yàn)資料的良好驗(yàn)證。運(yùn)用該模型研究了簡單坡面上降雨入滲產(chǎn)流的動(dòng)力學(xué)規(guī)律，分析了雨強(qiáng)、土壤初始含水量、浸透系數(shù)、坡面阻力，以及坡長、坡度等因素對坡面產(chǎn)流過程的影響規(guī)律，得出了一些有益的結(jié)論。關(guān)鍵詞：入滲產(chǎn)流坡面動(dòng)力學(xué)1概述雨水降落在坡面上將產(chǎn)生雨水的聚集并形成坡面水流。坡面水流是土壤水蝕過程的主要?jiǎng)恿?，搞清產(chǎn)流的動(dòng)力學(xué)特點(diǎn)是進(jìn)一步研究侵蝕過程規(guī)律的基矗坡面水流不同于一般明渠流動(dòng)，其水深極淺(一般只有幾毫米)，沿程不斷有質(zhì)量源和動(dòng)量源參加，使其隨時(shí)間和空間有

2、較大的變化。且坡面流的坡度較一般河渠陡得多，邊界條件也更為復(fù)雜。這些特點(diǎn)使得對坡面水流的研究有相當(dāng)?shù)碾y度。坡面產(chǎn)流研究已有很長歷史，但對它的數(shù)學(xué)求解還只有三十多年。60年代后期lhiser和Ligget(1967)將運(yùn)動(dòng)波模型引入坡面水流研究，大大簡化了計(jì)算工作，促進(jìn)了研究的開展。運(yùn)動(dòng)波模型是從一維圣維南方程簡化而來，其根本假設(shè)是水流的能坡和底坡相等，并借助hezy阻力公式得到流量和水深的關(guān)系。lhiser和Ligget的研究結(jié)果說明在運(yùn)動(dòng)波波數(shù)k10時(shí)，運(yùn)動(dòng)波模型可以很好地描繪坡面水流運(yùn)動(dòng)。而實(shí)際坡面流的運(yùn)動(dòng)波波數(shù)一般遠(yuǎn)大于10(沈冰等，1996)。因此，運(yùn)動(dòng)波近似是一種較好的數(shù)學(xué)描繪方式。

3、其后，又有對運(yùn)動(dòng)波理論的修正(Pne,1978,Gvindaraju,1988)，保存了水深的沿程變化項(xiàng)，相當(dāng)于壓力梯度，被稱為擴(kuò)散波模型。該模型擴(kuò)展了適用的參數(shù)范圍，但并無本質(zhì)性改良，因此實(shí)際應(yīng)用仍以運(yùn)動(dòng)波為主。也有使用完好圣維南方程求解實(shí)際問題的(戚隆溪，1997)。土壤入滲過程的研究也有很長歷史，從1911年提出概念明確形式簡單的Green-Apt積水入滲模型開場，相繼有Hrtn(1940),Philip(1957)等模型出現(xiàn)，但G-A模型仍以其簡單的形式，明晰的物理概念，良好的擴(kuò)展性和可信的應(yīng)用效果受到廣泛重視，特別是經(jīng)過EinLarsn(1973)和hu(1978)的兩次改良，使其可

4、應(yīng)用于不均勻降雨的入滲計(jì)算，更使它成為最有效和應(yīng)用最廣泛的模型。在國內(nèi)，G-A模型尚未受到重視，Hrtn模型曾得到相當(dāng)廣泛的運(yùn)用，但其參數(shù)的物理意義明顯不如G-A模型明晰。也有研究者使用更根本的土壤水分運(yùn)動(dòng)微分方程，但所需的參數(shù)更加難于獲取，計(jì)算也更為復(fù)雜。本文工作旨在建立物理概念明晰的降雨入滲產(chǎn)流綜合計(jì)算形式，并用以研究簡單坡面的產(chǎn)流過程，分析各主要因素的影響和各主要因素的影響和各主要參量的變化規(guī)律。以期對坡面產(chǎn)流的動(dòng)力學(xué)規(guī)律有清楚的認(rèn)識。2計(jì)算形式坡面流運(yùn)動(dòng)非常復(fù)雜，目前主要采用運(yùn)動(dòng)波理論、擴(kuò)散波或完好圣維南方程進(jìn)展描繪。正如前文所述，運(yùn)動(dòng)波近似理論在大多數(shù)情況下可以很好地描繪坡面流運(yùn)動(dòng)過

5、程，且計(jì)算簡單。因此本文仍采用一維運(yùn)動(dòng)波理論，即坡面流根本方程為(1)此處第二式直接使用了水力學(xué)中熟知的hezy公式和anning公式。其中，x為沿坡面向下的坐標(biāo)，t為時(shí)間(s)，h為水深(),q為單寬流量(2/s)，p為降雨強(qiáng)度(/s)，此處假設(shè)降雨方向垂直向下，i為入滲率(/s)，S0為坡面坡度，S0=sin，為坡面傾角，n為anning糙率系數(shù)。土壤的入滲過程對坡面流的形成和流動(dòng)過程影響很大，本文采用形式簡單、物理概念明晰的G-A入滲模型，其計(jì)算方程為i=dI/dt=K1+(S-i)S/II=Kt+S(S-i)ln(1+I/S(S-i)(2)其中K為土壤飽和導(dǎo)水率(浸透系數(shù))(/s)，S

6、為土壤飽和含水率，即有效孔隙率(%)，i為土壤初始含水率(%)，S為土壤吸力(),I為累積入滲量()。經(jīng)典的Green-Apt模型是干土積水入滲模型，其前提是在整個(gè)入滲過程中地表始終有積水。einLarsn1973年將其推廣應(yīng)用至降雨入滲的情況。設(shè)有穩(wěn)定的雨強(qiáng)p，只有p大于土壤的入滲才能時(shí)，地表才能形成積水。而在降雨的初始階段，全部降雨都滲入地下。由G-A模型知，入滲率是隨累積入滲量的增加而減小的。設(shè)想當(dāng)累積入滲量到達(dá)某一值時(shí)，i=p,此時(shí)開場積水，稱此累積入滲量為Ip。因此由G-A模型入滲公式可以導(dǎo)出開場積水時(shí)的Ip值Ip=(S-i)S/(p/K)-1(3)開場積水時(shí)間由tp=Ip/p給出。

7、因此整個(gè)過程的入滲率可表示為i=pttpi=K1+(S-i)S/Ittp(4)式中的I為積水開場后的累積入滲量(包含未積水時(shí)段的入滲量在內(nèi))。由于不是由t=0開場積水，I的計(jì)算須采用修正后的公式Kt-(tp-tS)=I-S(S-i)ln1+I/S(S-i)ttp(5)tS表示假設(shè)由t=0開場積水，到入滲量I=Ip(或i=p)時(shí)所需時(shí)間，可理解為一個(gè)虛擬時(shí)間，可計(jì)算如下KtS=Ip-S(S-i)ln1+Ip/S(S-i)(6)改良的主要思想是將整個(gè)過程假設(shè)為從一開場就是積水入滲，這樣該曲線在積水后局部相對于實(shí)際入滲曲線將向左平移tp-tS，將這條曲線向右平移tp-tS，再加上積水前的入滲強(qiáng)度等于

8、降雨強(qiáng)度的關(guān)系，就得到真實(shí)的入滲過程。但穩(wěn)定降雨在實(shí)際應(yīng)用中遠(yuǎn)不能滿足需要，hu(1978)將einLarsn(1973)改良的G-A模型再作推廣，使其可應(yīng)用于變化的降雨過程。根本作法是，對每個(gè)計(jì)算時(shí)段將地表狀態(tài)分為四種情況：1.開場無積水，完畢無積水2.開場無積水，完畢有積水3.開場有積水，完畢有積水4.開場有積水，完畢無積水在每一時(shí)段開場，降雨總量與入滲總量，剩余總量。根據(jù)兩個(gè)因子判斷時(shí)段完畢時(shí)是否有積水。假設(shè)時(shí)段開場無積水，使用因子u,假設(shè)時(shí)段開場有積水，使用因子p，其表達(dá)式為u=P(tn)-R(tn-1)-KS/(i-K)u0時(shí)段完畢將積水，u0仍無積水p=P(tn)-I(tn)-R(

9、tn-1)p0時(shí)段完畢仍有積水，u0積水消失(7)其中代表S-i,P(tn)代表tn時(shí)刻降雨總量，R(tn-1)代表tn-1時(shí)刻剩余總量?？梢宰C明，時(shí)段完畢時(shí)積水與否與此兩因子的正負(fù)等價(jià)。當(dāng)iK時(shí)，始終無積水，不用此兩因子判斷。3模型求解模型的求解也分為兩局部。運(yùn)動(dòng)波方程是一個(gè)非線性的對流型方程，求解采用了簡單的一階顯式迎風(fēng)格式。經(jīng)過比擬，一階顯式迎風(fēng)格式對此問題的計(jì)算結(jié)果在幅值、相位、守恒性幾個(gè)方面的綜合效果較其他一些格式為優(yōu)。我們曾使用過PrEissann四點(diǎn)偏心格式，它在計(jì)算此類坡面薄層水流運(yùn)動(dòng)波方程時(shí)耗散過大，表現(xiàn)為到達(dá)平衡產(chǎn)流的時(shí)間偏快。我們也曾試過中心差分格式和高階迎風(fēng)格式，結(jié)果發(fā)

10、現(xiàn)標(biāo)準(zhǔn)的中心差分格式仍然和在普通的線性對流型方程中的表現(xiàn)一樣是絕對不穩(wěn)定的。一些變形的中心差分格式能得到較好的結(jié)果，但在守恒性和穩(wěn)定性上不如一階迎風(fēng)格式。高階迎風(fēng)格式的結(jié)果也與此類似，其中二階迎風(fēng)格式實(shí)際計(jì)算時(shí)穩(wěn)定的urrant數(shù)比一階迎風(fēng)小一半。除此之外，一階迎風(fēng)格式的計(jì)算是最為簡單的，邊界條件也易于處理。將式(1)第一式寫為(8)其中F代表通量，F(xiàn)=F(h)=q=1/nh5/3S01/2,Sr為源項(xiàng)，Sr=ps-i。計(jì)算格式形式按有限體積法寫為hin+1-hni/t+Fi+1/2-Fi-1/2/x=Srni(9)一階迎風(fēng)格式的通量寫為Fi+1/2(h)=1/hhi5/3S01/2(10)t

11、n,tn+1時(shí)段入滲方程的求解首先根據(jù)兩個(gè)因子判斷tn+1時(shí)刻有無積水，假設(shè)無積水可直接由R(tn+1)=R(tn)得到Sr=0。假設(shè)有積水，按照方程(5)給出的關(guān)系直接用Netn法求解代數(shù)超越方程。假設(shè)是均勻降雨、坡度不變且土壤的物理參數(shù)也不變，對所有的網(wǎng)格可以采用同樣的積水時(shí)間，可以直接根據(jù)(5)式求解。假設(shè)由于各種參數(shù)的變化造成各網(wǎng)格的產(chǎn)流時(shí)間不同，那么還須考慮產(chǎn)流網(wǎng)格向相鄰未產(chǎn)流網(wǎng)格匯入造成該網(wǎng)格的既定入滲曲線的改變。這時(shí)u,p的公式中還須參加相應(yīng)的匯入量。4模型驗(yàn)證rgaliandLinsley(1965)曾在長為72feet(合22)，坡度S0=0.04的坡面上進(jìn)展了無入滲降雨產(chǎn)流

12、實(shí)驗(yàn)，降雨強(qiáng)度為3.66inhes/hur(1.55/in)。實(shí)驗(yàn)中分別采用了光滑和粗糙兩種外表?xiàng)l件。運(yùn)用本文模型對無浸透坡面上的降雨產(chǎn)流過程進(jìn)展模擬計(jì)算。幾種條件下單寬流量隨時(shí)間變化過程的數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果的比照方圖1，圖2，圖3所示。所取阻力系數(shù)值均按原作者進(jìn)展數(shù)值計(jì)算時(shí)所取值。結(jié)果說明，當(dāng)阻力系數(shù)取值較準(zhǔn)確時(shí)，數(shù)值解與實(shí)測數(shù)據(jù)符合較好。需要說明的是，實(shí)驗(yàn)中粗糙外表所設(shè)粗糙顆粒遠(yuǎn)大于在一般坡面上的尺度，因此其阻力系數(shù)也較一般無植被覆蓋的平坦坡面大得多。退水過程中由于已無雨滴打擊的影響，其阻力系數(shù)有所減校包括浸透過程的降雨產(chǎn)流過程計(jì)算我們采用了Lia(1992)的實(shí)驗(yàn)數(shù)據(jù)來驗(yàn)證。該實(shí)驗(yàn)在

13、長1,寬0.5，坡度S0=0.1的土質(zhì)坡面上進(jìn)展。降雨強(qiáng)度為0.03741/s。本文計(jì)算中涉及到的實(shí)測土壤參數(shù)和根據(jù)實(shí)測土壤參數(shù)與降雨總量率定的參數(shù)為：K=1.6710-6/s,S=0.506，i=0.0107，S=0.02。計(jì)算結(jié)果和實(shí)驗(yàn)結(jié)果見圖4。整個(gè)產(chǎn)流退水過程均符合較好。與原作者用土壤水分運(yùn)動(dòng)微分方程求得的結(jié)果也很接近。圖1浸透光滑坡面漲水過程Fig.1Rising-stageurvefrunffnnnpereablesthslpe圖2無浸透粗糙坡面漲水過程Fig.2Rising-stageurvefrunffnnnpereablearseslpe圖3無浸透粗糙坡面退水過程Fig.3F

14、alling-stageurvefrunffnnnpereablesthslpe圖4有浸透坡面產(chǎn)流過程Fig.4Runfpressnpereableslpe5坡面降雨產(chǎn)流規(guī)律上述模擬驗(yàn)證結(jié)果說明本文所建立的模型可以較好地模擬坡面降雨入滲產(chǎn)流的水動(dòng)力過程。因此，我們進(jìn)一步運(yùn)用該模型進(jìn)展數(shù)值模擬實(shí)驗(yàn)，分析各種因素對產(chǎn)流過程的影響。5.1降雨強(qiáng)度、土壤入滲率、初始含水率對產(chǎn)流過程的影響通過改變Lia實(shí)驗(yàn)中初始給定的參數(shù)值，分別對Lia實(shí)驗(yàn)在不同降雨強(qiáng)度，不同初始含水率(不同飽和度)和不同土壤浸透系數(shù)下的產(chǎn)流過程進(jìn)展數(shù)值預(yù)測，結(jié)果分別如圖5、圖6和圖7所示。圖例中的系數(shù)分別代表與原始實(shí)驗(yàn)參數(shù)的比值(后

15、文中圖例意義一樣)。限于篇幅，本文僅給出了單寬流量的變化過程，而未一一給出相應(yīng)各種情況的流速、水深和切應(yīng)力(=hS0)的變化曲線。除特別說明，它們的變化規(guī)律與單寬流量的規(guī)律一樣。各模擬結(jié)果的一個(gè)共同規(guī)律是，產(chǎn)流過程都有一個(gè)快速增長的初始階段，然后增長速度迅速降低，逐漸接近平衡產(chǎn)流狀態(tài)。到達(dá)平衡需要很長時(shí)間。而降雨停頓后的退水過程那么又變得很快。圖5不同降雨強(qiáng)度產(chǎn)流過程Fig.5Runffpressithvariusrainintensities圖5的計(jì)算結(jié)果說明，降雨強(qiáng)度增大，坡面流單寬流量隨之增大；坡面流的其它各水動(dòng)力因子(流速、水深和切應(yīng)力)均隨之增大。而產(chǎn)流開場的時(shí)間和產(chǎn)流的初始階段隨雨

16、強(qiáng)的增加逐漸縮短。退水過程也隨降雨強(qiáng)度的增加而略有延長。平衡時(shí)的流量幾乎和雨強(qiáng)成正比。根據(jù)模型，水深和流速與流量分別有3/5和2/5次冪的關(guān)系，因此它們與雨強(qiáng)也有近似3/5和2/5次冪的關(guān)系。切應(yīng)力那么與水深有同樣的規(guī)律。圖6不同初始含水率產(chǎn)流過程Fig.6Runffpressithvariusinitialisturentents圖7不同浸透系數(shù)產(chǎn)流過程Fig.7Runffpressithvariusinfiltratineffiients從圖6結(jié)果可以看到，土壤初始含水率越高，或者說土壤初始非飽和度越低，坡面上的產(chǎn)流量越大，各動(dòng)力學(xué)參量(流速、水深和切應(yīng)力)也相應(yīng)越大，且產(chǎn)流開場時(shí)間和到達(dá)

17、平衡的時(shí)間也有所提早。圖7所示的模擬計(jì)算結(jié)果說明，隨著土壤浸透系數(shù)的減小，產(chǎn)流過程的變化規(guī)律根本類似于其隨土壤含水率的變化規(guī)律。很顯然，土壤初始非飽和度低或浸透系數(shù)蝎使土壤在降雨產(chǎn)流的初始階段吸收的水量少，因此產(chǎn)流量較大。但在接近平衡時(shí)，土壤的初始非飽和度影響不大；浸透系數(shù)因?yàn)槭秋柡蜖顟B(tài)的入滲才能值，那么會(huì)影響到平衡產(chǎn)流量。另外，從圖1，2，3的計(jì)算結(jié)果還可以發(fā)現(xiàn)坡面阻力系數(shù)對產(chǎn)流過程也會(huì)產(chǎn)生一定影響。一般隨著阻力系數(shù)的增加，產(chǎn)流到達(dá)平衡的時(shí)間和退水時(shí)間均會(huì)延長，即起到一種延遲作用。而坡面阻力系數(shù)對于產(chǎn)流的開場時(shí)間，產(chǎn)流量，以及接近產(chǎn)流平衡時(shí)其他各動(dòng)力學(xué)參量那么沒有太大影響。5.2坡面長、坡度

18、對產(chǎn)流的影響圖8為不同坡長條件下的產(chǎn)流過程計(jì)算結(jié)果。說明隨著坡長的增加，出口處的產(chǎn)流量隨之增大，出口處坡面流的其他各水動(dòng)力因素也均增大。這與降雨強(qiáng)度增加有類似的作用，但對產(chǎn)流開場時(shí)間沒有影響。坡度變化對產(chǎn)流的影響情況那么比擬復(fù)雜。圖9為不同坡度下包含入滲過程的產(chǎn)流過程模擬計(jì)算結(jié)果。圖10、圖11中還給出了出口處水深、流速和切應(yīng)力隨坡度的變化。坡度對產(chǎn)流的影響是多方面的。圖9中給出了560范圍內(nèi)的8條產(chǎn)流過程曲線。說明隨著坡度的增大，產(chǎn)流有加快和減少的兩種趨勢。結(jié)合圖10可以看到，最大單寬流量(產(chǎn)流停頓時(shí)刻)僅在坡度很緩時(shí)有極小的增加，坡度大于5后均呈減小趨勢。這是因?yàn)樵谟晁怪庇诔潭让娼德涞募?/p>

19、設(shè)下，坡度增大使得實(shí)際承雨面積減小，或者可以理解為在一樣的坡長上有效降雨強(qiáng)度隨坡度增加而減校由于土壤的浸透系數(shù)也會(huì)隨坡度變化，我們的模型中采用了蔣定生(1998)關(guān)于浸透系數(shù)隨坡度增加而減小的關(guān)系，這會(huì)使得產(chǎn)流量有隨坡度增加的趨勢。但坡度大于5后浸透減少的作用不如承雨面積減少的作用大，其綜合的結(jié)果主要仍是使產(chǎn)流量減少。其次從圖10可知，坡度增大水深是減小的。主要的原因是產(chǎn)流量的減少。由運(yùn)動(dòng)波模型可知，在一樣流量下坡度增加將使水深減校因此，坡度增加的總體作用是使得水深減校圖11那么說明隨著坡度的增加，坡面流的流速和切應(yīng)力呈先增后減的變化趨勢，其間有一個(gè)到達(dá)最大值的臨界坡度。圖中給出了從160范圍

20、內(nèi)的出口處流速和水深的最大值(降雨停頓時(shí)刻)的計(jì)算結(jié)果。兩者的臨界坡度值均約在4050之間。由于流量和水深隨坡度增加同時(shí)減小，而在不同階段二者變化的速率不一樣，流量減小的速度逐漸增加，是上凸的曲線，水深的變化那么正相反，因此流速具有這種變化趨勢。切應(yīng)力的計(jì)算直接用了公式=hS0，因此是水深h和坡度sin值的變化規(guī)律共同決定了切應(yīng)力的變化趨勢。圖8不同坡長的產(chǎn)流過程Fig.8Runffpressithvariusslpelength圖9不同坡度的產(chǎn)流過程Fig.9Runffpressithvariusslpeangles圖11速度、切應(yīng)力隨坡度的變化Fig.11Variatinfvelityan

21、dshearstressattheutletithslpeangle上述計(jì)算是雨強(qiáng)相當(dāng)大(2.245/in,屬特大暴雨)，或雨強(qiáng)與浸透系數(shù)的比值很大(p/K=22.45)的情況。假設(shè)保持浸透系數(shù)不變，將雨強(qiáng)減小至1/in和0.5/in，所得結(jié)果與上述規(guī)律仍然一致，u,的臨界坡度范圍也無明顯變化，僅單寬流量略有增加的范圍有所擴(kuò)大，在p=0.5/in時(shí)為15以下，增幅仍然很校坡度在產(chǎn)流中的上述作用將有助于理解土壤侵蝕現(xiàn)象中侵蝕量的臨界坡度問題。6結(jié)論本文建立了坡面上降雨入滲產(chǎn)流的動(dòng)力學(xué)模型，用實(shí)驗(yàn)資料驗(yàn)證了模型的有效性。并通過模擬計(jì)算，分析了簡單坡面上降雨入滲產(chǎn)流的動(dòng)力學(xué)規(guī)律。得到以下幾點(diǎn)結(jié)論。(

22、1)僅降雨強(qiáng)度增大，所有參量(單寬流量，水深，流速，切應(yīng)力)均會(huì)增大。單寬流量與降雨強(qiáng)度根本呈正比關(guān)系。水深和流速與雨強(qiáng)分別有近似3/5和2/5次冪的關(guān)系，切應(yīng)力與水深有同樣的規(guī)律。(2)僅浸透率增大，所有參量減校(3)僅坡長增加，所有參量均增大。單寬流量與坡長也根本成正比。(4)僅土壤初始非飽和度(飽和含水量和前期含水量之差)增加，所有參量均減校(5)坡度增加，將會(huì)帶來幾個(gè)方面變化。坡長一樣時(shí)承雨面積將減小，所承受的雨量將減少；另一方面，坡度增大，浸透率將減小，導(dǎo)致產(chǎn)流增加。而在同樣的產(chǎn)流量下，坡度增大水流的水深減小流速將增大。綜合起來，坡長不變時(shí)流量隨坡度增大主要是減小的，僅在坡度較緩時(shí)有很少增加，水深隨坡度增大也減小，但流量變化是凸曲線，水深變化是凹曲線。流速和切應(yīng)力將先增加，在到達(dá)某一臨界坡度時(shí)變?yōu)闇p小，兩者的臨界坡度并不一樣，但都約在4050之間。參考文獻(xiàn)1J.L.P.deLia,delKINNIFfrverlandflnperviussurfae.InverlandFlA.J.ParsnsA.D.Abrahas(eds)69-88ULPress1992.2Gvindaraju,R.S.,nthediffusinavedelfrverlandflat