含单裂隙非饱和带中轻非水相流体修复的数值模拟

2020-09-27赵科锋王锦国曹慧群

水文地质工程地质 2020年5期

赵科锋，王锦国，曹慧群

(1.长江科学院流域水环境研究所，湖北武汉 430010；2.河海大学地球科学与工程学院，江苏南京 211100)

随着石油工业的不断发展，输油管线、储油罐、油槽车和油轮的泄漏事故不断增加，受石油污染的土壤和水体显著增加。为此，针对土壤和地下水中有机污染开展调查、评价，进行污染控制和治理已势在必行[1-3]。目前，土壤及地下水污染物理修复方法已有多种，包括：抽水处理法、污染土壤气提法、地下水曝气法和自然衰减法等。其中，多相抽提(Multi-Phase Extraction，MPE) 技术属于抽水处理法的一种土壤原位修复技术，主要通过使用真空抽提手段，同时抽取污染区域中的气体、地下水及浮油层，达到降低污染物浓度的目的[4-5]。相比传统的修复技术，多相抽提影响半径更大，可减少布井数量，能同时清除非饱和带及饱水带中的轻非水相流体(Light Non-Aqueous Phase Liquids, LNAPL)，尤其适合挥发性较强的污染物，同样适合渗透性较低的介质[6-8]。但MPE系统抽提设备复杂，启动及调试时间较长，修复成本高，限制了其使用及推广[9-10]。

国内外学者在油污染治理方面做了一些研究。王晓燕等[11]综合分析了SEAR技术修复土壤及地下水中NAPL污染的原理及发展。胡黎明等[1]通过离心模型试验研究了石油烃在非饱和土体中的迁移过程及SVE修复过程中浓度变化特征。王磊等[12]介绍了MPE技术的工艺特点和技术原理，总结了非饱和带、含水层性质和污染物特性的MPE适用性评估方法。张晶等[13]应用多相抽提技术回收土壤及地下水中的LNAPL，并评估了修复效果。白静[14]通过试验，利用表面活性剂强化GCW修复效果，为地下水中挥发性和半挥发性有机物的去除提供新途径。Leif Nelson等[15]通过试验，评估了在抽提井两侧注入水-CO2混合物，治理LNAPL污染的可行性。Jaehak Jeong等[16]建立数学模型预测了多层含水介质受到污染时LNAPL迁移与修复特征。Fatemeh Ebrahimi等[17]建立了一个分析模型，用于估算非均质含水层中LNAPL的分布和回收率。国外在油污染修复治理方面已有丰富的经验，相关研究成果已用于治理场地土壤及地下水污染。国内学者研究主要集中在污染物迁移规律及物理化学修复，相关学者虽已开展污染物控制与治理方面的研究[18]，但在原位冲洗修复技术方面的研究较少，鲜见原位冲洗技术相关机理方面的研究。在原位冲洗修复过程中，对污染物迁移规律认识不足，且限于技术与设备的因素，国内在应用相关技术修复土壤及地下水污染的工作较少。本文研究了土壤及地下水中LNAPL在单井抽提修复及原位冲洗修复时的运移规律，检验了LNAPL修复效果。

在评估修复可行性时，数值模拟相比物理试验具有技术成本低廉、可模拟复杂边界条件、模拟尺度大等优势。本文通过建立数值模型，分析了LNAPL在含裂隙基质中迁移与分布的规律。以此为基础，建立MPE抽提修复模型，分析抽提井的流量、位置及原位冲洗方法中注水井的位置等对修复效果的影响及LNAPL的时空分布规律。

1 模型网格剖分方案

TOUGH2(transport of unsaturated groundwater and heat)是由美国劳伦斯-伯克利国家实验室开发，用于模拟多维孔隙或裂隙介质中，多相、多组分混合水流及热运移的通用数值模拟程序[19]。T2VOC是TOUGH2中侧重于研究多孔介质或裂隙介质中水相、气相、非水相(NAPLs)三相流问题的子模块。主要用来模拟非水相流体(NAPLs)在饱水带及非饱水带中时空变化规律，包括非饱水带中有机污染物挥发、扩散及真空抽取模拟，向污染介质中注入空气及液体修复等模拟。

