原位非饱和黄土渗透系数研究

文 杰，韩金良，姚磊华，李论基

（1.中国地质大学（北京）工程技术学院，北京 100083；2.中国地质科学院 地质力学研究所，北京 100081）

1 引 言

我国黄土覆盖面积广、厚度大，黄土塬、梁边多发生降雨诱发的浅层滑坡。众多研究表明，非饱和渗透系数函数是分析降雨入渗乃至滑坡发生机制时至关重要的材料参数[1－3]。由于基质吸力的存在，非饱和渗透系数不同于饱和土，是随基质吸力或者体积含水率变化的函数，但具体为何种函数关系，尚难有定论[4－6]。

目前，常用的非饱和渗透系数的测定方法分为直接法和间接法[5]。间接法是以基质吸力-体积含水率关系曲线为基础，运用统计模型来确定非饱和渗透系数的一种方法[7－9]。Lu 等[5]指出，统计模型适用于砂土等孔径分布范围较窄的粗颗粒土，细颗粒土，如黄土则不适合采用统计模型来进行预测。

直接法主要为瞬态剖面法（室内、野外），瞬态剖面法测定非饱和渗透系数由Richard和Weeks 于1953年提出，如今已成为测定非饱和渗透系数最常用的试验方法。目前，对瞬态剖面法的研究主要集中于对室内瞬态剖面法试验装置的改进以及对数据的处理方法上[4－5]，关于原位瞬态剖面法测定非饱和渗透系数的研究较少，相比室内瞬态剖面法而言，原位瞬态剖面法更接近自然条件，计算结果更具实际应用价值。除此之外，Libardi[10]、邵明安[11－12]、张建丰[13]等提出了在田间测定条件下非饱和渗透系数的各种解法，上述田间测定方法在黄土边坡中的适用性尚需进一步的验证。

本文拟从自然条件下黄土降雨入渗监测数据出发，选取适用于瞬态剖面法、θ 法的数据计算得到黄土的非饱和渗透系数，分析黄土的非饱和渗透系数与体积含水率之间的关系，对比分析瞬态剖面法及θ 法计算结果，论证θ 法在非饱和黄土边坡计算中的适用性。同时，尝试根据实际情况对瞬态剖面法的数据处理以及θ 法的基本假定作出修正，以求得到更符合实际的结果。

2 降雨入渗监测

降雨入渗监测点位于宝鸡市渭滨区峪泉镇中岩山村7 组。2011年8月雨季曾经发生滑坡，现今一遇雨季，老滑坡体前缘仍会出现次级滑坡，为典型的降雨入渗诱发的浅层黄土滑坡。

据钻孔资料及露头分析，滑坡底界为三门系红黏土，滑体为覆盖于三门系红黏土之上的黄土，为揭示出降雨条件下诱发该滑坡的机制，在距离滑坡不远处埋置降雨入渗数据采集设备。

降雨入渗数据采集主要为含水率与负压监测，监测设备采用美国Decagon 公司研制的EC-5 土壤水分传感器、MPS-2 土壤水势传感器（见图1），本次监测共埋设了5 个EC-5 土壤水分传感器通道、5 个MPS-2 土壤水势传感器通道，每个土壤水分传感器通道和水势传感器通道埋置深度一致，水平距离为0.01 m。每个通道的埋置深度如表1 所示。

表1 传感器编号及埋置深度Table 1 Sensor number and embedment depth

2014年7月16日12︰00 开始使用监测仪器采集数据，7月20日以后采集间隔设置为30 min。

在监测时段内，根据宝鸡市气象局提供的2014年9月6～18日的降雨量数据，结合体积含水率与负压监测数据绘制成图2、3。

图2 9月6～18日体积含水率随降雨量变化曲线Fig.2 Change curves of volumetric water content changes with rainfall from September 6-18

图3 9月6～18日负压随降雨量变化曲线Fig.3 Change curves of negative pressure changes with rainfall from September 6-18

