Van Genuchten模型參數變化對土壤入滲特性的影響分析

范嚴偉，趙文舉，畢貴權(蘭州理工大學能源與動力工程學院，蘭州 730050)

非飽和土壤導水率、擴散率以及土壤水分特征曲線等是數值模擬非飽和土壤水分運動必不可少的重要參數，準確獲取這些參數是定量研究土壤水分和溶質運動的基礎[1]。多年來，國內外學者建立了許多經驗模型來描述土壤水分特征曲線和非飽和導水率曲線，其中，Van Genuchten模型應用最為廣泛[2]。基于Richards方程開發的HYDRUS-1D/2D/3D軟件[3]，均采用VG模型擬合。為此，許多土壤物理學家致力于確定VG模型中的參數。邵明安等[4]提出推求非飽和土壤水分運動參數的簡單入滲法，用以確定VG模型中的參數α和n；李春友[5]、馬英杰[6]、陳大春[7]、廖林仙[8]、郭向紅[9]等分別采用單純形調優法、阻尼最小二乘法、隨機粒子群算法、智能算法、混合遺傳算法擬合土壤水力參數；魏義長等[10]運用Matlab軟件的非線性擬合和非線性回歸函數獲得VG模型參數值；王薇等[11]利用RETC軟件依據土壤的顆粒組成和容重轉換出VG模型參數值；查元源等[12]根據瞬時剖面法，推求VG模型參數Ks、α和n。

模型參數的準確程度是數值計算精度高低的關鍵，而模型中各參數對模擬結果的影響程度因其物理意義的不同而有所差異。因此，有必要了解VG模型參數變化對土壤入滲特性的影響規律，以有效識別關鍵參數，減少率定參數的數量，提高模型的運行效率。李法虎等[13]對土壤水分運動參數的靈敏性能進行了理論分析；王志濤等[14]分析了VG模型參數變化對粉壤土入滲特性的影響。目前，尚未開展VG模型參數變化對不同質地土壤入滲特性的影響。為此，基于HYDRUS-1D模型，對VG模型中的θr、θs、α、n和Ks參數進行單因素擾動，模擬獲得土壤濕潤鋒運移距離和累積入滲量，分析VG模型參數變化對入滲特性的影響程度，以便為優化和改進模型結構提供依據。

1 數學模型建立

1.1 基本方程

HYDRUS-1D模型對土壤水分運動的模擬采用Richards方程[15]。其表達式為：

(1)

式中：z為垂向坐標，規定z向上為正；θ為體積含水量，cm3/cm3；h在飽和區和非飽和區分別為壓力水頭和基質勢，cm；t為入滲時間，min；K(h)為非飽和導水率，cm/min。

Richards方程中涉及的土壤水分特征曲線θ(h)和非飽和土壤導水率K(h)采用VG模型擬合。即

(2)

K(h)=KsSle[1-(1-S1/me)m]2

(3)

式中：Se=(θ-θr)/(θs-θr)；θs為殘余含水量，cm3/cm3；θs為飽和含水量，cm3/cm3；α、n和m為土壤物理特性有關的擬合參數，α單位為cm-1；n>1；m=1-1/n；l=0.5；Ks為飽和導水率，cm/min。

1.2 定解條件

在HYDRUS-1D模型中，土壤質地、容重是通過VG模型中的參數體現的，壓力水頭和土壤初始含水量對累積入滲量和濕潤鋒運移距離影響較小。初始條件和邊界條件不會影響VG模型參數，故在模擬時以定水頭作為HYDRUS-1D模型的上邊界條件，上邊界水頭值取為2 cm，入滲過程中，濕潤鋒未到達下界面。土壤初始含水量根據土壤質地取為一恒值，其中粉壤土為0.12 cm3/cm3，壤土為0.10 cm3/cm3，砂壤土為0.08 cm3/cm3，砂土為0.06 cm3/cm3。其定解條件為：

θ=θ0,0≤z≤L,t=0

(4)

h=h0,z=0,t>0

(5)

θ=θ0,z=L,t>0

(6)

式中：θ0為土壤的初始含水量，cm3/cm3；h0為初始壓力水頭，cm。

