基于Visual Modflow的某化工项目地下水污染预测

魏庆龙，杨金伟，靳廷朝

（河南省地矿局第一地质环境调查院，河南郑州450045）

在地下水资源的定量分析研究中，建立能较真实反映水文地质结构与流场特点的地下水数值模拟是地下水资源评价与管理的重要研究方法[1-2]，基于Visual Modflow软件的有限差分法可以满足复杂地下水的数值模拟[3]。

1 模拟区水文地质概况

模拟评价区处于豫北黄河冲积平原，不是一个完整的水文地质单元，地下水可以划分为松散岩类孔隙水、半胶结碎屑岩类裂隙水。其含水层组按埋藏条件，可进一步划分为浅层含水层组、中深层含水层组和深层含水层组。将浅层含水层和中深层含水层之间的粘土层当作模型的隔水底板。浅层含水层底板埋深90～120m，第一个连续稳定隔水层之上的上更新统与全新统含水层，浅层水为潜水—微承压水，含水介质为冲洪积及风积的粗、中、细、粉细砂组成，可见3层，单层厚度20～30m，总厚度为60～90m，由南向北，含水层厚度变薄。单井出水量为1800～2400m3/d，渗透系数13～20m/d，为中等富水区。

浅层地下水排泄主要以垂直蒸发为主，近些年来人工开采已成为重要的排泄方式，其次是向下游的径流排泄，以及向中深层地下水的越流补给。大气降水为本区浅层水的主要补给源。引黄渠道的渗漏补给也是本区浅层水的补给源。浅层水化学类型为HCO3·SO4-Na·Mg。

本区浅层地下水动态类型属气象-水文型。区内浅层地下水水位埋深差异较小，一般在3～7m之间，年变幅一般在0～4m左右。水位动态不全受降水的制约，而且又受引黄灌溉渠道水回渗因素影响，年内水位出现2个峰值，与降水和引黄灌溉时间相吻合。

2 地下水水动力场数值模拟

水是溶质运移的载体，地下水流场是溶质运移模拟的基础，在溶质运移模拟前，需先建立模拟区地下水流场模型，地下水流场预测模型为modflow。

2.1 地下水流数学模型

根据模拟区水文地质特征，可将模拟区浅层孔隙水含水层概化为非均质、各向异性、三维非稳定地下水流系统，其运动的数学模型为：

式中：Ω——渗流区域；

Kxx、Kyy、Kzz——渗透系数在x、y和z方向上的分量，假定渗透系数的主轴方向与坐标轴的方向一致，m/d；

h——水位（水头）标高，m；

ε——源汇项，m2/（d·m），ε＜0表示流出地下水系统，ε＞0表示流入地下水系统；

Ss——含水层给水性能参数，取重力给水度；

从图2可以看出，常规治理方案一般需将脱硫装置出口湿烟气从50～55 ℃加热至75～80 ℃(烟道及烟囱温降按5 ℃考虑)，需要较大的热量，现在电厂普遍采用除尘器前或脱硫塔前的高温烟气作为加热热源，高温烟气余热无法再进行节能综合利用，因此，治理方案应尽量减少饱和湿烟气的加热热量，并考虑烟气余热的综合利用效果。

t——时间，d；

h0——含水层的初始水位，m；

Γ2——流量边界；

Γ3——浅层含水层的底板；

ˉ——边界面的法线方向；

q(x，y，z，t)——二类边界的单宽流量，m3/（d·m），流入为正，流出为负。

2.2 地下水流数值模拟的建立

对整个区域模型采用不等距矩形网格剖分，其中厂区进行局部加密剖分，其他区域等距离剖分，剖分为94行133列，共剖分矩形网格单元12502个，其中有效单元11118个，无效单元1384个。模拟区网格剖分见图1，模型结构剖面图见图2。

图1 模拟区网格剖分图

2.3 模型识别

图2 模型结构剖面示意图

模型识别过程通常要在反复修改水文地质参数和调整某些源汇项基础上才能达到较为理想的拟合结果。此模型的识别过程采用的方法称为试估—校正法，属于反求参数的间接方法之一。运行计算程序，可得到这种水文地质概念模型在给定水文地质参数和各均衡项条件下的地下水位时空分布，通过拟合同时期的地下水水位，识别水文地质参数、边界值和其它均衡项，使建立的模型更加符合模拟区的水文地质条件。

