涌灌條件下非飽和土壤水分運動的數值模擬

張巧利，張 亮

(1.河南廣播電視大學機電工程學院,鄭州 450046；2.鄭州職業技術學院，鄭州 450000)

0 引 言

涌流灌溉是一種新型的適用于果樹及植樹造林的節水灌溉技術，它由滴灌演變而來，通過土壤根系層附近開挖一定深度灌水孔形成柱狀輻射水流入滲，形成類似橢球體濕潤體[1]。由于涌流灌溉入滲通常在地下進行，無污染，水分蒸發損失小，可以有針對性的濕潤作物根部區域土壤，作物能夠有效吸收并利用水分。涌流灌溉方式有效儲水空間大，毛管直接深入灌水孔，一旦堵塞，容易發現并更換，從而有效解決了滴灌、滲灌在山地經濟林灌溉時易堵塞和老化的問題，同時它對地形的適應性也很強。因此，涌流灌溉作為一種新型的節水方式，有一定的發展前景。研究不同灌水情況下土壤入滲濕潤體特性有重要的意義。通過室內外試驗可以獲得土壤不同分布界面上水分運動狀況，但由于試驗的局限性無法對各種狀況進行全面分析，且試驗過程存在較大誤差。近年來，以非飽和水運動理論為基礎的數值模擬方法被用來定量分析微灌條件下土壤水分分布狀況。計算機模擬方法簡便、靈活和快捷，在確定的初始、邊界條件下，可以模擬不同因素組合下土壤濕潤體的水分狀況，從而為合理確定灌水設計參數提供方便的研究手段。

國內外學者通過大量試驗手段和數值模擬方法研究了土壤水分及溶質運移[2-6]。Cote等[7]人運用HYDRUS-2D軟件模擬了滴灌入滲條件下土壤水分運移；許迪等[8]人通過建立點源土壤水分運動方程與溶質對流-彌散方程，并采用HYDRUS-2D求解了模型；周青云等[9]用有限差分法求解了蓄水坑灌條件下土壤水分入滲數學模型；依據非飽和土壤水分運動理論，李光永等[10-12]人建立了點源滴灌條件下土壤水分運動動力模型，并對模型進行求解；李明思[13]在圓柱體濕潤體模型的基礎上，建立了不同流量入滲下滴灌入滲的數學模型，并通過室內試驗研究了不同滴頭流量下各因子的變化過程，模型研究為滴灌系統的設計提供了重要的參考依據；本文依據非飽和土壤水動力學理論，試圖建立涌流灌灌溉條件下土壤入滲數學模型，并通過室內試驗對模型正確性予以驗證。

1 涌流灌溉條件下土壤水分運動的數學模型

1.1 模型條件假定

模型都是在理想條件下提出的，與實際有一定出入，須作如下假定：假定土壤為固相剛性介質，各項同性，土壤水分不可壓縮，忽略容質勢、溫度勢作用；忽略土壤水分運動滯后效應；土壤基質勢和土壤含水率的關系在土壤水分特征曲線中呈單值函數，土壤的密度(體積)和機械組成決定函數的形狀。

1.2 土壤水分運動方程

水分運動邊界為柱狀界面，屬于三維軸對稱入滲，依據達西定律及質量守恒定律，通過以下偏微分方程來描述土壤水分運動：

(1)

式中：h為土壤總水勢[L]；kx(θ)，ky(θ)，kz(θ)分別為土壤在x、y、z三個方向上的土壤非飽和導水率[LT-1]；θ為體積含水率[L3L-3]；t為計算時間[T]；x、y、z分別為三維坐標[L]。

1.3 初始及邊界

(1)初始條件[14]。

h0(x,y,z)=aθ-b0，t=0

(2)

式中：h0為土壤剖面初始負壓水頭，cm；θ0為土壤初始體積含水量，cm3/cm3；假定模型計算區土壤水分均勻分布，通過選定的土壤水分特征曲線將已知的土壤初始含水率轉化為需要模擬時初始輸入的土壤負壓剖面。其中a，b為土壤水分特征曲線常數。

(2)邊界條件。圖1為數值模擬計算區示意圖。 均質土壤水分運動可以視為軸對稱的三維入滲問題，假定計算區域較大，則水分在入滲過程中沒有到達邊界a′a″b″b′，a″c″b″無水量交換，為零通量界面。

