包气带作为连接大气圈、地表水圈和饱水带的重要关系纽带，其中的水分运移是整个水文循环过程中的关键环节。众多研究表明，包气带中的水分运移过程本质上是一个水气两相相互驱替的过程。在我国北部、西北地区，包气带厚度通常较大，且岩性结构较为复杂，这会对水气两相运移和污染物下渗等过程产生较大的影响。但在相关研究和一些实际工程项目如地下水环境影响评价、污染场地调查修复和污废水泄露模拟预测中，多忽略包气带的影响，采用简化的单相流模型，少数采用多相流模型的研究也是基于均质条件^[1-7]。对于野外常见的层状非均质岩性结构、不同入渗或泄露速率的情况，这种简化可能会引起较大的误差。

包气带中的水分运移涉及非常复杂的理化过程^[8-12]，而岩性结构是影响包气带水分运移、潜水蒸发以及地下水补给的重要因素之一。Colman等^[1]提出，在层状非均质结构条件下，水分的运移主要受细质土控制。王文焰等^[2]研究认为砂层在黄土中具有良好的阻水性和减渗性。解文艳等^[3]通过不同质地土壤的有压入渗试验发现，土壤质地由轻变重时，土壤入渗能力减小。张建丰^[4]研究了不同层状结构条件下土壤水分的入渗特性。郭会荣等^[5]通过土壤水流穿透试验指出，土壤水在粉细砂中通过活塞式入渗补给地下水，在亚黏土中则通过优势入渗补给地下水。王仕琴等^[6]以氢氧同位素为“示踪剂”，研究了不同土壤质地的降雨入渗过程。霍思远等^[7]利用Hydrus模拟了不同情况下包气带中黏土透镜体对降雨入渗补给饱和带地下水的影响。以上研究均探讨了岩性结构对包气带水分运移的影响，但忽略了入渗过程中包气带的气压变化，在两相流模型中，没有考虑层状非均质岩性结构对气体和水分运移的影响。

入渗速率方面，张杰等^[13]通过分析Hydrus-1D模拟结果，认为饱和导水率越大，入渗速率与变雨强时变曲线重合率越高。朱伟等^[14]通过土柱试验，探讨了土体初始水分状态等因素对入渗速率的影响。刘刚等^[15]对比了通气、封气条件下黏性土中水分入渗模拟试验数据，认为入渗速率、含水率和孔隙气压力在入渗过程中是协调变化的。刘秀花等^[16]综合畦灌试验数据和Hydrus-1D模拟结果发现，进水流量和历时是影响包气带水分补给饱和地下水的重要因素。上述研究均是在均质岩性结构单相流模型、两相流模型或非均质岩性结构单相流模型下，对包气带水分入渗过程中入渗速率的影响与变化进行分析，未涉及入渗速率对层状非均质岩性结构两相流模型中水气运移的影响。

本文利用数值模拟程序TOUGH2建立了层状非均质岩性结构和不同入渗速率条件下的水分入渗单相流和水气二相流数值模型，对比研究不同层状非均质岩性结构下气相对包气带水分运移的影响，分析包气带水气驱替规律。并针对包气带“上粗下细”结构模型，进一步讨论入渗速率对两相流计算结果的影响。旨在通过本文研究，拓展包气带水分运移机理及包气带对地下水补给的影响研究，丰富包气带水气二相流影响因素研究，同时为实际工程中饱和-非饱和渗流数值模型的选择提供决策依据。

1.   概念模型

由于施工开挖或填土、地面夯实或硬化、水平排水等工程措施和地层的天然沉积特征，实际污染场地、填埋场等的地层结构通常都具有水平层状非均质特征，相应设置模拟区域为垂向二维层状非均质剖面，下部饱和带岩性固定为粗砂，上部包气带分上、下两层且岩性不同，模型尺寸为100 m×60 m（60 m为垂向厚度），地下水埋深为40 m，即包气带厚度为40 m，饱水带厚度为20 m。