由图2、3 分析可知，体积含水率总的变化趋势是降雨时体积含水率增大，降雨停止后体积含水率减小，基质吸力的变化趋势为降雨入渗后基质吸力迅速降低；通道5 所在位置靠近地表，具有一定的波动性，变化范围大；通道1 位于土层较深处，体积含水率、基质吸力变化范围小，并且在降雨入渗发生较长一段时间后，体积含水率、基质吸力发生改变，反映了水分运移到土体深处需要一段时间，具有一定的滞后效应。

原位瞬态剖面法、θ 法等试验要求一个较长的雨水入渗过程后土体一定深度内含水率接近饱和并且处处相等，然后在不发生蒸发作用或者蒸发作用微弱并可忽略不计的情况下进行水分重分布过程[4－5]。

对2014年7月16日开始的监测数据进行分析，并结合中岩山村所在地区的实际气象情况，自9月28日开始的监测数据满足瞬态剖面法等的要求，相应的监测数据变化见图4、5。其中9月28日或更早的时间至10月2日为降雨入渗过程，土体内体积含水率接近饱和且处处相等，10月2日至11月21日为土体水分的重分布过程，该过程中无降雨过程发生，期间中岩山地区持续阴天，温度较低，场地表面无植被覆盖，可忽略蒸发作用效应，满足相关原位试验原理性要求。

图4 9月28日至11月21日含水率变化情况Fig.4 Variation of volumetric water content from September 28 to November 21

图5 9月28日至11月21日负压变化Fig.5 Variation of negative pressure from September 28 to November 21

3 瞬态剖面法

3.1 原理

瞬态剖面法是一种适用于室内或现场确定非饱和渗透系数的瞬态测量技术。该方法是指在瞬态流过程中量测一定深度内的负压和体积含水率剖面，负压水头加上重力水头即总水头，水头剖面用于计算不同深度处的特定时间的水力梯度[4－5]。某一点的流速是从体积含水率剖面上的变化计算出来。原位瞬态剖面法要求先将试验区土体用水浸湿，直到土体在水流作用下达到饱和状态，这时通过地面的水流停止。将试验区表面进行一定的处理以防止蒸发和入渗，以后水在土体内进行的过程即为水分重分布过程（即非稳态过程），然后开始量测重分布过程中含水率和吸力的变化数据用于计算。

由前述第2 节的说明及给出的图4 看出，2014年10月2日之前体积含水率持续上升接近饱和，10月2日之后开始下降，因此，以10月2日为瞬态剖面法计算当中的非稳态过程的开始，即取10月2日至11月21日体积含水率、负压数据作为瞬态剖面法计算用数据。

3.2 计算过程

排水过程中土层任意位置的总水头 hw等于该位置的位置水头z与测得的负压水头 hm之和（计算水头时，取地表面为基准面，向下为负），即

在特定的时间和特定的深度内，水力梯度可以从该深度水头剖面的坡度计算出来：

式中：iw为特定深度和特定时间内的水力梯度；为所考虑深度处水头剖面的坡度。

考虑向后差分格式的优越性，使用向后差分格式对式（2）进行求解，在某一计算时间段内，某深度的平均水力坡度计算如下：

式中：hi,t1、hi,t2分别为 t1和t2时刻位置 zi处的水头；hi+1,t1、hi+1,t2分别为 t1和t2时刻位置zi+1处的水头。

水的流速由不同深度和时间的体积含水率剖面得出，通过地面的水流通量为0，特定时间在地面(z=0)和地面下某一深度z 之间水的总体积可以按体积含水率剖面积分的方法计算：

式中：Vw为地面及深度z 之间的土中水体积；θw(z)为体积含水率，为深度z 的函数；A为试验区的平面面积。

使用式（4）计算一定深度范围内总含水率的变化时需假定体积含水率的分布函数，最常用的方法为分段线性拟合，该方法对数据的处理较为粗糙，张玉莲[14]、王红[15]等选用3 次样条方程拟合体积含水率分布函数，样条曲线有其自身的优越性，但其物理意义不太明显。本次试验所在土层为均匀的黄土层，其体积含水率随深度的分布应为连续变化，本次计算拟用连续性函数来拟合体积含水率分布函数。考虑以2014年10月2日为非稳态过程的开始，即t=0 d，绘制6、10、15、23、30、40、50 d 的体积含水率，对深度z 曲线，如图6 所示。拟合结果见表2。

