苏南沿江地块2，4-DCP迁移扩散的时空变化特征

尹 勇，薛 琦，张 华，苏 星，吴佳楠，张文艺

（1. 常州大学 环境与安全工程学院，江苏 常州 213164；2. 江苏龙环环境科技有限公司，江苏 常州 213031）

我国长江中下游沿江地区化工企业众多，这些企业在非正常工况下（如火灾爆炸、罐体泄漏、管线渗漏等）极易发生污染物泄漏，对周围环境造成破坏，这些污染物一旦扩散进入长江，则会对长江水体造成污染。2014年《全国土壤污染状况调查公报》显示：长江三角洲部分区域土壤污染问题较为突出［1］。近年来，随着化工园区污染预防、控制和管理措施的逐步完善，发生环境污染事故的风险大幅降低，但不可否认沿江化工园区内的污水处理站、危废库、储罐区等仍然是潜在污染源，存在污染物泄漏风险，必须予以高度重视。

地下水数值模拟是应用模拟软件建立数学模型来预测地下水及其中溶质的运移特征和未来发展趋势的一种方法，是当前预测、掌控沿江地区地下水污染问题的主要手段之一［2-3］。目前常用的地下水数值模拟软件有GMS、Visual MODFLOW、FEFLOW和TOUGH2等。其中，GMS软件具有良好的使用界面、强大的前处理和后处理功能以及优良的三维可视效果，应用最为广泛［4-6］。郑玉虎等［7］用GMS软件模拟了地下水中COD和氨氮的迁移特征，并在模拟过程中考虑了土壤吸附作用对溶质运移的影响，结果表明，研究区域土壤对COD和氨氮的吸附效果显著。高志鹏等［8］将HYDRUS-2D软件和GMS软件联用，模拟了氮素在河流-非饱和带-饱和含水层中的迁移特征，结果表明，氮素在地下水中的浓度逐渐升高，最终保持稳定。KHAYYUN等［9］利用GMS软件建立了MODFLOW和MT3DMS模型，研究了3种情景下工业废水对伊拉克底格里斯河左右两侧地下水水质的影响，结果表明底格里斯河左侧地下水水质受工业废水影响较大。RAHNAMA等［10］等调查了伊朗平原地区地下水的砷污染状况，并运用GMS软件模拟预测砷浓度的变化，发现黏土层厚度与砷浓度之间存在地理相关性。

本工作以苏南沿江某农药生产企业所在区域为研究对象，以该企业突发高浓度2，4-二氯苯酚（2，4-DCP）污水泄漏事故为情景，采用概念模型法，运用GMS软件建立地下水流场模型，以实测水位数据拟合模拟水位，通过模型识别验证获得研究区块的地下水流场模型，在此基础上建立三维地下水溶质运移模型，推演并分析地下水中2，4-DCP的迁移扩散和时空分布特征，判断其对周边环境及长江水体的影响，为沿江化工地块土壤及地下水污染管控提供技术支持。

1 研究区概况

某农药生产企业位于苏南沿江某工业园区（约10 km2），地貌单元属长江三角洲冲积平原。地面坡度小于0.5%，地形平坦，西北稍高，东南略低，平均地形高程约4.50 m，最高 5.80 m，部分地区仅2～3 m。地质构造属于扬子古陆东端的下扬子白褶带。第四系地层厚度160～200 m。

苏南沿江属北亚热带湿润性季风气候区，四季分明、雨量充沛、日照充足、无霜期较长，年平均气温约16.3 ℃，历年平均降雨量1 096 mm、最大降雨量1 888 mm、最小降雨量592 mm。

根据地质资料、现场调查和钻孔数据，确定该区域的地层范围属第四系全新统及上更新统长江下游三角洲冲积层，岩性主要为填土、淤泥质粉质黏土、粉质黏土、粉土和粉砂。水质检测结果显示该区域地下水化学类型主要为HCO3-Ca·Na型水，矿化度0.5～1.2 g/L，pH约6.2～11.4。地下水类型按潜藏条件分为孔隙潜水和承压水。其中，潜水主要埋藏于填土、淤泥质粉质黏土和粉质黏土中，其水位埋深一般为0.3～2.0 m；承压水主要埋藏于粉土和粉砂中，其埋深约10 m。苏南长江段5～10月为丰水期，11月之后进入枯水期，全年水位标高0.98～6.00 m，地下水水位会受其影响产生变化。潜水层水位变化主要受降水影响，年最大值大部分出现在7月，变化幅度0.55～1.28 m，年最小值常出现在枯水期，变化幅度0.06～0.30 m；承压水层年最高水位出现在7月，变化幅度0.61～1.64 m，年最低水位出现在1～2月，变化幅度0.03～0.40 m。

