地表滴灌條件下土壤濕潤體運移量化表征

陳 帥 毛曉敏

(中國農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院， 北京 100083)

0 引言

水資源短缺是干旱區(qū)農(nóng)業(yè)可持續(xù)發(fā)展的主要制約因素之一[1]。滴灌作為一種現(xiàn)代化高效節(jié)水灌溉技術(shù)，能根據(jù)作物需要對作物根區(qū)水分進行適時、適量的補充，從而有效提高水分利用效率[2]。地表滴灌條件下，在滴頭下方形成類似球狀的濕潤體，滴灌濕潤體的大小及水分分布特征會影響作物根系吸水及其伸展特性[3-4]，對作物生長和水分利用效率尤為重要。

國內(nèi)外學(xué)者對滴灌條件下土壤水分運動規(guī)律開展了諸多試驗研究。朱德蘭等[5]對粉質(zhì)粘土和重粉質(zhì)壤土中滴灌水分分布進行了試驗研究，提出了點源入滲中最大橫向、縱向濕潤直徑與滴頭流量、灌水量關(guān)系的數(shù)學(xué)模型；張振華等[6- 7]用土箱模擬研究了粘壤土點源入滲條件下土壤濕潤體特征值隨滴頭流量、土壤初始含水率和容重的變化規(guī)律，提出了地表滴灌入滲土壤濕潤體特征值的經(jīng)驗?zāi)Ｐ?；LI等[8]在砂土和壤土中研究滴頭流量對土壤水分分布的影響，建立了水平和垂直濕潤距離與灌水量的經(jīng)驗關(guān)系。這些研究成果簡便、實用，但多是針對特定的土壤質(zhì)地試驗得到的結(jié)果，由于土壤類型對滴灌濕潤體運移特征具有重要影響[9-10]，這些成果在其他土質(zhì)類型較難推廣應(yīng)用。

近年來，數(shù)值模擬方法在滴灌水分運動的研究中應(yīng)用較多。例如COTE等[11]運用HYDRUS-2D軟件定性分析了土壤滲透特性對地下滴灌中濕潤體形狀的影響；李光永等[12]提出了計算地表點源入滲土壤濕潤體特征值的通用數(shù)值方法，可用于常見土壤中滴灌濕潤體體積和入滲距離的計算，但需要量化土壤水分運動的水力特性參數(shù)；胡和平等[9]則在滴灌濕潤體運移的數(shù)值模擬中僅采用飽和導(dǎo)水率Ks來表征不同土質(zhì)對水分運動的影響。數(shù)值模擬對濕潤體的動態(tài)變化描述詳細，且在已知土壤水力特性參數(shù)的前提下具有很好的適用性。但也存在模擬過程復(fù)雜和模擬參數(shù)難以獲得等問題。

本文在搜集大量土質(zhì)水力特性參數(shù)的基礎(chǔ)上，采用HYDRUS-2D/3D進行不同質(zhì)地土壤中滴灌濕潤體運移特征的數(shù)值模擬。通過對數(shù)值模擬結(jié)果的分析，進一步構(gòu)建描述不同質(zhì)地土壤中濕潤體形狀特征的簡易經(jīng)驗公式，并采用針對性試驗結(jié)果對經(jīng)驗公式在不同土質(zhì)下預(yù)測結(jié)果的可靠性進行驗證。

1 模擬基本原理和模擬軟件

1.1 滴灌模擬基本原理

地表滴灌條件下土壤水分運動是三維流動，假設(shè)土壤各向同性，固相骨架不可壓縮，忽略地表土壤蒸發(fā)，單點源入滲條件下的土壤水分運動過程可用軸對稱二維Richards方程來描述[13]，即

(1)

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

h——土壤基質(zhì)勢或壓力水頭，cm

Kh——土壤導(dǎo)水率，cm/h

t——時間，h

r——徑向坐標(biāo)，cm

z——垂向坐標(biāo)，cm，規(guī)定z向上為正

