非饱和带排水及吸水过程中水分分布特征研究

张圣山，卫云波，顾 熙

1. 湖州中交投资发展有限公司，湖州 31300；2. 河海大学 地球科学与工程学院，南京 210013；3. 南京大学 地球科学与工程学院，南京 210023

1 引言

排水和吸水过程是非饱和带中两个相反的水气驱替过程。排水过程指的是多孔介质中含水量减少（即气驱水）的过程，而吸水过程指的是含水量增加（即水驱气）的过程（Lenormand and Zarcone,1984）。而根据吸水过程中水相压强（pw）与气相压强（pa）之间的大小关系，吸水过程又可以被细分为三类：自吸（即 free imbibition，此时pw=pa）、受迫吸水（即forced imbibition，此时pw>pa）以及受限吸水（即constrained imbibition，此时（pw

因为“瓶颈效应”的存在，孔隙尺度上的排水和吸水过程存在很大差别，如排水过程中存在不稳定的“海恩斯跳跃”现象（Haines, 1930）、临近孔隙干扰排水（Smiles et al., 1981）现象，以及吸水过程中存在的相邻孔隙的“协同充水”（cooperative pore filling）现象（Hughes and Blunt, 2000）、毛细阻断（snap off）现象（Berg et al., 2013）等等。这些孔隙尺度的排水/吸水现象上的差异是否会对达西尺度下的水分分布产生影响，文献中尚存在争议（Smiles et al., 1981; Ferrari and Lunati, 2014）。

本研究使用光透射法监测了二维砂箱中含水量的分布，并定性和定量地分析了排水 吸水过程中含水量分布的特征差异，目的是通过室内实验探讨孔隙尺度下的水/气驱替过程对于达西尺度下水分运移过程的可能影响，并加深人们对于非饱和带水气驱替过程的理解。

2 实验设计

2.1 实验材料与装置

实验中使用的介质，为美国尤尼明公司（Unimin Corporation, Le Sueur, MN）生产的高纯石英砂，商品名称Accusand。这种石英砂的SiO2含量大于92% 且有机质含量很低，因此透光性能良好，适合用作光透法实验的介质。实验所用的40/50目石英砂，颗粒密度2.66 g/cm3，粒径介于40~50目（0.300~0.425 mm）之间，平均粒径（d50）为0.359 mm，均一系数1.21（郑菲，2015）。实验中使用的二维砂箱，如图1a所示，内部尺寸为55 cm×45 cm×2.5 cm。砂箱结构由两块钢化玻璃板夹一个中心铝框构成。玻璃板与中心铝框之间使用硅酮密封胶黏合，以防止流体泄漏。中心铝框的左右边界完全封闭，上边界设置了六个进气口与大气相连通，下边界设置了六个进/出水口与蠕动泵相连，以控制砂箱中的水位（图1c）。为了尽可能避免出现水/气流动的优势通道，需采用湿法装填砂箱：将石英砂逐层装填进砂箱，使用小锤轻击砂面压实，同时调整水位，以保证水面始终高于砂面1 cm左右（若水位过高，则石英砂在沉降过程中有可能发生分选，出现分层的现象）。在装填过程中记录石英砂的总质量，结合石英砂颗粒密度及砂箱总体积可以计算得到实验所装填砂箱的平均孔隙度。

图1 实验装置示意图Fig. 1 The experimental setup (a) and side/front sketch view of the experimental setup (b, c)

2.2 光透法监测系统

实验中使用光透法监测系统（图1b）记录砂箱中的含水量变化。光透法监测系统由灯箱、砂箱和CCD（Charge Coupled Device）相机组成。灯箱作为光源放在砂箱一侧，其发出的光透射过砂箱以后，被放置在另一侧的CCD相机所接收。为了避免外界光源对于实验结果的干扰，整个实验需要放置在暗箱中进行。实验所用的暗箱尺寸约为2.5 m（长）×1.0 m（宽）1.0 m（高）。灯箱作为暗箱中的唯一光源，放置在暗箱前端，距砂箱约20 cm。CCD相机（AP2E, ApogeeInstruments, Auburn,CA)放置在暗箱后端，距砂箱约1.8 m。相机通过数据线与计算机相连，使用相机控制软件MaxIm DL（Diffraction limited, Ottwa）实时传输并记录透过砂箱的光强变化情况。CCD相机所记录的砂箱图片分辨率为1200像素（长）*980像素（宽），即实际砂箱中每1 cm长度在监测照片中约对应22个像素。

