饱和土壤中地下水渗流过程的研究在诸如地下水评价、水土污染防治等多个领域具有重要意义^[1-2]。地下水渗流过程特性通常是土壤等含水多孔介质属性的一种表征，并对流体中污染物运移行为有着决定性的影响，是联结介质属性与溶质运移的桥梁。含水介质中的很多溶质运移现象均可通过分析地下水渗流特性来探索其内在机制^[3]。如Bijeljic等^[4]的研究表明，流速的概率密度分布差异是导致具有不同复杂程度的碳酸盐岩介质中溶质不同运移特征的内在原因。Dou等^[5]则发现流场的非均质性增强会引起溶质的非菲克（non-Fickian）运移特征显著增强。地下水在土壤中渗流时常分布有大量的不流动区及优势流区，二者的空间分布及其在渗流区中所占比例差异对地下水中溶质运移具有显著影响，是地下水渗流特性的重要研究内容之一^[6-8]。

当分布有大量不流动区时，通常导致溶质滞留于其中，运移出相应区域的时间显著增加，最终引起拖尾现象的增强，这是非菲克运移的重要特征之一^[9]。而在地下水污染修复中最关键的问题就是何时能清除污染区域中的污染物，因此该类实践工程通常要充分考虑不流动区的影响。此外，不流动区的空间分布、大小直接决定了溶质滞留于其中的时间，进而与两区模型（TRM）的溶质质量交换系数、时间分数阶对流弥散方程（tFADE）的分数阶导数、连续时间随机行走模型（CTRW）的溶质等待时间概率分布等溶质运移模型参数密切相关^[10-11]。上述参数是相应模型的核心，对模型的模拟精度至关重要。通过分析地下水中不流动区分布特征及变化趋势，有助于推进不流动区与流动区之间质量交换的研究，进而为确定相应参数函数形式提供研究基础^[12-13]。

优势流区域通常引起流线的聚敛。流线聚敛就意味着周围的流体在短距离内迅速汇集流入渗透性高于其周围土壤介质的区域，即发生流体聚流现象^[14]。这对饱和土壤中溶质运移影响显著，主要表现在两方面。首先，由于优势流区流体流速明显高于流体平均流速，溶质进入优势渗流通道快速运移，通常导致早穿透现象增强，即首次穿透时间提前，这对于垃圾填埋场防护等至关重要。因此，在以上工程实践中必须考虑优势渗流区的大小、空间分布等特征。另一方面，溶质随地下水进入聚流区，溶质穿透给定数目的流线所需运移的横向距离明显缩短。而不同流线间的质量交换正是溶质稀释和不同溶质间混合的唯一机制^[15-16]。因此，优势流区的多少及其空间分布、流经优势流区的流量比例与溶质稀释作用、混合作用的强度密切相关，并最终影响非反应性污染物的浓度衰减及反应性污染物的反应速度与强度^[14，17]。如Rolle等^[18]的研究表明，优势渗流通道引起的聚流现象对溶质横向混合作用的增强效应远超过流场瞬变。Willingham等^[17]发现当存在优势渗流通道时，溶质反应速率的增强高达40%，且这种增强与优势渗流区的长度及优势渗流区与周围介质的渗透率差异程度相关。

土壤等多孔介质孔隙尺度结构对地下水渗流过程具有显著影响^[19-21]。例如含水介质孔隙尺度表面粗糙度增大会增加地下水流动中的水力迂曲度，导致介质渗透率降低。张学羿和窦智^[1]研究表明，黏土的孔隙尺度结构差异是引起其中地下水流速空间分布非均质的主要因素。魏赫鑫等^[22]的研究表明岩石比表面积、微孔孔径、孔容等介质微观结构特征对介质中流体分布具有显著影响。惠威等^[23]发现孔隙喉道的大小、分布与介质中流体赋存、渗流性能密切相关。多孔介质孔隙尺度结构差异引起的地下水渗流特性变化会进一步影响宏观尺度上地下水渗流及地下水中的溶质运移过程^[24-25]。然而，以往针对土壤等多孔介质中地下水流动的研究多是将其概化为连续介质后进行研究。这种方法降低了研究的复杂程度，却不能准确考虑孔隙尺度地下水渗流过程、流速分布差异及介质微观结构特征对其影响，而这正是深入理解宏观尺度上渗流、溶质运移现象的关键。近年来，计算机技术及微观成像技术快速发展，获取、构建土壤等含水介质微观结构变得较为容易，孔隙尺度模拟为在孔隙尺度上深入研究多孔介质中地下水渗流动力学、溶质运移提供了一个有力的研究手段^[26-28]。

