(北京林業(yè)大學(xué) 水土保持學(xué)院，北京 100083)

1 研究背景

土壤水分入滲是指水由土壤表面滲入土壤內(nèi)部的過(guò)程，是土壤水分運(yùn)動(dòng)中心環(huán)節(jié)[1]，研究土壤水分入滲規(guī)律對(duì)合理利用水資源，指導(dǎo)農(nóng)林業(yè)生產(chǎn)具有重要意義。相比于耗時(shí)長(zhǎng)、花費(fèi)大的傳統(tǒng)田間試驗(yàn)，精確的數(shù)值模型模擬能夠極大提高工作效率。基于Richards原理的HYDRUS-1D模型因其模擬結(jié)果的準(zhǔn)確性和廣泛的適用性，再加上便捷靈活的操作及集成參數(shù)優(yōu)化功能而被廣泛應(yīng)用于土壤水分運(yùn)動(dòng)研究中[2-4]。大量學(xué)者通過(guò)試驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)的對(duì)比發(fā)現(xiàn)此模型在土壤水分運(yùn)動(dòng)研究中具有很好的適應(yīng)性[5-9]。實(shí)際模擬過(guò)程中，模擬參數(shù)的確定往往是最關(guān)鍵的環(huán)節(jié)，參數(shù)的準(zhǔn)確性直接影響著模型率定效率以及模型輸出結(jié)果的可靠性[10-11]。因此，對(duì)土壤水力參數(shù)進(jìn)行敏感性分析，明確參數(shù)擾動(dòng)對(duì)入滲特性的影響規(guī)律，有助于提高模型模擬精度，便于更加高效的應(yīng)用。

模型的模擬精度除受到模型參數(shù)變化的影響外，外界輸入條件也不可忽略[12]。已有的相關(guān)研究主要集中于定水頭入滲條件下的參數(shù)變化對(duì)輸出結(jié)果的影響，如王志濤等[13]分析了Van Genuchten模型(簡(jiǎn)稱VG模型)參數(shù)變化對(duì)粉壤土入滲特性的影響；范嚴(yán)偉等[14]研究定水頭入滲情況下HYDRUS模型土壤水力參數(shù)變化對(duì)不同質(zhì)地土壤入滲特性的影響；劉鋼等[15]研究發(fā)現(xiàn)飽和含水量θs、飽和導(dǎo)水率Ks、進(jìn)氣吸力倒數(shù)α和孔隙尺寸分布指數(shù)n這4個(gè)參數(shù)對(duì)壓力水頭和溶質(zhì)運(yùn)移的影響較為顯著；霍思遠(yuǎn)等[16]研究表明針對(duì)不同的輸出結(jié)果參數(shù)敏感性不同，并提出參數(shù)敏感性分析及調(diào)參過(guò)程應(yīng)結(jié)合土壤參數(shù)的物理意義展開。以上研究多在定水頭入滲條件下進(jìn)行，未考慮到降雨入滲作為邊界條件的情況。相比于定水頭入滲，自然降雨下的入滲過(guò)程更符合實(shí)際情況，研究結(jié)果也更有助于實(shí)際應(yīng)用。張海闊等[17]以南京市2 a一遇設(shè)計(jì)暴雨過(guò)程作為上邊界輸入條件，主要在土層的垂直方向及水平方向研究參數(shù)擾動(dòng)對(duì)輸出結(jié)果的影響。其結(jié)果表明降雨入滲條件下土壤水力參數(shù)擾動(dòng)對(duì)入滲規(guī)律的影響與定水頭研究結(jié)果有較大差異。目前尚無(wú)降雨入滲條件下，參數(shù)擾動(dòng)對(duì)入滲隨時(shí)間變化過(guò)程的影響研究。

