土体非饱和渗透特性对降雨入渗的影响

陈旺旺，刘 刚，刘 畅，钟佳思，童富果

(三峡大学 水利与环境学院， 湖北 宜昌 443002)

连续性降雨、强降雨等是触发泥石流、滑坡等地质灾害的主要诱因[1-4]。非饱和土中雨水入渗，产生水气二相流的运动，继而改变边坡土壤的非饱和特性，引发边坡失稳[5-6]。土体的非饱和渗流特性包含了土水特征关系、水-气渗透关系，是土体微观结构在宏观上的综合反映。因此，进行不同非饱和渗流特性条件下的入渗过程进行数值模拟，对于揭示多个物理特性变量下土体的入渗规律具有重要意义。

非饱和渗流及其对边坡稳定性的影响有不少学者做了相关研究。朱伟等[7-8]通过设计一维降雨入渗试验，研究了非饱和渗流特性对降雨入渗水量的影响。丰光亮等[9]在室内模拟了人工降雨入渗非饱和土柱试验，提出鄂西恩施地区入渗影响区入渗前锋运移规律和降雨强度、降雨历时对入渗率的影响规律。Neuman[10]将有限元方法应用到求解饱和-非饱和渗流问题。陈善雄等[11]和陈守义[12]用积分有限差分方法模拟了降雨条件下土体中水分的运动情况，并对降雨条件下非饱和土坡等的稳定性的分析方法做了研究。上述的试验方法和数值模拟多建立在液相单相流理论上[13]。Lam等[14]假定土壤包含饱和和非饱和区域，在不饱和区域气体为连续体且与大气相通，利用二维有限元模型计算非饱和土体的应力状态。孙冬梅等[15]运用积分有限元差分方法模拟了土质边坡的降雨入渗过程，定量分析了空气阻力对水流入渗的影响。降雨入渗过程是一个涉及水气二相流耦合的非饱和渗流过程，以往涉及土体非饱和渗透特性的研究多集中于外界因素对降雨入渗规律的探讨，而与非饱和渗透特性直接相关的土水特征关系、水-气渗透关系对降雨入渗的影响也非常值得研究。

因此，本文基于非饱和土-水-气二相流理论，采用有限单元法对不同非饱和渗透特性条件下的入渗过程进行计算，探究了土水特征关系和水相对渗透曲线对稳定入渗强度的影响。对了解降雨诱发边坡失稳的机制，提高滑坡的预报和防治具有重要意义，也对城市防洪、农业节水灌溉、污染物质传输等诸多工程问题提供参考。

1 计算方程与数值模型

1.1 水气两相流控制方程

土体视为液、固、气三相混合的连续介质，各相物质在土中的运动可描述为关于时间和空间的偏微分方程组。土壤内部的液相流动主要受水压力梯度的驱动，水分通过土壤中的水流通道入渗，土体的断面面积大于实际水流的过水面积，因此Darcy定律计算的水流速度小于土中水流的真实流速。根据液相质量守恒方程得到孔隙水的非饱和流动方程[16]：

(1)

式中：ρw为水的密度，kg/m3；μw为液体黏滞系数，Pa·s；Krw为水相对渗透系数；K为本征渗透系数，m2；Sr为水饱和度；φ为土壤孔隙率；Pw为孔隙水压力；Qw为内源项。

气相在土壤中的流动受到重力和气体压力梯度的驱动，同时土中气体仅通过气流通道，土体的断面面积大于实际气流通过的面积，故由Darcy定律得来的气流速度小于土中气流的真实流速。根据气相质量守恒方程得到孔隙气的非饱和流动方程[16]：

(2)

式中：ρg为气的密度，kg/m3；μg为液体黏滞系数，Pa·s；Krg为气相对渗透系数；Pg为孔隙水压力。

在求解水、气二相流控制方程时，孔隙水压力和孔隙气压力之间恒有Pc≡Pg-Pw，式中Pc为基质吸力。方程(1)、方程(2)构成时间和空间的非线性偏微分方程，求解时对方程组的空间离散化采用Garlerkin有限单元法，时间域离散采用一维差分方法，即假设在时间域的时间段内变量随时间的线性变化。

1.2 本构模型关系与参数

方程(1)和(2)有Pw，Pg，Sr，Krw，Krg5个未知量，求解时还需要引入3个本构模型关系，即土水特征关系、水相对渗透关系、气相对渗透关系，反映了土体的含水率与基质吸力、渗水能力、渗气能力的关系。

土水特征关系选取适用于大多数土体的Van Genuchten模型[17]：

(3)

式中：Sr为水饱和度；Srw为残余水饱和度；P0为进气值；m为与土体孔径分布有关的参数。

水相对渗透关系选取为Van Genten-Mualem模型[18]：

(4)

式中：k为与材料特性有关的参数。

气相对渗透关系选取Brooks & Corey模型[19]：

(5)

式中：Srg为残余气饱和度。

方程涉及的其他计算参数为：ρw=1 t/m3，ρg=1.29 kg/m3，g=9.8 N/kg，φ=0.4，μw=1.0×10-3N·s/m2，μg=1.84×10-5N·s/m2，k=2.88×10-13m2。

1.3 几何模型

对于土质边坡而言，自然降雨的入渗方向主要为垂向，水分入渗也主要集中在坡体浅层区域，故本文在数值模拟时用到的几何模型为浅层均质土体。几何模型的大小由降雨入渗过程中水分入渗深度决定，而水分的入渗深度与土体特性、边界条件等因素相关的同时，还与降雨入渗时间相关。现选定几何模型高3 m，顶部单元厚度为0.1 m，底部单元厚度为0.2 m，共计42个节点20个单元。为研究土体非饱和渗透特性对入渗强度的影响，假定入渗强度远小于降雨强度，地表径流对入渗的影响可忽略不计，顶部及底部均为透水透气边界，且其边界气压力为大气压，侧边界为不透水不透气边界。

