近年来，随着我国城市化进程“退二进三”和“产业转移”政策落实步伐的加快，重污染行业大批关闭和搬迁，导致城市及其周边遗留大量污染地块。化工、农药厂是长三角地区一类典型的有机污染场地，具有数量多、污染程度高、污染类型复杂等特点。据文献报道，我国由工业企业搬迁而废弃遗留下来的“棕地”超过50多万块^[1]，其中包括大量的化工行业污染场地，进而导致大量化工污染物进入土壤和地下水，对当地生态和人居健康造成严重威胁^[2]。氯代有机污染物是此类场地地下水中一类代表性污染物，其易向下迁移并在含水层底部形成重质非水相液体(dense non-aqueous phase liquid，DNAPL)，例如氯苯类有机物能够在低渗透性的黏土介质中向下迁移到较深地层^[3]，导致地下水持续污染，对地下水环境产生潜在风险。因此，亟需开展我国场地地下水有机污染风险管控研究。

目前我国已经开展了一些基于保护人体健康的地下水风险评估研究^[4-7]，但是对保护水环境质量的重视不足。地下水污染具有隐蔽性、累积性和流动性特点^[8]，长期存在必然对场地周边地下水环境产生负面影响，因此应针对保护水环境质量开展地下水风险评估研究。地下水数值模型是一种预测地下水污染物溶质迁移的常见模型，国外利用对流弥散模型预测含水层污染物迁移的相关研究较多^[9]，也有研究者在此基础上考虑吸附及生物降解等自然衰减过程^[10-11]。国内对于污染物迁移模型的研究开展相对较晚^[12]，目前有研究者考虑污染物在含水层的吸附及生物降解对迁移模型进行了改进^[13-14]，但该模型对场地的特征参数和数据信息要求往往较高^[15-16]。国内外关于污染物在含水层中的迁移研究多集中于数值模型，关于解析模型的报道较少。英国水文地质风险评估(Hydrogeological Risk Assessment for Land Contamination)模型(简称P20模型)是一种求解污染物在含水层迁移过程浓度变化的解析模型，开发于1996年，发展较为成熟，特别对于小尺度范围污染场地的风险评估具有较强适用性，是定量评估地下水环境风险的重要工具^[17]。此外，美国RBCA模型(Risk-based Correction Action Tool Kit)对地下水溶质迁移风险的模拟同样借鉴了英国的P20模型^[18]。本研究以江苏省某溶剂化工场地为例，利用P20模型模拟地下水有机污染物的迁移趋势，以期为保护场地周边水环境质量提供理论依据。

1 材料与方法 1.1 研究区概况

研究场地为江苏某溶剂化工厂，占地面积约85 476 m³。如图 1所示，该厂区主要包括6个区域：制造车间、填料车间、办公楼、污水处理厂、杀虫剂生产车间、原料产品储存区。生产成品涉及的工艺原辅料包括苯、乙苯、氯苯、氯仿等，通过罐装密封储存在原料和产品储存区；而制造车间、填料车间及杀虫剂生产车间是农药生产活动的主要区域，主要生产包括甲胺磷类农药、氯碱、醚醛等化工产品。现场踏勘过程中，制造车间、生产车间及原料产品储存区内均闻到刺鼻的气味。在产品长期生产及储存过程中，原料及产品存在跑冒滴漏现象，可能造成场地地下水污染。场地的西南方向紧邻村庄，东北方向有一条河流。目前该场区周围无其他企业，主要人群分布在场地西南方的村庄内。

场地的水文地质概念模型如图 2所示，场地内地层结构自上而下分别为：回填土约1.5 m、黏土约7.5 m、粉质黏土约4.5 m、粉土约3.5 m。黏土层属于不透水层，粉质黏土层属微透水层，粉土层属弱透水层。场地内回填土与黏土层中存在孔隙潜水，但水量较少，水体不连续。

研究区内的地下水主要赋存于地表以下9 ~ 15 m的粉质黏土层和粉土层中，水位波动较小且水量相对充沛，抽水试验测得最大涌水量达到25.1 m³/d。场地东侧粉黏土层和粉土层进行抽水试验得出渗透系数为8.48×10^–4 ~ 1.00×10^–3 cm/s，导水系数2.638 ~ 3.11 m²/d。为确定地下水流向，用插值法作出等水位线如图 1所示，大致判断场区地下水整体流向为自西南向东北方向。场地周边地表水的水位未知，地下水与场区周围地表水的水力联系不明，无法确定补给关系。

