基于FlexPDE的膜孔灌土壤水分運動數值模擬

趙文剛，范嚴偉，劉曉群，石 林，邢旭光，馬孝義，宋 雯

(1.湖南省水利水電科學研究院，湖南 長沙 410007； 2. 蘭州理工大學能源與動力工程學院，甘肅 蘭州 730050；3.西北農林科技大學水利與建筑工程學院 旱區農業水土工程教育部重點實驗室，陜西 楊凌 712100)

土壤水分運動是陸地水循環的重要環節，也是土壤圈物質循環的重要組成部分，還是土壤中能量和溶質傳輸的主要載體[1]。掌握土壤水分運動規律對改進灌溉技術、合理利用水資源以及提高作物產量具有重要意義[2]。

目前，土壤水分運動規律研究主要依賴室內試驗、數值模擬。室內試驗因人工填土、夯實導致土壤各向同性受到影響；同時，對于多因素多水平試驗方案設置難以實現完全試驗，且完全試驗時無法保證因素絕對固定。數值模擬試驗依托Richards方程利用有限元法求解，已形成的軟件主要有HYDRUS和SWMS，其被廣泛應用于灌水方式[3-5]、水質狀況[6]、土壤分層情況[7-8]等對土壤水分運動的影響研究中，但缺乏實測資料支撐。基于二者優缺點，通常采用室內試驗、數值模擬相結合的方法。

SWMS軟件模擬土壤水分運動精度較高，但其計算前期需手動剖分網格、計算單元節點控制長度等工作，工作量大且易出錯；此外，其基于Fortran 語言開發導致用戶界面不友好。HYDRUS軟件雖已實現網格自動劃分，但網格劃分依賴經驗且網格密度單一，對于復雜計算區域，無法局部加密處理，影響計算精度。此外，為簡化計算，SWMS和HYDRUS等軟件模擬土壤水分運動規律時采用飽和導水率代替非飽和導水率，而研究發現非飽和導水率受孔隙狀況、溫度、含水率影響[9-10]。

FlexPDE是一款專門求解偏微分方程的軟件，具有網格自動劃分、網格數量和密度隨誤差自動調整等優勢，能通過編輯問題描述腳本文件將偏微分方程描述的系統轉化為有限元模型求解，用戶界面友好。Bader[11]利用FlexPDE對粘土薄膜的水分和溶質運移特性進行了模擬，且效果理想；Albrieu[12]采用FlexPDE模擬了定體積土壤中根系生長的吸水模型；Rühaak[13]通過FlexPDE模擬冰川效應下的三維水力機械耦合地下水流動；唐延貴[14]通過壓力板儀和非飽和土三軸儀開展試驗研究，結合FlexPDE分析了非飽和土-水特征和變形特征；此外，FlexPDE在溫度場、電場和熱傳導[15-17]等方面也有所應用。

相比常規地面灌溉而言，膜孔溝灌在重力影響下的垂向運移與水勢梯度作用下的側向運移共同作用對土壤進行局部濕潤的灌溉，使得灌水更加均勻、高效，也更節水。綜合上述，本文采用土壤擴散率與土壤比水容量的函數關系推求土壤非飽和導水率；建立土壤儲水量模型，分析土壤儲水量與累積入滲量的關系，驗證FlexPDE軟件在模擬膜孔溝灌土壤水分運動規律可靠性的同時也為合理確定灌溉技術要素提供參考依據，以便更好地指導農業生產。

1 材料與方法

1.1 試驗設計

供試土壤取自楊凌三道塬、渭河灘0～30 cm土層的粉粘壤土和砂壤土(顆粒組成見表1)。初始含水率取60%田間持水率[18]，膜孔間距一般為12～30 cm、膜孔直徑3～8 cm、灌溉水深2～10 cm[19-20]。據此，研究共設置4個處理，具體布設方案列于表2，兩種土壤的水力參數由RETC軟件擬合獲取(表3)。其中處理T1、T2采用A土箱進行試驗，處理T3、T4采用B土箱進行試驗。

表1 供試土壤顆粒組成

注：土壤質地分類依據國際制進行分類，下同。

Note: The soil texture was classified by international system, the same below.

表2 土箱試驗方案

表3 供試土壤水分特征參數擬合

注：ρ、θr、θs、α、n分別代表土壤容重、殘余體積含水率、飽和體積含水率、土壤進氣值的倒數、參數，下同。

Note:ρ、θr、θs、α、andnrepresent the bulk density, residual volumetric water rate, saturated volumetric water rate, the reverse of air entry value, and parameters, the same below.

1.2 試驗方法

