基于Visual Modflow的山区某再生铜地下水污染模拟

李 巍，卢玉东，卢阳春，郭 雯，梁少欢

(1.长安大学 环境科学与工程学院，西安 710054；2.旱区地下水文与生态效应教育部重点实验室，西安 710054)

0 引 言

我国人均水资源量很低，北方地区65%的生活用水都来自于地下水[1],但地下水污染问题却愈发严重,造成的水资源短缺给我国的经济发展带来了严重影响。铜工厂排放的废水既会造成资源损失，也会污染水环境，污染物中的重金属元素还会严重威胁到当地居民的饮水安全和人体健康[2]。因此探明地下水污染物迁移的规律，对地下水资源的保护、居民饮水质量的保证有着非常重大的意义。Visual Modflow是专门用于孔隙介质中的三维有限差分地下水流数值模拟软件[3]，国内外有很多人利用Visual Modflow对地下水溶质运移进行模拟[4]。马志飞[5]模拟了六价铬在某危险废物填埋场的特点，结果表明增加防渗墙和地表硬化能有效防止污染的扩散；尉鹏翔[6]研究了北京某污染厂区的苯和萘长时间渗漏污染羽的特点，发现污染物在迁移时其长轴与水流方向一致；马从安等[7]对某露天矿进行地下水质点示踪来判断矿坑涌水的水源及通道。国外的Laura K. Lautz[8]模拟了交错带地表水、地下水的汇流。山区地下水往往受地形地貌、岩性、地质构造条件等控制。表现为山岭与山坡地带主要为地下水补给径流区，河谷地带主要为排泄区，其特点为地下水位变化大，水流沿着沟道向下游流动。本文以陕南某山区再生铜为例，通过Visual Modflow软件建立地下水流模型、溶质运移模型，研究地下水溶质运移问题。模拟结果为地下水评价提供科学指导依据。

1 研究区概况

1.1 基本情况

研究区位于陕西省南部地区，在洛南县城西19 km处，属低中山区，高程在1 000～1 500 m之间，沟谷切深200～500 m，山势较低，斜坡在25°～40°之间，基岩大多被坡残积土覆盖，植被覆盖率30%～60%。

研究区处属暖温带半湿润季节性气候，年平均气温11.1 ℃，最热月份(7月)平均气温23.1 ℃，最低月份(1月)平均气温-1.9 ℃，年平均降水量为754.8 mm左右，蒸发量为1 034 mm左右。

1.2 水文地质

研究区丘陵山区是基岩(包括碳酸盐岩)强风化带潜水，顺地势往沟谷排泄，强风化带潜水以下受相对隔水的泥岩或裂隙不发育的岩层阻隔，为基岩裂隙承压水，径流迅速；沟谷区地势较为平坦，有薄层第四系分布，但厚度变化大，厚度小的地段仅下伏基岩强风化带赋存地下水，厚度较大的地段第四系松散层和下伏强风化带之间无稳定隔水层存在，共同构成具有双层含水介质的统一潜水含水岩组：

(1)松散岩类孔隙水。上部岩性为冲洪积物含碎石黏土含有少量碎石，透水性较强，不具备储水条件，仅零星含水，厚度4.3～7.3 m。

下部岩性为碎石土，碎石主要由全风化灰岩组成，充填物为粉质黏土和细砂，厚度4.2～6.1m。含水层富水性差，水量相对贫乏，且不稳定。

(2)灰岩风化带孔隙裂隙水。灰岩浅层风化带透水性强，厚度3.1～4.8 m，沟谷两侧含水层厚度大、渗透性能好，水量相对丰富。较大支沟含水层虽渗透性好，但由于含水层厚度薄，地下水赋存空间有限，地下水量相对较小，含水层富水性属低—中等。

浅层风化带基岩与未风化的深层基岩富水性有明显差别。风化带下部的灰岩，岩体完整性好，裂隙相对不太发育，且密闭性较好，富水性及导水性较差，构成了含水层底板。区内沟道无长年水流，降水以潜流形式补给第四系松散岩层地下水。地下潜水接受补给后，往沟谷下游排泄。

2 原理与方法

2.1 地下水流模型

