劉文光， 賈生海， 范嚴(yán)偉， 白貴林， 趙 彤

(1.甘肅農(nóng)業(yè)大學(xué) 水利水電工程學(xué)院， 甘肅 蘭州 730070； 2.蘭州理工大學(xué) 能源與動(dòng)力工程學(xué)院， 甘肅 蘭州 730050)

1 研究背景

有機(jī)融合農(nóng)業(yè)耕作措施與節(jié)水灌溉技術(shù)形成綜合的節(jié)水技術(shù)體系，是節(jié)水農(nóng)業(yè)關(guān)鍵技術(shù)的突破點(diǎn)，受到農(nóng)業(yè)和土壤相關(guān)科研工作者的重視[1]。研究表明，夾砂層的存在可減小土壤水分下滲速率，提高上層土壤的水分含量[2-5]?；趭A砂土壤中水分的入滲特點(diǎn)，在黃土中設(shè)置砂層的耕作方法得到認(rèn)可并應(yīng)用[6-7]。

目前，夾砂土壤入滲機(jī)理研究主要是針對(duì)畦灌和溝灌灌水模式，取得了豐富的成果，為夾砂土壤應(yīng)用提供了理論指導(dǎo)[8-13]。隨著地膜覆蓋技術(shù)的發(fā)展，膜孔灌逐漸發(fā)展起來，其利用地膜輸水并通過作物孔和專用孔入滲進(jìn)行灌溉，具有節(jié)水、保肥、灌水均勻度高等優(yōu)點(diǎn)，在中國(guó)西北干旱半干旱區(qū)得到應(yīng)用[14-15]。以往，膜孔入滲特性研究主要是針對(duì)均質(zhì)土壤，對(duì)夾砂土壤入滲特性研究相對(duì)較少。面對(duì)新形勢(shì)下農(nóng)田不同節(jié)水灌溉技術(shù)，研究膜孔灌水模式下夾砂層土壤水分運(yùn)動(dòng)規(guī)律具有重要的現(xiàn)實(shí)意義，同時(shí)可為水資源高效利用提供理論依據(jù)。

隨著計(jì)算機(jī)技術(shù)的發(fā)展，采用數(shù)值模擬方法進(jìn)行土壤水分運(yùn)動(dòng)規(guī)律的研究成為趨勢(shì)。由Simunek等[16]開發(fā)的SWMS軟件被廣泛用于模擬各種灌溉方式的土壤水分運(yùn)動(dòng)，并且取得了良好的效果。聶衛(wèi)波等[17]采用SWMS-1D軟件對(duì)均質(zhì)土一維土壤入滲特性進(jìn)行了數(shù)值模擬；范嚴(yán)偉等[18]利用SWMS-2D軟件對(duì)層狀土垂直一維入滲土壤水分運(yùn)動(dòng)規(guī)律進(jìn)行了數(shù)值模擬；Li等[19]利用SWMS-2D模型軟件對(duì)層狀土壤膜孔入滲量進(jìn)行了研究，但缺少必要的驗(yàn)證。

本文以室內(nèi)試驗(yàn)為基礎(chǔ)，針對(duì)土壤剖面存在砂質(zhì)夾層的特點(diǎn)，基于Richard方程建立膜孔灌條件下夾砂層土壤水分運(yùn)動(dòng)數(shù)學(xué)模型，利用SWMS-2D軟件求解，通過室內(nèi)試驗(yàn)對(duì)模擬結(jié)果進(jìn)行驗(yàn)證，分析數(shù)值模擬的夾砂土膜孔灌土壤水分運(yùn)動(dòng)規(guī)律的合理性及可靠性。以期借助數(shù)值模擬方法，進(jìn)一步認(rèn)識(shí)膜孔灌條件下夾砂土壤水分運(yùn)動(dòng)機(jī)理，為后期模擬分析土壤特性和灌水技術(shù)要素對(duì)夾砂土體中膜孔入滲特性的影響提供技術(shù)支撐。

2 材料與方法

2.1 室內(nèi)試驗(yàn)

