李欣麗,陳 力
(西安理工大學西北旱區生態水利國家重點實驗室,陜西 西安 710048)
非飽和土壤水是連接地表水和地下水的紐帶,也是農作物和自然植被所需水分的主要來源,研究土壤水分的運動規律可以揭示水文循環的本質,對提高農作物產量、保護生態環境等具有重要意義[1]。目前,有關土壤水分運動的研究手段主要以物理實驗和數值模擬為主。隨著計算機技術自20世紀80年代起開始普及,數值模擬方法因其操作簡便、效率高、受限小、耗時短、適用性廣泛等特點被廣泛應用于包氣帶水分運動的研究中,逐漸成為定量分析土壤水分入滲特性的重要手段[2-4]。近年來,國內外學者研發了許多適用于多孔介質中水分運動的數值模擬軟件,其中由美國農業部、美國鹽堿實驗室等機構開發的HYDRUS-1D軟件是模擬非飽和土壤水分運動的有力工具之一[5]。HYDRUS-1D模型中采用van Genuchten模型描述土壤水分特征曲線[6-7],其適用范圍廣且擬合效果較好,但其參數較多且形式較為復雜,并且不同的參數因其物理意義不同對模擬結果有著不同程度的影響,模型參數的準確性會直接影響著模型的精度,因此,模型參數的確定至關重要,而對模型參數進行敏感性分析則有助于參數的率定,從而有利于提高模型的模擬精度。
同時,外部輸入條件也是影響模型模擬精度的一個重要因素。目前,有關VG模型參數的敏感性分析研究多集中于定水頭入滲條件,如范嚴偉等[8]研究了van Genuchten模型各參數擾動對不同質地土壤入滲特性的影響;王志濤等[9]研究了定水頭入滲條件下van Genuchten模型參數變化對粉壤土入滲特性的影響;高志鵬等[10]分析了HYDRUS模型中土壤水力參數變化對水流和溶質運移變化的敏感性;霍思遠等[11]研究發現對于不同的模型輸出變量有著不同的參數敏感性。然而,變水頭入滲條件更符合實際入滲情況,有關該條件下VG模型參數的敏感性分析較為少見。因此,本研究采用HYDRUS-1D模型模擬土壤一維入滲,以成都市兩年一遇設計暴雨過程作為模型上邊界條件輸入,采用局部分析的方法,對單一模型參數進行敏感性分析,目的在于找出影響土壤入滲特性較為明顯的VG模型參數,明確參數敏感性的變化規律,以便于更精準更科學的模擬非飽和土壤水分運動,提高工作效率以及模擬結果的準確性,從而更好地服務實際應用。
HYDRUS-1D軟件采用Richards方程描述一維垂向非飽和土壤中的水分運動[5],其表達式為:
(1)
式中θ——土壤體積含水率,cm3/cm3;t——時間,d;z——垂向坐標,cm;h——土壤水勢,cm;K(h)——土壤非飽和導水率,cm/d。
本文采用van Genuchten模型對土壤水分特征曲線和土壤非飽和導水率進行擬合,其方程為[12]:
(2)
(3)
其中,
(4)