圖1 計算區示意圖

則a′a″b″b′邊界：

h0(x,y,z)=aθ-b0，t>0，x=xa，y=yb

(3)

a″c″b″邊界為下邊界，處理時考慮地下水埋深較大，看作第一類邊界：

h0(x,y,za″c″b″)=h0，t>0″

(4)

a′a″c″c′與c′b′b″c″為對稱邊界，看作不透水邊界，水通量為零。

(5)

a′b′c′邊界為上邊界，豎直方向水流通量為降雨強度或蒸發強度q0：

(6)

邊界o′o″d′d為柱狀入滲界面，由試驗觀測可知，有明顯的飽和積水區，為運動的邊界，較為復雜，一是灌水初期，灌溉出流速率小于土壤入滲速率時，灌水器所供應的水分很快滲入到土壤中，孔洞中未形成積水，土壤入滲速率等于邊界流量，在模型中根據所模擬計算結果中累積入滲量來推求邊界入滲速率；二是隨著入滲的進行，孔洞周圍土壤含水率達到飽和含水率，孔洞中形成積水，積水深度隨時間呈漸進增加趨勢，一定時間后達到穩定，此時邊界按動水頭邊界處理計算。ho′o″d′d=h′(h′通過試驗實測獲得，試驗時孔外安裝測壓管來觀測入滲水頭數據)。

1.4 模型求解

1.4.1 HYDRUS-3D

HYDRUS-3D是一款可以用來模擬土壤水分及溶質運移的有限元分析軟件。模型中的水流狀態為非飽和三維達西水流，計算過程忽略溫度和空氣對水流的影響，控制方程采用Richards方程，方程中可以通過嵌入匯源項來考慮作物根系對水分的吸收。程序可以靈活處理多種復雜邊界，有固定流量邊界、動水頭和定水頭邊界、大氣邊界及排水溝等。計算區域可以由不規則或各項異性非均質土壤構成。水流計算區域可以通過三棱柱網格劃分，采用伽遼金線狀有限元法求解控制方程。時間離散對飽和及非飽和均采用隱式差分法。離散化后的非線性控制方程組可采用迭代法將其線性化。

1.4.2 模型求解

利用HYDRUS-3D軟件求解以上所建模型[15]。模型計算區域為長方體，灌水孔置于中間，尺寸根據試驗決定，采用棱柱體(1.25×1.25×1.25×2.5)對計算區進行劃分，如圖2所示。由于孔洞周圍水勢梯度變化較大，為了提高計算的精確度，該處適當對網格體加密。模擬時間采用變時間步長計算，初始時間步長定為1 min，最大和最小時間步長定為10和3 min。

圖2 模擬計算流區的網格劃分

2 試驗驗證

2.1 土壤水動力學特性

試驗土壤土質為砂壤土，土壤容重為1.32 g/cm3，土壤穩定入滲率為1.5 mm/min，土壤平均含水率為10.4%。試驗土壤取自河南農科院試驗基地，物理性狀為：粒徑>0.02 mm占55.29%，0.02～0.002 mm占33.87%，粒徑<0.002 mm占10.84%。

根據經驗采用壓力膜儀來測定土壤水分特征曲線，定水頭法測定土壤飽和導水率。采用Van Genuchten模型對土壤水分特征曲線θ(h)和土壤導水率K(θ)進行擬合。

(8)

其中：

(9)

式中：以上方程中有5個獨立參數θs,θr,α,n,K，其中：θs為土壤飽和含水率， cm3/cm3；θr為土壤殘余含水率，cm3/cm3；m、n、α為模型的擬合參數，其中：m=1-1/n，α的取值與土壤物理性質有關，cm-1；Ks為土壤滲透系數，cm/min；Se為有效含水量，cm3/cm3；L為孔隙連通性參數，對大多數土壤一般取為0.5。供試土壤水分特性參數為：θr=0.028 3 cm3/cm3，θs=0.467 2 cm3/cm3，α=0.035 7 cm-1，h=1.48，Ks=0.125 cm/min。

2.2 測試方法