2.1.1 試驗(yàn)土樣 供試壤土取自蘭州市黃峪鎮(zhèn)祁家營(yíng)村蘋果園內(nèi)；砂土取自蘭州市秦王川盆地。取土深度為5～60 cm，將土樣帶回實(shí)驗(yàn)室，經(jīng)風(fēng)干、碾壓，過2 mm篩后制成試驗(yàn)土樣，并采用激光粒度儀(MS2000)測(cè)量其顆粒級(jí)配。

參照國(guó)際制土壤分類標(biāo)準(zhǔn)，試驗(yàn)土樣屬于粉壤土和砂土，土樣的基本物理特性見表1。

表1 試驗(yàn)土壤基本物理特性

2.1.2 試驗(yàn)方法 試驗(yàn)裝置主要由馬氏瓶和有機(jī)玻璃土箱組成，如圖1所示，在土箱一角安裝有膜孔室，并在側(cè)面設(shè)置取土孔和底部設(shè)通氣孔。將制備好的土樣依照試驗(yàn)方案(見表2)，以粉壤土容重為1.33 g/cm3和砂土為1.55 g/cm3分層(5 cm)裝入30 cm×30 cm×60 cm(長(zhǎng)×寬×高)的有機(jī)玻璃土箱。試驗(yàn)中，用馬氏瓶保持恒定水頭供水，按已確定的時(shí)間間隔記錄馬氏瓶讀數(shù)，并描濕潤(rùn)鋒運(yùn)移曲線，為防止水分蒸發(fā)在土箱上表面覆保鮮膜。試驗(yàn)結(jié)束后，從取土孔取土測(cè)定試驗(yàn)土體土壤含水量。

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

2.2.1 基本方程 假設(shè)每層試驗(yàn)土壤為理想的均質(zhì)且各向同性多孔介質(zhì)，不考慮蒸發(fā)、溫度以及土壤內(nèi)部空氣阻力對(duì)水分入滲的影響，在每層土壤具有相同的水分運(yùn)動(dòng)模型表達(dá)形式且有不同特性參數(shù)值的情況下，夾砂層土壤膜孔灌為充分供水條件下的三維入滲，可通過三維土壤水分運(yùn)動(dòng)基本方程的定解進(jìn)行數(shù)值模擬。

(1)

式中：x為橫向坐標(biāo)，cm；y為縱向坐標(biāo)，cm；z為垂向坐標(biāo)，向上為正，cm；C(φm)為比水容重，cm-1；φm為基質(zhì)勢(shì)，cm；φ為總水勢(shì)，cm；t為入滲時(shí)間，min；K(φm)為非飽和導(dǎo)水率，cm/min。

(2)

K(φm)=KsSe0.5[1-(1-Se1/m)m]2

(3)

式中：Se為土壤有效含水率，cm3/cm3；θ為土壤含水率，cm3/cm3；θr為土壤殘余含水率，cm3/cm3；θs為土壤飽和含水率，cm3/cm3；α為與進(jìn)氣值成反比的經(jīng)驗(yàn)參數(shù)，cm-1；n和m為影響土壤水分特征曲線形狀的經(jīng)驗(yàn)常數(shù)，m=1-1/n;Ks為土壤飽和導(dǎo)水率，cm/min。

表2 夾砂層土壤膜孔灌試驗(yàn)方案

2.2.2 定解條件 因膜孔具有軸對(duì)稱性，所以本試驗(yàn)可以簡(jiǎn)化為研究由1/4膜孔所控制的土壤水分入滲規(guī)律，簡(jiǎn)化后模型如圖2所示，OBCD-O′B′C′D′表示試驗(yàn)土體，其中，EFGH-E′F′G′H′表示夾砂層。由此確定的定解條件為：

(1)初始條件。

φ=φ0,t=0, 0≤x≤d, 0≤y≤d,

-h≤z≤0

(4)

式中：φ0為土壤的初始總水勢(shì)，cm；d為模擬土體的寬度，cm；h為模擬土體的深度，cm。

圖1 試驗(yàn)裝置圖