本文采用TOUGH2软件的T2VOC模块建立二维数值模型，模型网格按X、Z方向剖分，模拟区域有效网格长30 cm，高40 cm，模型网格剖分见图1。网格为不均匀剖分，X方向剖分规则为：-0.5～0和30.0～30.5为左右隔水边界；0～14和16～30分别平分为7个长度为2 cm的网格；14.00～14.85和15.15～16.00剖分为长度0.85 cm的网格，14.85～15.15剖分为中间长为0.3 cm的网格，有效网格总长度为 30 cm，网格不规则剖分的目的，主要是为了刻画中间0.3 cm宽的裂隙。Z方向剖分规则为：-0.5～0为底部隔水边界，边界中间处有0.3 cm宽的排泄孔；0～3为饱水带，水的模拟饱和度为1；3～14为毛细带，总高度为11 cm，平分为5组单元格，每格高2.2 cm；14～37为土壤带，此带与大气相通，14～15为1.0 cm厚的网格，15～37平分为11组，每组高为2.0 cm；37～40为LNAPL注入区，平均分为三组，每组高1.0 cm。根据以上规则，共剖分357个有效单元格，61个无效单元格。

图1 网格剖分图Fig.1 Schematic diagram of the grid design

2 定解条件及参数

建立模型时假定：

(1)模拟的基质中存在水相、气相、LNAPL三相，油、水不混溶，模型中的渗流为等温的三相渗流。

(2)基质和流体不可压缩，水和油的密度不变。

(3)不考虑水的蒸发，不考虑油的挥发。

(4)LNAPL渗入前，模型中不存在油相。

(5)考虑毛细力、重力及吸附作用，忽略降解作用。

2.1 边界条件

模型的上边界为注入区，无基质，与大气相通，压强为大气压，无水相。油相均匀从上边界注入，注油速率为1.44×10-4kg/s。下边界为饱水带，水的饱和度为1，底部有排泄孔；侧面边界，为隔水边界，水相与油相均不能透过。

2.2 初始条件

模拟初始时刻模型中没有油相，饱水带水面水平，水的饱和度为1，高度为3 cm；毛细带从上至下水的饱和度逐渐增加，高度为14 cm；土壤带只有气相，压强为大气压强。初始条件见表1。

表1 初始条件

3 模拟结果及分析

3.1 LNAPL注入后运移分布

LNAPL均匀注入裂隙-基质系统后，平铺在系统表面，在重力作用下向下迁移[21-22]，湿润性液体LNAPL进入干燥的土壤带后，优先进入裂隙中。LNAPL初始干-湿界面呈狭长形，污染区域饱和度相对较低；随着LNAPL的不断注入，污染面积不断增大，LNAPL到达裂隙底部后，在裂隙中聚集，以裂隙为核心向四周扩散，饱和度以裂隙为中心向外逐渐减小[20]。LNAPL迁移图像见图2。

从模拟结果可以看出，LNAPL在10 min到达毛细带，由于毛细力的作用，LNAPL主要以横向扩散为主，干-湿界面由狭长形逐渐变成椭圆形(图2中20～40 min)，最后填充整个系统，形状如矩形(LNAPL运移50 min以后)。相对于裂隙，两侧基质中LNAPL入渗明显滞后，初始锋面较为均匀、平行，随着时间的推移，由于受到裂隙的影响，裂隙入渗口处的基质受到越来越多的污染，相比距裂隙较远的基质，LNAPL的入渗深度较大，从图2图像变化可以看出，距离裂隙越远，基质受到裂隙的影响越小，入渗深度也相对较小。LNAPL入渗干-湿界面不是平行地到达饱水带，而是以裂隙为中心，呈圆滑的锥形迁移到饱水带。LNAPL入渗的过程就是一个油相驱替水相、气相，并占据其空间的过程，LNAPL在裂隙-基质系统中主要以裂隙为介质向下入渗至饱水带[23-24]，裂隙的存在，起到了“快速通道”的作用，降低了LNAPL在中砂基质顶部向下入渗的速率，在整体上缩短了LNAPL入渗至饱水带的时间，加速了基质污染的速率，加重了污染程度。

3.2 LNAPL停止注入后饱和度分布