1.3 模型參數

各類土壤及其Van Genuchten模型參數均采用HYDRUS-1D模型中典型土壤的設定值，結果見表1。模擬過程中，參數擾動幅度設置為5個水平，分別為-20%、-10%、0%、10%和20%。

表1 不同土質VG模型水力特性參數Tab.1 Hydraulic parameters in VG model of different soils

2 VG模型參數變化對入滲特性的影響分析

2.1 分析方法

為定量分析VG模型中各參數變化對土壤入滲特性的影響，采用靈敏度系數來表示參數變化對入滲過程的影響程度，即：

(7)

式中：RSi為靈敏度系數；ΔPi為某參數Pi的變化量；ΔQi為某參數Pi的變化而引起目標函數Qi的變化量。

2.2 參數θr變化對入滲特性的影響

定量分析參數θr對不同質地土壤濕潤鋒運移距離Z的影響，將不同擾動變幅下的θr值及T=100 min時的濕潤鋒運移距離值代入式(7)，得參數θr的靈敏度系數RSθr。如表2所示。

表2 參數θr對不同土質濕潤鋒運移距離的影響Tab.2 Influence of parameter θr to transport distance of wetting front for different soils

由表2可知：θr的擾動對濕潤鋒運移距離無影響。

同理可得出參數θr對不同擾動變幅下累積入滲量I的靈敏度。如表3所示。

表3 參數θr對不同土質累積入滲量的影響Tab.3 Influence of parameter θr to the cumulative infiltration for different soils

由表3可知：θr的擾動對土壤累積入滲量影響非常微弱。

綜合分析說明：土壤殘余含水量θr的變化對土壤入滲特性影響很小。考慮到土壤殘余含水量θr的獲取比較困難，因此，在沒有確定的或實驗測得的θr值的情況下，將穩定枯竭點上的基質勢(φ=-15 000 cm)所對應的含水量作為θr的估算值是可行的[16]。

2.3 參數θs變化對入滲特性的影響

模擬參數θs擾動下的濕潤鋒運移距離、累積入滲量與時間t的關系曲線圖，如圖1所示。

圖1 參數θs對不同土質入滲特性的影響Fig.1 Effect of parameter θs on infiltration characteristics for different soils

定量分析參數θs對不同質地土壤濕潤鋒運移距離的影響，將不同擾動變幅下的θs值及T=100 min時的濕潤鋒運移距離值代入式(7)，得參數θs的靈敏度系數RSθs。如表4所示。

表4 參數θs對不同土質濕潤鋒運移距離的影響Tab.4 Influence of parameter θs on transport distance of wetting front for different soils

由圖1及表4可知：θs的擾動與垂直濕潤鋒運移距離成負相關，對其影響較大，且隨土質由細變粗而影響逐漸增大；θs的負擾動對濕潤鋒運移距離的影響程度明顯強于正擾動。

同理可得出參數θs對不同擾動變幅下累積入滲量I的靈敏度。如表5所示。

由圖1及表5可知：θs的擾動與土壤累積入滲量成正相關，對其影響較大，且隨土質由細變粗而影響逐漸減小。

綜合分析說明：土壤飽和含水量θs的變化對土壤入滲特性影響較大。但θs可通過實驗準確測定。因此，在求解VG模型參數時，可將θs取為實測值，避免VG模型參數確定過程中θs的誤差。

表5 參數θs對不同土質累積入滲量的影響Tab.5 Influence of parameter θs to the cumulative infiltration for different soils

2.4 參數α變化對入滲特性的影響

模擬參數α擾動下的濕潤鋒運移距離、累積入滲量與時間t的關系曲線圖，如圖2所示。

圖2 參數α對不同土質入滲特性的影響Fig.2 Effect of parameters α on infiltration characteristics for different soils

定量分析參數α對不同質地土壤濕潤鋒運移距離的影響，將不同擾動變幅下的α值及T=100 min時的濕潤鋒運移距離值代入式(7)，得參數α的靈敏度系數RSα。如表6所示。

