(遼寧水利土木工程咨詢有限公司，遼寧 沈陽 110003)

入滲(又稱為下滲)，是水從地表垂直或水平滲入到土壤和地下的運動過程，這是自然界中水文循環(huán)必不可少的一部分。Hydrus是一個可用來模擬土壤水流及溶質(zhì)三維運動的有限元計算機模型，它可以靈活處理各類水流邊界，包括定水頭和變水頭邊界、給定流量邊界、滲水邊界、大氣邊界以及排水溝等。

入滲是半干旱地區(qū)水分循環(huán)的重要環(huán)節(jié)，以內(nèi)蒙古東部平原區(qū)這一半干旱地區(qū)土壤在膜下滴灌方式下的降雨入滲過程為例，本文研究了Hydrus-2D模型在土壤水分入滲過程模擬中的應(yīng)用，以期為模型在該地區(qū)的應(yīng)用及推廣提供借鑒。

1 試驗區(qū)概況

試驗區(qū)位于內(nèi)蒙古東部平原區(qū)，在通遼市開魯縣，為西遼河流域，屬大陸性溫帶半干旱季風(fēng)性氣候，年平均氣溫6.80℃，年平均降雨量340.5mm，無霜期148d，試驗區(qū)屬于典型的北方半干旱區(qū)，水資源匱乏，供水主要以地下水為主，而用水主要以農(nóng)牧業(yè)為主，隨著農(nóng)業(yè)生產(chǎn)規(guī)模的不斷擴大和氣候的變化，致使該地區(qū)的水資源短缺現(xiàn)象日益嚴(yán)重。

2 試驗區(qū)土壤條件

試驗點位于東經(jīng)121°23′，北緯43°37′，試驗點2m 深度范圍內(nèi)土壤分為3層，其中第1層0～30cm，第2層30～60cm，第3層60～200cm。用環(huán)刀在各層取土樣，同時，用自封袋和鋁盒在對應(yīng)各土層取土樣，每層3個土樣重復(fù)操作。將用環(huán)刀取的土樣帶回實驗室測定土壤容重，將用鋁盒取的土樣帶回實驗室測定土壤飽和含水率，將自封袋中的土樣風(fēng)干、碾壓、篩分、均勻混合后用篩分法測定土壤粒徑(各測定結(jié)果見表1)。

表1 實驗區(qū)土壤基本物理性質(zhì)

3 基于Hydrus-2D模型的膜下滴灌方式下的降雨入滲規(guī)律研究

Hydrus是用于模擬飽和與非飽和多孔介質(zhì)中水分、能量、溶質(zhì)運移的一維、二維和三維有限元計算模型軟件。該模型的水流狀態(tài)為飽和-非飽和達(dá)西水流，水流控制方程采用Richards方程，溶質(zhì)運移方程采用對流-彌散(CED)方程，模型方程求解采用Galerkin線性有限元法。軟件可以靈活處理各類水流邊界，包括定水頭和變水頭邊界、給定流邊界、滲水邊界、自由排水邊界、大氣邊界等。另外，該軟件也包括一個參數(shù)優(yōu)化算法，可用于各種土壤的水壓和溶質(zhì)運移參數(shù)的逆向估計，為使用者提供了反推土壤特性經(jīng)驗參數(shù)的途徑，從而提高了擬合效率和擬合效果。

3.1 建立模型

3.1.1 水分運動基本方程及模型

(1) 降雨入滲模型方程

因為膜下滴灌和地面灌的降雨入滲過程可以看成土壤水分在水平方向和垂直方向的運移過程，因此，本文基于Hydrus-2D軟件，建立膜下滴灌和地面灌二維降雨入滲模型。根據(jù)Richards方程，并考慮植物根系吸水的二維土壤水分運動方程為：

(1)

式中θ——土壤體積含水率，cm3/cm3；

t——時間，h；

k(θ)——土壤非飽和導(dǎo)水率，cm/h；