LNAPL在模拟的第50 min到达饱水带，之后停止注入，由于密度比水小，LNAPL漂浮在饱水带表面，不断聚集，LNAPL停止注入后的迁移图像见图3。从图中看出，LNAPL停止注入后，裂隙中LNAPL饱和度相比两侧基质中的LNAPL饱和度小，停止注入20 min时，裂隙中LNAPL的饱和度已经降到0.10，并且随着时间延长仍在继续减小，到停止注入120 min时，裂隙中LNAPL饱和度为0.07，饱和度幅度变化较小；基质中LNAPL饱和度很大，停止注入20 min时最大饱和度为0.40，到第120 min时，最大饱和度为0.20，饱和度随着时间变化幅度较大。这说明，孔隙越大，吸附性越差，残留的LNAPL越少。

图2 LNAPL注入时饱和度变化图像Fig.2 Migrate images during the LNAPL injection

图3 LNAPL停止注入后饱和度变化图像Fig.3 Migrate images after the LNAPL injection

LNAPL在饱水带聚集一定量后，由于重力作用，继续向下运移，饱水带水面向下凹陷，中间凹陷比两侧大，这主要是因为，LNAPL通过中间的裂隙首先达到饱水带，并大量聚集，中间聚集的LNAPL较两侧多，因此会形成图3中所示的中间凹陷较多的情形。LNAPL聚集一定量后，迫使水面下移，饱水带中的水从底部的排泄孔排出，多余的LNAPL从排泄孔排出，最终形成油-水平衡状态，饱水带中同时存在油相、水相。

LNPAL在迁移到饱水带的过程中，仍有部分LNAPL由于受到土壤颗粒的吸附作用长期残留在基质中，同时饱水带表面也会漂浮大量的LNAPL，这些微溶于水的污染物会长期存在，对土壤和地下水造成严重的污染，因此，对污染区域进行修复显得很有必要。

3.3 单井抽提修复

单井抽提修复是MPE修复的一种，抽提量是影响抽提修复的关键因素之一。抽提量过小，残留污染物难以被抽出且抽取时间长，效率低；抽提量过大，对设备要求高，成本费用高。抽提井真空泵位置同样影响抽提修复效果。因此，针对污染区域及LNAPL的性质，选择合适的抽提量及抽提井布置位置，将有助于提高污染区域的修复效果。

为了研究抽提量对受污染土壤修复的影响，设计井深30 cm，最大抽提量分别为1.0×10-4，3.0×10-4，5.0×10-4，7.0×10-4kg/s四种方案(图4a)；为了研究抽提井深对修复效果的影响，设计抽提流量为3.0×10-4kg/s，井深20，25，30，35 cm四种方案(图4b)。抽提压强均为370 Pa，模拟时长为2 h。

图4 抽提井抽提修复方案Fig.4 Extraction remediation scheme

3.3.1抽提量的影响

不同抽提量下LNAPL饱和度随时间变化见图5。从模拟结果可以看出，在无抽提井时，毛细水带基质中的LNAPL饱和度随着时间逐渐减小，饱和度最大值由0 min的0.20 减小到120 min的0.15，同时受污染的面积及污染程度也在减小。这表明，LNAPL在重力作用下仍向饱水带迁移，但迁移速率越来越小。从图5(a)可以看出，在第100～120 min，LNAPL饱和度几乎没有变化，油相、水相逐渐达到平衡。在有抽提井且抽提流量不为0的模型中，LNAPL饱和度随时间逐渐减小，不同抽提流量下，LNAPL饱和度变化程度不同。从图5(b)～(e)图像变化可以看出，同一时刻，模型抽提流量越大，LNAPL饱和度变化越大，说明从基质中抽取出的LNAPL越多，修复效果越好。从修复效果看，抽提量为1.0×10-4kg/s时，在第120 min，饱和度最大值为0.13，下降了0.07，仍有较大面积的污染，且污染区域的LNAPL饱和度大于0.10，修复效果最差。抽提量为3.0×10-4kg/s时，在第60 min，饱和度最大值为0.12，在120 min时，污染区域LNAPL饱和度在0.10以下，裂隙中饱和度最小，修复效果较好。修复效果最好的是抽提量为7.0×10-4kg/s的模型，在第20 min饱和度达到约0.12，在第120 min，饱和度最大值为0.07，下降了0.13，污染程度减小的最多，大部分污染区域经抽取LNAPL后，饱和度下降到约0.05。

