van Genuchten模型参数的物理意义

陈卫金，程东会，2，陶 伟

(1.长安大学环境科学与工程学院，陕西 西安 710054；2.长安大学旱区地下水文与生态效应教育部重点实验室，陕西 西安 710054)

vanGenuchten模型参数的物理意义

陈卫金1，程东会1，2，陶 伟1

(1.长安大学环境科学与工程学院，陕西 西安 710054；2.长安大学旱区地下水文与生态效应教育部重点实验室，陕西 西安 710054)

van Genuchten模型是目前拟合水土特征曲线应用最广泛的数学模型，但该模型参数多，且物理意义复杂，特别是参数α的物理意义至今尚未有统一的认识。文章通过理论推导，结合4种样品(粗砂、中砂、细砂和粉砂)实测的水土特征曲线，对van Genuchten模型的参数α、m和n的物理意义进行了探讨，重点研究了α与进气值和拐点处负压(hi)的关系。推导出的van Genuchten模型的1/α与hi的关系和实测数据均证明了1/α不仅与hi有关，还受参数m和n的影响。粗砂和中砂的1/α与hi大致相等，细砂和粉砂的1/α则大于hi。粗砂的进气值与1/α近似相等，其它介质中1/α与进气值之间存在更复杂的关系。另外，m和n作为独立变量在拟合水土特征曲线时精度更高，但在计算非饱和渗透系数时给定m=1-1/n的约束条件是必要的。

水土特征曲线；van Genuchten模型；进气值；拐点；非饱和渗透系数

非饱和介质基质势(负压，或吸力)是随介质含水率而变化的，表征这种性质的曲线称为非饱和介质水土特征曲线(SWCC)。水土特征曲线表示非饱和介质水分的能量和数量之间的关系，是反映非饱和介质基本性质、持水性能和水分运动规律的曲线[1]。根据水土特征曲线可以获得非饱和介质的水力学和土力学性质，如土壤的水分分布规律、渗透系数、滞后特性、体积变化等[2～5]，被广泛应用于非饱和介质的降雨入渗和水分蒸发、溶质的扩散和运移、非饱和土体或边坡的稳定性等多个方面[6～7]。从20世纪50年代至今，已有多达几十个数学模型被提出来表征水土特征曲线，包括van Genuchten模型(VG)、Brooks and Corey模型(BC)、Gardner模型、Fredlund and Xing模型、Mckee and Bumb模型和Williams模型等[8～14]。由于这些数学模型是经验或半经验的，模型中很多参数的物理意义不明确，难以直接测量，参数的获取方法是对实测的负压-含水率数据进行拟合而得到。而明确模型中参数的物理意义，对验证测试数据的可靠性十分重要，也可以为测定参数直接获取水土特征曲线数学模型的方法奠定基础。

VG模型是最为常用的水土特征曲线数学模型。该模型能够表示全负压范围内的水土特征曲线，对不同类型介质有很高的拟合精度，具有广泛的适用性，其表达形式如式(1)。

式中:θ——体积含水率/(cm3·cm-3)；

h——负压(cmH2O)，取正值；

θs、θr——分别为饱和含水率和残余含水率/(cm3·cm-3)；

α、m、n——模型参数。

本文的研究目标是以4种有代表性介质的实测水土特征曲线为基础，通过理论推导，讨论VG模型中参数α、m和n的几何学意义和物理意义，重点讨论α与进气值ha和水土特征曲线拐点处负压hi的关系，并分析m和n之间的相互关系及其与BC模型中参数λ的关系，以及m和n取值对非饱和渗透系数的影响。

1 试验方法

试验所用样品为河流沉积物中4种典型的介质，粉砂、细砂、中砂和粗砂。粉砂数据为查阅资料获取[19]，其它3种介质的样品为本次研究中在渭河陕西咸阳段河漫滩采集的河流沉积物。细砂、中砂和粗砂3种介质的粒径分析结果如图1(粉砂缺少粒径分析资料)。

