基于圓筒入滲儀的田間膜孔灌土壤入滲參數研究

范嚴偉 唐興鵬 史金紅 馬天花

(蘭州理工大學能源與動力工程學院，蘭州 730050)

0 引言

地膜覆蓋是中國干旱、半干旱地區常用的栽培技術[1-2]。膜孔灌是在地膜覆蓋基礎上發展起來的一種地面灌溉技術，其灌溉過程中，水在地膜上傳輸并通過作物孔或專用灌水孔(統稱為膜孔)滲入根區土壤，具有節水、增溫、保墑等優點[3-4]。該灌溉方式投資少、效益高且簡便易行，在農業節水灌溉領域具有良好的應用前景。

膜孔入滲規律研究是優化田間膜孔灌灌水技術要素的基礎，而土壤入滲參數的準確估算是量化表征膜孔入滲規律的關鍵[5]。目前，對膜孔灌土壤入滲參數的研究主要集中在室內試驗[6-8]和數值模擬[9-11]兩方面。室內試驗能夠較好地控制灌水技術要素，如土壤質地、容重和初始含水率以及膜孔直徑和灌溉水頭等，是揭示膜孔入滲規律的基本方法，但室內試驗所用土壤多為擾動土且質地類型較單一，導致試驗推求的土壤入滲參數普適性不強；數值模擬能夠快速全面地獲取不同灌水要素組合下土壤水分運動規律，是定量分析膜孔入滲特性的有效手段，但基于數值模擬結果構建的相關估算模型還需試驗數據的驗證。為此，一些學者采用數值模擬和試驗驗證相結合方法，開展了膜孔灌土壤濕潤體特性研究，如介飛龍等[12]基于HYDRUS模擬結果，建立了考慮土壤初始含水率和特征(水平和垂向)濕潤距離的累積入滲量模型，并利用室內試驗進行了驗證。由于模型中包含了特征濕潤距離這一待測值，使其在田間實際應用中受到限制，因為田間膜孔灌土壤表面覆有地膜，其特征濕潤距離難以觀測；FAN等[13]在數值分析影響濕潤體運移的主導因素的基礎上，建立了包含土壤飽和導水率、膜孔直徑和灌水時間的特征濕潤距離估算模型。該模型的建立是基于HYDRUS軟件中自帶的土壤質地類型，未考慮土壤容重的影響，另外將土壤飽和導水率考慮其中，而田間測定原狀土飽和導水率相對較麻煩[14-16]。因此，有必要發展一種能夠對田間膜孔灌濕潤體尺寸進行現場評估的新模型、新方法。

為體現土壤質地的廣泛性和研究成果的普適性，本文依據國際制土壤質地分類標準配置12種土壤質地類型，在測定大量土壤水力特性參數的基礎上，采用HYDRUS-2D/3D模擬研究土壤容重和膜孔直徑對12種質地土壤入滲特性的影響；基于數值模擬結果，對濕潤體運移模型的參數進行優化，進而提出一種新的膜孔灌濕潤體尺寸估算模型；設計一種反映田間膜孔灌入滲過程的試驗裝置，用于驗證估算模型的可靠性，實現估算模型在田間的實際應用，為設計運行管理人員提供簡便的實用工具。

1 材料與方法

1.1 技術路線

首先，設計一種真實反映田間膜孔灌入滲過程的試驗裝置——圓筒入滲儀；其次，利用圓筒入滲儀田間試驗驗證HYDRUS-2D/3D模擬結果的準確性，在此基礎上，模擬研究不同影響因素下膜孔灌土壤入滲率和濕潤體尺寸變化規律；然后，基于180組數值模擬結果，優化確定膜孔灌土壤濕潤體運移模型參數；最后，采用田間圓筒入滲儀試驗和室內土箱試驗評價模型的可靠性，實現膜孔灌濕潤體運移模型在田間的實際應用。具體技術路線如圖1所示。

圖1 技術路線圖

1.2 膜孔灌田間試驗設計

根據田間膜孔灌特點，設計了一種能夠真實反映膜孔入滲過程的圓筒入滲儀，如圖2所示。該儀器類似于圭爾夫(Guelph)入滲儀，利用馬氏瓶原理提供恒定水頭，但圓筒入滲儀將供水單元與滲水單元一體化，試驗前無需組裝調試，操作方便；試驗前圭爾夫入滲儀需鉆孔打洞，這無法反映土壤表面膜孔入滲情景，而圓筒入滲儀無需鉆孔，僅豎直放于待測土壤之上即可，可真實反映地膜覆蓋下膜孔入滲情景；另外圓筒入滲儀采用有機玻璃粘接而成，制作簡單，成本低廉。