1.2 样品采集与分析

场区的建筑及设备已拆除，而地下的管道及构筑物具体布局不明。地下水监测井位置如图 1所示，研究区内设置了16口地下水监测井，采样过程设置滤管的位置为地表下9 ~ 15 m，地下水取样深度为地表以下13 ~ 15 m(粉土层)，每次连续取2个平行样品。采集的地下水样品于4 ℃低温保存并于24 h内送往实验室检测，检测项目以挥发性有机物为主，分析方法参考GB/T 14848—2017^[19]《地下水质量标准》检测方法。污染物的检出浓度均为0.5 μg/L，实验室测样前进行了空白样品加标回收实验，污染物中苯系物和氯代有机物的加标回收率分别控制在85% ~ 95%，相对标准偏差低于15%。地下水中检出污染物主要包括苯、1, 2-二氯苯、乙苯、氯苯、1, 4-二氯苯、氯仿。

1.3 模拟分析方法 1.3.1 暴露概念模型

暴露概念模型表示场地污染源–暴露途径–受体之间的关系^[20]。如图 1所示场地的污染源很可能位于制造车间、污水处理厂、原料与产品储存区。如图 2概念模型所示，污染物入渗进入地下水后，可能会随着水流方向迁移到场地附近的水体中，进而对水环境产生风险。但目前场地边界地下水与场地外地表水的水力联系不明，因此基于保守原则，将场地边界地下水作为保护受体，沿地下水流方向，距离场地边界水环境最近的位置作为合规点。

1.3.2 风险评估模型

P20模型建立了多层次土壤及地下水修复目标值，适用于外部污染活动已经停止的污染土壤及地下水，即假设污染源为恒定源，且污染物浓度不随时间变化，不适用于无限污染源和污染物在源头处随时间发生变化的情景^[21]。根据污染源至受体的途径，模型分为4个层次，依据污染物的迁移途径，各层次的决定因子及对应的修复目标值不同(表 1)。土壤评估从第一层次考虑污染物的土壤–水相分配过程来推导修复目标；第二层次土壤修复目标值则考虑污染物在包气带土壤的淋溶过程，地下水修复目标值直接等于水环境质量标准；第三层次考虑污染物进入饱和带随地下水侧向迁移过程来推导修复目标；第四层次进一步考虑在受保护的水环境内稀释过程来推导修复目标。地下水风险评估因不涉及土壤，以保护水质标准为目标，从第二层次开始评估，修复目标值随评估层次的增加而逐渐升高。该模型的优势在于可以将修复集中于风险最严重的区域，以实现经济和环境的平衡。地下水评估模型假设条件为：①土壤不存在污染，即污染物已经绕过土壤区域或污染土壤已被清除；②地下水污染物在污染源和受体之间的迁移过程主要包括稀释、吸附、降解和弥散。

1) 第二次层地下水修复目标值推导。基于保护场地边界地下水修复目标值的推导从该层次开始，首先将修复目标值直接设定为地下水的某一标准浓度^[19]。然后，将污染羽中检测出的污染物最大浓度与修复目标值进行比较分析。如果浓度超过目标值，则需要采取修复措施或进行第三层次评估(包括可能的补充调查，我国标准未规定的指标可以参考欧美国家的水质标准)。

2) 第三层次地下水修复目标值推导。英国的Ogata Banks和Domenico经典解析模型简化了污染物的迁移过程，能够预测污染物的迁移变化，用于第三层次修复目标值的推导。该评估层次在第二评估层次的基础上，针对原场地地下水中风险较高的点位在场地边界建立不同的合规点，运用溶质迁移模型来预测场地边界地下水中污染物的浓度。预测污染物迁移的3种解析模型公式如式(1)~(3)所示。含水层对污染物阻滞作用可表现为有机碳对其吸附作用，因此模型中污染物的迁移速率是要小于含水层达西速率。式(4)~(6)给出了污染物的迁移速率及阻滞因子的计算模型。在综合式(1)~(9)基础上，由式(10)推导得出污染物迁移过程中的衰减因子，进而利用衰减因子来推导第三层次修复目标值。

