垂直線源灌土壤濕潤體尺寸預(yù)測模型研究

范嚴(yán)偉 邵曉霞 王 英 龔家國

(1.蘭州理工大學(xué)能源與動力工程學(xué)院， 蘭州 730050； 2.中國水利水電科學(xué)研究院流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室， 北京 100038)

0 引言

我國西北旱區(qū)土地資源豐富，日照充足，晝夜溫差大，是發(fā)展林果業(yè)的理想?yún)^(qū)域[1]。然而，這些地區(qū)降雨量有限，果樹生產(chǎn)在很大程度上取決于灌溉[2]。傳統(tǒng)灌溉方式耗水量大、水分利用效率低，不利于生態(tài)經(jīng)濟(jì)的可持續(xù)發(fā)展[3-4]。因此，灌溉技術(shù)和水資源管理的改進(jìn)將發(fā)揮重要作用。

垂直線源灌是一種適用于深根植物的節(jié)水灌溉方法，其灌水器垂直埋入土體，灌水過程中，水分直接進(jìn)入植物根部，濕潤體不易觀測[5]。濕潤體的形狀及大小影響著植物的生長與產(chǎn)量，了解垂直線源灌濕潤體動態(tài)變化特征，可確保灌水器在活性根區(qū)的精確放置，對設(shè)計經(jīng)濟(jì)高效的垂直線源灌溉系統(tǒng)至關(guān)重要[6-7]。土壤質(zhì)地是決定灌溉設(shè)計參數(shù)的重要因素，地下灌溉系統(tǒng)的設(shè)計應(yīng)考慮土壤質(zhì)地的影響[8-12]；相同土壤質(zhì)地條件下，土壤容重增加，孔隙度減小，土體入滲能力下降[13-16]；土壤初始含水率決定了滲透初期的土壤水勢，土壤初始含水率增加，濕潤體尺寸逐漸增大[5,11,17-20]。線源長度和直徑?jīng)Q定了灌水器滲水界面的大小，線源長度或直徑增大，滲水界面面積變大，意味著水分進(jìn)入土壤的通道增加，導(dǎo)致相同時間內(nèi)入滲水量增多，濕潤體隨之增大[21-24]。因此，從實(shí)用角度出發(fā)，線源直徑既要盡可能的大以加快其滲水速率，提高灌水均勻度，又要盡可能的小以減弱對作物根系生長的影響。灌水器埋深直接改變濕潤體水分分布位置，是實(shí)現(xiàn)作物根系與濕潤體有效匹配的關(guān)鍵因素，埋深過淺會增加地表水分無效蒸發(fā)，埋深過深又會引起深層滲漏和表土水分虧缺[25-29]。因此，埋深應(yīng)與土壤條件、根系分布及耕作要求等相適應(yīng)。灌溉必須適時適量，灌水時間過早或過晚、灌水定額過大或過小都是無益的。DU等[2]研究表明，中國西北干旱區(qū)蘋果根區(qū)土壤含水率低于田間持水率的50%～55%時，會對樹木生長和最終產(chǎn)量造成水分脅迫；周罕覓等[30]研究水肥耦合對3年生蘋果幼樹生長、產(chǎn)量、品質(zhì)及水肥利用的效應(yīng)，得出灌水下限為田間持水率的65%～75%；賈俊杰等[31]指出，SH矮砧蘋果幼樹滴灌條件下適宜灌水下限為60%的田間持水率。孫三民等[32]通過小區(qū)試驗(yàn)，確定13 L/(棵·次)的灌水量為適宜的新疆紅棗間接地下滴灌灌溉模式；張陸軍等[33]指出，陜北山地梨棗樹涌泉根灌時，每株2個灌水器及每個灌水器40 L/(株·次)的灌水量組合是適宜的布置方式；吳悠等[34]通過遮雨棚下可稱量式蒸滲桶試驗(yàn)得出，生育期內(nèi)8.4 L/(株·次)為柱狀蘋果樹相對節(jié)水的灌溉模式。