表6 參數α對不同土質濕潤鋒運移距離的影響Tab.6 Influence of parameter α on transport distance of wetting front for different soils

由圖2及表6可知：參數α的擾動與濕潤鋒運移距離成負相關，對其影響較小，且隨土質由細變粗而影響逐漸減小。

同理可得出參數α對不同擾動變幅下累積入滲量I的靈敏度。如表7所示。

表7 參數α對不同土質累積入滲量I的影響Tab.7 Influence of parameter α to the cumulative infiltration for different soils

由圖2及表7可知：參數α的擾動與土壤累積入滲量成負相關，對其影響較小，且隨土質由細變粗而影響逐漸減小。

綜合分析說明：參數α的變化對土壤入滲特性影響較小，且隨土質由細變粗而影響逐漸減小。在實際應用時，VG模型參數α一般是通過多次迭代擬合確定，但擬合過程中往往造成參數的不唯一性。因此，在擬合求解VG模型參數時，可給出參數α的取值范圍，以解決參數α不唯一問題，或采用簡單入滲法來直接推求，推求過程中應保證吸滲率的準確性[4]。

2.5 參數n變化對入滲特性的影響

模擬參數n擾動下的濕潤鋒運移距離、累積入滲量與時間t的關系曲線圖，如圖3所示。

圖3 參數n對不同土質入滲特性的影響Fig.3 Effect of parameter n on infiltration characteristics for different soils

定量分析參數n對不同質地土壤濕潤鋒運移距離的影響，將不同擾動變幅下的n值及T=100 min時的濕潤鋒運移距離值代入式(7)，得參數n的靈敏度系數RSn。如表8所示。

由圖3及表8可知：參數n的擾動與濕潤鋒運移距離成正相關，對粉壤土、壤土影響較大，隨土質由細變粗而影響迅速減小，表現為對砂壤土、砂土影響較小。

表8 參數n對不同土質濕潤鋒運移距離的影響Tab.8 Influence of parameter n on transport distance of wetting front for different soils

同理可得出參數n對不同擾動變幅下累積入滲量I的靈敏度。如表9所示。

表9 參數n對不同土質累積入滲量的影響Tab.9 Influence of parameter s to the cumulative infiltration for different soils

由圖3及表9可知：參數 的擾動與累積入滲量成正相關，對粉壤土、壤土影響較大，隨土質由細變粗而影響迅速減小，表現為對砂壤土、砂土影響較小。

綜合分析說明：參數n的變化對細質土壤入滲特性影響較大，對粗質土壤入滲特性影響較小。參數n與α相同，在實際應用時，一般是通過多次迭代擬合來確定，但擬合過程中容易造成參數n的不唯一性。因此，在擬合求解VG模型參數時，需給出參數n的取值范圍，以解決參數n不唯一問題。或采用簡單入滲法來直接推求，推求過程中應保證吸滲率的準確性[4]。

2.6 參數Ks變化對入滲特性的影響

模擬參數Ks擾動下的濕潤鋒運移距離、累積入滲量與時間t的關系曲線圖，如圖4所示。

圖4 參數Ks對不同土質入滲特性的影響Fig.4 Effect of parameter Ks on infiltration characteristics for different soils

定量分析參數Ks對不同質地土壤濕潤鋒運移距離的影響，將不同擾動變幅下的Ks值及T=100 min時的濕潤鋒運移距離值代入式(7)，得參數Ks的靈敏度系數RSKs。如表10所示。

表10 參數Ks對不同土質濕潤鋒運移距離的影響Tab.10 Influence of parameter Ks on transport distance of wetting front for different soils

由圖4及表10可知：參數Ks的擾動與濕潤鋒運移距離成正相關，對其影響較大，且隨土壤質地由細變粗而影響逐漸變大。

同理可得出參數Ks對不同擾動變幅下累積入滲量I的靈敏度。如表11所示。

表11 參數Ks對不同土質累積入滲量的影響Tab.11 Influence of parameter Ks to the cumulative infiltration for different soils

由圖4及表11可知：參數Ks的擾動與累積入滲量成正相關，對其影響較大，且隨土壤質地由細變粗而影響逐漸變大。

