1、引言

土壤持水曲线（soil-water retention curve,SWRC）是非饱和土基本的本构关系之一，它描述了饱和度与吸力之间的关系[1-3]。SWRC在干燥和湿润路径上有所不同，表现出明显的滞后特征[4]，显著影响非饱和土的渗流特性[5-6]、应力-应变关系[7-8]和抗剪强度[9-10]。理解土壤持水滞后机制可对SWRC滞后模型的构建提供重要参考。

为解释非饱和土持水行为的滞后特征，研究人员提出了不同的物理机制，如墨水瓶效应[11-12]、接触角滞后[13]以及土壤基质表面的水合作用[14]等，且基于上述物理机制提出了不同的SWRC滞后模型。在模拟复杂水力路径下土体的持水特性时，这些模型的预测能力通常优于传统模型[15-16]。如，Mualem[11]采用独立域假设提出了一个概念性的SWRC滞后模型，其可依据主边界曲线预测SWRC的扫描行为；Zhou[13]分析了接触角滞后机制，建立了SWRC滞后增量模型；Chen等[12]将土体孔隙空间简化为非均匀毛细管，并假设孔隙尺寸分布遵循分形定律，建立了一个可以描述持水滞后的SWRC分形模型。此外，Li[17]采用边界面塑性理论建立了可考虑水力滞后的SWRC增量模型。基于热力学理论，Wei等[18]通过引入不可逆含水率作为与水力滞后相关的内变量，建立了热力学水力滞后SWRC模型。刘艳等[19]综合考虑了可恢复和不可恢复含水率，改进了Wei等[18]的SWRC模型。徐炎兵[20]和Chen[21]等通过考虑残余含气量对土体水力滞后的影响分别建立了能模拟任意干湿循环路径下土壤持水滞后的SWRC模型。然而，试验研究表明[22-23]，土体初始孔隙比减小会导致SWRC滞回环显著收缩，而既有SWRC模型无法解释和描述这一现象背后的物理机制。

体积变形对土体持水特性的影响已受到广泛关注[24-27]。因此，研究人员提出了不同的SWRC改进模型来反映体积变形的影响[28-34]，其主要可分为两类：经验模型和理论模型。如陈正汉等[28]采用张力计研究了初始干密度对压实黄土持水特性的影响，并提出了双参数SWRC经验公式。该公式中的两个参数与孔隙率呈线性关系，能够较好地模拟不同初始孔隙率条件下压实黄土的SWRC。Huang等[29]在Brooks-Corey模型（BC模型）[30]中引入了孔隙比的单参数函数，提出了可描述孔隙比影响的BC模型；Gallipoli等[31]在Van Genuchten模型（VG模型）[15]中引入了幂函数以描述体积变形对饱和度的影响；基于Fredlund-Xing模型[16]，周葆春[32]和邹维列[33]等均将进气值表示为初始孔隙比的函数，提出了可考虑初始孔隙比变化的Fredlund-Xing扩展模型；Zhou等[34]基于饱和度和体积变形的内在约束，提出了一个与初始孔隙比相关的公式以量化初始干密度对SWRC的影响；陈正汉[1]详尽地阐述了初始孔隙比等因素对土体持水性质的影响。

对于基于特定物理机制的SWRC理论模型，如Simms等[35]基于网络方法研究了体积变化对非饱和土壤水力特性的影响；陶高梁等[36]假设土体孔隙尺寸分布（pore size distribution,PSD）服从分形分布，分析了孔隙比变化对进气值的影响，提出了可反映初始孔隙比影响的SWRC分形模型；胡冉等[37]假定变形状态下的PSD可通过缩放和移动参考状态下的PSD来获得，推导出了可描述体积变形对饱和度影响的SWRC改进模型。前述模型均通过不同方式反映体积变形对PSD和进气值的影响。但值得注意的是，Chen等[22]的试验研究表明，体积变形不仅会改变土体的PSD，也会改变孔隙结构的异质性（或孔隙结构的非均匀性）。土体压缩会显著降低孔隙结构异质性，弱化土体的墨水瓶效应，进而使得SWRC滞回环显著收缩。然而，既有SWRC模型重点关注了体积变形对PSD的影响，但忽略了体积变形对孔隙结构异质性的影响，因此难以解释和描述滞回环随土体孔隙比减小而收缩的现象。

