裂隙宽度空间变异性和泄漏条件对网络裂隙中DNAPLs运移影响研究

常 兴，骆乾坤，邓亚平，马 雷，钱家忠

（合肥工业大学资源与环境工程学院，安徽 合肥 230009）

地下水环境中存在大量难以消除的污染物，其中包括重非水相液体(dense nonaqueous-phase liquids，DNAPLs)。因为其与水不相容，在运移过程中由于毛管压力、黏滞力和重力的共同作用下会残存在运移路径中，长期存在且难以消除。裂隙介质是地下水三大赋存介质之一[1]。由于裂隙介质的各向异性和非均质性明显强于孔隙介质，所以裂隙中地下水受到DNAPLs污染后更难修复。为此，研究裂隙介质DNAPLs 的运移规律对此类污染物的修复尤为重要。

目前DNAPLs运移研究大多集中在孔隙介质[2-3]，主要是通过数值模拟[4-7]和室内试验[8-9]来研究孔隙介质的非均质性[10-12]以及地下水流流速[13]等对DNAPLs运移的影响。裂隙介质中的DNAPLs运移研究还处于起步阶段[14-15]，主要研究裂隙的相关参数以及裂隙的物理性质对DNAPLs运移及分布的影响。就目前来说，研究裂隙宽度空间的变异性主要围绕单裂隙介质[16-17]，裂隙隙宽的改变影响裂隙中水流的流态、流速。从裂隙水流态的研究出发延伸至单裂隙中裂隙隙宽的变异性对DNAPLs运移的影响，主要通过室内试验[18]、数值模拟[19-23]及求解对应数学模型[24]开展研究。裂隙隙宽变化，会改变DNAPLs 的流量、流速，最终影响DNAPLs 在裂隙介质中的运移范围和空间分布。泄漏速率对DNAPLs运移的研究集中在孔隙介质中[25]，研究结果表明DNAPLs 泄漏速率不同，泄漏总量不同，导致污染羽的空间展布不同。关于泄漏位置对DNAPLs 在裂隙中运移影响的研究很少[26]。

综上所述，国内外对DNAPLs运移规律的研究大多在单裂隙或多孔介质中进行，而对随机网络裂隙中DNAPLs运移规律及影响因素的研究不足，综合考虑裂隙宽度空间变异性和泄漏条件对DNAPLs 在裂隙中的运移影响的研究更少。因此本文结合真实网络裂隙数据，通过蒙特卡罗方法建立裂隙宽度空间变异性不同的随机网络裂隙，采用PetraSim软件模拟并分析裂隙宽度空间变异性、泄漏速率和泄漏位置对DNAPLs运移的影响，丰富裂隙介质中DNAPLs运移机理研究，为裂隙介质中DNAPLs 污染修复提供模型参考。

1 随机网络裂隙场的生成

为模拟天然裂隙网络中裂隙宽度空间变异性对DNAPLs运移的影响，本文基于已有统计资料[27]，使用蒙特卡罗方法(蒙特卡罗方法可以根据裂隙各几何参数所服从的分布函数，编程生成随机网络裂隙。此方法能最大程度模拟真实裂隙，是目前广泛使用的生成随机网络裂隙的方法)生成裂隙宽度空间变异性（方差不同，表明生成的裂隙场和裂隙宽度不同，此不同即为空间变异性）不同的3 组网络裂隙物理模型，主要参数见表1。

表1 不同裂隙场中的裂隙几何参数Table 1 Parameter setted for different fracture

将生成的随机网络裂隙坐标化。步骤如下：

第一步：按表1 中的裂隙几何参数生成随机线段，设定随机线段结果输出区域大小为100 cm×100 cm 的正方形区域，对生成的线段进行统计分析；

第二步：进行像素扫描，将上述区域离散为水平和垂直间距均为1 cm 的小网格共10 000 个，使用Python语言编程识别生成的线段所覆盖的网格块，并对其进行标记；

第三步：输出所有标记线段首尾网格的坐标及宽度信息，生成随机网络裂隙场。基于不同裂隙宽度空间变异性（σ2）的3 组随机网络裂隙场如图1所示。

图1 基于不同裂隙宽度空间变异性下的随机网络裂隙场Fig.1 Stochastic network fractures with spatial variability of different widths

2 算例介绍

