基于VG模型确定土水特征曲线基本参数*

潘登丽 倪万魁 苑康泽 张镇飞 王熙俊

(长安大学地质工程与测绘学院， 西安 710054， 中国)

0 引 言

土水特征曲线(SWCC)反映了土中含水量(质量或者体积含水量)与吸力之间的关系(Sillers et al.，2001a，2001b)，可以由其估计非饱和土的其他工程特性，如渗透系数，抗剪强度和体变特性(Cook， 1991; Rao et al.，2011)。以往大量的研究已经提出了很多封闭形式的经验模型来拟合SWCC数据点，提出的每一个经验模型都能很好地拟合SWCC的减湿和增湿曲线(Fredlund， 2006)。其中，Genuchten(1980)提出的经验公式是最常用的模型之一，如式(1)所示。

(1)

式(1)为饱和度形式表示的VG模型。式中，a，m，n为拟合参数；ψ为基质吸力(kPa)；S为饱和度。为了计算方便，可取m=1-1/n，此时VG模型演变为Mualem(1976)模型，只含有两个参数，形式简单，应用方便，是岩土工程文献中常用的模型。

土水特征曲线的基本参数有进气值ψaev，反弯点斜率Ki、残余饱和度Sr和残余吸力值ψr等，准确确定土水特征曲线基本参数是非饱和土强度、渗透性以及本构关系研究的基础(高游等， 2017)。近年来针对SWCC的影响因素方面的研究比较普遍(贾宝新等， 2018； 石振明等， 2018； 李同录等， 2019)，但是如何准确地确定SWCC基本参数的研究还较少。

目前对SWCC残余饱和度Sr较为广泛认可的定义是气连通而水不连通时的饱和度，它在非饱和土水力、力学特性研究中有十分重要的位置。如在非饱和土强度预测方面(Bishop， 1959)，Bishop单应力状态变量强度公式中，需要求解有效应力参数χ(0≤χ≤1)，而Vanapalli证明该参数可用有效饱和度Se(Se=(S-Sr)/(1-Sr))来代替，其中，Sr为残余饱和度(Vanapalli et al.，1996)； 在非饱和土渗透系数的预测方面，残余饱和度Sr直接或间接应用在很多非饱和渗透系数模型中，如CCG模型(Childs et al.，1950)、Burdine模型(Burdine， 1953)和Mualem模型(Mualem， 1976)。

残余基质吸力ψr的众多定义均是属于经验的定义，没有明确的物理意义和理论支撑。如Genuchten(1980)把残余基质吸力为1500,kPa时土体的状态定义为残余状态。Sillers et al. (2001a，2001b)认为残余含水量对应的基质吸力为3000,kPa，同时他认为残余状态是土中孔隙水从受毛细作用转变到受吸附力作用的状态，这种定义方法目前比较流行。

在传统图解法中，反弯点主要由人为主观确定，在此点基础上画一条切线，分别与边界区、残余区土水特征曲线的直线相交，这3条直线两两相交，得到交点A和R，如图1所示，分别为进气点和残余点，从而得到SWCC的基本参数。这种方法的精确度依赖于画图人员如何确定反弯点I，如图2所示，不同的反弯点会产生不同的ψaev，这些ψaev值有很大差异性(ψaev=17， 27， 43,kPa)，说明传统画图法由于主观的不确定性和个人的判断会产生很大的误差(Soltani et al.，2017)。

图1 单峰土水特征曲线基本参数定义Fig. 1 Definitions of unimodal SWCC variables

图2 不同的反弯点会产生不同的进气值Fig. 2 Different inflection points will lead to different air-entry values

如果土水特征曲线的基本参数用SWCC拟合参数来表示(如ψaev=f(a，m，n))，就会避免传统画图法带来的误差。Zhai et al. (2012)利用Fredlund和 Xing模型导出了一系列的公式，求得进气值和残余吸力值； 高游等(2017)基于修正Fredlund and Xing拟合方程，提出单峰和双峰土水特征曲线基本参数的确定方法； Soltani et al. (2017)利用4种常用的以体积含水率表示的VG模型(m=1;m=1-1/n;m=1-2/n;m≠1)导出进气值和残余状态吸力值的关系式。相比体积含水率，饱和度更能直观地反映非饱和土孔隙中充满水的程度。然而，目前还没有提出有关以饱和度形式表示的VG模型的基本参数的表达式，而且，如何利用VG模型求双峰土水特征曲线基本参数的研究还并未见到。