图1 试验样品的粒级分布特征Fig.1 Grain-size distribution for the samples

试验样品的水土特征曲线采用悬挂水柱法实测获得。悬挂水柱法的试验装置如图2。试验所用布氏漏斗为市面上购买的普通布氏漏斗，经测定漏斗的多孔板在吸力达到3～4 mH2O时穿透，即漏斗进气压力为3～4 mH2O，可以满足细砂、中砂和粗砂的测定。测定样品的水土特征曲线时，使布氏漏斗里的样品通过多孔板与U形管中的水柱建立水力联系，通过降低或升高U形管的水柱高度实现样品的排水或吸湿。U形管每次升降幅度约为5～20 cm，在容水率变化较大处升降幅度小，其它阶段升降幅度较大。样品含水率θ和相应的负压h(本文均取正值)通过压力平衡状态时U形管中水柱高度的变化ΔH获取。饱和样品通过多次降低U形管，获得多组含水率及相应负压的数据，测得样品实测的排干曲线；排干过程完成后，升高U形管中的水柱高度，则可获得实测的吸湿曲线[20]。该装置也可以测定从任意含水率开始排水或吸湿的扫描曲线。

图2 悬挂水柱法试验装置示意图Fig.2 Schematic diagram of the hanging water column method

2 试验结果

2.1实测的水土特征曲线

通过悬挂水柱法获得4种样品的水土特征曲线如图3。由于4个样品粒径的差异，水土特征曲线表现出的形状有明显差异。细砂和粉砂的曲线近似为“S”形，而粗砂和中砂则近似“J”形。

图3 试验样品的实测水土特征曲线Fig.3 Measured soil-water characteristics curves for the samples

2.2水土特征曲线的数学模型

利用VG模型对4个样品的实测水土特征曲线进行拟合求得模型参数，从而得到样品水土特征曲线的数学模型。最常用的VG模型为三参数模型，分别为α、n、m。当设定m=1-1/n的限制条件时，则变为二参数模型，2个参数分别为α和n。一般使用VG模型时，θr作为已知参数，其数值可以通过试验获取。但从本质上，随着吸力不断增大，θr的值不断小幅度地减小，并不存在固定的θr值，因此在一些涉及介质干燥过程的研究时，将θr也作为待求参数，此时VG模型为4参数模型，参数包括α、n、m和θr。

本文分别以二参数、三参数和四参数的VG模型对实测数据进行拟合，拟合结果如图4，获得的模型参数如表1。4种介质的VG模型参数表明，细砂的三参数模型图形上在θr附近略有偏离实测数据，参数α、n与四参数模型也有较大误差，粉砂、中砂、粗砂的三参数模型则均与四参数模型偏差不大；给定θr和m=1-1/n后，粗砂二参数拟合的n与四参数相比相差较大，而中砂、细砂和粉砂的二参数拟合的α、n与四参数均相差不大。因此，三参数模型对细砂拟合水土特征曲线误差较大，更适合于中砂、粗砂及更细粒的粉砂，而二参数模型对粗砂拟合水土特征曲线有较大误差，而对中砂、细砂和粉砂误差相对较小。总体上四参数模型能同时满足3种介质实测数据的拟合，拟合效果最好。

图4 试验样品不同参数VG模型的拟合曲线Fig.4 Fitting curves using the VG model with different parameter numbers for the samples(a)粉砂； (b)细砂； (c)中砂； (d)粗砂

3 讨论

3.1VG模型中的参数α

3.1.1参数α的几何学特征及其与进气值ha的关系

VG模型中参数α的值会影响拟合水土特征曲线的几何形态。从曲线形状而言，当其它参数不变时，α值越小，水土特征曲线中间平缓部分越高，即相同含水率所对应的负压越大(图5)。从物理含义而言，α值越小，介质的持水能力越强。细粒介质的α值小于粗粒介质的α值(表1)。