2.3 实验步骤

2.3.1 连续排水 吸水循环

实验开始前，将砂箱处于完全饱水状态。因此在实验开始前，CCD相机所记录的光强值，可以视作为完全饱水下的光强值Iws。完全饱水时，水所占的体积与孔隙体积相等，因此相应的饱和体积含水量θws就等于孔隙体积与总体积之比。此处体积含水量θw定义为：

式中，θw为体积含水量；V水为水的体积；V孔隙为孔隙的体积；V颗粒为颗粒的体积。饱水光强记录完成后，打开蠕动泵，控制砂箱经历三次排水 吸水的循环。首先是初次排水，以10 mL/min的流量从砂箱底部向外抽水，历时150分钟。初次排水结束以后静置50分钟，避免由于边界条件突然逆转造成的砂箱内含水量波动。静置结束后进行第一次吸水，使用蠕动泵以10 mL/min的速度从砂箱底部向内注水，历时100分钟。吸水结束以后，静置50分钟以完成第一次排水 吸水循环。随后，以相同的步骤又进行了两次排水 吸水的循环，流量与历时设置见表1。实验中，CCD相机的拍照间隔设置为1分钟，每次拍照的曝光时间设置为6 s。从完全饱水/干燥的条件下开始排水/吸水，称为初次排水/吸水。

表1 排水 吸水循环实验条件Table 1 Experimental conditions of drainage and imbibition

2.3.2 不同流量下的吸水 排水实验

为了研究流速对于含水量分布特征的影响，参照前人室内及野外实验的流速设置（Lehmann et al., 1998; Cartwright, 2014; Shoushtari et al., 2018），在同一个砂箱中使用三种其他流量（5 mL/min，100 mL/min，400 mL/min）又做了三次吸水 排水实验（表2）。在每一次吸水 排水实验开始之前，首先将砂箱中的水尽可能排干并静置 50分钟。接着使用设定的流量从砂箱底部向内注水至饱和，随后静置 50分钟并使用同样的流量从砂箱底部向外排水，直到将砂箱排干。对于5 mL/min 流量下的实验，CCD 相机的拍照间隔设置为1分钟，每张照片的曝光时间设置为 6 s。而对于100 mL/min和400 mL/min流量下的实验，CCD 相机的拍照间隔设置为3 s，每张照片的曝光时间设置为1 s。

表2 不同流量下吸水 排水实验条件Table 2 Experimental conditions of drainage and imbibition at different flow rates

2.4 分析方法及精度评价

CCD相机拍摄到的光强数据可以使用Niemet和Selker（2001）提出的方法转换成饱和度数据，其基本原理在文献中已有阐述（郑菲，2015；章艳红等，2014），本文在此就不做介绍。为了评估该光透法监测系统的测量精度，可以对比实际流入/流出砂箱的水分体积与光透法监测到的砂箱内水分体积的变化，如图2。这里使用均方根误差（RMSE）来表征实验过程中测量体积与实际体积的平均误差，即

图2 排水吸水循环实验中，实际流入/流出砂箱的水体积（actual volume）与光透法监测到的水体积（measured volume）变化对比Fig. 2 Comparison of the actual volume with measured volume of inflow/outflow water