濕潤體動態(tài)變化可通過濕潤鋒距離量化表征。國內(nèi)外學(xué)者開發(fā)了一些用于確定濕潤鋒距離的模型，其中，最常見的是分析模型[21,35-38]、數(shù)值模型[39-40]和經(jīng)驗(yàn)?zāi)Ｐ蚚11,41-45]。通常，通過求解特定初始和邊界條件的控制方程(Richards方程)來開發(fā)分析模型和數(shù)值模型，而使用實(shí)驗(yàn)或數(shù)值模擬的回歸分析來開發(fā)經(jīng)驗(yàn)?zāi)Ｐ?。文獻(xiàn)[46-48]對數(shù)值和經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行了比較和評估，研究表明，HYDRUS模型計算結(jié)果能較好地反映土壤水分運(yùn)動基本規(guī)律，但模型較復(fù)雜，需輸入大量參數(shù)才能模擬計算；另外，每個經(jīng)驗(yàn)?zāi)Ｐ投际峭寥浪μ匦院凸喔葏?shù)的函數(shù)方程，形式較簡單，但僅適用于具體的灌溉技術(shù)，如開發(fā)的滴灌或溝灌濕潤體預(yù)測模型并不適用于垂直線源灌。因此，有必要開發(fā)一種可以預(yù)測垂直線源灌土壤濕潤體尺寸的經(jīng)驗(yàn)?zāi)Ｐ?，為確定適宜的灌水技術(shù)參數(shù)和實(shí)現(xiàn)灌溉系統(tǒng)優(yōu)化運(yùn)行提供實(shí)用而方便的手段。

數(shù)值模擬方法可對不同土壤特性、不同灌水器規(guī)格和不同設(shè)計參數(shù)條件下的土壤水分運(yùn)動過程進(jìn)行模擬[49-51]。李淑芹等[52]、FAN等[53]通過試驗(yàn)驗(yàn)證了垂直線源灌HYDRUS-2D模擬結(jié)果的有效性。基于此，本文采用HYDRUS-2D軟件，模擬研究土壤質(zhì)地、初始含水率、線源長度、線源直徑和線源埋深對垂直線源灌濕潤體運(yùn)移特征值的影響；利用模擬數(shù)據(jù)篩選影響濕潤體運(yùn)移的主導(dǎo)因素，進(jìn)而構(gòu)建預(yù)測濕潤體尺寸的簡化經(jīng)驗(yàn)?zāi)Ｐ?；最后，通過土箱試驗(yàn)驗(yàn)證經(jīng)驗(yàn)?zāi)Ｐ偷目煽啃浴?/p>

1 材料與方法

1.1 試驗(yàn)設(shè)計

試驗(yàn)裝置由3部分組成：土箱、馬氏瓶和線源灌水器，如圖1所示。土箱由10 mm厚有機(jī)玻璃制成，長×寬×高為50 cm×50 cm×100 cm。土箱底部留有多個通氣孔(直徑2 mm)，以防氣阻發(fā)生。線源與土箱接觸面開取土孔(直徑2 cm，間距5 cm)，用于測量灌溉結(jié)束時的土壤含水率。線源采用1/4圓柱體，底端密封，管底向上l長度的柱面均勻開孔。馬氏瓶直徑為10 cm，高度為100 cm。試驗(yàn)前，將供試土樣按設(shè)定的初始含水率加水，均勻混合后，用塑料薄膜密閉靜置1 d，待土壤水分分布均勻后，按設(shè)計容重分層(5 cm)裝入土箱，以獲得均勻土壤剖面。為了便于觀察土壤濕潤體變化過程，將線源灌水器用紗布包裹，并置于土箱一角，確保灌水器管壁與土壤緊實(shí)接觸，待次日進(jìn)行入滲試驗(yàn)。試驗(yàn)中，馬氏瓶提供恒定水頭，按先密后疏的時間間隔記錄累積入滲量，并用馬克筆繪制濕潤鋒運(yùn)移圖。入滲達(dá)到設(shè)定灌水定額后停止供水，迅速從灌水器兩側(cè)預(yù)留孔取土，用干燥法(105℃干燥24 h)測定土壤含水率。為盡量消除試驗(yàn)誤差，每個試驗(yàn)重復(fù)3次。