表1 试验样品不同参数VG模型的拟合参数值

注：四参数为θr、α、n、m；三参数为α、n、m；二参数为α、n；R2为决定系数

图5 粗砂介质的VG模型中当α取不同值时的拟合曲线(n=9.95，m=0.42)Fig.5 Fitting curves of coarse sand using the VG model with different values of parameter α when n=9.95 and m=0.42

许多观点均把VG模型中的α与进气值ha相联系，甚至认为α值的倒数大致与ha相等。由于VG模型中，没有直接表征进气值的参数，而在Brooks and Corey模型(BC模型)中则直接定义了ha，因此ha通常利用BC模型拟合实测数据得到，BC模型的表达式为：[9]。

式中：ha——进气值/cmH2O；

λ——表征空隙大小分布特征的参数。

在以往几个研究中通过大量样品的统计分析或通过理论推导给出了几个α值与ha之间定量关系的经验表达式。Tinjum等利用24个黏土样品得到一个适用于细粒介质α值与ha关系的表达式，式(3)[22]：

另外Lenhard定义了一个有广泛适用性的α值与ha的关系，如式(4)[18]：

其中，Sx由一个经验的表达式确定，如式(5)：

理论上，VG模型中α值与ha之间的关系更复杂。Morel-Seytoux把VG模型和BC模型联系起来，定量地考察VG模型中α值与ha之间的这种复杂的关系，如式(6)[23]：

上述3个表达式也表明1/α值与ha没有直接相等关系，只有在某些特殊条件下二者才可能相近。为了清楚满足1/α值与ha相近的条件，本文将表1中二参数VG模型的拟合参数分别代入式(3)、式(4)和式(6)，可计算得各方程的1/α值，计算结果如表2。

表2 试验样品水土特征曲线在BC模型的进气值以及在VG模型的1/α

表2表明，从3个经验公式计算的1/α值均与ha存在一定的误差，但介质粒径越大，1/α与ha越接近，而细粒介质的1/α值与ha之间则相差较大，甚至达到几倍的误差。式(3)对4种介质的1/α计算结果误差很大，因为它是由黏土样品中得出的，所以式(3)只适用于极细粒的介质，对砂类介质不适用。式(4)和式(6)中，粗粒介质计算的1/α与拟合值接近，细粒介质计算的1/α值较拟合值偏大，因此式(4)和式(6)对中砂、粗砂适用性良好，而对粉砂和细砂误差较大。总之不论从理论上还是经验数据分析，1/α值与ha只有在粗粒介质中才可能近似相等。

3.1.2参数α与拐点处负压hi的关系

近年来，许多研究中将VG模型中的参数α与拐点处负压hi相联系，认为参数α的倒数与hi相近。事实上，α与hi的关系从理论上可以得到证明。在VG模型数学表达式中，如果求含水率的二阶导数，可得到式(7)：

由于拐点处含水率的二阶导数为零，将式(7)左端为零后化简，可得hi与α的表达式，如式(8)：

如果利用m=1-1/n的限制条件，则可得到hi与α的更简明的表达式，如式(9)：

式(9)表明，1/α与hi的关系还与参数m和n有关。更大的n值，或m接近于1时，才会使1/α更接近于hi。而较粗粒介质更易满足这些条件。本次研究将表1中拟合的相关参数值代入式(9)，可得不同介质的hi(表3)。表3表明式(9)计算得到的hi与中砂、粗砂拟合的1/α值均大致相等。细砂误差约为5 cmH2O，粉砂误差最大，大于40 cmH2O。

表3 样品水土特征曲线中hi与1/α的比较

3.2VG模型中的参数m和n

3.2.1参数m和n的几何学特征