由此可以看出，抽提流量越大，修复效果越好，但是水相、气相、油相被抽取出的越多，造成地下水面下降，可能形成新的环境问题；抽提量越大，成本费用越高。因此，在修复LNAPL造成的污染时，要充分考虑修复效果与成本两方面的因素，合理设置抽提流量及时间，在达到修复污染场地的同时，节约人力、物力、财力。

图5 不同抽提量下LNAPL饱和度随时间变化图像Fig.5 Images of the LNAPL saturation change under different extraction volumes

3.3.2抽提井位置的影响

图6显示了抽提井真空泵在不同深度下，LNAPL饱和度随时间的变化。从模拟结果可以看出，方案a与方案b修复效果差别不大，且效果较差，到模拟结束时，仍有较大面积的污染区域LNAPL饱和度在0.16以上。相比而言，方案c与方案d修复效果较好，方案d修复效果最好，在模型运行到第80 min时，饱和度降到了0.06以下，而方案c在模型运行结束时，仍有部分区域LNAPL饱和度约为0.11，修复效果稍差。对比图像发现，在毛细带及土壤带，方案c修复效果比方案d好，模型运行结束时，方案c毛细带及土壤带中LNAPL饱和度及污染区域明显比方案d的小。

由此看来，抽提井真空泵深度对修复效果影响较大。LNAPL在重力作用下主要聚集在饱水带表面及毛细带下部，这部分是LNAPL饱和度最大的区域，方案a与方案b抽提井真空泵埋深太浅，距LNAPL聚集的区域较远，由于抽提流量及抽提压强的限制，很难清除深部的污染物；方案d抽提井真空泵位于饱水带上部2 cm处，也是大部分LNAPL聚集的区域，因而抽提效果最好，由此可以看出合理的布设井位能弥补设备不足的缺点。方案d相比方案c，毛细带上部的修复效果稍差。抽提井真空泵埋深越大，对设备的要求越高。因此，在实际操作时，需充分了解LNAPL污染物的分布范围及污染程度，合理设计抽提井真空泵埋深，在设备允许的范围内，同时达到规定的修复效果及节约成本的目的。

图6 不同井深下LNAPL饱和度随时间变化图像Fig.6 Images of the LNAPL saturation change at different well depths

3.4 注水井位置对原位冲洗修复的影响

单井修复影响半径有限，针对大的污染区域，修复效果不理想，实际操作时，多采用群井的方式进行修复[25-27]。由于水泵对地下水的抽提作用，潜水位下降，原来被水侵没的土壤暴露出来，造成部分区域水位下降，引起新的环境问题。鉴于此，在抽提井两侧添加注水井，通过向土壤污染带或饱水带中注水的方式修复污染区域，这是单井抽水处理法的衍生方法，称为原位冲洗技术。

原位冲洗技术既弥补了单井抽提影响半径过小的问题，又补充地下水，防止水位下降过多的问题，是较为理想的清除LNAPL污染的方法。注水井的布设尤其重要，为此建立相关模型研究注水井布设位置对原位冲洗法修复效率的影响。设计方案如图7所示，抽提井真空泵长度均为2 cm，位于横向中间位置，竖直方向位于10～12 cm，抽提流量为3.0×10-4kg/s，抽提压强为370 Pa。两口注水井流量均为1.0×10-4kg/s，其中横向(X方向)注水井布设见图7(a)。竖直方向滤网长度为10 cm，位于10～20 cm处，抽提井两侧对称布设注水井。方案a～d注水井距抽提井距离分别确定为11，8，5，2 cm。竖直向(Z方向)布设见图7(b)。两口注水井在抽提井两侧，对称分布，距抽提井距离均为8 cm，竖直方向设计4种方案，方案A注水井滤网位于20～30 cm，方案B注水井滤网位于15～25 cm，方案C注水井滤网位于10～20 cm，方案D注水井滤网位于5～10 cm。通过不同布设方案探究注水井位置对修复效果的影响。

图7 注水井布设示意图Fig.7 Schematic diagram of the water injection well design

3.4.1注水井横向位置的影响

