基于Visual Modflow 软件的地下水污染物运移模拟研究

陈 峰，吴 涛

（中国煤炭地质总局 水文地质工程地质环境地质勘查院，河北 邯郸 056004）

0 引 言

水是人类生存不能缺少的资源，随着国民经济的迅速发展和人民生活水平的不断提高，人们对水资源的需求量也不断增长。地下水是工业、农业和人类等生活的重要部分，根据调查，我国大中城市中地下水资源均受到一定的污染。

近年来，人们对环境保护的意识不断增强，地下水环境质量也越来越受到重视。对一些可能会造成地下水污染的企业，防范措施要求也逐渐提高，但仍不能排除特殊情况下污染物泄露对地下水造成污染。应用科学有效的计算方法模拟污染物泄露的影响范围及程度，对于地下水污染的防治有着极其重要的意义。目前地下水系统数值模拟方法主要有有限差分法（FDM）、有限单元法（FEM）、边界元法（BEM） 和有限分析法（FAM） 等[1]，可以用来模拟地下水位下降、地下水位过量开采和地下水污染等问题。随着地下水科学与计算机技术的快速发展，地下水进行水流和水质运移采用数值模拟进行预测研究也得到了广泛的应用[2]。此次利用Visual Modflow软件来模拟研究区地下水污染物运移情况，为研究区地下水污染物的防治提供科学依据[3]。

1 Visual Modflow 软件

Visual Modflow 软件是加拿大Waterloo 水文地质公司在美国地质勘探局原MODFLOW 软件基础上应用现代可视化技术研发成的，该软件是能够建立三维地下水流动和模拟污染物迁移模型的软件，于1994 年8 月在国际上公开发行[12-13]，截止目前，该软件的三维地下水流模拟在国际上也得到认可。地下水一旦遭到破坏，由于地下水中的污染物降解和迁移比较缓慢，地下水一旦污染水化学就很难恢复，采用Visual Modflow 软件在地下水中建立污染物运移模型，来模拟地下水的运移过程。

2 研究区概况

2.1 研究区特征

研究区为河北省馆陶县工业园区，地貌特征属平原型，地势比较平坦，区域位于冲湖积平原区水文地质单元。研究区地下水赋存于第四系松散岩类孔隙含水层中，含水层岩性主要有细砂、粉砂、粉细砂。研究区范围内，浅层水底板一般埋深在40 ～55 m，含水层单位涌水量一般1.2 ～2.7 L/s·m，渗透系数5.12 m/d，为潜水。

该区地下水补给主要有大气降水的入渗补给、地下水侧向径流补给。地下水由西南流向东北，水力坡度在1/1 000 左右。地下水的排泄方式主要为地下水径流排泄和人工开采等。

2.2 地下含水层划分及特征

研究区的地下水主要赋存于第四系多层结构的松散岩层中，以第四系地层为基础，以水文地质条件为依据，自上而下将第四系划分4 个含水组，即Ⅰ、Ⅱ、Ⅲ、Ⅳ含水组。第Ⅰ含水组大致相当于全新统（Q）4、第Ⅱ含水组大致相当于上更新统（Q）3，第Ⅲ含水组大致相当于中更新统（Q）2，第Ⅳ含水组大致相当于下更新统（Q）1。

（1） 第Ⅰ含水组。

含水组底板埋深50 ～110 m，为一套冲积洪积—冲积湖积沉积物。含水层岩性以细砂为主，含水层厚度26 ～60 m，单层厚度一般5 ～l0 m，单位涌水量一般为5 ～17 m3/h。含水组上段为潜水，下段由微承压水变成承压水，水位埋深26.1 ～29.4 m。馆陶北县东北局部地区下部为咸水体外，其余地区矿化度0.9 ～2.0 g/L，为淡水。

（2） 第Ⅱ含水组。

含水组底界埋深200 ～230 m，为一套冲积湖积沉积物。岩性主要为细砂和粉细砂，东南部为中砂含砾石，颗粒稍粗且较厚，含水层厚度为28.2 ～78 m，单层厚度最大可达到20 m。为承压水，单位涌水量为6 ～l0 m3/h·m，水位埋深为29.5 ～30.1 m。主要为咸水，矿化度一般大于3 g/L。

（3） 第Ⅲ含水组。

