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

周 宏

（1.西北師范大學(xué)旅游學(xué)院，西北師范大學(xué)河西走廊研究院，甘肅 蘭州 730070；2.中國科學(xué)院西北生態(tài)環(huán)境資源研究院，中國生態(tài)系統(tǒng)研究網(wǎng)絡(luò)臨澤內(nèi)陸河流域研究站，甘肅 蘭州 730000）

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

圖1 試驗(yàn)區(qū)多年降雨?duì)顩r和干旱區(qū)沙丘-丘間低地地貌Fig.1 Location of experimental site,rainfall and image of a mature dune slack in arid region

1.2 原位試驗(yàn)與設(shè)計

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

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

1.3 土壤樣品收集與分析

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

表1 試驗(yàn)場地土壤水力學(xué)和物理特性Tab.1 Soil hydraulic and structural properties in experimental site

2 模型描述

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

三維水流運(yùn)動可以用如下Richards方程來描述:

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

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

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

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

2.1.2 初始和邊界條件

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

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

①初始條件：

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

②邊界條件：

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

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

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

①初始條件：

②邊界條件：

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

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

2.2 入滲模型理論

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

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

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

式中：D為土層厚度（cm）；下標(biāo)i、s、j、N分別表示初始狀態(tài)、飽和狀態(tài)、土層數(shù)和飽和層數(shù)。

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

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

式中：A是關(guān)于土壤性質(zhì)和水分導(dǎo)水系數(shù)函數(shù)（cm·min-1）；S為吸水率，是土壤基質(zhì)勢函數(shù)（cm·min-0.5）。

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

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

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

2.3 模擬效果評價

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

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

3 結(jié)果與分析

3.1 實(shí)測數(shù)據(jù)

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

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

3.2 觀測與模擬結(jié)果比較

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

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

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

表3 模擬結(jié)果擬合優(yōu)度參數(shù)Tab.3 Goodness-of-fit parameters for simulation results with models

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

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

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

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

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

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

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

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

4 討論

4.1 土壤水分入滲特性

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

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

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

4.2 土壤剖面水分狀況

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

5 結(jié)論

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

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

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

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