圖1 垂直線源灌試驗(yàn)裝置及灌水器細(xì)部結(jié)構(gòu)Fig.1 Experimental equipment for vertical line source irrigation and detailed structure of emitter1.調(diào)節(jié)水頭支架 2.馬氏瓶 3.灌水器 4.橡膠管 5.滲水孔6.土箱 7.取土孔 8.通氣孔 9.土壤表面 10.灌水器細(xì)部結(jié)構(gòu)

參照文獻(xiàn)[2,30-34]研究成果，取民勤地區(qū)砂壤土(容重γd=1.45 g/cm3，田間持水率θf=0.332 cm3/cm3，飽和導(dǎo)水率Ks=0.039 cm/min)和風(fēng)沙土(γd=1.56 g/cm3，θf=0.051 cm3/cm3，Ks=0.345 cm/min)，每種土壤采用2種處理(初始含水率θ0=60%θf、線源直徑d=4 cm、線源長度l=20 cm、線源埋深b=40 cm、灌水量V=40 L；θ0=70%θf、d=6 cm、l=30 cm、b=50 cm、V=40 L)進(jìn)行垂直線源灌土壤入滲試驗(yàn)。

1.2 數(shù)學(xué)建模

1.2.1基本方程

假設(shè)土壤是均勻和各向同性的，垂直線源灌可概念化為軸對稱的三維入滲過程。使用HYDRUS-2D模擬[54]。土壤水分運(yùn)動控制方程為Richards方程

(1)

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

z——垂向坐標(biāo)值，cm，向下為正

θ——土壤含水率，cm3/cm3

h——壓力水頭，cm

t——時間，min

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

采用van Genuchten-Mualem(VG-M)方程[55-56]描述土壤水分特征曲線和非飽和導(dǎo)水率

(2)

(3)

其中

m=1-1/n

式中Se——土壤相對飽和度

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

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

α——與進(jìn)氣值成反比的經(jīng)驗(yàn)參數(shù)，cm-1

n、m——影響土壤水分特征曲線形狀的經(jīng)驗(yàn)常數(shù)

1.2.2定解條件

圖2(圖中A、B、C分別為線源最高點(diǎn)、最低點(diǎn)和中心點(diǎn))為本研究中用于模擬不同建模情景的初始和邊界條件。

圖2 具有初始和邊界條件的計算域Fig.2 Computational domain with initial and boundary conditions

在所有的模擬情景中，土壤含水率按初始含水率設(shè)置；上邊界DE受大氣條件影響，考慮灌水過程中地表為干土層，蒸發(fā)量很小，為簡化計算，按零通量面設(shè)置；下邊界FG不受灌水影響，為自由排水，按零通量面設(shè)置；左邊界GH為灌水器中心軸，AD為塑料管壁，均無水量交換，按零通量面設(shè)置；右邊界EF灌溉水未到達(dá)，按零通量面設(shè)置；線源底部BH密封，為零通量邊界；滲水面邊界為充分供水方式，供水開始后很快達(dá)到飽和，可按定水頭邊界處理[20,22]。

綜上，初始條件可表述為

(4)

式中θ0——土壤初始含水率，cm3/cm3

Ω——計算域(圖2)

邊界條件可表述為

(5)

1.2.3模擬方案

采用單因素分析法，設(shè)置81種情景，模擬分析不同土壤質(zhì)地(表1)、θ0(50%θf、60%θf、70%θf)、d(2、4、6 cm)、l(10、20、30 cm)和b(30、40、50 cm)等因素對垂直線源灌濕潤體的影響。土壤質(zhì)地VG-M模型參數(shù)取自CARSEL等[57]資料以及文獻(xiàn)[52]，如表1所示。

