基于MODPATH 模型的单井抽出-回渗循环地下水水力控制程度的关键影响因素研究

王新港，杨 昱，王 磊，徐祥健，韩 旭，夏 甫，邓 圣，肖 瀚*，姜永海

1.河北工程大学能源与环境工程学院，河北 邯郸 056038

2.中国环境科学研究院，国家环境保护地下水污染模拟与控制重点实验室，北京 100012

针对高风险大型石化场地水文地质条件复杂、土壤和地下水中NAPL 相污染物种类多且浓度高等特点[1-4]，以及传统抽出处理异位修复存在污染羽捕获效果差、修复成本高等问题[5-7]，单井抽出-回渗循环原位修复技术创新性地利用抽水完整井将重度污染区污染地下水抽出后，先分离出NAPL 相污染物，再将污染地下水注入布设于地表的填料区处理，处理后的水在以抽水井为中心的一定范围等量回渗至包气带和含水层[8].该技术通过抽出-回渗实现重度污染区地下水流场的水力控制，防止污染物扩散蔓延的同时将受污染的地下水体与未受污染的清洁地下水体分隔，同时淋洗包气带中的污染物进入含水层并促进含水层中的污染物被持续有效抽出处理.

单井抽出-回渗循环原位修复技术的有效性主要取决于地下水流场的水力控制程度(抽水井对回渗水的捕获率，即回渗水量被抽水井抽出的比例)，足够的水力控制程度利于目标充分修复[9-11].影响地下水流场水力控制程度的因素包括抽水量[12-15]、渗透系数[16-19]、抽注水井布局[20]等.张帅等[21]研究了内循环区的水动力学特性得到回收率的理论公式和变化曲线，认为回收率取决于地下水侧向径流的方向和相对强度，随抽水流量的增大而增加.Shan[22]根据叠加原理推导了任意位置的两口井在承压含水层抽水时的捕获区解析解，发现捕获区的形状及相对位置随井向和井距的变化而变化；姜烈等[23]利用遗传算法和模拟退火法研究了不同井位布设，使得水力控制最大化，从而使硝酸盐污染地下水抽出处理成本最低；徐绍辉等[24]研究了地下水中石油类污染的时空分布规律，通过设置水力截获工程，发现水力截获可以导致地下水水流发生变化，形成新的降落漏斗，从而去除地下水中的部分污染物；汪成等[25]采用数值方法，发现水力控制可以使地下水中Cr(Ⅵ)的污染范围得到限制，从而达到风险管控的目的；郭飞等[26]在去除污染源的基础上采用水力截获技术，通过考虑污染场地的水文地质条件，调整抽水井的位置、间距、深度和抽水量得到了主要的治理策略；王新港等[8]通过实验室尺度砂柱物理模型和三维地下水流数值模型相结合的方法，三维刻画了单井抽出-回渗同步循环地下水修复技术实施时地下水流场的时空变化特征，发现影响该技术水力控制程度的关键因素为水文地质条件和水动力条件，即含水层渗透性强弱和抽出-回渗量大小.综上，已有研究仅筛选出水力控制程度的主要影响因素，但并未深入分析关键因素对水力控制程度的影响，未量化水力控制程度与关键因素之间的定量关系.

本研究构建实验室尺度砂柱物理模型，在验证单井抽出-回渗循环原位修复技术可行性的同时基于物理模型构建三维地下水流线示踪模型，通过校正的数值模型结合不同水文地质条件和水动力条件的情景设置，深入研究关键因素对抽出-回渗同步循环地下水水力控制程度的影响机制，精准识别并刻画回渗水质点的迁移轨迹，量化抽水井的捕获范围及回渗水的捕获率，明确抽出-回渗同步循环地下水水力控制程度与含水层渗透性强弱和抽出-回渗量大小之间的定量关系.本研究获取的水力控制程度与关键因素之间的定量关系对于指导该技术在实际污染场地的应用具有重要作用，可以有效指导参数设计，从而提高修复效率，缩短修复时间.

1 材料与方法

1.1 物理模型

1.1.1 装置搭建和测试方法