根据所掌握的资料，以2013年底流场作为初始流场（见图3），识别期选为2014年1月～2015年6月，2015年7～12月为验证期，应力期以月为单位，共划分为24个应力期。地下水位拟合效果见图4，可见地下水位模拟水位和实测水位误差均小于0.5m，所建立的数学模型达到了模型精度要求，可以利用该模型对研究区进行地下水污染情景预报。

3 地下水污染模拟预测

本次地下水污染模拟过程未考虑污染物在含水层中的吸附、挥发、生物化学反应，模型中各项参数予以保守性考虑。溶质运移模型采用MT3DMS进行求解。

3.1 溶质运移数学模型

地下水中溶质运移的数学模型为：

式中：ρb——介质密度，mg/(d·m)3；

θ——介质孔隙度，无量纲；

C——组分的浓度，mg/L；

t——时间，d；

x、y、z——空间位置坐标，m；

Dij——水动力弥散系数张量，m2/d；

Vi——地下水渗流速度张量，m/d；

Cs——组分的浓度，mg/L。

图3 初始流场

3.2 地下水污染预测情景设定

通过对厂区平面布置、污水收集排放系统等综合分析，确定规划区污水处理站为此次模拟泄露点。根据厂区规划特点，特征污染因子选为COD，地下水污染主要预测非正常工况防渗不正常情景

而非正常工况下，构筑物的最大渗水量取100倍的正常工况的渗水量进行计算。由于厂区有较完备的污染事故监测预警系统，考虑不正常工况持续的时间为2d，非正常工况排除后，污水池渗水量按2L/（m2·d）进行计算。

3.3 场区生产对地下水水质的影响评价

根据非正常工况防渗不正常情景下污染模拟预测结果，依据《地下水质量标准》（GB/T14848-1993）中CODMn的Ⅰ类和Ⅲ类水的要求，特征污染因子在1、5、10、20、30年对地下水的影响范围、超标范围见表1和图5。第1年CODMn浓度最大值分别为12.478mg/L；第30年CODMn浓度最大值分别为82.0544mg/L。第30年CODMn影响范围已到达厂区边界。此外，第27d，CODMn浓度超过地下水Ⅰ类标准，其浓度最大值为1.0288 mg/L，如果以厂区范围作为污染控制范围、以CODMn作为控制性污染因子，则非正常工况排除的时间不超过2d。

图4 地下水位拟合效果图

4 结论

根据研究水文地质资料和地下水观测资料，识别和验证了区域地下水运动水流运动模型，数值模拟模型具有较高的模拟精度。

根据厂区规划特点，特征污染因子选为COD。地下水污染预测情景设定为非正常工况，构筑物的最大渗水量取200L/（m2·d），不正常工况持续的时间为2d，非正常工况排除后，污水池渗水量按2L/（m2·d）进行计算。

若以《地下水质量标准》（GB/T14848-1993）中Ⅰ类和Ⅲ类水CODMn浓度标准值分别作为地下污染影响范围和超标范围，则非正常工况下，第1年CODMn浓度最大值为12.478mg/L；第30年CODMn浓度最大值为82.0544mg/L。第30年CODMn影响范围已到达厂区边界。此外，第27d，CODMn浓度超过地下水Ⅰ类标准，其浓度最大值为1.0288mg/L，如果以厂区范围作为污染控制范围、以CODMn作为控制性污染因子，则非正常工况排除的时间不超过2d。

图5 非正常工况下CODMn影响范围

表1 非正常工况下污染物对地下水的影响范围

[1]薛禹群，谢春红.地下水数值模拟[M].北京:科学出版社，2007.

[2]钱家忠，吴剑锋，朱学愚，等.地下水资源评价与管理数学模型的研究进展[J].科学通报，2011,46（2）:99-103.

[3]张洪霞，宋文.地下水数值模拟的研究现状与展望[J].水利科技与经济，2007,13（11）:794-796.

[4]郭瑞.兴平市地下水动态特征及数值模拟研究[D].西北农林科技大学，2008:5-8.

[5]徐娟花.地下水数值模拟现状[J].工程地质与计算机应用，2012(4):10-14.

[6]陈正国.数值模拟在地下水资源评价中的应用研究[J].黑龙江水利科技，2012,34(12):69-74.

[7]孙从军.地下水数值模拟的研究与应用进展[J].环境工程，2013,31(5):9-17.

[8]孙道林.典型垃圾堆放场对地下水污染模拟及控制技术研究[D].合肥工业大学，2008:23-35.

[9]刘娟，杨国勇，蒋兆华，等.FEFLOW在某市岩溶地下水四氯化碳污染中的应用[J].地下水，2011,33(4):87-91.