表1 HYDRUS模擬中9種典型土壤的VG-M模型參數(shù)Tab.1 VG-M model parameters of nine typical soils in HYDRUS simulation

注：*表示取自文獻(xiàn)[57]，** 表示取自文獻(xiàn)[52]。

不同質(zhì)地土壤田間持水率采用RAB等[58]建立的預(yù)測模型獲得，具體表達(dá)式為

(6)

式中θp——凋萎系數(shù)，可采用VG-M模型參數(shù)中的θr表示[11]，cm3/cm3

1.2.4求解方法

利用HYDRUS-2D進(jìn)行數(shù)值求解。求解過程中，采用隱式差分格式進(jìn)行時間離散，Galerkin有限元法對土壤剖面進(jìn)行空間離散。考慮到田間實(shí)際和計算精度的要求，確定有限單元計算域深度為100 cm，寬度為50 cm，空間步長為1 cm，時間步長為0.1 min，模擬歷時由灌水定額(40 L)決定。

1.3 分析方法

垂直線源灌濕潤體形狀近似為“梨”型[22,52]。選取5個特征值(A點(diǎn)水平方向、B點(diǎn)水平方向、C點(diǎn)水平方向、C點(diǎn)垂直向上和C點(diǎn)垂直向下)勾畫出濕潤體輪廓，點(diǎn)A、B和C見圖2。研究表明地下三維入滲土壤濕潤鋒運(yùn)移過程可采用冪函數(shù)描述，且具有很高的精度[22,59-61]。因此，采用冪函數(shù)定量分析垂直線源灌土壤濕潤鋒運(yùn)移過程，其具體表達(dá)式為

UC=l/2+b1ta1

(7)

DC=l/2+b2ta2

(8)

RA=d/2+b3ta3

(9)

RB=d/2+b4ta4

(10)

RC=d/2+b5ta5

(11)

式中UC、DC——C點(diǎn)垂直向上、垂直向下濕潤高度，cm

RA、RB、RC——A、B、C點(diǎn)濕潤半徑，cm

a1、a2、a3、a4、a5、b1、b2、b3、b4、b5——擬合參數(shù)

根據(jù)模擬結(jié)果分析濕潤體運(yùn)移規(guī)律，探討影響機(jī)理，篩選主導(dǎo)因素，采用式(7)～(11)擬合獲得a1、a2、a3、a4、a5、b1、b2、b3、b4和b5值，研究擬合參數(shù)與主導(dǎo)因素間的量化關(guān)系，進(jìn)而建立垂直線源灌土壤濕潤體運(yùn)移距離模型。

1.4 誤差分析

選取4個指標(biāo)，即平均絕對誤差MAE、均方根誤差RMSE、偏差百分比PBIAS和納什效率系數(shù)NSE，對濕潤體尺寸的實(shí)測值和預(yù)測值進(jìn)行誤差分析。指標(biāo)參數(shù)定義為

(12)

(13)

(14)

(15)

式中EMAE——平均絕對誤差

EPBIAS——偏差百分比

ENSE——納什效率系數(shù)

Mi——第i個實(shí)測值

Si——第i個模擬值

Mmean——實(shí)測值的平均值

N——數(shù)據(jù)總個數(shù)

MAE和RMSE的數(shù)值越接近0，PBIAS為[-10，10]，NSE越靠近1，表示模擬值與實(shí)測值差異越小，兩者吻合越好[62]。

學(xué)生又活躍起來，當(dāng)tan (α2-α1)不存在時，應(yīng)是右邊分母為零時，即1+k1k2=0，從而k1k2=-1．也可以這樣說，如果1+k1k2≠0，則tan (α2-α1)有意義，而tan (α2-α1)沒意義，所以必有1+k1k2=0．

2 結(jié)果與分析

2.1 擬合參數(shù)a1～a5影響因素分析

2.1.1擬合參數(shù)a1

