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

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

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

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

1 試驗(yàn)區(qū)概況

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

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

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

表1 實(shí)驗(yàn)區(qū)土壤基本物理性質(zhì)

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

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

3.1 建立模型

3.1.1 水分運(yùn)動(dòng)基本方程及模型

(1) 降雨入滲模型方程

因?yàn)槟は碌喂嗪偷孛婀嗟慕涤耆霛B過程可以看成土壤水分在水平方向和垂直方向的運(yùn)移過程，因此，本文基于Hydrus-2D軟件，建立膜下滴灌和地面灌二維降雨入滲模型。根據(jù)Richards方程，并考慮植物根系吸水的二維土壤水分運(yùn)動(dòng)方程為：

(1)

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

t——時(shí)間，h；

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

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

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

S——匯項(xiàng)，表示根系吸水速率，cm/h。

(2) 土壤水分運(yùn)動(dòng)模型

Hydrus模擬涉及土壤水分特征曲線和土壤導(dǎo)水率曲線兩類參數(shù)，本文對(duì)這兩類參數(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)驗(yàn)常數(shù)。其中m和n之間有一定的相關(guān)關(guān)系，一般認(rèn)為m=1-1/n時(shí)模型的簡化效果最好(Van-Genuchten 1980)。對(duì)于大多數(shù)土壤，l參數(shù)一般取0.5；

Se——土壤飽和度。

(3) 根系吸水模型

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

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

(4)

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

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

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

Tp——潛在蒸騰量。

3.1.2 模型概化

(1) 時(shí)間和空間離散化

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

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

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

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

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

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

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

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

表2 土壤水分運(yùn)動(dòng)參數(shù)

3.3 運(yùn)行模型及擬合效果分析

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

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

3.4 模型應(yīng)用

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

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

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

表3 試驗(yàn)因素水平

3.4.2 模型模擬與求解

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

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

3.4.3 模擬結(jié)果分析

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

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

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

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

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

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

4 結(jié) 語

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

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