基于Hydrus-2D模型的膜下滴灌對降雨入滲影響研究

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

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

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

1 試驗區概況

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

2 試驗區土壤條件

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

表1 實驗區土壤基本物理性質

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

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

3.1 建立模型

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

(1) 降雨入滲模型方程

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

(1)

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

t——時間，h；

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

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

z——觀測點垂直坐標，坐標向上為正，向下為負；

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

(2) 土壤水分運動模型

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

(2)

(3)

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

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

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

h——壓力水頭，cm；

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

Ks——飽和導水率；

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

Se——土壤飽和度。

(3) 根系吸水模型

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

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

(4)

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

α(h)——水分脅迫響應函數；

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

Tp——潛在蒸騰量。

3.1.2 模型概化

(1) 時間和空間離散化

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

(2) 作物根系概化

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

(3) 觀測點的設置

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

3.1.3 初始及邊界條件的設定

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

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

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

3.2 模型參數的確定

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

表2 土壤水分運動參數

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

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

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

3.4 模型應用

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

3.4.1 模擬試驗方案設置

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

表3 試驗因素水平

3.4.2 模型模擬與求解

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

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

3.4.3 模擬結果分析

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

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

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

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

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

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

4 結 語

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

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