图6 不同时刻含水率随深度的变化Fig.6 Variation of volumetric water content with depth at different times

表2 不同时刻含水率-深度的拟合公式Table 2 Fitting equations between volumetric water content and depth at different times

因此，选用对数函数来拟合体积含水率的分布情况（考虑对数函数参数必须大于0 的情况，深度为从0.000 1 m 开始；拟合常数为a、b），即

提取表2 中相关数据作拟合常数a、b 随时间变化曲线，如图7 所示，拟合公式见表3。

图7 拟合参数a、b 随时间变化Fig.7 Variation of a,b and time

表3 a、b 对时间拟合关系Table 3 Fitting relationship between a,b and time

根据不同时刻含水率对深度的拟合公式（5）及拟合常数a、b 对时间的拟合关系，可得非稳态重分布过程体积含水率与深度和时间的关系式为

因此，起始位置以下一定深度内总含水率Vw计算公式为

考虑曾经因为降雨入渗试验的需要，对以仪器埋置处为中心的4 m×4 m 范围内进行简单处理。

地面以下一定深度z 内，从两个相邻时间间隔dt 计算出来的水体积差dVw为该时间间隔内流过该深度点的水量，该点的流速 vw按下式计算：

该流速相应于两个相邻时间所得水力梯度的平均值。将流速 vw除以平均水力梯度 iave可算出渗透系数 kw为

对不同深度点和不同时间重复进行非饱和渗透系数的计算，因此，在一个试验中可计算出对应于不同含水率的渗透系数值。

由于通道5 所在位置接近地表，负压数据测量误差较大，因此，选取通道1～4 的数据进行计算。

选取4 个通道t 分别为6、10、23、40 d 数据用式（8）进行计算，计算结果列出于表4 中。

表4 瞬态剖面法计算结果Table 4 Calculated results of transient profile method

分析表4 中数据可知，存在两个渗透系数异常点，即平均体积含水率为39.2 %和38.7 %时，即1.32～2.22 m 深度范围内对应的平均含水率，其异常原因在于最深处不能使用向后差分法计算水力坡度，只能取与上一个深度范围相同的水力坡度，从而导致计算误差出现异常。

剔除异常点后，将计算所得渗透系数与平均含水率数据绘制成图8，拟合数据后得黄土的非饱和渗透系数与体积含水率变化曲线及关系式（10），相关系数为0.944 8。结果表明，黄土的非饱和渗透系数与体积含水率成指数型关系，当含水率较大并发生一定变化时渗透系数变化较大；含水率较小时，含水率发生同样的变化导致的渗透系数变化较小。

图8 瞬态剖面法计算结果及拟合结果Fig.8 Calculated results of transient profile method and fitting result

4 θ 法

1980年Libardi 等[10]提出了一种在田间测定非饱和渗透系数的方法，即θ 法，该法假设蒸发及入渗作用微弱，地下水的重分布只受到重力场的作用。他们将此法分别用于巴西的淋浴土、巴西的氧化土、美国加州的Panohe 黏壤土。任理[16]在河北省南宫试验场进行试验，用θ 法处理数据得到非饱和渗透系数，并用数值模拟验证了所得参数的可靠性。但该法在黄土中适用性尚未得到证明。

4.1 基本假定

（1）通过积水入渗，在地面下一定深度范围内（0

（2）t＞0时，土体表面（z=0处）的水通量为0；

（3）准饱和层内只有重力作用，各处吸力S 相等，即∂S/∂z=0，故 ∂φ/∂z=−1；

（4）土体内部在0～z 深度内平均含水率θ*与某一给定深度z 的体积含水率θ 之间存在线性关系，θ∗=mθ+n，m、n 均为拟合常数；

（5）非饱和渗透系数与体积含水率之间遵从如下关系式：