φ——土壤總水勢，cm，對于非飽和土，φ=φm±z，其中φm為土壤基質(zhì)勢；

z——觀測點垂直坐標(biāo)，坐標(biāo)向上為正，向下為負(fù)；

S——匯項，表示根系吸水速率，cm/h。

(2) 土壤水分運動模型

Hydrus模擬涉及土壤水分特征曲線和土壤導(dǎo)水率曲線兩類參數(shù)，本文對這兩類參數(shù)的刻畫釆用vanGenuchten-Mualem模型(Van-Genuchten 1980)：

(2)

(3)

式中θ(x)——土壤物理特征曲線；

θs——飽和含水率，cm3/cm3；

θr——剩余飽和率，cm3/cm3；

h——壓力水頭，cm；

K(h)——壓力水頭下的導(dǎo)水率；

Ks——飽和導(dǎo)水率；

α、m、n和l——經(jīng)驗常數(shù)。其中m和n之間有一定的相關(guān)關(guān)系，一般認(rèn)為m=1-1/n時模型的簡化效果最好(Van-Genuchten 1980)。對于大多數(shù)土壤，l參數(shù)一般取0.5；

Se——土壤飽和度。

(3) 根系吸水模型

在根系吸水模塊，Hydrus軟件共提供兩種模型，F(xiàn)eddes模型和S-Shape模型，本次模擬采用Feddes模型，計算公式為:

S(h)=α(h)b(x)Tp

(4)

式中S(h)——根系吸水速率，即從單位體積土體單位時間內(nèi)吸取的水量；

α(h)——水分脅迫響應(yīng)函數(shù)；

b(x)——根系吸水分配密度函數(shù)；

Tp——潛在蒸騰量。

3.1.2 模型概化

(1) 時間和空間離散化

對試驗區(qū)進(jìn)行原位觀測的時間內(nèi)，7月份降雨量最大，為81.10mm，其中21日及25日兩天降雨量分別為35.90mm、25.40mm，基于本文的研究目標(biāo)，選取這兩天的降雨量進(jìn)行模擬計算。土壤模擬時期分為兩段：第一段自2016年7月20日至2016年7月24日，共120個小時；第二段從 2016年7月24日至 2016年7月29日，共120個小時，各組模擬的時間都以h為單位，每隔2h輸出一個計算結(jié)果，用以與實測值進(jìn)行比較。模擬土層為0～2m的垂直剖面。將一維的土壤剖面平均剖分為201個節(jié)點，節(jié)點間距為1cm。在上述土層剖分的基礎(chǔ)上，根據(jù)土壤剖面特征和實測的土壤質(zhì)地、顆粒組成、容重等數(shù)據(jù)，將土層分成若干介質(zhì)層，以便賦予不同的土壤水分運動參數(shù)(同一介質(zhì)層的土壤水分運動參數(shù)相同，各層之間的土壤水分運動參數(shù)不同)。進(jìn)行土壤水分動態(tài)模擬時將土壤剖面分為0～30cm、30～60cm、60～200cm3層。

(2) 作物根系概化

玉米根系在土壤中的分布與生長方式可分為水平根和垂直根。根據(jù)一些學(xué)者研究發(fā)現(xiàn)玉米根系在水平方向延展距離一般為0.8m，在豎直方向延展距離一般為1.25m。玉米在不同質(zhì)地的土壤中根系分布各不相同，但總體上玉米大部分根系主要分布在0～40cm土層內(nèi)并在生育后期深層土壤根系量有所增加。在本次模擬研究中設(shè)定地表以下0～40cm深度范圍內(nèi)為主要根系吸水區(qū)。

(3) 觀測點的設(shè)置