研究区的含水层主要为第四系松散孔隙潜水、基岩裂隙水，整个研究区概化为单层潜水。将其概化为均质、连续分布的含水层。将研究区概化成非均质、各向异性、三维稳定流的地下水系统概念模型[9]。其运动方程可描述为：

(1)

H(x,y,z,t)|t=0=H0(x,y,z)(x,y,z)∈Ω

(2)

H(x,y,z,t)(x,y,z)∈A1=H1(x,y,z,t)(x,y,z)∈A1t>0

(3)

(4)

式中：H为水头；H0为含水层初始水头；H1为第一类边界水头；q1为第二类边界单位面积流量；Kxx、Kyy、Kzz分别为渗透系数在x、y、z上的分量；W为降水入渗补给强度；Ss为储水系数；t为时间；Ω为研究区范围；A1、A2分别为第一类、第二类边界条件。

2.2 溶质运移模型

不考虑污染物在含水层中的吸附、交换、挥发、生物化学反应，地下水中溶质运移的数学模型可表示为：

(5)

(6)

式中：αijmn为含水层的弥散度；Vm，Vn为分别为m和n方向上的速度分量；|v|为速度模；C为模拟污染质的浓度；ne为有效孔隙度；t为时间；W为源汇单位面积上的通量；Vi为渗流速度；C′为源汇的污染质浓度。通过工程分析，本文主要研究重金属污染物Pb2+、Cu2+地下水中的运移。

2.3 模拟范围及边界条件概化

本次模拟范围是洛南县卫东镇香沟沟谷内，面积约3.41 km2。厂区西、北、南以山为界，东南部约2 km处有一水库，水库东侧为下游村庄。

厂区西、北、南以山为界，将边界概化为山前侧向补给边界，东部以实测等水位线为界，将其概化为定水头边界。上边界：主要接受大气降水补给，概化为潜水面边界。下边界：等效定义为相对隔水边界[10,11]。

2.4 模型的识别验证

模型识别与验正是建立正确反映实际水文地质条件的数值模型的关键步骤之一，数值模拟工作的工作量主要集中在这一步骤。本次工作以第四系松散孔隙潜水、基岩裂隙水含水层的9观测孔(1个水文钻探孔和8个灌溉水井)的地下水动态资料对数值模型进行识别调试。

选择2017年9月观测水位代入模型，通过对水文地质参数的分区、数值进行反复调整，选取观测孔实际值与计算值拟合最好的情况，计算出模拟区地下水流场，再通过2017年12月观测水位进行了验证比较。调参后计算流场与实测流场吻合较好，等水位线基本一致，地下水动态过程曲线与实际监测曲线基本相似。其实测监测井的水位与计算流场水位对比如表1、图1所示。

表1 监测井的水位观测值和计算值比较 m

图1 模拟区含水层参数分区

其最大水位差值最大不超过2 m，一般在1 m左右结果表明地下水动态拟合情况较好。说明含水层结构、边界条件的概化、水文地质参数的选取是合理的，能够反映研究区水文地质条件的实际情况，可用来对未来地下水流场及溶质运移进行模拟预测。

根据当地的水文地质资料，溶质运移的纵向弥散度取10 m，降雨入渗系数取0.2，潜水蒸发系数取经验值0.15。通过模型的识别验证，最终得到模拟区模型的渗透系数的分区见表2、图2。

表2 渗透系数分区及取值

图2 模拟区含水层参数分区

3 结果与分析

通过工程分析，确定该工厂主要污染物为Pb2+、 Cu2+，Pb2+的最高浓度为1.5 mg/L，Cu2+的最高浓度为3 mg/L，位于工厂内初期雨水池处。由于初期雨水池对初期雨水进行收集处理，降雨、存水不连续，根据工程特点，地下水影响预测采用示踪剂短时注入模式，注入时间为30 d。

本项目初期雨水池尺寸为20 m×8 m×2 m，浸润面积为160(底面面积)+112×80%(侧面积的80%，蓄水设计高度为1.6 m)=249.6 m2。根据《地下工程防水技术规范》四级防水工程标准，验收时平均渗漏量2 L/(m2·d)为工程质量合格，则初期雨水池渗水量为499.2 L/d。由于下游敏感点距离污染源2 km，因此只对非正常工况下的渗漏进行模拟预测。非正常工况下，渗水量按正常工况下初期雨水池的10倍计算，渗水量为4 992 L/d，其中Pb2+、 Cu2+渗漏量见表3。根据己建立的地下水流及溶质运移模型，对此工厂的主要污染物Pb2+、Cu2+进行30 d、100 d、1 000 d的渗漏预测。