式中：β为拟合常数；θ0为水入渗一段时间（地下水重分布开始时）的准饱和层内体积含水率；k0为相应的渗透系数（饱和渗透系数）。

4.2 公式推导

描述地下水一维垂向运动的基本方程为Richards 方程：

对式（12）进行积分，得

式中：φ为总水势，包括吸力水头S和位置水头z两项。

由基本假定（2）可知，式（13）的最后一项为0，因此式（13）可改写为

式（14）左边项表示土壤深度z 以内贮存的水量随时间的变化率，可简化为

又根据基本假定（3）～（5），可将式（15）进一步简化为

在某一给定深度z 对式（16）从初始条件（t=0，θ=θ0）积分至任意时刻（t=t，θ=θ）可得

当t 较大时，式（17）可简化为

4.3 试验数据

由于自然降雨具有大面积分布的特点，故可视降雨入渗过程中地下水流具有垂向一维性。通过图4 曲线的变化趋势看出，整个过程可以分为两部分：9月29日或更早时间至10月2日为降雨入渗过程，此过程完成后整个土壤剖面内体积含水率基本相同且达到最大值，整个土壤剖面可视为基本假定（1）所定义的准饱和层；10月2日至11月21日中4 个通道体积含水率持续减小至稳定，表明期间无降雨渗入，该过程可视为在重力作用下的土壤水分重分布过程。因此，将10月2日作为土壤水分重分布过程的开始，对应于任理[16]所做野外试验的重分布试验的开始时刻，相应的体积含水率视为0θ，计算选用的数据同样为10月2日至11月21日数据。

4.4 计算步骤和结果

（1）由于通道5 临近地面，测量结果具有一定的波动性，同时深度越大，水分重分布过程所需时间过长，因此选用10月2日以后通道2～4 的测量数据进行计算。

（2）对通道2、3、4 的测量数据均进行平均含水率 θ*与θ 的线性回归（3 个计算深度至地面范围内，通过将该通道以上所有通道数据求算术平均值可得），得到的计算结果见图9和表5。

（3）对每一计算点进行θ0−θ与lnt 的线性回归（t 的单位为h）。由于假设t 较大，式（18）计算中选用时间靠后的数据进行拟合，其效果更佳，计算结果见图10和表6。

图9 θ*与θ 线性回归趋势Fig.9 Linear regression trend of θ* and θ

表5 θ*与θ 线性回归方程Table 5 Linear regression equations of θ* and θ

图10 θ0−θ与lnt 线性回归趋势Fig.10 Linear regression trend ofθ0−θ and lnt

表6 θ0−θ与lnt 线性回归方程Table 6 Linear regression equations ofθ0−θ and lnt

（4）根据式（18）、基本假定（4）及表5、6中各线性回归方程中的参数，由基本假定（5）即可求得3 个深度处黄土的非饱和渗透系数的函数表达式，将各公式中的系数取算术平均，并化成单位m/d，即可得到本方法所求结果，具体见表7。

4.5 对θ 法的修正

θ 法的基本假定（4）直接给出地表以下一定深度z 至地表范围内的平均含水率与该深度z 处的含水率呈线性关系，θ∗=mθ+n，m、n为拟合常数，并且随深度变化而变化，同时该假定中平均含水率直接通过若干个离散点数据取算术平均值求得，计算过程略显粗糙。因此，考虑3.2 节中给出的体积含水率随深度的分布函数对该基本假定进行重新设定。

根据3.2 节可得体积含水率随深度的分布函数如式（5），而平均含水率为含水率分布函数对深度进行积分并除以深度（同样考虑对数参数必须大于0，将积分起点设置为0.000 1），计算可得平均含水率分布函数如式（19）：

因此，对于任何深度处，平均含水率与该深度处含水率应满足以下条件：

将该式（20）代入式（16）重新进行积分得最终结果为