2 地下水数值模型建立

2.1 水文地质概念模型

研究区处于一个相对独立的水文地质单元，北部以长江为定水头边界；东西两侧河流与研究区地下水存在水力联系，故可概化为与河流平均状态相关的定水头边界；南部为自然的隔水边界。垂向上，研究区潜水层上部边界为水量交换边界，接受大气降雨入渗补给、河流入渗补给等，模型底部以不透水泥岩为隔水边界，具体模拟范围及边界条件见图1。根据研究区水文地质条件和地层岩性特征，污染物进入地下主要污染潜水层，从环境最不利角度出发，模拟层位设为潜水层和承压水层。

图1 研究区范围及边界条件示意图

本研究区内地下水水流运动遵循质量守恒原理，符合达西定律和渗流连续性方程。在空间维数上，以水平流为主，最终流入长江；垂向上，含水层之间存在水力流动，符合空间三维流特征。研究区水流系统的补径排特征、水位和流速随时间的变化产生差异，符合非稳定流特征。

2.2 研究区地下水数学模型

依据区域水文地质概念模型，本研究区的地下水流可概化为非均质、各向异性、空间三维结构、非稳定的地下水流系统，其对应的数学模型见式（1）。

式中：t为时间，d；Ω为模型模拟区；H为含水层水位，m；Kx、Ky、Kz分别为x、y、z方向上的渗透系数，m/d；μs为贮水率，m-1；W为含水层的源汇项，m3/d；H0为初始水位，m；Γ1为渗流区域的一类边界；Γ2为渗流区域的二类边界；n→为边界Γ2的外法线方向；K为三维空间上的渗透系数张量，m/d；q为二类边界上已知流量函数，流入为正、流出为负、隔水边界为0。

2.3 初始条件

将模拟区域内13个地下水监测孔的水位（监测初始时间为2021年1月）作为模拟预测的初始水位和模型识别验证的标准。

2.4 源汇项

区域地下水补给来源主要为垂向补给，来自大气降水入渗。研究区年均降雨量为1 096 mm，主要岩性为粉质黏土，因此，降水入渗补给系数取值0.18［11］。

区域地下水的主要排泄方式为蒸发，地下水蒸发量与其埋深有关。研究区年水面蒸发量为1 515.9 mm，地下水蒸发量比水面蒸发量小得多，研究区地下水位埋深约为2 m，主要岩性为粉质黏土，因此，蒸发系数取值0.06［11］。

2.5 水文地质参数选取

地下水模拟参数的取值见表1。参照《环境影响评价技术导则 地下水环境》（HJ 610—2016）［12］附录表B.1和B.2，潜水层水平渗透系数取0.1 m/d，给水度取0.02；由于苏南沿江为黏性土质，降雨入渗补给系数取0.000 52 m/d［12］。弱透水层与承压水层的渗透系数采用抽水试验（采用一个抽水井带两个观测井的非稳定流抽水试验方法）来确定，最终获得弱透水层、承压水层的水平渗透系数均为1.529 8 m/d。垂直渗透系数与水平渗透系数的比值按经验值取0.1。

表1 地下水模拟参数的取值

2.6 流场模型的建立与识别校正

在GMS软件中将模型在XY方向上以80×80进行剖分，垂直方向上分为3层，每个单元格长宽均为50 m，共计12 800个单元格，有效活动单元格为5 728个。将表1的参数输入GMS软件，建立研究区的地下水流场模型。

模型的识别和验证在地下水数值模拟过程中十分重要，是对所建立模型的一个校核过程，目的在于确保对该研究区的地下水数值模拟是正确的，所建立的模型是精确可靠的。主要原理是通过不断修正降雨入渗补给系数、渗透系数等水文地质参数，重复计算，直到模拟结果与水位观测点的实测值之间的误差在可接受范围内。GMS软件提供了一套辅助校核工具，可以帮助用户验证实测水位与模拟水位的吻合性。

在建立地下水流场模型后，根据所监测的水位资料对模型进行识别校正，最终得到校正后的地下水流场（见图2）。由图2可见：水流整体往东北方向流动，汇入长江，沿水流方向水位逐渐递减，但减幅不大。由于研究区的中部小块区域缺少部分实测点的水位数据，导致模拟值与实测值具有一定差别，但总体上看，校准后的地下水流模型能够较好地反映该区域地下水流场的特征。校核后的水文地质参数见表2。