土壤含水率、基質(zhì)勢和土壤非飽和導(dǎo)水率的關(guān)系采用van Genuchten模型來描述[14]，即

(2)

(3)

(4)

m=1-1/n

式中θr——土壤殘余含水率，cm3/cm3

θs——土壤飽和含水率，cm3/cm3

Ks——土壤飽和導(dǎo)水率，cm/h

Se——相對飽和度

α——土壤進氣值，cm-1

n、m——土壤形狀參數(shù)

l——孔隙關(guān)聯(lián)度參數(shù)，一般取0.5

單點源入滲的模擬求解區(qū)域是以AC為中心軸、ABCD為旋轉(zhuǎn)面的圓柱體，C點為滴頭位置(圖1)。

圖1 滴灌模擬計算區(qū)域示意圖 Fig.1 Sketch of flow domain in simulation of surface drip irrigation

CD為上邊界，灌水開始后，水分以較快的速度在土壤表面擴散，并在滴頭周圍形成土壤進水區(qū)域，進水區(qū)域先增大后趨于穩(wěn)定，BHATNAGAR等[15]研究表明，在滴灌模擬中采用穩(wěn)定的進水區(qū)域半徑對長時間的入滲結(jié)果影響較小，因此本研究采用穩(wěn)定進水區(qū)域半徑Rs[16-17]，在上邊界

(5)

式中q0——上邊界入滲流量，cm/h

不同土質(zhì)中Rs可近似表示為

(6)

式中q——滴灌流量，cm3/h

AC和BD分別為左右邊界，由于模擬區(qū)域的對稱性及模擬空間足夠大，左右邊界采用不透水邊界

(7)

下邊界采用自由排水邊界

(8)

假定模擬區(qū)域內(nèi)的土壤含水率均勻分布，因此初始條件

θz,r=θ0(0≤z≤Z, 0≤r≤R0,t=0)

(9)

1.2 數(shù)值模擬軟件與模擬情景設(shè)置

HYDRUS-2D/3D是用于模擬飽和-非飽和土壤中二維/三維水、熱、溶質(zhì)運移的模擬軟件[18]。其軸對稱垂向二維模塊可用于模擬點源入滲條件下土壤水分的運動和分布[19]。

模擬情景中選取不同質(zhì)地的土壤作為地表滴灌的影響因素。美國制土壤質(zhì)地分類系統(tǒng)根據(jù)土壤的機械組成將土壤劃分為12種不同的質(zhì)地，本研究以其中11種不同質(zhì)地的土壤作為滴灌模擬介質(zhì)(不包括粘土)，通過選取砂粒、粉粒和黏粒質(zhì)量分?jǐn)?shù)組合來確定代表土質(zhì)的樣本點，共選取33個樣本點(表1)。

如上所述研究中采用了van Genuchten函數(shù)來描述土壤水力特性。土壤水力特性參數(shù)受土壤顆粒組成和容重等因素的影響，直接通過試驗測定不同質(zhì)地土壤的水力特性參數(shù)是復(fù)雜的，常用做法是利用根據(jù)大量土壤數(shù)據(jù)構(gòu)建的轉(zhuǎn)換函數(shù)對土壤水力參數(shù)進行預(yù)測[20]。本研究根據(jù)33個樣本點的土壤機械組成(表1)、并假設(shè)土壤容重控制在1.5 g/cm3，利用HYDRUS軟件中的Rosetta Lite (v. 1.1)進行反演，反演得到的土壤水力特性參數(shù)見表1。SKAGGS等[21]研究表明Rosetta能夠?qū)Φ喂嗄M中的土壤水力特性參數(shù)進行快捷、可靠的預(yù)測。

表1 土壤機械組成及其水力特性參數(shù) Tab.1 Soil particle composition and parameters of soil hydraulic property