胶结度是土壤等多孔介质的基本属性，固体颗粒不同程度胶结是饱和多孔介质中最常见的现象之一，如土壤颗粒聚集、岩石压实等^[29-31]。固体颗粒胶结度增大时，常引起多孔介质的孔隙结构显著变化。如导致介质内部孔喉结构减少，孔隙表面积、迂曲度增加，介质非均质性增强等。但是，目前关于胶结度不同时土壤等多孔介质孔隙尺度结构差异对地下水渗流影响的研究尚较少。本文基于孔隙尺度模拟方法，研究胶结度不同的饱和土壤中地下水渗流特性差异；深入分析胶结度不同时，饱和土壤中地下水流速空间非均质性程度及流速概率密度分布，不流动区、优势流区的空间分布及比例差异，并探索相应渗流特性变化的内在机理。

1 材料与方法 1.1 土壤多孔介质的构建方法

目前，在孔隙尺度上重构多孔介质孔隙结构的方法主要包括扫描重构、等价网络模型以及模拟生成三类^[32-33]。其中，模拟生成具有能够较好反映介质孔隙结构、效率高、成本低等优点，被广泛应用于相关研究领域。

本文采用前期研究中基于整体重排算法编制的程序构建不同胶结度的土壤介质，该方法已被广泛应用于多孔介质中地下水流动、溶质运移等领域，能够反映多孔介质的孔隙结构特征^[34]。该方法中土壤颗粒由圆形粒子表示，构建不同胶结度的多孔介质主要包括三个步骤。首先，在给定区域内持续生成固体颗粒并随机确定其初始位置，直至初始孔隙率达到给定值。自然界中土壤颗粒粒径多数符合截断的对数正态分布^[35]。因此，生成过程中，固体颗粒粒径由截断的对数正态分布函数随机确定。然后，通过循环的搜索、重排相互重叠的固体颗粒来消除所有固体颗粒间的重叠，从而构建固体颗粒间完全不胶结的土壤介质。最后，基于完全不胶结的土壤介质，保持所有土壤颗粒位置不变，粒径增加相同的值，直至部分颗粒开始相互重叠，即土壤固体颗粒开始胶结。持续的增大固体颗粒粒径，从而构建固体颗粒不同程度胶结的土壤介质微观结构。胶结度由胶结的固体颗粒数（即与其他颗粒重叠的固体颗粒）占颗粒总数的比例（Pc）定量衡量。构建介质程序的详细信息可参考相关文献^[34-35]。

本文构建土壤多孔介质的区域范围为80 mm × 60 mm（x × y），控制固体颗粒粒径的截断对数正态分布函数的期望及方差分别为1.12 mm和0.25 mm。为防止生成过大或过小的固体颗粒，颗粒半径被限制在0.55 mm和1.80 mm之间。初始孔隙率设置为0.40。本文共构建了四组不同胶结度的土壤介质，孔隙率由大到小分别标记为PM40、PM35、PM30、PM25（分别表示孔隙度为0.40、0.35、0.30、0.25的多孔介质）。构建的介质如图 1所示，其中PM40为固体颗粒完全不胶结的土壤介质。本文统计了不同介质中胶结固体颗粒比例Pc，结果表明，从PM40到PM25，随着孔隙度从0.40减小至0.25，Pc从0增至60.20%。表明从PM40到PM25固体颗粒胶结度显著增大。土壤介质的属性如表 1所示，包括孔隙率n、固体颗粒平均半径r_ave、Pc等。