本研究以北京市2 a一遇設(shè)計(jì)暴雨過(guò)程為例作為上邊界輸入條件，選取HYDRUS-1D模型中的5個(gè)土壤水力參數(shù)(θs、Ks、α、n、殘余含水量θr)，研究前鋒雨型下單因素?cái)_動(dòng)對(duì)地表10 cm層土壤水分動(dòng)態(tài)分布、累積入滲量及濕潤(rùn)鋒運(yùn)移距離的影響，針對(duì)不同的輸出結(jié)果分析參數(shù)的敏感性。旨在完善降雨入滲條件下HYDRUS-1D模型參數(shù)敏感性研究規(guī)律，量化參數(shù)敏感性對(duì)于不同輸出結(jié)果的影響，明確參數(shù)敏感性隨時(shí)間變化規(guī)律，以便在實(shí)際應(yīng)用中注意高敏感性參數(shù)的精度對(duì)模擬結(jié)果的影響，提高模擬結(jié)果的準(zhǔn)確性及模型使用的效率。

2 模型原理及模型設(shè)計(jì)

2.1 模型原理

HYDRUS-1D模型僅考慮一維垂向運(yùn)移時(shí)的土壤水分運(yùn)動(dòng)[18]，采用Richards方程描述一維非飽和土壤水流的運(yùn)動(dòng),即

(1)

式中：t為時(shí)間(d);h為壓力水頭(cm);θ為土壤體積含水率;z為垂向空間坐標(biāo)(cm);K(h)為非飽和導(dǎo)水率(cm/d);s(z，t)為根系吸水速率(cm/d)。

該方程是由達(dá)西定律和質(zhì)量守恒定律推導(dǎo)而來(lái)，對(duì)于土壤水分運(yùn)動(dòng)過(guò)程的描述更加準(zhǔn)確，更適合分析土壤水分動(dòng)態(tài)變化特征。

本文使用應(yīng)用廣泛的VG模型擬合土壤水分特征曲線θ(h)和非飽和導(dǎo)水率K(h)，其方程為：

(2)

(3)

其中,

(4)

式中：m為經(jīng)驗(yàn)參數(shù)，n>1，m=1-1/n；l為經(jīng)驗(yàn)參數(shù)，取0.5。

2.2 數(shù)值模型設(shè)計(jì)

2.2.1 離散化

模擬土層為均質(zhì)土層，土壤類型選用軟件內(nèi)置壤土，土壤水力參數(shù)見圖1。模擬土柱高度為100 cm，采用等間隔剖分方式將包氣帶分成100份。模擬時(shí)長(zhǎng)共計(jì)24 h，模擬開始時(shí)設(shè)置2 h降雨條件。

圖1 模擬土柱示意圖

2.2.2 初始及邊界條件

