1 引言(Introduction)

地下水污染风险评价对地下水的管理与保护起着重要作用.英国地质学家Morris等(1987)指出, 地下水污染风险的定义是含水层中地下水由于其上的人类活动而遭受污染, 其程度到不可接受的水平的可能性.部分学者(Doerfliger et al., 1999)还考虑到地下水被污染后其价值功能可能发生的变化.国内外对地下水污染风险评价形成了各种体系, 从地下水污染风险评价的定义及其发展可以看出, 地下水污染风险评价主要是基于含水层脆弱性评价进行(Kumar et al., 2015;Hamza et al., 2015;邰托娅等, 2012), 而对含水层脆弱性评价使用最多的是叠置指数法(刘凡等, 2014).目前使用最为广泛的方法是基于美国环保局(EPA)1987年提出的DRASTIC模型或者其改进模型进行地下水脆弱性评价(Aller et al., 1987), 以及结合污染荷载、地下水价值进行地下水污染风险评价(申丽娜等, 2010; Ouedraogo et al., 2016).

传统的地下水污染风险评价体系以1988年世界卫生组织(WHO)与泛美卫生组织(PAHO)联合提出的方法(Forster et al., 1988)为代表, 简称WP体系, 评价结果受人为因素影响较大, 且评价结果不够精确.该体系将污染负荷与地下水防污性能分开考虑, 未能充分刻画出污染物在包气带中迁移转化的过程.

针对地下水污染风险评价体系存在的问题, 本文提出可以更准确刻画污染物在包气带中衰减过程和污染物输移过程的方法体系, 通过引进“折减系数”概念, 充分考虑包气带岩性、结构、厚度的差异、污染物属性特征等因素对污染物在包气带中迁移转化过程的影响, 对污染物在包气带中的衰减过程进行刻画.基于该体系方法, 使得污染物穿透包气带到达地下水这一过程的量化变得简洁、易操作.

2 方法(Methods) 2.1 评价方法体系

地下水的污染是一个复杂的过程, 首先需要明确地表污染源的类型.根据污染源类型, 对地表污染物排放量进行计算.地表污染物穿透包气带介质进入地下水的这一过程, 涉及到水文地质环境、污染物降解性、迁移性等诸多因素.本研究提出了“折减系数”来量化污染物在包气带中衰减的过程, 将地表污染物量与其对应的折减系数相乘得到其到达地下水面的量.再把到达地下水面的污染物量与该污染物毒性结合, 利用ArcGIS9.3软件得到基础地下水污染风险评价图.

地下水价值是其重要性的体现(孙才志等, 2015), 反映了地下水被污染后其价值功能可能发生的变化, 即一旦被污染, 地下水价值越大的地方后果越严重.将基础地下水污染风险评价和地下水价值按照叠置指数法叠加即可得到地下水污染风险分级, 该方法体系框架如图 1所示.

该体系分为4个层次, 首先, 根据包气带介质分区及每种污染物迁移属性得到每一种污染物对应每一个包气带分区的折减系数；然后, 由该污染物进入包气带的量与折减系数相乘得到该污染物到达地下水面的量；其次, 将到达地下水面的污染物量与污染物危害性叠加得到基础地下水污染风险；最后, 将基础地下水污染风险与地下水价值叠加得到地下水污染风险.该体系重在刻画污染物在包气带中的衰减过程, 是从污染物通过包气带时量的衰减和危害性不变的本质层面考虑的.

该评价方法体系与传统体系的不同主要体现在两方面.首先，在该体系中, 把污染物的量和属性分开来考虑, 把污染物属性分为毒性和影响污染物在包气带中迁移转化的属性.污染物的危害性与污染物到达地下水的量放在同一层次叠加计算, 这样考虑污染荷载更加合理, 把污染物对地下水产生的影响刻画的更清楚.其次, 折减系数实现了对污染物在包气带介质中迁移转化过程的量化, 既充分考虑了包气带介质的不同, 也考虑了每一种污染物性质的差异.与数值模拟法相比, 该方法更加简洁易操作；与叠加指数法(Gogu, 2000)相比, 该方法更精确.

