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

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

0 引言

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

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

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

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

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

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

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

(1)

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

t——時間，min

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

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

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

1.1.2土壤水力模型

HYDRUS-3D是一款用于模擬二維、三維土壤水分運(yùn)動、溶質(zhì)運(yùn)移、根系吸水及熱量傳輸?shù)挠邢拊?jì)算軟件。本文選取HYDRUS-3D中二維土壤水分運(yùn)動模擬模塊，該模塊在進(jìn)行數(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的二維積水入滲過程，根據(jù)模擬得到的數(shù)據(jù)分析并構(gòu)建相關(guān)數(shù)學(xué)模型，土壤水力參數(shù)選取HYDRUS-3D軟件中自帶參數(shù)[12]，見表1，表中θ0為土壤初始含水率。

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

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

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

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

1.1.3模型設(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)。模擬初始時間為0 min，結(jié)束時間為2 000 min。設(shè)置最大迭代次數(shù)(Maximum number of iterations)為10 000次，設(shè)置初始時間步長(Initial time)0.5 min，最小時間步長(Minimum time)0.1 min，最大時間步長(Maximum time)1 min，含水率允許最大誤差限(Water content tolerance)0.001 cm3/cm3，壓力水頭允許最大誤差限(Pressure head tolerance)1 cm。采用HYDRUS-3D中FE-Mesh模塊對所建的二維模型進(jìn)行網(wǎng)格劃分，網(wǎng)格劃分的類型為三角形，有限元柵格劃分參數(shù)中設(shè)置其有限元三角形外接圓最大直徑為1 cm，其他參數(shù)采用HYDRUS數(shù)據(jù)庫中自帶數(shù)據(jù)。為了提高模擬結(jié)果的精度，對積水周圍的部分區(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的積水面，在不入滲時視地表為大氣邊界，入滲開始后很快達(dá)到飽和，視為定水頭邊界；入滲面EF左右兩側(cè)邊界，以及邊界AD和BC假定為零通量邊界，故其水分通量為零；下邊界CD為不受地下水位影響的自由出流邊界。在入滲開始時，即在t=0時刻，土壤含水率為土壤初始含水率θ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土壤濕潤體形狀模型

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

(8)

式中xf——水平濕潤鋒長度，cm

zf——垂直濕潤鋒長度，cm

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

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

xf=a1tb1

(9)

zf=a2tb2

(10)

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

土壤水分?jǐn)U散率反映了土壤孔隙度、孔隙大小分布以及導(dǎo)水性能，并影響土壤中水分運(yùn)動狀況，它是研究水鹽運(yùn)動規(guī)律與調(diào)控以及表征土壤水動力學(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ù)R2、均方根誤差(RMSE)和相對誤差(RE)評價指標(biāo)進(jìn)行誤差分析，檢驗(yàn)相關(guān)性。

2 結(jié)果與分析

2.1 土壤濕潤體形狀變化特征

圖2 土壤入滲濕潤體形狀Fig.2 Soil infiltration wetting pattern shapes

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

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

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

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

表3 濕潤鋒比隨時間的變化
Tab.3 Wetting frontratio varied with time

土壤類型時間/min水平濕潤鋒長度/cm垂直濕潤鋒長度/cm濕潤鋒比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——水平濕潤鋒和垂直濕潤鋒長度的比值，即濕潤鋒比

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

wt=0.004hd+0.787 (R2=0.999)

(14)

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

(15)

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

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

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

圖4 水平濕潤鋒長度隨時間的變化曲線Fig.4 Variation curves of horizontal wetting front with time

表4 水平濕潤鋒長度擬合結(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

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

(16)

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

圖5 水平濕潤鋒長度變化率隨時間的變化曲線Fig.5 Changing curves of variation rate of horizontal wetting front with time

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

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

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

土壤水力參數(shù)是求解土壤水分運(yùn)動方程的重要參數(shù)，也是建立土壤水分運(yùn)動數(shù)學(xué)模型的重要基礎(chǔ)。在研究水平濕潤鋒運(yùn)移隨時間變化規(guī)律的基礎(chǔ)上，運(yùn)用HYDRUS軟件模擬的數(shù)據(jù)進(jìn)行擬合研究發(fā)現(xiàn)，不同土壤入滲時的冪函數(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%；同時，擬合發(fā)現(xiàn)參數(shù)b1是土壤飽和擴(kuò)散率Ds的二次函數(shù)，如圖6b所示，決定系數(shù)R2=0.981，RMSE為0.004，RE為0.011%，數(shù)值模擬效果較好。

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

a1=0.098Ds+1.175

(17)

(18)

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

(19)

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

2.5 濕潤體變化過程定量關(guān)系評價

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

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

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

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

由圖4、5已經(jīng)得出，從1 000 min左右開始水平濕潤鋒長度變化速率基本趨于穩(wěn)定，隨著時間的推移，濕潤鋒運(yùn)移變化規(guī)律無限趨于線性變化，也就意味著從此時開始計(jì)算值略小于實(shí)測值，而在實(shí)際研究中入滲時長并不會很長，因此只對水平濕潤鋒運(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是對細(xì)質(zhì)土壤濕潤體形狀的驗(yàn)證，圖中實(shí)測值表示水分入滲的實(shí)際濕潤體輪廓線，擬合值表示通過式(8)計(jì)算得到的標(biāo)準(zhǔn)橢圓曲線。從圖中可以看出，在入滲時刻50、100、200、500 min的實(shí)際濕潤體形狀和計(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)二維積水入滲的土壤濕潤體形狀可近似看成一個以垂直濕潤鋒和水平濕潤鋒分別為長半軸和短半軸的半橢圓面，同時這一研究結(jié)果更加證明了其他學(xué)者[14-16]對土壤濕潤體形狀研究結(jié)論的準(zhǔn)確性。

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

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

3 討論

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

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

已有學(xué)者通過試驗(yàn)得到，滴灌條件下濕潤鋒運(yùn)移過程在統(tǒng)計(jì)學(xué)上是時間的冪函數(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ù)增大。另外，張勇勇等[27]利用HYDRUS-2D模擬不同土壤物理性質(zhì)(土壤類型、初始含水率)、耕作技術(shù)參數(shù)和灌水技術(shù)參數(shù)組合的壟溝灌土壤水分入滲過程，采用空間矩分析方法定量分析不同因素對入滲濕潤體特征量的影響，發(fā)現(xiàn)初始含水率對濕潤體特征量的影響較其他因素小，土壤質(zhì)地對濕潤體特征量的影響差異明顯。白雪兒等[28]對不同初始土壤含水率和滴頭流量下滴灌土壤濕潤體特征研究發(fā)現(xiàn)，再分布后的濕潤體體積主要受灌水量的影響，可以選擇較小的初始含水率及較大的滴頭流量以提高濕潤體內(nèi)水分有效性。范嚴(yán)偉等[29]運(yùn)用HYDRUS-2D 對水平微潤灌濕潤體模擬進(jìn)行分析研究發(fā)現(xiàn)，土壤質(zhì)地對濕潤體特性有顯著影響，土壤質(zhì)地越黏重，濕潤鋒運(yùn)移速率越慢，濕潤體體積越小，其“圓心”越靠近微潤管，同時還發(fā)現(xiàn)在確定的土壤質(zhì)地條件下，初始含水率和壓力水頭對濕潤體特性有較大影響，濕潤鋒運(yùn)移距離及濕潤體體積均隨土壤初始含水率、壓力水頭的增大而增大。

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

4 結(jié)論

(1)所建數(shù)學(xué)模型濕潤鋒比與進(jìn)氣吸力呈線性遞增關(guā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)得到的濕潤鋒模型表征了水平濕潤鋒長度與飽和擴(kuò)散率間的關(guān)系，經(jīng)驗(yàn)證具有較好的效果，該模型揭示出土壤水平濕潤鋒長度與土壤水力參數(shù)具有密切關(guān)系。通過此模型發(fā)現(xiàn)，隨著時間的推移，水平濕潤鋒增加速率在1 000 min左右開始趨于緩慢，接近于一個穩(wěn)定值。

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