刘景兰，李立伟，朱明奕，葛菲媛，秦 磊

(1.天津市地质研究和海洋地质中心，天津 300170； 2.天津市生态环境科学研究院，天津 300191)

自1984年以来，我国石化园区经历了从有到无，再发展壮大的过程。1990年全国仅有省级以上园区13个，2014年则发展至近700家，且其中半数以上分布于我国沿海地区。随着我国石油化工行业的迅速发展，石化园区对周边水土环境构成的污染风险也成为了管理部门和学界关注的焦点。石化园区中的生产活动，往往会造成较为严重的水土污染问题。近些年，随着我国石化园区环境管理制度的不断完善，园区生产所产生的废水、废气和废渣逐步实现了统一管理，导致污水处理厂和固废暂存场地等区域成为了石化园区的潜在污染源。园区污水处理设施在非正常状况下可能发生渗漏事故，在设施检修的间隔期，污染物渗入地下水对周边环境所构成的风险尤其值得关注，针对非正常状况下污染物在地下水中迁移过程的预测对园区环境管理和园区周边敏感目标的保护具有重要的现实意义。

随着计算机模拟技术的发展，数值模拟方法被广泛用于地下水污染物迁移的模拟预测。目前国内外最为常用的数值模拟软件包括GMS(Groundwater Modeling System)、Visual MODFLOW和FEFLOW等，其中GMS和Visual MODFLOW是集成MODFLOW、MODPATH和MT3DMS等模拟程序开发的商业软件，这些模拟程序都是基于有限差分方法的，而FEFLOW软件则是基于有限元方法开发的，这些模拟软件为地下水污染物迁移转化的预测研究提供了有效工具。翟勇等[1]利用GMS软件对我国东北某化工场地地下水中挥发酚类物质在未来10年中的迁移和污染物浓度衰减过程进行了模拟预测。何华燕等[2]利用FEFLOW软件对氯代烃类污染物在地下水中的稀释扩散过程进行了模拟，发现卤代烃类污染物迁移能力较弱，扩散范围较小，对周边环境影响有限。魏传云等[3]利用Visual MODFLOW软件的MT3DMS模块，对某化工园区污染物在不同情景下的运移规律进行了模拟预测，并重点针对连续源强和瞬态源强条件下污染物迁移的差异性进行了对比分析。姜凤成等[4]采用数值模拟方法对某化工场地石油类污染物和重金属Zn在地下水中的时空演变规律进行了对比分析，发现石油类污染物泄漏后在含水层中可发生自然衰减，污染物浓度和范围随时间逐渐减小，而重金属Zn会随地下水不断向下游迁移。王平等[5]针对某有机污染场地地下水中污染物对周边拟建水库的潜在影响进行了数值模拟分析，发现水库建设会改变周边地下水流场，从而使得地下水污染羽范围有所扩大，但并不会对水库水质造成影响。张凤娟等[6]对淄博市张店化工园区建设对地下水环境的影响进行了模拟预测，发现在20年中石油类污染物在地下水中迁移了350 m的距离，影响面积达到3万m2，对地下水环境影响严重，有必要采取有效的防渗和监测措施加以防范。邓一荣等[7]利用GMS软件的MT3D模块对粤港澳大湾区某典型化工污染场地苯系污染物在地下水中的迁移过程进行了模拟分析，发现溶解态苯系污染物沿地下水流向往下游迁移，通过识别检验发现运移过程主要受对流弥散作用的控制。

本文以我国北方沿海某石化园区为研究区，以园区内污水处理厂作为潜在污染源，以石油烃为特征污染物，利用GMS软件的MODFLOW和MT3DMS模块建立石油烃三维迁移数值模型，对非正常状况下污染物渗入地下水后在未来30年中的迁移过程进行模拟预测，研究可为我国沿海地区石化园区环境管理提供科学依据。

1 材料与方法

1.1 研究区概况

研究区范围如图1所示，该石化产业园区建立于2003年7月，产业结构以精细化工和医药制造为主，占地面积约746.84 hm2。研究区地貌类型属海积冲积低平原，地势平坦，地面坡度小于1/10 000。该石化园区被河流围绕，北侧为城排河，南侧依次为荒地排河和独流减河，西侧城排明渠自北向南沟通城排河和独流减河，东侧是长青河。