为探讨包气带层状非均质岩性结构在单相流和两相流模型中对水分和气体运移的影响，固定包气带上层岩性为细砂，包气带下层岩性设置2 种情景：粗砂和亚黏土，分别模拟水在包气带中由低渗介质入渗至高渗介质（以下简称“上细下粗”岩性结构）和由高渗介质入渗至低渗介质（以下简称“上粗下细”岩性结构）的情形，见图1。介质参数取值参考相关数值模拟软件数据库和SUN等^[17]关于入渗引发准饱和带的研究，见表1。

图  1  包气带层状非均质岩性结构概念模型

Figure  1.  Conceptual model of the layered heterogeneous lithological structure in the vadose zone

下载: 全尺寸图片 幻灯片

表  1  非均质介质主要参数^[14]

Table  1.  Main hydrogeological parameters of heterogeneous media^[14]

介质	密度/（kg·m−3）	渗透率/（10−12m2）	孔隙度	相对渗透率-饱和度关系曲线 (VG-M模型)					毛细压力-饱和度关系曲线 (VG模型)
				λ	Slr	Sls	Sgr		λ	Slr	Sls	Pmax	1/P0
粗砂	2 650	30	0.31	0.416	0.06	1	0.05		0.416	0.06	1	5×106	3.0×10−4
细砂	2 650	8.0	0.35	0.430	0.10	1	0.05		0.430	0.10	1	5×106	2.5×10−4
亚黏土	2 650	0.16	0.40	0.447	0.15	1	0.05		0.447	0.15	1	5×106	1.0×10−4
注：λ为经验参数，${S}_{\mathrm{l}\mathrm{r}}$为残余水饱和度，${S}_{\mathrm{l}\mathrm{s}}$为饱和水饱和度，${S}_{\mathrm{g}\mathrm{r}}$为残余气饱和度，Pmax为最大吸力，1/${P}_{0}$为土壤进气值的倒数。

下载: 导出CSV 

| 显示表格

从二十世纪中叶至今，许多研究建立了毛细压力—饱和度关系曲线模型，如VG模型、BC模型、Gardner模型等，本次模拟选取较为常用的VG模型描述毛细压力与饱和度的关系。同时，采用VG-M模型描述相对渗透率与饱和度的关系。

网格剖分垂直与水平间距均为1 m，顶部为固定大气压力边界（1.013×10⁵ Pa），底部为固定压力边界（2.975×10⁵ Pa），左右两侧为无通量边界。在模型顶部第2 层网格中央设置30个网格长度（30 m，入渗范围）作为入渗补给源，模拟污染场地或填埋场中污废水面状泄露的情景，入渗速率设置为4.0×10⁻⁵ m/s。模型不考虑温度变化，设定为恒温25 ℃。

2.   研究方法

TOUGH2是模拟多相和多组分流体流动及热量运移的数值模拟程序^[18]，通过积分有限差分法（IFDM）对控制方程进行空间离散，包含多个状态方程（equation of state，简称EOS）模块，其中EOS9模块主要基于Richards方程模拟饱和-非饱和水流（即单相流）问题，因此其本质上只能模拟水分的运动，即入渗过程中水分运移到的区域气体会被自然排走。EOS3模块主要基于多组分质量和能量守恒方程模拟等温或非等温条件下的水气二相流，即可以模拟包气带入渗过程中水分与气体的运动及两者间的相互作用，因此与EOS9模块在入渗过程的模拟上存在核心差异^[19]。因为基于Richards方程的单相流模型相对为人所熟知，本文仅对EOS3模块水气二相流模型进行简单介绍。

任何单元体内水气二相流满足的质量和能量守恒方程可写成积分形式^[19]：

式中：${V}_{n}$——系统中的任意单元体；

${M}^{\kappa }$——各组分每单位体积的质能累计量；

${F}^{\kappa }$——各组分的质能通量；

${\Gamma }_{n}$——单元体${V}_{n}$的封闭边界；

${\boldsymbol{n}}$——指向单元体内部的单位法向量；