本文通过土体中三相比例的内在约束，探讨了常吸力下体积变形对土体饱和度影响的边界约束条件，提出了一个可描述体积变形对饱和度影响的增量公式。此外，基于墨水瓶效应分析了SWRC的滞后特性，假设吸力增量主要用于改变土体饱和度和克服墨水瓶效应，提出了一个物理模型以定量描述SWRC沿任意水力路径的水力滞后行为和滞回环随初孔隙比变化的演化特征。最后，采用11组试验数据对提出的模型进行了验证。

2、模型建立

饱和度对非饱和土的水力耦合行为有着显著影响[38-39]。在等温条件下，土体饱和度Sr受水力状态（基质吸力和水力路径）和力学状态（孔隙比）的影响，可表示为

式中：s为基质吸力；ζ为与水力路径相关的变量；e为孔隙比。忽略变量s、ζ和e之间的耦合关系，对式（1）微分可得

式（2）可进一步表示为

式（3）为非饱和土广义SWRC模型，其右侧第1项反映了吸力和水力路径的影响，第2项表示体积变形对土体饱和度的影响。

2.1常吸力下体积变形对饱和度的影响

在恒定吸力下，土体饱和度随孔隙比的变化规律是影响土体性质的重要边界条件之一。试验研究[40-42]表明，饱和度随着孔隙比的减小是增加还是减小，取决于土壤的性质和状态[43]。然而，现有的SWRC模型只能描述饱和度随孔隙比减小而增加的情况，而无法解释其随孔隙比减小而降低的现象。针对这一问题，依据Zhou等[34]提出的土体三相比例的内在约束，提出一个能够全面描述常吸力下饱和度随孔隙比变化的增量公式。

Zhou等[34]假设土体在变形过程中，孔隙体积与孔隙水体积存在如下关系，即

式中：d Vw为孔隙水体积增量；d Vv为孔隙体积增量；ω为比例系数（0≤ω≤1）。此外，根据饱和度与孔隙比定义可得

式中：Vs为土颗粒体积。由式（5）可得

将式（6）、（7）代入式（4）可得

结合式（5）和式（8）可得

由于0≤ω≤1，因此由式（9）可得

除体积变形外，影响饱和度的因素还有吸力和水力路径，因此，在常吸力下，∂Sr/∂e需满足式（10）的约束条件，即

式（11）即为常吸力下，孔隙比变化对饱和度影响的内在约束条件。

因此，根据式（11），提出如下公式来描述常吸力下孔隙比变化对土体饱和度的影响，即

式中：0<Sr≤1;β为拟合参数，其可通过给定吸力下的Sr-e数据拟合得到。式（12）可与任何既有SWRC模型相结合，以预测不同初始孔隙比下土体的持水特征。如图1所示为预测曲线与实测数的对比，可知式（12）既可描述∂Sr/∂e大于0的情况（见图1(a)），也可描述∂Sr/∂e小于0的情况（见图1(b)）。需要指出的是，图1(a)中的箭头表示饱和度随孔隙比的减小而减小的变化趋势。

图1 常吸力下孔隙比变化对土体饱和度的影响

2.2 SWRC水力滞后模型