因此，本文的目的是提出关于VG模型的单峰和双峰土水特征曲线基本参数(进气值，反弯点斜率和残余吸力值)的确定方法，再对滤纸法试验测得的数据点进行最优化拟合，获取拟合参数，验证此方法的可行性，再结合蒸发试验结果，分别得到单峰和双峰SWCC残余状态的确定方法。

1 由VG模型确定基本参数

图3 双峰土水特征曲线基本参数定义Fig. 3 Definitions of bimodal SWCC variables

图4 基于本文方法确定单峰SWCC的基本参数Fig. 4 The variables of unimodal SWCC are determined based on the method in this papera. L6； b. L7； c. L13

土体内部孔隙根据其孔隙分布可分为单峰孔隙分布，双峰孔隙分布和多峰孔隙分布。当土体中存在1组所占比重相对较高的孔隙(优势孔隙)时，土体即为单峰孔隙分布，其内部孔隙为颗粒间孔隙，其土水特征曲线称为单峰土水特征曲线(unimodal SWCC)，如图1所示，曲线只有1个陡降段。当土体中存在两组所占比重相对较高的孔隙(优势孔隙)，土体即为双峰孔隙分布，其内部孔隙主要为集聚体内孔隙和集聚体间孔隙，其土水特征曲线称为双峰土水特征曲线(bimodal SWCC)，如图3所示，曲线具有两个不同斜率的陡降段和一个水平过渡段(Burton et al.，2014; Li et al.，2014)。因此，土水特征曲线按孔隙分布类型可分为单峰土水特征曲线和双峰土水特征曲线。

1.1 单峰土水特征曲线基本参数的确定

要准确确定过渡区的切线方程，反弯点I(ψi，Si)的确定至关重要。在数学中，x是f(x)的拐点的一个必要条件是f(x)的二阶导数应该等于0。因此，VG模型的SWCC函数或式(1)的反弯点(图2中的点I)可以通过以下3个步骤得到：

(1)式(1)先对lg(ψ)求一阶导，可得到任意吸力值所对应的斜率：

(2)

(2)式(2)再对lg(ψ)求导：

(3)

(3)最后令式(3)等于0，可得反弯点I的坐标为:

(4)

把ψ=ψi带入式(2)中，可以得到过反弯点的切线斜率Ki:

(5)

由式(4)和式(5)，可知反弯点的坐标和斜率，据此可得过渡区直线的方程:

(6)

将边界区的水平线方程S=1带入式(6)，可得出进气值ψaev:

(7)

残余状态的确定方法有两种。

方法1：已有研究(Fleureau et al.，1993)表明，当土的含水量接近0时，各种类型土的吸力大致相同，都接近于106,kPa。也就是说，土水特征曲线要经过(106， 0)点。所以经过该点并与SWCC相切的切线(图1中的L3)和反弯点切线L2的交点即为残余点R1。

设L3与SWCC的切点为(ψt，St)，满足两个关系式：

(8)

Kr1由式(2)确定，但是计算发现，解不出ψt和St的解析解。因此只能由点(106， 0)引出一条直线，大致与拟合曲线相切，得到与L2的交点，即为R1(ψr1，Sr1)。

方法2：在曲线近似直线下降的高吸力段(残余区)任取两点(ψ1，S1)、(ψ2，S2)，连成直线L3′，L3′与L2的交点即为R2(ψr2，Sr2)。L3′的斜率Kr2为：

(9)

故直线L3′的方程为：

S=Kr2[lg(ψ)-lg(ψ1)]+S1

(10)

联立式(6)和式(10)，可得残余点的吸力值:

(11)

将式(11)带入式(1)，即解得残余饱和度Sr2。

本文按美国材料与试验协会标准(D18 Committee， 2010)对洛川原状黄土L6、L7、L13进行滤纸法试验，滤纸采用的是国产双圈牌滤纸。待吸力达到平衡后，由滤纸的含水率与滤纸率定公式计算试样的吸力，得到吸力-饱和度的数据点。3层试验黄土的基本物理指标见表1。

表1 洛川原状黄土的基本物理性质Table 1 Properties of Luochuan undisturbed loess

先利用VG模型(式1)对SWCC进行拟合，其拟合参数见表2，再根据上述确定单峰土水特征曲线基本参数的方法，在原拟合曲线的基础上分别画出边界区、过渡区和残余区的直线，见图4，最后确定出土水特征曲线的基本参数(表2)。