為了驗證模型的正確性以及所選參數的合理性，采用室內土箱入滲進行物理模型試驗。試驗裝置由土箱和供水系統兩部分組成，土箱由厚度10 mm的有機玻璃組成，規格為80 cm×80 cm×100 cm(長·寬·高)。試驗前將土壤風干過篩，按試驗要求密度分層裝土。供水采用馬氏瓶，土箱中部孔洞周圍用套管固定，管壁開孔，用紗網包裹，防止土壤進入孔洞堵塞出水孔，毛管伸入孔洞內進行灌水，灌水量通過開關進行調節。試驗結束后在以孔洞為中心一側對稱斷面不同位置取土測定含水率。

2.3 結果驗證

模擬計算了地表以下不同位置處(觀測點布置如圖3)涌流灌溉下土壤中心斷面含水率分布情況。圖4為模擬最終值與實測結果在灌水量14 L結束后距灌水孔不同位置處土壤含水量對比情況。由于模擬計算條件均為理想狀態，在初定的參數求得計算結果后與實測數據進行比對，利用實測數據反推模型中與土壤物理性質有關的參數，發現影響模擬結果的主要參數是m、α，通過修正參數計算結果，直到模擬值能較好地與實測值吻合為止。為了保證試驗數據的準確性，同一組試驗各做三次，取平均值作為最后的試驗結果。圖中對灌水量為14 L灌水420 min情況的模擬與實測結果對比分析表明，不管是水平方向或者垂直方向，灌水孔周圍水勢梯度變化大的點處土壤水含水率計算值與實測值比較接近，表明用數值模擬方法可以較好地表示濕潤體內水分分布狀況。

圖3 觀測點位置圖

圖4 土壤含水率模擬與實測值對比

表1為土壤入滲420 min后距離灌水孔不同位置觀測點處含水率分布狀況實測與模擬對比情況，觀測點布置如圖3所示。通過實測值與模擬值對比驗證，誤差定義為：(實測值-模擬值)/實測值。結果表明，相對誤差基本保持在±10%以內，兩者一致性較好。

3 模型的應用

通過模型對以上土壤參數條件下地下涌流灌溉土壤水分狀況進行初步模擬分析。

3.1 土壤濕潤體內水分分布狀況

作為一種新型的節水灌溉方式，掌握涌灌土壤濕潤體內含水量分布規律，對于科學合理施行作物水分管理，合理布置監測儀器有重要的意義。試驗實測獲得土壤水分在不同土質、流量、初始含水率下的動態分布狀況，耗時費力，且試驗過程繁瑣。通過數值模擬可以較好地解決以上問題，且模擬能夠為延伸試驗提供可靠的結果。圖5為灌水量為14 L，灌水420 min情況的入滲濕潤體內不同監測點處土壤含水率及負壓水頭隨時間的動態變化情況，從中可以得知，靠近灌水孔處土壤含水量最大，增加速率最快，當土壤含水量接近土壤飽和含水量時，增加速率明顯減緩，甚至不再增加。計算結果能連續反映土壤水分動態分布規律。

模擬土體具有相同的土壤初始含水率。從模擬得到的結果可以得出，灌水孔周圍的土壤水壓力幾乎接近正壓，含水率越接近飽和含水率。距離灌水孔水不同位置處對土壤水分的影響主要是濕潤鋒到達的時間不同，含水率增長速率隨著遠離灌水孔越小，且最終含水率也越小。

圖5 涌流灌溉下入滲剖面土壤模擬水分分布狀況

3.2 土壤濕潤體分布

在砂壤土條件下分別對灌水流量為2.0、3.6、5.4 L/h 3種情況進行模擬計算，分析灌水流量對土壤水分分布狀況的影響，其中，模擬的土體土壤初始含水率及灌水量均一致。圖6給出了3種不同流量下灌水結束后土壤水分等值線圖及土壤含水量分布，由圖6可知，灌水量相等，隨著流量的增加，濕潤圈形狀呈水平方面束窄，垂直方向伸長的變化趨勢，這是由于土壤水分在重力勢及基質勢的共同作用下運動，流量增大后，土壤濕潤體內部含水率迅速增大，灌水孔中積水深度也變大，所以濕潤體上半部分擴大較快；土壤含水率較大的范圍增加，在土壤濕潤體內部水分分布上，相同位置處土壤含水率隨流量的增大而增大。綜上可知，灌水流量較大時，土壤中靠近灌水孔附近局部土壤含水率在較短時間內達到較高值，從而水勢梯度變大，促進土壤水分擴散。但是流量過大會破壞周邊土壤結構，在來不及擴散的情況下孔洞內產生積水，向上入滲濕潤體擴散至地表，造成水分蒸發損失。所以選擇流量時，要充分考慮所選流量既不會破壞周圍工作條件亦不會發生深層滲漏和擴散到地表造成水分蒸發損失。