图1 研究区范围和水系分布示意Fig.1 Schematic of range and water system distribution of the study area

1.2 研究区水文地质条件

研究区场地位于我国北方海积低平原的咸水及盐卤水分布区。根据区域水文地质调查成果[8-10]，研究区400 m以浅的平原松散地层孔隙水可划分为4个含水岩组：第Ⅰ含水组属于浅层地下水系统；第Ⅱ—Ⅳ含水组属深层地下水系统。第Ⅰ含水组为咸水及盐卤水，底板埋深约100 m，其中顶部潜水含水层隔水底板埋深在18 m左右；第Ⅱ含水组在调查评价区内为咸水，咸水底界埋深在200 m以深；深部的Ⅲ、Ⅳ含水组均为承压淡水。第Ⅰ含水组属滨海平原冲海积层咸水及盐卤水区，该部分地下水无开发利用情况，不作为居民生活饮用水使用[11-12]。因此，研究将地下水环境保护目标设置为第Ⅰ含水岩组的顶部潜水含水层，以防止污染物进入该潜水含水层向下游迁移扩散对园区外部地下水和地表水体造成污染，并对下游敏感目标构成环境风险。

第Ⅰ含水岩组顶部含水层地下水水位埋深在0.24～1.29 m，平均水位埋深为0.70 m。包气带主要为人工填土层素填土、粉质黏土；含水层为第四系松散岩类，岩性以粉质黏土、淤泥质黏土为主。潜水含水层厚度为15.01～18.71 m，淤泥质黏土单层厚度2.90～9.60 m，粉质黏土单层厚度0.90～12.70 m。含水层渗透系数为0.05～0.20 m/d，单位涌水量为0.015～0.270 m3/(h·m)，富水性弱。地下水化学类型为Cl-Na型水，pH值为7.45～7.66，矿化度为15～54 g/L。

1.3 概念模型建立

1.3.1 含水层结构

模拟区内主要分布第四系松散岩类孔隙水，含水岩组由上至下分别为人工堆积层、全新统上组陆相冲积层和全新统中组浅海相沉积层，岩性由上至下依次为素填土、粉质黏土、粉土和淤泥质黏土。该潜水含水岩组富水性较差，其下伏粉质黏土厚度较大，分布连续，且水平渗透系数和垂直渗透系数在10-7～10-8cm/s，构成了地下水潜水良好的隔水底板。因此，本次将第四系松散岩类孔隙水概化为统一的潜水含水层，同时考虑垂向上的岩性变化，构成非均质各向异性的多层结构地下水潜水三维模型。

1.3.2 边界条件

根据模拟区水文地质条件及地下水流场特征，模拟区内地下水由区域中心分别向北西、南西和南东方向的河流进行径流排泄，几乎无地下水的侧向流入补给，因此将模型侧向边界均概化为定水头边界。模型的上边界为潜水面，整个含水层系统通过该边界接收大气降水补给和蒸发排泄，从而与外界进行垂向上的水量交换。

1.3.3 水文地质参数

由于区内潜水富水性较差，且水质普遍咸化，园区内企业几乎不对区内地下水潜水进行开发利用。因此，模拟过程中无需考虑人工开采地下水所导致的地下水潜水位的动态变化。基于上述原因，此次所构建地下水流数值模型为稳定流模型。

地下水稳定流模型所需率定的水文地质参数主要为渗透系数K，根据所收集的区内抽水试验计算结果，结合场地含水层结构特征，确定潜水含水层水平渗透系数Kh为0.05～0.20 m/d，将垂直渗透系数Kv设置比水平渗透系数小一个数量级。根据经验值，有效孔隙度设置为0.3。由于水动力弥散尺度效应的存在，难以通过野外或室内弥散试验获得真实的弥散度。因此，模拟参照经验值，按照保守的评价原则，取纵向弥散度值为15 m，横向弥散度值为1.5 m，垂向弥散度为0.15 m。

1.3.4 源汇项

模拟区内地下水主要接收降水补给，侧向流入补给几乎可以忽略。根据包气带岩性和地下水埋深特征，结合经验值，确定大气降水入渗系数为0.11。由于地下水潜水水位埋深较浅，一般为0.5～2.0 m，蒸腾作用较为强烈，区内地下水主要通过蒸发排泄，少部分地下水通过侧向流出进行排泄。此外，根据地下水潜水监测流场形貌判断，研究区内存在局部小规模地下水降落漏斗，推测该处可能存在少量地下水人工开采。