表2 单峰土水特征曲线的拟合参数和基本参数Table 2 Best fitting parameters and variables of the unimodal SWCCs

由表2可知，式(1)对试验点的拟合精度较高，均在0.970以上； L6、L7、L13的进气值ψaev，残余吸力值ψr和残余饱和度Sr均随着埋深增大而增大，而过渡区的直线斜率Ki并没有这一规律，而是L7的斜率最小。这可能是由于L7的颗粒组成中黏粒含量较高，细砂粒含量极少，导致土体在相同的吸力增量下排出的水量较少，表现为过渡区斜率最小。由图4可知，采用本文提出的方法可以比较精确地确定土水特征曲线的基本参数。传统图解法主要由操作人员人为确定过渡区切线，有很大的主观性，而本文方法的精度只依赖于试验数据点和拟合参数，精度大大提高。

1.2 双峰土水特征曲线基本参数的确定

双峰土水特征曲线可以看成是两个反“S”型SWCC曲线的叠加，其过渡区由两个斜率的下降段和一个水平段组成，如图3所示。低吸力段SWCC的下降段由集聚体间孔隙控制，ψaev、ψi、Si、Ki、ψr和Sr分别是该下降段的进气值、反弯点I的吸力值、饱和度和斜率、残余点R的吸力值和饱和度； 高吸力段SWCC的下降段由集聚体内孔隙控制，ψ′aev、ψ′i、S′i、K′i、ψ′r和S′r分别是该下降段的进气值、反弯点I′的吸力值、饱和度和斜率、残余点R′的吸力值和饱和度。

以集聚体间孔隙和集聚体内孔隙的孔径分界点D(ψd，Sd)为分界点，将双峰土水特征曲线分为10-1～ψd和ψd～106两部分，用式(1)进行分段拟合，获得拟合参数分别为a0，n0，m0和a1，n1，m1(m满足m=1-1/n)，这些参数将分别用于10-1～ψd段和ψd～106段的土水特征曲线基本参数的确定。

压汞试验得到的孔径分布曲线可获得集聚体间孔隙和集聚体内孔隙的分界孔径Dd，用Young-Laplace方程可换算出等价基质吸力：

(12)

式中，Ts为水-气交界面上的表面张力系数，当水与空气界面温度为20,℃时，Ts=0.072,75N·m-1；α为土颗粒与孔隙水之间的接触角，一般取0。

10-1～ψd段土水特征曲线基本参数的确定：该段反弯点坐标(ψi，Si)及反弯点处的斜率Ki和式(4)、式(5)相同，相应地过渡区直线方程和式(6)相同，将边界区的水平线方程S1=1带入式(6)，可得出进气值ψaev和式(7)相同。式(12)带入式(1)得到Sd，将过渡区水平线方程S2=Sd带入式(6)，可得到10-1～ψd段集聚体间孔隙所对应的残余点坐标R(ψr，Sr)：

(13)

ψd～106段土水特征曲线基本参数的确定：该段SWCC拟合公式需要进行修改为:

(14)

则相应地S′m(ψ)、S″m(ψ)为:

S′m(ψ)=SdS′(ψ)

(15)

S″m(ψ)=SdS″(ψ)

(16)

同样地，令式(16)等于0，可得反弯点I′的坐标为：

(17)

把ψ=ψ′i带入式(15)中，可以得到过反弯点的切线斜率K′i：

(18)

由式(17)和式(18)可得ψd～106段下降段切线方程为：

S=K′i[lg(ψ)-lg(ψ′i)]+S′i

(19)

将过渡区的水平线方程S2=Sd带入式(19)，可得ψd～106段的进气值点A′(ψ′aev，S′aev)，计算结果表明，ψd～106段ψ′aev的表达式与式(7)完全相同，只需要把a，n，m替换为a1，n1，m1即可。

ψd～106段集聚体内孔隙所对应的残余点坐标R′(ψ′r，S′r)的计算与单峰土水特征曲线相同，按1.1节中确定残余状态的两种方法来确定。

表3 双峰土水特征曲线的拟合参数Table 3 Bimodal SWCC best fitting parameters

表4 双峰土水特征曲线的基本参数Table 4 Bimodal SWCC variables

图5 洛川原状黄土L1的孔径分布图Fig. 5 Pore size distribution of Luochuan undisturbed loess L1

图6 基于本文方法确定双峰SWCC的基本参数Fig. 6 The variables of bimodal SWCC are determined based on the method in this paper