2.2 折减系数法的原理及推导

地表污染物随着下渗水进入包气带垂向运移过程中, 可在包气带介质发生各种衰减作用, 其中, 吸附和降解是最主要的作用(任剑飞等, 2015).这里将流出包气带介质的污染物质量与进入包气带介质的污染物质量之比定义为污染物穿透包气带时的折减系数.显然折减系数的大小不仅取决于污染物本身的物理化学性质, 还取决于包气带的岩性、结构、厚度等因素.为了获取折减系数, 可利用Multi-cell模型(Kong et al., 1996；Wang, 1999)结合推流的概念, 建立一个模型.按照介质结构顺序, 把包气带介质概化成多个连续的单元, 每一个单元作为一个独立的均质结构单元.首先假设研究区范围内的某一种污染物均匀地铺在地表, 其基本过程为：假设污染物随水垂向运移过程遵循推流模式, 在自地表向下进入每一个单元时, 水中的污染组分会与包气带介质产生吸附分配作用和降解作用, 在进入到下一个单元格后会以上一个单元格流出水溶液的浓度为起始浓度再进行分配和降解, 依次往复, 最终从最后一个单元中流出的液相浓度即为污染物穿过整个包气带之后的出水浓度.

依据质量守恒定理, 即：m_lr=m_lc+m_a+m_d(其中, m_lr为流入介质内的污染物质量, m_lc流出的污染物质量, m_a吸附在介质内的污染物质量, m_d降解反应降解的污染物质量), 可推导出相应的计算公式.

按照污染物本身的物理化学性质, 可将污染物在包气带中发生反应的类型分为4种模式：① 不发生吸附反应也不发生降解反应；② 只发生吸附反应；③ 只发生降解反应；④ 既发生吸附反应又发生降解反应.4种反应类型依次简称为模式1、模式2、模式3、模式4.结合包气带的岩性、厚度、颗粒大小, 量化出每一种污染物经过包气带后, 包气带对其的消减作用, 即“折减系数”.折减系数推导过程如下.

污染物进入包气带介质中, 依据质量守恒定理可得：

(1)

式中, M_in为在t时间内进入单元体的溶质质量；M_out为在t时间内流出单元体的溶质质量；M_s为t时间内由于吸附反应消耗的溶质质量, 模式1、模式3中为0；M_d为t时间内由于降解反应消耗的溶质质量, 模式1、模式2中均为0.

假设在理想条件下, 包气带介质是由均匀介质构成的, 溶质经过t时间从地表随水流下渗穿透包气带介质.对于模式2, 本研究中将污染物发生的吸附反应认为是瞬时平衡的可逆线性吸附, 符合线性吸附公式；对于模式3, 降解反应主要考虑一级生物降解反应；对于模式4, 假设包气带介质中溶质的吸附和解吸能够迅速地达到平衡, 属于线性等温吸附.并且, 在包气带土壤水介质中发生的生物降解反应符合一级衰减动力学方程.于是可以得到：

(2)

(3)

(4)

式中, n为土壤孔隙度；ρ为土壤介质容重；V为垂向单位土壤介质的体积；C₀为初始时刻液相浓度；C₁为t时间单元体的液相浓度；S₁为t时间内单元体的固相浓度, 模式1中为0；D₁为t时间内单元体消耗的浓度, 模式1、模式2中均为0；K_d为达到吸附平衡时固相和液相污染物的分配系数；K为降解速率常数.

折减系数W的计算公式如式(5)所示, 而污染物的阻滞系数R=1+ρK_d/n, 因此, 式(5)可以化简得到式(6).

(5)

(6)

