蘇李君 藺樹棟 王全九,3 王 康

(1.西安理工大學(xué)省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室， 西安 710048；2.西安理工大學(xué)理學(xué)院， 西安 710054；3.中國(guó)科學(xué)院水利部水土保持研究所黃土高原土壤侵蝕與旱地農(nóng)業(yè)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 陜西楊凌 712100)

0 引言

我國(guó)人均水資源量偏低，隨著人口的增加和經(jīng)濟(jì)的發(fā)展，愈來(lái)愈多的地區(qū)出現(xiàn)水資源供需緊張的狀況，北方地區(qū)的水資源已經(jīng)不能支撐其經(jīng)濟(jì)持續(xù)健康的發(fā)展，困境日趨明顯，降水更是稀少不均。面對(duì)現(xiàn)有水供應(yīng)日趨增加的態(tài)勢(shì)，特別是農(nóng)業(yè)灌水對(duì)精心設(shè)計(jì)和管理的灌溉系統(tǒng)的需求，合理設(shè)計(jì)灌溉系統(tǒng)、提高灌溉用水效率顯得尤為重要[1]。

滴灌技術(shù)以其顯著的節(jié)水、節(jié)肥和增產(chǎn)等優(yōu)點(diǎn)在我國(guó)已得到大面積的推廣和應(yīng)用，是目前最先進(jìn)的節(jié)水灌溉技術(shù)之一。滴灌濕潤(rùn)體和主要根系分布區(qū)域的一致性是提高灌溉水利用效率的基礎(chǔ)[2-3]。分析灌溉過(guò)程中不同時(shí)期的入滲趨勢(shì)、累積入滲率以及濕潤(rùn)鋒的遷移距離，可促進(jìn)灌溉系統(tǒng)的合理設(shè)計(jì)[1]。滴灌地表積水區(qū)(或稱飽和區(qū))對(duì)土壤濕潤(rùn)區(qū)的影響很早就被研究者所關(guān)注[4-7]，它對(duì)土壤濕潤(rùn)區(qū)有顯著的推動(dòng)作用。在滴灌技術(shù)設(shè)計(jì)中，土壤濕潤(rùn)區(qū)是確定土壤水平濕潤(rùn)鋒運(yùn)移距離與垂直濕潤(rùn)鋒運(yùn)移距離比值(簡(jiǎn)稱“濕潤(rùn)鋒比”)、灌水定額、以及滴頭間距等參數(shù)的基礎(chǔ)[8]，不少學(xué)者一直在探究滴灌條件下土壤水分運(yùn)動(dòng)的運(yùn)移機(jī)理。一般認(rèn)為，當(dāng)供水強(qiáng)度小于土壤入滲能力時(shí)，地表不形成積水；但當(dāng)供水強(qiáng)度超過(guò)土壤入滲能力后，地表形成積水，并同時(shí)產(chǎn)生徑流[9]。此外，F(xiàn)U等[1]還根據(jù)實(shí)驗(yàn)數(shù)據(jù)提出了累積入滲量和濕潤(rùn)鋒運(yùn)移距離的模型，研究結(jié)果可為灌溉系統(tǒng)的設(shè)計(jì)提供理論支持，提高灌溉用水效率。胡和平等[10]在地表滴灌條件下，運(yùn)用SWMS-2D 模型對(duì)滴頭流量、土壤初始含水率和飽和導(dǎo)水率等多種濕潤(rùn)體運(yùn)移情況進(jìn)行模擬分析，同時(shí)建立了土壤濕潤(rùn)體運(yùn)移的經(jīng)驗(yàn)方程。以上研究表明，滴灌水分入滲規(guī)律在灌溉系統(tǒng)中起到非常重要的作用，但是目前關(guān)于滴灌積水入滲模型仍需要進(jìn)行深入研究，尤其是能夠描述點(diǎn)源積水入滲機(jī)理的模型對(duì)研究滴灌土壤濕潤(rùn)區(qū)變化規(guī)律更為重要。因此，研究點(diǎn)源積水二維入滲條件下濕潤(rùn)鋒運(yùn)移規(guī)律對(duì)充分發(fā)揮滴灌的技術(shù)優(yōu)勢(shì)具有重要意義。

在前人研究的基礎(chǔ)上，本文以達(dá)西定律和質(zhì)量守恒原理為基礎(chǔ)，假設(shè)濕潤(rùn)體形狀和濕潤(rùn)鋒比不隨時(shí)間變化，僅濕潤(rùn)體大小隨時(shí)間的增加而發(fā)生變化，在此基礎(chǔ)上建立一個(gè)基于土壤水力參數(shù)的二維濕潤(rùn)鋒運(yùn)移模型。采用HYDRUS-3D軟件模擬二維積水入滲的水分運(yùn)動(dòng)過(guò)程。通過(guò)軟件模擬與數(shù)據(jù)分析，在研究濕潤(rùn)鋒運(yùn)移隨時(shí)間變化規(guī)律的基礎(chǔ)上，建立點(diǎn)源積水入滲情況下濕潤(rùn)鋒比以及水平濕潤(rùn)鋒與土壤水力參數(shù)間的模型，并對(duì)模型進(jìn)行驗(yàn)證。