綜合分析說明：飽和導水率Ks的變化對土壤入滲特性影響較大，且隨土壤質地由細變粗而影響逐漸變大。Ks可通過實驗測定，但空間變異較大，應考慮其尺度效應。或在求解VG模型參數時，同時測定土壤水分特征曲線θ(h)和非飽和土壤水擴散率D(θ)，采用同步擬合方法，得出VG模型參數θr、θs、α、n和Ks值，以消除單個參數引起的總體誤差。

3 結 語

基于HYDRUS-1D模型，采用單因素影響法，分析VG模型參數變化對土壤入滲特性的影響，結果表明：

(1)參數θr對濕潤鋒運移距離和累積入滲量影響很小，基本無影響。

(2)參數θs對濕潤鋒運移距離影響較大，成負相關關系，隨土質由細變粗而影響增大，負擾動強于正擾動；同時，θs與累積入滲量成正相關關系，對粉壤土、壤土影響較大，對砂壤土、砂土影響較小。

(3)參數α對濕潤鋒運移距離和累積入滲量影響都較小，成負相關關系，隨土質由細變粗而影響減小。

(4)參數n與濕潤鋒運移距離和累積入滲量成正相關，對粉壤土、壤土影響較大，對砂壤土、砂土影響較小。

(5)參數Ks對濕潤鋒運移距離和累積入滲量影響較大，成正相關關系，隨土質由細變粗而影響增大。

□

[1] 邵明安,王全九,黃明斌. 土壤物理學[M]. 北京:高等教育出版社, 2006.

[2] Van Genuchten M Th. A closed form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Sci.Soc.Am.J.,1980,44:892-898.

[3] Simunek J, Van Genuchten M Th, Sejna M. HYDRUS: Model Use, Calibration, and Validation[J]. Transactions of the ASABE,2012,55(4):1 261-1 274.

[4] 邵明安,王全九,Robert Horton. 推求土壤水分運動參數的簡單入滲法I.理論分析[J]. 土壤學報,2000,37(1):1-6.

[5] 李春友,任 理,李保國. 利用優化方法求算Van Genuehten方程參數[J]. 水科學進展,2001,12(4):473-478.

[6] 馬英杰,虎膽·吐馬爾拜,沈冰. 利用阻尼最小二乘法求解Van Genuchten方程參數[J]. 農業工程學報, 2005,21(8):179-181.

[7] 陳大春,馬英杰. 基于隨機粒子群算法的Van Genuchten方程參數優化求解[J]. 農業工程學報,2006,22(12):82-85.

[8] 廖林仙,邵孝侯,徐俊增. 基于智能算法推求Van Genuchten方程的參數[J]. 水利學報,2007,38(S1):696-700.

[9] 郭向紅,孫西歡,馬娟娟. 基于混合遺傳算法估計Van Genuchten方程參數[J]. 水科學進展,2009,20(5):677-682.

[10] 魏義長,劉作新,康玲玲,等. 土壤持水曲線Van Genuchten模型求參的Matlab實現[J]. 土壤學報,2004,41(3):380-386.

[11] 王 薇,孟 杰,虎膽·吐馬爾白. RETC推求土壤水動力學參數的室內試驗研究[J]. 河北農業大學學報2008,31(1):99-102,106.

[12] 查元源,周發超,楊金忠. 一種由土壤剖面含水率估算土壤水力參數的方法[J]. 水利學報,2011,,42(8):883-891.

[13] 李法虎,傅建平,孫雪峰. 土壤水分運動參數的確定及其靈敏性能分析[J].灌溉排水學報,1993,12(3):6-14.

[14] 王志濤,繳錫云,韓紅亮,等. 土壤垂直一維入滲對VG模型參數的敏感性分析[J]. 河海大學學報(自然科學版),2013,41(1):80-84.

[15] Simunek J, Van Genuchten M T, Sejna M. The HYDRUS-1D software package for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media[J]. University of California-Riverside Research Reports, 2005,3:1-240.

[16] 王金生,楊志峰, 陳家軍,等. 包氣帶土壤水分滯留特征研究[J]. 水利學報, 2000,31(2):1-6.