滴灌流量采用常用的1、2、3 L/h 3種。模擬中土壤初始含水率的取值略高于土壤殘余含水率，即0.1 cm3/cm3，這與實際農(nóng)田灌水前表層土壤的含水率接近。模擬區(qū)域徑向半徑R0為200 cm，垂向高度Z為200 cm，采用三角形網(wǎng)格對計算區(qū)域進行剖分，在滴頭位置對網(wǎng)格進行加密。模擬時段為1 500 min，采用變步長計算。

2 結(jié)果與討論

2.1 濕潤體運移規(guī)律

對11種不同的土質(zhì)(33個樣本點)，在不同滴頭流量(1、2、3 L/h)和不同灌水時刻(120、300、540、780、1 020、1 260、1 500 min)的濕潤體進行模擬，共得到693組濕潤體。對濕潤體運移規(guī)律與土質(zhì)關(guān)系進行了分析，圖2給出了滴頭流量為1 L/h，土壤容重為1.5 g/cm3，灌水時間為1 500 min的分析結(jié)果。

圖2 滴灌徑向和垂向遷移距離隨土壤顆粒質(zhì)量分?jǐn)?shù) 的分布變化 Fig.2 Variations of soil wetted volume with percentage of soil particles in radial and vertical directions

從圖2a中可以看出，濕潤體徑向遷移距離Ht隨土壤黏粒質(zhì)量分?jǐn)?shù)Cclay的增加而減小，在砂粒質(zhì)量分?jǐn)?shù)保持不變的條件下，徑向遷移距離與黏粒質(zhì)量分?jǐn)?shù)基本呈線性關(guān)系，即Ht∝Cclay；在相同的黏粒質(zhì)量分?jǐn)?shù)下，土壤砂粒質(zhì)量分?jǐn)?shù)的變化對徑向遷移距離影響較小，表明土壤黏粒質(zhì)量分?jǐn)?shù)對濕潤體徑向遷移距離具有重要影響。在圖2b中，濕潤體垂向遷移距離Dt隨土壤砂粒質(zhì)量分?jǐn)?shù)的增加而增大，但隨黏粒質(zhì)量分?jǐn)?shù)的增加而減小，因此垂向濕潤鋒在砂土中遷移較快，在粘土中發(fā)展較慢，且垂向遷移距離與土壤砂粒和黏粒質(zhì)量分?jǐn)?shù)分別滿足指數(shù)和冪函數(shù)關(guān)系。

土壤容重為1.5 g/cm3，灌水時間為1 500 min時，土壤滴灌濕潤體徑向和垂向遷移距離隨飽和導(dǎo)水率的變化規(guī)律如圖3所示。從圖中可以看出，不同滴灌流量下的徑向遷移距離總體表現(xiàn)為隨土壤飽和導(dǎo)水率的增大而增大，在土壤飽和導(dǎo)水率達到10 cm/h后呈現(xiàn)平穩(wěn)趨勢，其中粉質(zhì)土壤中徑向遷移距離在對應(yīng)的飽和導(dǎo)水率下表現(xiàn)較大增幅；而垂向遷移距離一直保持增加趨勢，且與飽和導(dǎo)水率具有較好的相關(guān)關(guān)系。這表明土壤飽和導(dǎo)水率增加到一定值后，濕潤體的徑向遷移受到限制，而垂向遷移則加速發(fā)展。

圖3 滴灌濕潤體徑向和垂向遷移距離隨土壤飽和 導(dǎo)水率的變化規(guī)律 Fig.3 Variations of soil wetted volume with soil saturated hydraulic conductivities in radial and vertical directions

2.2 濕潤體運移經(jīng)驗公式推導(dǎo)和檢驗

2.2.1經(jīng)驗公式推導(dǎo)

試驗表明，地表滴灌中土壤濕潤體形狀接近于半個橢球體[22]。通過模擬發(fā)現(xiàn)，在偏粘土中形成的濕潤體趨于平臥的半橢球體(徑向軸大于垂直軸)，而在砂土中易形成直立的半橢球體(徑向軸小于垂直軸)，因此不同土質(zhì)中濕潤體體積Vt可以用橢球體體積公式計算，即

