一維非均質條件下農田土壤水分運移模擬及應用

王佩浩，張 茜，王國芳，張吳平

（1.山西農業大學資源環境學院，山西太谷030801；2.山西農業大學軟件學院，山西太谷030801）

土壤水分作為四水轉化的中樞、植物生長發育的關鍵元素，無論是在自然降水條件還是人為的農事活動等的影響下，都有其重要的研究意義。國內外學者在設定均質土壤條件下已對不同初邊條件下土壤的入滲、再分配過程進行了大量的理論研究[1-5]。但在實際的田間土壤中，水力特征差異顯著的層狀土是田間土壤普遍存在的土壤質地結構，會導致土壤水分在均質土壤和層狀土壤中的運移具有很大的差異，進而導致模擬結果與真實的土壤水分運移過程存在較大差異，影響了土壤水分運移模擬模型在SPAC 系統和農業精準灌溉中的應用[6]。為此，國內外有學者采用HYDRUS 系列模型以及SWMS 等經典的模擬軟件模擬研究了非均質多層土壤中水分的入滲、再分配和排水等過程，得到了理想的模擬結果[7-8]。然而，在實際的土壤水分定量化運移過程中，需要與特定的研究目標結合起來，如研究犁底層特征對大氣降雨入滲和土壤再分配的影響，既需要把犁底層的存在對土壤水分運移過程實時地反映出來，又需要采用不同的關系模式。這種情景下，采用已存在的經典模型限制較多，不能靈活地改變不同變量及參數間的關系，限制了土壤水動力學模型在農業生產實踐中的應用。

本研究采用Richards 方程的基質勢形式描述了田間非均質土壤水分的運移過程，并采用Visual Studio 2017 面向對象的程序語言利用數值差分法求解了多層基質勢方程的數值解，在設定犁底層存在的條件下，模擬了土壤水分的運移過程，揭示出多層質地下農田土壤水分運移的規律特征，以期為研究區域農田土壤水分入滲變化提供依據。

1 材料和方法

1.1 研究區選取

大田土壤水分入滲試驗于2018 年8 月23 日至9 月8 日在太谷縣侯城鄉孟家莊村東北部試驗地進行。試驗地土壤類型為褐土性土，地理位置為東經112°35′，北緯37°26′。

1.2 模型初始值測定

試驗地土壤數據采集深度180 cm，每10 cm 為一層，分層用環刀取適量樣品進行室內試驗測定。測試指標包括土壤粒徑組成、容重、飽和含水量和分層土壤含水量（表1），其測定方法依據《土壤物理實驗指導》[9]。氣象數據取自太谷縣氣象局。

表1 各層土壤粒徑組成、容重及飽和含水量

1.3 模型原理與方法

1.3.1 土壤水分運動方程 本研究忽略土壤中的溫度與氣相對土壤水分的影響，借助一維飽和-非飽和Richards 控制方程描述多層質地下土壤水分運移[10]。

1.3.2 邊界條件的設定

1.3.2.1 初始條件 本研究設定時段初的土壤基質勢剖面為一維垂直土壤水分入滲的初始條件。

1.3.2.2 邊界條件 考慮到研究區域實際的降雨與蒸發過程，設定上邊界條件為二類邊界條件，在不同時段設定相應的土壤入滲速率，即在降雨時設定地表入滲速率，蒸發時設定蒸發速率[11]。

式中，f（t）在降雨時刻表示土壤入滲速率，蒸發時刻表示蒸發速率。

下邊界條件設定為常數邊界條件。

1.3.3 土壤水分特征曲線 通過粒徑組成并借助土壤傳遞函數得到本研究中各層土壤采用Van-Genuchten 模型表示的水分特征曲線初值[12]。

1.4 數值模型參數敏感性分析與確定

采用Visual Studio 2017 面向對象的程序語言，利用數值差分法求解多層基質勢方程的數值解。

1.4.1 敏感性參數的確定 將數值模型模擬結果與實測數據進行對比分析，通過對表征模型精度的主要參數進行敏感性分析，確定參數的調整范圍。土壤質地狀況的差異會引起土體孔隙狀況的不同，導致參數α 和n 發生變化，進而影響土壤水分特征曲線的變化[13-16]。因此，本研究確定的主要參數為α和n。

1.4.2 參數優化 將模型推求的各項土壤物理參數初值代入數值模型中進行模擬，不斷調整參數α和n 的大小，直至模型輸出含水量值與實測含水量（8 月23—30 日）達到最優擬合。參數最優調整依次為：0～10 cm 范圍n 值不變、α 值增加10%；10～50 cm 范圍n 值與α 值同時減少10%；50～120 cm范圍n 值與α 值同時增加10%；120～140 cm 范圍n 值與α 值同時增加30%；140～150 cm 范圍n 值與α值同時增加20%；150～170 cm 范圍n 值與α值同時增加30%；170～180 cm 范圍n 值與α值同時增加20%。參數α和n 的初值和優化值列于表2。

