滴灌雙點源土壤水分入滲特性數(shù)值模擬研究

李耀剛，蔣存仁

（1.甘肅天水市水務局，甘肅天水741000；2.甘肅省水利水電勘測設計研究院，蘭州730000）

滴灌作為一種現(xiàn)代高效節(jié)水灌溉技術在農(nóng)田節(jié)水灌溉中被廣泛使用，它可以根據(jù)田間不同作物需水要求進行灌溉。滴灌屬于局部微灌溉，通過頻繁、緩慢地將滴頭處消能后的少量水浸潤作物根系層土壤，在時間與空間上能夠比較精準地調控土壤含水量，同時，土壤地表水分無效蒸發(fā)和深層滲漏損失大幅減小，水分利用效率高達90%以上，相較其他農(nóng)田灌溉方式節(jié)水效益明顯[1]。作物根系經(jīng)由滴灌入滲形成的濕潤體獲取水分，根系分布狀況及生長狀態(tài)一定程度上受土壤濕潤體水分分布的影響。滴灌條件下土壤水分運動規(guī)律的研究可為滴灌系統(tǒng)設計中合理確定灌水持續(xù)時間、灌水流量及滴頭布置間距提供重要的理論依據(jù)。點源滴灌入滲土壤水分分布規(guī)律及濕潤體特性被國內外學者進行了大量試驗研究和數(shù)模模擬，得出了一些比較成熟的理論和方法[2-4]。張志剛[5]等研究了不同滴頭流量對地上滴灌土壤水分分布的影響，得出土壤濕潤鋒的入滲距離比值、運移速率與灌水時間呈冪函數(shù)關系。張松[6]等通過室內試驗研究了地下點源滴灌入滲土壤水分入滲規(guī)律，初步研究了地下點源滴灌在沙土中的適應性。數(shù)值計算方面以Richards方程為基礎的點源入滲數(shù)值模型被提出，裴青寶[7]等對多點源入滲條件下土壤水分溶質運移進行了模擬；李斯[8]運用Hydrus模型模擬了土壤有效水邊界對滴灌入滲土壤濕潤體特性的影響；許迪等[9]通過建立的非飽和土壤水分運動和溶質運移數(shù)值模型模擬了地下單點源滴灌入滲中水、肥分布規(guī)律；Skaggs[10]等對地上、地下滴灌土壤水分運動過程應用Hydrus-2d軟件進行了模擬。前人對點源滴灌入滲試驗及數(shù)值模擬進行了許多有益的探索，然而在農(nóng)田灌溉中，對于株間距較小的農(nóng)作物如棉花、玉米等，為使每株植物都能從土壤中獲取水分，滴灌入滲就需要在濕潤區(qū)形成濕潤帶，而多點源入滲在滴頭連線方向土壤濕潤體呈帶狀分布，因此多點源滴灌入滲濕潤體交匯情況下土壤水分運動的研究具有重要的現(xiàn)實意義。本文在上述研究的基礎上，以雙點源滴灌條件下的土壤水分分布為研究對象，根據(jù)非飽和土壤水動力學理論，利用軟件HYDRUS-3D建立了物理和數(shù)值模型，利用模型計算輸出結果分析了雙點源滴灌滴頭流量和灌水時間對土壤水分運動特征的影響，以期為進一步研究多點源滴灌條件下土壤水分運移規(guī)律提供理論依據(jù)。

1 數(shù)學模型

1.1 基本方程

忽略氣象因素、容質勢、溫度勢對土壤水分運動的影響，假定土壤各項同性且為剛性介質，土壤中的水分不可壓縮，不考慮土壤水分滯后效應。三維等溫均勻達西水流水分運移方程用修改過的Richards方程：

式中：ψ 為土壤總水勢， cm；θ 為土壤體積含水率，cm3/cm3； k(θ)為土壤非飽和導水率， cm/min；t 為計算時間，min，x、y、z分別為橫縱軸向坐標，cm。

1.2 定解條件

1.2.1 初始條件