(10)

大量的試驗及模擬研究表明，在土質(zhì)和滴頭流量確定的條件下，土壤中滴灌濕潤體體積Vt又與灌水量Vw，t之間存在顯著的正相關(guān)關(guān)系[7，10]

(11)

式中Δθ為滴灌過程中土壤濕潤體平均含水率的增量(cm3/cm3)，是與滴灌流量q(L/h)、土壤質(zhì)地和土壤初始含水率θ0有關(guān)的參數(shù)，對于特定土壤，在滴頭流量和土壤初始含水率確定的條件下，Δθ為一定值。通過模擬分析發(fā)現(xiàn)，在不同滴頭流量下，Δθ與滴頭流量成冪函數(shù)關(guān)系：Δθ∝q0.21，這與胡和平等[9]的研究結(jié)果一致。以土樣1(粉質(zhì)粘壤土)為例(圖4)，在不同滴頭流量下的濕潤體體積與灌水量滿足關(guān)系：Vt=4.88Vw，t/q0.21(R2=0.995)，因此，對于土樣1，不同滴頭流量下Δθ=Vw，t/Vt=q0.21/4.88=0.205q0.21。

圖4 土樣1(粉質(zhì)粘壤土)濕潤體體積與灌水量的關(guān)系 Fig.4 Relationship between wetted soil volume and irrigation amount for silty clay loam (No.1)

根據(jù)以上分析，在土壤質(zhì)地和初始含水率確定的條件下，濕潤體平均含水率的增量Δθ是關(guān)于滴頭流量q的函數(shù)。但對于不同質(zhì)地的土壤，Δθ不僅與滴頭流量q有關(guān)，而且與土壤性質(zhì)有關(guān)，若以飽和導(dǎo)水率Ks來代表土壤性質(zhì)，則Δθ是關(guān)于滴頭流量q和飽和導(dǎo)水率Ks的函數(shù)。通過對33個樣本點的Δθ與Ks進行分析(為消除q的影響，首先將不同流量下的Δθ除以q0.21，再與Ks進行回歸分析)，分析發(fā)現(xiàn)Δθ與Ks具有較好的負(fù)相關(guān)關(guān)系(圖5)。通過模擬結(jié)果發(fā)現(xiàn)，在僅用飽和導(dǎo)水率代表土壤性質(zhì)的情況下，不同土質(zhì)中Δθ與Ks滿足兩個冪函數(shù)關(guān)系

(12)

(13)

在應(yīng)用式(12)、(13)計算濕潤體平均含水率增量Δθ時，選擇公式的定量指標(biāo)是土壤的砂粒含量，土壤砂粒含量小于50%時，用式(12)求解，反之，用式(13)求解；定性指標(biāo)為是否為砂質(zhì)類土壤，若土壤偏砂質(zhì)(包括砂土、壤砂土、砂壤土、砂粘土和砂質(zhì)粘壤土)，用式(13)求解，否則用式(12)求解。

圖5 濕潤體平均含水率增量與飽和導(dǎo)水率的關(guān)系 Fig.5 Relationship between average water content increment in wetted volume and soil saturated hydraulic conductivity

若濕潤體平均含水率的增量Δθ已知時，根據(jù)式(10)和式(11)，濕潤體徑向遷移距離可表示為

(14)

由以上分析知，不同質(zhì)地土壤中的垂向遷移距離Dt與土壤飽和導(dǎo)水率Ks具有較好的相關(guān)關(guān)系(圖3b)。通過模擬分析，不同質(zhì)地土壤中Dt與Ks滿足經(jīng)驗關(guān)系(誤差10%的保證率為75%)

(15)

為提高經(jīng)驗公式的預(yù)測精度，根據(jù)土壤是否為偏砂質(zhì)，滴灌濕潤體垂向遷移距離Dt又可進一步表示為

(16)

(17)