表2 關鍵參數初值與優化值

1.4.3 模型精度的驗證 利用參數優化后的模型進行土壤水分運動變化模擬，模擬的土壤含水量值與實測值（9 月1—7 日）通過均方根誤差[17]（RMSE）和納什系數[18]（NSE）2 個指標進行驗證。二者統計分析結果如圖1 所示。

圖1 結果顯示，本次模擬結果表層0～10 cm模擬效果誤差大，其余各層的模擬效果都很好，表明模型模擬可對實測值進行預測，進而準確地反映出農田土壤水分的運移變化情況。

2 模型應用與分析

土壤的黏粒含量和容重是影響土壤水分入滲的重要因素，二者的變化會使土壤孔隙發生變化，進而引起參數α和n 的變化[13，19-20]。為了更好地檢驗模型關鍵參數α和n 對水分運移模擬的影響以及模型的適用性，本研究進行了如下模擬設定（圖2）。

據太谷縣多年降雨數據，設定模擬初期降雨強度為0.208 cm/h，降水時長24 h，后期設置為自然條件下的蒸發速率（此處輸入的是實測的蒸發速率數據）。下邊界條件設定為常數邊界條件，土壤水分自由下滲。代入數值模型進行模擬，模擬結果如圖3所示。由圖3 可知，總體來看，除第1 天外，在降水后的一段時間內，0～70 cm 土層土壤含水量隨時間增加基本呈逐漸降低趨勢，70～120 cm 土層含水量隨時間增加逐漸增大，120～180 cm 土層土壤含水量隨時間增加逐漸降低。120～180 cm 土層，2 種類型土壤各層含水量基本隨時間增加而降低，且變化趨勢基本相同。考慮該段土壤質地類型主要為粉黏壤和黏壤土，粉黏粒含量較大，表明土壤深度超過120 cm 土層處含水量變化受粒徑組成影響持水穩定，且受上層土壤影響較小，隨時間增加，下層土壤水分在重力和土壤本身吸力作用下入滲，因此，土壤含水量變化基本穩定。

圖3-a 與圖3-b 相比，0～50 cm 處含水量值圖3-a 各層含水量均低于圖3-b，該范圍內αa＞αb、na＞nb，使得圖3-a 的土壤進氣值小于圖3-b，因此，圖3-a 的土壤含水量較圖3-b 低。圖3-b 在30 cm處由于犁底層存在的影響，土壤含水量曲線較圖3-a 出現明顯的彎折，且較無犁底層的土壤含水量曲線分布更為緊密，此處na＜nb，因此，同一點處土水勢圖3-b 較圖3-a 降低，土壤吸力增大，受αa＞αb的影響，含水量圖3-b 高于圖3-a。50～70 cm 處模擬圖3-a 與圖3-b 的各層含水量均隨時間增加而降低，但圖3-a 的變化幅度明顯大于圖3-b，圖3-b 的變幅較小，這也表明具有較小空隙和較多黏粒組成的犁底層能夠很好地降低土壤水分入滲速率和流失速率，具有良好的吸持水能力。

3 結論與討論

本研究首先通過對研究區域進行實地采集與測量得到所需各項土壤物理數據，結合數值理論模型，搭建基于Visual Studio 的一維土壤水分運移模型。其次，借助實測數據對數值模型進行了參數敏感性分析及模擬精度驗證。最后，為更好地檢驗模型關鍵參數α 和n 對水分運移模擬的影響以及模型本身的適用性，通過分設不同性質土層進行模擬。研究結果表明，模型參數敏感性主要受土壤水力特征參數α 和n 的影響，參數值越大，土壤進氣值越低，含水量越小；n 值減小，土水勢降低，土壤吸力增大，含水量增加；在調整數值模型精度時，α 和n 的單個誤差一般為10%～30%；模擬值和測量值的變化過程比較吻合與一致，表明本研究構建的基質勢方程能夠實現對田間一維多層質地土壤水分運移過程的準確模擬。

本研究土壤水分運移模型可以較好地反映土壤水分變化特征，具有較高的應用價值。與其他模型相比，數據輸入與操作相對簡便，整體過程可控，能夠滿足對特定問題的專門研究。模型應用也較好地體現了土壤水力參數變化對模型精度的適用性。本研究后期將考慮在不同時空狀態下土壤水分的運移變化，同時結合植被生長、根系吸水等地上植物情況，完善土壤水分運移模型，以期為更好地研究區域農田土壤水分入滲狀況提供相應的參考依據。