式中，N为测量次数；Vactu为实际流入/流出体积；Vmeas为光透法测得的砂箱内水分体积变化。计算得到的RMSE值见表3，可见光透法监测数据的整体误差约在5%附近，最大不超过7%。这部分误差为光透法监测系统的系统误差，来源可能有很多方面，比如由于不透明边界的阻挡，贴近边界部分的介质饱和度无法被准确测量。另外，由于石英砂是自下而上逐层填入砂箱，这就导致砂箱下部的石英砂比砂箱上部更加密实一点。砂箱上下密实度的差异也会给光透法的测量引入系统误差(在计算饱和度及含水量时我们假设整个砂箱是均质的)。除了系统误差，光透法的测量结果也会受到随机误差的影响。光透法的随机误差主要由电噪声、光子散粒噪声、砂颗粒颜色差异等因素造成（Azooz et al., 1996）。与系统误差不同，由于随机误差会在求和时互相抵消，因此不会对整体的水量平衡造成影响。随机误差会造成单点测量结果的波动，因此影响到光透法测量结果的单点精度。这里使用残余水区域以及完全饱水区域的含水量波动标准差 (±2%S.D.)来评估光透法测量结果的随机误差。在残余水区域以及完全饱水区域，含水量一般被认为是常数，分别等于θwr以及θws。因此这两个区域中含水量的波动可以被归因于光透法测量的随机误差。实验过程中，光透法测得残余水以及完全饱水区域的含水量波动标准差见表3。可以看到，在实验过程中随机误差造成的测量结果波动约在2%附近，最大不超过3%，因此其对测量结果的影响基本可以被忽略。

表3 光透法监测结果的均方根误差及随机误差Table 3 Comparison of the RMSEs and random errors during the experiments

3 实验结果与分析讨论

3.1 排水 吸水过程含水量分布对比

图3显示了排水 吸水循环实验中六个阶段的含水量分布曲面。从图3a（初次排水）中可以看到，在砂箱的不同高度处，含水量的分布呈现出不同的规律。根据含水量的分布规律，可以将曲面分为三个部分。在图中的I号区域（即湿润锋以上的区域），含水量值小于 0.10，区域内的含水量基本没有横向波动，曲面相对平滑。而当含水量介于 0.10~0.25之间时（区域II，即湿润锋所在的区域），含水量出现了比较明显的横向波动，曲面变得不再平滑，并能观察到“指状”的起伏。当含水量增加到0.25以上时（区域III，即湿润锋以下的区域），含水量的横向波动又变得不明显起来。在图3b、c（第二和第三次排水）的湿润锋上也可以观察到“指状”起伏；但是对于吸水过程，在图3b、d、f 中，类似的“指状”起伏却不是特别明显。排水过程和吸水过程中湿润锋上的这种含水量分布上的差异可以被归结于以下原因：排水过程是一个相对不稳定的过程，由于孔隙尺度下海恩斯跳跃以及临近孔隙干扰排水这些瞬变现象的存在，排水过程介质中含水量的波动会大于相对稳定的吸水过程。

图3 含水量分布曲面Fig. 3 The moisture content distribution for each period

非稳定的排水过程，即孔隙尺度上的海恩斯跳跃现象，很早就有人研究（Klimentos and Mccann,1990）。海恩斯跳跃一般发生在气体驱替水分从孔喉进入孔腹之时。在排水过程中，由于毛细力的作用，水气交界的弯液面往往不会均匀地流过孔隙，而是直接从所处的孔喉迅速地跃变到下一个孔喉处。这种跃变往往也伴随着毛细压力的瞬间下降（Reynolds et al., 2009），造成相邻孔隙里水分分布的不稳定波动。Rücker等（2009）指出，在达西尺度下的排水过程，可以视作是一系列孔隙尺度下“雪崩式”不停发生的海恩斯跳跃事件的集合。从图3中排水与吸水过程湿润锋的对比可以推测，孔隙尺度下的不稳定排水过程影响到了湿润锋上的含水量分布，并且这类影响能够在达西尺度下被观察到。

为了更进一步研究排水与吸水过程中含水量波动的差异，将湿润锋上含水量的横向分布提取出来加以对比。图4显示的是第二次排水与吸水过程中，湿润锋上含水量横向分布的对比。从图4可以看出来，不管是排水过程还是吸水过程，其含水量的横向分布都存在一定的波动。在排水过程中含水量的波动范围更大（0.145~0.175），而吸水过程中含水量的波动范围（0.150~0.165）仅为排水过程的 1/2。排水过程中湿润锋上含水量的大范围波动，可以视作是导致图3中“指状”起伏出现的直接因素。