1 模型構(gòu)建和研究方法

1.1 數(shù)學(xué)模型

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

在膜下滴灌條件下，土壤水分運(yùn)動(dòng)屬于三維土壤水分入滲問(wèn)題，但可簡(jiǎn)化為二維問(wèn)題來(lái)考慮[11]。根據(jù)達(dá)西定律、質(zhì)量守恒定律，在忽略溫度、空氣阻力以及土壤水滯后效應(yīng)等因素對(duì)土壤水分入滲影響的同時(shí)，并假設(shè)土壤為各向同性，入滲模型可簡(jiǎn)化為以積水中心為對(duì)稱點(diǎn)的二維入滲，此時(shí)土壤水分運(yùn)動(dòng)可用二維Richards方程表示為

(1)

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

t——時(shí)間，min

x、z——空間坐標(biāo)值，cm

D(θ)——非飽和擴(kuò)散率，cm2/d

K(θ)——非飽和導(dǎo)水率，cm/min

1.1.2土壤水力模型

HYDRUS-3D是一款用于模擬二維、三維土壤水分運(yùn)動(dòng)、溶質(zhì)運(yùn)移、根系吸水及熱量傳輸?shù)挠邢拊?jì)算軟件。本文選取HYDRUS-3D中二維土壤水分運(yùn)動(dòng)模擬模塊，該模塊在進(jìn)行數(shù)值模擬時(shí)，需要土壤水分特征曲線的參數(shù)，主要包括土壤滯留含水率、土壤飽和含水率、土壤形狀參數(shù)、土壤進(jìn)氣吸力、飽和導(dǎo)水率以及土壤孔隙連通性系數(shù)。土壤水力特征模型采用Brooks-Corey模型，即

(2)

(3)

其中

m=3n+2

式中Ks——飽和導(dǎo)水率，cm/min

n——形狀系數(shù)m——參數(shù)

S——有效飽和度

θr——滯留含水率，cm3/cm3

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

h——土壤吸力，cm

hd——進(jìn)氣吸力，cm

選取表1中的土壤模擬積水半徑為1 cm的二維積水入滲過(guò)程，根據(jù)模擬得到的數(shù)據(jù)分析并構(gòu)建相關(guān)數(shù)學(xué)模型，土壤水力參數(shù)選取HYDRUS-3D軟件中自帶參數(shù)[12]，見(jiàn)表1，表中θ0為土壤初始含水率。

表1 土壤水力參數(shù)Tab.1 Soil hydraulic parameters

為了驗(yàn)證本文構(gòu)建模型的準(zhǔn)確性，另采用不同的土壤水力參數(shù)對(duì)所建模型進(jìn)行驗(yàn)證，選取土壤參數(shù)方法是在表1中相應(yīng)土壤水力參數(shù)的基礎(chǔ)上通過(guò)合理調(diào)整n和hd獲得，見(jiàn)表2。

表2 驗(yàn)證采用的土壤水力參數(shù)Tab.2 Soil hydraulic parameters verified

注：*表示調(diào)整參數(shù)后的土壤類型。

1.1.3模型設(shè)定與時(shí)空離散

HYDRUS-3D工作空間的幾何模型設(shè)為X-Z平面，所建模型區(qū)域?yàn)?00 cm×100 cm的正方形平面，入滲源為半徑為1 cm的積水區(qū)域，位于上邊界中心點(diǎn)位置處。如圖1a建立直角坐標(biāo)系。

建立幾何模型后將Brooks-Corey模型作為土壤水力特征模型，不考慮水力滯后效應(yīng)。模擬初始時(shí)間為0 min，結(jié)束時(shí)間為2 000 min。設(shè)置最大迭代次數(shù)(Maximum number of iterations)為10 000次，設(shè)置初始時(shí)間步長(zhǎng)(Initial time)0.5 min，最小時(shí)間步長(zhǎng)(Minimum time)0.1 min，最大時(shí)間步長(zhǎng)(Maximum time)1 min，含水率允許最大誤差限(Water content tolerance)0.001 cm3/cm3，壓力水頭允許最大誤差限(Pressure head tolerance)1 cm。采用HYDRUS-3D中FE-Mesh模塊對(duì)所建的二維模型進(jìn)行網(wǎng)格劃分，網(wǎng)格劃分的類型為三角形，有限元柵格劃分參數(shù)中設(shè)置其有限元三角形外接圓最大直徑為1 cm，其他參數(shù)采用HYDRUS數(shù)據(jù)庫(kù)中自帶數(shù)據(jù)。為了提高模擬結(jié)果的精度，對(duì)積水周圍的部分區(qū)域網(wǎng)格進(jìn)行加密處理，有限元計(jì)算網(wǎng)格如圖1b所示。

