基于未冻水含量曲线的冻土渗透系数模型研究∗

2022-08-02张致龙关陈晨

新疆大学学报(自然科学版)（中英文） 2022年4期

张致龙，关陈晨

(新疆大学建筑工程学院，新疆乌鲁木齐 830017)

0 引言

水分迁移在冻土的冻胀和融沉中起着重要作用[1]，精确评估冻土中的水分迁移需要确定渗透系数.在冻土孔隙中，液态水和孔隙冰共同存在，孔隙冰的产生可降低冻土的渗透系数，土壤中冰含量随温度的变化而变化，因此，冻土的渗透系数也随温度变化.表征冻土中水分迁移的一个主要挑战是如何在冻结和融化过程中确定温度、初始含水率和渗透系数之间的变化关系[2−3].

由于冻土的渗透系数测试比较困难，在这方面的研究有限.冻土渗透系数的测量一般采用两种方法：直接测量法和间接测量法.冻土渗透系数的直接测量首先由Burt等[4]提出，其设计了一种等温条件下液体直接通过冻土的渗透仪，并在测试中采用低浓度的乳糖和聚乙二醇替代水分.Horiguchi等[5]以油为流体测量了0∼-0.35 ℃之间土壤冻结过程中的渗透系数，发现冻土渗透系数随温度的降低而急剧减小.张虎等[6]以二甲基硅油作为流体，对含水量为62%的青藏高原粉质黏土在0∼-0.6 ℃范围的渗透系数进行了测量.间接测量法是假设冻土中的水分迁移与不饱和融土中的水分迁移相同，这种假设是基于融土中干燥和湿润过程与冻土中冻结和融化过程的液态水变化相似[7]，认为土水特性曲线可以代替未冻水含量曲线[8].

理论计算方面，都假设冻土中的未冻水含量和融土中的含水量变化与渗透系数的函数关系是相同的.Flerchinger等[9]采用这一假设计算冻土渗透系数，发现当液态水含量小于某一阈值时，必须人为设置冻结土壤的渗透系数为零.Lundin[10]认为应用融土的渗透系数方程计算冻土的水流量会使计算值过高，因此加入了一个阻抗参数.Watanabe等[11]建立了一种毛细管束模型来描述冻土中的水流，但此模型同样采用土水特性曲线（SWCC）代替未冻水含量曲线（SFCC），并且忽略了孔隙壁对水流的影响.然而，粉质黏土中3∼100 nm孔隙的累积体积约占总孔隙体积的一半，试验和理论研究发现，由于小孔隙（纳米孔）壁对水滑移和有效粘度的影响，在相同压力梯度作用下，小孔隙中的水通量相对体积水显著变化[12−14]，比通过不同直径碳纳米管膜且使用无滑移的Hagen-Poiseuille方程计算的水通量大1到5个数量级，在接触角较小的纳米孔中，流速还会降低[15−16].此外，测试土水特性曲线比较困难，并且花费时间较长，但随着核磁共振（NMR）和时域反射计（TDR）技术的发展，使得未冻水含量曲线的测量更为容易，并且能更准确地反应未冻水含量的变化情况.因此，一些研究人员认为应用直接测量的未冻水含量曲线能更准确地预测冻土的渗透系数[17].

考虑到上述使用SWCC的局限性，并且以上理论模型都没有考虑孔隙壁对水的滑移和有效粘度的影响.为了更准确、更容易地预测冻土的渗透系数，本文提出一种利用未冻水含量曲线、考虑孔隙壁对水的滑移和有效粘度影响的新模型来描述冻土中水的迁移，并讨论了土壤冻结过程中渗透系数的变化，将计算结果与试验值进行比较，验证了模型的合理性.

1 理论推导

此模型基于Hagen-Poiseuille方程，采用冻土的未冻水含量曲线，并考虑了冻土中孔隙壁对水流的影响，推导出饱和冻土渗透系数与初始含水率、温度的函数关系.

1.1 冻土中的孔隙

当温度低于土的冻结温度时，冰先从土中的大孔隙中心生长，但由于孔隙壁表面力（范德瓦耳斯力）的作用，在孔隙壁附近和小孔隙中的水仍然保持液体状态.如图1所示，已冻大孔隙中的未冻水膜厚度d可以表示为[18]：