以粗砂排干过程为例，讨论参数m和n取值对拟合水土特征曲线几何学特征的影响(图6)。图6表明，参数n的取值影响拟合水土特征曲线的弯曲程度，n越小，拟合的水土特征曲线越平缓，且曲线绕着一个固定轴点随n变化。在轴点右侧较低“膝盖”处，n越小，“膝盖”越低，而轴点左侧相反。此外，n的取值还与介质孔隙大小分布有关，n越小孔隙粒径分布范围越广，n越大孔隙粒径分布越均匀(图1)。参数m的取值影响拟合水土特征曲线左侧拐弯处的弯曲程度，m越小，曲线越平缓，左侧拐弯弯曲程度越小，拐弯越靠右，且除接近θs部分外，m越小，相同含水率对应的负压越大。

图6 参数m、n变化的VG模型拟合曲线Fig.6 Fitting curves using the VG model with variable parameters m and n(a)α=0.058,m=0.40; (b)α=0.058, n=9.95

3.2.2参数m和n之间的关系

三参数和四参数的VG模型中，参数m和n均为独立变量，二参数模型中设定了m=1-1/n的限制条件。从样品拟合结果可以看出(表2和图4)，给定m=1-1/n后，中、细粒介质拟合得到的参数α和n的值均与四参数的拟合值接近，误差很小；粗粒介质拟合参数α的值与四参数的拟合值接近，但n值的差距较大。因此可以初步判断m=1-1/n的限制条件更适合粒径较细的介质，而对粗粒介质，该限制条件会使水土特征曲线的数学模型出现明显的误差，应将m和n作为独立变量来拟合求参。

尽管VG模型中的参数m和n作为独立变量可以更准确的表达水土特征曲线，但是当涉及到非饱和渗透系数时，建立参数m和n之间的关系却是必要的。van Genuchten以VG模型为基础，设定m=1-1/n的限制条件，得到了非饱和渗透系数的表达式，如式(10)。

式中：K——非饱和渗透系数/(m·d-1)；

Ks——饱和渗透系数/(m·d-1)。

当m和n自由变化时，m的取值对K的影响显著。以细砂为例，给定n值而取不同的m值，分析m取值对求K的影响(图7)。当mlt;1-1/n时，K-h曲线几何形态为先下降再上升，且上升部分随着m取值越接近m=1-1/n越靠近x轴(图7a)；当m≥1-1/n时，K-h曲线均呈“L”形下降趋势，且m值越接近m=1-1/n平缓部分越靠近x轴(图7b)；当m=1-1/n时，曲线平缓部分最终趋向于x轴。因此，当m值自由变化时，求得的K值在高负压段可能会有较大误差，在给定限制条件m=1-1/n时求得的K才比较准确。

图7 参数m不同取值对非饱和渗透系数的影响(m=1-1/n时，n=2.929，m=0.658 6；α=0.019)Fig.7 Effect of the values of m on unsaturated hydraulic conductivity when n=2.929, m=0.658 6 and α=0.019 with m=1-1/n(a) mlt;1-1/n; (b)m≥1-1/n

3.2.3VG模型m与BC模型λ的关系

Lenhard 通过容水率dSe/dh将BC模型和VG模型统一，得到了BC模型中λ和VG模型中的m(m=1-1/n)的关系，如式(11)[18]：

将拟合的m值代入式(11)，结果如表4。由表4比较可知，Se取0.5时，除粉砂误差较大外，其余3种介质计算的λ值与拟合λ值均大致相等。说明式(11)对粉砂不适用。

表4 拟合的m、λ值与计算的λ值

4 结论

(1)van Genuchten模型中参数α不仅与拐点处负压有关，还受参数m和n的影响。更大的n值，或m接近于1时，才会使1/α更接近于hi。而粗粒介质更易满足这些条件。

(2)van Genuchten模型中参数α与进气值之间的关系较为复杂，也受其它模型参数的影响。只有介质较粗时，如粗砂或更粗粒的介质，1/α才与进气值近似相等，细粒介质中1/α值均大于进气值ha。