圖2 圓筒入滲儀田間試驗示意圖

圓筒入滲儀供水單元由儲水室和進氣管組成，進氣管一端通過進氣閥與大氣相通，另一端伸入儲水管直至底部，用于控制入滲室保持恒定水頭；滲水單元由入滲室、排氣管、出水管和透水墊組成，入滲室相當于膜孔裝置，排氣管用于試驗初期入滲室空氣的排出，透水墊可防止出水管水流對地面的擊濺。田間試驗具體操作步驟為：①選擇待測區域，確定待測地點，然后去除待測點地表的植被及石塊等大體積雜物，將地面整平。②關閉進氣閥和出水閥，通過橡膠塞孔往儲水室加滿水。③封閉橡膠塞，打開出水閥，直至儲水管水位穩定，將圓筒入滲儀豎直放入待測地點。④記錄儲水管初始水位，打開進氣閥，記錄不同時刻水位變化，待變化速率基本穩定時，停止試驗。通過以上4步即可確定膜孔穩滲率，進而推求土壤入滲參數，整個試驗過程中土壤無擾動。需要說明的是，文中所建濕潤體運移模型的驗證部分，需將水分傳感器埋入土中并回填土，為了減少土壤擾動產生的試驗誤差，待回填土與原土物理性質基本相同后再開展試驗，并增列步驟⑤，即：⑤在步驟④過程中，觀測不同時刻地表濕潤半徑，而垂向濕潤深度通過圓筒入滲儀下方水分傳感器獲得，記錄水分達到傳感器的時間。

1.3 膜孔灌數值模擬

1.3.1基本方程

假設土壤是均質且各向同性的，膜孔灌土壤水分運動可簡化為軸對稱二維入滲問題。采用HYDRUS-2D/3D軟件[17]進行模擬，其基本方程為Richards方程，即

(1)

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

t——時間，min

r——徑向坐標，cm

h——壓力水頭，cm

z——垂向坐標，規定向下為正，cm

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

通過VAN GENUCHTEN-MUALEM(VG-M)模型[18-19]來描述土壤水分特征曲線和土壤非飽和導水率，計算式為

(2)

(3)

式中Ks——土壤飽和導水率，cm/min

Se——土壤相對飽和度

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

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

α——與進氣值成反比的經驗參數，cm-1

β——影響土壤水分特征曲線形狀的經驗常數

1.3.2定解條件

充分考慮田間膜孔灌實際情況，設定不同建模情景的初始和邊界條件，如圖3所示，圖中S為膜孔間距，cm；θ0為土壤初始含水率，cm3/cm3；θc為田間持水率，cm3/cm3。

圖3 膜孔灌土壤水分運動計算域示意圖

為確保計算域邊界設置與田間實際相符，考慮膜孔的對稱性，以膜孔圓心為起點設置上邊界和左邊界；采用1/2膜孔間距作為右邊界，需要說明的是，對于不同的田間植物，其膜孔間距不等，差異較大，文中邊界設置主要是針對較大膜孔間距的自由入滲情景；對于下邊界可設置垂向距離足夠大(本文取100 cm)，以滿足模擬結束時刻不受灌水影響。在所有的模擬情景中，計算域內土壤水分均按初始含水率設置，根據文獻[20-23]的田間適宜灌水下限范圍(50%～70%田間持水率θc)，本文按60%田間持水率設置；上邊界AB為膜孔外圍覆膜地表，不受蒸發和降水影響，按零通量邊界設置；下邊界CE不受灌水影響，按自由排水邊界設置；左邊界EF為膜孔中心線，無水量交換，按零通量邊界設置；右邊界BC灌溉水分未達到，按零通量邊界設置；邊界AF為上方有水頭的膜孔區，考慮田間膜孔上方水頭(2～10 cm)難以控制且對土壤入滲影響較小，為簡化計算，本文按水頭6 cm定邊界設置。

綜上所述，初始條件可描述為

(4)

式中Ω——計算域(圖3)

最終邊界條件可總結為

(5)

式中H——膜孔上方水頭，cm

1.3.3模擬方案

