基于多入滲模型的荒漠砂質土壤積水入滲模擬對比

周 宏

（1.西北師范大學旅游學院，西北師范大學河西走廊研究院，甘肅 蘭州 730070；2.中國科學院西北生態環境資源研究院，中國生態系統研究網絡臨澤內陸河流域研究站，甘肅 蘭州 730000）

入滲是指水分進入土壤的過程，而土壤水分入滲作為水文循環中關鍵環節，對連接地表水資源和降雨以及調控水分分配至關重要［1-2］，許多學者基于不同尺度、不同土壤和氣候區域進行了深入研究［3］。在非均質土壤中提出了用以表征入滲速率、累積入滲量等參數入滲模型，包括Kostiakov 和Horton經驗模型，Philip近似物理模型和Green-Ampt、基于Richard 方程算法的Hydrus-1D/2D/3D 等物理模型。但不同土壤條件中模型應用范圍存在差異，Kostiakov方程簡單，能夠很好地擬合實測入滲數據，但其參數無物理意義，無法解析入滲過程詳細信息［4-5］，Philip是估算最終入滲速率較常用的模型，已在砂質土壤［6］、黏土-壤土農田中得到了驗證［7-8］。Green-Ampt模型主要用于均勻土壤介質中的積水入滲，借助其簡便性和物理基礎，擴展模擬了一系列條件下的入滲、降雨和徑流過程［9-12］。非飽和水分運移過程模擬主要采用基于質量守恒定律和達西定律的Richards 方程來表達［13］，但Richards 方程通常用非線性描述，無法求解析解，導致計算過程較復雜，特別在復雜的初始和邊界條件中。由美國農業部鹽漬土實驗室開發的Hydrus-1D 軟件可廣泛用于模擬飽和-非飽和介質中一維土壤水分運移，其Hydrus-2D/3D 軟件應用于模擬二維或三維土壤水分運移和空間分布，并已獲得大量的驗證和應用［14-16］。

干旱區綠洲荒漠區有許多代表性的景觀單元類型，如荒漠、沙丘和農田等［17］。其中沙丘及其丘間低地在荒漠生態系統水力連接和水循環中扮演著重要的角色［18］，事關荒漠生態系統的結構功能演化和綠洲穩定。在年均降水量小于200 mm干旱區，降水具有降水稀少、變率大、歷時短等特征［19］，盡管降水量少，但確是干旱生態系統中唯一水分補給來源，尤其是儲存在淺層土壤水分調控地表徑流和降雨入滲及其再分配過程［20］。水分入滲作為土壤水文過程重要組成部分，許多學者已針對不同尺度的干旱荒漠生態系統土壤包氣帶水分入滲方面開展了廣泛研究［21］。積水入滲是由于入滲強度小于降雨強度導致一種壓力入滲［22］，尤其在土壤結皮較厚的干旱區，水分停留時間延長，相應入滲量會降低，易導致積水入滲現象發生，而在丘間低地入滲速率降低，會導致下墊面土壤含水率提高［23］。

國內外有關土壤入滲過程研究方法較多，如土柱回填法、單雙環入滲、Hood土壤入滲儀等，其中采用單環定水頭對水分入滲進行原位測量，可用于估算土壤水分入滲特性［2,24-25］。然而，由于野外試驗需要花費大量的時間和費用，而簡單、可靠的滲透試驗數據處理方法顯得尤為必要。但由于土壤異質性，導致大量入滲模型在模擬土壤水分變化時仍然有很大局限性。因此，不同時空條件下，土壤水分運移過程模型的識別和選擇仍需不斷嘗試，基于此，本研究通過野外原位入滲試驗，分析砂質土壤中水流入滲過程并評價不同模型在土壤水分入滲模擬中的性能，旨在尋找降雨條件下適宜模型，以便合理預測沙丘積水入滲過程，提升對荒漠沙丘局地土壤包氣帶水分運移及交換過程認知。

