(1.山東電力工程咨詢院有限公司,山東 濟南 250100; 2.山東黃河勘測設計研究院,山東 濟南 250013； 3.山東正元建設有限責任公司，山東 濟南 250100)

研究降雨條件下的非飽和土滲流過程對許多巖土工程或環(huán)境巖土工程的實踐有重要意義，例如，自然或人工邊坡穩(wěn)定性分析[1-4]、礦床或城市垃圾填埋場設施的覆蓋層設計[5-6]、涉及膨脹土和濕陷性土的工程問題[7-8]，等等。解決這些問題的關鍵和難點在于計算降雨入滲過程在非飽和土體中引起的滲流場[9]。

Richards方程是描述非飽和土體中瞬態(tài)滲流的控制方程。針對該方程的求解，以往的研究多采用數(shù)值模擬手段給出數(shù)值解[10-12]，而解析方法的研究相對較少。實際上，與數(shù)值解相比，解析解更清晰簡潔，能更直觀地顯示定解問題對邊界條件的響應，且便于對各影響因素進行參數(shù)敏感性分析，并能對數(shù)值解的結果進行驗證[13]。Srivastava 等推導了均質、非均質土體降雨條件的一維瞬態(tài)入滲解析解[14]，基于此模型，詹良通[13]和李寧[9]等人分別推導了考慮坡度影響和降雨形式的滲流解析解。Chen 等[15-16]提出了非穩(wěn)恒降雨條件下的入滲解析模型。Huang等[17]和Ghotbi等[18]也分別基于Fourier積分變換和同倫分析方法(HAM)提出了降雨入滲過程的解析解。在初始狀態(tài)的處理方面，上述研究或將存在前期降雨并達到穩(wěn)態(tài)入滲時刻作為初始條件[9,13-14,16]，或假定初始含水率(水頭)分布為線性或特定函數(shù)[15,17-18]，而根據(jù)實際情況，初始狀態(tài)是可以選擇在任意時刻，相應含水率(水頭)分布往往是不規(guī)則的。誠如Zhan 等[19]所述，前期降雨所形成的初始狀態(tài)對水頭分布的影響甚至大于后續(xù)降雨。朱偉等指出在非穩(wěn)定滲流中，能不能正確地考慮初始水分量將會對計算結果產(chǎn)生非常大的影響[20]。在采用的數(shù)學算法方面，前人多借鑒Ozisik (1980)關于偏微分方程的解法[21]，在解析表達式中存在復雜的積分變換項[9，15-17]。該類方法在理論上是先進的，但數(shù)學算法稍顯復雜，不利于直接應用。因此，需要進一步研究任意初始狀態(tài)的入滲過程解析解，并考量初始狀態(tài)的影響。

本文采用指數(shù)函數(shù)描述非飽和土體的土水特征曲線及其滲透系數(shù)隨水頭的變化，對任意初始狀態(tài)的非飽和土入滲過程解析解進行了推導。最后基于所推導模型，分析初始狀態(tài)對含水率分布的影響。

1 理論推導

1.1 基本假設

(1) 土體為均質，非飽和滲流由基質吸力差(水頭差)決定，不考慮溶質吸力的影響。

(2) 不考慮非飽和土的回滯特性。

(3) 地下水位始終穩(wěn)定，不受土表降雨入滲過程的影響。

(4) 降雨過程雨強穩(wěn)定，即為均勻降雨。

1.2 控制方程

Richards為非線性偏微分方程，目前它已衍生出3種基本格式，即壓力水頭格式(h-form)、混合格式(mixed form)和體積含水率格式(θ-form)，其中混合格式被證實具有嚴格遵守質量守恒的特性，可表示為

(1)

式中，h為孔隙水壓力水頭；θ為體積含水率；k為滲透系數(shù)；t為時間；z為縱坐標，正方向向下(如圖1所示)。

圖1 入滲過程示意Fig.1 Schematic illustration of infiltration process

Yanful和Mousavi[22]指出，在非飽和土滲流中相較于大的的壓力水頭差，式(1)中的重力項可以忽略。同時，忽略重力項可以極大簡化推導過程，避免解析解中復雜的積分或級數(shù)項[23]。因此，Richards公式可以簡化為

(2)

參照文獻中常用的指數(shù)函數(shù)形式[13-17,24]，土水特征曲線和滲透系數(shù)方程可描述為

θ(h)=θr+(θs-θr)eαh

(3a)

k(h)=kseαh

(3b)

式中，θs和θr分別為飽和和殘余體積含水率；ks為飽和滲透系數(shù)；α為減飽和系數(shù)。將式(3)代入式(2)可得

(4)

為了簡便計算，定義土的擴散系數(shù)D和無量綱的體積含水率Θ分別為

(5a)

Θ=(θ-θr)/(θs-θr)

(5b)

將式(3a)、(5a)和(5b)代入式(4)，并進一步化簡，可得：

(6)

1.3 入滲過程解析解

降雨強度與土體飽和滲透系數(shù)的相互關系會對入滲過程產(chǎn)生影響。對于降雨強度小于飽和滲透系數(shù)的情況，即q≤ks，降雨將全部入滲，土表不會形成積水，也即整個過程均處于定流量入滲階段。而對于降雨強度大于飽和滲透系數(shù)的情況，即q>ks，其起始為定流量入滲階段，降雨一段時間后，表層土體達到飽和，降雨不能完全入滲，土表存在積水，定流量入滲階段轉化為土表積水階段。

1.3.1定流量入滲階段

在定流量入滲階段，上邊界為恒定的流量邊界，也即降雨強度。此時的初始和邊界條件可表示為

Θ(z,t)|z=L=1

(7a)

(7b)

Θ(z,t)|t=0=f(z)