图4 第二次排水与吸水过程中，湿润锋上含水量横向分布的对比Fig. 4 Comparison of the lateral moisture distribution during main drainage and main imbibition. The insets show the locations in the sand chamber where the moisture content data is extracted

3.2 含水量空间分布相关性分析

图5中显示的是第二、第三次排水/吸水过程中湿润锋上含水量横向分布的对比。由图5a可以发现，在第二次排水与第三次排水过程中，含水量波动的趋势非常相似。同样的，在图5b中，第二次吸水与第三次吸水的含水量波动趋势也比较一致。根据图5中的数据，计算了这四个阶段中湿润锋上含水量横向分布的斯皮尔曼相关系数r （表4）。从表4可知，第二与第三次排水过程的含水量分布显示出极显著的强相关性（r=0.71，p<0.01），同样的，第二与第三次吸水过程的含水量分布也显示极显著的强相关性（r=0.79，p<0.01）。

表4 不同排水/吸水过程，湿润锋上含水量分布的斯皮尔曼相关系数Table 4 Spearman correlation coefficients (r) for the moisture distribution during drainage and imbibition.

图5 (a) 第二、第三次排水过程及 (b) 第二、第三次吸水过程中，湿润锋上含水量分布对比Fig. 5 Comparison of the moisture distribution on the wetting front at the 50th min of (a) main drainage and secondary drainage,(b) main imbibition and secondary imbibition

图5和表4中不同排水(或吸水)过程含水量分布的强相关性可以被归因于石英砂介质本身的空间分布特征。Suekane等（2009）使用X射线层析成像法研究了砂岩岩芯以及玻璃珠介质中压实程度与含水量分布的关系，并发现在平衡状态下压实度越大的地方含水量越大。Zhou等（2010）在随后的研究中指出，平衡状态下的含水量分布与介质孔隙度分布呈现较强的相关性，并且这种相关性不随实验流速的改变而改变。虽然本次实验砂箱内的局部孔隙度分布难以测量，但是从不同排水（或吸水）过程含水量分布的强相关性可以推测，实验中观察到的含水量波动并不是随机出现的，而是受到了介质本身的空间分布特征（局部孔隙度分布/局部压实度差异）的控制。此外，从表4中还可以看到，排水过程含水量分布与吸水过程含水量分布之间的相关性较弱，说明含水量的分布不仅与介质本身的空间分布特征有关，还会受到水分运动方向的影响。

3.3 含水量分布的Hurst指数

为了检验砂箱中含水量的横向分布是否符合分数布朗运动（fractional Brownian motion，fBm）模型分布，在砂箱的不同高度处提取出含水量的横向分布，并分别使用差值法和V/S分析法作双对数散点图。图6显示了三个不同位置的分析示例。从图6中可以发现，不管是差值法还是V/S分析法得到的散点图中，窗口长度（∆x）的对数与该窗口上统计量（|∆θw|或V/S2）的对数均呈较好的线性关系，并且两种方法估计得到的H值非常接近。由此可知，fBm可以很好地描述砂箱中含水量的横向波动。

图6 (a) 差值法及 (b) V/S分析法估计的Hurst指数值的对比Fig. 6 Comparison of Hurst coefficients estimated by (a) Intensity Difference method and (b) V/S analyses

图7显示了每个排水/吸水过程，砂箱中不同高度处含水量横向分布的Hurst指数值。由图7可见，不管是排水还是吸水过程，Hurst指数的峰值都出现在湿润锋上，排水过程的峰值H在0.5附近，而吸水过程的峰值H只有0.2左右。H值越小意味着含水量更倾向于局限在均值附近波动，因此波动范围越小；相对的，更大的H值意味着含水量的波动范围也更大。H值在湿润锋上达到峰值，意味着湿润锋上的含水量分布波动范围最大，这一点与我们在图3中观察到的现象一致。另外，吸水过程中0.2的H峰值显示了在吸水过程中含水量的波动受到了较强的抑制，其分布有强烈的回归均值的趋势；而排水过程中0.5的H峰值则说明湿润锋上含水量的波动更接近随机波动，相比吸水过程受到的抑制较少，更接近一个不稳定的过程。Moebius和Or（2012）在他们的孔隙尺度排水实验中指出，达西尺度下看似连续运动的湿润锋实际上在孔隙尺度上看来是非常不规则和不稳定的。从图7中可以看到，使用fBm模型及Hurst指数可以较好刻画这种不稳定排水现象造成湿润锋上含水量波动的精细结构。