$\kappa$——物质组分（例如水、空气等）；

${q}^{\kappa }$——组分$\kappa$的源汇项；

t——时间。

所有组分质量表示为：

式中：$\phi$——孔隙度；

β——气相或液相，$\beta =1,2,\cdots ,{N}_{\mathrm{P}\mathrm{H}}$；

${S}_{\beta }$——各个相的饱和度；

${\rho }_{\beta }$——各个相的密度；

${X}_{\beta }^{\kappa }$——各个相组分的质量分数。

${F}^{\kappa }$的对流通量为各相中该组分的质量流量之和：

${F}^{\kappa }$的扩散通量可表示为：

式中：${\bar{{\boldsymbol{D}}}}_{{\boldsymbol{\beta}} }^{{\boldsymbol{\kappa}} }$——水动力扩散张量。

基于EOS9、EOS3模块分别构建层状非均质岩性结构条件下水分入渗单相流和水气二相流数值模型，通过模拟结果对比分析气相对包气带水分运移的影响，探究包气带水气驱替规律。

3.   水分入渗单相流与水气二相流模拟结果对比与讨论

通过模拟得到了不同的层状非均质岩性结构下单相流和水气二相流的水分入渗状况，分析了不同岩性结构下水分入渗单相流与水气二相流的特点和区别，以及“上粗下细”岩性结构下包气带渗透性和入渗速率对单相流与水气二相流入渗过程的影响。

3.1   岩性结构

3.1.1   “上细下粗”岩性结构

“上细下粗”岩性结构模型模拟水分在包气带由低渗透介质向高渗透介质入渗，针对包气带水分入渗非稳定过程，考虑到我国北部、西北地区相对干旱的气候条件，假设污废水泄漏前场地未接受降水入渗补给，并参考相关研究^[17]的模型设置方式，以无入渗时毛细作用下的稳定水分分布作为初始条件，在实际应用时应对此加以注意。

稳定状态下单相流模型和水气二相流模型的水相饱和度和孔隙气压力分布随土层厚度变化基本一致：水相饱和度在潜水面以上整体变化较小，在包气带上、下层介质交界面处略有减小，从距离潜水面约5 m处，随深度增加，水相饱和度迅速增大，到达潜水面后，水相完全充满介质孔隙，水相饱和度恒定为1；孔隙气压力在包气带均为大气压（1.013×10⁵ Pa），在饱水带随深度线性增加至底部边界设定压力（2.975×10⁵ Pa）。

“上细下粗”岩性结构水分入渗非稳定过程的模拟时长为80 d，单相流与水气二相流模型计算得到水分入渗52 d水相饱和度分布见图2。二者水相饱和度的分布基本一致，入渗的水在穿过包气带双层介质界面时，由于下层粗砂的渗透性较强，水相饱和度明显降低。

图  2  “上细下粗”岩性结构条件下水分入渗52 d水相饱和度分布

Figure  2.  Distribution of water saturation on the 52th day of infiltration under the “upper fine and lower coarse” lithological structure

下载: 全尺寸图片 幻灯片

在水气二相流模型的包气带双层介质交界面上、下分别取2 个观测点D、E，见图1（a）。D、E孔隙气压力与水相饱和度随时间的变化曲线见图3。 D、E孔隙气压力仅在水分入渗初始阶段达到最大峰值，在湿润锋推进阶段波动幅度较小，气压变化对水流的阻碍可以忽略不计。计算结果表明，无论是单相流模型，还是水气二相流模型，入渗水均在56 d补给到潜水面，因此在“上细下粗”岩性结构下，可以忽略气相对水分入渗过程的影响。

图  3  “上细下粗”岩性结构交界面上、下观测点气压与水相饱和度变化曲线

Figure  3.  Variations of pressure and water saturation at points D and E under the “upper fine and lower coarse” lithological structure

下载: 全尺寸图片 幻灯片