1.2 渗流模拟与数据统计分析

本文运用有限元模拟软件COMSOL Multiphysics 5.1中的层流模块求解Navier-Stokes方程，模拟多孔介质中孔隙部分中地下水流动：

$ \rho \nabla \cdot \left( {\bf{u}} \right) = 0 $

(1)

$ \rho \left( {{\bf{u}} \cdot \nabla } \right){\bf{u}} = - \nabla {{p}} + \mu {\nabla ^2}{\bf{u}} $

(2)

式中，u为地下水的流速矢量，m·d^–1；p为压强向量，Pa；为水的黏度，Pa·s；为水的密度，kg·m^–3。

模型中地下水在水力梯度的作用下，由左向右流动。左边界、右边界设置为第一类边界，介质的上边界、下边界、内部流固边界设置为不流动边界，如图 1所示。通过调整左边界、右边界的水头差使地下水平均流速达到目标值。本文主要模拟了平均流速分别为0.1 m·d^–1、1 m·d^–1、10 m·d^–1的三种地下水典型流速，相应渗流模型入口处压强如表 2所示。图 1为平均流速为1 m·d^–1时模拟得到的流速空间分布。

为定量分析地下水流速空间分布非均质性的变化，本文计算了地下水流速u及其沿地下水主要流动方向、垂直于地下水主要流动方向上的流速矢量u_x、u_y的变异系数CV_u、CV_ux、CV_uy，变异系数越大表明流速的非均质越强，计算公式为：

$ {\text{C}}{{\text{V}}_{\text{u}}} = \frac{{{\sigma _{\text{u}}}}}{{{u_{{\text{ave}}}}}} = \frac{{\sqrt {{{\int_{{V_{\text{P}}}} {{{\left( {u - {u_{{\text{ave}}}}} \right)}^2}d{V_{\text{p}}}} } \mathord{\left/ {\vphantom {{\int_{{V_{\text{P}}}} {{{\left( {u - {u_{{\text{ave}}}}} \right)}^2}d{V_{\text{p}}}} } {{V_{\text{P}}}}}} \right. } {{V_{\text{P}}}}}} }}{{{{\int_{{V_{\text{P}}}} {ud{V_{\text{p}}}} } \mathord{\left/ {\vphantom {{\int_{{V_{\text{P}}}} {ud{V_{\text{p}}}} } {{V_{\text{P}}}}}} \right. } {{V_{\text{P}}}}}}} $

(3)

式中，u为地下水流速，m·d^–1；u_ave、_u分别为流速的平均值和标准差，单位均为m·d^–1。V_P为被流体占据的孔隙区域的面积，单位为m²，随着胶结度增大，V_P显著减小，且孔隙区域V_P与固体颗粒间的流固界面不规则程度显著增加。

为定量研究不同介质中流速分布的差异，本文统计了x方向、y方向上局部流速值与相应方向上流速平均值比值的概率密度分布函数：

$ p\left( {{l_{\text{1}}},{l_{\text{2}}}} \right){\text{ = }}\frac{{m\left( {{l_{\text{1}}},{l_{\text{2}}}} \right)}}{M} $

(4)

式中，l₁、l₂为局部流速值与相应方向上流速平均值的比值，其中l₁＜l₂；p（l₁，l₂）、m（l₁，l₂）分别为局部流速值与相应方向流速平均值的比值处于[l₁，l₂]区间的概率及格点数，M为孔隙空间中流速统计格点总数。

2 结果与讨论 2.1 流速空间非均质性差异

