应用Visual Modflow软件对宿州市某制药厂地下水污染物运移过程的模拟及污染预测

马国明

（安徽省地质实验研究所，安徽合肥230041）

0 引言

随着社会经济的高速发展，人口快速增长和城市规模的扩大对浅层地下水资源的污染现象日趋严重，因此，迫切需要有效防治地下水污染的方法和技术监管手段，以便更好地保护和利用地下水资源。目前地下水污染数值模拟技术正在迅速发展[1～3]，多种地下水数值模拟软件(GMS、Modflow 等)被应用到地下水调查评价工作当中。目前国内地下水模拟更加注重于客观水文地质条件认识的真实性，数值模型结构概化的合理性[4～5]。本文旨在使用Visual Modflow 软件对宿州市某制药厂附近浅层地下水的溶质迁移污染问题进行数值模拟，探讨分析泄漏污染物随地下水的运移对周边地下水环境造成的影响，探索监测、防治地下水污染的新方法。

1 研究区概况

1.1 地质条件

研究区地处淮北冲积平原的中部，地貌条件简单，地面标高在23.00～26.70m，地形平坦，总体趋势北和西北部略高，向东南和南缓倾。区域地层属华北地层区淮河地层分区淮北地层小区。前第四纪地层除缺失太古界、元古界、上奥陶统至下石炭统、三叠系、侏罗系中-下统、白垩系外，自寒武系至第四系其它地层均有发育。本研究区地表出露地层为第四系，松散层包括新近系(N)和第四系(Q)，厚度在100～300m，由北向南逐渐变薄。下伏基岩主要为奥陶系下统萧县组、马家沟组和中统老虎山组，石炭系上统本溪组、太原组，二叠系上统石千峰组，古近系古新统双浮组。

1.2 水文地质条件

根据区内地下水含水介质的不同和赋存条件的差异，本区地下水类型可划分为松散岩类孔隙水、碎屑岩类孔隙-裂隙水和碳酸盐岩类岩溶裂隙水三种。其中松散岩类孔隙含水层自上而下可划分为三个含水岩组：①松散岩类孔隙含水岩组：含水层主要由中晚更新世地层所组成，岩性为粉细砂、亚砂土、亚黏土、黏土等，一般埋深20～30m，含水层厚度5～25m，底板埋深小于50m。②中深层松散岩类孔隙含水岩组：分布广泛，地层主要为中、下更新统。含水层埋深一般在35～110m之间，发育4～5层砂，岩性以分选性较好的细砂为主，少量中细砂。单层厚5～15m，累计厚度35～55m。砂层底板埋深一般为50～100m。水力性质为承压水，水质优良，是城乡主要供水水源。③深层松散岩类孔隙含水岩组：广泛分布，地层多为下更新统和新近系，岩性以细砂、中砂为主，少量粗砂，多为半胶结状。埋深一般大于100m，底板埋深150～200m，砂层单层厚度一般在15～20m，累计厚度一般10～40m。该层水量较贫乏，开采量较小。

研究区浅层松散岩类孔隙含水岩组补给方式主要为降水入渗补给、地表水下渗补给和灌溉回渗补给，其中，降水入渗补给是研究区地下水的主要补给来源，除此之外，还有少量河流和农田灌溉下渗补给。

区域上浅层松散岩类孔隙水径流总体上是由西北流向东南。有人工开采的地域，浅层松散岩类孔隙水流场受地形和开采双重影响，局部会形成漏斗。浅层松散岩类孔隙水的排泄主要为蒸发排泄，农村人畜生活用水和灌溉用水的开采。中深层水除了侧向径流排泄外，还有城市供水的集中开采。

2 地下水数值模拟

2.1 模拟软件的选取

Visual Modflow 软件是近几年成功应用于水文地质学领域的一款计算机软件，技术上成熟，功能强大，使用方便，正在被越来越多的水文地质、环境地质工作者所接受。