圖1 計(jì)算區(qū)域及網(wǎng)格劃分Fig.1 Simulation areas and meshing

計(jì)算區(qū)域入滲面EF是半徑為1 cm的積水面，在不入滲時(shí)視地表為大氣邊界，入滲開始后很快達(dá)到飽和，視為定水頭邊界；入滲面EF左右兩側(cè)邊界，以及邊界AD和BC假定為零通量邊界，故其水分通量為零；下邊界CD為不受地下水位影響的自由出流邊界。在入滲開始時(shí)，即在t=0時(shí)刻，土壤含水率為土壤初始含水率θ0。綜上所述，其初始條件和邊界條件可總結(jié)為：

(1)初始條件

初始條件為

θ(x,z,0)=θ0
(-50 cm≤x≤50 cm，0 cm≤z≤100 cm)

(4)

(2)邊界條件

入滲面EF邊界

θ(x,z,t)=θs(x0≤x≤xl，t>0)

(5)

CD自由排水邊界

(6)

AE、FB、AD、BC零通量邊界

(7)

式中x0——積水區(qū)域左邊界，cm

xl——積水區(qū)域右邊界，cm

1.1.4土壤濕潤(rùn)體形狀模型

濕潤(rùn)體的大小與形狀對(duì)作物生長(zhǎng)起關(guān)鍵性作用，土壤水分入滲濕潤(rùn)體的研究是滴灌以及涌泉根灌技術(shù)要素和水分運(yùn)動(dòng)規(guī)律研究的基礎(chǔ)[13]。研究表明積水入滲的濕潤(rùn)體形狀可近似看成一個(gè)橢球或橢圓面[14-16]。針對(duì)二維土壤積水入滲問(wèn)題，本文假設(shè)土壤濕潤(rùn)體是以垂直濕潤(rùn)鋒和水平濕潤(rùn)鋒為長(zhǎng)半軸和短半軸的半橢圓面

(8)

式中xf——水平濕潤(rùn)鋒長(zhǎng)度，cm

zf——垂直濕潤(rùn)鋒長(zhǎng)度，cm

1.1.5濕潤(rùn)鋒運(yùn)移模型

土壤濕潤(rùn)鋒分布一方面可以反映灌溉水量是否滿足灌溉要求，另一方面可以反映灌溉質(zhì)量，其運(yùn)動(dòng)規(guī)律是灌溉管理的理論基礎(chǔ)。不同灌溉技術(shù)因濕潤(rùn)土壤的方式不同，濕潤(rùn)鋒運(yùn)移規(guī)律存在明顯差異，而同一灌溉技術(shù)條件下，又因土壤理化性質(zhì)、灌溉參數(shù)以及土壤水力參數(shù)的不同，濕潤(rùn)鋒運(yùn)移也不一致[3]。研究表明積水入滲的濕潤(rùn)鋒運(yùn)移過(guò)程可采用入滲時(shí)間的冪函數(shù)表示[14,17-18]，即

xf=a1tb1

(9)

zf=a2tb2

(10)

式中a1、a2、b1、b2——擬合參數(shù)

土壤水分?jǐn)U散率反映了土壤孔隙度、孔隙大小分布以及導(dǎo)水性能，并影響土壤中水分運(yùn)動(dòng)狀況，它是研究水鹽運(yùn)動(dòng)規(guī)律與調(diào)控以及表征土壤水動(dòng)力學(xué)的重要參數(shù)之一[19]。因此本文考慮構(gòu)建飽和擴(kuò)散率與冪函數(shù)參數(shù)之間的關(guān)系。土壤非飽和擴(kuò)散率的表達(dá)式[20]為

(11)

令θ=θs，得到土壤飽和擴(kuò)散率的表達(dá)式為

(12)

1.2 數(shù)據(jù)處理及誤差分析

數(shù)據(jù)運(yùn)用Excel進(jìn)行處理，并用AutoCAD和OriginPro進(jìn)行繪圖以及Matlab進(jìn)行模型驗(yàn)證，同時(shí)采用決定系數(shù)R2、均方根誤差(RMSE)和相對(duì)誤差(RE)評(píng)價(jià)指標(biāo)進(jìn)行誤差分析，檢驗(yàn)相關(guān)性。

2 結(jié)果與分析