1 材料與方法

1.1 研究區概況

研究在巴丹吉林沙漠北部邊緣的中國科學院臨澤內陸河流域研究站進行，該站屬于中國生態系統研究網絡（CERN），地理位置位于100°9′30″E，39°24′41″N，海拔1381 m。年均氣溫7.5 ℃，最高氣溫39 ℃，最低氣溫-27 ℃，年均潛在蒸發量變化為1900～2088 mm，全年濕度變化為7.3%～80.9%。多年平均降雨量為120 mm，高峰出現在7—9 月。降雨是研究區域最敏感的氣象因子，2006—2017年該區域共接收到降雨事件約500次（圖1），地下水埋深約為4.9 m，地下水的飽和毛管上升不會影響表層土壤水分。

圖1 試驗區多年降雨狀況和干旱區沙丘-丘間低地地貌Fig.1 Location of experimental site,rainfall and image of a mature dune slack in arid region

1.2 原位試驗與設計

在選定的試驗區開展單環入滲試驗，首先去除土壤表面可能阻礙入滲環插入的殘留物，然后用落錘將鐵環夯實，嵌入深度約為5 cm（鐵環尺寸：直徑20 cm，高度為30 cm）。此外，要注意保持鐵環邊緣與土壤表面垂直，并盡量減少對土壤擾動，試驗開始前在環內墊上濾紙，以防止水流破壞土壤表面，造成水頭不均勻（圖2），并沿著入滲環下端中心位置插入土壤水分觀測探頭，探頭插入土壤剖面深度依次為10 cm、20 cm、40 cm、60 cm、80 cm，數據采集器連續收集土壤水分變化值，時間間隔1 min（圖2）。

圖2 單環試驗示意圖及土壤包氣帶剖面和探頭深度Fig.2 The profiles of soil vadose zone and schematic diagram of single ring experiment

1.3 土壤樣品收集與分析

在研究區0～80 cm內用土鉆分層采集原狀土壤樣品，土樣均來自入滲剖面中心，然后烘干法求土壤容重，浸泡法測飽和含水量，激光粒度儀獲取土壤粒徑組成，定水頭法求飽和導水率KS，土壤水力參數是水分運移模擬的關鍵因素，主要特征數據見表1。

表1 試驗場地土壤水力學和物理特性Tab.1 Soil hydraulic and structural properties in experimental site

2 模型描述

2.1 Hydrus-1D/3D模型理論

三維水流運動可以用如下Richards方程來描述:

式中：θ為體積含水率（cm3·cm-3）；h為壓力水頭（cm）；t為時間（h）；K（h）為非飽和導水率函數（cm·d-1）；z為向上正的空間坐標（m）；S為匯源項，通常表示根吸水率（S-1），本試驗在裸露沙丘進行，不考慮根系吸水。其中簡化的一維水流運動方程如下：

2.1.1 土壤水力性質 用Van Genuchten-Mualem（VG）模型描述土壤水分特征曲線，其土壤水力函數如下：

式中：θr和θs分別為土壤殘余含水量和土壤飽和含水量（cm3·cm-3）；Se為飽和度（無量綱）；Ks為飽和導水率（cm·d-1）；n為孔徑分布指數（無量綱）；m=1-1/n為土壤水分特征曲線參數；φ為土壤基質勢；l為孔隙連通性參數（無量綱）?；赩G模型土壤水力參數見表2。

表2 試驗地不同剖面土壤水力參數組成Tab.2 Soil parameters for different soil layers in the experiment site

2.1.2 初始和邊界條件

（1）基于Hydrus-2D二維計算模擬

模擬區域入滲面以上DA、BC為不透水邊界；地表DE、FC為大氣邊界，下底面AB為自由排水邊界；EF 在入滲開始后很快達到飽和，為定水頭邊界，其邊界和初始條件如下：

①初始條件：