1.3.5 源强概化

根据该化工产业园区污水处理厂设计进水水质指标限值以及各污染因子进水浓度和标准限值所计算的标准指数，发现石油类标准指数最高，对地下水所构成的潜在威胁最大。因此，将研究中污水处理设施渗漏情境下对地下水可能造成影响的预测因子确定为石油类。

根据项目工程设计，污水处理各构筑物的池体均为钢筋混凝土。池体的泄漏量参考《给水排水构筑物工程施工及验收规范》(GB 50141—2008)中关于满水试验验收的要求，钢筋混凝土池体满水试验验收标准为2.0 L/(m2·d)，假设项目在非正常状况下池底由于地面不均匀沉降或地下水对池体的腐蚀等多种因素影响下，出现防渗层破裂情况，破裂程度引起的地下水渗漏量按照验收标准的10倍计算[13-17]，即20 L/(m2·d)。

沉淀池池体长9.7 m、宽8.5 m、高5.0 m，本次采用池体的表面积作为泄漏面积，即264.45 m2，结合单位面积渗漏量20 L/(m2·d)，可计算得到非正常状况下沉淀池池体渗漏量为5.289 m3/d。由工程分析可知石油类污染物的进水浓度为20 mg/L，假定污水处理设施的地下水防控或检漏周期为120 d，即发生非正常状况后120 d发现并进行修复切断渗漏源，因此进入含水层中石油类污染物渗漏量为12 213.6 g。

1.4 数学模型和软件选取

1.4.1 地下水流模型

针对本次所建立的非均质各向异性三维稳定流地下水流模型，选取下列数学模型和定解条件进行求解：

(1)

式中，h为地下水潜水位；Ω为地下水渗流区域；Γ1和Γ2分别为模型第一类边界和第二类边界；Kxx、Kyy和Kzz分别为x、y和z主方向的渗透系数；λ为通过单位水平面积潜水面注入含水层的水量；ω为源汇项，包括降水入渗补给和井抽水量等；h0(x，y，z)为第一类边界的地下水水头函数；q(x，y，z)为第二类边界单位面积流量函数。

1.4.2 溶质运移模型

地下水水质模型是定量评价和研究地下水水质问题的数学模型。建立地下水水质模型时，将溶质运移过程中的对流、弥散和吸附等反应项考虑在内，用于描述地下水中溶质的运移状态。模型的偏微分方程式及定解条件如下：

(2)

式中，方程左侧前2项为弥散项，之后2项为对流项，最后一项为由化学反应或吸附解吸所产生的溶质增量；Dxx和Dyy分别为x和y两个主方向的弥散系数；μx和μy为x和y方向上的实际水流速度；c和c0分别为溶质浓度和初始浓度；Ω为溶质运移的区域；Γ为三类边界；φ为边界溶质通量；v为渗流速度；gradc为浓度梯度。

1.4.3 软件选取

研究所选用的GMS软件综合了MODFLOW、MODPATH和MT3DMS等已有的地下水模型基础上研发的。GMS软件相较于Visual MODFLOW等地下水模拟软件最大的特点是提出了概念模型建模的方法，这为GMS与GIS软件间提供了良好的数据接口，使得建模过程更为灵活简便。除此之外，GMS因其强大的前处理和后处理功能以及优良的三维可视化功能，目前已经成为了国际上最受欢迎的地下水模拟软件之一。

1.5 网格剖分和应力期划分

1.5.1 网格剖分

空间网格剖分是建立地下水水流和溶质运移数值模型的基础。考虑到需要利用所建立的模型对重点区域污染物迁移过程进行模拟预测，因而在相应区域进行了网格的加密剖分。东西向剖分网格64个，南北向剖分网格74个，垂向上剖分网格5层，共剖分网格23 680个，其中有效计算单元15 960个。具体空间网格剖分如图2所示。

图2 数值模型空间网格剖分示意Fig.2 Schematic of the spatial grid profiles of the numerical model

1.5.2 应力期划分

此次模拟将地下水非稳定流模型和溶质运移模型模拟期划分为2个应力期：①第1个应力期设置为120 d，对应污染物入渗阶段；②第2个应力期设置为11 000 d(30年)，对应污染物扩散阶段，最长时间为污水处理厂最长服务年限(30年)。