因為膜下滴灌存在著覆膜和起壟，而它們的存在影響了農(nóng)田土壤水的入滲、蒸發(fā)、運移等過程，因此為了更好地研究該灌溉方式下降雨入滲過程，可以把膜下滴灌分為3個土壤剖面進(jìn)行研究，分別為膜中剖面、膜邊剖面和溝中剖面，同時，在各個剖面設(shè)置觀測點，共設(shè)置30個觀測點，平均分布在膜中、膜邊和溝中3個剖面上，即每個斷面分配10個觀測點，其中在土壤深度0～40cm范圍內(nèi)每隔10cm設(shè)置1個觀測點，在40cm深度一下每隔20cm設(shè)置1個觀測點。

3.1.3 初始及邊界條件的設(shè)定

a.初始條件：Hydrus軟件提供了負(fù)壓水頭和土壤含水率兩種初始條件形式，本次模擬選取了土壤含水率初始條件形式。因為模擬時期分為兩段，所以在進(jìn)行模擬時以降雨前一天作為每一次模擬的初始時間，即 2016年7月 20日和 2016年7月24日，給定當(dāng)天的整個剖面各個位置的實測土壤含水率作為初始條件，相鄰測點間的含水率按線性插值程序自動給出。每一次模擬又分成膜下滴灌組和地面灌組并分別給出初始條件。

b.邊界條件：根據(jù)實際情況邊界條件設(shè)定如下：?上邊界，將膜下滴灌覆膜的部分設(shè)定為無水流邊界，而無膜部分土壤與大氣直接相通，既有降雨的入滲又受到蒸發(fā)作用，所以將其設(shè)定為大氣邊界。地面灌方式下的土壤與大氣直接相通，所以其上邊界也設(shè)定為大氣邊界；?下邊界，研究區(qū)域內(nèi)年平均地下水埋深為9m，大于所模擬的2m深度，故膜下滴灌和地面灌下邊界都設(shè)定為自由排水邊界。

c.大氣邊界為時變邊界條件，主要包括降雨量、蒸發(fā)量、蒸騰量等信息，其中降雨數(shù)據(jù)根據(jù)實測降雨資料輸入，蒸發(fā)狀況則由模型根據(jù)輸入的氣象條件及植被情況自行計算。

3.2 模型參數(shù)的確定

在此模型模擬中將土層劃分為3層進(jìn)行模擬。根據(jù)土壤的顆粒組成和容重，利用Hydrus軟件自帶的神經(jīng)網(wǎng)絡(luò)預(yù)測模塊，首先對土壤水分運動參數(shù)設(shè)定初始值，然后求解正問題，并將結(jié)果與實測值進(jìn)行比較，反饋并修改土壤水分運動參數(shù)，直至模擬值與實測值差距最小，這時得到的土壤水分運動參數(shù)認(rèn)為是最優(yōu)的結(jié)果(各參數(shù)最終結(jié)果見表2)。

表2 土壤水分運動參數(shù)

3.3 運行模型及擬合效果分析

為直觀評價模型的模擬效果，下面給出了模擬時間段內(nèi)各層土壤含水率實際觀測值和對應(yīng)模擬值的對比關(guān)系圖(如圖1所示)。

由圖1可看出：?在觀測深度下土壤含水率的實際值和模擬值的變化趨勢大體相同；?在土壤含水率較高時的模擬值與實測值偏差相對較小；?膜下滴灌三個剖面的模擬效果比膜下滴灌整體層面上的模擬效果相對較好。總的來說，模擬值能很好地反映出土壤含水率隨土層深度變化時的變化趨勢，所建模型對土壤含水率的模擬效果較好，可用于試驗區(qū)田間土壤水分運動的模擬。

3.4 模型應(yīng)用

降雨入滲過程主要受到降雨強度、初始含水率、下墊面方式、土壤特性等因素的影響。對于相同的研究對象其土壤特性一般保持不變，所以在研究同一區(qū)域降雨入滲時主要考慮另外三個影響因素。

3.4.1 模擬試驗方案設(shè)置

