基于降雨入滲全過程的非飽和濕潤峰模型

徐 旭，席 越，姚文娟

（1.上海大學(xué) 土木工程系，上海 200444；2.上海大學(xué) 力學(xué)與工程科學(xué)學(xué)院，上海 200444）

1 研究背景

自改革開放以來，中國的城鎮(zhèn)化建設(shè)進(jìn)程不斷加快。但由于粗放的發(fā)展模式，城市自然水文循環(huán)過程遭到破壞，從而導(dǎo)致城市內(nèi)澇災(zāi)害的頻繁發(fā)生［1］。因此，緩解城市管網(wǎng)的排水壓力，優(yōu)化生態(tài)環(huán)境，最大限度地實現(xiàn)雨水在城市區(qū)域的積存、滲透和凈化，已成為人類生存需要解決的重大問題［2］。為此，許多學(xué)者進(jìn)行了降雨入滲過程的研究［3-6］，以便充分發(fā)揮土壤的蓄水能力，準(zhǔn)確預(yù)測徑流。

目前，具有代表性的降雨入滲計算方法為Mein-Larson降雨入滲模型［7］和數(shù)值模擬方法［8-10］。Mein-Larson模型由于基于物理假設(shè)且計算簡便，從而被廣泛使用。但由于該模型采用水平濕潤峰假設(shè)，假定土壤濕潤區(qū)均為飽和含水率，土壤入滲過程中的水力傳導(dǎo)度取為飽和滲透系數(shù)，這會導(dǎo)致土壤蓄水量被高估。從而學(xué)者們不斷地改進(jìn)該模型。李寧等［11］引入Van Genuchten水力傳導(dǎo)函數(shù)，忽略土壤吸力梯度對入滲過程的影響，建立了可以反映低強(qiáng)度降雨條件下非飽和土入滲過程的改進(jìn)模型。于寧宇等［12］引入Van Genuchten基質(zhì)吸力模型考慮土壤基質(zhì)吸力對入滲過程的影響，建立了適用于吸力顯著土壤的入滲模型，并且可體現(xiàn)土壤含水率隨入滲時間的變化。唐揚等［13］在模型中加入了土壤初始含水率的分布情況對降雨入滲的影響。以上修正模型均采用水平濕潤峰假設(shè)，而實際入滲過程中土壤水分剖面呈明顯的層狀分布［14］。由于該假定過于簡化，從而限制了模型的計算精度。同時，也會影響降雨入滲過程中各時間節(jié)點的確定。

本文參考前人的實驗以及數(shù)值模擬結(jié)果，對降雨入滲全過程進(jìn)行分段，假設(shè)各階段土壤水分剖面的分布狀態(tài)，提出修正的濕潤峰模型。同時，結(jié)合達(dá)西定律，引入土體非飽和參數(shù)，使修正后的模型可以更加真實地反映實際入滲過程，并且能夠更加精確地預(yù)測土壤累計入滲量和積水時間。

2 模型建立

降雨入滲按照降雨強(qiáng)度與土壤飽和滲透系數(shù)的相對大小主要可分為兩類。本文分別針對以下兩種降雨條件建立修正模型。

2.1 降雨強(qiáng)度小于土壤飽和滲透系數(shù)當(dāng)降雨強(qiáng)度小于土壤飽和滲透系數(shù)時，整個降雨入滲過程為非飽和入滲，所有雨水均滲入土壤，故不產(chǎn)生積水，又可稱為無壓入滲。設(shè)降雨強(qiáng)度為r，土壤入滲率為q，降雨入滲可簡化為一維入滲情況，由達(dá)西定律可得：

式中：K（θ）為土壤非飽和滲透系數(shù)；zf為濕潤區(qū)浸潤深度；h1、h2分別為土壤表層和濕潤峰處基質(zhì)勢吸力水頭。