式中：θ0為初始含水量（cm3·cm-3）。

②邊界條件：

式中：E(t)為入滲速率；H(t)為恒定水頭。

（2）基于Hydrus-3D三維計算模擬

三維模擬區域入滲面OPNM 為大氣邊界；側面OPGH、GPMJ、MNIJ 和HION 為不透水邊界，下底面GHIJ為自由排水邊界；以直徑R=20 cm入滲圓面開始后很快達到飽和，為定水頭邊界，其邊界和初始條件如下：

①初始條件：

②邊界條件：

二維與三維模擬均采用1 cm 壓力水頭作為單環入滲下的水流邊界條件（圖3），并將模擬區域離散成徑向網格間距為2 cm、垂直網格間距為1 cm的有限元，總模擬時長12 h。

圖3 Hydrus-2D/3D模擬區域及邊界條件Fig.3 The Hydrus-2D/3D simulation domain and boundary conditions

2.2 入滲模型理論

2.2.1 Green-Ampt Green-Ampt 模型最初是利用達西定律對定水頭條件下土柱積水入滲進行分析，入滲方程如下：

式中：I為累積入滲量（cm）；i為入滲率（cm·min-1）；Ks為飽和導水率；Sf為濕潤鋒面吸力；Zf為概化的濕潤鋒深度（cm）；θi為初始土壤含水量（cm3·cm-3）。

分層均質土壤入滲時，通過整合Green-Ampt模型推導出累積入滲量和濕潤鋒隨時間變化的隱函數，分別表示如下：

式中：D為土層厚度（cm）；下標i、s、j、N分別表示初始狀態、飽和狀態、土層數和飽和層數。

式中：tN為濕潤鋒到達第N層界面所需時間（min）。

2.2.2 Philip Philip模型最初由Philip提出，累積入滲量和入滲率表示為：

式中：A是關于土壤性質和水分導水系數函數（cm·min-1）；S為吸水率，是土壤基質勢函數（cm·min-0.5）。

式中：Δθ為飽和含水量和初始含水量的差值。

2.2.3 Kostiakov-Lewis 修正后長周期的Kostiakov 模型描述為［26］：

式中：B、C分別為方程參數（C＞0,0＜B＜1）；i（t）、if分別為入滲率和穩定入滲率。

2.3 模擬效果評價

模擬值與觀測值比較用以評價模型性能，本研究以決定系數（R2）和均方根誤差（RMSE）作為最終評價指標。

式中：Xobs,i為觀測值；Xmod,i為模擬值；為觀測樣本均值，obs,i代表第i（1,2,…,n）個觀測值；n為觀測數據總個數。RMSE 指標值越小，表明模擬誤差越小，而R2的指標值越接近1，表明模擬精度越高，擬合度越好。

3 結果與分析

3.1 實測數據

結果表明，試驗開始后入滲速率持續下降，300 min 左右達到穩定入滲速率0.57 cm·min-1（圖4a），累積入滲量隨時間增加而遞增，觀測時段內累積入滲量達到398 cm（圖4b），入滲開始后濕潤鋒同時沿豎直和水平兩個方向快速推進，200 min后推進速度趨緩，入滲結束后，水平方向濕潤鋒距離較垂直方向濕潤鋒距離增加了45%（圖4c）?？梢园l現冪函數可以較好的擬合實測數據（R2=0.85）。然而，冪函數只是估算土壤水分入滲的一個簡單的經驗模型，無具體物理參數意義。此外，結果表明，在10～80 cm深度內土壤各個剖面含水量差異明顯（圖4d），達到穩態土壤含水量隨土層深度增加而降低，但20 cm剖面的最大含水率明顯低于其他剖面，且均低于飽和含水量，這可能與容重、粒徑構成等土壤性質差異有關。