模型模擬0~100 cm深度范圍內的均質土壤,土壤類型選用HYDRUS-1D軟件中內置砂壤土、壤土以及黏壤土,土壤水流模型選用VG模型,不考慮滯后效應,模擬時長為24 h,上邊界為大氣降水,下邊界為自由排水。芝加哥雨型在國內外短歷時暴雨雨型的推求中應用較為廣泛[13],其雨型多為單峰雨型。根據趙劉偉等[14]的研究,成都市短歷時降雨以單峰型為主,故本文使用成都市新暴雨公式推求成都市區2年一遇設計暴雨量,利用芝加哥雨型推求設計暴雨過程,并作為模型上邊界條件進行輸入,選取2 h作為降雨歷時,給定雨峰相對位置r=0.33(0 (5) 式中i——降雨強度,mm/min;t——降雨歷時,min;P——重現期,a。 利用芝加哥雨型所推求的成都市區兩年一遇設計暴雨過程見圖1。 圖1 成都市區兩年一遇暴雨過程 選取模型輸出變量中的壓力水頭、濕潤鋒運移距離、累積入滲量作為研究對象,分析VG模型各參數對輸出變量的影響。 局部分析法因其操作簡單且高效被廣泛應用于模型參數敏感性分析中,故本文采用局部分析的方法對VG模型各個參數進行敏感性分析[15]。模擬時,將其中一個參數進行擾動而不改變其他參數值,重新運行模型,依次得到各情景下的模型輸出變量值進行敏感性分析。 本文采用相對敏感值對參數進行歸一化處理,計算敏感性指數I[16],進行參數之間的敏感性對比。 (6) 式中O——模型輸出結果;Fi——影響模型輸出結果的參數值;ΔO——模型輸出結果的改變量;ΔFi——參數的改變量。 根據敏感性指數值,對敏感性進行分類,具體見表1。 表1 敏感性分類 VG模型主要包含殘余含水量(θr)、飽和導水率(Ks)、飽和含水量(θs)、與進氣吸力相關的參數(α)、孔隙尺寸分布指數(n)5個土壤水力參數,本文以HYDRUS-1D軟件中提供的砂壤土、壤土、黏壤土的經驗參數值為基準,將各參數分別上下擾動10%、20%,得到各情景下的參數取值,參數基準值見表2。根據相關研究[17-18]可知20%內的參數增減變化處于合理范圍內。 表2 不同質地土壤VG模型水力參數 圖2為Ks擾動下不同質地土壤3個不同時刻(120、360、1 440 min)壓力水頭隨深度的變化規律。從圖2可以看出,Ks的擾動對壓力水頭的影響較大,通過觀察同一時刻壓力水頭曲線之間所圍面積,面積越大,其影響越大。隨著時間的推移,參數Ks的影響區域逐漸變深。從圖3可以看出,Ks擾動對壓力水頭的影響隨土壤質地變細而減小;在降雨過程中,Ks的擾動與壤土、黏壤土壓力水頭呈正相關,即敏感性指數為負值(由于計算采用負壓,故敏感性指數為負時呈正相關),壓力水頭隨Ks增大而增大,降雨結束后,Ks與土壤表層10 cm內壓力水頭呈負相關,而與其他深度呈正相關。該結論與定水頭入滲條件下研究結論并不相同,如高志鵬等[10]研究表明Ks的擾動對壓力水頭的影響呈非線性正相關,其原因可能是定水頭入滲條件下研究水分入滲貫穿整個模擬期且其敏感指數采用平均值,而本文中降雨歷時僅為2 h且其敏感性受降雨影響較大。同一質地土壤的Ks敏感性在降雨結束時達到最大,之后隨時間增加而減小,且不同時刻所影響的土壤深度范圍也不同,隨濕潤鋒運移而不斷變化,例如T=120 min時,Ks對砂壤土的影響范圍為0~40 cm,而T=360 min時,Ks對砂壤土的影響范圍為40~60 cm。Ks的敏感性隨土壤質地變細而減小,Ks的擾動對砂壤土的影響最大,其敏感性指數最高可達-103,并且影響范圍為整個土壤剖面。 a)砂壤土 a)砂壤土 Ks擾動下濕潤鋒運移距離隨時間的變化規律曲線見圖4。計算T=1 440 min時刻各質地土壤濕潤鋒運移距離在Ks不同擾動下的敏感指數值,見表3(敏感類別按照參數對輸出變量影響最大的敏感性指數絕對值進行排序)。 圖4 Ks擾動下不同質地土壤濕潤鋒運移距離的變化規律 表3 Ks擾動對不同質地土壤濕潤鋒運移距離的影響 由圖4、表3參數的敏感指數可知,隨著時間的推移,濕潤鋒逐漸向深處推進,其增加速度呈先快后慢的趨勢,且增速在降雨結束后逐漸減緩。砂壤土在Ks擾動下的濕潤鋒運移距離最遠,黏壤土次之,壤土最淺。Ks擾動對濕潤鋒運移距離影響較大,Ks的擾動與濕潤鋒運移距離呈正相關,且其敏感性隨土壤質地變細而減小。Ks的正負擾動對砂壤土的濕潤鋒運移距離影響程度基本相同,而對于壤土和黏壤土來說,Ks的正擾動影響程度略大于負擾動。Ks與土壤滲透性能有關,質地較粗的土壤Ks較大,滲透能力較強,故濕潤鋒運移距離遠,由于其孔隙較大且排列較為疏松,故Ks擾動對粗質地土壤的影響較大。 由表4可知,3種質地壤敏感類別都為III,且Ks的擾動與累積入滲量呈正相關關系,且Ks的敏感性隨土壤質地變細而減小。Ks的擾動對3個輸出變量都具有較大影響,且其影響程度均隨土壤質地變細而減小。Ks空間變異顯著,可由實驗測定,但應考慮其尺度效應。有相關研究[19]表明,Ks可根據土壤理化參數(土壤質地、粒徑分析等)建立土壤轉換函數,估測其空間分布。 表4 Ks擾動對不同質地土壤累積入滲量的影響 θr的擾動對壓力水頭的影響規律與Ks相似,但影響程度遠小于Ks。以壤土為例(圖5),降雨結束時,θr擾動的影響范圍為0~20 cm,并且其最大敏感指數值僅為-2.65(圖6),而同情況下Ks擾動的敏感指數可達14.12。與其他參數相比,θr的擾動對壓力水頭的影響較小,并且θr對壓力水頭的影響程度于同一時刻呈先增加后減少的趨勢。 a)砂壤土 a)砂壤土 由圖7、表5參數的敏感指數可知,在θr擾動下,土壤質地越粗,濕潤鋒運移距離越遠;θr的擾動對砂壤土的濕潤鋒入滲深度影響較小,且其影響隨土壤質地由粗至細逐漸增加。θr對VG模型低含水率段影響顯著,質地較細的土壤有著較高的θr,其持水能力強,孔隙較細且連通性較差,故θr擾動對細質地土壤的影響較大。 圖7 θr擾動下不同質地土壤濕潤鋒運移距離的變化規律 表5 θr擾動對不同質地土壤濕潤鋒運移距離的影響 由表6可知,3種質地土壤的敏感類別分別為I、II、III,且θr的擾動對砂壤土、壤土的累積入滲量影響甚微,其敏感性指數值均低于±0.1,可忽略不計;而對于黏壤土來說,在θr擾動為-20%時,也僅為-0.21,影響較小;θr的敏感性隨土壤質地變細而增加,但與其他參數相比,其影響很小。累積入滲量一般與初始含水量有關,而與θr關系較小。總的來看,θr的擾動對不同輸出變量的影響規律不同,但影響程度均較低。由于實際殘余含水量θr的獲取較為困難,在θr未確定時,可將壓頭為-15 000 cm時所對應的含水量作為θr的初估值[12]。 表6 θr擾動對不同質地土壤累積入滲量的影響 由圖8、9可知,θs的變化對壓力水頭的影響較大,其敏感性隨土壤質地變細而減小,并且θs的正擾動對壓力水頭的影響大于負擾動。以砂壤土為例,T=120 min時,+20%擾動情景下θs的最大敏感指數為106.5,是-20%擾動的2.13倍。降雨過程中,θs的擾動與壤土、黏壤土壓力水頭呈負相關,降雨結束后,θs與土壤表層20 cm內壓力水頭呈正相關。該結論與定水頭入滲條件下研究結論并不相同,如高志鵬等[10]研究表明θs的擾動對壓力水頭的影響呈非線性負相關,定水頭入滲會使得表層土壤在短時間內達到飽和,其水分入滲貫穿整個模擬期且敏感指數采用平均值,而本文中降雨歷時僅為2 h且其敏感性受降雨影響較大。 a)砂壤土 a)砂壤土 由圖10和表7參數的敏感指數可知,θs擾動對濕潤鋒運移距離的影響較大,θs的擾動與砂壤土、壤土濕潤鋒運移距離呈負相關,且其負擾動影響大于正擾動。濕潤鋒運移距離在θs擾動下隨土壤質地變細而減小。由表8可知,3種質地土壤的敏感類別分別為II、III、III,且θs的擾動與累積入滲量呈正相關,且其敏感性隨土壤質地由粗到細逐漸增加。對于砂壤土和壤土來說,θs的正負擾動對累積入滲量影響程度基本相同。θs越小,越易飽和,濕潤鋒運移距離越遠,累積入滲量越大,而粗質地土壤飽和導水率大,其孔隙較大且排列較為疏松,更易入滲,故θs在同擾動情況下對細質地土壤累積入滲量的影響更大。θs擾動對3個輸出變量的影響規律不同,但影響程度均較大。實際上,θs可通過實驗準確獲得,故應盡可能取實測值以減少誤差。 圖10 θs擾動下不同質地土壤濕潤鋒運移距離的變化規律 表7 θs擾動對不同質地土壤濕潤鋒運移距離的影響 表8 θs擾動對不同質地土壤累積入滲量的影響 由圖11、12可知,α對地表0~20 cm附近壓力水頭的敏感性較高,且其敏感性隨土壤質地變細而減小;α的影響區域隨時間推移而逐漸擴大,以壤土為例,T=120 min時,α的影響范圍為0~20 cm,而T=360 min時,其影響范圍增至0~40 cm。 a)砂壤土 a)砂壤土 由圖13及表9可知,α的負擾動越大,濕潤鋒的運移距離越深;α的擾動對壤土、黏壤土濕潤鋒入滲深度的影響較大,敏感類別為IV;α是對壤土濕潤鋒運移影響最大的參數。由表10可知,3種質地土壤的敏感類別分別為II、II、III,且α的變化對累積入滲量的影響較小,其敏感指數范圍僅為-0.26~0.14,影響程度比θr略高。α一般被認為進氣吸力的倒數,粗質地土壤進氣吸力較小,更易排水,其α值更大,故α在相同擾動情況下對粗質地土壤濕潤鋒運移的影響較小。累積入滲量一般與初始含水量及飽和含水量有關,與α關聯較小。α的擾動對濕潤鋒運移距離影響較大,對壓力水頭和累積入滲量影響較小,但略大于θr。α是與進氣值和孔徑分布相關的參數,通常經由土壤水分特征曲線多次迭代擬合確定,有研究[20]表明,利用簡單入滲法直接推求VG模型參數α具有較高的測定精度,但需保證吸滲率的準確度。 圖13 α擾動下不同質地土壤濕潤鋒運移距離的變化規律 表9 α擾動對不同質地土壤濕潤鋒運移距離的影響 表10 α擾動對不同質地土壤累積入滲量的影響 由圖14、15可以看出,降雨結束時,n的擾動對0~10 cm土壤層壓力水頭的影響最大,以砂壤土為例,其最大敏感性指數值可達-350;并且n擾動的影響范圍隨土壤質地變細逐漸減小,n的擾動對壓力水頭的影響規律與θs相似,但其影響程度大于θs,其敏感性隨土壤質地變細而減小。以壤土為例,-10%擾動下n的最大敏感值為75.3,是同情況下θs的12.86倍。降雨結束之后,n的擾動對3種質地土壤壓力水頭的影響較小,均不超過±3.5。 a)砂壤土 a)砂壤土 由圖16及表11可知,n負擾動對濕潤鋒運移距離的影響大于正擾動,n擾動對細質地土壤影響較大。n擾動與壤土的濕潤鋒運移距離呈負相關,該結論與定水頭入滲條件下研究結論并不相同,如范嚴偉等[8]研究發現n擾動與濕潤鋒運移距離呈正相關關系,其原因可能是定水頭入滲條件下供水速率均勻,入滲取決于土壤滲透能力,而變水頭入滲條件供水速率不恒定,入滲取決于供水速率或者土壤滲透能力,且n不確定性強,易受土壤結構、質地等因素以及外部條件的影響。由表12可知,3種質地土壤的敏感類別分別為III、IV、IV,且n的擾動與累積入滲量呈正相關關系,其敏感性隨土壤質地變細而顯著增加,質地較細的土壤n較小,其結合水能力較強,水分不易排出,在n擾動相同的情況下,細質地土壤累積入滲量所受到的影響較大,n的正負擾動對累積入滲量的影響程度基本相同;相比于其他參數,n是對壤土、黏壤土累積入滲量影響最大的參數。總之,n擾動對3個輸出變量都具有較大影響,其確定與參數α類似。 圖16 n擾動下不同質地土壤濕潤鋒運移距離的變化規律 表11 n擾動對不同質地土壤濕潤鋒運移距離的影響 表12 n擾動對不同質地土壤累積入滲量的影響 結合前人研究[8-9,21]對比,邊界條件不同也會導致參數敏感性的改變。變水頭入滲條件與定水頭入滲對模型輸出結果的影響存在差異,定水頭入滲條件下的累積入滲量大于變水頭入滲,而濕潤鋒運移距離則小于后者。入滲水頭對入滲系數有著較大的影響[22],定水頭入滲條件下入滲水頭不變,入滲主要取決于土壤滲透條件,而變水頭入滲條件下入滲水頭隨時間變化,入滲取決于土壤滲透條件與供水條件兩者的相互作用。對于輸出變量濕潤鋒運移距離來說,定水頭入滲條件下參數θr、Ks的敏感性大小與變水頭入滲基本相同,參數θs、α敏感性小于變水頭入滲,而參數n大于變水頭入滲;對于累積入滲量來說,定水頭入滲條件下參數Ks的敏感性大小與變水頭入滲基本相同,參數θr、θs、α大于變水頭入滲,而參數n小于變水頭入滲。 土壤質地對不同輸出結果參數敏感性有著一定的影響[23-24],土壤的持水能力與孔隙大小及其分布狀況有關,細質地的土壤有著較高的θr和θs,持水能力較強,孔隙較細且連通透性較差,濕潤鋒運移距離較小,累積入滲量較小,故在θs擾動下,細質地土壤所受影響較大;粗質地土壤Ks較大,滲透能力強,在Ks擾動下,粗質地土壤所受影響較大;土壤質地與α、n存在某種聯系,土壤質地越粗,α、n越大,在相同的α、n擾動下,細質地土壤濕潤鋒運移距離的變化更大;累積入滲量一般與初始含水量及θs有關,與θr、α、n關聯較小,故參數敏感性也較小。 VG模型參數應依據自身物理意義進行確定。不同質地土壤條件下的θr數值較小,且變化范圍較小,具有較強的確定性,并且θr對輸出變量的影響均較小,故可采用PTFs法進行推測或者使用壓頭為-15 000 cm時所對應的含水量;θs的變化范圍較小,確定性較強,實驗測定的θs值較為準確,故實際應用時盡量選取實測值;α、n具有很強的不確定性,隨土壤質地和結構不同而變化,土壤持水能力隨α增加而減小,有研究[25]表明,α不僅與進氣吸力有關,與n也存在非線性關系,只有在土壤質地較粗時,1/α與進氣值近似相等,細質地土壤的1/α均大于進氣值;n的取值影響著土壤水分特征曲線的形狀,并與土壤孔隙大小分布有關,采用土壤水分特征曲線多次迭代擬合確定;Ks變化范圍較大,不確定性強,空間變異性較大,Ks可由試驗測定,但應根據所選研究對象選擇適合的方法,例如,采用變水頭法測定較細土壤的Ks較為合適,并可基于實測值結合土壤水分運動數據進一步優化擬合,以提高估計效率,減小不確定性。除此之外,也可根據土壤理化參數建立土壤轉換函數,估算Ks空間分布。 變水頭入滲條件下,不同質地土壤VG模型各參數對模型輸出結果的敏感性具有以下規律。 a)VG模型各參數對壓力水頭的影響均隨土壤質地變細而減小,并且影響區域隨時間增加逐漸變深。同質地土壤各參數敏感性在降雨結束時達到最大,之后隨時間增加而減小。各參數敏感性排序為:n>θs>Ks>α>θr(砂壤土、壤土);n>θs>Ks>θr>α(黏壤土)。其中,n、Ks、θs是對壓力水頭影響較大的參數。 b)VG模型各參數對不同質地土壤濕潤鋒運移距離的影響具有較大差異,各參數敏感性排序為:θs>Ks>α>n>θr(砂壤土);α>θs>Ks>θr>n(壤土);n>α>θr>θs>Ks(黏壤土)。 c)同一參數的正負擾動對累積入滲量的影響程度基本相同,各參數敏感性排序為:Ks>n>θs>α>θr(砂壤土);n>Ks>θs>α>θr(壤土);n>θs>Ks>α>θr(黏壤土)。其中,n、Ks、θs是對累積入滲量影響較大的參數,并且其擾動與累積入滲量呈正相關關系。 d)邊界條件的不同也會導致參數敏感性發生改變,變水頭入滲下VG模型各參數對模型輸出結果的敏感性與定水頭入滲不同,故進行模擬時應確保邊界條件的真實性。
1.3 輸出變量選取
2 參數敏感性分析
2.1 參數敏感性分析理論依據

2.2 參數取值

3 結果與分析
3.1 Ks敏感性分析





3.2 θr敏感性分析





3.3 θs敏感性分析





3.4 α敏感性分析





3.5 n敏感性分析





4 討論
5 結論