滴灌模擬區(qū)域幾何邊界如圖1 所示。尺寸為：100 cm×100 cm×60 cm（長×寬×高)，模擬土壤介質選擇沙壤土，容重為1.35 g/cm3，假定土壤初始剖面含水率（負壓水頭）均勻分布，為方便計算，本文將土壤含水率通過土壤水分特征曲線轉換為負壓水頭，其值為0.085 cm3/cm3(負壓水頭為-1 200 cm水柱）。求解方程的初始條件為：

式中：h0為土壤初始負壓水頭，cm。X、Y、Z分別為模擬計算區(qū)域的橫向、縱向及垂向最大距離。本文中，X=100 cm，Y=100 cm，Z=60 cm。

1.2.2 邊界條件

上邊界：

（1）滴頭處：

式中：q為滴頭流量，L/h；E為土壤表面蒸發(fā)強度，cm/min。

（2）滴頭之外：

（3）地表積水，飽和區(qū)h=0；非飽和區(qū)：

下邊界：h=h0，z=Z，0≤t≤T。

左右邊界：假定邊界對稱，考慮邊界上水分通量為0：

1.3 模型求解

采用隱式差分法離散求解方程式（1）：

式中：i、j 分別為距離、時間節(jié)點編號；j+1、j 為當前和前一計算時刻；k為迭代次數(shù)。

應用HYDRUS-3D對建立的數(shù)值模型進行求解運算。計算區(qū)域采用有限元網(wǎng)格進行剖分，通過GALERKIN 有限元法求解。滴頭處水勢梯度大，通過增加計算節(jié)點加密網(wǎng)格，滿足計算精度要求。初始、最小和最大時間步長分別按0.01、0.001、10 min 計；運算結果每5 min間隔輸出一次；最大迭代次數(shù)為10 次；壓力水頭、含水率容差分別為1、0.001。初始條件采用土壤初始負壓水頭，假定計算區(qū)域土壤初始含水率均勻分布，計算區(qū)內負壓水頭均為-1 200 cm。邊界條件采用連續(xù)流輸入，計算區(qū)下邊界設為自由排水面。為分析濕潤體內土壤水分變化情況，在計算區(qū)距離滴頭水平、豎直間隔5 cm處設置20個觀測點。

2 試驗驗證

2.1 土壤水動力學特性

采用沙壤土作為供試土樣，土壤容重為1.35 g/cm3，土壤平均含水率為0.08 cm3/cm3。供試土樣物理性狀：粒徑大于0.02 mm 的沙粒含量占比55.29%，粒徑介于0.02～0.002 mm 的粉粒含量占比33.87%，粒徑小于0.002 mm 的黏粒含量占比10.84%。土壤水分特征曲線通過壓力膜儀測定擬合，土壤飽和導水率通過定水頭法測定。土壤導水率和土壤水分特征曲線通過van Genuchten模型進行擬合：

式中：5 個參數(shù)θs，θr，α，n，Ks獨立；θr為土壤殘余含水率，cm3/cm3；θs為土壤飽和含水率，cm3/cm3；α、m、n 為擬合參數(shù)，其中α 的大小與土壤物理性質密切相關，cm-1；L 為孔隙連通性參數(shù)；Se為土壤有效含水量，cm3/cm3；Ks為土壤滲透系數(shù)，cm/min。

本文中，土壤水分特性參數(shù)擬合值分別為：θs=0.45 cm3/cm3，θr=0.018 cm3/cm3，α=0.009 cm-1，n=1.6，Ks=0.006 cm/min。

2.2 結果驗證

圖2 給出了沙壤土在滴頭流量為1.2 L/h，灌水時間為300 min時距滴頭水平和垂直方向不同位置處土壤含水率的實測與模擬結果對比。從含水率柱狀圖可以看出，實測與模擬結果比較接近，水平距離土壤含水率相對誤差分別為5.74%、3.38%、4.42%、4.21%；垂直距離土壤含水率相對誤差分別為0.68%、3.36%、5.09%、8.45%、4.94%，均在10%以內。因此，濕潤體內土壤水分分布的數(shù)值模擬能較好地反映土壤水分分布的實際情況。

3 雙點源滴灌濕潤體內土壤水分分布模擬結果與分析

利用建立的雙點源入滲模型模擬了灌水時間、滴頭流量對滴灌土壤水分分布規(guī)律的影響，對濕潤體內土壤水分分布規(guī)律進行了初步探索研究。

3.1 灌水時間對濕潤體內土壤水分分布的影響

土壤水分入滲濕潤體大小的變化規(guī)律，是正確進行雙點源滴灌系統(tǒng)田間布置和灌溉水分管理的重要依據(jù)。沿滴頭布置方向的橫斷面可以反映雙點源灌水條件下的土壤水分交匯過程，因此本文主要對該斷面上的土壤水分分布規(guī)律進行分析。圖3反映了滴頭流量為1.2 L/s，灌水持續(xù)時間分別為250、300、350、400、500、600 min 時，沙壤土土壤水分分布的模擬過程。坐標（0、0）、（30、0）為滴頭所在位置。從圖3 可以看出，入滲開始后，濕潤鋒未交匯前，土壤水分入滲遵循點源入滲，濕潤體內土壤含水率等值線剖面形狀近似呈橢圓形。土壤水分隨時間推移由滴頭向外不斷擴散，土壤含水率等值線圖離滴頭由近及遠的位置上呈由疏到密分布。灌水時間為250 min 時，點源入滲濕潤鋒水平距離小于15 cm，垂直入滲深度大于15 cm，濕潤鋒沒有交匯，滴頭處土壤含水量達到土壤飽和含水量。當灌水時間超過268 min 后，濕潤體開始融合交匯。隨著入滲時間推移，濕潤體形狀逐漸由2個分離的半橢圓變成花生殼形狀，除靠近滴頭處土壤含水率較高的區(qū)域土壤含水率等值線圖近似于1/4 橢圓形外，遠離滴頭的濕潤區(qū)內土壤等值線圖融合為一條曲線，并且隨著時間的推移，曲線向下移動，等值線形狀從屋脊形過渡為水平。交匯區(qū)內下部土壤含水率等值線分布平緩，上部分布彎曲。灌水時間增加，滴頭處土壤水分飽和區(qū)范圍逐漸增大，交匯面土壤含水率增加。

分析雙點源滴灌入滲條件下濕潤體內的土壤含水率變化對于科學合理地施行作物水分管理，制定合理的滴灌灌溉制度有一定的指導意義。圖4 所示滴頭流量為1.2 L/h 灌水600 min 情況下的距離滴頭不同位置觀測點（測點布置見圖5）的模擬土壤含水率變化情況，由于滴頭對稱布置，濕潤體對稱分布，故選取一側的觀測點模擬數(shù)據(jù)進行分析。由圖4 可和，灌水開始后，滴頭處（N1）土壤含水率迅速升高，10 min 內由0.08 cm3/cm3增加到0.42 cm3/cm3，隨后增長速率趨于平緩。由測點N4、N3、N2可知，當時間分別為25、110、268 min時水分運移至距滴頭水平方向5、10、15 cm 處；當時間分別為25、90、250、480 min 時水分運移至距滴頭正下方豎直方向5、10、15、20 cm 處。土壤水分平均水平運移速率分別為0.2、0.09、0.05 cm/min，平均垂直運移速率分別為0.2、0.11、0.06、0.04 cm/min。垂直方向土壤入滲速率高于水平方向，這是由于土壤水分在垂直方向入滲受到基質勢和重力勢共同作用。隨著時間推移水勢梯度減小，土壤水分在水平和垂直方向平均入滲速率均減小。當濕潤鋒開始交匯時，N2 的土壤含水率迅速升高，灌水結束后土壤含水率由0.08 cm3/cm3提高到0.3 cm3/cm3。灌水結束時地表N4點的含水率小于0.35 cm3/cm3，地表以下N5的含水率大于0.35 cm3/cm3，N6的含水率大于N3，可以看出，滴灌點源入滲過程中，距滴頭水平和垂直相同距離測點含水率有差距，垂直方向含水率高于水平方向，隨著水分擴散，水平和垂直方向土壤含水率的差距逐漸減小。這是因為在滴灌入滲過程中基質勢是土壤水分水平運動的主要動力，而土壤水分向下運動的動力除基質勢外還有重力勢。N7 的含水率小于N2，因為濕潤鋒的交匯先從地表開始，接著逐漸向下進行交匯，交匯時，交匯面處土壤含水率明顯大于下層土壤。由測點N9、N11、N14可以看出，土壤水分運移到該點的時間分別為60、225、480 min。

3.2 流量對土壤濕潤體內水分分布的影響

對雙點源土壤入滲不同滴頭流量供水情況下的土壤濕潤體內含水率分布狀況，采用試驗實測的方法來獲得土壤含水率分布剖面相對比較困難，而采用數(shù)值模擬的方法可以得到可靠的計算結果。圖6 為2 種滴頭流量（1.2、2.4 L/h）、灌水時間（250、600 min）組合情況下的濕潤體內土壤含水率等值線剖面。由圖6 可知，土壤含水率等值線圖類似于濕潤體形狀，土壤水分由滴頭向外擴散，水勢梯度減小，等值線圖由疏變密。灌水時間相同，滴頭附近各節(jié)點處土壤含水量，土壤濕潤體剖面面積，高含水區(qū)面積及同一位置濕潤區(qū)的土壤平均含水率均隨灌水流量的增大而變大。其它條件不變，流量增加，濕潤鋒運移速率加快，水平和垂直距離增加，交匯時間變短，反之亦然。地表因滴頭流量增加而形成積水區(qū)，濕潤鋒水平運移速率隨著積水區(qū)半徑增大而加快，滴頭流量增加使得濕潤鋒與積水區(qū)邊界間水分梯度減小。說明滴頭流量對濕潤體形狀和土壤含水率的分布均有影響。

交匯面含水率分布狀況的研究對合理制定灌溉制度有一定指導意義。圖7所示為灌水結束后，地表滴頭連線中心交匯面縱剖面含水率變化情況。q=2.4 L/h 時，距地表面豎直向下5 cm 處含水率為0.38 cm3/cm3；0～15 cm 范圍內剖面含水率較高，高于0.35 cm3/cm3，未達到土壤飽和含水率；土壤表層向下20 cm 以內土壤含水率緩慢減小，由0.38 cm3/cm3減小到0.35 cm3/cm3；距土壤表層向下20 cm 以下，含水率減小速率增大，豎直向下20 ～30 cm 范圍內，含水率由0.35 cm3/cm3減小到0.08 cm3/cm3。q=1.2 L/h 時，距地表面豎直向下3 cm 處土壤含水率最大，為0.31 cm3/cm3；0～7 cm 范圍內含水率超過0.3 cm3/cm3，距土壤表層向下7 cm 范圍內土壤含水率緩慢減小，由0.31 cm3/cm3減小到0.30 cm3/cm3；距土壤表層向下7 cm 以下，含水率減小速率增加，豎直向下7～18 cm 范圍內，含水率由0.3 cm3/cm3減小到0.08 cm3/cm3。滴頭流量對濕潤體內含水率的變化有一定影響。

4 結 語

（1）本文在非飽和土壤水動力學理論的基礎上試圖建立了雙點源滴灌土壤水分運動數(shù)值模型，應用流體建模有限元分析軟件HYDRUS-3D 進行模擬，其結果與實測數(shù)據(jù)進行對比，一致性較好。在理想狀態(tài)下，可以通過數(shù)值計算方法來探索研究雙點源土壤水分入滲影響因素。

（2）滴灌雙點源入滲土壤含水率的分布同時受灌水時間和滴頭流量的影響。當?shù)晤^流量q=1.2 L/h，灌水時間為268 min時，濕潤鋒開始由點源入滲過渡為雙點源交匯入滲，濕潤體形狀逐漸由2個分離的半橢圓變成花生殼形狀，除靠近滴頭處土壤含水率較高的區(qū)域土壤含水率等值線剖面近似于1/4 橢圓形外，遠離滴頭的濕潤區(qū)內土壤等值線圖融合為一條曲線，并且隨著時間的推移，曲線向下移動，形狀從屋脊形過渡為水平。隨著灌水時間延長，滴頭處土壤水分飽和區(qū)范圍逐漸增大，交匯處土壤含水率增加。灌水時間相同，濕潤鋒運移速率、滴頭附近各節(jié)點處土壤含水量，土壤濕潤體剖面面積，高含水區(qū)面積及同一位置濕潤區(qū)的土壤平均含水率均隨灌水流量的增大而變大。滴頭流量、灌水時間對雙點源滴灌土壤水分入滲均有一定影響。