表2 校核后的水文地质参数

图2 校正后的地下水流场

选择研究区域内13个监测井的实测水位与模拟水位进行拟合，设置实测值与计算值的极差为0.5 m，置信区间为95%，结果见图3。由图3可见：算点都均匀分布在45°线两侧，R2为0.973 6，表明拟合程度较高，误差较小。

图3 水位实测值与计算值的拟合曲线

图4为水位拟合误差情况图，根据误差棒颜色的不同，可以判断误差的大小：绿色表示实测值与计算值的误差在置信区间范围内；橙色表示误差超出置信区间，但不超过置信区间的200%；红色表示误差超出置信区间200%。由图4可见，只有一个点的误差棒指示为橙色，其余均为绿色，没有出现红色。模拟结果中水位绝对误差小于0.5 m的观测点占已知水位观测点的92%，远高于《地下水资源管理模型工作要求》（GB/T 14497—1993）［13］的规定：预测水位值绝对误差在0.5 m内的观测点须占所有已知水位观测点的70%以上。这表明所建立的模型精度符合要求。

图4 水位拟合误差情况图

综上，从地下水流场和观测孔水位对比两方面看，水位观测孔拟合结果良好，地下水流场宏观拟合效果良好，模型识别的水文地质参数与水文地质条件基本符合，能够反映地下水流动特征，具有较高的精度。因此，可基于该地下水流模型模拟2，4-DCP的迁移转化特征。

3 2，4-DCP的迁移转化

3.1 2，4-DCP运移模型的控制方程

模拟2，4-DCP的迁移时只考虑对流、弥散和吸附作用，描述三维地下水流系统中某种污染物迁移运动的控制方程见式（2）。

式中：R为滞后因子；θ为地下介质的孔隙度；ρ为溶质的质量浓度，mg/m3；xi，xj分别为沿X轴和Y轴的距离，m；Dij为水动力弥散系数张量，m2/d；vi为渗流或线性孔隙水流速度，m/d，与单位流量或达西流量qi有关，vi=qi/θ；qs为单位体积含水层源（正值）和汇（负值）的体积流量，1/d；ρs为源汇流中溶质的质量浓度，mg/m3；λs为吸附相的反应速率常数，1/d；λ为溶解相的反应速率常数，1/d；ρb为多孔介质密度，mg/m3。

3.2 污染情景及参数设置

本次模拟情景设置为该农药厂地块内污水收集管破损，导致工业废水通过土壤包气带进入地下水造成污染。假设裂口面积约为0.01 m2，污水渗漏量约为0.7 m3/d，污染类型为连续入渗。经调查取样分析，该工业废水中污染物以苯系物为主，其中2，4-DCP质量浓度较高，为107 mg/L。因此，选择2，4-DCP作为典型污染因子进行污染物迁移规律研究。

模拟期为7 300 d（20 a），模拟输出时间段为100 d，1 000 d，1 825 d（5 a），3 650 d（10 a），5 475 d（15 a），7 300 d（20 a）。本次模拟不考虑化学反应、生物降解以及外部源汇项（溶质通过源进入或者通过汇离开研究区域），即忽略扩散过程中2，4-DCP的生物降解和化学分解等因素。

模拟过程中相关参数的设置见表3。其中：2，4-DCP的干预值参照荷兰干预值标准［14］；吸附常数经实验室等温吸附实验测得；纵向弥散度通过实验室土柱实验和现场试验获得，横向弥散度通常比纵向弥散度低一个数量级，垂向弥散度比横向弥散度低一个数量级［15-16］，横纵弥散度比值取0.2，垂纵弥散度比值取0.05；孔隙度采用沿江地区的经验值［17］。

表3 模拟过程中相关参数的设置

3.3 2，4-DCP迁移扩散过程中的时空变化特征

初始化三维地下水溶质运移模型，将相关参数导入已建立的地下水流场模型，运行后得到2，4-DCP在不同模拟时间条件下的时空变化特征，结果见图5和表4。如图5可见：2，4-DCP沿地下水水流朝东北方向呈椭圆形扩散，污染面积逐渐扩大（图5a）；当模拟时间为5 475 d（15 a）时，2，4-DCP在水平方向上已经扩散至长江江岸（图5e）；当模拟时间为7 300 d（20 a）时，2，4-DCP在水平方向上已经扩散进入长江（图5f）。

