一維非均質(zhì)條件下農(nóng)田土壤水分運移模擬及應(yīng)用

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

（1.山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院，山西太谷030801；2.山西農(nóng)業(yè)大學(xué)軟件學(xué)院，山西太谷030801）

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

本研究采用Richards 方程的基質(zhì)勢形式描述了田間非均質(zhì)土壤水分的運移過程，并采用Visual Studio 2017 面向?qū)ο蟮某绦蛘Z言利用數(shù)值差分法求解了多層基質(zhì)勢方程的數(shù)值解，在設(shè)定犁底層存在的條件下，模擬了土壤水分的運移過程，揭示出多層質(zhì)地下農(nóng)田土壤水分運移的規(guī)律特征，以期為研究區(qū)域農(nóng)田土壤水分入滲變化提供依據(jù)。

1 材料和方法

1.1 研究區(qū)選取

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

1.2 模型初始值測定

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

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

1.3 模型原理與方法

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

1.3.2 邊界條件的設(shè)定

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

1.3.2.2 邊界條件 考慮到研究區(qū)域?qū)嶋H的降雨與蒸發(fā)過程，設(shè)定上邊界條件為二類邊界條件，在不同時段設(shè)定相應(yīng)的土壤入滲速率，即在降雨時設(shè)定地表入滲速率，蒸發(fā)時設(shè)定蒸發(fā)速率[11]。

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

下邊界條件設(shè)定為常數(shù)邊界條件。

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

1.4 數(shù)值模型參數(shù)敏感性分析與確定

采用Visual Studio 2017 面向?qū)ο蟮某绦蛘Z言，利用數(shù)值差分法求解多層基質(zhì)勢方程的數(shù)值解。

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

1.4.2 參數(shù)優(yōu)化 將模型推求的各項土壤物理參數(shù)初值代入數(shù)值模型中進行模擬，不斷調(diào)整參數(shù)α和n 的大小，直至模型輸出含水量值與實測含水量（8 月23—30 日）達到最優(yōu)擬合。參數(shù)最優(yōu)調(diào)整依次為：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%。參數(shù)α和n 的初值和優(yōu)化值列于表2。

表2 關(guān)鍵參數(shù)初值與優(yōu)化值

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

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

2 模型應(yīng)用與分析

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

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

圖3-a 與圖3-b 相比，0～50 cm 處含水量值圖3-a 各層含水量均低于圖3-b，該范圍內(nèi)αa＞αb、na＞nb，使得圖3-a 的土壤進氣值小于圖3-b，因此，圖3-a 的土壤含水量較圖3-b 低。圖3-b 在30 cm處由于犁底層存在的影響，土壤含水量曲線較圖3-a 出現(xiàn)明顯的彎折，且較無犁底層的土壤含水量曲線分布更為緊密，此處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 結(jié)論與討論

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

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