(2)邊界條件。在試驗(yàn)中，膜孔OAA′處水頭保持恒定；采用保鮮膜將上邊界AA′DCB覆蓋，因此沒有水分蒸發(fā)，所以為零通量面；由于膜孔的對(duì)稱性，1/4膜孔邊界面的OBB′O′和ODD′O′為零通量面；由于模擬土體體積較大，故可以認(rèn)為邊界面DCC′D′、BCC′B′和O′B′C′D′水分不會(huì)到達(dá)，即邊界無水量交換，總水勢(shì)會(huì)保持恒定。依上述分析，邊界條件為：

邊界OAA′：

(5)

邊界AA′DCB：

(6)

邊界OBB′O′：

y=0,t≥0

(7)

邊界ODD′O′：

x=0,t≥0

(8)

邊界DCC′D′：

φ=φ0, -h≤z≤0, 0≤x≤d,y=d,t≥0

(9)

邊界BCC′B′：

φ=φ0, -h≤z≤0, 0≤y≤d,x=d,t≥0

(10)

邊界O′B′C′D′：

φ=φ0, 0≤x≤d, 0≤y≤d,z=-h,t≥0

(11)

式中：h0為灌溉水深，cm；r為膜孔半徑，cm；T為灌水時(shí)間，min。

圖2 計(jì)算簡(jiǎn)化圖

2.2.3 數(shù)值求解方法 應(yīng)用軟件SWMS-2D對(duì)上述模型進(jìn)行求解。將模擬土體剖面劃分為方形計(jì)算單元，z軸方向間隔1 cm，x軸方向先密后疏，計(jì)算區(qū)域的深度為30 cm，半徑為30 cm。試驗(yàn)土壤的VG模型參數(shù)見表3。

2.3 誤差分析

采用以下4個(gè)指標(biāo)，用于分析模擬值和實(shí)測(cè)值之間的誤差，分別為：

(12)

(13)

(14)

(15)

式中：MAE為平均絕對(duì)誤差；RMSE為均方根誤差；PBIAS為偏差百分比；NSE為納什效率系數(shù)；Yiobs為第i個(gè)實(shí)測(cè)值；Yisim第i個(gè)模擬值；Ymean為實(shí)測(cè)值的平均值；n為數(shù)據(jù)總個(gè)數(shù)。MAE、RMSE和PBIAS的數(shù)值越接近0，NSE的數(shù)值越靠近1，表示模擬值與實(shí)測(cè)值差值越小，兩者吻合越好。

表3 試驗(yàn)土壤的Van Genuchten模型參數(shù)

3 結(jié)果分析

用SWMS-2D軟件求解膜孔灌條件下夾砂層土壤水分運(yùn)動(dòng)數(shù)學(xué)方程，獲得累積入滲量、濕潤(rùn)鋒運(yùn)移距離和土壤剖面含水量的模擬結(jié)果，將模擬結(jié)果與實(shí)測(cè)結(jié)果進(jìn)行對(duì)比驗(yàn)證。

3.1 累積入滲量的數(shù)值模擬與驗(yàn)證

圖3為不同處理下夾砂層土壤膜孔灌累積入滲量的模擬值與實(shí)測(cè)值的比較。

圖3 累積入滲量的模擬值與實(shí)測(cè)值對(duì)比

由圖3可以看出，膜孔直徑對(duì)累積入滲量有影響，膜孔直徑大者對(duì)應(yīng)的累積入滲量較大，這是因?yàn)槟た字睆皆酱?，水分進(jìn)入土壤的通道越大，導(dǎo)致其累積入滲量越大；夾砂層埋深對(duì)累積入滲量有一定的影響，夾砂層的埋深大者對(duì)應(yīng)的累積入滲量較大，主要是因?yàn)槿霛B水分達(dá)到砂-土界面向下運(yùn)動(dòng)的條件是砂土的進(jìn)水吸力等于界面處粉壤土的基質(zhì)勢(shì)，夾砂層埋深小者的入滲水分先到達(dá)砂-土界面，隨著入滲時(shí)間的增加，夾砂層埋深小者粉壤土的基質(zhì)勢(shì)增加程度較大，對(duì)應(yīng)的水分入滲速率較慢，累積入滲量較小，從而表現(xiàn)出一定的減滲現(xiàn)象。