圖6為滴頭流量為1～3 L/h，土壤容重為1.5 g/cm3，灌水時間為120～1 500 min時濕潤體徑向和垂向遷移距離的HYDRUS模擬值與經(jīng)驗公式計算值(式(14)、(15))的對比圖，從圖中可以看出，經(jīng)驗公式計算值與HYDRUS模擬值基本分布在1∶1線附近，且R2分別為0.937和0.959，表明經(jīng)驗公式具有一定的可靠性。

圖6 濕潤體徑向和垂向遷移距離的經(jīng)驗公式計算值 與HYDRUS模擬值對比 Fig.6 Comparison of soil wetted volume development between results of empirical method and HYDRUS simulation

2.2.2經(jīng)驗公式檢驗

在給定土壤飽和導(dǎo)水率Ks和滴頭流量q的前提下(q=1～3 L/h，初始土壤體積含水率為0.1 cm3/cm3左右)，應(yīng)用式(15)或式(16)、(17)可以估算出滴灌濕潤體垂向遷移距離Dt，將計算出來的Dt代入式(14)，即可計算出徑向遷移距離Ht。

圖7 砂壤土中徑向和垂向遷移距離隨時間變化情況的 經(jīng)驗公式計算結(jié)果與通用算法解(文獻[12])對比 (飽和導(dǎo)水率Ks=80.64 cm/d) Fig.7 Comparison of soil wetted volume development between results of empirical method and empirical algorithm solution for sandy loam (Ref.12, Ks=80.64 cm/d)

分別利用文獻[6、8、12、19、23]中的試驗數(shù)據(jù)或經(jīng)驗?zāi)Ｐ?，對本研究中建立的滴灌濕潤體運移經(jīng)驗公式(式(14)～(17))進行檢驗。圖7中濕潤體徑向和垂向遷移距離的通用算法解是文獻[12]中根據(jù)土壤水力特性參數(shù)(θr、θs、α、n和Ks)通過無量綱計算和有量綱轉(zhuǎn)化得到的結(jié)果，利用本研究中的經(jīng)驗公式僅根據(jù)飽和導(dǎo)水率Ks也能得到可靠的計算值。圖8～11分別給出了濕潤體運移的試驗實測值與經(jīng)驗公式計算值的對比，從圖中可以看出在砂土、砂壤土條件下經(jīng)驗公式較好地預(yù)測了徑向和垂向的遷移過程，而對粘壤土中徑向濕潤鋒遷移過程的預(yù)測效果稍差，主要是由于土壤質(zhì)地偏粘，滴灌過程中容易形成地表積水。另外，在經(jīng)驗公式與文獻中經(jīng)驗?zāi)Ｐ秃驮囼灁?shù)據(jù)的統(tǒng)計分析中，砂土、砂壤土和壤土條件下經(jīng)驗公式計算的濕潤鋒遷移距離的均方根誤差RMSE均小于3 cm，決定系數(shù)R2都在0.75以上，粘壤土條件下由于地表積水，經(jīng)驗公式計算的均方根誤差稍大(4.83 cm)。通過將本研究中經(jīng)驗公式與文獻中經(jīng)驗?zāi)Ｐ秃驮囼灁?shù)據(jù)的對比與統(tǒng)計分析，表明不同土質(zhì)中滴灌濕潤體經(jīng)驗公式具有較高的預(yù)測精度。

圖8 砂土中徑向和垂向遷移距離隨時間變化情況的 經(jīng)驗公式計算結(jié)果與實測(文獻[8])對比 (飽和導(dǎo)水率Ks=57.6 cm/d) Fig.8 Comparison of soil wetted volume development between results of empirical method and laboratory measurements for sandy soil (Ref. 8, Ks=57.6 cm/d)

圖9 砂壤土中徑向和垂向遷移距離隨時間變化情況的 經(jīng)驗公式計算結(jié)果與實測(文獻[19])對比 (飽和導(dǎo)水率Ks=47.04 cm/d) Fig.9 Comparison of soil wetted volume development between results of empirical method and laboratory measurements for sandy loam (Ref. 19, Ks=47.04 cm/d)