降雨入滲初期，土壤中形成的濕潤區(qū)域不斷擴(kuò)大，其表層含水率及浸潤深度同時增加，土壤含水率沿深度方向逐漸減小。王文焰等［15］根據(jù)實驗結(jié)果，假定入滲過程中非飽和區(qū)含水率剖面為1/4橢圓曲線。彭振陽等［16］同樣采取該假定，結(jié)合Richards模型計算結(jié)果對Green-Ampt模型進(jìn)行改進(jìn)，取得了良好的計算效果。本文同樣假設(shè)在降雨入滲初期，土壤水分剖面呈橢圓形分布。設(shè)土壤表層含水率為θi，初始含水率為θ0，其濕潤區(qū)平均含水率為：

本文引入Van Genuchten模型［17］，如式（3）、式（4）所示。以濕潤區(qū)的平均含水率所對應(yīng)的土壤基質(zhì)勢和非飽和滲透系數(shù)等效替代實際的土壤基質(zhì)勢以及非飽和滲透系數(shù)。

由于降雨強(qiáng)度小于土壤飽和滲透率，其入滲率q等于降雨強(qiáng)度r，將式（2）、式（3）和式（4）帶入式（1）中可得：

設(shè)降雨歷時為T，根據(jù)水量守恒定律可知，土壤累計入滲量等于其含水量增量，即得：

將式（6）帶入式（5）可消去zf，即得：

上式即為降雨初期土壤入滲的控制方程，可通過迭代的方法計算出給定T時刻所對應(yīng)的土壤水分剖面分布。

表層土壤最終的入滲率為表層土體含水率的函數(shù)［18］，由式（1）可知其數(shù)值上等于表層土壤含水率所對應(yīng)的非飽和滲透系數(shù)。對于恒定雨強(qiáng)的降雨，隨著入滲過程的推進(jìn)，土壤表層含水率逐漸增加至一定值后不再發(fā)生變化［19］，將該值稱之為臨界含水率（θc）。此時，降雨強(qiáng)度與表層土壤的入滲率相等。故臨界含水率可由式（8）計算得出。

另外，有試驗結(jié)果顯示，在降雨入滲后期，土壤的浸潤深度隨時間呈線性增長［20］，且土壤過渡層所占總濕潤層的比例線性下降［16］。而且根據(jù)數(shù)值模擬結(jié)果可知，對于恒定降雨強(qiáng)度的入滲，過渡層的厚度將逐漸趨于穩(wěn)定。基于以上原因，本文假設(shè)已形成的橢圓形濕潤區(qū)形狀不再發(fā)生變化，并稱之為無壓入滲過渡區(qū)，隨著入滲進(jìn)程均速向下推進(jìn)。過渡區(qū)以上的濕潤區(qū)含水率均為臨界含水率，土壤水分剖面開始出現(xiàn)分層現(xiàn)象。

為方便敘述，稱出現(xiàn)分層前的階段為無壓入滲過渡階段，分層后為無壓入滲穩(wěn)定階段。對于無壓入滲穩(wěn)定階段，根據(jù)水量守恒定律可得：

式中：Tc為土壤表層含水率由初始含水量增長至臨界含水率所用的時間，稱為臨界時間。zfc為過渡區(qū)的深度范圍。由式（9）可知，在無壓入滲穩(wěn)定階段土壤濕潤區(qū)深度隨入滲時間線性增長，其變化規(guī)律與試驗結(jié)果相吻合。

設(shè)累計入滲量為I，基于以上的假定與推導(dǎo)，累計入滲量可以用下列表達(dá)式計算：

2.2 降雨強(qiáng)度大于土壤飽和滲透系數(shù)當(dāng)降雨強(qiáng)度大于土壤飽和滲透系數(shù)時，整個入滲過程將經(jīng)歷三個階段：降雨強(qiáng)度控制入滲階段，非飽和狀態(tài)控制入滲階段和飽和狀態(tài)控制入滲階段。

降雨初期土壤的入滲能力大于降雨強(qiáng)度，此時土壤入滲率取決于降雨強(qiáng)度，故稱為降雨強(qiáng)度控制入滲階段。該階段內(nèi)所有雨水均滲入土壤，不產(chǎn)生地表徑流。與降雨強(qiáng)度小于土壤飽和滲透系數(shù)時的無壓入滲過渡階段相同，可同樣采用式（7）作為該階段的控制方程。