模型通過(guò)循環(huán)迭代計(jì)算獲取穩(wěn)定剖面含水率分布作為初始條件。下邊界選擇自由排水邊界。在短歷時(shí)暴雨雨型設(shè)計(jì)方面，芝加哥雨型在國(guó)內(nèi)外應(yīng)用廣泛[19]。利用北京市II區(qū)新暴雨公式推求北京市2 a一遇設(shè)計(jì)暴雨量，降雨歷時(shí)選用120 min，根據(jù)王彬雁等[20]的研究，北京市80%以上的雨型為單峰型，且在北京市單峰雨型中，雨峰位于前中后的比例差別不大。故本研究只針對(duì)單峰型雨型進(jìn)行模擬，給定雨峰相對(duì)位置r(0

(5)

式中：q為設(shè)計(jì)暴雨強(qiáng)度(L/(s·hm2))；t為降雨歷時(shí)，取120 min;P為重現(xiàn)期，本研究以2 a一遇設(shè)計(jì)暴雨為例。北京市II區(qū)的P和t滿足2 a≤P≤100 a，5

北京市II 區(qū)2 a一遇芝加哥雨型推求的暴雨過(guò)程如圖2所示。

圖2 北京市II 區(qū)2 a一遇芝加哥雨型過(guò)程

2.2.3 擾動(dòng)參數(shù)設(shè)置

由模型原理可知，HYDRUS-1D模型的土壤水力參數(shù)主要有殘余含水量θr、飽和含水量θs、飽和導(dǎo)水率Ks、進(jìn)氣吸力倒數(shù)α和孔隙尺寸分布指數(shù)n，同時(shí)這5個(gè)參數(shù)也是應(yīng)用模型時(shí)需要輸入的主要的土壤參數(shù)。采用HYDRUS-1D模型中壤土的設(shè)定值作為標(biāo)準(zhǔn)參數(shù)，各參數(shù)均以10%為步長(zhǎng)，共分5種情況進(jìn)行上下擾動(dòng)，各情況對(duì)應(yīng)的參數(shù)設(shè)置見表1，根據(jù)前人相關(guān)研究可知各參數(shù)20%的增減變化幅度處于合理范圍內(nèi)。

表1 土壤水力參數(shù)擾動(dòng)結(jié)果

2.2.4 輸出變量

分別選擇HYDRUS-1D模型中t=120 min時(shí)的累積入滲量、地面10 cm層的土壤水分動(dòng)態(tài)分布及濕潤(rùn)鋒運(yùn)移距離為研究對(duì)象，分析參數(shù)擾動(dòng)對(duì)降雨結(jié)束后累積入滲量、地面10 cm深度的土壤水分動(dòng)態(tài)分布及濕潤(rùn)鋒運(yùn)移距離的影響。

3 土壤水力參數(shù)對(duì)計(jì)算結(jié)果的影響分析

3.1 分析方法

Morris[21]提出的篩選法是一種運(yùn)算簡(jiǎn)單且效率較高的敏感性分析方法，被廣泛應(yīng)用于各種水文模型中，本文采用修正的Morris篩選法對(duì)HYDRUS-1D模型參數(shù)進(jìn)行敏感性分析，參數(shù)靈敏度取多次擾動(dòng)所計(jì)算出的Morris系數(shù)的平均值，即

(6)

式中：S為參數(shù)靈敏度；N為模型運(yùn)行次數(shù)；Yi為第i次運(yùn)行后模型的輸出值；Yi+1為第i+1次運(yùn)行后模型的輸出值；Y0為模型計(jì)算結(jié)果的初始值；Pi為模型第i次運(yùn)行時(shí)參數(shù)取值相對(duì)于初始值變化的百分率；Pi+1為模型第i+1次運(yùn)行時(shí)參數(shù)取值相對(duì)于初始值變化的百分率;P0為參數(shù)初始取值。

參數(shù)敏感度分級(jí)如表2所示。

表2 參數(shù)敏感度分級(jí)

3.2 土壤水力參數(shù)對(duì)10 cm層土壤含水量動(dòng)態(tài)分布的影響

土壤含水量的變化與降雨過(guò)程相關(guān)程度較高，尤其是表層土壤含水量的變化[1]。本研究借助10 cm處土壤含水量的變化來(lái)反映表層土壤水分動(dòng)態(tài)變化。從圖3可以看出，表層土壤含水量隨降雨時(shí)間的增加迅速增加，降雨結(jié)束時(shí)達(dá)到最大。隨后，表層水分向深層遷移，表層含水量逐漸減少。各土壤水力參數(shù)對(duì)水分特征曲線及非飽和導(dǎo)水率影響程度不同是導(dǎo)致參數(shù)敏感性差異的主要原因。

圖3 10 cm層土壤含水量動(dòng)態(tài)分布

圖4 參數(shù)擾動(dòng)對(duì)10 cm層土壤含水量分布的影響程度

從圖4可知，α、Ks、n擾動(dòng)的敏感性均在雨峰后(t=50 min)達(dá)到最大，θr在t=40 min時(shí)達(dá)到最大，降雨結(jié)束后其敏感性均降低。α、Ks在降雨結(jié)束后，其|S|值均<0.05，敏感性級(jí)別為IV;θr在降雨結(jié)束后其敏感性也顯著降低，|S|值略>0.05。這說(shuō)明降雨結(jié)束后θr、α、Ks擾動(dòng)對(duì)10 cm層土壤水分動(dòng)態(tài)分布基本無(wú)太大影響。而θs擾動(dòng)在降雨開始后達(dá)到一個(gè)次峰值(|S|值為0.83)，經(jīng)歷一個(gè)短暫的下降后增加至峰值(|S|值為0.95)并保持微弱的浮動(dòng),說(shuō)明θs擾動(dòng)在整個(gè)10 cm層土壤水分動(dòng)態(tài)分布過(guò)程中，均具有較大影響。參數(shù)n的敏感性在降雨50 min時(shí)達(dá)到最大，其|S|值為1.37，敏感性級(jí)別為Ⅰ，隨后其敏感性急速下降，并隨后續(xù)的入滲過(guò)程呈現(xiàn)一個(gè)較緩的增長(zhǎng)。5個(gè)參數(shù)對(duì)10 cm層土壤水分動(dòng)態(tài)分布的敏感性最值排序?yàn)閚>θs>α>Ks>θr，其敏感度級(jí)別分別為Ⅰ、Ⅱ、Ⅱ、Ⅱ、Ⅱ。此外，從圖3還可以看出，θr和θs擾動(dòng)與10 cm層土壤含水量動(dòng)態(tài)分布的變化呈正相關(guān)關(guān)系，而α、n、Ks擾動(dòng)與10 cm層土壤含水量動(dòng)態(tài)分布的變化呈負(fù)相關(guān)關(guān)系。