结合水气二相流模型包气带气流分布模拟结果可知，水入渗至渗透性较强的粗砂层后，水相饱和度降低，气流通道变宽，气相流动变快，湿润锋前端气体受到压缩后，在侧向逸出的同时还会形成向上的气体出流通道，被压缩的气体迅速得到释放。故粗砂层内的观测点E气压变化幅度较小，气相对水分入渗的阻滞较弱，可以近似忽略。因此在“上细下粗”岩性结构条件下，基本可以直接采用单相流模型对包气带水分运移进行模拟。

3.1.2   “上粗下细”岩性结构

在“上粗下细”岩性结构条件下，先计算一个稳定场作为模型初始条件。与“上细下粗”岩性结构类似，稳定状态下单相流与水气二相流模型计算得到的水相饱和度和孔隙气压力分布随土层厚度的变化基本一致，由于介质差异，2 种模型计算得到的水相饱和度在上、下层介质交界面处均存在拐点。

“上粗下细”岩性结构下水分入渗非稳定渗流过程模拟时长为100 d，单相流与水气二相流模型计算得到80 d水相饱和度分布见图4。两者在下层低渗透介质中均出现水相近饱和区，但二相流模型湿润锋明显滞后于单相流，且二相流计算介质交界面处的水分有明显的水平运移，说明气相对入渗水起到了明显的阻滞作用。

图  4  “上粗下细”岩性结构条件下水分入渗80 d水相饱和度分布

Figure  4.  Distribution of water saturation on the 80th day of infiltration under the “upper coarse and lower fine” lithological structure

下载: 全尺寸图片 幻灯片

入渗水由包气带细砂层进入包气带亚黏土层时，由于亚黏土层的渗透性比上覆介质差，对水相的运移起到一定的阻滞作用，导致水相饱和度增大并接近饱和。为进一步说明气相对水流入渗的阻滞作用，在水气二相流模型的包气带双层介质交界面上、下分别取2 个虚拟观测点D′（50,38）、E′（50,31）见图1（b），得到D′、E′孔隙气压力与水相饱和度随时间的变化曲线，见图5。

图  5  “上粗下细”岩性结构交界面上、下观测点气压与水相饱和度变化曲线

Figure  5.  Variations of pressure and water saturation at points D′ and E′ under the “upper coarse and lower fine” lithological structure

下载: 全尺寸图片 幻灯片

从图5中可以看出，双层介质交界面下方E'点孔隙气压力在湿润锋到达前后开始增大，最大增加31677 Pa。相较而言，双层介质交界面上方D'点孔隙气压力波动无论是在水分入渗初始阶段还是湿润锋推进阶段，都显得微不足道。亚黏土层的气压变化已经明显阻滞了水相的入渗，并形成顶托，致使入渗水流在交界面处发生明显的侧向运移。计算结果表明，单相流水分入渗模型中水分在80 d补给到潜水面，而二相流水分入渗模型在98 d补给到潜水面。由此可见，此时忽略气相对水流的阻滞作用会造成较大的误差，需要采用水气二相流模型进行模拟。

水气二相流模型包气带气流分布模拟结果与“上细下粗”岩性结构的模拟结果存在明显不同：“上粗下细”岩性结构模型中，水分入渗至包气带双层介质界面后，由于亚黏土层的低渗透性，水相饱和度迅速增大，气流通道变窄，气相流动变慢，湿润锋前端气体受到压缩后，难以形成向上的出流通道，只能由侧向逸出，被压缩的气体得不到充分释放，因此黏土层湿润锋前端的孔隙气压力持续增大，从而阻挡了水相的入渗。同时，气压变化也会引起潜水面处水压的升高。

选取潜水面处F'得到孔隙水压力随时间变化曲线见图6。水分入渗至亚黏土层后，孔隙水压力升高明显，待水分入渗补给到潜水面，压力增加了4 283 Pa，相当于水头上升了0.44 m，在入渗水补给到潜水面之前，地下水水位由于气相的压缩已经明显升高，即发生了Lisse效应。Lisse效应指湿润锋前端气体受到压缩，致使观测井中水位上升的现象。Thal Larsen根据荷兰利塞村的水位数据最早描述了这种现象^[20]。

