基于HYDRUS-1D的不同質(zhì)地土壤入滲過(guò)程數(shù)值模擬

劉彬彬，劉堯兵，張科鋒

(1. 浙江大學(xué)建筑工程學(xué)院，浙江 杭州 310058; 2. 浙江大學(xué)寧波理工學(xué)院，浙江 寧波 315100)

世界上40%的食物是由灌溉農(nóng)業(yè)生產(chǎn)的，消耗全世界70%的淡水資源[1]。然而，目前我國(guó)的農(nóng)業(yè)灌溉效率為 45%,僅為發(fā)達(dá)國(guó)家的一半左右, 農(nóng)業(yè)節(jié)水潛力巨大[2-3]。農(nóng)業(yè)用水效率的提高對(duì)緩解我國(guó)水資源短缺現(xiàn)狀, 解決供水危機(jī)以及保護(hù)生態(tài)環(huán)境具有重要作用。

優(yōu)化灌溉速率和時(shí)間對(duì)于節(jié)省農(nóng)業(yè)用水至關(guān)重要。然而，制定灌溉方案時(shí)涉及的因素眾多，包括：土壤質(zhì)地及與之相關(guān)的土壤水力特性、作物的根深、根區(qū)土壤的含水量以及天氣情況等。目前，在優(yōu)化灌溉的研究方面，土壤水分傳感器的應(yīng)用已日益普遍[4-7]，它可以實(shí)時(shí)監(jiān)測(cè)土壤剖面不同深度的含水量或水勢(shì)，有助于做出精確的灌溉策略。但由于埋設(shè)的土壤水分傳感器數(shù)目有限，它們不能精確地確定灌溉后土壤剖面水分分布和最終的潤(rùn)濕深度。此外，傳感器的維護(hù)和觀測(cè)均需要一定的技術(shù)含量,只適用于土地規(guī)模經(jīng)營(yíng)的農(nóng)場(chǎng)。隨著土壤科學(xué)和計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬方法越來(lái)越多地應(yīng)用于土壤水分運(yùn)動(dòng)的研究中，已成為定量分析土壤水分入滲特性的有效工具[8-10]。美國(guó)國(guó)家鹽分實(shí)驗(yàn)室(US-Salinity laboratory)研發(fā)的 HYDRUS-1D 軟件，可對(duì)一維飽和-非飽和土壤中水、熱及溶質(zhì)運(yùn)移進(jìn)行數(shù)值模擬[11]。國(guó)內(nèi)外已有諸多文獻(xiàn)報(bào)道使用HYDRUS軟件開展在不同灌溉模式的水分運(yùn)動(dòng)和作物根系吸水等方面的研究，取得了可喜的成果[12-15]。但研究?jī)H限于特定的條件，尚未見使用該軟件對(duì)不同土壤的濕潤(rùn)深度和灌溉量之間關(guān)系的系統(tǒng)研究，對(duì)指導(dǎo)優(yōu)化灌溉、節(jié)約農(nóng)業(yè)用水局限性明顯。

本研究基于HYDRUS-1D軟件，針對(duì)不同的土質(zhì)，在不產(chǎn)生地表徑流的情況下，系統(tǒng)地模擬灌溉結(jié)束時(shí)和灌溉結(jié)束24 h之后的土壤含水量分布，期望能獲得在不同土壤條件下灌溉水分重分布后土壤濕潤(rùn)深度和灌溉量之間的定量關(guān)系，為不同根深的作物制定合理的灌溉方案提供參考。

1 模型原理

1.1 土壤水分運(yùn)動(dòng)方程

假定土壤為均質(zhì)、各向同性的多孔介質(zhì)，忽略溫度與土壤中的氣相對(duì)土壤水分的影響，不考慮根系吸水與源匯項(xiàng)，用以研究飽和-非飽和流動(dòng)問題的一維Richards控制方程[16]：

(1)

式中,θ是單位體積土壤含水量(cm3·cm-3)；t為時(shí)間(h)；z是垂直空間坐標(biāo)，也可認(rèn)為是土壤深度(cm)；K是土壤導(dǎo)水率(cm·h-1)；h是土壤壓力水頭(cm)。

1.2 土壤水力特征參數(shù)方程

土壤水力特征參數(shù)的描述采用van-Genuchten模型，其表達(dá)形式為[17]：

(2)

(3)

(4)

式中,θs是土壤飽和含水率(cm3·cm-3)；θr是土壤殘余含水率(cm3·cm-3)；Ks是飽和導(dǎo)水率(cm·h-1)；α為與進(jìn)氣吸力相關(guān)的參數(shù)(1·cm-1)；n和m是形狀系數(shù)；n的大小直接決定著對(duì)土壤水分特征曲線的坡度影響；Se為有效含水率；l為經(jīng)驗(yàn)擬合參數(shù),一般為0.5。