(3)虽然m和n作为独立变量在拟合水土特征曲线时精度更高，但在计算非饱和渗透系数时，给定m=1-1/n的约束条件是必要的，否则计算的渗透系数会在高负压段出现较大的误差。

[1] 雷志栋，杨诗秀，谢森传，等.土壤水动力学[M].北京：清华大学出版社，1988：18-24. [LEI Z D, YANG S X, XIE S C,etal. Soil water dynamics[M]. Beijing: Tsinghua University Press, 1988: 18-24. (in Chinese)]

[2] Kovacs G. Seepage Hydraulics[M]. Amsterdam: Elsevier Science Publishers, 1981.

[3] 刘晓东，施建勇.基于土水特征曲线预测城市固体废弃物(MSW)非饱和渗透系数研究[J].岩土工程学报，2012，34(5): 855-862. [LIU X D, SHI J Y. Unsaturated conductivity of MSW based on soil-water characteristic curve[J]. Chinese Journal of Geotechnical Engineering, 2012, 34(5): 855-862. (in Chinese)]

[4] 冯立, 张茂省, 孙萍萍, 等. 非饱和土脱湿与吸湿水力特性对比研究[J]. 水文地质工程地质, 2016, 43(2): 134-139. [FENG L, ZHANG M S, SUN P P,etal. A study of the hydro-mechanical properties of unsaturated loess under the drying and wetting path[J]. Hydrogeology amp; Engineering Geology, 2016, 43(2): 134-139. (in Chinese)]

[5] Fredlund D G, Rahardjo H. Soil Mechanics for Unsaturated Soils[M]. Hoboken:Wiley, 1993.

[6] 王勇. 压实非饱和黏土层生活源污染质运移及其效应[D]. 徐州：中国矿业大学， 2014. [WANG Y. Life source pollutant migration and its inducing effects in unsaturated compacted clay layer[D]. Xuzhou:China University of Mining and Technology, 2014. (in Chinese)]

[7] 尹平保，贺炜，张建仁，等. 洞庭湖大桥锚碇基础土围堰饱和-非饱和稳定性分析[J]. 水文地质工程地质，2016，43(2): 62-68. [YIN P B, HE W, ZHANG J R,etal. Stability of soil cofferdam for anchorage foundation of the Dongting Lake bridge considering unsaturated soil characteristics[J]. Hydrogeology amp; Engineering Geology, 2016, 43(2): 62-68. (in Chinese)]

[8] Genuchten M T V. A close-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(44): 892-898.

[9] Brooks R J, Corey A T. Hydraulic Properties of Porous Media[J]. Colorado State University Hydrology paper, 1964(3):22-27.

[10] Gardner W R. Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[J]. Soil Science, 1958, 85(4): 228-232.

[11] McKee C R, Bumb A C. The importance of unsaturated flow parameters in designing a monitoring system for hazardous wastes and environmental emergencies[C]// Proceedings of Hazardous Materials Control Research Institute National Conference. Houston:1984:50-58.

[12] Mckee C R, Bumb A C. Flow-Testing Coalbed Methane Production Wells in the Presence of Water and Gas[J]. Spe Formation Evaluation, 1987, 2(4): 599-608.

[13] Fredlund D G, Xing A. Equations for soil-water characteristic curve[J]. Canadian Geotechnical Journal, 1994, 31: 521-532.

[14] Williams J, Prebble R E, Williams W T,etal. The influence of texture, structure and clay mineralogy on the soil moisture characteristic[J]. Australian Journal of Soil Research, 1983, 21(1): 15-19.

[15] Genuchten M T V, Nielsen D R. On describing and predicting the hydraulic properties of unsaturated soil[J]. Annales Geophysicae, 1985, 3(5): 615-628.