试验装置如图1 所示，有机玻璃柱半径为50 cm，高度为100 cm，内部填充高度为70 cm 的80～120 目(100 目=0.15 mm)的石英砂，石英砂每5 cm 进行均匀填充及饱和夯实，搭建各向同性的均质含水层.抽水管半径1 cm，长度100 cm，位于有机玻璃柱中心处，为了抽出含水层所有深度的水，抽水管全部开筛并包裹160 目筛网以防止含水层石英砂进入抽水管内.

图1 单井抽出-回渗同步循环地下水水力控制技术试验装置Fig.1 Conceptual diagram of single well pumpingrecharge synchronous cyclical groundwater hydraulic control technique device

为实时监测模拟含水层不同位置、不同深度的水位，在装置内部布设48 个水位监测点(见图2)，从上到下共计8 层，采用上密下疏方式埋深(分别为1、3、5、7、14、28、48 和70 cm)，每层采用内密外疏方式设置6 个监测点位(与井轴的距离分别为2、4、8、14、22 和32 cm).通过抽水试验、补水试验和排水试验，测定出含水层渗透系数(K)为0.01 cm/s，孔隙度为0.4，给水度为0.09，持水度为0.31.抽水管顶部通过橡胶管连接蠕动泵(保定兰格BT100-1L)，蠕动泵的抽水流量范围为0.005～12 cm3/s，本研究选取1、2.5、5、10 cm3/s 四个抽出-回渗量(Q)运行砂柱物理模型.

图2 单井抽出-回渗同步循环地下水水力控制技术试验装置Fig.2 Conceptual diagram of single well pumpingrecharge synchronous cyclical groundwater hydraulic control technique device

1.1.2 回渗装置搭建

回渗装置布设于含水层顶部以抽水井井轴为中心的范围(见图1)，有机玻璃柱半径为32 cm，铺设厚度为15 cm 的圆柱体填料层(微米零价铁和石英砂混合装填，用于去除污染物)，填料层底部与含水层顶部交界处为两层160 目砂网(两层砂网之间铺设一层直径为0.4 cm 的玻璃珠)，防止填料层中的微米零价铁和石英砂随回渗水进入含水层.抽水管抽出的水通过橡胶管输送至回渗装置中，形成一定高度的液面，通过填料层缓慢均匀回渗至含水层中.

相较于王新港等[8]构建的原砂柱物理模型采用的注水孔分散回渗方式，本研究构建的回渗装置可以大幅度增加均匀回渗能力，并保证一定的回渗范围〔回渗半径(R)为32 cm〕.

1.2 数值模型

参照实验室尺度砂柱物理模型，王新港等[8]利用GMS 模拟软件等比例构建了基于MODFLOW 模块的地下水流数值模型，用于三维刻画单井抽出-回渗同步循环过程中地下水流场的时空变化情况[27-30].本研究基于该地下水流数值模型，利用GMS 模拟软件构建了地下水流线示踪MODPATH 模型，精准识别并刻画回渗水质点的迁移轨迹，量化抽水井的捕获范围及回渗水的捕获率，明确水力控制程度[28].

1.2.1 地下水流数值模型构建

概念模型结合砂柱物理模型构建(见图1)，模拟区侧壁和底部由砂柱边界阻挡，设定为隔水边界.模拟区顶部有回渗装置，顶部设定为变流量边界(接受回渗补给).抽水井概化为定流量边界.因装置顶部设置密封盖防止水分蒸发，且抽水井抽出水同步等量回渗，因此装置内水量始终保持平衡稳定.模拟含水层水文地质参数与砂柱物理模型的水文地质参数一致.

数值模型基于MODFLOW 模块，在概念模型的基础上，结合实验室尺度砂柱物理模型等比例构建，采用有限差分的方法将模拟区域进行空间离散，剖分为100 行、100 列和4 层共计40 000 个网格，其中有效网格数为31 400.

1.2.2 地下水流数值模型验证

评价模型模拟精度采用纳什效率系数(Nash-Sutcliffe efficiency coefficient,NSE)，其中NSE 的取值范围为(-∞，1]，当NSE 为正值时说明结果可信，越接近1表明模型模拟精度越高[31].通过模拟水位与实测水位的对比校验模型，发现NSE 为0.88(见图3)，表明模型效果较好，可以较好地模拟地下水流场[32-33].