表4 模拟过程中2，4-DCP的迁移情况

图5 不同模拟时间条件下2，4-DCP在第一层含水层中的时空分布特征（埋深0.2～1.5 m）

由表4可见：在0～100 d内，2，4-DCP的水平迁移距离为121.15 m，扩散速率最快，平均扩散速率为1.211 5 m/d；100～1 000 d时，2，4-DCP的扩散速率有所降低，水平迁移距离为64.40 m，平均扩散速率为0.071 6 m/d；1 000～1 825 d时，水平迁移距离为53.19 m，平均扩散速率为0.064 5 m/d；1 825～3 650 d时，水平迁移距离为94.41 m，平均扩散速率为0.051 7 m/d；3 650～5 475 d时，水平迁移距离为76.17 m，平均扩散速率为0.049 3 m/d；5 475～7 300 d时，2，4-DCP的扩散速率最低，水平迁移距离为57.35 m，平均扩散速率降低至0.037 1 m/d。综上，随着模拟时间的延长，2，4-DCP的迁移扩散速率在不断降低。如果不采取任何污染控制措施，到7 300 d（20 a）后，2，4-DCP在水平方向上的最大迁移距离为466.67 m；污染面积将达到110 982 m2，而且，随着时间的推移，污染面积还会持续增大。由表4还可看出：在污染前期（0～3 650 d），随着模拟时间的延长，污染羽中心2，4-DCP的质量浓度逐渐增大；当模拟时间为7 300 d时，污染羽中心的2，4-DCP质量浓度达到最大，为41.453 mg/L，超过荷兰干预值标准（0.03 mg/L）1 375倍。而且，随着时间的推移，污染羽中心2，4-DCP的浓度还会持续增大。

图6为模拟时间7 300 d时垂直方向上地下水中2，4-DCP的分布情况。由图6可见：在垂直方向上，2，4-DCP由潜水层先后穿透弱透水层和承压水层，扩散至模型底部，最大污染深度可达40 m；随着地下水深度增加，污染物扩散速度明显减慢，由初始的-5 m处的41.453 mg/L扩散至-25 m处的4.607 mg/L大约经历7 300 d（约20 a），表明2，4-DCP的垂直扩散速率相对较低，这与江水丰水期的顶托作用有关。

图6 模拟时间为7 300 d（20 a）时垂直方向上地下水中2，4-DCP的分布

以上分析表明：土壤吸附对2，4-DCP的滞留作用有限，2，4-DCP很容易穿过粉质黏土层，对地下水造成污染；2，4-DCP的净化主要依靠含水层的稀释作用，这意味着该区域地下水一旦受到污染，只依靠地下水的循环是难以实现自净的。由此可见，该厂工业废水的泄漏对周边地区的土壤地下水环境安全以及长江生态存在一定威胁，应及时采取有效措施防止2，4-DCP的扩散。

4 结论

由于沿江地块场地土壤及水文地质条件的复杂性和不确定性，很难建立起科学、准确反映污染物弥散、吸附和降解等环境行为的数学模型。基于此，本工作采用一种间接模拟方式，即先利用GMS软件构建能反映研究区域的地下水流场模型，然后利用地下水观测井实测数据检验、校对地下水流场模型，从而建立起能准确反映本研究区域地下水流场的模型。在此基础上，利用建立的流场模型构建溶质运移模型，揭示污染物扩散规律，即推演地下水中2，4-DCP在迁移扩散过程中的时空分布特征，判断其是否会对周边环境及江水造成污染影响。

a）模拟结果表明：2，4-DCP将会随着地下水流往长江方向迁移，100 d后，2，4-DCP的水平迁移距离为121.15 m，1 000 d后为185.55 m，5 a后为238.74 m，10 a后为333.15 m，15 a后为409.32 m（届时2，4-DCP将扩散至长江江岸），20 a后为466.67 m；20 a后污染羽中心2，4-DCP的质量浓度将达到41.453 mg/L，最大污染深度为40 m，污染面积为110 982 m2，2，4-DCP将扩散至长江水体。可见，该厂工业废水的泄漏对周边地区的土壤地下水环境安全以及长江生态存在一定威胁，应及时采取有效措施防止2，4-DCP的扩散。

b）本文所建立的沿江地区污染物扩散与预测方法，对于沿江化工地块土壤及地下水污染管控有一定的参考价值。