以上結(jié)果說(shuō)明參數(shù)的敏感性并不是一個(gè)恒定不變的，不同的研究時(shí)段下參數(shù)敏感性不同，根據(jù)研究需要確定參數(shù)敏感性變化區(qū)間對(duì)簡(jiǎn)化試驗(yàn)工作具有重要意義。如在模擬分析10 cm層土壤水分動(dòng)態(tài)分布規(guī)律時(shí)，應(yīng)盡可能確保n、α、θs及Ks參數(shù)的準(zhǔn)確性，在分析降雨前后10 cm層土壤水分動(dòng)態(tài)分布規(guī)律時(shí)，應(yīng)尤其注意θs參數(shù)的準(zhǔn)確性。因此，參數(shù)敏感性隨時(shí)間的變化也應(yīng)該作為模型模擬時(shí)需要考慮的因素之一。

3.3 土壤水力參數(shù)對(duì)累積入滲量的影響

根據(jù)表3可以看出，5個(gè)參數(shù)擾動(dòng)累積入滲量的敏感性排序?yàn)閚>Ks>θs>α>θr，其敏感性級(jí)別分別為Ⅰ、Ⅱ、Ⅱ、Ⅲ、Ⅲ。其中θr擾動(dòng)對(duì)累積入滲量的影響很小，其|S|值僅為0.07。而θs、n、Ks的擾動(dòng)對(duì)累積入滲量的影響很大，尤其是參數(shù)n，其|S|值為1.27。因此，在模擬分析累積入滲量時(shí)，應(yīng)尤其注意θs、n、Ks這3個(gè)參數(shù)的準(zhǔn)確性。

表3 參數(shù)擾動(dòng)對(duì)累積入滲量的影響

此外，還可看出α與θr擾動(dòng)與累積入滲量的變化呈負(fù)相關(guān)關(guān)系，θs、n、Ks擾動(dòng)與累積入滲量的變化呈正相關(guān)關(guān)系。此結(jié)論與前人定水頭研究結(jié)果不同，如范嚴(yán)偉等[14]認(rèn)為θs擾動(dòng)與累計(jì)入滲量變化呈負(fù)相關(guān)關(guān)系，而本研究發(fā)現(xiàn)為正相關(guān)關(guān)系，本研究結(jié)果與張海闊等[17]的結(jié)論一致。

3.4 土壤水力參數(shù)對(duì)濕潤(rùn)鋒運(yùn)移距離的影響

定量分析t=1 440 min時(shí)，各參數(shù)擾動(dòng)對(duì)濕潤(rùn)鋒運(yùn)移距離的影響，并計(jì)算相應(yīng)的靈敏度，如表4所示。

