采用敏感性分析和最差情况分析进行地下水污染风险评价

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