2.1 土壤濕潤(rùn)體形狀變化特征

圖2 土壤入滲濕潤(rùn)體形狀Fig.2 Soil infiltration wetting pattern shapes

以積水中心為坐標(biāo)原點(diǎn)，用HYDRUS-3D模擬土壤水分在水平方向和垂直方向上的二維入滲，不同時(shí)刻濕潤(rùn)鋒輪廓如圖2所示。從圖中可以看出，在整個(gè)入滲過(guò)程中水平擴(kuò)散距離和垂直入滲深度均隨著入滲時(shí)間的增加而增加，但濕潤(rùn)體運(yùn)移速率在后期明顯逐漸減小。入滲開始后土壤表層很快達(dá)到飽和，視為定水頭，隨著入滲時(shí)間的增加，同一土壤的濕潤(rùn)鋒運(yùn)移垂直入滲深度略快于水平擴(kuò)散距離。這是由于隨著時(shí)間的增長(zhǎng)，土水勢(shì)梯度減小，在基質(zhì)勢(shì)和重力勢(shì)的共同作用下，導(dǎo)致土壤水向下運(yùn)動(dòng)能力增強(qiáng)，使?jié)駶?rùn)體的形狀近似可以看作一個(gè)橢圓面[21]，細(xì)質(zhì)土壤(hd>11 cm)的濕潤(rùn)體形狀變化過(guò)程可以假設(shè)為以積水中心處為圓心的橢圓形，但粗質(zhì)土壤(壤砂土和砂土，hd<8.7 cm)的濕潤(rùn)體形狀變化過(guò)程不是以積水中心為圓心的橢圓形，隨著入滲時(shí)間的增加橢圓中心逐漸向積水中心下方偏移(如圖2k、2l中壤砂土和砂土濕潤(rùn)體形狀)，因此本文研究所得模型適用于細(xì)質(zhì)土壤。

2.2 水平濕潤(rùn)鋒與垂直濕潤(rùn)鋒的關(guān)系

通過(guò)數(shù)值模擬可知，土壤在不同入滲時(shí)間的濕潤(rùn)鋒形狀基本相似，但水平和垂直濕潤(rùn)鋒長(zhǎng)度隨時(shí)間發(fā)生變化。選取砂土、壤砂土、粘土和砂質(zhì)粘壤土為例，表3給出了這4種土壤濕潤(rùn)鋒長(zhǎng)度及濕潤(rùn)鋒比隨時(shí)間的變化過(guò)程。由表可看出，細(xì)質(zhì)土壤(hd>11 cm)的濕潤(rùn)鋒比隨時(shí)間變化幅度較小，垂直濕潤(rùn)鋒與水平濕潤(rùn)鋒變化速率基本保持一致，但粗質(zhì)土壤(砂土和壤砂土，hd<8.7 cm)的濕潤(rùn)鋒比隨時(shí)間變化幅度較大，隨著入滲時(shí)間的增加垂直濕潤(rùn)鋒變化速率比水平濕潤(rùn)鋒變化速率較快，因此濕潤(rùn)鋒比隨時(shí)間的增加而減小，同時(shí)可以結(jié)合圖2發(fā)現(xiàn)，粗質(zhì)土壤濕潤(rùn)體橢圓中心逐漸向積水中心下方偏移，所以本研究主要考慮細(xì)質(zhì)土壤，并假設(shè)細(xì)質(zhì)土壤水平濕潤(rùn)鋒和垂直濕潤(rùn)鋒的比值不隨入滲時(shí)間變化，令

圖3 濕潤(rùn)鋒比與土壤水力參數(shù)關(guān)系Fig.3 Relationship between wetting frontratios and soil hydraulic parameters

表3 濕潤(rùn)鋒比隨時(shí)間的變化
Tab.3 Wetting frontratio varied with time

土壤類型時(shí)間/min水平濕潤(rùn)鋒長(zhǎng)度/cm垂直濕潤(rùn)鋒長(zhǎng)度/cm濕潤(rùn)鋒比5022.46029.5300.76110025.86034.7310.745砂土20030.19552.5230.57550040.36189.0110.453100049.001134.5000.364200058.474215.0100.272509.65012.3850.77910013.20517.5880.751壤砂土20017.11523.1980.73850025.10236.5440.687100033.04854.5230.606200042.55080.5560.528504.6664.9670.9391005.5755.9220.941粘土2006.8227.2550.9405009.2299.8170.940100011.87412.6100.942200015.53616.5120.941507.0127.7620.9031008.95010.0500.891砂質(zhì)粘壤土20011.60813.1100.88550016.73518.8110.890100022.36025.1820.888200030.13733.8610.890

(13)