该软件主要由MODFLOW、MT3DMS、MODPATH 和ZONEBUDGET 四个功能模块组成，可分别用来模拟地下水流的运动情况；模拟污染物粒子的运动轨迹；以及计算给定区域的总水量及其与周围区域的水量交换情况，特别是对于分析特定的水量变化情况非常有用。本次实际案例即是应用Visual Modflow 软件对某制药厂附近地下水的溶质迁移问题进行数值模拟。本次模拟以氨水储罐泄漏导致部分氨氮通过防渗破损层进入含水层为案例。

2.2 水文地质概念模型的建立

数值模拟结果的正确与否取决于水文地质概念模型的合理程度。建立概念模型之前，需要充分了解实际的水文地质条件，认真分析地下水系统，统计组织相关的测量试验数据。通过对实地水文地质条件的概化，确定模型的基本要素，为建立地下水数值模型奠定基础[6～8]。

研究区地下水渗流系统按岩层的渗透性自上而下划分为第一含水层组（浅层）、第一隔水层、第二含水层组（中深层）和第二隔水层（图1）。以黏土层作为模型隔水层。研究区调查显示，对地下水资源的开采量是相对稳定的，本次假设其为稳定流状态，因此概化成非均质各向同性条件下的三维稳定流地下水系统。

图1 水文地质概化模型示意图Figure 1. Schematic diagram of hydrogeological conceptionalization model

2.3 数学模型

在水文地质概化模型的基础上，综合以上水文地质基本要素及特征，通过对研究区的水均衡分析，利用已有资料，将研究区含水层概化为三维稳定、非均质各向同性的地下水流系统，数学模型如下：

式中：

Ω—为研究区地下水渗流范围；

H—为含水层水头值(m)；

S1—为地下水渗流区第一类边界；

S2—为地下水渗流区第二类边界；

kxx，kyy，kzz—分别表示纵向、横向和垂向的渗透系数(m/d)；

w—为地下水的源汇项(m3/d)；

H0(x，y，z)—为初始的水头函数(m)；

H1(x，y，z)—为已知地下水第一类水头边界函数(m)；

q(x，y，z)—为第二类单位面积流量边界函数(m3/d)；

n—为第二类S2边界的外法线方向。

2.4 模型检验与参数选取

模型建立之后，需要确定各项水文地质参数来检验模型的正确性。这些参数包括作为源汇项部分的用于计算各种地下水补排资源量的参数，如降水渗入系数、灌溉回渗系数和河流侧向渗漏系数，以及含水层的水文地质参数：渗透系数、给水度、导系数、越流系数及释水系数等。

根据实际抽水试验观测资料、水化学分析测试资料、经验值资料，收集的历史资料以及试验获得的各层水文地质参数，对模型含水层赋值反复模拟计算，通过计算水位和实际水位拟合分析，反复调参，最终选取最佳拟合度方案确定含水层各项参数（表1）。

表1 含水层水文地质参数分区表Table 1. Hydrogeological parameters of aquifer

2.5 模型识别

根据水文地质模型所建立的数学模型，能否确切地反映计算区的实际水文地质条件，需要进行验证，这个过程就是模型识别。本次模拟识别先是选取研究区已有的抽水试验结果作为本次模型的水文地质参数初始值，然后用试错法调整水文地质参数，拟合研究区地下水各含水层的流场，直到能反映实际流场特征为止。

通过调整水文地质参数和拟合研究区流场，模型校正工作取得了比较好的效果，各含水层的实际流场和拟合流场对比见图2和图3所示。

图2 浅层含水层流场拟合图Figure 2. Shallow aquifer flow field fitting

图3 中深层含水层流场拟合图Figure 3. Middle and deep aquifer flow field fitting

2.6 地下水溶质运移模拟

2.6.1 溶质运移数学模型

（1）控制方程

根据本区实际水文地质条件，本次采用的地下水溶质迁移模型是一个基于三维水流模型下的三维溶质弥散模型，假设溶液密度保持不变，不考虑吸附作用和动力反应过程，吸附相和溶解相的速率相等。根据上述假设，本次研究的地下水中溶质迁移的数学模型表述如下：

式中：

C—地下水中的各组分溶解相浓度，ML-3；

θ—含水层孔隙度，无量纲；

t—时间，T；

xi—沿坐标系轴向的距离，L；

