基于Visual MODFLOW的某工业园区地下水污染溶质迁移预测

华兴国,龙 泉,郑 颖,王可可,林 清,罗曼琳

(四川省生态环境科学研究院,四川 成都 610000)

地下水是水资源的重要组成部分,地下水的合理利用和有效保护,对于支撑经济社会高质量发展和维系良好生态环境具有重要作用。工业园区作为重要的地下水污染源类型,查明地下水污染的迁移情况,对开展地下水污染源头防治,保护地下水环境具有重要意义。

地下水溶质运移模拟是用于模拟和计算地下水中溶质运移的动态特征及其演变方向。溶质运移研究始于上世纪50年代,首先由Amundson和Lapidus对溶质运移模型进行了初步的尝试,提出了一个类似于对流弥散方程(GDE)的数学模型[1],我国于上世纪80年代才开始对溶质运移的研究,起步较晚,初期主要为农业环保部门对农业生产过程中产生的氮、磷和农药等的迁移转化研究。目前溶质运移理论研究已经比较完整的理论体系[2],地下水模型的求解方法主要有解析法、物理模拟法及数值模拟法三种。其中数值模拟法是现今主要的求解方法[3],本文选取某典型工业园区,对其生产过程中,污染物泄漏后在地下水中迁移的情况进行模拟预测,为园区地下水污染防治和源头管控提供科学依据。

1 研究区概况

1.1 地理环境

研究区位于金沙江河谷西侧的谷坡上,地貌上属于剥蚀构造中的台状中山。金沙江从北向南流经园区,整体地势西高东低,倾向金沙江。园区内最低点945 m,最高点1 650 m,最大相对高差705 m。

1.2 地层岩性

研究区出露的地层由新到老主要有第四系全新统(Qh)砂卵石、第三系上新统(N2)灰色页岩、细砂岩、侏罗系下统冯家河组(J1f)紫红色泥岩夹灰色石英砂岩、三叠系上统大菁组(T3bd)灰色砂岩、页岩和砾岩以及岩浆岩(δo2)石英闪长岩。

1.3 水文地质条件

研究区主要包含三大地下水类型,各类地下水分别赋存于不同岩层和岩体中,主要包括松散堆积层孔隙水,碎屑岩孔隙裂隙层间水和基岩裂隙水。松散堆积层孔隙水地层岩性为粉细砂岩、泥岩在研究区北部广泛出露,碎屑岩孔隙裂隙层间水地层岩性为石英砂岩,在研究区南部广泛出露,基岩裂隙水地层岩性为石英闪长岩,主要分布于研究区东侧临近金沙江区域。

研究区以浅部的孔隙、裂隙地下水为主,主要由高处向低处运移,在沟谷两侧或斜坡带上以动态极不稳定的裂隙泉排泄,受地形地貌和岩性的控制,仅经过短途由地势高的丘包向地势低的冲沟、河谷径流,受裂隙展布规律控制,最终向东侧金沙江方向排泄。

2 研究方法

2.1 预测方法

基于资料收集和现场调查,分析并掌握研究区的环境和水文地质特征,建立地下水流动的污染物迁移的数学模型,根据污染源特征确定污染源强及预测参数,建立以 Visual MODFLOW 数值计算水质预测模型,对研究区地下水污染溶质迁移情况进行预测。

2.2 地下水流数学模型

综合研究区的地层岩性、地下水类型、地下水补径排特征、地下水动态变化等水文地质条件及研究区水均衡分析等,在现有资料的基础上,将研究区地下流系统概化成非均质各向异性准三维非稳定地下水流系统,用1式的数学模型表述[5]:

(1)

式中:H(x,y,z,t)标识模拟区任一点(x,y,z)任一时刻t的水头值(m);Ω表示地下水渗流区域;S1为模型的第一类边界;S2为模型的第二类边界;Kxx,Kyy,Kzz分别表示 x,y,z主方向的渗透系数(m/d);w表示源汇项,包括降水入渗补给、蒸发、井的抽水量和泉的排泄量(d-1);μs表示单位贮水率(1/m);H0(x,y,z)表示初始地下水水头函数(m);H1(x,y,z)为第一类边界己知地下水水头函数(m);q(x,y,z,t)为第二类边界己知单位面积流量或单宽流量函数(m3/d·m2),零流量边界或隔水边界 q=0。

2.3 溶质运移数学模型

溶质运移的数学模型[6],如下所示:

(2)

(3)

式中:αijmn为含水层弥散度(m),Vm、Vn分别为m和n方向上的速度分量(m/d);C为含水层中污染物的浓度(mg/L);n为含水层的有效孔隙率;t为时间(d);C’为源汇项中污染物的浓度(mg/L);W为面状源汇项强度(m3/(d·m2));Vi为地下水渗流速度(m/d)。

3 结果与分析

3.1 地下水流数值模型

3.1.1 水文地质概念模型