圖6 不同流量濕潤體含水率等勢線圖(單位：cm3/cm3)

4 結論與討論

(1)本文依據非飽和土壤水動力學理論，建立了涌流灌溉條件下土壤水分運動數學模型，通過采用HYDRUS-3D對所建模型進行模擬，模擬結果采用室內試驗實測數據進行了驗證，通過實測數據對模型參數進行修正，發現參數α、m對模型計算影響較大。最后經試驗數據修正以后的參數代入模型計算，結果與實測數據吻合較好。

(2)應用所建模型分析了地下涌流灌溉土壤水分運動的一些規律，從而可以為涌流灌溉系統的合理設計和運行提供理論指導。

(3)由于涌流灌溉系統流量一般較大，根系吸收水分對土壤水分運動有一定影響，在今后的研究中有待進一步增加根系吸收模型。

[1] 吳普特，朱德蘭.涌泉根灌技術研究與應用[J].灌排機械工程學報，2010，28(4)：355-357.

[2] Thorburn P J,Cook F J,Bristow K L.Soil-dependent wetting from trickle emitters:implications for system design and management[J].Irrigation Science,2003,22(3/4):121-127.

[3] Provenzano G.Using HYDRUS-2D simulation model to evaluate wetted soil volume in subsurface drip irrigation systems[J].Journal of Irrigation and Drainage Engineering,2007,133(4):342-349.

[4] Zhou Q Y,Kang S Z,Zhang L,et al.Comparison of APRI and Hydrus-2D models to simulate soil water dynamics in a vineyard under alternate partial root zone drip irrigation[J].Plant Soil,2007,291(1/2):211-223.

[5] 趙偉霞，蔡煥杰，陳新明，等.無壓灌溉土壤濕潤體含水率分布規律與模擬模型研究[J].農業工程學報，2007，23(3)：7-12.

[6] 康銀紅，馬孝義，李 娟，等.黃土高原重力式地下滴灌水分運動模型與分區參數研究[J].農業機械學報，2008，39(3)：90-95.

[7] Cite CM.Bristow K L,Charlesworth P B.et al.Analysis of soil wetting and solute transport in subsurface trickle irrigation[J].Irrigation Science.2003，22(3/4)：143-156.

[8] 許 迪，程先軍. 地下滴灌土壤水運動和溶質運移數學模型的應用[J]. 農業工程學報，2002，18(1)：27-31.

[9] 李光永，鄭耀泉，曾德超，等. 地埋點源非飽和土壤水運動的數值模擬[J]. 水利學報，1996，(11)：47-56.

[10] 周青云，孫西歡，康紹忠. 蓄水坑灌條件下土壤水分運動的數值模擬[J]. 水利學報，2006，37(3)：342-353.

[11] 李光永，曾德超，段中鎖，等. 地埋點源滴灌土壤水運動規律的研究[J]. 農業工程學報，1996，12(3)：66-71.

[12] 李光永，曾德超，鄭耀泉. 地表點源滴灌土壤水運動的動力學模型與數值模擬[J]. 水利學報，1998，(11)：1-25.

[13] 李明思，康紹忠，孫海燕. 點源滴灌滴頭流量與濕潤體關系研究[J]. 農業工程學報，2006，22(4)： 32-35.

[14] Warrick A W. Soil water dynamics[M]. Oxford: Oxford University Press, 2003.

[15] Simunek J，Sejna M, van Genuchten M T. The HYDRUS-2D Software Package for Simulating the Two-Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Media[M]. California: U.S. salinity laboratory, Agricultural Research Service, U.S.Department of Agriculture, Riverside, 1999.