鄄城县工业园区污染物下渗及迁移预测研究

毕桂真 孟令杰 韦爱菊

（山东省菏泽市水文局 菏泽 274008）

1 地下水流数值模型的建立及求解

1.1 水文地质概念模型

从含水层结构上看，研究区为多层结构区。据实测资料，园区内潜水埋深一般为1.51～4.73m，含水层厚度为14.6m，潜水层之间夹有平均厚度在5～10m的粘质砂土、粉质粘土弱透水层。此次研究把潜水作为计算目的层，通过分析钻孔资料和水文地质调查资料，将含水层概化为非均质、各向同性含水层。所以模拟区地下水的水动力条件概化为稳定的三维流。

1.2 地下水流数学模型

根据前述的水文地质概念模型，研究区地下水运动数学模型如下：

式中：H—地下水位标高（m）；H0—第一类边界水位标高（m）；t—时间（t）；A—模拟区四周边界；K—渗透系数（m/d）；Z—潜水面水头分布值（m）；x，y，z—坐标变量（m）；rw—抽水井的半径（m）；SS—含水层的释水率（L-1）；Qi—第 i口井的抽水流量（m3/d）；W—研究区三维流的源汇项（L/d）；μ—给水度；U—渗流研究区域；f（x，y，z）—初始水位分布（m）；（r，θ）—辅助柱坐标变量。

1.3 模型离散化及基础资料的给定

1.3.1 空间的离散化

研究范围是一个规则的矩形区域。将研究区垂向上剖分为4层，其中潜水含水层1层，弱透水层1层，微承压含水层1层和隔水层1层，由于研究区地处平原，地形平坦，面积较小，地形高程和各层厚度以平面的形式输入到模型中。

1.3.2 时间的离散化

选取2017年4月1日为模拟的起始时间，7300日后为模拟的终止时间。最大时间步长为200d，严格控制每次迭代的误差。在每个应力期保持含水层补给和排泄强度不变。

1.3.3 初始条件的输入

初始水位采用2017年4月统测的分层观测水位，将实测水位输入到模型中从而得到研究区数值模拟的初始流场，其中弱透水层的初始水位取上下两个含水层水位的算术平均值。

1.3.4 边界条件的输入

四周边界概化为第一类边界条件，每月1日水位值作为水位观测值输入到Visual MODFLOW中。

1.3.5 水文地质参数的初值

根据2017年3月16日稳定流抽水试验资料和水文地质勘察室内土工试验资料，结合研究区地质、水文地质条件，将以上资料所得参数输入到地下水流模型中。

1.3.6 源汇项的处理

潜水的补给来源包括降雨入渗补给、灌溉入渗补给和渠系渗漏补给，其中降水入渗补给占主导地位。地下水排泄方式包括侧向径流排泄、潜水蒸发排泄和人工开采。补给和排泄按补给强度和开采强度处理，人工开采按双井开采流量计算。降水采用鄄城县多年平均降水量609.7mm，蒸发强度采用魏楼闸水文站多年平均水面蒸发量，潜水蒸发极限埋深参考邓集试验站成果取4.0m。

1.4 数学模型的识别、验证

1.4.1 模型识别

把各种水文地质资料代入模型，以迭代残差和最大水头改变量最小为目标。通过调整分区参数值使二者尽量小，并据此来判断所用水文地质参数及分区是否合理。经反复调整参数，获得了较为满意的水文地质参数。

1.4.2 模型验证

选用2016年（年降水量611.0mm）作为验证年份，以周边11#、15B#、32A#多年地下水动态监测井作为参照井进行验证。经验证，各监测井实测水位与计算水位差值的绝对值小于0.05m，表明各观测孔的水位计算结果与实测结果吻合很好，验证了所取参数的合理性。同时对比2017年4月潜水的实测水位与模型计算水位，模拟流场与实际流场的变化趋势基本一致，可以运用到地下水水质模型中。

1.5 地下水流数学模型预测

在考虑各项补给项的条件下，使用建立的地下水流数学模型，假定园区工业取用水户双井稳定流开采地下水，单井出水量1200m3/d。连续抽水48h后实际观测记录与预测结果见表1。

2 地下水水质模型及预测

通过Visual MODFLOW中的MT3DMS模块计算污染物质的运移情况，可以求出污染物在地下水系统中的变化规律，预测研究区污染物质在不同时刻、不同的情况下所导致的地下水污染程度。

2.1 溶质运移数学模型

根据研究区的具体条件，采用下述的溶质运移模型：