式中wt——水平濕潤(rùn)鋒和垂直濕潤(rùn)鋒長(zhǎng)度的比值，即濕潤(rùn)鋒比

濕潤(rùn)鋒比可以作為一個(gè)確定滴灌灌水參數(shù)的指標(biāo)[22]，不同土壤的濕潤(rùn)鋒比差異較大，分別分析濕潤(rùn)鋒比與土壤水力參數(shù)(n、hd、Ks)之間的關(guān)系，如圖3所示。由圖可知，相對(duì)于n和Ks，wt與hd具有較好的線性關(guān)系，通過(guò)擬合得到wt與hd的關(guān)系為

wt=0.004hd+0.787 (R2=0.999)

(14)

結(jié)合式(13)、(14)，可得垂直濕潤(rùn)鋒與水平濕潤(rùn)鋒之間的關(guān)系為

(15)

由式(15)可知，在二維積水入滲條件下已知土壤某一入滲時(shí)刻水平濕潤(rùn)鋒長(zhǎng)度和土壤進(jìn)氣吸力時(shí)，就可以通過(guò)濕潤(rùn)鋒比計(jì)算得到此時(shí)刻垂直方向上的濕潤(rùn)鋒運(yùn)移距離。

2.3 水平濕潤(rùn)鋒運(yùn)移特征

由式(15)可知，若水平濕潤(rùn)鋒長(zhǎng)度已知，則可以由土壤進(jìn)氣吸力估計(jì)二維積水入滲的土壤濕潤(rùn)體范圍。雖然水平濕潤(rùn)鋒長(zhǎng)度可以直接在室內(nèi)或田間試驗(yàn)中觀測(cè)得到，但如果水平濕潤(rùn)鋒長(zhǎng)度與入滲時(shí)間之間的關(guān)系未知，則本文方法就不能直接預(yù)測(cè)濕潤(rùn)體范圍。因此有必要進(jìn)一步分析土壤水平濕潤(rùn)鋒長(zhǎng)度與時(shí)間的變化關(guān)系。不同土壤濕潤(rùn)體水平濕潤(rùn)鋒長(zhǎng)度隨時(shí)間變化如圖4所示?？梢钥闯?，積水入滲條件下濕潤(rùn)鋒長(zhǎng)度與時(shí)間存在顯著的冪函數(shù)關(guān)系，這與其他學(xué)者的一些研究結(jié)果相同[14,17-18]。對(duì)表1中細(xì)質(zhì)土壤的水平濕潤(rùn)鋒長(zhǎng)度運(yùn)用式(9)進(jìn)行擬合，結(jié)果如表4所示，決定系數(shù)R2均大于0.990，擬合效果較好。

圖4 水平濕潤(rùn)鋒長(zhǎng)度隨時(shí)間的變化曲線Fig.4 Variation curves of horizontal wetting front with time

表4 水平濕潤(rùn)鋒長(zhǎng)度擬合結(jié)果
Tab.4 Horizontal wetting front fitting results

土壤類型擬合結(jié)果R2粘土xf=1.260t0.3300.996壤土xf=1.442t0.3880.996粉砂土xf=1.416t0.3810.999粉壤土xf=1.389t0.3770.998粘壤土xf=1.329t0.3580.997砂壤土xf=1.679t0.4240.999砂質(zhì)粘土xf=1.330t0.3570.998粉質(zhì)粘土xf=1.205t0.3270.996砂質(zhì)粘壤土xf=1.490t0.3930.999粉質(zhì)粘壤土xf=1.316t0.3440.997

對(duì)式(9)求導(dǎo)可得

(16)

計(jì)算得到水平濕潤(rùn)鋒長(zhǎng)度變化率隨時(shí)間變化的關(guān)系曲線，如圖5所示。

圖5 水平濕潤(rùn)鋒長(zhǎng)度變化率隨時(shí)間的變化曲線Fig.5 Changing curves of variation rate of horizontal wetting front with time

式(9)表明水平濕潤(rùn)鋒長(zhǎng)度隨著入滲時(shí)間的增加而呈冪函數(shù)增大，由圖5可以看出，在入滲初始階段約200 min內(nèi)水平濕潤(rùn)鋒長(zhǎng)度迅速增大，約200 min后，濕潤(rùn)鋒長(zhǎng)度緩慢增大。在入滲過(guò)程中水平濕潤(rùn)鋒長(zhǎng)度變化速率隨入滲時(shí)間的增加而持續(xù)減小，在入滲初始的500 min內(nèi)，下降迅速，約1 000 min后下降緩慢，基本趨于穩(wěn)定。

2.4 濕潤(rùn)鋒運(yùn)移定量分析

圖6 水平濕潤(rùn)鋒運(yùn)移模型參數(shù)與飽和擴(kuò)散率的關(guān)系Fig.6 Relationship between parameters of horizontal wetting front migration model and saturation diffusivity