图3 模拟结果与实测结果的精度分析Fig.3 Conceptual diagram of single well pumping-rech

1.2.3 地下水流线示踪数值模型构建

基于构建并验证的地下水流MODFLOW 数值模型，构建地下水流线示踪MODPATH 模型，选取回渗范围内表层地下水质点作为模拟因子，向前(后)跟踪模拟地下水质点的运行轨迹，刻画装置内部地下水流线和抽水井的捕获范围，判断回渗范围边界处的水质点是否可以被抽出井捕获并抽出，明确抽出-回渗水力控制程度.

由于地下水流MODFLOW 数值模型已通过砂柱物理模型的实测水位校准，可以较好地模拟砂柱物理模型中不同位置、不同深度地下水水位，并较好地刻画地下水流场特征，因此未校准地下水流线示踪MODPATH 模型.

2 结果与讨论

2.1 抽水井捕获范围和水力控制程度

基于MODPATH 模型的模拟结果(见图4)，不同抽出-回渗量(1、2.5、5、10 cm3/s)条件下，砂柱物理模型抽出-回渗平衡时，抽水井捕获半径分别为31.29、31.63、31.67、31.71 cm，抽水井捕获范围分别为3 076、3 143、3 151 和3 162 cm2.由于该砂柱物理模型的回渗范围为3 217 cm2(回渗半径为32 cm)，抽水井捕获范围均小于回渗范围，因此回渗水无法全部被抽水井捕获并抽出，抽出-回渗水力控制程度分别为95.63%、97.69%、97.93%和98.17%.

图4 不同流量下回渗水质点的迁移轨迹Fig.4 Migration trajectory of re-infiltration water mass at different flow rates

2.2 关键因素对水力控制程度的影响机制研究

为了深入研究水文地质条件和水动力条件对抽出-回渗同步循环地下水水力控制程度的影响，根据不同水文地质条件和水动力条件的情景设置，构建地下水流线示踪MODPATH 模型，模拟不同情景下回渗水质点的迁移轨迹，分析抽水井的捕获范围及回渗水的捕获率，量化捕获率与含水层渗透性强弱和抽出-回渗量大小之间的定量关系.

不同水文地质条件和水动力条件的情景设置如表1 所示，其中R代表回渗半径，K代表含水层渗透系数，Q代表抽出-回渗量.为避免回渗水流靠近装置侧壁引起的边壁效应的影响，同时避免回渗水流过于集中在抽水井周围，回渗范围不应过大或过小，因此设置了25、30、32 和35 cm 四种回渗半径情景(装置半径为50 cm)；设置了粉砂(K=0.009 cm/s)、细砂(K=0.02 cm/s)、中砂(K=0.04 cm/s)和粗砂(K=0.09 cm/s) 4 种具有代表性的含水层岩性情景[34]；依据砂柱物理模型设置1、2.5、5、10 cm3/s 四种抽出-回渗量情景.

表1 不同水文地质条件和水动力条件的情景设置Table 1 Scenarios with different hydrogeological and hydrodynamic conditions

64 种水文地质条件和水动力条件情景下模拟得到的抽水井的捕获半径和捕获率如图5 和表2 所示.从图5 可以看出，不同水文地质条件和水动力条件下，当回渗半径为25 cm 时，抽水井的捕获半径为24.59～24.98 cm，捕获范围为1 899.62～1 960.36 cm2，捕获率为96.77%～99.88%；当回渗半径为30 cm 时，抽水井的捕获半径为29.28～29.89 cm，捕获范围为2 693.34～2 806.74 cm2，捕获率为95.26%～99.63%；当回渗半径为32 cm 时，抽水井的捕获半径为31.22～31.71 cm，捕获范围为3 062.07～3 158.95 cm2，捕获率为95.18%～98.17%；当回渗半径为35 cm 时，抽水井的捕获半径为33.41～34.57 cm，捕获范围为3 506.73～3 754.47 cm2，捕获率为91.15%～97.55%.可以发现，当回渗范围和含水层渗透系数一定时，随着抽出-回渗量的增加，捕获范围和捕获率均增大；当回渗范围和抽出-回渗量一定时，随着含水层渗透系数的增大，捕获范围和捕获率均减小；当含水层渗透系数和抽出-回渗量一定时，随着回渗范围的增大，捕获范围和捕获率均减小.