然而, 实际情况中包气带介质垂向结构十分复杂, 由不同类型的土壤介质组成, 每一种土壤介质都有不同的阻滞系数, 而且每一种土壤介质的厚度及渗透系数都不同, 所以, 污染物穿透每一层土壤介质的时间t也不同.按照土壤介质结构顺序, 把包气带介质分成多个部分, 每一种介质作为一个单独的均质结构, 使用上述均质公式(4)进行计算, 将在包气带介质中发生吸附作用(或者吸附作用和降解作用)后流出的该部分的污染物液相浓度作为进入下一个均质部分的初始浓度进行计算, 其他各个土壤介质遵循同样的方法, 最终从最后一个土壤介质类型中流出的液相浓度即为污染物穿透整个包气带厚的出水浓度.由于每个水文地质区包气带介质的厚度也不同, 因此, 本研究中将计算出的折减系数进行归一化处理, 污染物在每一个水文地质分区内的归一值为该污染物穿透包气带介质后的量与水流穿透该包气带介质厚度的时间t的比值.

在包气带介质中既发生吸附反应也发生降解反应的污染物穿透包气带介质时的折减系数W推导公式为：

(7)

式中, 对于模式1, R_i、K始终为0；对于模式2, K始终为0.各模式的折减系数公式具体见表 1.

2.3 包气带介质分区概化

包气带介质的岩性、结构、厚度均会影响污染物在其中迁移转化的过程, 而实际中包气带介质结构复杂, 为了便于计算, 以及实现量化过程的可操作性, 需要对其进行分区概化.基于所收集资料, 通过研究区已有的钻孔剖面资料对研究区进行包气带介质分区, 利用钻孔剖面与研究区水文地质图结合研究区包气带形成过程, 考虑研究区地形地貌特征, 推断不同包气带介质中最具典型的包气带岩性与厚度, 从而实现对研究区包气带水平分区和纵向分层.

3 材料(Materials) 3.1 研究区概况

石家庄地处河北省中南部, 位于北纬37°26′~38°46′, 东经113°31′~115°29′之间.本次研究区主要包括石家庄市区、鹿泉市、正定县、灵寿县与平山县区域, 研究区位置与范围如图 2所示.区内最大河流为滹沱河, 所形成的河床及河漫滩一般宽度为2000~4000 m.由于上游修建了黄壁庄水库, 河水受到人工控制, 所以通常滹沱河河道内只有汛期才有水通过.研究区内有黄壁庄水库与岗南水库两大水库, 岗南水库是河北省管辖的库容最大的水库, 岗南水库与黄壁庄水库主要作用以防洪为主, 并结合发电、农业灌溉、工业供水和省会石家庄城市生活、环境供水等为一体, 为石家庄最重要的水源地.

3.2 地下水污染风险评价 3.2.1 污染荷载

研究区主要污染源类型有工业污染源、生活污染源、农业污染源、垃圾填埋场及排污河5类, 地下水中主要的污染物在这5类污染源类型中都会产生, 因此, 认为地表每一种污染物的排放量是这5类污染源排放的共同结果.实测代表性污染源中各种污染物的排放浓度, 根据每种污染源的情况对该污染源污染物的排放量化, 其计算公式为：

(8)

式中, Q_ij为实测污染源j污染物i浓度(mg·L^-1), λ为地表入渗系数, S_j污染源j区域面积(m²)；L_ij表示污染源j单位面积上污染物i浓度(mg·m^-2·a^-1), Q_j表示污染源j污水排放量(L·a^-1), 具体对每一种污染源而言计算方法不同(李广贺等, 2015).这5类污染源产生的特征污染物共有21种, 分别为：COD、氨氮、亚硝酸根、氯离子、氟、硝酸根、硫酸根、钠离子、铁、锰、砷、苯并(a)芘、荧蒽、芘、芴、苊、萘、甲苯、三氯甲烷、六氯苯、1, 2二氯苯.利用ArcGIS9.3软件分别对每一种污染物在5类污染源类型中的排放量进行矢量图绘制, 并利用软件中的栅格叠加功能将该种污染物在每一种污染源类型中的排放量进行叠加计算, 最终得到市研究区污染物地表分布, 结果见图 3.

