1、 采用敏感性分析和最差情况分析进行地下水污染风险评价 比利时 Marijke Huysmans 李 烨 译；冯翠娥、魏国强 校译 本文阐述了如何利用敏感性分析和最差情况分析来对地下水污染进行风险评价，并将这一方法应用于匈牙利境内的一个研究区。在该研究区存在几个地下水污染源，而且位于饮用水井附近。主要关注的是污染源是否对饮用水井造成威胁。由于资料有限，模拟的结果具有很大的不确定性。应用敏感性分析来评价这一不确定性。 一、概 述 确定与地下水污染有关的环境风险是一个常见的研究问题。通常需要研究受到污染的 地下水是否会到达饮用水井、河流、生态脆弱区或动植物敏感区（ Calow， 1998）。计算机模
2、型是常用的进行地下水流动和污染预测的工具。由于资料的缺乏和模型参数的异质性，往往会造成模拟结果的不确定性。尽管与预测相关的不确定性可能非常关键，但是通常都会被忽略，特别是在资料相当缺乏的情况下（ Levy 等， 1998）。 有几种解决参数不确定性的方法。常用的随机分析方法是蒙特卡罗（ Monte Carlo）模拟（ Asante-Duah， 1998）。应用这一方法，可以在所能输入的概率分布范围内自由选择输入值，并计算出每种情况 下的输出值。通过重复计算可以确定输出的分布范围。尽管蒙特卡罗模拟方法的功能强大，而且渐近收敛，但是计算效率却很差。另外，还需要确定每一个模拟参数的概率函数，这是在输
3、入信息比较缺乏的情况下，该方法的一个严重的缺陷。蒙特卡罗模拟通常与地质统计学方法相结合，然而，地质统计学方法需要通过大量的数据来准确描述每一个参数的空间变化，而在实际工作中，通常不可能获得如此充足的资料。解决不确定性问题最直接的方法是敏感性分析和（或）最差情况分析。通过敏感性分析，可以确定由于输入变量和参数造成的输出变化，这是一种检验模型中输出变 量对输入变量的灵敏性的分析方法。在最差情况分析中，需要给定每一个变量和参数最差的可能值，这样可以得出模型的最差输出结果。 本文根据敏感性分析和最差情况分析进行风险评价，研究区 Mtzalka 市拥有 25400人口，位于匈牙利东部，靠近罗马尼亚和乌克
4、兰边界。 Mtzalka 沿 Kraszn 河分布， Kraszna河与 Tisza 河相汇，最终排泄到多瑙河。 Mtzalka 周围有几个潜在的污染源。其中一个地下水污染源是城市垃圾处理场，占地面积 80 万 m2，没有加适当的防渗层，在高水位期地下水水位会到达垃圾场的底部边界 。另一个地下水污染源是在 1971 1997 年间利用的污水氧化池（ Nauner， 2000）。现在，该池被土壤和植被覆盖，但是在地下土壤中仍存留大量的污水和污泥。第三个地下水污染源是污水处理厂，一级处理是通过滤网将木头、纸和塑料等去除，并通过大的沉降池将固体和水分离开来，大部分固体沉到池的底部，在这一阶段， 70
5、%左右的固体是污泥，约 1 万 m3 左右 。与污染处理系统没有联系的住宅区是第四个地下水污染源，这些住宅区的污染坑没有用混凝土加固，因此，污水容易到达地下，特别是在地下水位较高时。工业活动是第五个地下水污染源。主 要关注的问题是这些污染源是否会对饮用水井造成威胁，要研究这一问题相当复杂，但是资料却很有限。 二、地质条件 Mtzalka 位于匈牙利平原，是帕诺尼亚（ Pannonian）山间盆地的一部分。研究区内更新统沉积物厚度约为 260m，下更新统厚度为 110m，中更新统厚度为 90m，上更新统厚度为 60m。下更新统的沉积物主要由砂砾组成，是该区的含水层系统，是重要的饮用水源。中更新统
6、是区域隔水层，由渗透性差的冲积相或湖相粉砂和粘土组成。上更新统由中细砂和粉砂组成，也可以作为含水层，但渗透性比上更新统要差。下部的粘土层是隔水层，可以作为地下水流动和运移模型的隔水底板。 根据 23 口水井的钻孔资料可以对不同地层单元进行评价，通过将更新统划分为 6个水文地层单元可以将复杂的地质条件进行简化（表 1）。第一层是各向异性的含水层，渗透性较强，由许多较薄的砂层、粉砂层和粘土层组成；第一层和第二层是主要由粘土层和砂质粘土层组成的半透水层，渗透性较差；第三层是连续的砂层，有几口水井将该层作为水源；第五层和第六层是最好的含水层或含水单元，水力传导系数最大，所有的饮用水都是从这两层获取。
7、表 1 更新统的 6 个水文地质单元描述 水文地质单元 平均厚 度（ m） 描述 地层 第一层 65 薄砂、粉砂和粘土层 上更新统 第二层 25 粘土或粘质砂层 中更新统 第三层 7 砂层 第四层 40 粘土或粘质砂层 第五层 20 粗砂砾层和粘土互层的各向异性单元 下更新统 第六层 100 厚砂砾和粘土互层 三、地下水流动模型 根据 MODFLOW 来求解微分方程。 （一）边界条件 水文地质模型是一个 9km 10km 260 km 的模型。上部的更新统粘土沉积物表示模型的不透水边界；根据测压图和剖面图获得测压水头条件，将测压水头条件作为含水层一、三、五和六 的边界；第二层和第四层的垂直边界
8、是零通量边界，这两层主要由粘土组成，地下水流动相当慢，而且水流主要是在垂直方向上运动，因此，假定沿垂直边界的水平通量为 0。模型的东部边界是一条河流，这条河流的坡度为 13cm/km，河流底部沉积物厚度约为 0.7m，水力传导系数约为 10 6m/s。按照每月的抽水资料，确定模型中从井场抽取的地下水量。共有 15 口抽水井，总抽水量为 15000m3/天。根据 Thorntwaite方法，估计含水层的补给量为 30mm/年。 （二）栅格 设置 6 层 104 行和 112 列的栅格，基本单元是 100100m2，在抽水井附近栅格大小一般是 5050m2，每个单元一般者不得超过其邻近单元的 1.
9、5 倍。为了进行计算，每个单元的长宽比不超过 10。 （三）水力传导系数 利用抽水试验、排泄和水位下降资料以及粒度分布资料，确定研究区不同地层的水力传导系数。在第一层进行了 12 次试验，在第五层进行了 3 次试验，在第六层进行了 8次试验。采用 Thiem-Dupuit 等式来分析排泄和水位下降资料，在第一层进行了 13 次分析，在第五层进行了 1 次，在第六层进行了 22 次。采用 Beyer 公式和 Zamarin 公式对第六层的 6 个土样进行了粒度分析 ，这两个公式是分析粒度与传导系数之间关系的经验公式。在第二层和第四层没有分析水力传导系数，取以前的研究数据。 计算每一层水力传导系数
10、的平均值。在水平层的沉积物中，水平水力传导系数要高于垂直水力传导系数，假定 Kh与 Kv 的比等于 10（表 2）。第五层和第六层是渗透性最强的更新统地层；第二层和第四层是渗透性最差的更新统地层，它们形成了地下水流动的天然屏障。 表 2 每一层的平均水力传导系数测量值 层 水平水力传导系数 Kh(m/s) 垂直水力传导系数 Kv(m/s) 一 5.810 5 5.810 6 二 1.310 7 1.310 8 三 1.310 5 1.310 6 四 1.310 7 1.310 8 五 7.010 4 7.010 5 六 3.710 4 3.710 5 （四）校 正 通过 32 个测压计测量地下
11、水水位来进行校正，其中 17 个测压计位于上更新统（第一层），另外还有 15 个位于下更新统（第六层）。根据稳定态条件对模型进行校正，首先利用 “试错 ”法对水力传导系数和补给进行校正，之后采用 PEST 进行自动校正。 根据 “试错法 ”获得的水力传导系数见表 3。 表 3 采用“试错法” 校正前后的水力传导系数 初始值 校正值 Kh1 5.810 5 m/s 区 1： 1.510 4 m/s 区 2： 2.010 5 m/s 区 3： 8.010 5 m/s 区 4： 1.510 5 m/s Kv1 5.810 6 m/s 区 1： 1.510 5 m/s 区 2： 2.010 6 m/s
12、 区 3： 8.010 6 m/s 区 4： 1.510 6 m/s Kh2 1.310 7 m/s 9.010 8 m/s Kv2 1.310 8 m/s 9.010 9 m/s Kh3 1.310 5 m/s 1.310 5 m/s Kv3 1.310 6 m/s 1.310 6 m/s Kh4 1.310 7 m/s 9.010 8 m/s Kv4 1.310 8 m/s 9.010 9 m/s Kh5 7.010 4 m/s 7.010 4 m/s Kv5 7.010 5 m/s 7.010 5 m/s Kh6 3.710 4 m/s 9.010 5 m/s Kv6 3.710 5 m
13、/s 9.010 6 m/s 有效入渗率 30 mm/年 25 mm/年 采用 PEST 进行自动校正，这是一个参数估计程序。采用 Gauss-Marquardt-Levenberg运算法则，通过 PEST 将均方差之和最小化，采用了两个约束条件。第一个约束条件是每一层水平方向上的水力传导系数都相等；第二个约束条件是选择每个参数的最大值和最小值。将初始值除以 10 作为下限，将初始值乘以 10 作为上限。在本研究过程中，通过自动参数估计与试错校正获得的平均绝对误差相同。自动参数估计并没有减小误差，因此，在进一步分析工作中，采用的是试错校正。 （五）结 果 计算东西向剖面的测压水位，结果表明：第
14、一层在河流西边的位置上，地下水补给河流。在河流东边，没有明显的地下水径流。应当注意到，出于计算的原因，模型中 Kraszna河的深度要大于实际观测值。如果减小模型中的河流深度，污染物就会在河流上游方向流入 Kraszna 河，一些污染物会沉积在河流底部。在第一层与第五、六层之间，垂向地下水径流量不大。在第六层，抽水井在地下水流动过程中起着重要作用。 采用 MODFLOW 的区域预算模块计算稳定态的水均衡，可以了解该区不同地层、河流、入渗、抽水和边界之间能量的相互作用。水均衡误差要 控制在 1%之内。在第一层，输入的水量主要是源自侧边界，入渗占总输入量的 17%，大部分输出的水量补给了河流，另一
15、部分输出量是沿侧边界排泄，从第一层通过侧边界流向第二层的水量为 8208m3/天。第二层是粘土层，没有侧边界流量，从第一层进入的水量沿垂直方向直接进入第三层。第三层是细砂层，从第二层获得的水量几乎全部流入第四层，沿侧边界输入的水量和输出的水量相差不大。第四层也是粘土层 ,没有侧边界流量，从第三层获得的水量直接流向第五层。第五层获得的水量通过侧边界和垂直方向流向第六层，第五层中的抽水井也占一部分输出量。在 第六层，通过侧向流获得的输入量为 4218m3/天，从第五层获得的输入量为 5252m3 /天，主要是通过抽水井输出水量，抽水量为 6181m3/天。这说明通过抽水井获取的水量中，有一部分是来
16、自垂直补给，仅通过侧边界补给，不能满足 6181m3/天的抽水量。 四、运移模型 采用两种不同的方法，对运移情况进行模拟，这两种方法分别是采用 MODPATH 进行颗粒示踪（ Pollock， 1994）和采用 MT3D 进行运移模拟（ Zheng and Wang， 1999），包括水平对流和弥散运移。没有考虑由于扩散造成的运移，因为在高渗透 性的环境中，扩散引起的运移量可以忽略不计（ Carges and Baehr， 1998）。由于没有关于污染物阻滞因数的资料，假定阻滞因数为 1，这是比较安全的假定。 （一）边界条件 在运移模型的边界，假定浓度梯度。因为没有 3 种主要污染源（即城市垃
17、圾处理场、污水氧化池和污水处理厂）不同污染物浓度的资料，因此，任意选择 1000 个连续的浓度值。 （二）运移参数 每层的主要输入性质是有效孔隙度以及纵向和横向弥散性。根据前人研究的文献，将有效孔隙度统一选择为 0.10（ Anderson and Woesner， 1992）。 确定弥散性时要复杂一些，弥散值与测试或观测范围有关（ Zheng and Bennett， 1995）。在运移模型中选择的栅格大小为 50 和 100m。根据 Gelhar 等（ 1992）的研究结果，对于 50m 的栅格，纵向弥散度约为 0.3m，对于 100m 的栅格，纵向弥散度约为 5m。为了简化，选择整个研究
18、区的纵向弥散度为 5m。由于资料不足，选择横向弥散度比纵向弥散度小一个数量级，垂直方向的横向弥散度比纵向弥散度小两个数量级。这样在运移模型中，水平横向弥散度为 0.5m，垂直横向弥散度为 0.05m。 （三）结 果 通过模拟， 可以看出，颗粒（即污染物）迁移的深度不大，到达的最深处只有 11m，但是沿水平方向迁移量较大，到达河流附近。根据计算结果，最先到达河流的颗粒来自污水处理厂，经过 10 年的时间到达河流；最后到达河流的颗粒来自城市垃圾处理场，需要 18 年的时间到达河流。颗粒不会达到深井，因此不会污染饮用水源。 MT3D 运移模型表明，经过 8 年左右的时间，污水处理厂的污染物会到达河流
19、；经过 9 年时间，污水氧化池的污染物到达河流；经过 13 年左右的时间，城市垃圾处理场的污染物到达河流。计算主要污染源下 5、 15、 25 和 35m 深度的污染物浓度，结 果表明，在主要污染源以下 15m，污染物浓度低于污染源以下 5m 处污染物浓度的几倍，深度为25m 和 35m，污染物的浓度相当低。根据运移模型也证实了颗粒示踪的结果。根据可利用的数据进行研究，污染物似乎不会到达较深的地区，因此，位于下更新统的饮用水井也不会受到这些污染物的威胁。 五、敏感性分析 在本研究中，由于资料有限，所以对边界条件进行了简化和假定。这样当然会影响结果的准确性，而且污染物不会对饮用水井造成威胁的结论
20、也会遭到质疑。为了检验结论是否准确，选择最差情况进行敏感性分析。 首先，对第一层到第二层和第五层到第 六层的边界条件、水力传导系数、河流参数以及沿垂直方向的水流量进行分析。不同层之间的垂直地下水流动对于溶解污染物向晚更新统地层迁移起着重要的作用。根据敏感性分析结果可以看出，从第一层到第二层的地下水流对所有的边界条件都非常敏感；如果将第一层的水头增加 2m，就会使第一层到第二层的水流量增加 24%；如果将第五层和第六层的水头增加 2m，就会使从第一层到 2层的水流量减少 24%。从第五层到第六层的水流量取决于第五层和第六层的边界条件。如果将第五层和第六层的边界水头降低 2m，就会使从第五层到第六
21、层的水流量增加11%；如果 将第五层和第六层的边界水头增加 2m，就会使地下水流量减少 19%。通过垂直流量与水力传导系数、电导率和入渗的敏感性分析，可以看出，改变粘土层的水力传导系数 K2和 K4，会对水流量造成很大的影响。如果 K2增加 10 倍，从第一层到第二层的水流量就会增加 135%；如果 K4 增加 10 倍，就会使第五层到第六层的水流量增加 61%。 其次，需要分析模拟出来的运移时间和浓度对第一层的水力传导系数、有效孔隙度和弥散程度的敏感性。期间用到了以下缩写： t1：污水处理厂的污染物到达河流的时间；t2：污水氧化池的污染物到达河流的时间； t3： 垃圾处理场的污染物到达河流的
22、时间； c1：污水处理厂 15m 以下的污染物浓度； c2：污水氧化池 15m 以下的污染物浓度； c3：垃圾处理场 15m 以下的污染物浓度。 计算结果表明，随着水力传导系数增加，有效孔隙度减少和弥散度增加，污染物向河流的运移时间减小，水力传导系数较大会造成污染源以下的溶质浓度减小，这是由于如果孔隙介质渗透性较强，污染物更容易在水平和垂直方向上运移。有效孔隙度对浓度分布情况影响不大，而弥散度增加会使污染源以下污染物的浓度增加，而且会造成污染物沿与水流垂直的方向运移，使污染物到达较深的地方。 需要 注意的是，在敏感性分析当中，参数值变化时，无需重新进行校正。但是，这样会得出一些有趣的结果，也许
23、这就是将来需要进一步研究的课题。例如，假设天然补给量是 50mm/年，而不是 25mm/年，重新校正后第一层的水力传导系数可能会变大，根据这一校正值重新评价水流量分布和污染物运移可能是一种新的尝试。 六、最差情况分析 根据敏感性分析，大多数参数对污染物向下迁移的影响都是已知的。利用给定输入参数值（实际值或可能值），通过将污染物向下游迁移速度最快，来确定最差情况，参数值见表 4。如前所述，第一层边界的测压水头增加了 2m， 第五层和第六层的边界条件减小了 2m，粘土层的水力传导系数增加了 10 倍，这就意味着水平和垂直水力传导系数分别为 110 6m/s 和 110 7m/s。所有纵向弥散度都增
24、加了 5 倍，这样，该值现在为 25。所有其它输入参数和变量都为初始值，因为这些值受污染物向下迁移的影响都不大。 在最差情况中，溶质仍会向河流迁移，但是不会迁移到晚更新统的饮用水井中。然而，污染物会到达较深的地方，一部分污染物会到达地下 77m 的地层，这样，就会穿透第二层数米。通过最差情况说明，污水处理厂、污水氧化池和垃圾处理场的污染物不可能到达晚更新 统的生产井中，污染源没有处在生产井的截获区中。 表 4 最差情况的参数值 参数 相对于初始值的校正值 第一层的边界条件 2m 第三层的边界条件 初始值 第五层和第六层的边界条件 2m K1 初始值 K2 10 K3 初始值 K4 10 K5
25、初始值 K6 初始值 弥散度 5 七、结果和讨论 本研究的主要目标是确定污水处理厂、污水氧化池和垃圾处理场的溶解性物质是否可以到达研究区晚更新统的饮用水井中。通过地下水流动和运移模型的模拟，结果表明，生产井不会受这些污染源的威胁。然而，由于 资料有限，模型中的边界条件和参数具有不确定性，这样当然会导致模拟结果的不确定性。为了使研究的结果更为可靠，同时还进行了敏感性分析和最差情况分析。通过敏感性分析，研究了边界条件和参数值对污染物向下游迁移的影响；通过最差情况分析，确定污染物向下所能达到的最大深度。即使是最差情况中，来自污染源的污染物也不会到达饮用水井中。 这一研究表明，敏感性分析和最差情况分析是处理水文地质模型不确定性的较为合适的方法，因为这样可以使根据有限资料模拟的结果变得更为可靠。 译自 Environ Geol 2006, 50