土壤水力參數(shù)是求解土壤水分運(yùn)動(dòng)方程的重要參數(shù)，也是建立土壤水分運(yùn)動(dòng)數(shù)學(xué)模型的重要基礎(chǔ)。在研究水平濕潤(rùn)鋒運(yùn)移隨時(shí)間變化規(guī)律的基礎(chǔ)上，運(yùn)用HYDRUS軟件模擬的數(shù)據(jù)進(jìn)行擬合研究發(fā)現(xiàn)，不同土壤入滲時(shí)的冪函數(shù)參數(shù)a1與土壤飽和擴(kuò)散率Ds間呈線性關(guān)系，其中土壤飽和擴(kuò)散率Ds可將表1中的參數(shù)代入式(12)計(jì)算得到，變化關(guān)系如圖6a所示，決定系數(shù)R2=0.984，RMSE為0.017，RE為0.013%；同時(shí)，擬合發(fā)現(xiàn)參數(shù)b1是土壤飽和擴(kuò)散率Ds的二次函數(shù)，如圖6b所示，決定系數(shù)R2=0.981，RMSE為0.004，RE為0.011%，數(shù)值模擬效果較好。

經(jīng)過(guò)擬合得到參數(shù)a1、b1與Ds間滿足

a1=0.098Ds+1.175

(17)

(18)

聯(lián)立式(9)、(17)、(18)可以推導(dǎo)得出水平濕潤(rùn)鋒長(zhǎng)度與土壤飽和擴(kuò)散率間的關(guān)系，即

(19)

綜合式(15)和式(19)可知，在給定土壤水力參數(shù)(包括n、Ks、hd、θs、θr)時(shí)，就可以估算得到某一時(shí)刻的水平濕潤(rùn)鋒長(zhǎng)度和垂直濕潤(rùn)鋒長(zhǎng)度。

2.5 濕潤(rùn)體變化過(guò)程定量關(guān)系評(píng)價(jià)

采用表2中的水力參數(shù)分別對(duì)濕潤(rùn)鋒比與進(jìn)氣吸力的關(guān)系、冪函數(shù)參數(shù)a1和b1與飽和擴(kuò)散率的關(guān)系、水平濕潤(rùn)鋒模型以及土壤濕潤(rùn)體形狀進(jìn)行驗(yàn)證。實(shí)測(cè)值和擬合值的吻合程度運(yùn)用決定系數(shù)R2、RMSE和RE評(píng)價(jià)指標(biāo)來(lái)進(jìn)行誤差分析，檢驗(yàn)?zāi)Ｐ偷木_性。

圖7是對(duì)濕潤(rùn)鋒比和冪函數(shù)參數(shù)a1、b1的驗(yàn)證，由圖7可以看出，濕潤(rùn)鋒比、冪函數(shù)參數(shù)a1、b1的實(shí)測(cè)值與擬合值之間有較好的吻合度，決定系數(shù)R2均大于0.950，RMSE均不大于0.027，RE不大于0.124%，表明式(14)可以很好地用來(lái)描述濕潤(rùn)鋒比與土壤進(jìn)氣吸力間的線性變化關(guān)系；式(17)、(18)可以較好地表示冪函數(shù)參數(shù)a1、b1與土壤飽和擴(kuò)散率間的變化關(guān)系。

圖7 模型驗(yàn)證Fig.7 Model validation

圖8 水平濕潤(rùn)鋒驗(yàn)證Fig.8 Horizontal wetting front validation

由圖4、5已經(jīng)得出，從1 000 min左右開始水平濕潤(rùn)鋒長(zhǎng)度變化速率基本趨于穩(wěn)定，隨著時(shí)間的推移，濕潤(rùn)鋒運(yùn)移變化規(guī)律無(wú)限趨于線性變化，也就意味著從此時(shí)開始計(jì)算值略小于實(shí)測(cè)值，而在實(shí)際研究中入滲時(shí)長(zhǎng)并不會(huì)很長(zhǎng)，因此只對(duì)水平濕潤(rùn)鋒運(yùn)移數(shù)學(xué)模型(式(19))前1 000 min進(jìn)行驗(yàn)證，如圖8所示。驗(yàn)證結(jié)果顯示，其決定系數(shù)R2在0.941～0.999之間，RMSE在0.219～1.150 cm之間，RE在0.015%～0.954%之間，在一定誤差允許范圍內(nèi)擬合效果較好。