为了研究注水井与抽提井之间的距离对修复效果的影响，设计了4种横向(X方向)注水井方案，研究LNAPL饱和度随时间变化(图8)。从模拟结果可以看出， LNAPL饱和度随时间减小，模型运行结束时，注水井附近LNAPL饱和度均保持在约0.05，比抽提井附近LNAPL饱和度小，且以注水井为中心，向外LNAPL饱和度有变大的趋势。方案a在模拟的第40 min，注水井附近的LNAPL饱和度明显降低，低饱和度的分布区有向上部及下部扩散的趋势，这是由于注水井距隔水边界较近，注入大量的水受到阻挡后，向竖直方向扩散的缘故。方案a中注水井距抽提井较远，大于两者共同作用下的影响半径；方案b位于抽提井影响半径边缘。方案b修复效果明显优于方案a和方案c，图8(c)LNAPL饱和度变化图像显示，注水井与隔水边界之间的污染物浓度变化较小，模型结束时，LNAPL饱和度最小的区域集中在注水井与抽提井之间，这是因为注水井距抽提井太近，注水井注入的水很快运移到抽提井中，被抽取出，由于运移距离过短而未能携带大量的污染物，导致只有注水井与抽提井附近污染物被抽取出，边界附近的污染物浓度依然很大。方案d修复效果最差，LNAPL饱和度只在注水井与抽提井附近的小范围内比较小，相比图5(c)同等状态下(不存在注水井的情况)，修复效果差。

由此看出，注水井与抽提井的距离不宜过大或过小。距离过大或超出影响半径，注入的水运移速率较慢，修复时间长，效果差；距离过小，注入的水很快进入抽提井被抽出，携带的污染物量小，修复效率低。实际布井时，需根据场地土壤的渗透性、颗粒大小及污染物的分布范围及污染物的性质进行合理布井[28-30]。

图8 注水井距抽提井不同距离时LNAPL饱和度随时间变化图像Fig.8 Images of the LNAPL saturation at different distances between the injection well and extraction well

3.4.2注水井埋深的影响

从图8看出，方案b是最优的修复方案，但注水井上部与下部LNAPL饱和度依然很大，污染物依旧大量存在。可见，注水井滤网埋深的布设存在不合理之处。为此，设计了图7(b)所示的布井方案，探究注水井埋深对修复效果的影响。LNAPL饱和度随时间变化图像见图9。

抽提井的位置直接影响了污染物的抽取效率，注水井的位置影响到污染物除去的范围以及修复程度。从模拟结果看，方案D至方案A修复效果依次变好。方案D的修复效果最差，LNAPL去除程度最低，残留的污染物主要集中在远离注水井的非饱和带上部，饱和度约为0.10。对比方案A与方案D，方案A饱水带以上区域LNAPL饱和度已降到0.07以下，在第80 min 时LNAPL饱和度变化已较缓慢，基本达到稳定状态，总体修复效果最好；方案D只有饱水带附近污染物浓度较低，土壤带及毛细带污染物浓度依然较高。

注水井在向受污染区域注水的过程中，除了向两侧扩散外，还在重力作用下向下扩散，注入的水不仅降低了污染物浓度，而且起到了“冲洗”吸附在土壤颗粒上污染物的作用。注水井周围的污染物随着水的注入不断被稀释，并被携带到非饱和带下部，被抽取至地表。

从图9看出，方案B、方案C、方案D修复效果不明显，主要是因为，注水井滤网距离污染范围的上边界较远，而水受重力作用只会向下运移，故而起不到稀释与“冲洗”的作用，导致修复效果不理想。在采用注水井与抽提井同时工作时，需充分掌握污染物的分布范围，注水井埋深不宜过大，滤网上部与修复范围的上边界齐平为最佳。

图9 注水井筛管距底部边界不同距离时LNAPL饱和度随时间变化图像Fig.9 Images of the LNAPL saturation at different distances from the bottom border of the injection well screen

3.5 原位冲洗修复下LNAPL饱和度变化规律

原位冲洗技术设计主要有布井位置、数量、滤网位置、运行参数等。井位布设建立在充分了解污染物分布的基础上，在设计时要充分考虑设备的可行性、成本费用及环境保护等问题[31-32]。

3.5.1布井方案

