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

劉彬彬，劉堯兵，張科鋒

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

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

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

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

1 模型原理

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

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

(1)

式中,θ是單位體積土壤含水量(cm3·cm-3)；t為時間(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的大小直接決定著對土壤水分特征曲線的坡度影響；Se為有效含水率；l為經(jīng)驗擬合參數(shù),一般為0.5。

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

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

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

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

1.4 數(shù)值實驗方案

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

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

2 結(jié)果與分析

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

2.1 灌溉后土壤含水量實測值和模擬值之間的對比

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

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

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

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

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

為便于表示，將模擬方案統(tǒng)一命名為RxTy，其中，x表示灌溉速率(cm·h-1)，y表示灌溉時間(h)。就淤泥、粉砂壤土、砂質(zhì)粘壤土三種不同土質(zhì)，采取了多種灌溉速度和時間的組合進(jìn)行數(shù)值實驗。對于淤泥，采取以下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；對于粉砂壤土，9種灌溉組合方案分別是：R0.7T1、R0.7T1.5、R0.7T2、R0.7T2.5、R0.7T3、R0.8T1、R0.8T1.5、R0.9T1、R1.2T0.5；對于砂質(zhì)粘壤土，灌溉組合方案設(shè)計為：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種不同的灌溉速度和時間組合，模擬結(jié)果較多。作為比較分析，只選取每種土質(zhì)的兩種灌溉方案(灌溉量較小和較大的各一種)，并予以詳細(xì)地分析在灌溉過程中和灌溉結(jié)束后土壤含水量分布的動態(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可以很明顯地看出，不同的模擬方案中的上層濕潤土在灌溉結(jié)束時均接近飽和狀態(tài)，隨著深度的增加，從接近飽和區(qū)的底部過渡到初始含水量區(qū)，含水量驟降，以粉砂壤土R0.7T3方案為例，土壤剖面含水量從0.432 cm3·cm-3下降至初始含水量0.172 cm3·cm-3，土壤深度變化僅為4 cm。另一方面，由于水分重分布的作用，灌溉結(jié)束24 h后和灌溉結(jié)束時的土壤剖面含水量有很大的差異，在灌溉結(jié)束時，灌溉水儲存在距表層9.2 cm深的土層中，但經(jīng)過24 h后，濕潤鋒向下滲透至20.6 cm，約有55.3%的灌溉水參加了水分重分布。從圖2中還可發(fā)現(xiàn)，不管何種土質(zhì)，灌溉結(jié)束后土壤水分重分布現(xiàn)象明顯，重分布后的土壤含水量更趨均勻，尤其是對于排水性能較好的土質(zhì)，土壤水分重分布尤為明顯。因此灌溉時必須考慮土壤水分重分布的情況，以免造成水資源浪費(fèi)。

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

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

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

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

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

圖4 土壤濕潤深度隨時間的變化Fig. 4 Variation of soil wetting depth with time

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

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

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

3 結(jié) 論

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

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

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

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

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