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.