由孔隙结构异质性引起的墨水瓶效应是导致SWRC具有滞后效应的主要原因。土体的非均匀孔隙可以用图2所示的几何形状来简化描述。如图所示，半径较大（R1）的部分称为孔体，半径较小（R2）的部分称为孔喉。孔喉与孔体之比被定义为异质性因子η(0<η≤1）。η越小，则表明孔隙结构的非均匀程度越大，土体的墨水瓶效应也越强。

孔隙半径和吸力之间的关系可以用毛细管定律来描述[16]，即

式中：在平衡条件下，s=ua‒uw，其中ua为孔隙气压，uw为孔隙水压；c=2Tscosθ，其中Ts为孔隙水的表面张力，θ为固-液接触角，假定其为常数。因此，当水气接触面位于孔体和孔喉处时，相应的吸力分别为c/R1和c/R2。

假设图2所示的毛细管处于一个封闭系统中，在该系统中可以应用轴平移技术（水压保持不变，吸力的变化等于气压的变化）。需要注意的是，当基质吸力大于水气压差时（s>ua‒uw），孔隙水会向上流动；当s<ua‒uw时，孔隙水会从毛细管中排出。

图2 非均匀孔隙中的持水行为

首先，毛细管中充满水，处于干燥路径。当封闭系统中的气压从0开始逐渐增大，但ua‒uw值小于c/R1时，毛细管不会排水。虽然ua‒uw增大，但毛细管的饱和度不会发生变化（Sr=1），因此SWRC是一条平行于吸力轴的直线（见图2(c)中的AB段）。当ua‒uw达到c/R1时，毛细管上部孔体中的孔隙水被排干，水气接触面位于孔喉顶部。此时，毛细管的饱和度为Sr,1。如果ua‒uw继续增大但小于c/R2，则毛细管饱和度在此过程中保持不变，如图2(c)中CD段所示。当ua‒uw>c/R2时，毛细管中的水全部排出，饱和度为0。

随后，土体处于湿化路径，如图2(b)所示。当ua‒uw逐渐减小到c/R2时，毛细管中的半月板上升到孔喉和孔体的交界处（见图2(b)）。此时毛细管的饱和度为Sr,2。随着ua‒uw继续减小但大于c/R1，毛细管饱和度不会发生变化（见图2(c)中的EF部分）。当ua‒uw减小到c/R1时，毛细管中又会充满水。从图2(c)可以看出，在这种干湿循环中，SWRC具有明显的滞后特性。根据图2(c)可知，在给定的饱和度下，干燥曲线和润湿曲线之间的距离可表示为

根据式（14）可知，土体孔隙结构越均匀（即异质性程度越小），η越大，进而使得SWRC的滞后程度越小。若孔隙为横截面均匀的光滑毛细管（η=1），则SWRC的滞后效应将消失。由此可见，SWRC的水力滞后是由孔隙结构的异质性引起的，因此式（3）中的变量ζ可用η代替，即

Chen等[22]的试验研究表明，孔隙比变化不仅会改变PSD，还会改变孔隙结构的异质性。由于半径较小的孔喉比半径较大的孔体更难被压缩，孔隙空间的异质性在压缩过程中将减弱。因此，土体墨水瓶效应会随着初始孔隙比的减小而降低，进而导致滞回环随着孔隙比减小而收缩的现象。

如图3所示[22]，土体的初始孔隙比减小时，孔体的体积发生了显著减小，而孔喉体积变化很小。因此，若忽略孔喉体积的变化，根据式（14）则有

式中：η'和为土体在变形后的异质性因子和孔体半径。

图3 孔体与孔喉体积变化特征

试验证据[45]表明，对于黏土或者粉土，等效孔径ln R与孔隙比e呈线性关系，即ln R∝e。一方面孔体体积显著大于孔喉体积，土体孔径R主要是孔体半径贡献的；另一方面，土体在压缩过程中，由于孔喉的半径变化很小，而孔体半径发生了很大的变化，因此孔径R与孔隙比的关系可近似假设为孔体半径R1与孔隙比e的变化关系，即ln R1∝e，故有

式中：κ为比例系数；e'为土体在变形后的孔隙比。

结合式（16）和（17）可得

基于前述滞后机制，下面将构建一个SWRC模型，用于定量描述非饱和土壤沿任意水力路径的持水行为。由于主干燥曲线数据相对容易通过室内试验获得，下文将采用主干燥曲线作为参考曲线来构建SWRC模型。在参考状态下的主干燥曲线可用VG模型[15]描述，即

式中：sd为主干燥曲线上的吸力；md和nd为与孔径分布有关的参数；αd为与参考状态下土体进气值有关的参数。

式中：λ为孔径变化对异质性因子影响程度的参数；ηe=ηminexp[κ(e-e')],ηmin为湿润路径上异质性因子的最小值。

在湿润路径，孔隙结构异质性随土体变形和吸力变化而演化的规律则可表示为

结合式（12）、（17）、（19）、（20）和（21）可得

式（22）即为提出的可模拟水力滞后与体积变形耦合的非饱和土SWRC增量模型。该模型共包含7个参数，其中4个（αd、nd、md和β）来源于主干燥曲线，可以过拟合两条不同初始孔隙比下的主干燥曲线试验数据确定；其余参数（κe、λ和ηmin）是孔隙结构参数，可以通过拟合扫描曲线试验数据来确定。这些参数可用于描述土体在不同初始孔隙比下滞回环的演化特性。但若只描述给定初始孔隙比下土体的水力滞后性质时，仅需其中5个参数（αd、nd、md、λ和ηmin）；若只描述初始孔隙比对土体主干燥或主湿润曲线的影响时，仅需其中4个参数（αd、nd、md和β）。

3、模型验证与讨论

本节采用了不同土壤（包含砂土、粉土和黏土，且经历了复杂的应力和水力路径）的11组试验数据对提出的模型进行验证。

3.1常吸力下土体变形对饱和度的影响

图4为3种不同土壤在常吸力下的饱和度Sr与孔隙比e的关系。由于只描述常吸力下土体的孔隙比与饱和度的关系，因此仅需联合式（12）和式（9），模型参数为αd、nd、md和β。

图4 恒定吸力下体积变形对土壤饱和度的影响

对于Speswhite高岭土[42]，参数αd、nd、md和β可通过吸力为150 k Pa和1 k Pa下的试验数据进行校准（αd=0.2 k Pa-1、nd=1.13、md=0.12、β=10），然后用于预测其他常吸力下的Sr-e关系。如图4(a)所示，预测曲线与试验数据吻合得很好。图4(b)和图4(c)分别为高岭土[46](αd=0.015 k Pa-1、nd=1.23、md=0.19、β=16）和松散压实土[47](αd=0.005 k Pa-1、nd=1.5、md=0.34、β=3）在不同吸力下Sr-e试验数据与预测曲线的对比。如图4(a)～4(c)所示，在所有示例中，预测曲线均与试验数据吻合良好。另一方面，如图1(a)所示，所提出的模型还可描述土体在压缩过程中饱和度随孔隙比减小而降低的情况，进一步表明了提出模型的可行性。

3.2变初始孔隙比下的SWRC特性

图5为Speswhite高岭土[48]和Municipal膨胀土[49]在不同初始孔隙比下的持水特征，因此亦仅需联合式（12）和式（9）进行预测，包含参数为αd、nd、md和β。其中，参数αd、nd、md可通过拟合任意初始孔隙比下的SWRC数据确定，β可通过拟合任意额外一条SWRC数据确定。

图5 初始孔隙比对SWRC的影响

如图5所示，初始孔隙比的减小会导致进气值增大，且在给定吸力下饱和度也会相应增加。对于Speswhite高岭土，αd、md、nd和β的值分别为0.009 k Pa-1、0.18、1.22和1.9，而对于Municipal膨胀土，这些参数分别为0.001 k Pa-1、0.23、1.3和2.0。预测曲线与试验数据吻合良好，表明所提出的SWRC模型可有效描述初始孔隙比变化对土体持水特征的影响。

3.3给定初始孔隙比下SWRC的水力滞后特性

Gillham等[50]通过伽马射线衰减研究了不同水力路径下Dune砂的持水特征，其试验结果如图6所示。由于只描述给定初始孔隙比下土体的水力滞后性质，因此，本小节仅需5个参数（αd、nd、md、λ和ηmin），未用到式（17）。

图6 Dune砂的持水滞后行为

模型参数αd、nd和md可由VG模型根据主干燥曲线校准，参数ηmin和λ可由任意扫描曲线确定（αd=0.3 k Pa-1、md=0.87、nd=7.5、ηmin=0.55、λ=0.6）。只需确定水力路径的反转点，即可用校准的参数预测任意扫描干燥曲线、扫描润湿曲线及主润湿曲线。如图6所示，预测曲线与试验数据吻合良好。

3.4不同初始孔隙比下的滞回环演化特性

对于吸力较高的情况（>1 MPa），所提出的模型仍能较好地描述SWRC的滞后行为。图7展示了吸力范围较大的Viadana粉土的SWRC数据[51]。如图7(a)所示，初始孔隙比e0=0.66的Viadana粉土先沿主湿润路径加载至350 k Pa，再沿扫描湿润路径卸载至10 k Pa，然后沿扫描干燥路径加载至105 k Pa，最后沿主湿润路径卸载至1 k Pa。在这两个干湿循环（ABCDA）过程中，SWRC表现出明显的滞后特性。同样，参数αd、nd、md、ηmin和λ可通过主干燥曲线和任意扫描曲线确定（αd=0.003 5 k Pa-1、md=0.65、nd=0.78、ηminref=0.15、λ=0.40和κ=3.39），所获得的参数可用于预测任意其他水力路径下的SWRC。如图7(a)所示，预测曲线与试验数据吻合良好。

图7 Viadana粉土在不同初始孔隙比下的水力滞后行为

另一方面，如前所述，孔隙比的变化不仅会改变土壤的孔径分布，还会改变孔隙空间结构的异质性，因此会显著影响SWRC的滞后行为。图7(b)为初始空隙比e0=0.93的Viadana粉土的SWRC数据。比较图7(a)和图7(b),e0=0.66的Viadana粉土的滞回环尺寸明显小于e0=0.93下的滞回环尺寸。结合从e0=0.66下Viadana粉土的SWRC中获得的参数，可通过e0=0.93的Viadana粉土的单一扫描曲线来校准参数β (β=3）。然后，所获得的参数可用于预测Viadana粉土在任何其他孔隙比和水力路径下的扫描行为。如图7(b)所示，预测曲线与试验数据吻合得很好。

Chen等[22]对不同孔隙比下的高岭土和南阳黏土进行了干湿循环试验，研究了墨水瓶效应弱化对土体持水滞后行为的影响，其试验结果如图8所示。对于高岭土而言（αd=0.005 3 k Pa-1、md=0.38、nd=0.95、ηmin=0.05、λ=2、κ=8.3和β=0.3），如图8(a)所示，初始孔隙比e0=0.88的滞回环尺寸显著小于初始孔隙比为1.10的滞回环尺寸。若任意取某一饱和度，如Sr=0.75,e0=1.10对应的ln(1/η)=1.73，显著大于e0=0.88对应的ln(1/η)=0.90，这表明滞回环随土体初始孔隙比减小而显著收缩。图8(b)所示的南阳黏土也表现出了相似的特征（αd=0.003 1 k Pa-1、md=0.37、nd=0.6、ηmin=0.01、λ=1、κ=8.5和β=0.9）。由图8(a)和图8(b)可知，预测曲线与试验数据吻合得很好，表明提出的模型可较好地描述滞回环随土体体积变形而演化的特征。

图8 滞回环随初始孔隙比减小而收缩的特性

3.5复杂水力耦合路径下的SWRC特性

Raveendiraraj[42]进行了一系列吸力控制的非饱和土持水试验，研究了Speswhite高岭土（Test-B6）在如图9(a)所示复杂应力和水力路径下的持水行为，SWRC试验结果如图9(b)和9(c)所示。模型参数确定步骤如下：首先确定水力路径初始点C，采用式（21）拟合图9(b)中的湿润扫描曲线CD，标定参数αd=0.009 5 k Pa-1、md=0.12、nd=2.54、ηmin=0.25和λ=0.34；采用式（12）拟合图9(c)中AB段，标定参数β=10。吸力引起孔隙比变化的试验数据可直接用于计算土体饱和度。如图9所示，预测曲线与试验数据吻合良好。

图9 Speswhite高岭土（Test-B6）在水力耦合下的SWRC

Fig.9 SWRC of Speswhite kaolin (Test-B6) under hydro-mechanical coupling

图10为高岭土（Test-A9）在复杂水力和应力路径下的持水特性试验结果。与Test-B6中的模型参数确定方法相同，采用式（21）拟合图10(b)中的湿润扫描曲线CD，标定参数αd=0.033 k Pa-1、md=0.057、nd=1.06、ηmin=0.5和λ=0.9；采用式（12）拟合图10(c)中AB段，标定参数β=4。如图10所示，预测曲线与试验数据吻合良好。

图1 0 Speswhite高岭土（Test-A9）在水力耦合下的SWRC

4、结论

从体积变形和水力路径对土壤持水特性的影响出发，建立了能够考虑水力滞后与变形耦合的SWRC模型，并通过试验数据验证了该模型的可行性，主要工作和结论如下：

(1）基于试验数据，根据土体中三相比例的内在关系，提出了一个可用于描述体积变形对SWRC影响的经验增量方程，该方程能够较为全面地描述饱和度随着孔隙比的减小而增大或减小的现象。

(2）采用墨水瓶效应揭示了土体持水滞后的物理机制，并假设吸力增量主要消耗在改变土体饱和度和克服墨水瓶效应两个方面，构建了可描述任意水力路径下土体持水特征的理论模型。

(3）孔径较小的孔喉比孔径较大的孔体更难被压缩，将导致土体孔隙结构随着孔隙比的减小而更加均匀，进而弱化墨水瓶效应，引起土体持水滞后程度减小。因此，通过将恒吸力下饱和度随孔隙比变化的增量方程与基于墨水瓶效应的理论模型相耦合，建立了可定量描述和解释滞回环随初始孔隙比变化的SWRC模型。最后，采用不同水力和应力路径下的试验数据验证了模型的可行性。

参考文献:

[1]陈正汉.非饱和土与特殊土力学[M].北京:中国建筑工业出版社, 2022.

[3]陈正汉,苗强强,郭楠,等.关于持水特性曲线研究的几个问题[J].岩土工程学报, 2023, 45(4):671-679.

[4]蔡国庆,田京京,李舰,等.考虑变形及滞回效应影响的三维土-水特征曲面模型[J].土木工程学报, 2019,52(11):97-107.

[5]周葆春,陈志.压实膨胀土非饱和渗水系数函数的密度与水力滞回效应[J].岩土工程学报, 2019, 41(10):1800-1808.

[8]孙德安.非饱和土力学特性及本构模型[J].岩土工程学报, 2023, 45(1):1-23.

[19]刘艳,赵成刚.土水特征曲线滞后模型的研究[J].岩土工程学报, 2008, 30(3):399-405.

[20]徐炎兵,韦昌富,陈辉,等.任意干湿路径下非饱和岩土介质的土水特征关系模型[J].岩石力学与工程学报,2008, 27(5):1046-1052.

[24]高游,李泽,孙德安,等.考虑初始孔隙比影响的单/双峰土-水特征曲线模型研究[J].岩土力学, 2022, 43(6):1441-1452.

[25]周葆春,孔令伟,梁维云,等.压缩过程中非饱和膨胀土体变特征与持水特性的水力耦合效应[J].岩土工程学报, 2015, 37(4):629-640.

[26]翟钱,戴国亮,赵学亮.土-水特征曲线对非饱和砂土抗剪强度的影响[J].岩土工程学报, 2020, 42(7):1341-1349.

[27]林志强,钱建固,时振昊.毛细-吸附作用下考虑孔隙比影响的单/双峰土体持水曲线模型[J].岩土力学,2021, 42(9):2499-2506.

[28]陈正汉,谢定义,王永胜.非饱和土的水气运动规律及其工程性质研究[J].岩土工程学报, 1993, 15(3):9-20.

[32]周葆春,孔令伟.考虑体积变化的非饱和膨胀土土水特征[J].水利学报, 2011, 42(10):1152-1160.

[33]邹维列,张俊峰,王协群.脱湿路径下重塑膨胀土的体变修正与土水特征[J].岩土工程学报, 2012, 34(12):2213-2219.

基金资助:国家重点研发计划(No.2022YFC3803000,No.2016YFC0800200)~~;

文章来源:陈可,王琛,梁发云,等.考虑水力滞后与变形耦合的非饱和土持水曲线模型[J].岩土力学,2024,45(12):3694-3704+3716.