设计注水井、抽提井及观测点位置，如图10 所示。考虑到效率因素，抽提井流量取3.0×10-4kg/s，抽提真空泵位于5～7 cm处(图6中最优方案)，注水井综合考虑横向与竖直方向的影响，选取最优方案(横向距抽提井8 cm，竖直方向滤网位于20～30 cm处)。观测点监测相应点位LNAPL饱和度的变化，因LNAPL饱和度分布以中间裂隙为轴左右对称，故点位布设裂隙左侧。观测点Ob1、Ob2、Ob3探究不同深度LNAPL饱和度的变化，观测点Ob4、Ob5、Ob6分别监测邻近裂隙、裂隙中及抽提真空泵处的LNAPL饱和度变化。从模型运行开始到结束，6 个观测点共采集数据180 组。

3.5.2不同位置LNAPL饱和度变化规律

图11为各观测点LNAPL饱和度变化情况。观测点Ob1、Ob2、Ob3、Ob6处LNAPL饱和度均是随时间逐渐减小，四个观测点饱和度降幅分别为0.108，0.196，0.122，0.128。同一时刻，各观测点LNAPL饱和度：Ob1

从观测井Ob1、Ob2、Ob3处饱和度对比看出，深度越大，LNAPL饱和度越大；从饱和度变化速率看，Ob2处饱和度变化最大(Ob1、Ob2、Ob3处前4 000 s LNAPL饱和度下降平均速率分别为2.61×10-5，2.90×10-5，2.73×10-5s-1)。对比Ob4、Ob5、Ob6处饱和度，前1 300 s，Ob5Ob6。

观测点Ob4位于注水井中下部分靠近裂隙的地方，当注水井开始注水，会迫使水-油锋面向Ob4迁移，增大LNAPL的浓度，当水-油锋面过后，注入的水开始运移到Ob4处，此时LNAPL浓度逐渐降低。图11中观测点Ob4处LNAPL饱和度在615 s之前增大，在之后缓慢降低。随后水-油锋面到达裂隙中的观测点Ob5处，导致LNAPL饱和度缓慢增大，这是由于裂隙渗透性很大，模型未运行之前LNAPL饱和度已近乎最低，处于平衡状态，随着注水井的注水，水-油混合物迁移到裂隙中，增加了裂隙中污染物浓度。当注水井继续工作时，大量的水运移到裂隙中，稀释了裂隙中的混合物，同时裂隙中的水-油混合物在重力作用下向下运移，进入到抽提井中抽出至地表，于是在模型运行4 800 s 后，Ob5处再次出现LNAPL饱和度下降。初始时刻观测点Ob6处LNAPL饱和度最大(饱和度为0.2)，相比观测点Ob3初始值0.197略大，但随模型运行饱和度很快变小，并在同时刻小于Ob3处的值。

观测点Ob6位于抽提井真空泵处，随着抽提作用，大量的水相、气相、油相混合物被迅速抽出系统，因而LNAPL饱和度会在短时间内下降；观测点Ob3距抽提井相对较远，受抽提作用影响较小，在359 s之前饱和度降速较慢，在之后抽提作用影响扩展到观测点Ob3，故而在359 s之后，LNAPL饱和度迅速下降。对比Ob4与Ob6，Ob6处LNAPL下降平均速率为2.10×10-5s-1，远离注水井与抽提井的Ob4处LNAPL饱和度变化相对缓慢很多，速率为5.79×10-6s-1。由于远离注水井，LNAPL得不到有效“冲洗”，无法在短时间内被携带至抽提井处抽出，同时由于没有有效的稀释，故而变化速率较为缓慢。

3.6 原位冲洗修复效果检验

根据原位冲洗修复时LNAPL饱和度变化规律及井位布设规律，原位冲洗修复最优模式为：抽提井埋深到达饱水带(35 cm处)，注水井横向位于抽提井影响半径边缘(距抽提井8 cm)，竖向位于污染区域上边缘。

图12为原位冲洗修复最优模拟方案下，LNAPL饱和度随时间变化图像。从模拟结果可以看出，在第80 min时，污染物已得到有效控制，修复面积达到96%，修复率达到75%。在第80～120 min污染物浓度变化不大，但仍在继续减小。到第120 min时，修复面积接近100%，LNAPL饱和度控制在0.05。由此看出，原位冲洗方法能较好地清除LNAPL。相比单井抽提模式，原位冲洗模式修复效率高、效果好，在某种程度上能弥补对设备要求较高的限制，而且不会引起大面积的地下水位下降，避免在修复污染物时造成新的环境问题。