利用HYDRUS-2D模擬結(jié)果，采用式(7)，擬合獲得不同影響因素下a1，如圖3所示。

圖3 擬合參數(shù)a1隨飽和導(dǎo)水率Ks的變化規(guī)律Fig.3 Variations of fitting parameter a1 with saturated hydraulic conductivity Ks

由圖3可見，相同Ks時，擬合參數(shù)a1主要受線源長度和埋深的影響，而土壤初始含水率和線源直徑對其影響較小。主要是線源長度和埋深影響濕潤鋒到達(dá)地表的時間，線源越長或埋深越淺，濕潤鋒到達(dá)地表時間越短，導(dǎo)致擬合參數(shù)a1時數(shù)據(jù)量減少，產(chǎn)生了擬合誤差。不同Ks情況下，擬合參數(shù)a1隨Ks先增大后減小，但增減幅度不大(0.172～0.355)，為簡化計算，可取平均值，即a1=0.3。

2.1.2擬合參數(shù)a2

利用HYDRUS-2D模擬結(jié)果，采用式(8)，擬合獲得不同影響因素下a2，如圖4所示。

圖4 擬合參數(shù)a2隨飽和導(dǎo)水率Ks的變化規(guī)律Fig.4 Variations of fitting parameter a2 with saturated hydraulic conductivity Ks

圖5 擬合參數(shù)a3隨飽和導(dǎo)水率Ks的變化規(guī)律Fig.5 Variations of fitting parameter a3 with saturated hydraulic conductivity Ks

2.1.3擬合參數(shù)a3

利用HYDRUS-2D模擬結(jié)果，采用式(9)，擬合獲得不同影響因素下a3，如圖5所示。

由圖5可見，土壤初始含水率、線源直徑、線源長度和埋深對擬合參數(shù)a3影響較小。不同Ks情況下，擬合參數(shù)a3隨Ks先增大后減小，但增減幅度不大(0.269～0.381)，為簡化計算，可取平均值，即a3=0.348。

2.1.4擬合參數(shù)a4

利用HYDRUS-2D模擬結(jié)果，采用式(10)，擬合獲得不同影響因素下a4，如圖6所示。

圖6 擬合參數(shù)a4隨飽和導(dǎo)水率Ks的變化規(guī)律Fig.6 Variations of fitting parameter a4 with saturated hydraulic conductivity Ks

由圖6可見，相同Ks時，擬合參數(shù)a4主要受線源長度的影響，而土壤初始含水率、線源直徑和埋深對其影響較小。線源越長，其滲水速率越快，導(dǎo)致濕潤體水分疊加效應(yīng)增強(qiáng)，同時，由于重力勢的作用，下部濕潤體的疊加效應(yīng)強(qiáng)于上部。不同Ks情況下，擬合參數(shù)a4隨Ks先增大后減小，但增減幅度不大(0.298～0.414)，為簡化計算，可取平均值，即a4=0.374。

2.1.5擬合參數(shù)a5

利用HYDRUS-2D模擬結(jié)果，采用式(11)，擬合獲得不同影響因素下a5，如圖7所示。

由圖7可見，類似于擬合參數(shù)a4，相同Ks條件下，擬合參數(shù)a5主要受線源長度的影響，不同Ks情況下，擬合參數(shù)a5隨Ks先增大后減小，但增減幅度不大(0.278～0.379)，為簡化計算，可取平均值，即a5=0.346。

圖7 擬合參數(shù)a5隨飽和導(dǎo)水率Ks的變化規(guī)律Fig.7 Variations of fitting parameter a5 with saturated hydraulic conductivity Ks

綜上所述，擬合參數(shù)a1、a3、a4和a5隨Ks、θ0、d、l和b的變化規(guī)律不明顯，且變化幅度較小。為簡化計算，分別取其平均值。a2隨θ0、d、l和b的增減而稍有變化，但變化較小，而隨Ks的增大而增大，兩者具有較好的冪函數(shù)關(guān)系。