圖9是對(duì)細(xì)質(zhì)土壤濕潤(rùn)體形狀的驗(yàn)證，圖中實(shí)測(cè)值表示水分入滲的實(shí)際濕潤(rùn)體輪廓線，擬合值表示通過(guò)式(8)計(jì)算得到的標(biāo)準(zhǔn)橢圓曲線。從圖中可以看出，在入滲時(shí)刻50、100、200、500 min的實(shí)際濕潤(rùn)體形狀和計(jì)算得到的標(biāo)準(zhǔn)橢圓曲線基本吻合，驗(yàn)證結(jié)果顯示決定系數(shù)R2在0.920～0.972之間，RMSE在0.096～2.542 cm之間，RE在0.330%～5.970%之間，表明在一定的誤差范圍內(nèi)二維積水入滲的土壤濕潤(rùn)體形狀可近似看成一個(gè)以垂直濕潤(rùn)鋒和水平濕潤(rùn)鋒分別為長(zhǎng)半軸和短半軸的半橢圓面，同時(shí)這一研究結(jié)果更加證明了其他學(xué)者[14-16]對(duì)土壤濕潤(rùn)體形狀研究結(jié)論的準(zhǔn)確性。

綜上所述，在誤差允許范圍內(nèi)，本研究所建立的模型可以為分析點(diǎn)源積水水分運(yùn)移、滴灌灌水器的選擇和布設(shè)以及為合理設(shè)計(jì)滴灌系統(tǒng)等相關(guān)問(wèn)題提供科學(xué)依據(jù)。

圖9 濕潤(rùn)體形狀驗(yàn)證Fig.9 Wetting pattern shape validation

3 討論

本文在Richards方程的基礎(chǔ)上采用Brooks-Corey土壤水力特征模型，并假設(shè)濕潤(rùn)體形狀和濕潤(rùn)鋒比不隨時(shí)間變化、積水中心為橢圓圓心的情況下推導(dǎo)出點(diǎn)源積水半徑為1 cm的二維入滲濕潤(rùn)鋒比及濕潤(rùn)鋒長(zhǎng)度與土壤水力參數(shù)間的數(shù)學(xué)模型。

通過(guò)數(shù)值模擬發(fā)現(xiàn)進(jìn)氣吸力小于8.7 cm的粗質(zhì)土壤(砂土和壤砂土)濕潤(rùn)鋒比隨時(shí)間的增加呈遞減趨勢(shì)，這是因?yàn)殡S著入滲時(shí)間的增加，重力勢(shì)作用明顯，粗質(zhì)土壤垂直濕潤(rùn)鋒運(yùn)移速率大于水平濕潤(rùn)鋒運(yùn)移速率。因此本研究只針對(duì)進(jìn)氣吸力大于11 cm的細(xì)質(zhì)土壤進(jìn)行了詳細(xì)的研究。同時(shí)從圖8可以看出，粉砂土*和粘壤土*的實(shí)測(cè)值和擬合值吻合效果相對(duì)較差，決定系數(shù)R2分別為0.941、0.961；砂質(zhì)粘壤土*和壤土*效果較好，R2分別為0.975、0.966；粘土*、砂質(zhì)粘土*、粉壤土*、粉質(zhì)黏土*、粉砂粘壤土*和砂壤土*實(shí)測(cè)值和擬合值的吻合效果很好，R2分別為0.994、0.982、0.981、0.992、0.996、0.999。總體看來(lái)10種土壤實(shí)測(cè)值和擬合值吻合效果相對(duì)較好，引起偏差的主要因素可能是受冪函數(shù)系數(shù)a1的影響較大，對(duì)參數(shù)a1驗(yàn)證發(fā)現(xiàn)，實(shí)測(cè)值與擬合值并不完全一致，部分驗(yàn)證點(diǎn)分布在直線y=x兩側(cè)而并沒(méi)有落在直線上。