图7 排水 吸水过程中砂箱不同高度处含水量分布的Hurst指数值Fig. 7 Comparison of Hurst coefficients in the sand chamber during each period

3.4 不同流量条件下含水量分布的Hurst指数

前文提到，较大的Hurst指数与孔隙尺度不稳定排水过程有关。因此，如果增大排水流量，加快排水速度，应该也能观察到 Hurst 指数的升高。图8中显示了不同流量条件下的Hurst指数对比。从图8中可以看到，排水过程中随着流量从5 mL/min升高到400 mL/min，湿润锋区域（区域II）的H值变化很小，而湿润锋以上区域（区域I）的H值变化较为明显，尤其是当流量增加到400 mL/min时，区域I的H值从0升高到了0.3附近。区域II的H值变化不明显是因为该区域H峰值本来就已经很高，在0.5左右的H值意味着含水量的波动已经接近随机波动，因此排水流量的进一步增加无法造成区域II中H值的增加。而区域I中H值的明显升高可以被归因于快速排水造成的不稳定现象。随着排水流量的增加，越来越多的水分被滞留于湿润锋上方来不及及时排出。由于这部分滞留的水分无法在湿润锋上方稳定存在，会随着排水过程的进行逐步排出。正是这些仍然在缓慢排出的滞留水分造成了区域I中Hurst指数的升高。另一方面，观察吸水过程的Hurst指数可以发现，随着吸水流量增大，吸水速度加快，计算得到的Hurst指数反而出现了些许下降。从图8中可以看到，随着吸水流量从5 mL/min上升到400 mL/min，砂箱中Hurst指数的峰值从0.2左右下降到了0.1左右。H峰值的下降意味着随着吸水流量升高，湿润锋上含水量的波动有更强烈的均值回归的倾向，其波动范围受到了抑制。关于吸水速度对于湿润锋形态的影响，前人也做了不少研究（Berg et al., 2013, Nguyen et al.,2006）。Hughes和Blunt指出，在吸水速度较快的时候，粘性力占主导而毛细力受到抑制，因此水分更倾向于优先充填临近的孔隙，而不是在毛细力的作用下上升形成指状毛细水（Berg et al., 2013）。此时相邻孔隙的“协同充水”现象发生得很频繁，在孔隙尺度下来看，水气交界面往往是连续且平滑的。而当吸水速度较低时，毛细力占据了主导地位，在毛细力的作用下水分沿着孔隙直径较小的路径上升，形成了一处处不规则的指状毛细水。根据图8中的结果可以推测，这种吸水速度对于水气交界面的影响同样能在湿润锋上体现出来。由于介质本身存在一定的非均质性，低流量下湿润锋上的水分在毛细力的作用下不规则上升，体现为含水量在横向上的波动。然而随着吸水速度加快，毛细力的作用受到抑制，湿润锋上的水分分布变得更加规则，所以此时含水量的波动也受到了抑制，H值也随之下降。

图8 不同流量条件下吸水/排水过程含水量分布Hurst指数对比Fig. 8 Comparison of Hurst coefficients under different flow rates

4 结论

通过二维砂箱实验研究了排水及吸水过程中含水量分布的差异，并使用分数布朗运动模型定量分析和讨论了含水量的波动特征。在排水过程中，由于“海恩斯跳跃”等孔隙尺度不稳定排水现象的存在，导致含水量在横向上出现较大范围的波动，并且在含水量曲面上能观察到含水量波动造成的指状起伏。而在吸水过程中，含水量的波动范围更小，含水量曲面也更为平整。含水量分布在排水过程和吸水过程之间的差异说明了孔隙尺度的一些常见过程能够对达西尺度下的水分运动过程造成影响，而不是在提升尺度的过程中被简单抵消掉。