Dij—地下水动力弥散系数，L2T-1；

vi—地下水的平均实际流速，LT-1；

qs—单位体积的地下水流量，L3T-1；

Cs—源汇项中的组分浓度，ML-3；

∑Rn—化学反应条件浓度随时间变化，ML-3T-1。

在模型拟合验证后，本次地下水稳定流条件下三维溶质运移模型直接利用Visual Modflow 中的MT3DMS 模块进行不同时段的预测计算，模型的边界及初始条件设置如下：

（2）初始条件

本次溶质模型的污染源概化为浓度边界，初始浓度设定为C0，具体表述为：

式中：Γ2—为通量边界；

fi(x,y,t)—代表边界弥散通量的已知函数。

2.6.2 溶质运移参数的确定

根据相关经验结果，本次研究的含水层为孔隙水，对应的弥散度位于1～10之间，因此本着最不利原则，本次模拟弥散度取10(图4)。

图4 孔隙介质数值模型的LgaL-LgLs图Figure 4. LgαL-LgLs graph for a numerical model of porous media

2.6.3 模型条件的概化

根据实地客观条件，污染源位置按实际设计概化。模型将上述情形的污染源以面源形式设定浓度边界。由于资料所限，在模拟污染物扩散时，主要考虑对流、弥散的作用，不考虑吸附作用、化学反应等其它因素。

为了分析该工厂在氨水储罐泄露的情形下对厂区以及下游地下水的影响范围和程度，利用上述拟合好的流场水流模型，预测不同的时间段污染物的迁移扩散情况。

图5 预测氨氮泄漏污染羽范围变化图Figure 5. Variation diagram of the predicted range of ammonia and nitrogen leakage plum

表2 污水渗漏事故氨氮对地下水水质的影响情况Table 2. The influence of ammonia nitrogen on groundwater quality after a sewage leakage accident

3 预测结果分析

本次模拟是基于非正常工况下氨氮泄漏进行的模拟预测，本次模拟污染物迁移扩散过程未考虑土壤的吸附作用，以及污染物与土壤和地下水中物质的相互化学反应，在泄漏量和泄露浓度一定的情况下，模型中污染物扩散在时间和空间上显示一定规律。

(1)泄漏时间：假设非正常工况下防渗层破损，氨短时间泄漏。

(2)预测结果：将污染源输入模型，模拟预测发生渗漏事故后100 天、1000 天、10 年和20 年污染羽的变化情况(表2)。参照《地下水质量标准》氨氮浓度，即0.5mg/L。

由表2和图5可见，非正常工况下情景发生后，如果罐区发生泄漏，渗漏的污染物会对第一层(浅层)含水层的水质造成一定影响，第二层(中深层)含水层基本无影响。污染物泄漏时间较短，泄漏过程中，污染物在水动力条件下不断由东南向西北方向运移，并且污染羽范围逐渐缩小，泄漏第100天后，污染晕水平运移34.3m；在运行第1000 天后，污染晕水平运移108.0m；泄漏后3650 天，污染晕水平运移175.5m；泄漏后7300 天，污染晕水平运移265.0m。说明发生事故后，如果防渗区破损，泄漏的氨水未及时清理，泄漏渗入含水层，则污染物会对下游浅层含水层造成一定影响。因此，非正常工况下，污染物有限时间的泄漏之后，需要及时采取处理措施，防止污染物泄漏下渗，污染浅层地下水。

4 结语

根据Visual Modflow软件的模拟分析结果，研究区在不合理排放或意外事故状态泄漏的污染物，一旦泄漏渗入含水层，将对浅层地下水水质造成影响，甚至污染，其污染程度取决于污染源的浓度。因此营运期建设单位需要加强罐体、围堰以及防渗层的维护保养，一旦发生泄漏事故，应及时清理并检查泄漏区防渗层，防止污染物进入包气带和含水层，保护好有限的地下水资源。

该模型的建立和模拟分析方法可对污染物意外泄漏后随地下水流运移对周边地下水环境造成的影响进行预测和实时评估，对探索监测浅层地下水污染的方法具有一定意义。