圖4 觀測時段的入滲速率、累積入滲量、濕潤鋒距離和土壤含水量變化Fig.4 Dynamic of observed infiltration rate,cumulative infiltration,depth of wetting front and soil water content during the study period

3.2 觀測與模擬結果比較

3.2.1 入滲速率 4 種水文模型的入滲率模擬值與觀測值對比如圖5 所示。結果表明，入滲率的R2系數均值在0.75～0.95，Philip 的R2系數最高（0.95），Kostiakov 的模擬結果略大于觀測值，相反Green-Ampt 的模擬結果略低于觀測值，尤其在入滲后期。Hydrus-1D 模擬結果與實測匹配較低，決定系數R2較?。?.79）?？傊?，荒漠砂質土壤積水入滲條件下，Philip 較其他入滲模型能更好地描述土壤水分入滲過程。

圖5 入滲速率模擬值與實測值比較Fig.5 Comparison of simulated infiltration rate from models with observed result

3.2.2 累積入滲量 模擬與觀測累積入滲量關系如圖6 所示，結果表明，4 種模型累積入滲量模擬值均與實測數據有較好的相關性。而Philip 模型R2最高，RMSE 系數最低，Kostiakov 模型RMSE 次之（表3），但Hydrus-1D模型下RMSE值系數相對較高，表明高估了累積入滲量，尤其在入滲初期，但4種入滲模型R2系數均達到了0.98以上，因此，4種入滲模型均是預測砂土累積入滲量較為有效的模型，尤其是Philip模型。

表3 模擬結果擬合優度參數Tab.3 Goodness-of-fit parameters for simulation results with models

圖6 累積入滲量模擬值與實測值比較Fig.6 Comparison of simulated cumulative infiltration from models with observed result

3.2.3 濕潤鋒距離 模擬與觀測濕潤鋒距離之間的關系如圖7 所示，結果表明，Green-Ampt 和Hydrus-1D模型的濕潤鋒距離模擬值與觀測值差異較大，尤其是Hydrus-1D 明顯高估了濕潤鋒推進距離，其均值R2小于Philip 模型。一種可能的解釋是Hydrus-1D 忽略了土壤水側滲和氣流，以往研究表明，Hydrus-1D高估了濕潤區的蓄水能力，但綜合RMSE和R2均值結果可知，Philip模型對砂土濕潤鋒推進距離的預測效果較好。

圖7 濕潤鋒距離模擬值與觀測值比較Fig.7 Comparison of simulated wetting front depth from models with observed result

3.2.4 土壤含水量與Hydrus-2D 模擬 基于Hydrus-2D模擬土壤含水量如圖8所示。結果表明，土壤水分模擬值與觀測值擬合度較低，R2均值僅為0.82，尤其在表層0～10 cm，R2僅為0.65，綜合0～80 cm的RMSE和R2均值結果可知，Hydrus-2D模型并不能很好地預測積水入滲條件下土壤水分變化過程。

圖8 基于Hydrus-2D模擬土壤含水量與實測值比較Fig.8 Comparison of simulated soil water content by Hydrus-2D model with observed result

3.2.5 土壤含水量與Hydrus-3D 模擬 基于Hydrus-3D 模擬土壤含水量如圖9 所示，結果表明，土壤含水量模擬值與觀測值基本一致，0～80 cm 土壤剖面的R2均值為0.93，RMSE 均值為0.02 cm3·cm-3（表4）。盡管較Hydrus-2D 相比，Hydrus-3D 對入滲后期土壤含水量預測較低，但是綜合考慮R2和RMSE值，Hydrus-3D是描述和模擬砂土積水條件土壤含水量較為理想的選擇。

表4 H ydrus-2D和Hydrus-3D模擬結果擬合度參數Tab.4 Goodness-of-fit parameters for simulation results between Hydrus-2D and Hydrus-3D

圖9 基于Hydrus-3D模擬土壤含水量與實測值比較Fig.9 Comparison of simulated soil water content by Hydrus-3D model with observed result