研究区主要含水层为第三系上统昔格达组(N2x),三叠系大菁组(T3dq)和石英闪长岩裂隙水。研究区内较大的河流为金沙江,位于研究区东侧,地表水由北向南径流。园区位于区域斜坡之上,园区东侧以金沙江为界,其余区域由地表分水岭构成园区所在水文地质单元的边界。地下水顺地形由北西向南东径流,最终排泄至金沙江,金沙江为区域最低基准排泄面。

3.1.2 模拟区边界条件

模型考虑地下水主径流方向、水文地质边界及该含水介质中污染因子的迁移性能,本项目渗流场模拟范围:东侧以项目区最低排泄基准面金沙江为界,北侧、西侧和南侧以地表分水岭为界。模型概化范围共计78.5 km2。模拟区东～西方向作为模型的 x 轴方向,长10 850 m,每108.5 m 划分一个网格,南～北方向作为模型y 轴方向,长9 770 m,每97.7 m 划分一个网格,垂直于 xy 平面向上为模型 z 轴正方向,模拟范围 975～2 040 m。根据含水介质的性质及模型计算需求,垂向上概化为3层。

本次模拟区以评价区最低排泄面金沙江设置为River边界;项目北侧、西侧、南侧地表分水岭设置为Constant Head 边界,金沙江东岸非本次模拟区设置为无效单元格;其余网格为计算单元格。模型网格划分网格见图1。

图1 模型边界设置

3.1.3 模型识别与验证

按照前述建立的数值模型,边界条件和计算参数,以稳定流运行得到的流场作为初始渗流场,见图2。根据模拟结果,受区内地形地貌、含水岩组等条件控制,项目北西侧地下水位较高,由北西侧往南东侧,区域内地下水位呈下降趋势,这与研究区所在区域水文地质条件相符,故利用模型模拟结果所得流场作为评价区初始渗流场基本合理。

图2 初设流场模拟结果(单位：m)

根据现场调查结果等资料,选取沿地下水主径流方向分布的监 测井JC31及JC09对模型地下水渗流场进行校验。各水文点现场水位实测值介于1 134.9～1 163.4.0 m,模拟水位介于1 122～1 222 m,结果位于95%置信区间之内,故利用模型计算所得流场作为项目区初始渗流场基本合理(见图3)。

图3 实测水位与模型计算水位对比图(单位：m)

3.2 溶质运移模拟及结果

3.2.1 模拟预测情景

假定埋地污水管线的连接交汇部位(如焊缝)出现开裂或破损情况,那么泄漏的污染物将会污染地下水。假设泄漏的污染物质全部进入地下水,破裂泄漏的等效孔径按50 mm 计算,将渗漏检测+修复时间共定为150 d,泄漏量计算公式如下。

Q=S×V

(4)

式中:Q为渗漏量(m3/d),S为截面积(m2),V为流速(m/s)经计算,渗漏量为12 780 m3。

泄漏污染物浓度参考《镁、钛工业污染物排放标准》(GB25468-2010)中水污染物排放控制要求中间接排放污染物浓度要求,初始浓度以园区地下水现状浓度为准,选取氨氮和COD为评价因子(见表1)。

表1 废水中各污染物浓度

3.2.2 预测结果

模拟结果显示,管道出现破损后150 d,泄漏点中心氨氮浓度为3.5 mg/L,影响范围19 890.2 m2,运移距离为92.8 m。150 d后泄漏停止,泄漏1 000 d后,即泄漏停止750 d后,污染晕中心浓度为1.6 mg/L,影响范围17 801.4 m2,运移距离为102.2 m。泄漏3 650 d后,污染晕中心浓度为0.6 mg/L,周边地下水氨氮不在出现超标。泄漏点中心COD浓度为16 mg/L,影响范围18 556.9 m2,运移距离为87.5 m。150 d后泄漏停止,泄漏1 000 d后,即泄漏停止750 d后,污染晕中心浓度为9 mg/L,影响范围15 970.7 m2,运移距离为96.4 m。泄漏3 650 d后,污染晕中心浓度为3.5 mg/L,影响范围423.1 m2,运移距离为15.2 m(见图4～图9)。

图4 管道泄漏150 d后氨氮浓度分布图

图6 管道泄漏1 000 d后氨氮浓度分布图

图8 管道泄漏3 650 d后氨氮浓度分布图

4 结语

(1)园区企业管道发生破损后,管道废水通过进入地下水,致使地下水中氨氮和COD超过地下水III类质量标准,氨氮最大超标倍数为6倍,COD最大超标倍数为4.3倍。

(2)废水泄漏后,地下水污染面积最大为19 890.2 m2,最大迁移距离为102.2 m,对应污染因子为氨氮。

(3)管道废水停止泄漏后,地下水中氨氮和COD浓度逐步降低,泄漏停止750 d后,地下水中氨氮最大浓度下降至1.6 mg/L,COD浓度下降至9 mg/L,泄漏停止3 500 d后,地下水中氨氮最大浓度下降至0.6 mg/L,周边地下水氨氮不在出现超标。COD最大浓度下降至3.5 mg/L,影响范围423.1 m2,运移距离为15.2 m。