對(duì)累積入滲量的模擬值與實(shí)測(cè)值進(jìn)行統(tǒng)計(jì)特征分析：MAE為0.05 L、RMSE為0.08 L、PBIAS為-0.81%、NSE為0.99。說明各時(shí)段數(shù)值模擬結(jié)果與實(shí)測(cè)值非常吻合，表明SWMS-2D能較好地模擬夾砂土壤膜孔灌累積入滲量的變化過程。個(gè)別處理實(shí)測(cè)值與模擬值相差較大，主要是為保持水頭恒定，在試驗(yàn)開始前，需用凡士林將塑料薄膜粘貼于膜孔處，然后向膜孔室內(nèi)注水至指定水位線，試驗(yàn)開始時(shí)再將薄膜去除，在這一操作過程中會(huì)產(chǎn)生水量損失(薄膜帶出水)，從而造成開始時(shí)間段內(nèi)的讀數(shù)不準(zhǔn)確；還有裝土不均勻和讀數(shù)誤差等也會(huì)對(duì)試驗(yàn)數(shù)據(jù)的準(zhǔn)確性有一定的影響。

3.2 濕潤(rùn)鋒運(yùn)移距離的模擬與驗(yàn)證

根據(jù)夾砂層土壤膜孔灌濕潤(rùn)鋒運(yùn)移的特殊性，分別進(jìn)行不同處理下濕潤(rùn)鋒水平向和垂直向運(yùn)移距離的實(shí)測(cè)值與模擬值的對(duì)比分析，如圖4所示。

分析圖4中濕潤(rùn)鋒的運(yùn)動(dòng)規(guī)律，可知，在其他條件相同時(shí)，膜孔直徑大者水平向濕潤(rùn)鋒運(yùn)移距離較大，夾砂層埋深小者水平向濕潤(rùn)鋒距離較大；每一種處理下水平向濕潤(rùn)鋒運(yùn)移距離明顯大于垂直向，這是因?yàn)閵A砂層的存在，在入滲水分到達(dá)砂-土交界面時(shí)，上層粉壤土的基質(zhì)勢(shì)小于夾砂層的進(jìn)水吸力，只有當(dāng)兩者相等時(shí)水分才會(huì)開始下滲，于是水分會(huì)在上層土壤中蓄積并繼續(xù)向水平方向運(yùn)移，而垂直向則運(yùn)移緩慢。

對(duì)濕潤(rùn)鋒運(yùn)移距離的模擬值與實(shí)測(cè)值進(jìn)行統(tǒng)計(jì)特征分析，模擬值與實(shí)測(cè)值的MAE為0.28 cm、RMSE為0.36 cm、PBIAS為-1.56%、NSE為0.99。說明各時(shí)段數(shù)值模擬結(jié)果與實(shí)測(cè)值一致性良好，表明SWMS-2D能較好地模擬夾砂土壤膜孔灌各個(gè)方向濕潤(rùn)鋒運(yùn)移的變化過程。由于在試驗(yàn)中，只有當(dāng)試驗(yàn)土體含水量達(dá)到一定量時(shí)，才可以觀測(cè)到濕潤(rùn)鋒，而且土壤中的禁錮氣體產(chǎn)生的壓力對(duì)入滲水流會(huì)有阻礙作用，所以，對(duì)濕潤(rùn)鋒運(yùn)移距離的模擬值進(jìn)行適當(dāng)減小。

3.3 土體剖面含水量模擬與驗(yàn)證

圖5為不同處理下夾砂層土壤膜孔灌濕潤(rùn)體剖面含水量的模擬值與實(shí)測(cè)值的比較，實(shí)線為surfer軟件利用模擬值繪制的土體剖面含水量等值線，圓點(diǎn)為實(shí)測(cè)值，原點(diǎn)0為膜孔中心，測(cè)點(diǎn)坐標(biāo)為(x，0，z)，單位為cm。