參照國際制土壤質地分類標準，配置4個類別12種土壤[24]。將每種土壤設置3組容重(γ)，通過離心機法測定土壤水分特征曲線，RETC擬合確定θs、θr、α和β，定/變水頭法測定Ks，環刀法測定θc，相關參數見表1。在H=6 cm、θ0=60%θc情況下，設置180種情景(12種土質、3組容重、5個膜孔直徑)，模擬分析土壤容重和膜孔直徑對12種質地土壤入滲率和濕潤體尺寸的影響。

表1 不同容重下12種土壤質地VG-M模型參數和田間持水率

1.4 膜孔灌濕潤體運移過程描述

膜孔入滲過程中土壤濕潤體形狀近似于半橢球體，其特征軸為水平向濕潤半徑R和垂直向濕潤深度Z，兩者與時間具有良好的冪函數關系[5,7,13]，表達式為

R=atb+0.5D

(6)

Z=ctd

(7)

式中a、b、c、d——入滲參數

1.5 誤差分析

為了驗證估算模型可靠性，采用均方根誤差(RMSE)和納什效率系數(NSE)對估算值和實測值進行誤差分析。當NSE越趨近于1，RMSE越靠近于0，說明估算值越接近于實測值，模型精度越高。

2 結果與分析

2.1 HYDRUS-2D/3D數值模擬驗證

采用HYDRUS-2D/3D軟件模擬獲得粉壤土(D=8 cm、γ=1.33 g/cm3)不同時刻土壤入滲率和濕潤距離，將模擬結果與圓筒入滲儀田間試驗數據進行對比，如表2所示。

由表2可知，不同時刻土壤入滲率和濕潤距離模擬值與實測值相對誤差絕對值介于1.56%～9.09%之間，平均值為4.45%；入滲率的均方根誤差為0.196 mL/min，納什效率系數為0.992；濕潤距離(濕潤半徑、濕潤深度)的均方根誤差為0.181 cm，納什效率系數為0.990。說明文中所建數學模型是正確的，采用HYDRUS-2D/3D軟件模擬田間膜孔灌入滲過程可靠。

表2 土壤入滲率和濕潤距離模擬值與實測值對比

2.2 不同影響因素對入滲率的影響

在相同土壤初始含水率(θ0=60%θc)和灌溉水頭(H=6 cm)情況下，選定4種具有代表性的質地土壤(砂土、壤土、黏壤土、黏土)進行數值分析。圖4為不同土壤容重和膜孔直徑組合下4種質地土壤入滲率變化規律。

圖4 不同土壤容重和膜孔直徑組合下4種質地土壤入滲率變化曲線

由圖4可知，不同影響因素組合下土壤入滲率隨時間的變化規律相同，20組入滲率變化曲線都存在快速入滲和穩定入滲兩個階段，土壤質地越黏重，達到穩滲階段的用時越長，但基本在120 min后進入穩滲狀態。土壤質地對入滲率影響非常顯著，相同條件下，砂土的穩滲率明顯大于壤土、黏壤土和黏土，表現出土壤質地越黏重，土壤穩滲率越小的現象；土壤容重和膜孔直徑均對入滲率有較大影響，相同條件下，穩滲率隨膜孔直徑增大而增大，隨土壤容重增大而減小。這是因為膜孔直徑越大，土壤受水面積越大，單位時間內入滲水量越多，使穩滲率隨之增大；土壤容重越大，土壤滲透能力越弱，單位時間內入滲水量越少，導致穩滲率降低。

2.3 不同影響因素對濕潤體運移的影響

2.3.1土壤容重對濕潤體運移的影響

在θ0=60%θc、H=6 cm和D=6 cm情況下，對4種質地土壤(砂土、壤土、黏壤土和黏土)容重變化下的濕潤體特征值進行數值分析，其中砂土容重為1.55、1.60、1.65 g/cm3，其他土質取1.30、1.35、1.40 g/cm3。圖5為不同土壤容重條件下4種質地土壤的濕潤體運移曲線。

圖5 不同容重下4種質地土壤濕潤體運移曲線

由圖5可知，濕潤半徑和濕潤深度均隨時間的延長而逐漸增大，土壤質地越黏重，濕潤體運移越緩慢。如灌水結束(240 min)時，容重為1.35 g/cm3的壤土、黏壤土和黏土的濕潤半徑分別是砂土(γ=1.60 g/cm3)的0.49、0.39、0.26倍，而濕潤深度僅為0.35、0.27、0.17倍，說明土壤質地對濕潤體運移有顯著影響，在灌溉系統設計中土壤質地是必須要考慮的。對于同一土壤質地，濕潤半徑和深度均隨容重的增加而減小。以壤土為例，灌水結束時容重為1.35、1.40 g/cm3的濕潤半徑比1.30 g/cm3的分別減小7.3%和13.7%，而濕潤深度分別減小7.7%和14.6%，說明土壤容重變化會影響濕潤體運移過程，在灌溉系統優化設計中也是應該考慮的影響因素。