Ogata Banks模型(非稳态)：

$ \begin{gathered} {C_{{\text{ED}}}} = \frac{{{C_0}}}{8}\exp \left\{ {\frac{x}{{2{a_x}}}\left[ {1 - \sqrt {\left( {1 + \frac{{4\lambda {a_x}}}{u}} \right)} } \right]} \right\}{\text{erfc}} \hfill \\ \left\{ {\frac{1}{{2\sqrt {{a_x}ut} }}\left[ {x - ut\sqrt {\left( {1 + \frac{{4\lambda {a_x}}}{u}} \right)} } \right]} \right\} \hfill \\ + \exp \left\{ {\frac{x}{{2{a_x}}}\left[ {1 + \sqrt {\left( {1 + \frac{{4\lambda {a_x}}}{u}} \right)} } \right]} \right\} \cdot \hfill \\ {\text{erfc}}\left\{ {\frac{1}{{2\sqrt {{a_x}ut} }} \cdot \left[ {x + ut\sqrt {\left( {1 + \frac{{4\lambda {a_x}}}{u}} \right)} } \right]} \right\} \hfill \\ \cdot \left[ {{\text{erf}}\left\{ {\frac{{z + \frac{{{S_z}}}{2}}}{{2\sqrt {{a_z}x} }}} \right\} - {\text{erf}}\left\{ {\frac{{z - \frac{{{S_z}}}{2}}}{{2\sqrt {{a_z}x} }}} \right\}} \right] \cdot \hfill \\ \left[ {{\text{erf}}\left\{ {\frac{{y + {S_y}}}{{2\sqrt {{a_y}x} }}} \right\} - {\text{erf}}\left\{ {\frac{{y - {S_y}}}{{2\sqrt {{a_y}x} }}} \right\}} \right] \hfill \\ \end{gathered} $

(1)

非稳态Domenico模型：

$ \begin{gathered} {C_{{\text{ED}}}} = \frac{{{C_0}}}{2}\exp \left\{ {\frac{x}{{2{a_x}}}\left( {1 - \sqrt {\left[ {1 + \frac{{4\lambda {a_x}}}{u}} \right]} } \right)} \right\} \cdot \hfill \\ {\text{erfc}}\left\{ {\frac{1}{{2\sqrt {{a_x}ut} }}\left( {x - ut\sqrt {\left[ {1 + \frac{{4\lambda {a_x}}}{u}} \right]} } \right)} \right\} \cdot \hfill \\ {\text{erf}}\left\{ {\frac{{{S_y}}}{{2\sqrt {{a_y}x} }}} \right\}{\text{erf}}\left\{ {\frac{{{S_z}}}{{4\sqrt {{a_z}x} }}} \right\} \hfill \\ \end{gathered} $

(2)

稳态Domenico模型：

$ \begin{gathered} {C_{{\text{ED}}}} = {C_0}\exp \left\{ {\frac{x}{{2{a_x}}}\left( {1 - \sqrt {\left[ {1 + \frac{{4\lambda {a_x}}}{u}} \right]} } \right)} \right\} \cdot \hfill \\ \;\;\;\;\;\;\;\;{\text{erf}}\left\{ {\frac{{{S_y}}}{{2\sqrt {{a_y}x} }}} \right\} \cdot {\text{erf}}\left\{ {\frac{{{S_z}}}{{4\sqrt {{a_z}x} }}} \right\} \hfill \\ \end{gathered} $

(3)

$ u = \frac{{Ki}}{{n{R_{\text{f}}}}} $

(4)

$ {R_{\text{f}}} = 1 + \frac{{{K_{\text{d}}}\rho }}{n} $

(5)

$ {K_{\text{d}}} = {K_{{\text{oc}}}} \times {f_{{\text{oc}}}} $

(6)

$ {a_x} = {\text{10% }}x $

(7)

$ {a_{\text{y}}} = {\text{1% }}x $

(8)

$ {a_{\text{z}}} = {\text{1% }}x $

(9)

衰减因子的计算：