圖4 濕潤(rùn)鋒運(yùn)移距離的模擬值與實(shí)測(cè)值對(duì)比

圖5 濕潤(rùn)體剖面含水量的模擬值與實(shí)測(cè)值對(duì)比

分析圖5中濕潤(rùn)體水分分布結(jié)果不難看出，濕潤(rùn)體的濕潤(rùn)鋒面形狀與均質(zhì)土的半球體或半橢圓體不同，而是在粉壤土層與砂層的交界處出現(xiàn)了明顯的不連續(xù)性，表現(xiàn)出了一定的阻水性。在其他條件相同時(shí)，膜孔直徑對(duì)濕潤(rùn)體水分分布有影響，膜孔直徑大者對(duì)應(yīng)的濕潤(rùn)體體積較大；夾砂層埋深對(duì)濕潤(rùn)體水分分布有較大的影響，埋深小者入滲水分會(huì)穿過夾砂層到達(dá)下層土壤，水分會(huì)在每層土壤中重新分布，不與上層土壤的濕潤(rùn)鋒面連續(xù)。

取相應(yīng)點(diǎn)位的含水率模擬值與實(shí)測(cè)值進(jìn)行統(tǒng)計(jì)特征分析，模擬值與實(shí)測(cè)值的MAE為0.11 cm3/cm3、RMSE為0.05 cm3/cm3、PBIAS為-6.36%、NSE為0.97?？傮w上，模型計(jì)算出來的濕潤(rùn)體含水量的模擬值與實(shí)測(cè)值相近，表明SWMS-2D能較好地模擬獲得夾砂土壤膜孔灌濕潤(rùn)體水分分布狀況。個(gè)別試驗(yàn)點(diǎn)土壤含水量誤差較大，可能是由于土壤夯填不均勻，土鉆取土有誤差，且因取土量少造成稱量不準(zhǔn)確等誤差產(chǎn)生。

4 結(jié) 論

(1)累積入滲量、濕潤(rùn)鋒運(yùn)移距離和濕潤(rùn)體剖面含水量的模擬值與實(shí)測(cè)值的MAE分別為0.05 L、0.28 cm和0.11 cm3/cm3；RMSE分別為0.08 L、0.36 cm和0.05 cm3/cm3，PBIAS分別為-0.81%、-1.56%和-6.36%；NSE分別為0.99、0.99和0.97，說明，數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果具有較好的一致性。所建模型能較準(zhǔn)確地反映夾砂層土壤膜孔灌水分運(yùn)動(dòng)規(guī)律。

(2)數(shù)值模擬結(jié)果與試驗(yàn)結(jié)果均表明，夾砂層具有一定的阻水和減滲作用，能夠提高上層土壤的持水能力，而且使得不同處理下水平向濕潤(rùn)鋒運(yùn)移距離明顯大于垂直向。膜孔直徑和夾砂層埋深對(duì)累積入滲量、濕潤(rùn)鋒運(yùn)移距離和土體剖面含水量分布均有不同程度的影響：在其他條件相同時(shí)，膜孔直徑和夾砂層的埋深大者對(duì)應(yīng)的累積入滲量較大；膜孔直徑大者和夾砂層埋深小者水平向濕潤(rùn)鋒運(yùn)移距離較大；膜孔直徑大者對(duì)應(yīng)的濕潤(rùn)體體積較大；夾砂層埋深對(duì)濕潤(rùn)體水分分布有較大的影響，會(huì)使入滲水分不連續(xù)分布。

(3)本文利用SWMS-2D軟件對(duì)夾砂層土壤膜孔灌土壤水分入滲方程求解是可行的，模擬結(jié)果是可靠的。因此，SWMS-2D軟件可為進(jìn)一步研究不同夾砂層質(zhì)地、入滲水頭、初始含水率和夾砂層厚度等條件下的夾砂層土壤膜孔灌水分運(yùn)動(dòng)規(guī)律提供理論依據(jù)和技術(shù)手段。