式（20）表明平均含水率与某一深度处含水率呈线性关系，且其斜率恒定为1，因此，在计算渗透系数过程中无须进行平均含水率计算以及某一深度处含水率的拟合分析。式（21）表明，θ0−θ与lnt仍呈线性关系，且斜率不变，但截距发生改变。应用式（21）重新计算得到修正的θ 法所求渗透系数的结果，见表8。

表8 修正θ 法计算结果Table 8 Calculated results of modified θ method

修正的θ法与θ法的计算结果差异体现在中间结果 k0（饱和渗透系数）上，将两种方法在不同通道处计算出的 k0值列于表9 中。

表9 修正θ 法与θ 法计算的 k0值Table 9 Calculated k0 of modified θ and θ methods

计算结果对比显示，修正θ 法的计算中间值 k0具有较小的离散性，从一定程度上说明修正θ 法的可靠性。

6 结果对比

将瞬态剖面法结果数据及相应体积含水率下θ法、修正的θ 法得到结果绘制于图11 中。

图11 3 种计算方法结果对比Fig.11 Comparison of three methods

由瞬态剖面法计算结果（图8）分析可知，黄土的非饱和渗透系数与体积含水率呈指数关系，证明θ 法的假定（1）的合理性。从图11 中可看出，瞬态剖面法与θ 法、修正的θ 法计算结果在数量级上具有一致性。同一体积含水率下，瞬态剖面法计算结果普遍大于θ 法，这是由于θ 法中对式（17）进行简化时需依赖于t 较大的情况，因此，依据式（18）进行拟合时，需选用t 较大时的数据，但实际当t 较大时水分重分布作用接近完成，渗流缓慢，因此，相比于瞬态剖面法，θ 法计算的结果偏小。由于对θ 法中平均含水率的求取进行了接近实际情况的修正，使得平均含水率与体积含水率拟合关系中的斜率减小，从而导致所求的黄土的非饱和渗透系数较θ 法小，而计算的中间结果饱和渗透系数 k0值具有更小的离散性。

由于θ 法作了t 较大的假定，求出的黄土的非饱和渗透系数较瞬态剖面法偏小，但结果与瞬态剖面法相比具有数量级上的一致性，同时θ 法计算过程中不需要负压数据，减少了野外测试的工作量，并且计算过程极为简单，可用于黄土的非饱和渗透系数的测定。在使用过程中需根据实际体积含水率的分布情况对基本假定（4）作出相应的修正，以求得到更为合理的结果。

7 结 论

（1）使用瞬态剖面法处理黄土降雨入渗监测数据，得出黄土的非饱和渗透系数与体积含水率之间的关系呈指数关系。

（2）θ 法求取非饱和渗透系数无需测定基质吸力，试验方法、求解过程简单，对比瞬态剖面法计算结果，发现θ 法亦适用于黄土的非饱和渗透系数的求解，值得推广。

（3）通过分析实际情况，尝试用对数函数来拟合不同时间体积含水率与深度的分布关系，并将此关系应用到瞬态剖面法计算过程以及对θ 法的改进中，提高了计算结果可靠度。

（4）由于θ 法及相应修正方法作了较多假定，瞬态剖面法计算结果大于θ 法以及修正的θ 法；基于含水率实际分布的修正θ 法计算结果与θ 法计算结果相差不大，表明含水率的分布对θ 法测定结果影响很小，但修正θ 法的计算中间结果饱和渗透系数 k0值具有更小的离散性。

[1]吴礼舟,黄润秋.非饱和土渗流及其参数影响的数值分析[J].水文地质工程地质,2011,38(1):94－98.WU Li-zhou,HUANG Run-qiu.A numerical analysis of the infiltration and parameter effects in unsaturated soil[J].Hydrogeology &Engineering Geology,2011,38(1):94－98.

[2]简文星,许强,童龙云.三峡库区黄土坡滑坡降雨入渗模型研究[J].岩土力学,2013,34(12):3527－3533.JIAN Wen-xing,XU Qiang,TONG Long-yun.Rainfall infiltration model of Huangtupo landslide in Three Gorges Reservoir Area[J].Rock and Soil Mechanics,2013,34(12):3527－3533.