將a1～a5分別代入式(7)～(11)，得

UC=l/2+b1t0.300

(16)

(17)

RA=d/2+b3t0.348

(18)

RB=d/2+b4t0.374

(19)

RC=d/2+b5t0.346

(20)

2.2 擬合參數(shù)b1～b5影響因素分析

利用HYDRUS-2D模擬結(jié)果，采用式(16)～(20)，再次擬合獲得不同θ0、d、l和b組合下的b1～b5，如圖8所示。

由圖8可見，擬合參數(shù)b1、b2、b3、b4和b5均隨Ks的增大而增大，θ0、d、l和b對其影響相對較小。進(jìn)一步分析發(fā)現(xiàn)b1、b3、b4和b5與Ks具有很好的冪函數(shù)關(guān)系，而b2與Ks呈線性關(guān)系?；诖耍傻么怪本€源灌濕潤體尺寸簡化預(yù)測模型，即

(21)

(22)

(23)

(24)

(25)

2.3 模型驗(yàn)證

利用試驗(yàn)數(shù)據(jù)對簡化預(yù)測模型進(jìn)行驗(yàn)證。將砂壤土和風(fēng)沙土的Ks分別代入式(21)～(25)，得2種土壤的垂直線源灌濕潤體尺寸簡化預(yù)測模型為：

民勤砂壤土

(26)

民勤風(fēng)沙土

(27)

將簡化預(yù)測模型計算值與試驗(yàn)特征值(2個處理，3個重復(fù))進(jìn)行對比分析，如圖9所示。

圖8 擬合參數(shù)b1～b5隨飽和導(dǎo)水率Ks的變化規(guī)律Fig.8 Variations of fitting parameter b1～b5 with saturated hydraulic conductivity Ks

圖9 2種土壤垂直線源灌土壤濕潤體尺寸計算值與實(shí)測值對比Fig.9 Comparisons of calculated and measured values of wetted soil dimensions for two soil types undervertical line source irrigation

采用式(12)～(15)，對計算值與實(shí)測值進(jìn)行統(tǒng)計分析，計算結(jié)果見表2。

由表2可知，MAE和RMSE接近0，PBIAS為-4%～9%之間，NSE靠近1(NSE不小于0.929)，說明簡化預(yù)測模型計算值與實(shí)測值一致性良好，但仍存在一定誤差，究其原因可能是垂直線源灌土壤濕潤鋒運(yùn)移距離受土壤質(zhì)地影響最大，而土壤初始含水率以及線源直徑、長度和埋深對其尚有一些影響，為了簡化計算，僅考慮了土壤質(zhì)地的影響，建立了單變量模型，從而在一定程度上影響了計算結(jié)果的準(zhǔn)確性。另外，僅采用飽和導(dǎo)水率Ks來表征不同土壤質(zhì)地濕潤鋒運(yùn)移規(guī)律也是存在部分誤差的原因之一。

表2 計算值與實(shí)測值統(tǒng)計分析Tab.2 Statistical analysis of calculated and measured values

3 結(jié)論

(1)垂直線源灌濕潤體尺寸主要受土壤質(zhì)地影響，土壤質(zhì)地越粗(Ks越大)，濕潤鋒運(yùn)移越快；線源長度、線源直徑和埋深對其影響較小。

(2)土壤濕潤鋒運(yùn)移過程符合冪函數(shù)關(guān)系，冪函數(shù)系數(shù)隨Ks的增大而增大，冪函數(shù)指數(shù)在垂直向上和水平方向上變化較小，而在垂直向下方向上隨Ks的增大而增大。

(3)提出了包含Ks的垂直線源灌土壤濕潤體尺寸預(yù)測模型，利用試驗(yàn)驗(yàn)證了預(yù)測模型的有效性，初步實(shí)現(xiàn)了由土壤物理參數(shù)預(yù)測垂直線源灌土壤濕潤鋒運(yùn)移距離的可能。