基于FlexPDE的膜孔灌土壤水分運(yùn)動(dòng)數(shù)值模擬

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

(1.湖南省水利水電科學(xué)研究院，湖南 長沙 410007； 2. 蘭州理工大學(xué)能源與動(dòng)力工程學(xué)院，甘肅 蘭州 730050；3.西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院 旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室，陜西 楊凌 712100)

土壤水分運(yùn)動(dòng)是陸地水循環(huán)的重要環(huán)節(jié)，也是土壤圈物質(zhì)循環(huán)的重要組成部分，還是土壤中能量和溶質(zhì)傳輸?shù)闹饕d體[1]。掌握土壤水分運(yùn)動(dòng)規(guī)律對(duì)改進(jìn)灌溉技術(shù)、合理利用水資源以及提高作物產(chǎn)量具有重要意義[2]。

目前，土壤水分運(yùn)動(dòng)規(guī)律研究主要依賴室內(nèi)試驗(yàn)、數(shù)值模擬。室內(nèi)試驗(yàn)因人工填土、夯實(shí)導(dǎo)致土壤各向同性受到影響；同時(shí)，對(duì)于多因素多水平試驗(yàn)方案設(shè)置難以實(shí)現(xiàn)完全試驗(yàn)，且完全試驗(yàn)時(shí)無法保證因素絕對(duì)固定。數(shù)值模擬試驗(yàn)依托Richards方程利用有限元法求解，已形成的軟件主要有HYDRUS和SWMS，其被廣泛應(yīng)用于灌水方式[3-5]、水質(zhì)狀況[6]、土壤分層情況[7-8]等對(duì)土壤水分運(yùn)動(dòng)的影響研究中，但缺乏實(shí)測資料支撐?；诙邇?yōu)缺點(diǎn)，通常采用室內(nèi)試驗(yàn)、數(shù)值模擬相結(jié)合的方法。

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

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

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

1 材料與方法

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

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

表1 供試土壤顆粒組成

注：土壤質(zhì)地分類依據(jù)國際制進(jìn)行分類，下同。

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

表2 土箱試驗(yàn)方案

表3 供試土壤水分特征參數(shù)擬合

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

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 試驗(yàn)方法

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

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

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

2 模型建立

2.1 土壤水分運(yùn)動(dòng)模型

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

(1)

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

土壤基質(zhì)勢(shì)和擴(kuò)散率與含水率的關(guān)系可分別采用van-Genuchten方程[1]和常用的經(jīng)驗(yàn)公式[21]來擬合，即式(2)、(3)；土壤非飽和導(dǎo)水率與土壤基質(zhì)勢(shì)的關(guān)系可采用式(4)擬合。

(2)

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

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

(3)

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

(4)

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

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

(5)

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

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

初始條件：

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

上邊界條件：

下邊界條件：

左邊界條件：

右邊界條件：

圖2 膜孔灌土壤水分運(yùn)動(dòng)簡化模型示意圖Fig.2 The diagram of simplified model for water flow of film-hole irrigation

2.3 土壤儲(chǔ)水量模型

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

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——輸入描述標(biāo)簽項(xiàng)；SELECT——控制FlexPDE參數(shù)；VARIABLES——偏微分方程變量；DEFINITIONS——定義偏微分方程中出現(xiàn)的參數(shù)；EQUATIONS——偏微分方程；BOUNDARIES——邊界條件；MONITORS and PLOTS——數(shù)值解的圖形；END——腳本結(jié)束。

3 結(jié)果與分析

3.1 土壤水分動(dòng)態(tài)模擬模型有效性驗(yàn)證

土壤剖面含水率分布特征可以反映膜孔灌溉的入滲規(guī)律[2]，為了驗(yàn)證膜孔灌溉條件下土壤水分運(yùn)動(dòng)模型的有效性，對(duì)土壤剖面含水率的實(shí)測值與模擬值進(jìn)行對(duì)比分析。

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

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

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

3.2 土壤儲(chǔ)水量與累積入滲量的關(guān)系

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

綜合圖3、4可以看出，基于實(shí)測土壤非飽和擴(kuò)散率和土壤水分特征曲線，F(xiàn)lexPDE對(duì)粉粘壤土水分運(yùn)動(dòng)模擬效果優(yōu)于砂壤土，但均滿足精度要求。這可能因?yàn)樯叭劳链罂紫遁^多，飽和含水率相對(duì)較低，相同入滲量使得擴(kuò)散率變化量增大，土壤濕潤鋒運(yùn)移速度加快，含水率分布變化加快，進(jìn)而導(dǎo)致模擬比較困難。

4 模型應(yīng)用

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

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

圖4 儲(chǔ)水量和累積入滲量與入滲時(shí)間的關(guān)系Fig.4 The relationship among soil water storage, accumulative infiltration, and infiltration time

4.1 模擬方案

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

4.2 不同處理交匯時(shí)間、水分分布和儲(chǔ)水量分析

4.2.1 膜孔灌單向交匯時(shí)間 由表7可知，對(duì)于膜孔間距、直徑相同但水深不同的處理，交匯時(shí)刻基本相同，這表明入滲水深對(duì)膜孔灌單向交匯時(shí)間的影響不顯著，這可能因?yàn)槿肽に铍m增大，但只增加了初滲時(shí)刻土壤表面膜孔處的水勢(shì)梯度，起到瞬時(shí)加速濕潤土壤的作用，一旦膜孔處土壤水分達(dá)到飽和時(shí)，土壤水勢(shì)梯度就變?yōu)樗诛柡屯寥琅c土壤初始土水勢(shì)之間的差異，入滲水深失去作用。此外，膜孔直徑、水深相同條件下，膜孔間距越小，交匯所需時(shí)間縮短，這主要因?yàn)殚g距減小，縮短了水分水平運(yùn)動(dòng)距離；膜孔間距、入滲水深相同條件下，膜孔直徑越大，達(dá)到交匯時(shí)間越短，這主要由于膜孔孔徑增大一方面增大了初始時(shí)刻的濕潤面，另一方面縮短了水分水平運(yùn)動(dòng)距離。

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

表5 土壤水力特性參數(shù)擬合

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

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

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

表7 各處理單向交匯時(shí)刻

圖5 儲(chǔ)水量與入滲時(shí)間的關(guān)系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 結(jié) 論

基于對(duì)60 cm土層深度的土壤水分分布監(jiān)測，同時(shí)測定了土壤擴(kuò)散率，并結(jié)合FlexPDE軟件對(duì)膜孔灌土壤水分以及土壤儲(chǔ)水量進(jìn)行模擬，主要得到以下結(jié)論：

1)通過Richards方程對(duì)土壤水分運(yùn)動(dòng)進(jìn)行模擬可知FlexPDE軟件適合于粉粘壤土的土壤水分入滲模擬，滿足精度要求；

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

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

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