3.2.2 包气带介质分区概化

研究区地下水主要赋存于多层交叠第4系松散沉积物中, 水文地质条件具有明显的水平和垂直分带性.自西向东含水层介质粒径变小, 含水层层次逐渐增多, 富水性由强变弱；在垂向上, 含水层的上部及下部砂层粒度较细, 厚度较小, 中部砂层粒度较粗, 厚度较大.岩性结构也由最初的单一砂、卵砾石层过渡至顶部有薄层粘性土覆盖的砂、卵砾石层, 最后变成粘性土夹杂着多层砂卵砾石层.地下水径流条件由强变弱, 富水性由好变差.污染物在包气带中迁移转化的程度与进入地下水中量的多少取决于包气带的岩性、结构、厚度.选择2条有代表性的剖面, 剖面岩性如图 4a、4b所示(注：均根据地质矿产部河北水文工程地质勘察院(2004)简化绘制).研究区主要由滹沱河的冲洪积扇形成, 滹沱河两侧包气带岩性存在一定相似性, 在岩性剖面没有覆盖的区域参考滹沱河对称区域包气带.研究区包气带厚度分布特征由研究区地下水位埋深得到.

对石家庄平原区包气带介质进行分区概化, 共概化为9个不同类型包气带介质(图 4c), 分别为A、B、C、D、E、F、G、H和I区, 对于每个地质区包气带岩性剖面参考相应区域整体包气带属性来代替整个水文地质区的岩性.岩性剖面用典型的钻孔岩性来代替整个水文地质区的岩性, 包气带岩性分区概括如图 4d所示.

3.2.3 折减系数计算

研究区的包气带介质概化为9个不同的水文地质分区, 依据折减系数的定义, 对于每一种污染物, 在9个包气带介质分区对应着9个不同的折减系数.包气带介质的土壤容重参考各岩性土壤密度经验值(唐大雄等, 1999)赋值；渗透系数参考松散岩石渗透系数, 孔隙度参考文献(王大纯等, 1995)赋值；由于分配系数K_d随包气带介质中有机碳含量的变化而变化, 实际情况中介质中的有机碳含量较难直接获取, 考虑用土壤介质中粘粒的百分含量代替(何江涛等, 2012).为了便于说明, 在本研究考虑的污染物中选择COD为例计算折减系数, 结果如表 2所示.

COD在土壤水环境中会发生吸附、沉淀反应, 且容易被生物降解, 其中, 生物降解是导致COD降低的主要原因.根据文献调研, 可以得到COD在包气带介质中的降解系数为0.0453(罗定贵等, 2009).认为COD在包气带介质中发生吸附反应和降解反应, 按照穿透包气带介质时发生吸附和降解反应对应的折减系数推导公式, 计算得到COD在A、B、C、D、E、F、G、H和I区的折减系数依次分别为0.0274、0.0682、0.1890、0.0731、0.0679、0.0621、0.0801、0.0678和0.0908.因此, COD穿透包气带介质后的量及归一化值按式(9)、(10)计算.

(9)

(10)

式中, Q₁为COD穿透包气带介质后的量(mg·m^-2·a^-1), Q₂为COD地表量(mg·m^-2·a^-1), a为折减系数, Q_t为COD归一化值(mg·m^-2·a^-1), t为水流穿透该包气带介质厚度的时间.

地下水受到污染源污染, 归根结底表现为污染源排放的一定数量的特征污染物透过包气带, 对地下水环境造成破坏.进入含水层特征污染物的数量与污染源的排放数量有关, 而特征污染物对地下水环境的破坏则需考虑其自身属性(Wang et al., 2012).污染物的迁移性与降解性等性质已经包含在污染源量化过程中, 这里考虑特征污染物的毒性, 参考DZ/T 0290-2015《地下水质量标准修订》和GB5749-2006《生活饮用水卫生标准》.其中, 污染物浓度限值越大, 说明其毒性越小, 根据污染物浓度限值从大到小排序, 依次赋值.基础污染风险的计算公式为：