图  6  “上粗下细”岩性结构下水分入渗导致潜水面孔隙水压力变化历时曲线

Figure  6.  Variation of pore water pressure on the water table caused by infiltration under the “upper coarse and lower fine” lithological structure

下载: 全尺寸图片 幻灯片

3.2   包气带渗透性

以包气带上层介质渗透率为基准，通过保持其余参数不变只改变“上粗下细”岩性结构中包气带下层介质渗透率，增设5种模拟情景，进一步探讨包气带层状非均质岩性结构在单相流和两相流模型中对水分和气体运移的影响。模拟情景设置及相应参数取值见表2。需要说明的是，在进行不同岩性结构水分入渗单相流与水气二相流模拟时，相对渗透率函数和毛管压力函数均根据介质岩性进行了相应的设置。在研究包气带上、下层介质渗透性差异对单相流与水气二相流入渗过程的影响时，为了更加直观地揭示渗透性差异的影响，假设不同情景间仅渗透率存在区别^[21]（表2），VG-M和VG模型参数设置与“上粗下细”岩性结构水分入渗单相流与水气二相流研究保持一致。

表  2  不同模拟情景参数设置

Table  2.  Parameter settings of different simulation cases

模拟情景设置		渗透率/（10−12m2）	上下层介质 渗透率比值
上层介质（细砂）		8.0	/
情景1	下层介质1 （亚黏土）	0.16	50
情景2	下层介质2	0.20	40
情景3	下层介质3	0.40	20
情景4	下层介质4	0.50	16
情景5	下层介质5	0.80	10
情景6	下层介质6	1.6	5

下载: 导出CSV 

| 显示表格

通过对比各情景模拟结果可知，包气带下层介质渗透率的增加会导致下层介质内水分运移速率增大，使气相对水相的阻滞能力明显减弱，不同情景下单相流与水气二相流模型中水分入渗补给到潜水面的时间变化见图7。

图  7  不同情景下单相流与水气二相流模型水分入渗补给至潜水面所需时间变化曲线

Figure  7.  Time variation curve of infiltration recharge to the water table in the single-phase flow and water-gas two-phase flow models under different simulation cases

下载: 全尺寸图片 幻灯片

从图7中可以看出，当包气带上下层介质渗透率比值为50倍（情景1）时，单相流模型中水分80 d补给到潜水面，二相流模型中水分98 d补给到潜水面，两者相差约18 d。随着包气带下层介质渗透率的增加，这一差距逐渐缩小。当渗透率比值缩小至20倍（情景3）时，单相流和二相流模型中水分补给至潜水面所需时间分别为72，74 d，表明气相对水流的入渗仍有一定程度的阻滞作用。当下层介质渗透率为5×10⁻¹³ m²，上下层介质渗透率比值为16倍（情景4）时，单相流与水气二相流模型中下渗水70 d均补给到潜水面，且水相饱和度分布一致，即此时单相流与水气二相流模型已无明显区别。此后继续增大下层介质渗透率（情景5、6），两者模拟结果仍保持一致。

同时，包气带双层介质交界面下方E'点的气压波动幅度随上下层介质渗透率比值的减小而逐渐减小（图8），这是由于下层介质渗透率增加促进了水相和气相的流动，减缓了气体蓄积作用。因此，包气带下层介质的渗透性越差、上下层介质渗透率比值越大，气相对水流入渗的阻滞作用越明显，此时包气带水分运移模拟应采用水气两相流模型。

图  8  不同情景下“上粗下细”岩性结构界面下方观测点气压与水相饱和度变化曲线

Figure  8.  Variations of pressure and water saturation at point E′ under different simulation cases

下载: 全尺寸图片 幻灯片

在模型设定的底部压力保持不变及污废水泄漏前场地未接受降水入渗补给等条件下，当包气带上下层介质渗透率比值小于16倍时，气相对水流的阻滞作用可以忽略，此时可以采用单相流模型进行模拟。