式中：C—溶解于水中的污染物浓度；n—孔隙度；xi—空间坐标；qw—源（正值）或汇（负值）的单位流量；Dij—水动力弥散系数张量；t—时间；F—固相表面的溶质浓度；C0—源汇项的浓度；Ω—空间区域；Vi—地下水渗透流速；（x，y，z）—空间位置；r—研究区边界。

2.2 水质模型的建立

水质模型以水流模型为基础建立，水质模型的概化与所建立的水流概念模型相符。水质模拟区范围、含水层结构、边界类型划分、源汇项的概化均与水流概念模型相同，流体概化为不可压缩的均质流体，粘度和密度均为常数。

2.2.1 模拟因子的选择

此次研究为稳定流情况下预测研究区污染物的运移情况，考虑污染物质在含水层中的线性等温吸附（平衡），无反应项，选择化学需氧量（COD）为示踪剂。

2.2.2 边界条件、初始条件

研究区水质模型选定2017年4月1日作为初始时刻，初始时刻含水层中污染物的浓度为0。

表1 园区抽水试验观测记录与预测结果对照表

2.2.3 模型参数

溶质运移模型涉及的参数中含水介质的有效孔隙度（n）由试验所得n=0.19，其他参数的取值如下：

（1）地下水渗流速度：按照2017年3月16日实际抽水试验结果，含水层地下水垂向渗透系数速度为1.0521×10-4cm/s，由2017年3月26日室内土工试验，弱透水层地下水垂向渗透系数确定为2.73×10-6cm/s，隔水层为9.94×10-8cm/s。

（2）纵向弥散系数的确定：溶质运移模型需要的弥散系数是在结合抽水试验进行的野外弥散实验，同时取研究区的土样，采用一维土柱弥散实验法，进行室内土柱弥散实验而确定的。经过计算，评价区纵向弥散系数为1.46m/d。

（3）其他参数选用经验值。

3 预测结果

假定园区发生污染事故，污染物下渗污染当地地下水。以地面表层受到一定浓度COD污染作为初始条件。在假定的园区工业取用水户双井稳定流开采地下水的条件下，通过Visual MODFLOW中的MT3DMS模块和已建立的溶质运移模型，预测不同时间段内污染物垂直和水平运移的情况。设定污染羽外围污染物的浓度达到0.5mg/L时为该污染物的影响范围。

3.1 不同时段污染物的垂直运移情况

以地面表层受到浓度为5000mg/L的COD污染为初始条件，经过365d（1年）的垂直运移，地下水已经受到了轻度的污染。此时，在潜水的表层有仍有少许积累外，污染物已经到达潜水含水层的中下部，中心位置的浓度达到100mg/L。经过5475d（15年）的垂直运移，污染物在垂向上已经能够到达隔水层的下部，实现了对隔水层的渗透性运移。运移7300d（20年）后，垂向上的影响范围没有显著的增加，但潜水表层的污染物浓度明显增加，污染物在横向和纵向上得到了积累，地下水中的污染物运移范围进一步扩大，浓度递增显著，中心位置处不但浓度高达1000mg/L，而且范围扩大了近10倍。

3.2 不同时段污染物的水平运移情况

为对比不同浓度下污染物的水平运移情况，分别对初始浓度5000mg/L、500mg/L的COD做污染物的运移模拟，对比在7300d后的运移等值线，5000mg/L浓度COD的10mg/L等值线运移了582m，污染羽（浓度为0.5mg/L）运移了680m；500mg/L浓度COD的10mg/L等值线运移了426m，污染羽运移了606m。说明污染物浓度越高，污染羽在相同时段运移距离越远速度越快。

4 结论与建议

根据鄄城县工业园区水文地质条件的调查，运用Visual MODFLOW模型对工业园区受到污染物污染后污染物垂直和水平运移的模拟计算，得出以下结论：

（1）当园区受到污染物污染后，随着时间的推移，园区地下水从受到轻度污染逐步演变到重度污染。从初期污染物在垂向上到达隔水层的下部并渗透性运移过隔水层，到后期地下水中的污染物运移范围进一步扩大，浓度递增显著，污染物浓度逐年线性增加，若不及时治理，污染范围会进一步扩大。随着抽水时间的延长，地下水埋深逐年增加，地下水位逐年降低，污染羽的运移速度在第二含水层显著增加。

（2）园区场地基础之下第一岩土层不能满足天然防渗层的要求，园区内事故污染下渗会造成当地地下水的污染。鉴于园区地层不能满足最低渗透性标准，建议采用在园区四周和底部铺设土工膜的方法，进行人工水平防渗和垂直防渗，防止由于园区受到污染物污染对地下水水质造成影响