根据64 种情景的模拟结果，发现水力控制程度与含水层渗透性强弱、抽出-回渗量大小和回渗半径之间存在定量关系，因此对数据进行拟合，其中捕获率是水力控制程度的评价标准.作为优化目标，选择线性拟合、指数拟合和对数拟合处理数据.拟合效果显示，对数拟合效果优于另外两种模拟结果，可以较好地拟合模型输出数据，此外，影响捕获率的关键参数为含水层渗透性强弱、抽出-回渗量大小和回渗半径，采用对数拟合后所得的公式相对于另外两种拟合公式更加简单且方便应用.拟合结果如图6 所示，拟合获取的定量关系表达式如式(1)所示.

图6 捕获率与含水层渗透性强弱和抽出-回渗量大小之间的定量关系的数据拟合结果Fig.6 Data fitting of quantitative relationship between capture rate and aquifer permeability strength and magnitude of extraction-return seepage

式中：η为水力控制程度，%；Q为抽出-回渗量，cm3/s；K为含水层渗透系数，cm/s；R为回渗半径，cm.

推导出的定量关系表达式，可以指导不同模拟情景参数范围内污染场地实现水力控制的关键参数的设计工作.但针对实际污染场地应用存在一定局限性，由于单井抽出-回渗循环原位修复技术涉及回渗环节，需要通过水文地质勘察确定含水层渗透系数(K)，对于粉砂、细砂、中砂和粗砂范围内的岩土类型更适用，对于含水层渗透情况差的，如砂质粉土、黏土等岩土类型应用该技术时存在一定局限.此外需要确定地块最大允许抽水量(Q)(保证含水层可以连续抽水而不会疏干，并且回渗水可以使地下水水位保持动态平衡)，在确定了水力控制程度(η)的基础上，可利用该定量关系表达式确定回渗半径(R)和回渗范围.

结果显示，64 种模拟情景的捕获率均大于90%，含水层渗透性越好，抽出-回渗量越大的时候捕获率进一步提升，可以较好地完成水力控制.主要是由于模拟情景为较为稳定的地下水环境，不存在地下水的侧向径流，可以达到较好的水力控制效果.对于天然状态下的地下水环境，侧向径流会受到抽水井的扰动，导致影响区域内存在水头的降低，使原本单向均匀的流场变为不均匀流场，从而影响捕获率.此外模型设置为均质含水层，而实际场地渗透性常为非均质，回渗一般是不均匀的，对捕获率也存在一定影响.因此通过模型结果得到的定量关系表达式存在一定的不确定性，对该问题的解决需要进一步设置具有不同水力梯度和非均质水文地质条件的情景，通过数值模型进一步优化参数关系.

3 结论

a) 不同抽出-回渗量(1、2.5、5、10 cm3/s)条件下，砂柱物理模型抽出-回渗平衡时，抽水井捕获范围分别为3 076、3 143、3 151 和3 162 cm2，均小于回渗范围(3 217.0 cm2)，水力控制程度分别达到95.63%、97.69%、97.93%和98.17%.

b) 影响水力控制程度的关键因素包括水动力条件和水文地质条件，对应参数为渗透系数、抽出-回渗量和回渗范围.当回渗范围和含水层渗透系数一定时，随着抽出-回渗量的增大，捕获范围和捕获率均增大；当回渗范围和抽出-回渗量一定时，随着含水层渗透系数的增大，捕获范围和捕获率均减小；当含水层渗透系数和抽出-回渗量一定时，随着回渗范围的增大，捕获范围和捕获率均减小.

c) 水力控制程度与含水层渗透性强弱、抽出-回渗量大小和回渗半径之间存在定量关系，定量关系表达式可用于指导不同情景和条件下实际污染场地实现水力控制的关键参数的设计.