2 土体非饱和渗流特性与入渗强度的关系

2.1 土水特征关系对入渗强度的影响

土水特征关系表述了土壤基质吸力和土壤饱和度的关系，本文采用的Van Genuchten模型中，参数m和参数P0值的变化会直接影响关系曲线的结果形式，本节将对这两个参数分别研究。为便于分析比较，采用相对入渗强度之比来衡量渗流特性变化对入渗强度的影响强弱。入渗强度之比是一个无量纲变量，定义为不同渗透特性下的相对入渗强度与某一渗透特性下的相对入渗强度的比值。

(1) 参数m变化对入渗强度的影响。参数m是反映土体孔隙结构特征的重要参数之一，与孔径分布指数相关。通过模拟不同参数m条件下的入渗过程，探究孔径分布对土壤入渗强度的影响。根据相关资料可知参数m的常见取值范围[20]为0.3～0.5，计算过程中，参数m取0.30、0.35、0.40、0.45、0.50一共5组(见图1)，饱和度取值从0.15～1.00共18组，参数P0=0.278，k=0.9，相互组合共计模拟90组不同入渗过程。基质吸力随饱和度的变化率因参数m的改变而改变，参数m越小，基质吸力在饱和度越低时的影响越大，而当饱和度较高时，参数m对基质吸力的影响较小。

图1不同参数m条件下土水特征曲线

如图2所示，参数m的改变只影响饱和度较低时的基质吸力，随着饱和度的增大，参数m的改变对基质吸力影响不再明显，故参数m主要影响饱和度较低时的土体稳定入渗强度，当饱和度较高趋近于1时，土壤内部的入渗通道是水分入渗的主要通道，几乎不受基质吸力的影响，故土体稳定入渗强度几乎不受参数m影响，且其极小值的大小和位置不随参数m的改变而改变。

图2相对入渗强度与饱和度的关系曲线

(2) 参数P0变化对入渗强度的影响。参数P0为进气值，通过模拟不同参数P0条件下的入渗过程，反映土体的最大孔隙尺寸对入渗强度的影响。参数P0的取值范围为0.2～2.0，计算过程中，参数P0为0.2、0.6、1.0、1.4、1.8一共5组(见图3)，饱和度取值从0.15～1.00共18组，参数m=0.359，k=0.9，相互组合共计模拟90组不同入渗过程。基质吸力随饱和度的变化率因参数P0的改变而改变，参数P0越大，基质吸力在饱和度越低时的影响越大，随着饱和度的增大，水相对渗透系数对入渗起主导作用，参数P0对基质吸力的影响逐渐减小直至为0。

图3不同参数P0条件下土水特征曲线

当饱和度增大时，入渗强度先减小后增大，相同饱和度的条件下，入渗强度随参数P0的增大而增大(见图4)。因为P0增大使得饱和度较高时的基质吸力增大，入渗强度受基质吸力的影响明显，故入渗强度极小值点向饱和度较高处偏移。当饱和度趋近于1时，水分主要是通过土壤孔隙通道流动，水相对渗透系数是主要影响因素，基质吸力此时的作用并不明显。

图4相对入渗强度与饱和度的关系曲线

2.2 水相对渗透曲线对入渗强度的影响

水相对渗透曲线反映了水相对渗透系数与饱和度的关系，是影响入渗强度的关键因素。本文采用Van Genten-Mualem模型，模型中参数k反映了水相对渗透系数随饱和度的变化规律。本文参数k的取值范围为0.5～0.9，计算过程中，参数k为0.5、0.6、0.7、0.8、0.9一共5组(见图5)，饱和度取值从0.15～1.00共18组，参数P0=1.33，m=0.4，相互组合共计模拟90组不同入渗过程。水相对渗透系数随着饱和度的增大而增大，土壤残余含水和完全饱和时的水相对渗透系数分别为0和1，相同饱和度下水相对渗透系数随着参数k的减小而减小。

图5不同参数k条件下的水相对渗透曲线

由图6可以看出稳定入渗强度是取决于基质吸力和水相对渗透系数两者结果，水相对渗透系数减小时入渗强度减小，则参数k越小稳定入渗强度也越小。随着饱和度的增大，水相渗透系数的变化会影响孔隙通道的过水能力，所以饱和度较高时入渗强度受参数k的变化影响明显。土体稳定入渗强度的极小值位于饱和度较高时，参数k的减小会导致水相对渗透系数的降低进而导致入渗强度的极小值减小。

图6相对入渗强度与饱和度的关系曲线

3 结 论

(1) 在非饱和入渗过程中，基质吸力的变化直接影响水分入渗过程中的水压力梯度。当饱和度较低、参数m较小时，基质吸力对土体的入渗强度影响比较明显；随着饱和度的增大并趋近于1时，土体入渗强度几乎不受参数m的影响。

(2) 当饱和度较低时，基质吸力驱动水分入渗，入渗强度随着参数P0的增大而变化明显；当饱和度较大时，水相对渗透系数较大，孔隙通道作为水分主要下渗通道，参数P0增大对入渗强度的影响微弱。

(3) 土体稳定入渗强度的极值点不仅受基质吸力和相对渗透系数的影响，还与土体本征渗透性关系密切，土体本征渗透性决定了土体稳定入渗强度极小值的位置。随着本征渗透性的增大，土体渗水能力增强，极小值的位置向低饱和度运动，此时基质吸力的影响可忽略不计；当本征渗透性足够小时，入渗受基质吸力的影响，此时稳定入渗强度随饱和度的增大而减小。