本文首先对地下水平均流速相同时，不同胶结度的饱和土壤介质中地下水渗流特性进行分析研究。平均流速为1 m·d^–1时的地下水渗流模拟结果如图 1所示。从图中可看出，当胶结度不同时，流速空间分布特征差异显著。如在土壤颗粒完全不胶结的介质PM40中，不同位置处流速的差异较弱，流速空间分布较为均一，无明显的流速较小、较大的区域。而随着胶结度的增大，流速空间分布的非均质性逐渐增强，局部流速远小于平均流速的区域、远大于平均流速的区域均明显增多。如在胶结度最大的介质PM25中，流速远大于平均流速的区域显著增多，且相互连接形成了连续的优势渗流通道。

变异函数CV_u、CV_ux、CV_uy计算结果如表 3所示，对于PM40、PM35、PM30、PM25四种介质，流速变异函数CV_u分别为1.233、1.478、1.921、2.098。从PM40到PM25，随着胶结固体颗粒的比例由0增至60.20%，CV_u增大了70.15%。这表明流速空间非均质性显著增强。此外，CV_ux、CV_uy结果显示从PM40到PM25，x、y方向上流速分量的变异函数分别增大了1.98倍（从1.172到8.237）和2.24倍（从1.985到4.446）。这表明随着胶结度增大，沿地下水主要流动方向、垂直于地下水主要流动方向上的流速分量空间非均质性均显著增强。

2.2 流速概率密度分布差异

x方向、y方向上局部流速值与相应方向上流速平均值比值的概率密度分布函数统计结果如图 2所示。分别绘制了常数坐标图和半对数坐标图，以更好地对比流速较大区域、流速较小区域的概率密度差异。

与上述分析结果一致，土壤颗粒胶结程度对两个方向上的流速概率密度分布均有影响。在PM40中，x方向、y方向上流速分量的概率密度函数分布均较为集中，流速数值多数分布于接近平均流速的范围内，流速较大、流速较小的区域所占比例均较低。而随着胶结度增大，流速概率密度分布函数越来越发散，局部流速接近平均流速的区域所占比例大幅度降低，流速远大于平均流速的区域、远小于平均流速的区域均迅速增加。例如在PM40、PM35、PM30、PM25中，|u_x|/u_x，ave值为1的峰值处对应的比例分别为6.59%、5.18%、3.61%、2.44%（表 4，Pc_peak）。即随着胶结土壤颗粒比例从0增至60.20%，局部流速值等于平均流速的区域的比例降低了62.97%。此外，本文发现，沿x方向和y方向上流速较小区域的流速值均集中分布于平均流速值的10^–13倍至10^–10倍这个较小的区间内。在PM25中，|u_x|/u_x，ave在[10^–13，10^–10]区间处形成的流速分布峰的峰值处的比例为1.74%，占|u_x|/u_x，ave值为1的峰值处所对应比例（3.29%）的52.89%。这表明局部流速值处于[10^–13，10^–10]间的区域在地下水渗流区中占有的比例已经非常大。

同时本文统计了平均流速分别为0.1 m·d^–1、10 m·d^–1时，不同胶结度的介质中沿x、y方向上的流速分量概率密度分布函数，并与平均流速为1 m·d^–1时的速度概率密度函数对比，结果如图 3所示。从图中可看出，在胶结度由小到大的四种介质中，当平均流速变化时，速度的概率密度分布函数基本不变。这也说明土壤介质孔隙结构对地下水流速概率密度分布具有决定性的影响。此外，虽然已有研究表明通过模拟固体颗粒膨胀作用构建的土壤介质能够反映胶结度不同时介质微观结构特征及其差异^[5]，但是，本文未考虑胶结物透水性及内部孔隙等自身特性对地下水渗流的影响。因此关于胶结物自身特性对地下水渗流的影响尚需进一步深入研究。

2.3 不流动区及优势流区域差异