膜孔灌溉裝置由馬氏瓶、有機玻璃土箱和膜孔裝置組成(圖1)：馬氏瓶用于維持恒定水頭；土箱規格(長×寬×高)為10 cm × 10 cm × 60 cm(A)和15 cm × 15 cm × 60 cm(B)，土箱A水平向從距邊壁2.5 cm間隔2.5 cm、垂向從距上邊壁5 cm間隔5 cm設置取土孔，土箱B水平向從距邊壁3 cm間隔3 cm、垂向從距上邊壁5 cm間隔5 cm設置取土孔，在土箱相鄰兩面對應位置均設置取土孔(可視為重復)；對于膜孔裝置，將1/4膜孔(即土箱一角)作為研究對象。土樣經風干、過2 mm篩，按設定容重每5 cm厚度裝填，土壤裝填完成后，除膜孔處均采用薄膜覆蓋以防止水分散失；試驗過程中監測馬氏瓶水位獲得入滲量；實驗結束時，從取土孔取土，采用干燥法測定不同深度土壤含水率及土壤剖面水分分布特征。

土壤擴散率采用水平土柱法測定，有機玻璃土柱材質，直徑(內徑)5 cm、長36 cm。土樣經風干、過2 mm篩，按設定容重將含水率為風干含水率土壤均勻裝入水平土柱中，試驗過程中忽略重力作用，馬氏瓶供水，記錄馬氏瓶水位變化及濕潤鋒每過1 cm所用時間；試驗結束時，從濕潤鋒處開始取土，采用干燥法測定土壤含水率。

a.膜孔裝置；b.取土孔；c.通氣孔；d.土面；e.馬氏瓶 a. film-hole settings; b. soil collecting hole; c. air vents; d. soil surface; e. Markova bottle圖1 膜孔灌實驗裝置Fig.1 The experimental equipment of film hole

2 模型建立

2.1 土壤水分運動模型

膜孔灌溉條件的Richards方程如下[20]：

(1)

式中，h為土壤基質勢(cm)；k(h)為非飽和土壤導水率(cm·min-1)；t為時間(min)；x為水平距離(cm)；z為垂直距離(cm)。

土壤基質勢和擴散率與含水率的關系可分別采用van-Genuchten方程[1]和常用的經驗公式[21]來擬合，即式(2)、(3)；土壤非飽和導水率與土壤基質勢的關系可采用式(4)擬合。

(2)

式中，θr為殘余體積含水率(cm3·cm-3)；θs為飽和體積含水率(cm3·cm-3)；n、α、m為經驗常數，其中m=1-1/n。

D(h)=D0×exp(-β×(θ0-θ(h)))

(3)

k(h)=C(h)×D(h)

(4)

式中，D0、β、θ為經驗常數，θ0為某一特征含水率，本文中取θ0=0；C(h)為土壤比水容量(cm-1)；D(h) 為土壤擴散率(cm2·min-1)。

本研究中，利用FlexPDE軟件對土壤水分運動特性進行模擬，且模型可通過式(1)、(4)獲得二維條件下的Richards方程：

(5)

2.2 模型初始條件和邊界條件

將膜孔直徑所在平面作為研究對象(圖2)，由于關于z軸對稱，故可取OABC圍成的區域作為研究對象，s為膜孔半徑，L為土箱的長度，H為土箱的高度，則初始條件和邊界條件可表述為(T為入滲時間)：

初始條件：

h(x,z,0)=h0(x,z) 0≤x≤L,0≤z≤H

上邊界條件：

下邊界條件：

左邊界條件：

右邊界條件：

圖2 膜孔灌土壤水分運動簡化模型示意圖Fig.2 The diagram of simplified model for water flow of film-hole irrigation

2.3 土壤儲水量模型

取OABC面作為研究對象(圖2)，t和t+Δt時刻的土壤儲水量按式(6)、(7)計算，W(t)為t時刻土壤剖面整體含水量，則Δt時間間隔內土壤儲水量ΔW可按式(8)計算：

W(t)=?θ(x,z,t)dxdy

(6)

W(t+Δt)=?θ(x,z,t+Δt)dxdy

(7)

ΔW=?θ(x,z,t+Δt)dxdy-?θ(x,z,t)dxdy

(8)

2.4 FlexPDE軟件簡介

FlexPDE軟件求解偏微分方程問題的可讀腳本文件主要由下面幾部分組成[22]：TITLE——輸入描述標簽項；SELECT——控制FlexPDE參數；VARIABLES——偏微分方程變量；DEFINITIONS——定義偏微分方程中出現的參數；EQUATIONS——偏微分方程；BOUNDARIES——邊界條件；MONITORS and PLOTS——數值解的圖形；END——腳本結束。

3 結果與分析

3.1 土壤水分動態模擬模型有效性驗證

土壤剖面含水率分布特征可以反映膜孔灌溉的入滲規律[2]，為了驗證膜孔灌溉條件下土壤水分運動模型的有效性，對土壤剖面含水率的實測值與模擬值進行對比分析。