由于降雨強(qiáng)度大于飽和滲透系數(shù)，表層土壤的滲透系數(shù)無法增長至與降雨強(qiáng)度相等，表層土壤的含水率先達(dá)到飽和含水率（θs）后保持不變。由式（5）可知，當(dāng)土壤濕潤區(qū)平均含水率不變時，隨著土壤濕潤區(qū)深度增加，土壤的入滲率不斷減小，導(dǎo)致此后的降雨強(qiáng)度大于土壤的入滲率，部分雨水無法滲入土壤，從而形成地表徑流。降雨入滲過程開始進(jìn)入積水入滲階段。

在積水入滲初期，土壤的入滲能力小于降雨強(qiáng)度，但大于土壤飽和狀態(tài)下的入滲率（即飽和滲透系數(shù)）。該階段土壤的入滲率取決于非飽和土的入滲能力。假定該階段土壤含水率

剖面仍呈橢圓形分布，土壤表層含水率為定值θs，故土壤濕潤區(qū)的平均含水率為：

由水量守恒定律可得：

將式（5）、式（11）代入式（12），整理可得：

其中，zfp=4·r·T［/π·（θs-θ0）］，為積水點所對應(yīng)的土壤浸潤深度。式（14）即為非飽和狀態(tài)控制入滲階段的控制方程。

當(dāng)土壤入滲能力減小至等于土壤飽和滲透系數(shù)時（即q=Ks時），由于此后的入滲速率恒定，與無壓入滲穩(wěn)定階段類似，故同樣假設(shè)已形成的橢圓形濕潤區(qū)形狀不再發(fā)生變化，稱之為積水入滲過渡區(qū)，隨著入滲進(jìn)程向下均速推進(jìn)。過渡區(qū)以上的濕潤區(qū)含水率均為飽和含水率，土壤水分剖面開始出現(xiàn)分層現(xiàn)象。降雨入滲過程開始進(jìn)入飽和狀態(tài)控制入滲階段。根據(jù)水量守恒定律可得：

式（15）即為飽和狀態(tài)控制入滲階段的控制方程。其中，zfs為積水入滲過渡區(qū)的深度范圍，可由式（16）計算得出。Ts為達(dá)到飽和狀態(tài)控制入滲階段所需的時間，可將zfs代入式（14）計算得出。

為方便后文敘述，類比前一種降雨條件下的入滲過程，將積水入滲非飽和狀態(tài)控制階段與飽和狀態(tài)控制階段稱之為積水入滲過渡階段與穩(wěn)定階段。

設(shè)累計入滲量為I，基于以上假定與推導(dǎo)，各階段的累計入滲量可以用下列表達(dá)式計算：

3 結(jié)果分析與驗證

3.1 計算結(jié)果分析本節(jié)首先利用上節(jié)提出的修正模型計算土壤在不同降雨強(qiáng)度下的入滲情況。修正模型的計算流程如圖1所示，本文采用Matlab軟件實現(xiàn)整個計算過程。計算過程中的土質(zhì)參數(shù)取自文獻(xiàn)［13］，土質(zhì)類型為粉質(zhì)黏土，考慮土壤含水率沿深度方向均勻分布的情況，不考慮地下水埋深的影響，設(shè)土壤初始含水率為0.2。具體水力參數(shù)見表1。

表1 土體水力參數(shù)

圖2為當(dāng)降雨強(qiáng)度分別為10 mm/h和30 mm/h時，土壤水分剖面在不同時刻的分布情況。兩個降雨強(qiáng)度分別對應(yīng)上節(jié)中提到的兩種降雨入滲情況。從圖中可見，隨著入滲時間的增加，土壤的濕潤區(qū)范圍不斷擴(kuò)大。本模型可以反映在降雨入滲過程中，土壤含水率分布沿深度方向并不均勻。地表含水率較高，土體內(nèi)部含水率隨土壤深度逐漸減小。另外，圖中可以明顯的顯示出土壤含水率剖面的分層現(xiàn)象。相較于傳統(tǒng)的水平濕潤峰假設(shè)模型，本模型更符合實際的入滲狀況。