式中：A为Hamaker常数，A=-10−19.5[17]；ρi为冰密度，ρi=917 kg/m3；Lf为冰的融化潜热，Lf=3.34×105J/kg；Ti为土的冻结温度.由Gibbs-Thomson方程，未冻小孔隙的临界直径（未冻孔隙最大直径）γGT可以表示为[18−19]：

式中：σls为冰水界面能，σls=2.90×10−2J/m2.

冻土中的孔隙包含两种：含冰孔隙和无冰孔隙，含冰孔隙所占总孔隙的比例随温度的变化而变化.未冻孔隙的临界直径和已冻孔隙的未冻水膜厚度与温度变化关系可以由式（1）和式（2）计算得到.如图2所示，当土温度低于-1 ℃时，未冻孔隙的直径小于100 nm，已冻孔隙的未冻水膜厚度小于5 nm，均为纳米孔隙，因此，必须考虑孔隙壁对水流的影响.由于未冻孔隙的临界直径比已冻孔隙的未冻水膜厚度大十几倍，所以，冻土中的液态水主要是通过未冻孔隙迁移的.图3显示了不同初始含水量Nanticoke黏土[20]和青藏高原粉质黏土[21]的孔径累积体积分布曲线，由图3可知，3∼100 nm的孔隙累积体积约占低含水土壤孔隙总体积的一半.研究表明，孔隙壁附近的水粘度比体积水的粘度高，水与小孔隙壁相互作用使水的物理性质表现出剧烈的变化[22−23].因此，研究水在冻土中迁移时必须考虑水在小孔隙壁上的滑移和有效粘度的变化.

1.2 冻土中的水流

在压力梯度作用下，水流过一个给定直径的孔隙，不考虑水在孔隙壁上的滑移、层流，由力平衡得出Hagen-Poiseuille方程：

式中：J为水通量；µ∞为体积水的粘度，µ∞=9.62×10−7Pa·s；R为孔隙直径；∂p/∂z为压力梯度.当温度低于-1℃时，由于冻土中未冻孔隙的直径都小于200 nm，考虑水在孔壁上的滑移和有效粘度的变化，水通量由式（3）变为[24]：

式中：ls,t为实际滑移长度；µ(d)为孔隙水的有效粘度.水的实际滑动长度作为土-水接触角的函数可以描述为[25]：

式中：θ为土-水接触角，θ=25˚[26]；C为常数，C=0.41[24].

为了得到冻土孔隙中水的有效粘度，采用孔隙中的界面和体积水粘度的加权平均值[27]：

式中：µi为水的粘度；Ai为孔隙中土-水接触面积；At为孔隙横截面积.µi/µ∞与接触角θ为线性关系[24]：

为了理解孔隙壁对水通量的影响，引入了增强因子，增强因子定义为水通量的测量值与J的比值[12−13]，如果公式（4）可以描述实测水通量，则增强因子可以表示为：

如图4所示，由公式（8）计算出的增强因子与冻土孔径的关系.计算结果表明：由于孔隙壁对水流的影响，冻土中水的流量随着孔隙直径的减小而降低，最小可以降低为原来的50%，降低的速率也随着孔隙直径的减小逐渐变大.

1.3 冻土渗透系数

假设在圆柱形的冻土样品中形成一个网状的孔隙结构，每一个孔隙的长度为L，土样的长度为H（H

式中：∂z=0 ∼H；ξ为土柱长度H与水流过实际路径L的比值（H/L），ξ=0.6[11]；M为单位面积上某一孔径的孔隙数目；nJ为单位面积上孔径RJ的孔隙数目.

利用未冻水含量曲线（未冻水含量wu与温度T）建立一个实际土样和模型之间的对应关系，即将未冻水含量曲线中的未冻水含量分成若干宽度为△wu的等分，冻土中的未冻水含量wu(T)与温度的变化相关，类比达西定律，由式（9）可以得到冻土的渗透系数为：

式中：k=0,1,2,···，未冻水含量每下降△wu对应的k增加1.依据式（10），当k=0时，K为相同条件下融土的渗透系数，因为融土中水分以在大孔隙中流动为主，忽略孔隙壁对水流的影响，融土的渗透系数可以表示为：