如引言中所述，地下水中不流动区和优势流区的空间分布、所占比例为地下水渗流特性的重要组成部分，是土壤等多孔介质中地下水渗流特性的重要研究内容。根据已有研究，将流速值小于平均流速100倍的区域定义为不流动区，将流速值大于平均流速5倍的区域定义为优势流区^[6]。四种不同胶结度的介质中不流动区、优势流区的空间分布如图 4所示。从图中可以看出，PM40中几乎不存在不流动区和优势流区。在PM35、PM30等土壤颗粒轻微胶结的介质中，开始出现零星分布的不流动区和优势流区。这些区域分散分布于整个孔隙空间中，相互之间连通性较差。此外，随着胶结度增大，不流动区、优势流区同时扩大。在胶结度最大的土壤介质PM25中，不流动区不再是分散分布的零星区域，而是相互连接形成面积较大的不流动区。同时，优势流区相互连通，逐渐形成较为连续的优势渗流通道。而已有研究表明优势渗流通道的长度对地下水中溶质运移有显著影响^[19]。

为定量衡量不同胶结度的土壤介质中不流动区和优势流区的差异。本文统计了平均流速为1 m·d^–1时，四种介质中不流动区的比例Pc_stag、优势流区的比例Pc_pref，如表 3所示。结果表明，从PM40到PM25，地下水中不流动区比例Pc_stag从2.06%增至48.31%，增大了23倍多。而优势流区的比例Pc_pref从0.27%增至2.41%，增大了近9倍。表明随着胶结度增大，地下水中的不流动区、优势流区同时显著增多。

2.4 渗流特性差异内在机理

从速度概率密度分布函数可看出流速大小为平均流速的10^–13倍到10^–10倍的区域是不流动区的主要部分。本文深入分析了该类区域所处的空间部位，进而探讨不同胶结度的土壤中地下水渗流特性变化的内在机制。图 5为PM25中相应区域的空间分布。

从图 5中可看出不流动区主要分布在三类部位。首先，土壤颗粒胶结通常导致部分原本相互连通的孔隙形成典型的死端孔隙结构。流体进入死端孔隙后流动缓慢形成不流动区域，如图 5a所示。此外，土壤颗粒相互胶结形成的胶结团，并不简单地等效于较大的土壤颗粒，胶结团表面通常形成不规则的流固界面，这在一定程度上增加了固体介质的表面粗糙度。在形状不规则的固体表面分布有大量的不流动区，如图 5b所示。最后，土壤颗粒胶结增大到一定程度，常导致流体不能在靠近上游的一侧直接流入部分区域，而是通过下游与流动区域连通的孔隙反向流入该区域，从而形成面积较大的不流动区，如图 5c所示。地下水流经上述形成不流动区域部位的阻力增大，从而汇入连通性相对较好的孔隙部分，引起优势流区域扩大。在不流动区、优势流区同时增大的综合影响下，最终导致流速空间非均质性增强。

3 结论

本文基于孔隙尺度模拟，研究了土壤颗粒不同程度胶结的饱和土壤介质中地下水渗流特性的差异，并探讨了渗流特性变化的内在机制。得到如下结论：不同胶结度的土壤中地下水渗流特性有显著差异。随着胶结度增大，流速空间非均质性显著增强，流速概率密度分布函数越来越发散，不流动区、优势流区比例同时显著增大。从PM40到PM25，随着胶结的土壤颗粒百分数Pc从0增大到60.20%，地下水渗流速度空间变异函数增大了70.15%，不流动区比例增大了23倍多，而优势流区的比例Pc_pref增大了近9倍。此外，研究表明地下水渗流特性的以上变化主要是由胶结度增大时土壤介质孔隙尺度结构差异引起的。在胶结形成的死端孔隙结构、胶结团不规则边界处、上游连通性差下游与流动区相连通的区域通常形成不流动区。不流动区增加，流体汇入连通性较好的孔隙中快速流动，导致优势流区也随胶结度增大而增多，甚至形成连续的优势渗流通道。在上述二者的共同作用下，流速空间非均质增强。本文的研究可丰富土壤等多孔介质孔隙尺度结构差异对地下水渗流特性影响领域的研究，并为研究多孔介质中溶质运移行为及其内在机理提供研究基础。