2 结果与讨论

2.1 地下水流模型识别检验

2.1.1 地下水位拟合

通过模型参数识别，得出模型计算地下水流场和实际监测地下水流场拟合情况如图3所示。

图3 计算水位与实际监测水位拟合Fig.3 Fitting of calculated water level and actual monitoring water level

从图3中可以看到，模型计算地下水流场和实际监测地下水流场拟合情况较好，二者之间的水位误差整体上小于±0.1 m，说明所建立的地下水流模型能够较为准确地反映评价区地下水潜水的真实流动状态，所建立的模型可用于评价区地下水中污染物迁移过程的模拟预测。

2.1.2 地下水均衡分析

通过模型参数识别，得出模型验证期地下水水量均衡结果见表1。从表1中可以看出，模拟区地下水潜水主要接受大气降水入渗补给，并主要通过潜水蒸腾方式进行排泄，并有少量地下水通过侧向流出进行排泄。

表1 模拟区地下水均衡分析Tab.1 Groundwater balance analysis in simulated area

2.2 污水处理设施渗漏模拟预测

2.2.1 地下水流场动态模拟

渗漏事故情景下，污水渗入地下会构成新的地下水潜水补给来源，除此之外，模拟区地下水系统的其他源汇项基本保持不变。因此，非正常状况下污染物在地下水中的迁移预测，有必要在已有地下水流模型的源汇项条件和含水层特征的基础上，增加入渗补给地下水潜水的污水水量，重新进行地下水潜水非稳定流模型计算。污水入渗时间设置为120 d，计算得到的污水处理厂周边地下水流场变化情况如图4所示。

图4 非正常状况下污水处理厂周边地下水流场变化情况Fig.4 Changes of groundwater flow field around sewage treatment plant under abnormal conditions

从图4中可以看到，随着渗漏时间延长，在污水处理设施周边开始形成地下水水丘，且规模不断增大，在120 d时达到最大，且在模拟1 000 d时也未能完全恢复。地下水水丘的形成使得泄漏点周边地下水水力梯度增大，这会加快污染物向周边的扩散速率。

2.2.2 地下水溶质运移模拟

污水处理设施泄漏对地下水污染预测结果如图5所示。

图5 污水处理设施泄漏对地下水污染预测结果Fig.5 Prediction results of groundwater pollution caused by leakage of sewage treatment facilities

从图5中可以看出，在模拟期内石油类对潜水含水层将会造成污染，并存在超标现象。石油类在潜水含水层中不同时间的超标范围及最大运移距离见表2。

表2 非正常状况下污水处理设施泄漏对含水层影响范围Tab.2 List of influence range of aquifer caused by leakage of sewage treatment facilities under abnormal conditions

从表2中可以看到，1 000 d污染物最远迁移41 m左右，超标范围为4 879 m2左右；30年(11 000 d)污染物最远迁移56 m左右，超标范围为6 482 m2左右。根据模拟结果，在污染物发生泄漏30年后，地下水中污染物运移距离有限，污染羽前缘距离下游地表水体还有120 m左右的距离，这主要与研究区地下水水力梯度较小、同时含水层渗透系数较低有关。模拟结果表明，污水处理设施泄漏对周边地表水体影响较小，但由于周边土壤和地下水遭受污染，对泄漏点周边敏感目标可能构成环境风险，有必要在污水处理设施下游进行地下水水质定期监测。

3 结论

(1)利用GMS构建了研究区三维非稳定地下水流和地下水溶质运移数值模型，通过对污水处理设施非正常状况下废水的泄漏进行模拟预测可知：池体泄漏条件下，泄漏点周边会形成地下水水丘，在120 d时规模达到最大，加速了污染物向周边的迁移速率；研究区属低渗透性地层，预测期内污染物运移距离较小，尚未对下游地表水体产生影响。

(2)为防控地下水污染，避免对泄漏点周边地下水环境造成影响，结合预测结果，建议在污水处理设施下游30～50 m处布设跟踪监测井，定期开展地下水监测工作，以便及时发现污染并截断污染源。

(3)此次模拟遵从保守性原则，未考虑污染物的吸附和降解等生物地球化学过程，后期可根据管理需要，通过补充相关反应动力学和热力学参数，进一步提高模型模拟预测精度。