表4 參數(shù)擾動(dòng)對(duì)濕潤(rùn)鋒運(yùn)移距離的影響

從表4可知，5個(gè)參數(shù)對(duì)濕潤(rùn)鋒運(yùn)移距離的敏感性排序?yàn)棣?θs>Ks>n>θr，其敏感性級(jí)別分別為Ⅱ、Ⅱ、Ⅱ、Ⅱ、Ⅲ。θr擾動(dòng)對(duì)濕潤(rùn)鋒運(yùn)移距離的影響較小。而θs、α、n及Ks擾動(dòng)對(duì)濕潤(rùn)鋒運(yùn)移距離的影響均較大，其|S|值分別為0.65、0.78、0.41、0.50。因此在模擬分析濕潤(rùn)鋒運(yùn)移距離時(shí)應(yīng)盡可能確保θs、α、n及Ks這4個(gè)參數(shù)的準(zhǔn)確性。此外，θs、α、n的擾動(dòng)與濕潤(rùn)鋒運(yùn)移距離呈負(fù)相關(guān)關(guān)系，Ks及θr擾動(dòng)與濕潤(rùn)鋒運(yùn)移距離呈正相關(guān)關(guān)系。該結(jié)論與定水頭條件下的研究結(jié)果并不一致。如范嚴(yán)偉等[14]研究發(fā)現(xiàn)在定水頭輸入條件下n與濕潤(rùn)鋒運(yùn)移距離呈正相關(guān)，本研究為負(fù)相關(guān)。

從圖5可以看出，不同參數(shù)擾動(dòng)下的土壤濕潤(rùn)深度隨時(shí)間的變化規(guī)律類似。土壤的濕潤(rùn)深度隨降雨過(guò)程時(shí)間增加基本呈線性增加，隨后在土壤水分重力和毛管孔隙的作用下，土壤濕潤(rùn)深度隨濕潤(rùn)鋒向深層運(yùn)移而增加，增速較降雨期間減緩。

圖5 土壤濕潤(rùn)鋒運(yùn)移距離隨時(shí)間的變化

結(jié)合累積入滲量及濕潤(rùn)鋒運(yùn)移距離的研究結(jié)果，對(duì)比前人定水頭的研究結(jié)果可知，不同的邊界條件下模型輸出結(jié)果不同。即邊界條件也是影響模型模擬結(jié)果的重要因素之一，因此在模型模擬時(shí)應(yīng)根據(jù)實(shí)際情況選擇合適的邊界條件。

4 討 論

水文模型的使用可極大的提高工作效率，但一般而言，水文模型通常會(huì)涉及到很多參數(shù)的輸入，由于對(duì)水文模型機(jī)理的認(rèn)知不足或參數(shù)本身難以通過(guò)試驗(yàn)確定等原因，通常很難做到同時(shí)提高很多參數(shù)的精度，輸入?yún)?shù)的不確定性往往會(huì)對(duì)模擬結(jié)果產(chǎn)生較大影響。參數(shù)敏感性分析的目的就是確定輸入?yún)?shù)對(duì)結(jié)果的影響程度，從而可以降低參數(shù)率定過(guò)程中的計(jì)算消耗，提高模型的量化能力，降低參數(shù)維數(shù)等以提高模型使用效率。

在參數(shù)敏感性分析方面有大量的研究方法[22-23]，本文采用最常用的Morris篩選法進(jìn)行分析。在應(yīng)用HYDRUS-1D模型進(jìn)行模擬時(shí)，模型參數(shù)的準(zhǔn)確性及邊界條件會(huì)對(duì)輸出結(jié)果產(chǎn)生較大影響。因此對(duì)不同的受影響研究對(duì)象開展參數(shù)的敏感性研究有助于簡(jiǎn)化后續(xù)的參數(shù)率定工作，提高模型應(yīng)用效率及模擬精度。