圖10 壤土中徑向遷移距離隨時間變化情況的 經(jīng)驗公式計算結(jié)果與實測(文獻[23])對比 (飽和導(dǎo)水率Ks=47.52 cm/d) Fig.10 Comparison of soil wetted volume development between results of empirical method and laboratory measurements for loam (Ref. 23， Ks=47.52 cm/d)

圖11 粘壤土中徑向和垂向遷移距離隨時間變化 情況的經(jīng)驗公式計算結(jié)果與實測(文獻[6])對比 (飽和導(dǎo)水率Ks=11.59 cm/d) Fig.11 Comparison of soil wetted volume development between results of empirical method and laboratory measurements for clay loam (Ref. 6, Ks=11.59 cm/d)

2.3 經(jīng)驗公式適用性

土壤初始含水率對滴灌濕潤體的運移具有一定的影響，為分析本研究中濕潤體運移經(jīng)驗公式對于不同初始含水率的適用性，將經(jīng)驗公式在不同土質(zhì)和滴灌流量下的計算值與在相應(yīng)條件下的HYDRUS模擬值進行比較，其中HYDRUS模型分別采用凋萎點含水率(基質(zhì)勢為-1 500 kPa)[24]和0.15 cm3/cm3兩種初始含水率。從圖12可以看出，垂向濕潤鋒的經(jīng)驗公式計算值與不同初始含水率下的HYDRUS模擬值比較接近，平均相對誤差MRE均在10%左右；結(jié)果表明徑向濕潤鋒計算值的平均相對誤差MRE也均在10%以內(nèi)，表明經(jīng)驗公式適用于土壤含水率比較低的情況下滴灌濕潤體遷移過程的預(yù)測。

本研究中計算滴灌濕潤體遷移的經(jīng)驗公式是關(guān)于滴灌流量、滴灌時間和飽和導(dǎo)水率的函數(shù)，具有計算過程簡便、預(yù)測精度較高的特點，能為滴灌系統(tǒng)設(shè)計提供參考依據(jù)。但應(yīng)用本經(jīng)驗公式進行滴灌設(shè)計又受到一定的條件限制：①土壤的初始含水率不宜太高。②在預(yù)測滴灌濕潤體在不同土質(zhì)中的遷移狀況時，對土壤飽和導(dǎo)水率的可靠性要求較高。

圖12 垂向濕潤鋒的經(jīng)驗公式計算值與不同初始 含水率下HYDRUS模擬值對比 Fig.12 Comparison of soil wetted volume development in vertical direction between results of empirical method and HYDRUS simulation under different soil initial water contents

3 結(jié)束語

基于非飽和土壤水分運動基本原理，首先，采用HYDRUS-2D/3D軟件模擬了多種土質(zhì)地表滴灌條件下的土壤水分運動過程，根據(jù)數(shù)值模擬結(jié)果分析了土壤的顆粒含量和飽和導(dǎo)水率對濕潤體運移的影響規(guī)律。進而以飽和導(dǎo)水率作為土壤質(zhì)地的代表參數(shù)，建立了描述不同土質(zhì)中地表滴灌濕潤體運移的經(jīng)驗公式。經(jīng)驗公式中滴灌濕潤體的動態(tài)變化僅是滴灌流量、滴灌時間和土壤飽和導(dǎo)水率的函數(shù)，簡單、實用。與數(shù)值模擬結(jié)果、文獻試驗結(jié)果的對比表明，該經(jīng)驗公式對不同土質(zhì)下的濕潤體運移規(guī)律預(yù)測效果較好，且在初始含水率偏小、土壤飽和導(dǎo)水率可靠的情況下，模擬精度更高。此經(jīng)驗公式可為農(nóng)業(yè)生產(chǎn)中地表滴灌的設(shè)計和應(yīng)用提供指導(dǎo)。