不同降雨強(qiáng)度下，土壤表面的含水率隨時間的變化趨勢如圖3所示。降雨初期，土壤的表層含水率呈非線性增加。當(dāng)入滲持續(xù)一段時間后，土壤表層含水率將會逐漸趨于一定值，即臨界含水率。當(dāng)降雨強(qiáng)度小于土壤飽和滲透系數(shù)時，臨界含水率小于土壤飽和含水率，并且與降雨強(qiáng)度成正比。當(dāng)降雨強(qiáng)度大于土壤飽和滲透系數(shù)時，臨界含水率等于土壤飽和含水率。當(dāng)降雨強(qiáng)度越大，表層土壤達(dá)到臨界含水率所需的時間越短。

圖4為本文提出的修正模型在不同降雨強(qiáng)度下，土壤浸潤深度隨時間的變化。由圖可知，降雨前期，各降雨強(qiáng)度下土壤浸潤深度隨時間呈非線性增長。降雨后期，土壤浸潤深度隨時間呈線性增長。對于降雨強(qiáng)度小于土壤飽和滲透系數(shù)的入滲過程，降雨強(qiáng)度越大，其增長速率越大。當(dāng)降雨強(qiáng)度小于土壤飽和滲透系數(shù)時，土壤浸潤深度曲線在降雨后期基本平行，且其數(shù)值相差不大。該結(jié)果與實際入滲規(guī)律更加接近。

圖1 修正模型計算流程圖

圖2 不同時間不同降雨條件下土壤含水率分布

圖3 土壤表層含水率隨時間的變化

不同降雨強(qiáng)度下土壤的累計入滲量如圖5所示。兩種降雨條件下，土壤累計入滲量隨時間的變化規(guī)律有所不同。當(dāng)降雨強(qiáng)度小于土壤飽和滲透系數(shù)時，土壤的累計入滲量呈線性增長，且其增長速率等于降雨強(qiáng)度。當(dāng)降雨強(qiáng)度大于土壤含水率時，土壤的累計入滲量的增長分為三個階段。降雨初期呈線性增長，增長速率等于降雨強(qiáng)度。降雨中期呈非線性增長，此時的增長速率由土壤的非飽和性質(zhì)控制。降雨后期呈線性增長，增長速率等于土壤飽和滲透系數(shù)。故從圖中可以看出，在降雨后期，當(dāng)降雨強(qiáng)度為30 mm/h和40 mm/h時，兩者的累計入滲量的增長曲線平行。另外，兩者的累計入滲量也十分接近，說明對于高強(qiáng)度降雨條件下的降雨入滲，最終的入滲量主要取決于土壤飽和滲透性質(zhì)。

圖4 土壤浸潤深度隨時間的變化

圖5 累計入滲量隨時間的變化

3.2 模型驗證為驗證修正模型計算結(jié)果的準(zhǔn)確性，將修正模型和數(shù)值模型的計算結(jié)果進(jìn)行對比。另外，由于Mein-Larson模型為采用“水平濕潤峰”假設(shè)的經(jīng)典降雨入滲模型，故將其計算結(jié)果加入比較，對本文修正方法的效果進(jìn)行評價。將傳統(tǒng)模型與修正模型計算結(jié)果稱為累計入滲量的預(yù)測值，作為縱軸。將有限元軟件的模擬結(jié)果稱為累計入滲量的計算值，作為橫軸。當(dāng)數(shù)據(jù)點越靠近對角線，則說明兩者吻合愈好。本文選擇土壤浸潤深度與累計入滲量的計算結(jié)果作為驗證項目。對于累計入滲量的對比，僅選擇降降雨強(qiáng)度大于土壤飽和滲透系數(shù)的情況，即降雨強(qiáng)度為30 mm/h和40 mm/h的計算結(jié)果。具體對比結(jié)果見圖6及圖7。

圖6 浸潤深度對比

圖7 累計入滲量對比

