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

張巧利，張 亮

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

0 引 言

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

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

1 涌流灌溉條件下土壤水分運動的數(shù)學(xué)模型

1.1 模型條件假定

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

1.2 土壤水分運動方程

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

(1)

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

1.3 初始及邊界

(1)初始條件[14]。

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

(2)

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

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

圖1 計算區(qū)示意圖

則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′邊界為上邊界，豎直方向水流通量為降雨強度或蒸發(fā)強度q0：

(6)

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

1.4 模型求解

1.4.1 HYDRUS-3D

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

1.4.2 模型求解

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

圖2 模擬計算流區(qū)的網(wǎng)格劃分

2 試驗驗證

2.1 土壤水動力學(xué)特性

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

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

(8)

其中：

(9)

式中：以上方程中有5個獨立參數(shù)θs,θr,α,n,K，其中：θs為土壤飽和含水率， cm3/cm3；θr為土壤殘余含水率，cm3/cm3；m、n、α為模型的擬合參數(shù)，其中：m=1-1/n，α的取值與土壤物理性質(zhì)有關(guān)，cm-1；Ks為土壤滲透系數(shù)，cm/min；Se為有效含水量，cm3/cm3；L為孔隙連通性參數(shù)，對大多數(shù)土壤一般取為0.5。供試土壤水分特性參數(shù)為：θ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 測試方法

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

2.3 結(jié)果驗證

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

圖3 觀測點位置圖

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

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

3 模型的應(yīng)用

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

3.1 土壤濕潤體內(nèi)水分分布狀況

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

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

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

3.2 土壤濕潤體分布

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

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

4 結(jié)論與討論

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

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

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

[1] 吳普特，朱德蘭.涌泉根灌技術(shù)研究與應(yīng)用[J].灌排機械工程學(xué)報，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] 趙偉霞，蔡煥杰，陳新明，等.無壓灌溉土壤濕潤體含水率分布規(guī)律與模擬模型研究[J].農(nóng)業(yè)工程學(xué)報，2007，23(3)：7-12.

[6] 康銀紅，馬孝義，李 娟，等.黃土高原重力式地下滴灌水分運動模型與分區(qū)參數(shù)研究[J].農(nóng)業(yè)機械學(xué)報，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] 許 迪，程先軍. 地下滴灌土壤水運動和溶質(zhì)運移數(shù)學(xué)模型的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報，2002，18(1)：27-31.

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

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

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

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

[13] 李明思，康紹忠，孫海燕. 點源滴灌滴頭流量與濕潤體關(guān)系研究[J]. 農(nóng)業(yè)工程學(xué)報，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.