$ {\text{AF}} = \frac{{{C_{{\text{ED}}}}}}{{{C_0}}} $

(10)

式中：C_ED为污染物迁移至场地边界合规点处的浓度，mg/L；x为污染物平行于地下水流向的迁移距离，m；y、z分别为污染源至地下水污染羽中心线的横向和垂向距离，m；S_y、S_z分别为地下水污染源宽度和厚度，m；λ为一阶衰减常数，1/d；α_x、α_y、α_z分别为地下水纵向、横向和垂向弥散度，m；u为污染物迁移速率，m/d。n为含水层有效孔隙度，无量纲；K为水力传导系数，m/d；i为水力梯度，无量纲；R_f为污染物阻滞因子，无量纲；K_d为土壤–水分配系数，L/kg；其中K_oc为含水层有机碳分配系数，L/kg；f_oc为含水层有机碳含量，以百分数计。

1.3.3 模型参数取值

基于自然衰减下的污染物在含水层中侧向迁移过程的模拟，弥散作用是描述污染物在多孔介质中受分子扩散及对流扩散的影响结果，常用弥散系数来表征^[22]。弥散系数是模型中较难确定的一个参数，可进一步转化为弥散度，目前弥散度在模型中的取值颇具争议性，关于弥散度的取值缺乏依据。在场地调查阶段缺少对弥散度的实测数据，本研究根据迁移距离结合式(7)~(9)计算得出，该模型设定横向弥散度为纵向弥散度的十分之一，此外将垂向弥散度设定为纵向弥散度的百分之一^[20]。模型主要考虑了污染物沿水流方向上的弥散，使迁移的模拟结果相对保守。

第三层次评估的场地特征参数主要来自现场实际测量和当地水文地质勘查资料，主要如表 2所示，此外检出污染物的主要理化性质如表 3所示。衰减速率应该基于对场地地下水中污染物长期的观测数据计算，本案例中无地下水长期监测数据，需要借助模型推算目标污染物的衰减速率。由于污染物中不涉及放射性衰变，因此该层次评估利用了仅考虑污染物溶解相的一阶衰减模型来计算一阶衰减速率。

1.3.4 敏感性分析

第三层次模型主要考虑污染物在含水层的自然衰减，进而模拟侧向迁移过程中的浓度分布。参数本身的不确定性会影响模型对污染物自然衰减过程模拟结果的可信度^[23-26]。第三层次涉及的参数较为复杂且模型的选择性较多，关键是分析各参数对于衰减因子的计算的影响。P20模型是利用衰减因子AF来表征污染物在含水层的自然衰减，衰减因子受到诸多参数的影响，特别是如含水层容重ρ、水力传导系数i、含水层有机碳含量f_oc等现场采集受人为影响较大的参数^[18]，场地调查工作采用加密采样来降低场地调查过程造成的参数不确定性。本研究采用局部敏感性分析的方式确定不同参数对修复目标值的影响^[27]。修正的Morris筛选法是目前常见的一种敏感性分析方法，利用灵敏度系数S衡量参数的敏感性，如式(11)所示^[28]。

$ S = \left[ {\sum\limits_{j = 0}^{m - 1} {\frac{{\left( {{Y_{j + 1}} - {Y_j}} \right)/{Y_0}}}{{\left( {{P_{j + 1}} - {P_j}} \right)}}} } \right]/m $

(11)

式中，P_j是敏感性分析参数，Y_j表示修复目标，m为模型迭代次数。本研究选取了含水层容重ρ(g/cm³)、含水层孔隙度n、含水层水力传导系数k(d/m)、含水层水力梯度i、含水层有机碳含量f_oc、含水层有机碳分配系数k_oc(L/kg)、含水层迁移距离x(m)作为敏感性分析参数，以5% 的振幅进行计算。

2 结果与讨论 2.1 地下水污染特征

地下水样中主要超标污染物为挥发性有机污染物，包括苯、乙苯、氯苯、1, 2-二氯苯、1, 4-二氯苯、氯仿。污染物的检出浓度如图 3所示，污染物浓度的最大值、最小值、平均值及检出率如表 4所示。场地水资源主要用途为工业用水，因此需要以Ⅲ类水质标准来分析检出污染物浓度是否超标。其中苯、氯苯和氯仿的检出率高于50%，最高浓度分别达到了10.37、130.36和670 mg/L，分别超出GBT14848—2017地下水质量Ⅲ类标准的1 037倍、435倍、11 166倍。