4 討論

4.1 土壤水分入滲特性

單環積水入滲儀是測量土壤水分入滲速率的最常用方法之一，已被廣泛用于評估和監測降雨入滲進入土壤過程。研究表明，隨著時間的延長，入滲速率隨著入滲量的增加而降低，最終達到穩定入滲率［27］。本研究發現，砂土的穩定入滲率為0.55 cm·min-1，這與前人研究砂土入滲速率范圍為0.2～0.77 cm·min-1結果相一致［28］，結果差異可能由于土壤質地、初始土壤含水量和地表覆蓋條件導致。土壤水分運移過程中，非飽和區水分濕潤鋒推進動力學特性研究具有重要的意義，研究指出，單環入滲受環下水分橫向運動影響，其基質勢通量在濕潤鋒處占主導地位［29］，本研究發現濕潤鋒距離水平方向較垂直方向大。此外，在入滲發生的前100 min內，濕潤鋒推進速度隨時間增加而減小，之后保持恒定不變（圖10），一種可能的解釋是砂土引發的漏斗流在水平方向上受基質勢驅動，而在垂直方向上受重力勢驅動。

圖10 垂直和水平兩個方向濕潤鋒推進速度Fig.10 Comparison of observed the velocity of the wetting front with two directions

先前研究表明，Philip 模型可用于砂質土壤入滲速率與入滲時間關系預測［30］，Zolfaghari 等也做了類似報道［31］，然而，Duan等［28］和Ogbe等［32］通過研究入滲模型對砂土水分入滲適應性發現，Hortons模型較適宜累積入滲量預測，Wang等［13］通過模擬與實測值比較，證實Hydrus-1D 能夠預測層狀土柱砂土積水條件下的入滲速率。然而，由于模型參數和性能隨土壤條件和時間的不同而存在差異，需要大量試驗以確定入滲模型適用范圍和條件。

4.2 土壤剖面水分狀況

土壤水分分布狀況是影響干旱區荒漠植物根系吸水有效性的重要因素，本研究結果表明，土壤含水量表層高于深層，積水入滲導致土壤水分在整個濕潤區呈不均勻非飽和狀態分布，這可能與土層結構對土壤的持水性能影響有關，盡管整個剖面土壤質地為砂土，但由于各層之間粒徑組成比例有差異，導致持水性略有不同，出現土壤含水量空間變異［33-34］。此外，Hydrus-3D較Hydrus-2D對土壤含水量數值模擬結果效果最佳，其優點是考慮了水分在土壤中的三維變飽和多孔介質運動。

5 結論

以自然沙丘中的丘間低地為研究對象，通過單環入滲試驗，以確定土壤水力參數，并分析了砂土在固定水頭滲流區的土壤水分變化對入滲響應，以驗證幾種入滲模型在積水入滲條件下的適應性狀況，取得以下主要結論：

（1）入滲初期入滲速率隨時間急劇下降，之后趨于穩定，整個入滲過程中，水平濕潤鋒的推進速度較垂直方向快，入滲速率、累積入滲量和濕潤鋒深度與時間呈冪函數關系。

（2）受邊界和初始條件以及土壤水力參數設置等因素影響，Hydrus-1D 水文模型在砂土水分入滲模擬中表現欠佳，但當前試驗條件下，通過性能評價指標系數比較，特定土壤環境中Philip 較Green-Ampt、Kostiakov-Lewis 和Hydrus-1D 模型能夠合理描述砂土入滲參數，反映較真實的土壤水動力過程。

（3）單環入滲試驗需要同時考慮環內一維和環外三維水流過程，而Hydrus-3D 模型基本能反映砂土不同剖面三維土壤水分變化?？傊?，結合Philip和Hydrus-3D 水文模型確定砂土入滲是可行的，但此結果只是基于砂土，能否應用其他類型土壤水分入滲過程的評估有待進一步研究。