(7c)

式中，L為地下水位深度；q為降雨強度；f(z)為任意初始狀態(tài)含水率分布函數(shù)。式(7a)表示在地下水位處，下邊界土體始終飽和，Θ始終為1；式(7b)表示上邊界為恒定流量邊界條件；式(7c)表示初始狀態(tài)。

首先求穩(wěn)態(tài)入滲時的含水率分布，設為w(z),則w(z)滿足：

(8a)

w(z)|z=L=1

(8b)

-D(θs-θr)w′(z)|z=0=q

(8c)

聯(lián)立式(8)可以求得w(z)的表達式為

(9)

由于原問題的邊界條件為非齊次(如式(7))，因此為方便求解，引入一變量使其齊次化，令：

(10)

將式(9)和式(10)代入式(6)和式(7)，可得：

(11a)

(11b)

(11c)

(11d)

這里，式(11a)為控制方程，式(11b)和式(11c)為邊界條件，式(11d)為初始條件。根據(jù)附錄中的推導過程，式(11)的解可以寫為

(12)

由于式(12)中存在積分項，較難定量求得。為簡便應用，用一個二次函數(shù)來擬合初始含水率分布，以消除積分變量，簡化計算。由于z=L時，Θ始終為1，因此令：

(13)

式中，a和b為擬合參數(shù)。將式(9)和式(13)代入式(12)，并進一步化簡，可得定流量入滲階段土體含水率瞬態(tài)分布的解析解表達式：

(14)

1.3.2土表積水階段

對于q>ks的情況，假設土表積水開始時間為tp。不考慮積水造成的土表水頭堆積的影響，則當t>tp時，上邊界條件為恒定含水率邊界條件。

首先對于積水時刻tp，可令式(14)等于1求得。設此時的含水率分布為p(z)，則p(z)的表達式為

(15)

當t>tp時，定義一個新的時間變量T=t-tp，將T代入式(6)可得：

(16)

在土表積水階段，初始和邊界條件可以表示為

Θ(z,T)|z=L=1

(17a)

Θ(z,T)|z=0=1

(17b)

Θ(z,T)|T=0=p(z)

(17c)

根據(jù)Carslaw和Jaeger(1959)[25]的解法，式(16)在邊界條件式(17)下的解為

(18)

由于式(18)含p(z)函數(shù)未解積分項，若直接將式(15)代入求解，算法復雜，且不能得到其解析表達式，限制其直接應用。因此，我們用二次函數(shù)來代替式(15)，以實現(xiàn)簡化。因為當T=0時，在z=0和z=L處，Θ=1，所以二次函數(shù)的表達式可寫為

(19)

其中，c為擬合參數(shù)，利用式(19)擬合式(15)來可以求得其值。將式(19)代入式(18)可求得存在土表積水階段的含水率分布解析解：

(20)

2 算例分析和討論

本文針對降雨強度和飽和滲透系數(shù)的關系，分別進行如下兩個算例討論。如表1所示，在算例1中降雨強度q小于飽和滲透系數(shù)ks；而在算例2中，q大于ks。對算例1和算例2，均進行兩種初始條件(線性和非線性含水率分布)的分析。

表1 邊界條件和土性參數(shù)Tab.1 The boundary conditions and soil properties

圖2為降雨強度小于飽和滲透系數(shù)時，不同時刻土體內(nèi)含水率的分布，其值均由式(14)計算得出。圖2(a)中，a和b分別賦值0和1來擬合線性初始含水率分布，而圖2(b)中，a和b分別為2和-1.5。比較圖2(a)和(b)，不難發(fā)現(xiàn)初始狀態(tài)對含水率分布產(chǎn)生非常大的影響。另外，圖2(b)的初始狀態(tài)可認為存在較強的前期降雨，由于后期降雨強度小于飽和滲透系數(shù)，土體內(nèi)部滲流較快，使得表層土體含水率在t=5 h小于初始含水率。

圖3展示了降雨強度大于飽和滲透系數(shù)時，不同時刻土體內(nèi)含水率的分布。其初始狀態(tài)設定與圖2情況相同。對圖3(a)和(b)，其積水時刻tp分別為9.69 h和9.39 h，土表邊界在ttp時上邊界始終維持在飽和含水率，含水率分布由式(20)計算求得，且以tp時刻的含水率分布為初始條件。對應圖3中的兩種

圖2 算例1中不同時刻土體內(nèi)含水率分布Fig.2 Water content redistribution during infiltration for Case 1

圖3 算例2中不同時刻土體含水率分布Fig.3 Water content redistribution during infiltration for Case 2

情況，式(20)中的擬合參數(shù)c分別為1.53和1.68。隨著時間推移，含水率分布的變化幅度逐漸變小，并最后趨于穩(wěn)態(tài)。由圖3可見，初始條件不僅影響含水率的大小，而且影響含水率分布形態(tài)。

3 結 論

(1) 基于指數(shù)函數(shù)來描述土水特征曲線和滲透系數(shù)隨水頭的變化，得到降雨強度大于或小于飽和滲透系數(shù)的入滲過程解析解。該模型可以考慮任意初始狀態(tài)，并且避免了復雜的數(shù)學算法表達，無需借助編程來實現(xiàn)數(shù)學計算，易于應用。計算結果表明初始狀態(tài)對入滲過程的計算結果有顯著影響。

(2) 提出的模型能夠揭示特定情況下的入滲過程機理，形式簡單，能克服數(shù)值模擬的復雜數(shù)學算法及計算收斂問題，便于考察入滲過程的機理。

(3) 但同時需要考慮到，解析解的得出是基于一定的假設和簡化。實際應用中，應將解析解手段和數(shù)值模擬方法結合應用。