1.3 模型邊界條件和初始條件

本模擬是為了研究水分重分布后土壤濕潤(rùn)深度和灌溉水量之間的關(guān)系。對(duì)淤泥、粉砂壤土、砂質(zhì)粘壤土三種土質(zhì)進(jìn)行了以下的系統(tǒng)模擬：在不產(chǎn)生徑流量的情況下，設(shè)計(jì)在3 h內(nèi)不同的灌溉速率和時(shí)間的組合，得到模擬土層在灌溉結(jié)束時(shí)和灌溉結(jié)束24 h后的水分分布情況，并確定灌溉水分重分布后的土壤濕潤(rùn)深度。

模擬區(qū)域?yàn)?00 cm厚的土壤深度，等分為500層，每層厚度為0.2 cm。模型區(qū)域上邊界為大氣邊界，下邊界為自由排水邊界。由于模擬的持續(xù)時(shí)間較短，故忽略土壤表面蒸發(fā)，同時(shí)也不考慮根系吸水作用。

計(jì)算區(qū)域的初始含水量θ0為土壤有效含水量的中間值，即為田間持水量θFC和永久枯萎點(diǎn)含水量θPWP的中間值[18]。對(duì)大部分作物而言，當(dāng)根區(qū)土壤含水量低于θ0時(shí)，作物無(wú)法獲取生長(zhǎng)所需的潛在蒸騰量，因而影響作物生長(zhǎng)[18]。

1.4 數(shù)值實(shí)驗(yàn)方案

本數(shù)值實(shí)驗(yàn)采用的土壤水力特征參數(shù)如表1所示[19]，同時(shí)表1還列出了田間持水量θFC和永久萎蔫點(diǎn)含水量θPWP。

表1 土壤的水力特性參數(shù)Table 1 Soil hydraulic characteristics

2 結(jié)果與分析

本研究的結(jié)果分兩部分，第一部分是驗(yàn)證灌溉后土壤含水量的實(shí)測(cè)值和模擬值之間的吻合度，以確定用本文提出的數(shù)學(xué)方法能否較準(zhǔn)確地模擬實(shí)驗(yàn)結(jié)果；另一部分是通過(guò)對(duì)淤泥、粉砂壤土、砂質(zhì)粘壤土三種土質(zhì)的系統(tǒng)模擬，量化不同土質(zhì)的土壤水分分布與灌溉水量和時(shí)間之間的關(guān)系。

2.1 灌溉后土壤含水量實(shí)測(cè)值和模擬值之間的對(duì)比

本驗(yàn)證示例取自文獻(xiàn)[20]。灌溉實(shí)驗(yàn)在直徑5 cm、長(zhǎng)160 cm的柱體里進(jìn)行，土質(zhì)為細(xì)砂壤土，密度為1.47 g·cm-3，灌溉量為5 cm和10 cm兩種，觀測(cè)灌溉后不同時(shí)間的柱內(nèi)土壤含水量分布，土壤的水力特征參數(shù)和土壤含水量觀測(cè)方法見文獻(xiàn)[20]。本研究取10 cm灌溉量、灌溉1天后的土壤含水量分布作為驗(yàn)證示例。圖1顯示了土壤含水量實(shí)測(cè)值和模擬值之間的對(duì)比，不難看出土壤含水量實(shí)測(cè)值和模擬值之間吻合良好，模擬值和實(shí)測(cè)值之間的R2值達(dá)0.9以上，表明兩者之間高度相關(guān)。

圖1 灌溉后土壤含水量實(shí)測(cè)值和模擬值之間的對(duì)比Fig. 1 Comparison of soil water content between measurement and simulation

為客觀評(píng)價(jià)模擬結(jié)果，對(duì)模擬的土壤含水量值進(jìn)行了統(tǒng)計(jì)分析，計(jì)算的統(tǒng)計(jì)指標(biāo)包括：均方根誤差(RMSE)、相對(duì)均方根誤差(RRMSE)、納什系數(shù)[21](NSE)和平均誤差(AE)，統(tǒng)計(jì)分析結(jié)果見表2。從表2中可以發(fā)現(xiàn)，RMSE、RRMSE和AE的值均低，說(shuō)明模擬值和實(shí)測(cè)值接近；而NSE的值較高，一般認(rèn)為若NSE的值超過(guò)0.75，模型可以被判定為“非常好”。這些統(tǒng)計(jì)數(shù)據(jù)表明了采用本文的數(shù)值方法能滿意地模擬灌溉后土壤含水量分布。大量的文獻(xiàn)也表明HYDRUS-1D模型能較好預(yù)測(cè)土壤含水量的變化[22-24]。

表2 模擬結(jié)果的統(tǒng)計(jì)分析Table 2 Statistical analysis of simulated results