[16] Benson C H, Chiang I, Chalermyanont T,etal. Estimating van Genuchten parametersαandnfor clean sands from particle size distribution data[J]. Geotechnical Special Publication, 2014(233): 410-427.

[17] Schaap, Marcel G, van Genuchten, Martinus Th. A modified mualem-van Genuchten formulation for improved description of the hydraulic conductivity near saturation[J]. Vadose Zone Journal, 2006, 5(1): 27-34.

[18] Lenhard R J, Parker J C, Mishra S. On the correspondence between Brooks-Corey and van Genuchten models[J]. Journal of Irrigation amp; Drainage Engineering, 1989, 115(4): 610-611.

[19] 赵雅琼. 非饱和带土壤水分特征曲线的测定与预测[D]. 西安：长安大学, 2015. [ZHAO Y Q. Determination and prediction of soil water characteristic curve in unsaturated zone[D]. Xi’an: Chang’an University, 2015. (in Chinese)]

[20] Sharma R S, Mohamed M H A. An experimental investigation of LNAPL migration in an unsaturated/saturated sand[J]. Engineering Geology, 2003, 70(3-4): 305-313.

[21] Milly, P C D. Estimation of Brooks-Corey parameters from water retention data[J]. Water Resources Research, 1987, 23(6): 1085-1089.

[22] Tinjum J M, Benson C H, Blotz L R. Soil-water characteristic curves for compacted clays”[J]. Journal of Geotechnical amp; Geoenvironmental Engineering, 1999,123(11):629-630.

[23] Morel-Seytoux H J, Meyer P D, Nachabe M,etal. Parameter equivalence for the Brooks-Corey and van Genuchten soil characteristics: preserving the effective capillary drive[J]. Water Resources Research, 1996, 32(5): 1251-1258.

责任编辑

：汪美华

PhysicalsignificanceoftheparametersinthevanGenuchtenmodel

CHEN Weijin1, CHENG Donghui1,2, TAO Wei1

(1.SchoolofEnvironmentalScienceandEngineering,Chang’anUniversity,Xi’an,Shaanxi710054,China; 2.KeyLaboratoryofSubsurfaceHydrologyandEcologicalEffectsinAridRegion(Chang’anUniversity),MinistryofEducation,Xi’an,Shaanxi710054,China)

The van Genuchten model is widely used in fitting soil-water characteristic curves. However, it has numerous parameters with the complex physical significance, and especially the parameterαhas not been well understood so far. This paper discusses the physical significance of the parametersα,mandnin the van Genuchten model on the basis of the results of both theoretical derivation and the measured values of the soil-water characteristic curves of the samples (silt sand, fine sand, medium sand, coarse sand), and focuses on the relationships between parameterαand air entry values or matrix potential at the inflection point (hi). The results illustrate that the relationship between the values of 1/αandhiis related to not onlyhibut also the parametersmandn, and the value of 1/αis approximately equal tohifor the medium sand and coarse sand, while it is greater thanhifor the fine sand and silt sand. The air entry value is equal to value of 1/αfor the coarse sand, and the complex relationships occur for the other media. In addition, there is a higher degree of accuracy when the water-soil characteristic curves are fitted asmandnare independent variables. It was necessary to give a constraint condition, i.e.,m=1-1/n, when the unsaturated hydraulic conductivity is calculated.

soil-water characteristic curve; van Genuchten model; air entry value; inflection point; unsaturated hydraulic conductivity

10.16030/j.cnki.issn.1000-3665.2017.06.22

S152.7

A

1000-3665(2017)06-0147-07

2017-01-20;

2017-03-05

国家自然科学基金资助项目(41472220)；中央高校基本科研业务费专项资金(310829162015)

陈卫金(1991-)，男，硕士研究生，主要从事地下水环境方面的研究。E-mail:578677316@qq.com

程东会(1969-)，男，教授，主要从事地下水资源与环境方面的教学和研究工作。E-mail: chdhbsh@chd.edu.cn