含水组底板埋深300 ～320 m，为冲积湖积沉积物，岩性主要为细砂，厚度20 ～50 m，北薄南厚，为承压水。基本全为淡水，矿化度小于2 g/L。

（4） 第Ⅳ含水组。

含水组岩性冲积湖积和冰水沉积物，为承压水，岩性主要为细中砂，单位涌水量一般5 m3/h·m左右。全为淡水，矿化度小于1 g/L。

2.3 地下含动态特征

从浅层淡水长期观测资料分析，水位变化属降水入渗—人工开采型，年内水位具有明显的季节性变化，影响水位变化的主要因素是降水量大小，每年l ～3 月份为水位平缓期，是年内最高水位；3 ～6 月中旬为农业开采期，且降雨量少，水位一般下降最大，6 月中旬为年内最低水位，7 ～9 月份为降雨期，水位逐渐回升，至下一年1 月份达到年最高水位，年变幅0.2 ～3.4 m。深层淡水较浅层淡水具有滞后效应，滞后期l ～2 个月，年变幅2 ～4 m。

3 地下水污染物运移模型建立

3.1 地下水流数值模型的建立

3.1.1 水文地质概念模型

水文地质概念模型是对水文地质条件的简化，对地下水系统的科学概化，是把含水体的边界性质、地层内部的结构、地下水的渗透性能、地下水的水力特征、地下水的补给情况和排泄情况等条件概化，为了更方便建立数学和物理模拟的基本模式[4-7]。建立研究区的水文地质概念模型是进行预测评价的第一步。

首先对研究区水文地质条件进行概化处理，研究区内水文地质条件简单，为潜水含水层。由于研究区范围只是水文地质单元的一部分，模型边界为人为划定边界，研究区西南边界和东北边界以地下水等水位线为界，按一类边界处理；西北边界、东南边界基本与地下水等水位线垂直，按隔水边界处理；在垂向上，浅层含水层自由水面为系统的上边界，通过该边界，浅层水与系统外发生垂向交换，比如大气降水的入渗补给；以含水层底板作为模型的下边界。

3.1.2 数学模型的建立

通过对水文地质概念模型的分析，依据渗流连续性方程和达西定律，建立研究区与地下水系统水文地质概念模型相对应的二维稳定流数学模型[8-11]：

式中：H为地下水水头，m；K为渗透系数，m/d；H0（x、y） 为开始地下水水头函数，m；H1（x、y）为第一类边界地下水水头函数，m；q（x、y） 为含水层二类边界单位面积过水断面补给流量函数，m2/d；ε 为源汇项强度（包括开采强度等），m/d；Ω 为渗流区域；B1为水头已知边界，第一类边界；B2为流量已知边界，第二类边界；n为渗流区边界的单位外法线方向。

3.1.3 网格剖分

对研究区进行网络剖分，网格剖分在水平方向上用正交网格剖分法，将模型剖分成80×80 的单元格，在重点研究区域进行了细化剖分，在原剖分基础上进行了2 倍加密，如图1 所示。

图1 模型剖分及边界条件示意Fig.1 Model division and boundary conditions

3.1.4 污染物质的确定

根据现场调查分析得出，研究区的污染不是一个明显的污染源，而是有存在多个次要或者一个主要污染源，分析得出评价区内污染物COD 负荷比73.56%、氨氮负荷比26.44%，COD 为主要污染物。此次利用Visual Modflow 软件对这两种污染物进行数字模拟。

3.1.5 水文地质参数的选取

为了更加准确地表达出该区水文地质情况，模型中参数的确定主要依据研究区已有的水文地质勘查成果，并结合研究区的区域水文地质资料和多年的经验值，确定了含水层参数值，见表1。

表1 水文地质模型参数Table 1 Parameters of hydrogeological model

3.1.6 源汇项处理

（1） 大气降雨入渗补给量。

降雨入渗补给量是指大气降雨渗入到土壤并在重力作用下渗透补给地下水的水量。年降雨入渗补给量的计算公式为：

式中：Pr为年降雨入渗补给量；P为年降雨量；a为年降雨入渗补给系数；F为计算区面积。

研究区多年平均降雨量为445 mm/a，降雨入渗系数参照研究区内水文地质资料，取值为0.18。

（2） 侧向补给。

侧向补给量用达西公式计算，公式如下：

式中：Q为侧向补给量，m3/d；K为渗透系数，m/d；D为剖面宽度，m；M为含水层厚度，m；I为垂直于剖面的水力坡度。在模型中，根据补给边界的以上参数自动计算边界的流入量。