2.2 不同土質(zhì)的灌溉模擬分析

為便于表示，將模擬方案統(tǒng)一命名為RxTy，其中，x表示灌溉速率(cm·h-1)，y表示灌溉時(shí)間(h)。就淤泥、粉砂壤土、砂質(zhì)粘壤土三種不同土質(zhì)，采取了多種灌溉速度和時(shí)間的組合進(jìn)行數(shù)值實(shí)驗(yàn)。對(duì)于淤泥，采取以下9種灌溉組合方案： R0.7T0.5、R0.9T0.5、R0.5T1、R0.6T1、R0.5T1.5、R0.6T1.5、R0.5T2、R0.5T2.3、R0.5T2.6；對(duì)于粉砂壤土，9種灌溉組合方案分別是：R0.7T1、R0.7T1.5、R0.7T2、R0.7T2.5、R0.7T3、R0.8T1、R0.8T1.5、R0.9T1、R1.2T0.5；對(duì)于砂質(zhì)粘壤土，灌溉組合方案設(shè)計(jì)為：R0.9T1、R1T1.5、R1T2、R1T2.5、R1.1T2.8、R1.25T1、R1.25T1.5、R1.25T1.8、R1.45T0.5。

2.2.1 土壤含水量重分布 如上所述，每種土質(zhì)有9種不同的灌溉速度和時(shí)間組合，模擬結(jié)果較多。作為比較分析，只選取每種土質(zhì)的兩種灌溉方案(灌溉量較小和較大的各一種)，并予以詳細(xì)地分析在灌溉過(guò)程中和灌溉結(jié)束后土壤含水量分布的動(dòng)態(tài)變化，這六種方案分別是：淤泥R0.9T1、0.4T1.5、粉砂壤土R0.9T1、粉砂壤土R0.7T3、砂質(zhì)粘壤土R0.9T1和砂質(zhì)粘壤土R1.1T2.8(圖2)，灌溉量的變化范圍0.6～3.08 cm。

從圖2可以很明顯地看出，不同的模擬方案中的上層濕潤(rùn)土在灌溉結(jié)束時(shí)均接近飽和狀態(tài)，隨著深度的增加，從接近飽和區(qū)的底部過(guò)渡到初始含水量區(qū)，含水量驟降，以粉砂壤土R0.7T3方案為例，土壤剖面含水量從0.432 cm3·cm-3下降至初始含水量0.172 cm3·cm-3，土壤深度變化僅為4 cm。另一方面，由于水分重分布的作用，灌溉結(jié)束24 h后和灌溉結(jié)束時(shí)的土壤剖面含水量有很大的差異，在灌溉結(jié)束時(shí)，灌溉水儲(chǔ)存在距表層9.2 cm深的土層中，但經(jīng)過(guò)24 h后，濕潤(rùn)鋒向下滲透至20.6 cm，約有55.3%的灌溉水參加了水分重分布。從圖2中還可發(fā)現(xiàn)，不管何種土質(zhì)，灌溉結(jié)束后土壤水分重分布現(xiàn)象明顯，重分布后的土壤含水量更趨均勻，尤其是對(duì)于排水性能較好的土質(zhì)，土壤水分重分布尤為明顯。因此灌溉時(shí)必須考慮土壤水分重分布的情況，以免造成水資源浪費(fèi)。

圖3顯示了地表以下5 cm處土壤含水量隨時(shí)間的變化，灌溉期間土壤含水量隨時(shí)間的增加急劇增加，灌溉結(jié)束后由于水分重分布的作用，含水量隨時(shí)間的增長(zhǎng)而減少，初期的減少速率大于后期，24 h后基本趨于穩(wěn)定。以粉砂壤土R0.7T3為例，在3～10 h期間土壤含水量從0.431 cm3·cm-3下降到0.346 cm3·cm-3，而在其后的14 h內(nèi)，含水量?jī)H下降了0.03 cm3·cm-3。對(duì)于同一種土質(zhì)，灌溉量小的土壤含水量變化較灌溉量大的更為平緩。

此外，還發(fā)現(xiàn)對(duì)于同一種土質(zhì)，當(dāng)灌溉量相同時(shí)，不同的灌溉速度在灌溉結(jié)束24 h之后土壤剖面的含水量分布相同。這說(shuō)明了重分布后的土壤含水量只與灌溉量有關(guān)，與灌溉速度和時(shí)間關(guān)系不大。

圖2 不同質(zhì)地土壤含水量的動(dòng)態(tài)分布Fig. 2 Dynamic distribution of soil water content in various soils

圖3 土壤含水量隨時(shí)間的變化Fig. 3 Variation of soil water content with time