[3]简文星,许强,吴韩,等.三峡库区黄土坡滑坡非饱和水力参数研究[J].岩土力学,2014,35(12):3517－3522.JIAN Wen-xing,XU Qiang,WU Han,et al.Study of unsaturated hydraulic parameters of Huangtupo landslide in Three Gorges Reservoir Area[J].Rock and Soil Mechanics,2014,35(12):3517－3522.

[4]FREDLUND D G,RAHARDJO H.Soil mechanics for unsaturated soils[M].Ottawa:Wiley-Interscience,1993.

[5]LU N,WILLIAM J.Unsaturated soil mechanics[M].Hoboken:John Wiley and Sons Inc.,2004.

[6]VAN GENUCHTEN M T,LEIJ F J.The RETC code for quantifying the hydraulic functions of unsaturated soils[R].Riverside:U.S.Department of Agriculture,Agricultural Research Service,Report IAG-DW12933934,1991.

[7]李萍,李同录,王红,等.非饱和黄土土-水特征曲线与渗透系数Childs和Collis-Geroge 模型预测[J].岩土力学,2013,34(增刊):184－189.LI Ping,LI Tong-lu,WANG Hong,et al.Soil-watercharacteristic curve and permeability prediction on Childs&Collis-Geroge model of unsaturated loess[J].Rock and Soil Mechanics,2013,34(Supp.):184－189.

[8]刘海宁,姜彤,刘汉东.非饱和土渗透函数方程的间接确定[J].岩土力学,2004,25(11):1795－1799.LIU Hai-ning,JIANG Tong,LIU Han-dong.Indirect determination of permeability function equation of unsaturated soils[J].Rock and Soil Mechanics,2004,25(11):1795－1799.

[9]叶为民,钱丽鑫,白云,等.由土-水特征曲线预测上海非饱和软土渗透系数[J].岩土工程学报,2005,27(11):1262－1265.YE Wei-min,QIAN Li-xin,BAI Yun,et al.Predicting coefficient of permeability from soil-water characteristic curve for Shanghai soft soil[J].Chinese Journal of Geotechnical Engineering,2005,27(11):1262－1265.

[10]LIBARDI P L,REICHARDT K,NIELSEN D R,et al.Simple field methods for estimating soil hydraulic conductivity[J].Soil Science Society of America Journal,1980,44:3－7.

[11]邵明安,王全九.推求土壤水分运动参数的简单入渗法:I 理论分析[J].土壤学报,2000,37(1):1－8.SHAO Min-an,WANG Quan-jiu.A simple infiltration method for estimating soil hydraulic properties of unsaturated soils,I:Theoretical analysis[J].Acta Pedologica Sincia,2000,37(1):1－8.

[12]邵明安,王全九.推求土壤水分运动参数的简单入渗法:II 实例验证[J].土壤学报,2000,37(2):217－224.SHAO Ming-an,WANG Quan-jiu.A simple infiltration method for estimating soil hydraulic properties of unsaturated soils,II:Experimental results[J].Acta Pedologica Sincia,2000,37(1):217－224.

[13]张建丰,王文焰.野外一维垂向入渗试验确定土壤水分运动参数[J].水土保持学报,1994,8(1):69－72.ZHANG Jian-feng,WANG Wen-yan.Parameters determination of soil water movement by field one-dimensional experiment of vertical infiltration[J].Journal of Soil and Water Conservation,1994,8(1):69－72.

[14]张玉莲.非饱和土渗透系数瞬态剖面测量方法及仪器的改进[D].哈尔滨:哈尔滨工业大学,2010.

[15]王红,李同录,付昱凯.利用瞬态剖面法测定非饱和黄土的渗透性曲线[J].水利学报,2014,45(8):997－1003.WANG Hong,LI Tong-lu,FU Yu-kai.Determining permeability function of unsaturated loess by using instantaneous profile method[J].Journal of Hydraulic Engineering,2014,45(8):997－1003.

[16]任理.野外条件下非饱和土壤水力传导度的确定[J].水利学报,1989,(11):49－55.