由圖6可知，傳統(tǒng)模型和修正模型的浸潤深度計算值根據(jù)降雨強(qiáng)度不同，其相對大小有所差異。當(dāng)降雨強(qiáng)度小于土壤飽和滲透系數(shù)時，由于兩種模型入滲量計算值相等，本文所提出的模型考慮了該情況下土壤處于非飽和狀態(tài)并且含水率分布不均，故修正模型的計算結(jié)果大于傳統(tǒng)模型的計算結(jié)果。

當(dāng)降雨強(qiáng)度大于土壤飽和滲透系數(shù)時，如圖7所示，由于本本文模型考慮前期非飽和入滲階段，故其入滲量計算值小于傳統(tǒng)模型計算值，并使其浸潤深度計算值將小于傳統(tǒng)模型計算值。另外，本文模型的計算結(jié)果與模擬結(jié)果的最大誤差小于5 cm，相對誤差均小于5%。根據(jù)以上分析結(jié)果可知，本文提出的修正模型在傳統(tǒng)模型的基礎(chǔ)上，修正效果較好，對于計算結(jié)果的準(zhǔn)確性方面有了很大的提高。

積水時間也是降雨入滲過程的重要指標(biāo)之一。表2即為不同計算方法所得的積水時間的對比。從表中的數(shù)據(jù)可以看出，修正模型相較于傳統(tǒng)模型更加接近數(shù)值解，說明本文的修正模型對于節(jié)點的假定合理，當(dāng)土壤表層含水率達(dá)到飽和含水率時開始出現(xiàn)積水。

表2 積水時間的模型解及數(shù)值解的比較

為了檢驗本文提出的模型對于其他類型的土壤的適用性，本文另外選取3組不同土質(zhì)土壤進(jìn)行計算，并與有限元結(jié)果進(jìn)行比較。具體土體參數(shù)［16，21］見表3。4組土壤的土-水特征曲線如圖8所示。

圖9即為4種土質(zhì)在不同降雨強(qiáng)度下，兩種計算方法所得的累計入滲量的對比圖。雖然圖中顯示有數(shù)據(jù)點偏離較大，但實際誤差均小于2 cm。砂土相差最大，基本差值范圍在1～2 cm之間。粉土與壤土相差范圍較小，其范圍在0.1～0.5 cm之間。由此表明本文提出的修正模型對于不同土質(zhì)土壤適應(yīng)性良好。

綜上所述，本文所建立的修正模型假定合理，對于兩種降雨入滲情況計算結(jié)果的準(zhǔn)確性較好，其計算結(jié)果滿足誤差范圍，并且更加接近真實的入滲情況。另外，本文模型也可以適用于其他類型的土壤。

表3 土體水力參數(shù)

圖8 不同土質(zhì)土壤的土-水特征曲線

圖9 不同類型土的模型解與數(shù)值解的對比

4 結(jié)論

（1）與現(xiàn)有的入滲模型相比，本文提出的入滲模型符合實際入滲過程中的土壤含水率分布，能夠更加真實地體現(xiàn)出土壤表層含水率和浸潤深度隨時間的變化趨勢。

（2）對于降雨強(qiáng)度小于飽和滲透系數(shù)和降雨強(qiáng)度大于土壤飽和滲透系數(shù)的兩種降雨入滲情況，本文模型所得的計算結(jié)果與有限元軟件模擬結(jié)果吻合，明顯優(yōu)于傳統(tǒng)模型。

（3）當(dāng)降雨強(qiáng)度大于土壤飽和滲透系數(shù)時，當(dāng)表層土壤飽和后，即出現(xiàn)積水。土壤累計入滲量最終取決于土壤的飽和滲透性質(zhì)，降雨強(qiáng)度對其影響有限。

（4）本文明確了降雨入滲過程中各階段的狀態(tài)和節(jié)點發(fā)生時間，與有限元方法相比更加簡單，便于計算并且計算效率得到提高，且具有良好的計算精度，因此可以利用本模型對積水時間與土壤蓄水量進(jìn)行計算和預(yù)測。