宋孝玉,白 鵬,王 娟,史文娟,王全九,4

(1 西安理工大學(xué) 西北旱區(qū)生態(tài)水利工程國家重點實驗室培育基地,陜西 西安 710048；2 中國科學(xué)院 地理科學(xué)與資源研究所,北京 100101；3 河北省水利科學(xué)研究院,河北 石家莊 050051；4 中國科學(xué)院 西北農(nóng)林科技大學(xué) 水土保持研究所 黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點實驗室,陜西 楊凌 712100)

降雨入滲過程是土壤水循環(huán)的重要組成部分，是農(nóng)業(yè)灌溉、工程設(shè)計、水文過程模擬的基礎(chǔ)。國內(nèi)外學(xué)者對土壤水分入滲過程進(jìn)行了大量的研究[1-2]，建立了一系列物理基礎(chǔ)模型和經(jīng)驗公式模型。國外自1911年建立概念明確、形式簡單的Green-Ampt積水入滲模型以來，相繼建立了Richards方程、Horton 模型、Philip模型、Smith模型等，而國內(nèi)學(xué)者的研究大多是對國外經(jīng)典入滲模型的改進(jìn)和應(yīng)用[3-4]。相對于經(jīng)驗入滲公式，物理模型具有明確的物理意義，便于建立其特征參數(shù)與土壤物理特征間的關(guān)系，因而得到了廣泛應(yīng)用[5]。在物理模型中，無論是基于Darcy定律簡化的Green-Ampt模型，還是基于土壤水分運動方程解析法的Philip模型、Parlange模型等，都只適用于簡單定解條件下的土壤水分運動問題，對于復(fù)雜條件的土壤水分運動問題，數(shù)值解法仍然是最客觀、有效的方法[6]。

目前，Richards方程在國外廣泛應(yīng)用于各種水文、農(nóng)業(yè)、工程施工模擬軟件中，如HYDRUS[7]、SWAP[8]、SWMS-2D[9]等，國內(nèi)對Richards入滲方程研究還不夠重視，有限的研究大多集中在室內(nèi)土柱入滲過程或點源滴灌的模擬應(yīng)用方面[10-13]，在野外降雨條件下的應(yīng)用研究尚比較少[14]。本研究在陜北黃土坡地人工降雨試驗的基礎(chǔ)上，采用有限差分法對Richards方程進(jìn)行數(shù)值模擬求解，分析其在降雨入滲、產(chǎn)流過程中的應(yīng)用情況，以期為該方程在該地區(qū)墑情、徑流預(yù)報以及坡面侵蝕方面的應(yīng)用提供參考。

1 材料與方法

1.1 試驗區(qū)概況

野外降雨試驗于2010-03-2011-09在陜西榆林市米脂縣銀州鎮(zhèn)孟岔村紅棗示范園進(jìn)行，數(shù)據(jù)測試和分析于2012-12完成。米脂縣位于陜西省北部，榆林市中部偏東，屬無定河中游，地理位置為東經(jīng)109°49′～110°29′，北緯37°39′～38°5′；總?cè)丝?1.6萬，其中農(nóng)業(yè)人口占89.4%；總土地面積1 212 km2，川地主要種植玉米，坡地主要種植果樹以及雜糧作物，在國家退耕還林的政策措施下，當(dāng)?shù)剞r(nóng)戶在坡地大面積種植果樹，其中棗樹種植面積最廣。米脂縣屬典型的黃土高原丘陵溝壑區(qū)，屬于溫帶半干旱性氣候區(qū)，以無定河為分水嶺，地勢總體東西高中間低。試驗區(qū)年降雨較少，晝夜溫差較大，多年平均降雨量413 mm，主要集中在7-9月，其中最大年降雨量691 mm，最小年降雨量268.3 mm。區(qū)域土壤以黃綿土為主，土質(zhì)松軟，水土流失嚴(yán)重。試驗區(qū)內(nèi)種植棗樹品種主要為梨棗，樹齡為7～8年，株高2 m左右，株行距為2 m×3 m，栽培模式為矮化密植。

1.2 試驗設(shè)計與測定方法