（3） 侧向排泄量。

侧向排泄量用达西公式计算，公式如下：

式中：Q为侧向排泄量，m3/d；K为渗透系数，m/d；D为剖面宽度，m；M为含水层厚度，m；I为垂直于剖面的水力坡度。

（4） 蒸发。

蒸发是指潜水（其中埋深小于4 m 时） 在毛细管力的作用下由下向上运动，最终以参加陆面蒸散发形式散逸到大气中的水分损失量。此次研究区内浅层地下水的埋深均超过了4 m，因此此次地下水蒸发量按零计。

3.1.7 模型的识别及验证

该阶段是模拟过程中的重要部分，先对模型进行识别，给出每个参数的范围值，采用软件自动进行模拟，期间也经过手动进行调整，如果实际与计算水位相差较大，然后根据变化情况在进行调试，直到两者拟合较好为止。

根据研究区地下水分布规律、地下水动态特征以及收集的地下水动态资料，将识别期定为2018年5 月。地下水流场拟合如图2 所示。由图2 可知，在识别期实际地下水流场与模拟的地下水流场基本一致，可以认为实际值与模拟值拟合情况较好。通过以上分析，此次在研究区建立的模拟地下水模型基本符合研究区的水文地质条件。

图2 地下水流场拟合图Fig.2 Groundwater flowfield fitting

3.2 地下水溶质运移数值模型建立

（1） 数学模型的建立。

溶质运移的二维水动力弥散方程的数学模型如下[8-10]：

式中：右端前二项为弥散项，后二项为对流项，最后一项为由于化学反应或吸附解析所产生的溶质的增量；D为弥散系数；μx，μy为x、y方向的实际水流速度；c为溶质浓度：ML-3；Ω 为溶质渗流的区域，量纲：L2；c0为初始浓度：ML-3。

（2） 弥散度的确定。

研究污染物在地下水中迁移转化规律的最重要参数之一是弥散度，反映渗流系统弥散特征的一个综合参数是弥散系数D，在忽略其分子的扩散时，它是介质弥散度仅和孔隙流速V的函数。在地下水溶质运移方程中，表征含水层介质弥散特征的参数是水动力弥散系数，它可表示为：

式中：αL，αT分别为纵向和横向孔隙尺度弥散度，是仅与介质特性有关的参数。结合区域水文地质条件特征，确定潜水含水层纵向弥散度取值0.2 m。

（3） 模拟时段设定。

此次预测模拟工作以5 500 d 为模拟时间，研究污染物浓度时空变化过程与规律，更好的评价厂区污染物泄露后对地下水环境可能造成的直接影响和间接危害提供依据。

4 地下水污染物运移模拟预测

模拟预测根据污染风险分析的情景设计，在选定优先控制污染物的基础上，分别对地下水污染物在不同时段的运移距离和影响范围进行模拟预测。

根据对研究区模拟计算地下水的水质情况和污染源的分布及类型，选取对地下水环境质量影响较大的污染物作为污染物溶质运移对象，通过标准指数法排序，确定了本次污染物，研究区污染物收集池中未200 mg/L。非正常工况下，假设厂区污水收集池防渗层发生破损，则确定厂区污水收集池为模拟泄漏点。研究区地下水现状监测指数0.8 mg/L，以此数据为边界进行计算。污染物发生泄露后，运移模拟情况如图3 所示。预测结果表明，在5 000 d 内污染物羽状污染范围最大为0.2 km2，最大运移距离436 m，在5 000 d时污染物浓度达到边界值，对下游村庄不会造成影响。

图3 污染物运移模拟图Fig.3 Simulation of pollutant transport

5 结 论

此次模拟计算是依据研究区内地下水的水化学现状、污染源的分布及类型，选择对地下水环境质量影响较大的污染组分；在同样的浓度和同样体积的污水注入含水层的条件下，如果污染物指数含量不超标，而其余污染物则更不会超标。模拟结果表明，污染物在地下水运移过程中，污染范围由最初的近圆形逐渐变成长轴方向与水流方向一致的近椭圆形。

污染物在地下水中随水流运移，随着时间的变化，污染羽中心向下游运移，运移范围逐步扩大，但浓度降低，5 000 d 时污染物浓度接近临界值，最大运移436 m，对下游的村庄没有影响。