表3 非正常工况下Pb2+、Cu2+的渗漏量

根据初期雨水池的功能特点及下游设置的检漏井监测频率，超标范围为《地下水环境质量标准》(GB/T 14848-2017)Ⅲ类标准限值，影响范围为Pb2+、Cu2+污染因子的最低检出限值，对Pb2+、Cu2+进行不同时限渗漏情况的预测。

初期雨水池渗漏Pb2+对地下水的污染预测结果见表4及图3、4。

图3 Pb2+渗漏30 d后

图4 Pb2+渗漏40 d后

根据Pb2+的预测结果显示，渗漏30 d时，影响范围为下游112 m，超标面积为3 192 m2，下游的检漏井监测到污染物超标，清理初期雨水池，检查并修复可能的裂缝，此后仅已经渗入地下水体中的污染物继续扩散，渗漏40 d时，影响范围为下游62 m，超标面积为1 327 m2；按照污染物浓度不超过《地下水环境质量标准》(GB/T 14848-2017)Ⅲ类标准限值要求，按照短时注入模式，可见污染晕沿着水流方向运移，浓度降低并不断减小，在58 d时运移至距水池30 m处，污染物呈点状，58 d以后未检出。

表4 初期雨水池预测因子Pb2+渗漏对地下水影响范围

初期雨水池渗漏Cu2+对地下水的污染预测结果见表5及图5-7。

表5 初期雨水池预测因子Cu2+渗漏对地下水影响范围

图5 Cu2+渗漏30 d后

图6 Cu2+渗漏100 d后

图7 Cu2+渗漏200 d后

根据Cu2+的预测结果显示，渗漏30 d时，影响范围为下游88 m，两侧47 m，超标面积为0 m2；渗漏100 d时，影响范围为下游168 m，两侧51 m，超标面积为0；渗漏200 d时，下游影响距离为87 m，两侧40 m，超标范围为0；污染物未超标。按照污染物浓度不超过《地下水环境质量标准》(GB/T 14848-2017)Ⅲ类标准限值要求， Cu2+污染物在切断污染源后影响范围继续扩大，污染晕沿着水流向下游移动，但浓度降低，污染晕迁移在282天时运移至距水池76 m处，污染物呈点状，282 d以后未检出。

Pb2+和Cu2+类污染物由于其浓度不同，检出、超标限值不同，在时间上表现出了污染程度的不同。Pb2+在切断污染源后，在100 d时已经未检出；Cu2+在切断污染源后，影响范围扩大，浓度降低，在1 000 d时未检出。

4 结 论

本文通过Visual Modflow对地形复杂的山区进行地质建模，且污染模式采用更符合实际情况的30 d短期注入模式对地下水污染进行模拟，有如下结论：

(1)利用Visual Modflow对陕南某山区建立水文地质模型，溶质运移模型。通过观测水位及参数校正，使地下水流模型拟合的较好。通过溶质运移模型来研究污染物的迁移规律，为该区域地下水资源评价提供依据。

(2)污染物泄露后，其运移轨迹沿着地下水流西北—东南方向，由于模拟区处于沟谷内，水力坡度较大，西北方向的扩散范围远小东南方向。且污染晕的形状在时间上由一开始的圆形逐渐变化为羽状直至消失。

(3)Pb2+污染物在地下水而迁移中，30 d短时注入时，超标污染晕逐渐增大，在切断污染源后，污染晕沿着水流方向运移并逐渐消失；Cu2+污染物由于超标限值较大，检出限值较低，在切断污染源后，其影响范围还会继续增大并且向下游运移，方向基本与水流方向一致，浓度会随着时间而减小，直至消失。

(4)由于场区距离下游水库1.87 km，距下游水井2.32～2.94 km，按Pb2+污染物短时注入30 d进行计算，污染晕超标范围距离水库为1.75 km，距离下游水井2.20～2.82 km，故预测时段内，污染物质不会对下游水库、水井产生影响。

□