降雨徑流試驗在徑流小區(qū)內(nèi)進(jìn)行，徑流小區(qū)選擇在地形、坡向、土壤、植被情況等有代表性的地段上，平整土地，去除雜草，使坡面、坡向均勻一致?？紤]到人工降雨用水的便利，小區(qū)選擇在離供水水源地較近的地段。根據(jù)坡面降雨入滲試驗研究需要，共修建10個徑流小區(qū)，坡向為陽坡，坡度測量利用全站儀進(jìn)行測量。徑流小區(qū)寬度均為3 m，四周用石棉瓦做防側(cè)滲處理，埋深30～40 cm，底部修建V字形出水口，通過徑流桶收集徑流。供水水源位于山頂蓄水池，通過管道將用水引至徑流小區(qū)，山頂蓄水池與徑流小區(qū)之間垂直高度為40 m左右，人工降雨利用重力水頭差進(jìn)行。降雨噴頭采用美國雨鳥牌1804散射型噴頭，噴射半徑為4 m，單噴頭最大噴水量為0.8 m3/h，可通過調(diào)節(jié)閥門開度以及噴頭組合方式調(diào)節(jié)降雨強度大小，經(jīng)試驗驗證，降雨均勻度在0.8以上。試驗裝置見圖1。

降雨徑流試驗在無風(fēng)條件下進(jìn)行，均勻布設(shè)雨量筒測定雨強。降雨歷時30 min，降雨開始后，記錄產(chǎn)流時刻。地面產(chǎn)流后，每間隔2 min讀取1次徑流筒體積，并取水樣。

降雨入滲試驗前，利用土鉆法取土，烘干法測定土壤含水量，每隔10 cm取1次土樣，每層3個重復(fù)，測定深度為1 m，取土后回填鉆孔。由于坡長較短，不考慮匯流過程，入滲量根據(jù)水量平衡原理由降雨量和徑流量的差值得到。

圖1 降雨入滲-產(chǎn)流試驗裝置示意圖

土壤水分特征用離心機測定；土壤顆粒分析用篩析法和吸管法相結(jié)合的方法測定；土壤飽和導(dǎo)水率在田間測定，用直徑13 cm有機玻璃取田間原狀土柱，取樣深度約30 cm，土柱下端用細(xì)紗布固封，用馬氏瓶定水頭供水，待下端出水后，每隔10 min測定1次出流水量，直至出流水量恒定為止，試驗重復(fù)3次。

1.3 降雨入滲產(chǎn)流模型及其求解

1.3.1 降雨入滲產(chǎn)流模型 Richards方程以達(dá)西定律和液體連續(xù)方程為基礎(chǔ)描述非飽和土壤水分運動，對于一維垂直入滲情況，其方程形式如下：

(1)

式中：θ為土壤體積含水量，cm3/cm3；t為時間，min；z為垂直方向坐標(biāo)(向上為正)；k(θ)為土壤非飽和導(dǎo)水率，cm/min。

1)初始條件。對于降雨入滲過程而言，初始條件為：

θ(z,0)=θ(z)。

(2)

2)上邊界條件。降雨入滲初期，地表入滲能力大于降雨強度，入滲率等于降雨強度，邊界為第2類邊界，即有：

(3)

式中：D為土壤水分?jǐn)U散率，D=KdΨ/dθ，其中Ψ為土壤基質(zhì)勢，cm；f為土壤入滲率。

當(dāng)表層土壤水分達(dá)到飽和時，入滲率小于降雨強度，地表出現(xiàn)積水，上邊界條件轉(zhuǎn)化為第1類邊界條件，即：

θ(0,t)=θ0(t)或Ψ(0,t)=Ψ0(t)。

(4)

式中：θ0(t)和Ψ0(t)分別為地表土壤含水量和基質(zhì)勢。

3)下邊界條件。對于黃土高原坡地，地下水位極深，下邊界條件為重力排水邊界，相應(yīng)的水分通量為下邊界處的導(dǎo)水率，此時下邊界處的含水率或基質(zhì)勢梯度為0，即：

(5)

式中：L為下邊界深度，cm。

產(chǎn)流量(R)可根據(jù)水量平衡原理由降雨量(P)與累積入滲量(I)之差得到，即：

R=P-I。

(6)

非飽和導(dǎo)水率k(θ)一般通過Van Genuchten[15]和Mualem[16]提出的土壤水分特征曲線與非飽和導(dǎo)水率間的關(guān)系計算，即：

(7)

(8)

(9)

式中：Ks為土壤飽和導(dǎo)水率；Se為有效飽和度；θr為土壤滯留水量；θs為土壤飽和含水量；n和m為形狀系數(shù)，m=1-1/n；h為土壤吸力，cm。