3.3   入渗速率

入渗速率的改变通常会影响包气带水分运移速率进而导致包气带的气压变化^[22]，为进一步探讨包气带“上粗下细”岩性结构水分入渗模型选择，比较不同入渗速率对水气两相流计算结果的影响，在“上粗下细”岩性结构水分入渗单相流与水气二相流模型（下层介质渗透率k=1.5×10⁻¹³ m²，入渗速率q=4×10⁻⁵ m/s）的基础上增设入渗速率为2×10⁻⁵ ，8×10⁻⁶ m/s的对照组。模拟结果表明，随着入渗速率的减小，水分补给到潜水面所需时间逐渐变长。

经过计算，在q=2×10⁻⁵ m/s时，单相流模型中水分在143 d补给到潜水面，而相应的水气二相流模型水分在150 d补给到潜水面，表明包气带黏土层中气相对水流仍有一定的阻滞作用。图9为水分入渗130 d水相饱和度的分布，可以看到单相流和两相流模拟结果仍有差异。

图  9  入渗速率q=2×10⁻⁵ m/s时，“上粗下细”岩性结构水分入渗130 d水相饱和度分布

Figure  9.  Distribution of water saturation on the 130th day of infiltration under the “upper coarse and lower fine” lithological structure when q=2×10⁻⁵ m/s

下载: 全尺寸图片 幻灯片

当q=8×10⁻⁶ m/s时，单相流与水气二相流模型计算结果基本一致，水分均在291 d补给到潜水面。这说明随着入渗速率的减小，气相对水流入渗的阻滞作用逐渐减弱。

对于水气二相流模型，双层介质交界面下方E′点处水相饱和度和孔隙气压力随时间的变化曲线见图10。通过对比可以发现，入渗速率越小，包气带水分运移速率随之变小，气体被压缩后能够在短时间内得到释放，气压增幅较小，对水流的阻滞能力减弱。因此，入渗速率越大，气相对水流入渗的阻滞作用越明显，此时包气带水分运移模拟应采用水气两相流模型。经模型参数敏感性分析，在本模型设定条件下，当q<1×10⁻⁵ m/s时，气相的影响可以忽略，此时可以采用单相流模型进行模拟。

图  10  入渗速率q=2×10⁻⁵ ，8×10⁻⁶ m/s时，“上粗下细”岩性结构界面下方观测点气压与水相饱和度变化曲线

Figure  10.  Variations of pressure and water saturation at point E′ under the “upper coarse and lower fine” lithological structure when q=2×10⁻⁵ m/s and 8×10⁻⁶ m/s

下载: 全尺寸图片 幻灯片

4.   结论

本文利用TOUGH2模拟对比了不同包气带层状非均质岩性结构和入渗速率对包气带单相流与水气二相流水分入渗过程的影响，选取VG模型描述毛细压力与饱和度关系以及设定的模型初始条件和底部边界等条件下，得到以下结论：

（1）在包气带“上细下粗”岩性结构条件下，水分入渗单相流与水气二相流模拟结果基本一致，包气带气体被压缩后能够迅速释放。从宏观上看，气相的影响可以忽略，可直接采用单相流模型对包气带水分运移进行模拟。

（2）在包气带“上粗下细”岩性结构条件下，气相对水分运移的阻滞作用显著，在潜水面出现明显的Lisse效应。且包气带下层介质的渗透率越小、上下层介质渗透率比值越大，气相对水流入渗的阻滞作用越明显，单相流与两相流的运移结果差别越大。在本文模型设定的底部压力保持不变及污废水泄漏前场地未接受降水入渗补给等条件下，包气带上下层介质渗透率比值大于16倍左右时，忽略气相影响会造成较大的误差，需要采用水气二相流模型模拟包气带水分运移。

（3）入渗速率对包气带气压影响较大，在包气带“上粗下细”岩性结构条件下，入渗速率越大，气相对水流入渗的阻滞作用越明显，此时包气带水分运移模拟应采用水气两相流模型。