2.3.2膜孔直徑對濕潤體運移的影響

選取容重為1.60 g/cm3的砂土以及容重為1.35 g/cm3的壤土、黏壤土和黏土，在θ0=60%θc和H=6 cm情況下，對4種質地土壤膜孔直徑變化(D=4、6、8 cm)下的濕潤體特征值進行數值分析。圖6為不同膜孔直徑條件下4種質地土壤的濕潤體運移曲線。

圖6 不同膜孔直徑下4種質地土壤濕潤體運移曲線

由圖6可知，土壤質地差異對濕潤半徑和深度的影響規律與圖5相似，兩個方向上濕潤體運移距離由大到小表現為砂土、壤土、黏壤土、黏土。相同土壤質地和容重條件下，濕潤半徑和深度均隨膜孔直徑的增大而增大。以D=4 cm為基準，灌水結束時壤土的濕潤半徑和濕潤深度分別為15.48 cm和15.36 cm，當膜孔直徑增至6 cm和8 cm時，濕潤半徑分別增加14.5%(2.25 cm)和26.8%(4.15 cm)，濕潤深度分別增加12.8%(1.96 cm)和23.0%(3.53 cm)，說明膜孔直徑對濕潤體運移存在較大影響，是膜孔灌系統中非常重要的技術參數。

2.4 膜孔灌濕潤體運移模型建立

通過數值分析膜孔灌濕潤體運移影響因素可知：土壤質地、膜孔直徑、土壤容重和灌水時間均會對濕潤半徑和濕潤深度產生影響，是不容忽視的主導因素。因此，基于180組模擬結果，利用式(6)和式(7)擬合獲得入滲參數(a、b、c和d)。

對于180組膜孔灌土壤濕潤體運移過程的擬合，其決定系數R2均大于0.968，說明采用式(6)和式(7)描述濕潤半徑和濕潤深度與時間的關系合適。分析發現12種土壤質地入滲參數b和d的數值變化較小且無明顯規律，b介于0.28～0.42之間，d介于0.25～0.41之間，為簡化計算，分別取180組b和d的平均數，即b=0.35和d=0.32。將確定的b和d代入式(6)和式(7)，原式可進一步變換為

R=at0.35+0.5D

(8)

Z=ct0.32

(9)

利用式(8)和式(9)再次對HYDRUS-2D/3D模擬結果進行擬合，獲得a和c。發現R2≥0.912，說明將b和d取平均值來簡化模型可行。分析發現不同土壤質地下入滲參數a和c的數值差異很大且規律性不強，采用求平均值法來確定a和c已不現實。進一步分析發現，a和c與穩滲率f0存在良好的冪函數關系(R2≥0.972)，如圖7所示。

圖7 入滲參數a和c與穩滲率f0關系曲線

由圖7可知，a和c均隨f0的增大而增大，其冪函數關系可表示為

(10)

(11)

將式(10)和式(11)分別代入式(8)和式(9)，進一步得到

(12)

式(12)包含了f0這一待定參數，因此，如何確定f0是準確估算濕潤半徑和濕潤深度的關鍵。通過數值分析穩滲率影響因素可知：相同土壤質地條件下，穩滲率隨膜孔直徑增大而增大，而隨土壤容重增大而減小。選取4種代表性土壤質地，繪制穩滲率與土壤容重和膜孔直徑關系曲線，如圖8所示(圖中砂土f0左軸讀數，其他土壤f0右軸讀數，下同)。

圖8 4種土壤質地下f0與γ和D關系曲線

由圖8可知，穩滲率與土壤容重和膜孔直徑之間均呈良好的冪函數關系，據此，提出f0與γ和D的函數關系式為

f0=kγuDv

(13)

式中k——待定入滲參數

u、v——擬合參數

對12種土壤質地不同容重和膜孔直徑條件下濕潤體運移模擬結果進行多元非線性回歸，確定u=-6.3和v=1.1。k與土壤質地有關，數值差異較大，無顯著性規律，如圖9所示。

圖9 4種土壤質地f0與γ-6.3D1.1關系曲線