1.3.2 模型求解 模型采用有限差分法求解，根據(jù)土壤水分運動參數(shù)的表達(dá)式、初始條件和邊界條件，按照一定的時間、空間步長進(jìn)行離散，根據(jù)時段初土壤含水率預(yù)報時段末值，利用追趕法求解土壤水分運動參數(shù)的差分方程，得到時段末土壤含水率，當(dāng)時段末的土壤含水量預(yù)報值與計算值的誤差滿足要求，則認(rèn)為迭代收斂，可進(jìn)入下一階段求解。否則，以時段末的土壤含水率計算值為新的預(yù)報值，繼續(xù)迭代，直至滿足誤差要求為止。

2 結(jié)果與分析

2.1 不同降雨強度條件下的降雨入滲產(chǎn)流模擬

試驗在經(jīng)過平整、去除雜草，且坡面、坡向均勻一致的棗樹徑流小區(qū)內(nèi)進(jìn)行，由于5月初試驗地降雨量較少，棗樹剛剛發(fā)芽，因此不考慮降雨截留的影響，選取初始含水率相同的3個徑流小區(qū)開展不同降雨強度的人工降雨入滲、產(chǎn)流試驗，降雨歷時30 min，各場次降雨入滲、產(chǎn)流情況的統(tǒng)計結(jié)果見表1。

用1.3.1的降雨入滲、產(chǎn)流模型對表1不同降雨強度徑流小區(qū)的試驗過程進(jìn)行模擬，用F1-Obs、F2-Obs、F3-Obs和Q1-Obs、Q2-Obs、Q3-Obs分別代表降雨強度為1.3,1.0,0.8 mm/min時的累積入滲量和徑流實測值；以F1-Sim、F2-Sim、F3-Sim和Q1-Sim、Q2-Sim、Q3-Sim分別表示其模擬值。不同降雨強度下各徑流小區(qū)的累積入滲量、徑流量的實測值與模擬值見圖2和圖3。

表1 不同降雨強度下棗樹徑流小區(qū)的降雨入滲、產(chǎn)流情況

圖2 不同降雨強度下徑流小區(qū)累積入滲量實測值與模擬值的比較

降雨強度對土壤入滲的影響主要反映在降雨強度與土壤入滲能力的大小關(guān)系上，降雨過程中，土壤入滲能力隨土壤含水率的增加而減小。降雨初始階段，降雨強度小于土壤入滲能力，實際土壤入滲率就等于降雨強度，降雨全部入滲。隨著土壤入滲能力的減小，某一時刻以后，降雨強度大于土壤入滲能力，地面開始積水。由圖2，3可知，相同條件下，降雨強度越大產(chǎn)流時間越短，入滲達(dá)到穩(wěn)定的時間也越短，累積入滲量和徑流量也越大，降雨歷時為30 min時，Q1-Obs是Q2-Obs的3.34倍，是Q3-Obs的6.06倍。由圖2、3還可知，利用本研究建立的降雨入滲、產(chǎn)流模型對不同降雨強度徑流小區(qū)降雨入滲、產(chǎn)流的模擬結(jié)果與實測值比較吻合，符合土壤實際入滲、產(chǎn)流規(guī)律。

2.2 不同土壤含水率條件下的降雨入滲產(chǎn)流模擬

為了檢驗?zāi)Ｐ蛯Σ煌蕳l件下降雨入滲、產(chǎn)流的模擬效果，試驗繼續(xù)于5月初在棗樹徑流小區(qū)內(nèi)進(jìn)行，選取初始含水率不同的3個徑流小區(qū)開展相同降雨強度(1.0 mm/min)的人工降雨入滲、產(chǎn)流試驗，降雨歷時30 min，各場次降雨入滲、產(chǎn)流情況見表2。

表2 不同土壤含水率下棗樹徑流小區(qū)的降雨入滲、產(chǎn)流情況