式中：K0为相同含水量融土的渗透系数.因此，冻土的渗透系数也可以表示为：

依据式（1）和式（2），当孔隙直径RJ>γGT（冻结孔隙）时，RJ=d.冻土中的未冻水含量等于各孔隙中的液态水之和，可以表示为：

式中：孔径为RJ的孔隙数目nJ可以由式（13）和未冻水含量曲线计算得到.本冻土渗透系数模型仅需要准确测量土的未冻水含量曲线和在相同条件下融土的渗透系数，因此，可以更容易和准确地计算冻土的渗透系数.

2 理论模型的验证

2.1 试验材料和测试方法

试验在冻土工程国家重点实验室进行，选用的土样为青藏高原粉土和粉质黏土，与文献[6]的土样相同.将土样在实验室内风干、碾碎、过2 mm筛后分成若干份备用，土样的粒径分布如图5所示.为了制备不同含水率的饱和土壤样品，在土样中加入预设含量去离子水，搅拌均匀，放入密封的塑料袋中24 h，使水分充分均匀，然后进行24 h真空饱和.采用干燥法测定土壤样品的含水率，实测体积含水量分别为46.3%、51.2%、53.5%、62.2%、50.4%（粉土），并计算干密度，结果见表1，并将干燥法测得土样含水率作为标准值.为减小因试样非均质带来的误差，将土样分别逐层均匀地填埋于圆柱形塑料容器内，制成直径为6.18 cm、高为12 cm的圆柱形试样.

表1 土壤样品的干密度和初始含水率

试验中选用的5TM土壤水分传感器，是由Decagon公司生产的ECH20土壤水分传感器系列之一，体积含水量测量范围为0∼100%，工作温度范围为-40∼50 ℃.5TM土壤水分传感器可以测量土壤或介质中的介电常数，由于土壤中水的介电常数远大于空气与其它基质的介电常数，因此，土壤的介电常数大小主要取决于土壤中水分的含量，通过测量土壤的介电常数就可以达到测量土壤中水分含量的目的.试验时将温度设置为阶梯型下降，温度分别为0.5、-0.5、-1、-1.5、-2、-3、-5、-10 ℃的5TM探头可以测得冻结过程中土壤介电常数的变化.

核磁共振方法被广泛应用于测量未冻水含量，利用纽曼公司生产的冻土核磁共振仪（NMR），对体积含水量为51.2%土样在冻结过程中的未冻水含量变化进行测量.NMR核磁共振法可以快速精确地在实验室获得冻土中的未冻水含量，NMR法测量的是磁场中质子的自由感应衰变值FID，而样品中水分含量与水的FID值成一定比例关系，所以通过测量负温下相应介质的FID值就可以根据相应的换算公式得到准确的未冻水含量.通过对应5TM测得的相同土样的介电常数值，得到在冻结过程中未冻水含量与介电常数的关系曲线，得出转化公式，如图6所示，从而依据测出的介电常数，计算出其它的未冻水含量值.

2.2 未冻水含量

图7为未冻水含量（体积含水率）随温度的变化关系.如图7所示，在0∼-0.5 ℃时，未冻水含量显著降低，且相同土质下，初始含水量较大的土样下降幅度更大，如初始含水量为46.3%和62.2%的粉质黏土，未冻水含量分别从46.3%减小到约20%以及从62.2%减小到约10%，这是因为相同土质下初始含水量越大，土中大孔隙总体积占比越高，在温度低于冰点且大孔隙冻结后，未冻水含量下降幅度就越大；当温度在-0.5∼-2.5 ℃时，较大孔隙开始冻结，未冻水含量下降幅度减缓；当温度在-2.5∼-5 ℃时，土体中的未冻水主要存在于约50 nm（图2）的孔隙中，因此未冻水含量变化幅度很小；同时，在相同温度下，初始含水量越小的土样，未冻水含量越高，这是因为初始含水量越小，未冻小孔隙体积占比越大，而未冻水主要存在于未冻小孔隙中，因此，随着初始含水量的减小，未冻水含量逐渐增加.