由圖9可知，對于任一土壤質地，u和v統一取值-6.3和1.1是可靠的(R2≥0.988)，而k是與土壤質地相關的待定系數。需要說明的是，土壤質地是根據顆粒機械組成劃分，這就導致了農田土壤質地類型保持不變而不同地塊差異較大的客觀事實，因此，在灌溉系統設計時土壤質地應分情況予以考慮?；诖?，將u=-6.3和v=1.1代入式(13)，并與式(12)聯立，最終可構建出膜孔灌土壤濕潤體運移模型

(14)

式(14)僅含一個待定入滲參數k，其值與土壤質地有關，僅需一組室內土箱試驗或田間圓筒入滲儀試驗獲得穩滲率f0，進而反推求得。

2.5 膜孔灌濕潤體運移模型評價

為了評價本文所建濕潤體運移模型的適用性和可靠性，分別采用文獻[5]和文獻[7]土箱試驗以及圓筒入滲儀田間試驗進行驗證。

文獻[5]采用榆林壤土(γ=1.40 g/cm3)，開展膜孔直徑為2～6 cm的單點膜孔入滲試驗?；贒=6 cm的累積入滲量數據，獲取f0=8.89 mL/min，推求k=10.34，進而得到榆林壤土膜孔灌濕潤體運移模型

(15)

采用文獻[5]中D=3 cm和D=5 cm濕潤體半徑和深度實測資料對式(15)進行驗證，如圖10a所示。

圖10 4種土壤膜孔灌土壤濕潤體尺寸計算值與實測值對比

文獻[7]通過西安粉砂土(γ=1.30 g/cm3)研究了膜孔直徑為4～8 cm的渾水膜孔灌入滲特性。本次驗證采用D=4 cm的入滲量數據推求f0和k，獲取西安粉砂土膜孔灌濕潤體運移模型

(16)

利用文獻[7]中D=6 cm和D=8 cm的濕潤半徑和深度實測資料對式(16)進行驗證，如圖10b所示。

文中模型建立的目的是為了能夠對田間膜孔灌土壤濕潤體尺寸進行快速的現場評估。因此，模型的可靠性和適用性也理應接受田間現場試驗資料的驗證。為此，分別選取蘭州市魏嶺鄉狗牙山村和臨夏回族自治州西河鎮紅城寺村基本農田，開展圓筒入滲儀田間試驗。為了體現土壤容重的影響，在相近地塊選取兩處天然容重存在差異的試驗點，分別用于參數f0和k的推求和濕潤體尺寸的驗證。

采用容重為1.33 g/cm3砂壤土的圓筒入滲儀(D=5 cm)水量變化計算f0，確定k=11.01。將k代入式(14)，得蘭州砂壤土膜孔灌濕潤體運移模型

(17)

利用容重為1.41 g/cm3砂壤土的圓筒入滲儀(D為4、6 cm)濕潤體運移實測資料對式(17)進行驗證，如圖10c所示。

采用容重為1.31 g/cm3粉壤土的圓筒入滲儀(D=5 cm)水量變化計算f0，確定k=8.25。將k代入式(14)，得臨夏粉壤土膜孔灌濕潤體運移模型

(18)

利用容重為1.38 g/cm3粉壤土的圓筒入滲儀(D為4、6 cm)濕潤體運移實測資料對式(18)進行驗證，如圖10d所示。

由圖10可知，模型計算值與試驗實測值在濕潤半徑和濕潤深度上均具有高度的一致性，采用T-Test檢驗，P>0.05，說明兩者之間均無顯著性差異。統計指標RMSE接近于0，介于0.020～0.170 cm之間，NSE不小于0.995，趨近于1，表明模型可靠性良好，適用于田間不同質地類型土壤濕潤體運移過程的預測。

3 結論

(1)設計了一種能夠模擬田間膜孔入滲過程的圓筒入滲儀，其供水單元與滲水單元一體化，結構簡單，操作方便，可用于測定田間原狀土膜孔灌溉模式下的穩滲率。

(2)相同土壤質地條件下，穩滲率與土壤容重和膜孔直徑呈良好的冪函數關系，冪函數指數分別為-6.3和1.1，冪函數系數只需1組圓筒入滲儀田間試驗即可推求。

(3)相同土壤質地條件下，濕潤半徑和濕潤深度與穩定入滲率和灌溉時間呈良好的冪函數關系，建立了一種適用于不同土壤質地類型的濕潤體運移模型。

(4)田間和室內試驗均證實模型可靠性良好和適用性較強，實現了膜孔灌濕潤體運移模型在田間的實際應用，方便管理人員進行快速的現場評估。