關(guān)于HYDRUS-1D模型參數(shù)的敏感性，前人的研究多集中在定水頭輸入條件下[15-18]，本研究則以降雨入滲作為邊界條件，以10 cm層土壤水分動(dòng)態(tài)分布、累積入滲量和濕潤(rùn)鋒運(yùn)移距離為受影響研究對(duì)象，研究了其隨5個(gè)主要土壤水力參數(shù)變化的敏感性規(guī)律，研究結(jié)果豐富了降雨入滲條件下HYDRUS-1D模型參數(shù)敏感性結(jié)果，對(duì)于模型的使用效率及模擬結(jié)果準(zhǔn)確性的提高具有實(shí)踐指導(dǎo)意義。

模型參數(shù)的確定應(yīng)根據(jù)其物理意義具體問(wèn)題具體分析。θs可通過(guò)實(shí)驗(yàn)準(zhǔn)確測(cè)定，因此可將θs取為實(shí)測(cè)值，從而避免VG模型參數(shù)確定過(guò)程中θs的誤差。有研究表明飽和導(dǎo)水率Ks與土壤基本性質(zhì)之間存在一定關(guān)系，因此可通過(guò)土壤質(zhì)地或土壤粒徑分析等土壤理化參數(shù)，建立土壤轉(zhuǎn)換函數(shù)，估算或預(yù)測(cè)Ks的空間分布。α是土壤進(jìn)氣值的倒數(shù)，在實(shí)際應(yīng)用中可提前給出α的取值范圍，根據(jù)多次迭代擬合確定。n與土壤粒徑分布有關(guān)，n值的確定方法與α類似。

筆者建議，在應(yīng)用HYDRUS-1D模型進(jìn)行模擬前可針對(duì)擬輸出結(jié)果進(jìn)行參數(shù)敏感性分析，明確相應(yīng)參數(shù)的敏感性程度或高敏感性區(qū)間，進(jìn)而減少敏感性低的參數(shù)的試驗(yàn)確定工作，提高工作效率。本研究的土壤類型為軟件內(nèi)置壤土，土層為均質(zhì)土層，未考慮土壤大孔隙及層狀土體結(jié)構(gòu)對(duì)土壤水分運(yùn)移的影響。且模型上邊界條件為2 a一遇芝加哥雨型，與實(shí)際降雨過(guò)程可能有一定差異。后續(xù)研究者可結(jié)合實(shí)測(cè)土壤入滲數(shù)據(jù)及土層條件進(jìn)行深入研究，從實(shí)測(cè)數(shù)據(jù)與模擬數(shù)據(jù)的擬合角度去進(jìn)一步分析參數(shù)敏感性對(duì)于模型應(yīng)用的具體影響。

5 結(jié) 論

降雨入滲條件下，5個(gè)參數(shù)擾動(dòng)對(duì)HYDRUS-1D模型輸出結(jié)果的敏感性存在如下規(guī)律：

(1)模型參數(shù)的敏感性并非一成不變，針對(duì)不同的輸出結(jié)果，參數(shù)敏感性隨時(shí)間變化規(guī)律不同。根據(jù)研究需要確定參數(shù)敏感性區(qū)間對(duì)簡(jiǎn)化試驗(yàn)工作具有重要意義。

(2)邊界條件也是影響模型參數(shù)敏感性的重要因素之一，不同的邊界條件會(huì)對(duì)輸出結(jié)果產(chǎn)生較大影響，因此在模型模擬時(shí)應(yīng)力求邊界條件的真實(shí)性。

(3)對(duì)于不同的受影響研究對(duì)象，5個(gè)參數(shù)的敏感性不同。如n、α、θs及Ks這4個(gè)參數(shù)對(duì)10 cm層水分動(dòng)態(tài)分布影響較大，θs、n、Ks對(duì)累積入滲量影響顯著。θs、α、n及Ks對(duì)濕潤(rùn)鋒運(yùn)移距離的影響較大。模擬時(shí)應(yīng)確保高敏感性參數(shù)的準(zhǔn)確性。