已有學(xué)者通過(guò)試驗(yàn)得到，滴灌條件下濕潤(rùn)鋒運(yùn)移過(guò)程在統(tǒng)計(jì)學(xué)上是時(shí)間的冪函數(shù)[5,22-24]。也有不少學(xué)者根據(jù)水量平衡原理以及地表積水的原理，建立了以滴頭流量與土壤入滲率比為基礎(chǔ)的積水半徑數(shù)學(xué)模型[7,25]，該模型認(rèn)為滴灌條件下積水半徑與滴頭流量和土壤入滲率之間比值的0.5次方成正比。趙曄等[26]以水量平衡原理為基礎(chǔ)，引用考斯加科夫土壤入滲公式建立點(diǎn)源積水半徑運(yùn)移的數(shù)學(xué)模型，發(fā)現(xiàn)積水半徑隨積水時(shí)間呈冪函數(shù)增大。另外，張勇勇等[27]利用HYDRUS-2D模擬不同土壤物理性質(zhì)(土壤類型、初始含水率)、耕作技術(shù)參數(shù)和灌水技術(shù)參數(shù)組合的壟溝灌土壤水分入滲過(guò)程，采用空間矩分析方法定量分析不同因素對(duì)入滲濕潤(rùn)體特征量的影響，發(fā)現(xiàn)初始含水率對(duì)濕潤(rùn)體特征量的影響較其他因素小，土壤質(zhì)地對(duì)濕潤(rùn)體特征量的影響差異明顯。白雪兒等[28]對(duì)不同初始土壤含水率和滴頭流量下滴灌土壤濕潤(rùn)體特征研究發(fā)現(xiàn)，再分布后的濕潤(rùn)體體積主要受灌水量的影響，可以選擇較小的初始含水率及較大的滴頭流量以提高濕潤(rùn)體內(nèi)水分有效性。范嚴(yán)偉等[29]運(yùn)用HYDRUS-2D 對(duì)水平微潤(rùn)灌濕潤(rùn)體模擬進(jìn)行分析研究發(fā)現(xiàn)，土壤質(zhì)地對(duì)濕潤(rùn)體特性有顯著影響，土壤質(zhì)地越黏重，濕潤(rùn)鋒運(yùn)移速率越慢，濕潤(rùn)體體積越小，其“圓心”越靠近微潤(rùn)管，同時(shí)還發(fā)現(xiàn)在確定的土壤質(zhì)地條件下，初始含水率和壓力水頭對(duì)濕潤(rùn)體特性有較大影響，濕潤(rùn)鋒運(yùn)移距離及濕潤(rùn)體體積均隨土壤初始含水率、壓力水頭的增大而增大。

綜上所述，積水半徑、土壤質(zhì)地、土壤初始含水率和滴頭流量是影響濕潤(rùn)體的主要因素，本文基于積水半徑為1 cm的數(shù)值模擬結(jié)果較好，建立了水平濕潤(rùn)鋒長(zhǎng)度、濕潤(rùn)鋒比與土壤水力參數(shù)間的關(guān)系，通過(guò)此關(guān)系可以估算土壤濕潤(rùn)體范圍，在一定程度上可以為分析點(diǎn)源積水運(yùn)移問(wèn)題提供科學(xué)依據(jù)。但由于數(shù)值模擬時(shí)假定土壤質(zhì)地均勻，初始含水率等在一定理想條件下進(jìn)行，與實(shí)際自然環(huán)境狀態(tài)下的土壤結(jié)構(gòu)存在一定程度的差異，入滲受到土壤空間異質(zhì)性、作物生長(zhǎng)以及氣象因素的影響比較復(fù)雜，因此所得結(jié)論還需通過(guò)室內(nèi)和田間試驗(yàn)進(jìn)一步驗(yàn)證完善，同時(shí)在后期工作中會(huì)對(duì)積水半徑、初始含水率的影響以及同時(shí)適用于細(xì)質(zhì)土壤和粗質(zhì)土壤的幾何模型等進(jìn)一步深入研究，以便于對(duì)模型進(jìn)行優(yōu)化和完善。

4 結(jié)論

(1)所建數(shù)學(xué)模型濕潤(rùn)鋒比與進(jìn)氣吸力呈線性遞增關(guān)系，模型簡(jiǎn)單、便于應(yīng)用，且冪函數(shù)系數(shù)和指數(shù)可以分別用土壤飽和擴(kuò)散率的一次多項(xiàng)式和二次多項(xiàng)式進(jìn)行估計(jì)，決定系數(shù)R2均大于0.980，RMSE在0.004～0.017之間，RE在0.011%～0.013%之間，推導(dǎo)得到的濕潤(rùn)鋒模型表征了水平濕潤(rùn)鋒長(zhǎng)度與飽和擴(kuò)散率間的關(guān)系，經(jīng)驗(yàn)證具有較好的效果，該模型揭示出土壤水平濕潤(rùn)鋒長(zhǎng)度與土壤水力參數(shù)具有密切關(guān)系。通過(guò)此模型發(fā)現(xiàn)，隨著時(shí)間的推移，水平濕潤(rùn)鋒增加速率在1 000 min左右開始趨于緩慢，接近于一個(gè)穩(wěn)定值。

(2)建立的模型能較好地表征地表二維積水入滲規(guī)律，細(xì)質(zhì)土壤濕潤(rùn)體是以積水中心為圓心、以垂直濕潤(rùn)鋒和水平濕潤(rùn)鋒分別為長(zhǎng)半軸和短半軸的半橢圓面。粗質(zhì)土壤濕潤(rùn)體也可以近似看成橢圓面，但粗質(zhì)土壤濕潤(rùn)體的圓心位于積水中心下方，因此，短半軸大于水平濕潤(rùn)鋒，長(zhǎng)半軸小于垂直濕潤(rùn)鋒。