圖3為通過FlexPDE軟件模擬膜孔溝灌灌水結束時刻的土壤剖面水分分布狀況。從圖中可以看出：靠近膜孔位置，土壤水分以膜孔為中心呈橢圓形分布；離膜孔位置越遠，膜孔對土壤水分分布特征影響越小，呈水平形態分布。

4種處理按照設計的灌水定額所需的灌水時間分別為120、642、209 min和568 min，灌水結束后，通過對比土壤剖面含水率實測與模擬值可以發現，除T3中土壤水分模擬值小于實測值之外，其余處理中土壤水分模擬值均高于實測值，處理T1～T4土壤含水率的平均相對誤差分別為9.600%、5.930%、8.340%和14.040%。進一步分析發現，粉粘壤土FlexPDE的模擬效果較砂壤土理想，即對于粉粘壤土，土壤含水率相對誤差小于5%、10%和20%的比例分別為40%、68%和84%；砂壤土，其相對誤差小于5%、10%和20%的比例分別為12%、40%和72%。而范嚴偉[23]基于SWMS對此研究得出，各處理平均相對誤差分別為10.230%、5.830%、7.660%和8.540%；其中對于粉粘壤土，土壤含水率相對誤差小于5%、10%和20%的比例分別為48%、72%和96%，與本研究結果對比發現，FlexPDE的模擬誤差較SWMS略高，但二者差異很小，亦滿足精度要求，同時鑒于FlexPDE的優點(前已述及)，因此該軟件更適合對粉粘壤土的水分入滲進行模擬。

土壤含水率模擬與實測值之間存在一定差異，分析其原因認為：(1)土壤水分特征曲線測定依賴有限點的土壤基質勢與含水率擬合，且用土壤脫濕過程來代替土壤吸濕過程，存在一定滯后性；(2)非飽和土壤水分運動參數的測定導致誤差；(3)缺少實時監測設備，灌水結束時土壤含水率測定有所延遲；(4)假設土壤各向同性，但實際土壤各向異性；(5)人工裝土與設計容重存在差異。

3.2 土壤儲水量與累積入滲量的關系

由圖4可以看出，各處理累積入滲量、儲水量與入滲時間均呈冪函數關系，累積入滲量與土壤儲水量的關系可以用線性模型描述，且所有擬合模型均具有較高精度(R2>0.900)。其中，累積入滲量隨時間延長而增大，但入滲速率隨時間延長而減小；儲水量隨累積入滲量增大以恒定速率增大。基于儲水量與累積入滲量的線性關系，可通過儲水量的變化反映累積入滲量的變化情況。

綜合圖3、4可以看出，基于實測土壤非飽和擴散率和土壤水分特征曲線，FlexPDE對粉粘壤土水分運動模擬效果優于砂壤土，但均滿足精度要求。這可能因為砂壤土大孔隙較多，飽和含水率相對較低，相同入滲量使得擴散率變化量增大，土壤濕潤鋒運移速度加快，含水率分布變化加快，進而導致模擬比較困難。

4 模型應用

膜孔灌，就單個膜孔而言，屬于局部不交匯灌溉；但在實際田間管理中，隨灌溉水量、歷時增加，灌溉水發生交匯入滲，從而導致土壤水分再分布。水分運動及再分布特性受入滲影響顯著，而膜孔灌溉水分入滲特性受膜孔直徑和間距等因素影響，因此，為提高田間灌溉水利用率，有必要深入研究膜孔直徑、間距和灌溉水深等因素對土壤水分運動特性的影響，從而實現田間水的有效管理。

圖3 不同處理土壤剖面水分分布Fig.3 The moisture distribution in soil profile under different treatments

圖4 儲水量和累積入滲量與入滲時間的關系Fig.4 The relationship among soil water storage, accumulative infiltration, and infiltration time

4.1 模擬方案

試驗土壤為楊凌一級階地的粉粘壤土(表4)，利用RETC擬合的水力參數性質見表5。膜孔交匯灌溉水分模擬，重點研究膜孔間距、入滲水深和膜孔直徑3個因素對入滲的影響，初始含水率取60%田間持水率，入滲時間設置為300 min，具體模擬試驗方案見表6。

4.2 不同處理交匯時間、水分分布和儲水量分析