3 结果与讨论

3.1 理论模型的验证

为了验证模型的可靠性，将模型的计算结果与其它试验结果进行对比.如图8所示，依据式（1）、（2）、（12）和（13）对不同土壤和含水量的渗透系数进行计算.粉质黏土的渗透系数计算采用与文献[4]和[6]相同的土样和初始含水率测得的未冻水含量曲线，粉土的初始含水量与文献[8]接近.由图8可知，本模型计算的冻土渗透系数与试验结果吻合较好，由于初始含水率和土壤类型（黏土、高岭土）的不同，其它测试结果与计算值略有不同，但也在初始含水率51.2%∼62.2%的计算结果范围内，这说明计算结果是可靠的.本模型能更容易、更准确地计算冻土渗透系数，提高预测冻胀和融沉的准确性.这是因为土的冻胀由两部分组成，原位水冻胀和迁移水冻胀，水分迁移量大小主要由驱动力和渗透系数决定，试验测得冻土渗透系数非常困难，一般采用替代溶液法，且需冻土温度大于-1.5 ℃，因此，准确给出不同含水率与土质在负温下的渗透系数是精确预测冻胀融沉的关键.

3.2 冻土的渗透系数

图9为不同初始含水量饱和土样的渗透系数随温度的变化关系.在0∼-0.5 ℃时，冻土的渗透系数随着温度的降低而快速减小，减小幅度受初始含水量和土质类型的影响，如：粉土从10−6cm/s减小到10−11cm/s，含水率62.2%的粉质黏土从约10−6cm/s减小到10−10cm/s，这是因为当温度低于冻结温度时，大孔隙中的液态水首先相变成冰，在初始含水量较高的粉质黏土和粉土中，大孔隙的体积占总孔隙体积的比例更高，未冻的小孔隙体积占的比例较小，而未冻水主要在小孔隙中流动，使得高含水量的粉质黏土和粉土的渗透系数更小；结合图2可知，当温度达到-0.5 ℃时，孔径约200 nm以上的孔隙全部冻结，使得未冻孔隙的体积占比急剧减小，因此，在冻结初期，土壤的渗透系数下降剧烈；在-0.5∼-2.5 ℃时，各土样渗透系数的降低速率减缓，减小幅度基本相同，约减小了2个数量级，水分迁移困难；在-2.5∼-5 ℃时，由于各土样渗透系数都小于10−10cm/s，粉土已小于10−13cm/s，虽然各土样渗透系数最大相差2个数量级，但渗透系数都很小，水分迁移非常困难；小于-5 ℃的冻土渗透系数更小，一般计算认为水分不再迁移，因此，不进行更低温度渗透系数的计算[30−31].这是因为温度低于-2.5℃之后，孔径约50 nm以上的孔隙全部冻结，不仅未冻孔隙占比更少，而且水分迁移还受到孔隙壁作用.

图10为-0.5 ℃时饱和冻土的未冻水含量、渗透系数与初始含水量之间的关系.由图10可知，未冻水含量和渗透系数随初始含水量的增大而减小，且变化趋势基本相同.这是因为在-0.5 ℃时，直径小于216 nm的孔隙中液态水并没有冻结（图2），随着饱和土初始含水量增加，孔隙中所容纳的水分体积增加，土颗粒的间距增大，部分小孔隙变为较大孔隙（γGT>216 nm）而冻结，未冻水主要处于小孔隙中（γGT<216 nm），而土中的小孔隙数量和总体积减小，在相同条件下，低含水量土样中的未冻小孔隙的总体积大于高含水量土样.因此，随着饱和冻土的初始含水量增加，未冻水含量和渗透系数逐渐减小.

3.3 冻融过程中冻土的渗透系数

为了避免重复分析，我们只测量了初始含水量为46.3%的冻结粉质黏土在0∼-5 ℃冻融过程中的未冻含水量，如图11所示，融化过程的未冻含水量小于相同温度下的冻结过程.计算对应的渗透系数如图12所示，在相同温度下，融化过程的渗透系数大于冻结过程，差值超过1个数量级，在温度低于-3 ℃后，渗透系数小于10−11cm/s，水分迁移困难.