繼續(xù)用1.3.1中的降雨入滲、產(chǎn)流模型對不同土壤含水率的試驗過程進(jìn)行模擬，用F1-Obs、F2-Obs、F3-Obs和Q1-Obs、Q2-Obs、Q3-Obs分別代表初始土壤含水量為0.08,0.15,0.24 cm3/cm3時的累積入滲量和累積徑流量實測值；用F1-Sim、F2-Sim、F3-Sim和Q1-Sim、Q2-Sim、Q3-Sim分別代表其模擬值。得到不同土壤初始含水率條件下，3個徑流小區(qū)累積入滲量、徑流量實測值與模擬值，結(jié)果見圖4和圖5。從圖4、5可以看出，土壤初始含水率對降雨入滲、產(chǎn)流過程有重要影響。土壤初始含水率越大，產(chǎn)流開始時間越早，平均入滲率越小，趨于穩(wěn)定入滲的時間也越短，相同時間內(nèi)，累積產(chǎn)流量越大，根據(jù)實測徑流值可知，Q3-Obs是Q2-Obs的1.65倍，是Q1-Obs的3.35倍。由圖4、5還可知，本研究的降雨入滲、產(chǎn)流模型對不同初始土壤含水率徑流小區(qū)降雨入滲、產(chǎn)流的模擬結(jié)果與實測值比較吻合，符合土壤實際入滲、產(chǎn)流的規(guī)律。

圖4 不同初始土壤含水率條件下徑流小區(qū)累積入滲量實測值與模擬值的比較

2.3 模型誤差分析

由圖2～5可知，本研究1.3.1中的降雨入滲產(chǎn)流模型對徑流小區(qū)降雨入滲、產(chǎn)流數(shù)值的模擬結(jié)果與實測值比較吻合，均符合土壤實際入滲、產(chǎn)流的規(guī)律。為了定量分析該模型的模擬精度，對不同時刻入滲產(chǎn)流模擬值和實測值的誤差進(jìn)行分析，結(jié)果見表3。

表 3 徑流小區(qū)降雨入滲產(chǎn)流模擬值與實測值的誤差分析

由表3可知，累積徑流量模擬結(jié)果除試驗1的相對誤差為12.5%外，其余相對誤差均在10%以內(nèi)，所有試驗的平方根誤差均小于0.4 mm；累積入滲量相對誤差均小于3%，平方根誤差均小于0.3 mm，模擬結(jié)果比較理想。由表3還可以看出，產(chǎn)流起始時間的模擬值普遍滯后于實測值，滯后時間最少0.8 min，最多2.55 min，這是由坡面飽和導(dǎo)水率的空間變異性和降雨的不均勻性所致。由于自然界包氣帶土壤結(jié)構(gòu)的復(fù)雜性，土壤表層體積質(zhì)量以及大孔隙的空間差異都可能引起飽和導(dǎo)水率的不同，從而導(dǎo)致局部提前產(chǎn)流。另外，由于降雨的不均勻性，地表在雨強比較大的區(qū)域首先開始產(chǎn)流，這也是導(dǎo)致模擬產(chǎn)流起始時間滯后的原因之一。

3 小 結(jié)

本研究運用Richards入滲方程結(jié)合水量平衡方程模擬了黃土坡地不同降雨強度和不同土壤初始含水率條件下徑流小區(qū)的降雨入滲、產(chǎn)流過程，通過與實測值的對比，表明該模型模擬結(jié)果合理，可作為該地區(qū)降雨入滲、產(chǎn)流的計算公式。

模擬產(chǎn)流起始時間普遍滯后于實測值，滯后時間最少0.8 min，最多2.55 min。這是由于坡面飽和導(dǎo)水率的空間變異性以及降雨的不均勻性引起的，關(guān)于這2種原因?qū)Ξa(chǎn)流起始時刻的具體影響規(guī)律還有待于進(jìn)一步研究。

[參考文獻(xiàn)]

[1] 王文焰，汪志榮，王全九，等．黃土中Green-Ampt入滲模型的改進(jìn)和驗證 [J].水利學(xué)報，2003，34(5)：30-34．

Wang W Y,Wang Z R,Wang Q J,et al.Improvement and evaluation of the Green-Ampt model in loess soil [J].Journal of Hydraulic Engineering,2003,34(5):30-34.(in Chinese)

[2] 毛麗麗，張心平，雷廷武，等．用水平土柱與Green-Ampt模型方法測量土壤入滲性能的原理與誤差 [J].農(nóng)業(yè)工程學(xué)報，2007，23(12)：64-68．

Mao L L,Zhang X P,Lei T W,et al.Principles and errors of measuring in-filtrability with horizontal soil column and Green-Ampt model [J].Transactions of the Chinese society of Agricultural Engineering,2007,23(12):64-68.(in Chinese)

[3] 李 毅，王全九，邵明安，等．Green-Ampt入滲模型及其應(yīng)用 [J].西北農(nóng)林科技大學(xué)學(xué)報：自然科學(xué)版，2007，35(2)：225-230．