2.2.2 土壤濕潤(rùn)深度與時(shí)間和灌溉量的關(guān)系 土壤濕潤(rùn)深度與時(shí)間的變化如圖4所示。不難發(fā)現(xiàn)，不同土質(zhì)間的土壤濕潤(rùn)深度隨時(shí)間的變化規(guī)律類似。灌溉期間土壤的濕潤(rùn)深度隨時(shí)間的增長(zhǎng)而迅速增加，且兩者基本呈線性關(guān)系，隨后在土壤水分重力和毛細(xì)孔的作用下，土壤的濕潤(rùn)深度隨時(shí)間的增加而繼續(xù)增加，但增速下降。以砂質(zhì)粘壤土R1.1T2.8為例，在灌溉期間，土壤的濕潤(rùn)深度增加量為14.2 cm，灌溉結(jié)束后的前5 h，土壤的濕潤(rùn)深度增加了5.6 cm，但在其后的16 h內(nèi)，土壤的濕潤(rùn)深度只增加了5.2 cm，尤其在20～24 h期間，濕潤(rùn)深度僅增加了1.0 cm，說(shuō)明土壤水分重分布在24 h內(nèi)基本完成。對(duì)于同一種土質(zhì)，在灌溉量小的情況下，土壤含水量的重分布更為平緩。

圖4 土壤濕潤(rùn)深度隨時(shí)間的變化Fig. 4 Variation of soil wetting depth with time

圖5顯示了三種土質(zhì)模擬得到的在灌溉結(jié)束24 h后灌溉總量和土壤濕潤(rùn)深度的關(guān)系。對(duì)于不同的土質(zhì)，其濕潤(rùn)深度與總灌溉量有顯著的線性關(guān)系，擬合優(yōu)度R2均在0.99以上，擬合的斜率相差不明顯，介于5.15(淤泥)和5.95(砂質(zhì)粘壤土)之間(見表3)。表3給出了灌溉總量和土壤濕潤(rùn)深度的簡(jiǎn)單定量關(guān)系，該關(guān)系對(duì)作物生產(chǎn)中的優(yōu)化灌溉具有實(shí)際指導(dǎo)意義。已有研究表明，作物從根深以外的土層中獲取水分的能力較弱，因此優(yōu)化灌溉的目標(biāo)之一就是要避免過(guò)量灌溉導(dǎo)致灌溉水在重分布后流出根區(qū)范圍。對(duì)于不同的作物，已有大量有關(guān)作物根系分布方面的研究[25-26]，作物的最大根深可以估計(jì)，其根深在生長(zhǎng)過(guò)程中與累計(jì)氣溫呈線性關(guān)系[27-28]，這就允許在作物生長(zhǎng)階段動(dòng)態(tài)地估計(jì)作物的根深。在作物根深已知的條件下，通過(guò)借助表3中的方程，可以推求出濕潤(rùn)根區(qū)土壤所需的水量，為優(yōu)化灌溉提供科學(xué)依據(jù)。

圖5 不同質(zhì)地土壤灌溉結(jié)束24 h后的土壤濕潤(rùn)深度與總灌溉量的關(guān)系Fig. 5 Relationship between soil wetting depth and total irrigation amount after 24 hours of irrigation

表3 土壤濕潤(rùn)深度與灌溉量間的擬合方程Table 3 Fitted equations relating soil wetting depth and irrigation amount

3 結(jié) 論

本文運(yùn)用HYDRUS-1D軟件，采用不同的灌溉速率和灌溉時(shí)間組合，對(duì)灌溉后不同土質(zhì)的土壤的土水動(dòng)力學(xué)特性進(jìn)行了系統(tǒng)的數(shù)值模擬試驗(yàn)，量化了在灌溉結(jié)束24 h后土壤濕潤(rùn)深度與灌溉總量的關(guān)系。基于以上結(jié)果，得出了以下幾點(diǎn)結(jié)論：

1)土壤的含水量分布和濕潤(rùn)深度只與土壤種類和灌溉量有關(guān)。對(duì)同一種土壤而言，當(dāng)灌溉量相等時(shí)，灌溉結(jié)束24 h后土壤含水量分布相同。

2)在灌溉結(jié)束后，上層土壤含水量接近飽和，灌溉水儲(chǔ)存在上層土壤中。水分重分布后，土壤透水性不同，土壤濕潤(rùn)深度增加不同。因此，對(duì)透水性好的土壤在制定灌溉方案時(shí)須考慮水分重分布的影響。

3)水分重分布后土壤的濕潤(rùn)深度與灌溉量之間存在著顯著的線性關(guān)系。該關(guān)系可根據(jù)灌溉量估計(jì)土壤濕潤(rùn)深度，避免灌溉水量過(guò)多而產(chǎn)生的水量損失。

本文的研究結(jié)果可為不同根深的作物制定合理的灌溉方案提供參考，對(duì)作物生產(chǎn)中的水分優(yōu)化利用有一定的指導(dǎo)意義。