基于生成的随机网络裂隙场，建立DNAPLs运移数值模拟模型。假设模拟区域为100 cm×100 cm 的二维XZ剖面，X轴方向和Z轴方向剖分步长均为1 cm，共离散10 000 个网格。地下水面位于模型顶部。左右边界均为定水头边界，上下边界为无通量边界。模型首先达到水-气平衡，初始水相饱和度Sw为1.0。模型区域内水力梯度J设置为0.01，水流方向为从左向右（图2）。因为四氯乙烯（PCE）是目前地下水中普遍存在的污染物，且已经成为国内外专家研究地下水生物修复的热点问题，所以选取PCE 作为DNAPLs 示例。由于PCE 很难溶解于水相中(溶解度约为150 mg/L)，故本文不考虑PCE 的溶解及吸附，假设其完全以非水相形式存在。由于本次模拟温度设定为20℃，所以不考虑热物理相关参数。PCE 以线源方式（共3 个网格）从模型顶部3 个可能位置（图2）以恒定的速率泄漏。综合考虑模型尺寸和泄漏速率，设定总模拟时间为600 s。模型主要参数设置见表2。

图2 模拟区域概念模型图Fig.2 Conceptual model of simulated area

表2 模型参数值Table 2 Model calculation parameters

3 试验方法

本次采用多相流数值模拟软件PetraSim 进行模拟。将生成的随机网络裂隙读入PetraSim 模拟软件，进行不同裂隙宽度空间变异性、泄漏速率和泄漏位置下PCE 的运移模拟。通过空间矩评估PCE 污染羽的空间变化特征，以此分析PCE 的运移规律，其中分别用零阶矩、一阶矩和二阶矩来表征PCE 的总质量、质心位置和空间展布。由于本例不考虑污染物运移过程中的化学反应，所以污染物的总质量（零阶矩）在运移过程中保持不变。

3.1 裂隙宽度空间变异性

选择3 种不同裂隙宽度方差（σ2=0.25，0.50，0.75）的网络裂隙场进行模拟，泄漏位置为图2 中位置2，并对模拟结果进行对比分析。

3.2 泄露速率

在同一网络裂隙中，选取3 种泄漏速率在裂隙场中进行模拟，泄漏位置为图2 中泄漏位置2。高、中、低泄露速率分别为1.4×10-4，1.4×10-5，1.4×10-6kg/s[25]。3 种裂隙场中的模拟结果规律相似，受文章篇幅所限，仅选择 σ2=0.25 的裂隙场作为代表，并对模拟结果进行对比分析。

3.3 泄漏位置

在同一网络裂隙场中，设置3 个泄漏位置（图2），进行不同泄漏位置PCE 运移模拟。3 种裂隙场中的模拟结果规律相似，受文章篇幅所限，仅选择 σ2=0.75的裂隙场作为代表，并对模拟结果进行对比分析。

4 结果

图3 显示了不同裂隙宽度空间变异性的网络裂隙中PCE 在模拟末时刻的饱和度分布，裂隙宽度空间变异性的改变使得PCE 空间展布发生改变，PCE 到达模型底部的时间不同，PCE 饱和度不同。

图3 网络裂隙不同空间变异性下末时刻PCE 的饱和度分布Fig.3 Corresponding saturation distribution of DNAPLs for different spatial variability

图4 为3 种不同泄漏速率下在模拟末时刻的PCE 饱和度分布图。由于泄漏速率不同，相同时间内释放的PCE 总量不同，空间展布不同，最终PCE 饱和度不同。

图4 不同泄漏速率条件下末时刻PCE 的饱和度分布Fig.4 Corresponding saturation distribution of DNAPLs for different leakage rate

图5 为同一模型不同泄漏位置末时刻的PCE 饱和度分布图。在同一模型中，泄漏位置的不同导致相同时间内PCE 污染羽分布范围不同，饱和度不同。

图5 不同泄漏位置条件下末时刻PCE 的饱和度分布Fig.5 Corresponding saturation distribution of DNAPLs for different leakage location

5 分析与讨论

5.1 裂隙宽度空间变异性对PCE 运移范围的影响

裂隙宽度空间变异性不同，PCE 运移范围不同。分析原因，由于宽度较大裂隙内毛管压力较小，构成优势通道，控制PCE 运移路径，PCE 能够很快到达模型底部边界，所以PCE 多集中在宽度较大的裂隙中。同时，由于裂隙宽度在垂向和横向上的分布不同，致使PCE 饱和度空间展布不同。图6（a）显示了具有不同裂隙空间变异性下的PCE 在Z轴的质心位置。空间变异性大的网络裂隙（σ2=0.75），PCE饱和之后通过宽度较大的裂隙通道继续运移至模型底部其他横向分布的裂隙中，因而PCE 质心位置随时间降低减缓。裂隙宽度空间变异性小( σ2=0.50，0.25)时，泄漏的PCE 到达模型底部后大多在运移通道内堆积。另外，由于裂隙的宽度突变和交叉角度的改变使得局部裂隙通道内蓄积的PCE 的饱和度增大，减小了PCE 的运移速率，因而PCE 质心位置在模拟后期呈上升趋势。图6（b）显示了裂隙宽度不同空间变异性的网络裂隙中污染羽围绕质心的空间分布，横轴表示不同方差的裂隙场，纵轴表示不同裂隙场对应的横向，纵向的空间分布。不同裂隙宽度方差条件下污染羽绕质心的垂向展布与横向展布均存在差异。裂隙空间变异性的改变，使得PCE运移路径和空间展布方向不同。裂隙空间变异性小，裂隙宽度离散程度较小，隙宽相对集中，PCE 通过裂隙通道向周围运移，PCE 饱和度的空间变异性较小。随着裂隙空间变异性的增大，裂隙隙宽离散程度增大，相对较窄和相对较宽的裂隙占比增大，相对较宽的裂隙通道构成优势通道，PCE 通过优势通道快速向周围运移，使得PCE 饱和度的空间变异性增大。