Li Y,Wang Q J,Shao M A,et al.Green-Ampt model and its application [J].Journal of Northwest A&F University:Natural Science Edition,2007,35(2):225-230.(in Chinese)

[4] 梁海軍，劉作新，王振營，等．地下滲灌土壤水分運動數(shù)值模擬 [J].農(nóng)業(yè)工程學(xué)報，2008，24(10)：56-60．

Liang H J,Liu Z X,Wang Z Y,et al.Numerical simulation of soil water movement under infiltration irrigation [J].Transactions of the Chinese society of Agricultural Engineering,2008,24(10):56-60.(in Chinese).

[5] 王全九，來劍斌，李 毅．Green-Ampt模型與Philip入滲模型的對比分析 [J].農(nóng)業(yè)工程學(xué)報，2002，18(2)：13-16．

Wang Q J,Lai J B,Li Y.Comparison of Green-Ampt model with Philip infiltration model [J].Transactions of the Chinese society of Agricultural Engineering,2002,18(2):13-16.(in Chinese).

[6] 邵明安，王全九，黃明斌．土壤物理學(xué) [M].北京：高等教育出版社，2006：114-115．

Shao M A,Wang Q J,Huang M B. Soil physics [M].Beijing:Higher Education Press,2006:114-115.(in Chinese).

[7] Simunek J,van Genuchten M T,Sejna M.The HYDRUS-1D software package for simulating the movement of water,heat,and multiple solutes in variably saturated media,version 3.0,HYDRUS software series 1 [M].California,USA：[s.n.],2005:15-189.

[8] Van Dam J C,Groenendijk P,Hendriks R F A,et al.SWAP version 3.2：Theory description and user manual [M].Wageningen,Netherlands:Alterra,2008：11-203.

[9] Jacquesa D,J,Mallantsa D,et al.Modeling coupled water flow,solute transport and geochemical reactions affecting heavy metal migration in a podzol soil [J].Geoderma,2008,145(3/4):449-461.

[10] 周云成，何 娜，王 展，等．變水頭邊界條件下土壤水分運動數(shù)值模擬 [J].沈陽農(nóng)業(yè)大學(xué)學(xué)報，2005,36(4):454-457．

Zhou Y C,He N,Wang Z,et al.Numerical simulation of movement of soil moisture under variable water head boundary [J].Journal of Shenyang Agricultural University,2005,36(4):454-457.(in Chinese)

[11] 李 毅，王全九，王文焰，等．入滲、再分布和蒸發(fā)條件下一維土壤水運動的數(shù)值模擬 [J].灌溉排水學(xué)報，2007，26(1)：5-8．

Li Y,Wang Q J,Wang W Y,et al.Mathematical simulation of soil water movement under infiltration,redistribution and evaporation [J].Journal of Irrigation and Drainage,2007,26(1):5-8.(in Chinese)

[12] 聶衛(wèi)波，馬孝義，王術(shù)禮．溝灌土壤水分運動數(shù)值模擬與入滲模型 [J].水科學(xué)進(jìn)展，2009，20(5)：668-676．

Nie W B,Ma X Y, Wang S L.Infiltration model and numerical simulation of the soil water movement in furrow irrigation [J].Advances in Water Science,2009,20(5):668-676.(in Chinese)

[13] 吳夢喜．飽和-非飽和土中滲流Richards方程有限元算法 [J].水利學(xué)報，2009，40(10)：1274-1280．

Wu M X.Finite-element algorithm for Richards equation for saturated-unsaturated seepage flow [J].Journal of Hydraulic engineering,2009,40(10):1274-1280.(in Chinese)

[14] 郭向紅，孫西歡，馬娟娟．降雨灌溉蒸發(fā)條件下蘋果園土壤水分動態(tài)數(shù)值模擬 [J].農(nóng)業(yè)機械學(xué)報，2009，40(11)：68-73．

Guo X H,Sun X H,Ma J J.Numerical simulation for root zone soil moisture movement of apple orchard under rainfall-irrigation-evaporation [J].Transactions of the Chinese Society for Agricultural Machinery,2009,40(11):68-73.(in Chinese)

[15] Van Genuchten.A close-form equation for predicting the hydraulic conductivity of unsaturated soils [J].Soil Science,1980,44:892-898.

[16] Mualem Y.A new model for predicting the hydraulic conductivity of unsaturated porous media [J].Water Resource Research,1976,12(3):513-522.