2 蒸发试验及残余状态的确定

第1节里提出了两种确定残余状态的方法，由表2和表4可知，这两种方法确定的残余状态相差较大，残余饱和度相差最大的L7达到了11%，而L1的残余吸力值ψr1和ψr2竟相差了2×104,kPa，因此，究竟哪种方法可以有效地确定残余状态？

陶高梁等(2018)通过进行模型拟合、水分蒸发试验以及核磁共振试验，得出这样的结论：水分蒸发试验是确定残余含水率有效可行的直接方法。因此，本文将采用蒸发试验获取残余含水率，从而判断第1节的两种方法的准确性。

取L1、L6、L7、L13的原状环刀样，每层2个平行试验，利用真空饱和法进行饱和，置于温度为25,℃的恒温箱里，自然风干，隔一定时间称量环刀样总质量，得出土中含水量和蒸发时间的关系(图7)。

图7 洛川原状黄土蒸发试验结果Fig. 7 Evaporation test results of Luochuan undisturbed loess

由图7中土中含水量随时间的变化规律可知，不同层原状黄土的蒸发过程可明显分为3个阶段：常速率阶段、降速阶段和残余阶段(唐朝生等， 2011)，符合一般土体水分蒸发的3个阶段。

按陶高梁等(2018)中定义分界时间点的方法，由图7可确定降速阶段和残余阶段的分界时间点(临残时间tr)，tr对应的含水率即为残余含水率ωr。利用图7的试验结果就可以确定残余含水率ωr，最后结合质量含水率与饱和度之间的关系式(20)得到蒸发残余饱和度Srzf，结果如表5所示。分别计算Sr1、Sr2与Srzf之间的相对误差，列于表5。

(20)

式中，ωr为残余质量含水率；Gs为土粒比重；e为土体的孔隙比。

表5 Srzf值以及Sr1、Sr2与Srzf的相对误差计算结果Table 5 Srzf values and the relative errors between Sr1、Sr2 and Srzf

通过表5可以看出，蒸发试验得到的残余含水率在4.52%～7.22%之间，彼此相差不大。对于L1，其Sr1与Srzf的相对误差e1小于Sr2与Srzf的相对误差e2，故建议用方法1确定双峰土水特征曲线的残余状态； 对于其余3层黄土，L7、L13的e2远小于e1，尽管L6的e2大于e1，但这两种方法确定的残余点R1(1002， 0.179)和R2(1136， 0.156)很接近，故建议用方法2来确定单峰土水特征曲线的残余状态。

黄土是由多种矿物组成，其吸附结合水(薄膜水)的能力主要是由其亲水性黏土矿物决定，非黏土矿物对其水合能力的影响很小(王铁行等，2014)。也就是说，黄土残余含水率的大小由其亲水性矿物的含量多少决定。已有研究表明，残余含水率与矿物成分特性有密切关系(Bates et al.，2013)。土样的矿物成分特性不同，其残余含水率就会不同。采用X射线衍射仪(XRD)分析实验黄土的矿物成分，其主要矿物组成成分如表6所示。

表6 试验黄土的主要矿物成分和含量Table 6 The main mineral composition and content of the loess

从表6可以看出试验所用的4层洛川黄土的主要矿物成分为弱亲水性的石英、方解石和斜长石，三者的总量为74%～83.1%。试验黄土中的亲水性黏土矿物为绿泥石和伊利石，这两种亲水性黏土矿物的含量在4层黄土中的分布较为均匀，其含量分别相差为3%、2%，可以认为4层试验黄土中的亲水性黏土矿物的含量相差很小，这可以解释4层黄土的残余含水率值相差不大的原因。

3 结 论

(1)基于VG模型提出了确定单峰和双峰土水特征曲线基本参数的方法。首先由式(1)进行拟合确定其拟合参数，再根据拟合方程的几何关系求出过渡区和残余区的直线方程，最后由直线交点确定土水特征曲线基本参数。双峰土水特征曲线较单峰土水特征曲线的不同是，要根据孔径分界点将双峰土水特征曲线分成两段，进行分段拟合，分段求基本参数。

(2)用洛川标准剖面4层原状黄土的SWCC试验数据分别验证了求解单峰和双峰土水特征曲线基本参数方法的可行性。利用本文方法可较准确地确定土水特征曲线的基本参数。

(3)分析蒸发试验结果，建议用方法1确定双峰土水特征曲线的残余状态，用方法2来确定单峰土水特征曲线的残余状态。