图6 不同空间变异性下污染羽空间分布特征Fig.6 Spatial moments of DNAPLs contaminant plume for different spatial variability

5.2 泄漏速率对PCE 运移速率和空间分布的影响

图7（a）显示了3 种不同泄漏速率情形下PCE 在Z轴的质心位置。从图中可得，高泄漏速率由于流量大，相同时间内PCE 泄露总量最多，受重力影响最大，因此能最快克服毛管压力，集中向下运移，因此PCE能最快到达模型底部并在裂隙通道内堆积，导致质心位置快速降低。相较于高泄漏速率，中泄漏速率相同时间内PCE 泄漏总量减少，导致其重力不易克服毛管压力，虽然污染羽能到达模型底部并在模型底部裂隙通道内堆积，但是质心位置较高泄漏速率降低较慢。低泄漏速率由于流量最小，PCE 在运移过程中受毛管压力影响最大，导致运移速率最慢，运移深度最小，污染羽大多集中在注入点下方，质心位置降低最慢。图7（b）显示了3 种不同泄漏速率条件下PCE 污染羽围绕质心的空间分布。由于相同时间内高泄漏速率释放的总量最多，PCE 在裂隙相对较宽的通道最先饱和，待饱和之后，PCE 会进入相对较窄的裂隙通道中，所以空间展布范围较中、低泄漏速率大。相同时间内，中泄漏速率情形下PCE 与高泄漏速率均运移至模型底部并形成高饱和度区；但由于入渗速率较低，到达模型底部的时间较长，所以相同时间内中泄漏速率情形下PCE 空间展布范围变化较大。低泄漏速率情形下由于流量小，相同时间内释放总量最小，难以克服毛管压力，因此PCE 集中在泄漏点下方且空间展布范围最小。

图7 不同泄漏速率下空间分布特征Fig.7 Spatial moments of DNAPLs contaminant plume for different leakage rate

5.3 泄漏位置对PCE 质心位置和空间分布的影响

泄露位置不同，PCE 污染羽的质心位置和空间分布不同。图8（a）为同一模型不同泄漏位置PCE 在Z轴质心位置图。从图中可得，由于较宽的裂隙通道内毛管压力较小，利于优势通道的形成，从而控制PCE 的运移路径。PCE 通过优势通道快速运移，使得污染羽在垂直方向上的分布范围明显增大，所以PCE能够很快运移至模型底部。因此，即使同一模型中，由于泄漏位置不同，污染物运移的重力方向的空间变异性也不相同，导致PCE 质心位置也不相同。图8（b）为同一模型不同泄漏位置的PCE 围绕质心的空间分布。PCE 受重力作用向下运移，同一模型不同泄漏位置其重力方向网络裂隙分布不同，污染物到达模型底部时间不同。同样由于网络裂隙空间变异性的不同，PCE 空间分布也不相同。

图8 不同泄漏位置空间分布特征Fig.8 Spatial moments of DNAPLs contaminant plume for different leakage location

6 结论

（1）随机网络裂隙宽度的空间变异性是影响DNAPLs 饱和度分布的一个重要因素。网络裂隙宽度空间变异性的增大，增大了网络裂隙中相对较宽的裂隙的比例，相对较宽的裂隙会形成优势通道，DNAPLs 通过这些优势通道快速向周边运移，导致水平方向和垂直方向上DNAPLs 的一阶矩和二阶矩发生改变。

（2）随机网络裂隙中DNAPLs运移范围受泄漏速率影响。随着泄漏速率的增大，相同时间内，污染物泄漏总量增大，运移速率加快，空间分布范围增大，一阶距和二阶矩变化增大。

（3）同一模型中，泄漏位置不同，导致DNAPLs运移重力方向的裂隙介质的空间变异性不同，DNAPLs运移路径和污染羽分布范围也不同。