2.2 污染物迁移趋势

根据污染物检测分析结果，GW1处检出苯的浓度最高，GW6处检出乙苯、氯苯、1, 2-二氯苯、1, 4-二氯苯的浓度最高，GW14处检出氯仿的浓度最高(图 3)。因此假设GW1、GW6、GW14监测井为地下水污染源的位置，以场地边界地下水环境为保护受体，设置污染源下游沿地下水流动方向与下游场边界交汇处为本层次评估的合规点，分别记为CP1、CP2、CP3(图 1)，迁移距离分别为256、155、132 m，运用Ogata Banks、Domenico稳态和非稳态3种解析模型预测污染物迁移至合规点过程中浓度的变化。

模型预测结果如图 4 ~ 图 6所示，Ogata Banks、Domenico稳态和非稳态3种解析模型预测污染物浓度的变化趋势近乎相同，因为污染物浓度在短距离内迅速下降，导致随后迁移中污染物浓度变化不明显。同时Ogata Banks模型与Domenico稳态和非稳态模型针对氯仿的模拟结果差异较小，相比Ogata Banks模型，Domenico稳态模型由于仅考虑在污染羽流在沿水力梯度方向的弥散度，进一步简化了污染物的迁移过程。通过3种模型迁移模拟分析，污染物浓度沿水流方向在30 m内会显著下降，30 m之后的迁移过程中浓度衰减变化较为缓慢。苯、氯苯和氯仿均难以被生物降解^[29]，且氯仿浓度较高，同时有机碳分配系数相比其他污染物较低，含水层阻滞因子较小，因此在相同的水力条件下更容易迁移^[30]。如图 4 ~ 图 6，各模型的预测结果判断多数污染物在迁移6 ~ 20 m达到检出限，氯仿则迁移至75 m达到检出限，各模型预测污染物在合规点处的浓度均远低于检出限，因此场地内地下水污染物对场地边界处地下水的环境风险均为可接受，无需开展基于保护场地边界地下水的修复工作。

2.3 敏感性分析

运用Ogata Banks、Domenico稳态和非稳态模型模拟了关注污染物在含水层自然衰减下的迁移过程，利用敏感性分析比较模型参数对迁移模拟结果的影响。本研究选取7个参数进行敏感性分析，以±5% 振幅在参数均值附近取值，利用式(11)计算各参数的敏感性判别因子(S)值，并对其绝对值(|S|)进行排序，结果如表 5所示。其中对6种污染物迁移模拟结果影响最大的参数为迁移距离x，说明合规点的选取对迁移模拟影响最为显著；含水层容重ρ、含水层有机碳含量f_oc和有机碳分配系数K_oc对自然衰减迁移模拟影响也较大，表明含水层的容重及有机碳对污染物的阻滞作用是影响迁移过程的主要因素之一。如张大定等^[31]报道污染物在不同含水层中有机碳分配系数K_oc相差极大，目前可以采取利用正辛醇-水分配系数来预测其值。同时发现含水层孔隙度n、含水层水力传导系数k和含水层水力梯度i等参数对迁移模拟均有影响，这些参数在条件允许的情况下也应该进行补充调查。

3 结论

该场地地下水污染区域分布较为离散，制造车间、原料与产品储存区及污水处理厂房位置检出6种有机污染物浓度均超地下水质量Ⅲ类标准，利用Ogata Banks和Domenico模型，通过设置不同合规点位置预测地下水污染物向场地北侧边界迁移的趋势，发现污染物对场地边界水环境的风险基本可以忽略，需要综合考虑健康风险进一步确定场地是否需要修复。敏感性分析表明影响污染物迁移模拟的因素较多，后期可以在监测点和合规点之间增加长期地下水监测，用以校准衰减因子提高地下水污染侧向迁移模拟的可信度。另外有机物在含水层形成的非水相是影响侧向迁移过程的主要因素之一，今后有必要开展对于含水层中有机物两相流的侧向迁移模拟的研究，增加模型的使用范围。