(11)

式中, R为地下水基础污染风险, C_i为第i种污染物毒性评分, Q_i为第i种污染物到达地下水面的量的评分.利用ArcGIS软件的栅格叠加功能将所有污染物对地下水的污染风险叠加起来得到地下水污染风险评分, 利用Natural Breaks进行分级即得到基础地下水污染风险分级(图 5a).

由于资料缺少, 本研究中仅用研究区富水性强弱来表示地下水价值大小, 在本文中, 不影响对折减系数的讨论.由研究区水文地质图得到该区的富水性分级(图 5b).岗南水库到黄璧庄水库下游河道及周围富水性最强, 周围区域较强.

3.3 地下水污染风险评价结果与验证

地下水污染风险评价依据公式Z=R×a₁+V×a₂计算, 其中, Z为地下水污染风险, R、V依次分别为基础地下水污染风险、地下水价值；a₁、a₂分别为地下水基础污染风险、地下水价值的权重值.本研究认为地下水基础污染风险比地下水价值稍微重要, 计算得到的a₁、a₂分别为0.6、0.4.在得到研究区地下水污染风险分值的基础上, 利用Natural Breaks方法实现等级划分, 得到的石家庄市平原区地下水污染风险分级如图 6a所示.

石家庄市平原区地下水污染风险评价结果表明, 地下水污染风险最高的区域位于石家庄市区东北侧, 该区域为老工业基地及其他工业聚集区, 很多重型污染企业分布在此区域, 不同类型污染源排放的相应污水可能是造成该地区地下水风险较高的一个原因.其次为滹沱河河道及周边区域, 受上游水库截流的影响, 研究区内大多数河流常年断流, 形成“有河皆干、有水皆污”的状况, 且河道包气带基本以砂性土为主, 这种包气带结构和岩性有利于地下水的补给, 且该带地下水极富水, 故其脆弱性强.加上滹沱河河道长期受居民区污水影响, 综合各类因素导致该地区地下水污染风险较高.地下水污染风险最低的地区为农田区.整体来说, 包气带介质岩性和厚度的差异、污染物属性特征与分布污染源类型的不同是造成地下水风险评价存在差异的主要原因.其中, 工业污染类型是影响地下水污染的主要污染源, 其次为居民区污水的排放, 农业污染源对地下水的影响最小.

基于对研究区地下水采样测试结果, 采用单指标最大类别法评价每一个采样点处地下水水质, 利用ArcGIS的空间分析功能得到研究区地下水采样点水质分级图如图 6b所示.可以观察到, 地下水水质分布与污染风险分布大致一致, 污染风险较高区域对应点地下水水质也较差, 污染风险较低区域地下水水质较好.

4 结论(Conclusions)

1) 引入“折减系数”实现了污染物由地表到达地下的衰减量化过程刻画, 使得污染物穿透包气带到达地下水这一过程的量化变得简洁、易操作.折减系数推导过程涉及到包气带介质的岩性、厚度、孔隙度、土壤容重, 污染物降解速率等, 可以较为充分地概括包气带水文地质条件及污染物的性质.

2) 该评价体系运用于石家庄滹沱河平原区的结果表明, 地下水污染风险最高的地区多为工业区, 与研究区工业污染源分布区域相符合；其次为滹沱河河道及其附近区域, 河道包气带以砂性土为主, 河道及其附近地区常年有排污河河水下渗, 且该带含水层富水性极好, 故其污染风险较高；大部分农田区污染风险较低, 农田区污染源排放量较小.

致谢(Acknowledgement): 感谢中国地质科学院水文地质环境地质研究所孙继朝教授团队, 以及中国环境科学院刘琰研究员团队在本研究中给予的大力支持.