模擬試驗設(shè)置為兩因素多水平試驗，兩因素分別為：降雨強度、初始含水率。根據(jù)試驗區(qū)多年降雨資料，降雨強度分別設(shè)置為2mm/h、4mm/h，8mm/h3個水平；所觀測的土壤深度范圍內(nèi)共有3種土質(zhì)，共分3層，結(jié)合2016年試驗所測得的含水率數(shù)據(jù)，3層土壤初始含水率分別設(shè)置為對應(yīng)層土壤飽和含水率的40%、60%和80%。試驗因素水平見表3。

表3 試驗因素水平

3.4.2 模型模擬與求解

本次模擬試驗是在已建立好的膜下滴灌和地面灌模型基礎(chǔ)上進(jìn)行的，在模擬時僅對降雨條件和初始含水率條件進(jìn)行改變，其他信息不做改變。所模擬的一次累積降雨量為32mm，當(dāng)雨強為2mm/h時設(shè)置降雨歷時為16h，連續(xù)降雨，中間無間隔；當(dāng)雨強為4mm/h時設(shè)置降雨歷時為8h，連續(xù)降雨，中間無間隔；當(dāng)雨強為8mm/h時設(shè)置降雨歷時為4h，連續(xù)降雨，中間無間隔。在進(jìn)行各次模擬試驗時設(shè)置模擬總時長都為120h，初始時刻設(shè)置為降雨開始前4h。

模型計算中，初始步長為0.5h，最小步長為0.001h，最大步長為1h，每隔2h輸出一個計算結(jié)果。

3.4.3 模擬結(jié)果分析

因為降雨入滲基本上在雨后1—2d內(nèi)就能完成，所以在分析時取雨后48h輸出的模擬結(jié)果，通過對模擬結(jié)果進(jìn)行整理，繪制出不同降雨強度、不同初始含水率條件下，膜下滴灌在雨后48h時的土壤各層含水率狀況圖(如圖2所示)。

圖2 不同初始含水率條件下雨后48h時的土壤含水率分布圖

由圖2可知：?在降雨強度不變的情況下，初始含水率越低，雨后經(jīng)過一段時間其土壤含水率就越大，這種現(xiàn)象在淺層尤其明顯，原因是初始含水率越低其土壤入滲能力越強，在同等的降雨強度下，入滲量就越多，從而含水率也就越大；?在初始含水率不變的情況下降雨強度越低，降雨后土壤含水率在整體上就越大，并且這種現(xiàn)象在淺層最為明顯，由此可見土壤入滲能力隨著降雨強度的增加而降低。

為進(jìn)一步了解初始含水率和降雨強度對土壤入滲的影響，做出了模擬試驗雨后48h的累積入滲量分布圖(如圖3所示)。

圖3 不同雨強、不同初始含水率條件下雨后48h累積入滲量分布圖

由圖3可知：初始含水率和降雨強度對降雨累積入滲量都有影響，具體表現(xiàn)為：在降雨強度一定的情況下，初始含水率越高降雨累積入滲量越低；在初始含水率一定的情況下，降雨強度越大降雨累積入滲量相對來說就越低。

4 結(jié) 語

本文以開魯縣田間降雨入滲過程為例，介紹了 Hydrus-2D 軟件組成，基于室內(nèi)試驗取得的各項土壤物理特征參數(shù)，模擬了膜下滴灌方式下土壤含水量變化過程，并以部分實測數(shù)據(jù)對模擬的入滲過程進(jìn)行了驗證，從而實踐了Hydrus-2D 模型在土壤水入滲過程中的應(yīng)用，模型能很好地模擬出膜下滴灌方式下降雨入滲過程，可用于實驗區(qū)田間土壤水分運動的模擬計算。

試驗結(jié)果表明：在降雨強度不變情況下，初始含水率越低，雨后經(jīng)過一段時間土壤含水率就越大；在初始含水率不變的情況下降雨強度越低，降雨后土壤含水率在整體上就越大；在降雨強度一定的情況下初始含水率越高，降雨累積入滲量越低；在初始含水率一定的情況下降雨強度越大，降雨累積入滲量相對來說就越低。