4.2.1 膜孔灌單向交匯時間 由表7可知，對于膜孔間距、直徑相同但水深不同的處理，交匯時刻基本相同，這表明入滲水深對膜孔灌單向交匯時間的影響不顯著，這可能因為入膜水深雖增大，但只增加了初滲時刻土壤表面膜孔處的水勢梯度，起到瞬時加速濕潤土壤的作用，一旦膜孔處土壤水分達到飽和時，土壤水勢梯度就變為水分飽和土壤與土壤初始土水勢之間的差異，入滲水深失去作用。此外，膜孔直徑、水深相同條件下，膜孔間距越小，交匯所需時間縮短，這主要因為間距減小，縮短了水分水平運動距離；膜孔間距、入滲水深相同條件下，膜孔直徑越大，達到交匯時間越短，這主要由于膜孔孔徑增大一方面增大了初始時刻的濕潤面，另一方面縮短了水分水平運動距離。

表4 模擬方案供試土壤顆粒組成

表5 土壤水力特性參數擬合

4.2.2 土壤儲水量變化 從圖5可以看出，儲水量隨時間延長而增大；相同膜孔直徑、膜間距條件下，入膜水深對土壤儲水量影響微弱；相同膜孔直徑、入膜水深條件下，膜孔間距對于土壤儲水量的影響也較微弱。分析其原因認為，儲水量變化主要與土壤含水率變化相關，含水率變化與進入土體水量相關，膜孔直徑相同的情況下，相同時間周期內水深或膜間距變化使得進入土體水量差異不大，進而導致土壤儲水量變化不大。相同膜間距、入膜水深條件下，膜孔直徑對儲水量影響較顯著，且隨膜孔直徑增大而增大，這主要因為隨膜孔直徑增大，水分進入土壤通道增大，使得入滲水量增加，進而導致儲水量增大。通過對土壤儲水量和入滲時間擬合發現，各處理均呈現出良好的冪函數關系(表8)。

4.2.3 土壤水分分布狀況 圖6顯示，在土壤入滲過程中當濕潤鋒未發生交匯時為自由入滲，水分分布呈橢圓形；當濕潤鋒發生交匯時，隨時間推移土壤水分運動逐漸趨向水平向下推進。對比M1和M2、M3和M4可以看出，膜孔間距較大時，濕潤鋒發生交匯時間延后，土壤水分分布不均勻；而膜孔間距較小時，土壤水分分布均勻程度高。通過對其它處理對比分析也均反映出與上述類似的特性。由此可知，膜孔間距對土壤水分分布影響顯著：膜孔間距越小，土壤水分分布不均勻程度越小；而其他條件相同時，入滲水深越大，相同入滲時間，垂向濕潤鋒會有所增大，但增大幅度不明顯；膜孔間距相同條件下，灌溉水深、膜孔直徑對水分入滲的影響比較微弱，這與范嚴偉[3]、費良軍[6]等的研究結果一致。

表6 膜孔交匯灌溉模擬方案

表7 各處理單向交匯時刻

圖5 儲水量與入滲時間的關系Fig.5 The relationship between soil water storage and infiltration time

處理Treatment擬合公式Fitting functionR2處理Treatment擬合公式Fitting functionR2M1W=0.381t0.7170.999M10W=0.468t0.6940.999M2W=0.339t0.7310.992M11W=0.467t0.6950.999M3W=0.354t0.7300.999M12W=0.457t0.6980.999M4W=0.400t0.7030.999M13W=0.456t0.7180.993M5W= 0.349t0.7310.998M14W=0.453t0.7180.993M6W=0.397t0.7030.999M15W=0.543t0.6860.999M7W=0.461t0.7000.999M16W=0.460t0.7140.992M8W=0.459t0.6980.999M17W=0.541t0.6850.999M9W=0.468t0.6960.999M18W=0.531t0.6880.998

注：鑒于篇幅所限，只列舉M1～6，8。 Note: Due to space limitation, only treatments 1 to 6 and 8 are listed.圖6 膜孔單向交匯水分分布Fig.6 The moisture distribution of film-hole single-line convergence infiltration

5 結 論

基于對60 cm土層深度的土壤水分分布監測，同時測定了土壤擴散率，并結合FlexPDE軟件對膜孔灌土壤水分以及土壤儲水量進行模擬，主要得到以下結論：

1)通過Richards方程對土壤水分運動進行模擬可知FlexPDE軟件適合于粉粘壤土的土壤水分入滲模擬，滿足精度要求；

2)基于建立的土壤儲水量模型，儲水量與累積入滲量呈現較好的線性關系，且二者均隨入滲時間的延長而增加，但增長速率逐漸減小；

3)基于粉粘壤土建立的土壤水分、儲水量模型，對膜孔單向交匯入滲進行預測，結果表明，膜孔間距對于交匯初期土壤的水分分布影響較大，膜孔直徑對于土壤交匯時間以及儲水量影響顯著。

本文利用FlexPDE對于膜孔灌溉土壤水分狀況進行初步模擬，對于實際膜孔灌溉技術要素的確定具有指導意義。在本研究中，只是對于土壤初始含水率較高的情況進行了模擬，對于土壤初始含水率較低的情況需要進一步研究。