京津冀典型区域地下水污染风险评价方法研究

李 翔, 汪 洋, 鹿豪杰, 洪 慧, 肖 超, 李 娟, 席北斗*

1.中国环境科学研究院，国家环境保护地下水污染模拟与控制重点实验室，北京 100012 2.生态环境部土壤与农业农村生态环境监管技术中心，北京 100012

在京津冀地区，75%以上城镇生活饮用水均来自地下水，地下水环境质量是京津冀地区的饮水安全与社会经济发展的重要保障[1]. 随着对地下水的持续开采，京津冀局部地区地下水位持续下降，由此引发的各种环境问题日益突显. 由于地下水系统的复杂性与隐蔽性，地下水污染问题难以发现，且一旦污染很难治理[2]. 因此，通过开展地下水污染风险评价方法研究，寻求适合京津冀地区的地下水污染风险评价方法，识别地下水易遭受污染的区域，以分区域、分级别防范地下水污染，可为管理者合理利用地下水资源、有针对性地制定地下水保护措施提供理论支撑.

地下水污染风险也被称为“地下水的易污性”“地下水系统防污性能”“地下水污染潜力”等[3-4]. 地下水污染风险问题最早由法国水文地质学家Margat于1968年提出[5]，提出伊始，人们仅关注到土壤-含水层系统本身，即仅研究自然因素对地下水脆弱性大小的影响，认为地下水脆弱性为一种含水层内在的属性，主要与包气带的保护能力及饱和带的净化能力有关. 1987年在“土壤与地下水脆弱性国际会议”上，首次提出地下水的脆弱性与人为因素有关，认为不同的污染源，地下含水层对其敏感性大小也不同，因此其地下水污染风险也不同[6]. 目前国内外常见的地下水污染风险评价的方法主要有迭置指数法、过程模拟法、统计方法和模糊数学方法等[7].

迭置指数法由于具有评价方法简单、评价参数易于获取等优点，是目前最简单且应用最广泛的地下水污染风险评价方法. 但是由于因子的评分和权重体系多基于经验方法获取，客观性和科学性较差. 典型的评价模型有DRASTIC、GOD、SINTACS等[8-14]，其中DRASTIC模型的使用较为普遍[15]. 统计方法虽然能够客观地筛选出地下水污染的主要影响因素，并通过回归方程给出适当的权重值，避免了专家评判的主观性，但其没有涉及污染发生的基本过程并且评价时需要大量的监测资料和相关信息作为基础，应用有限[16]. 常用的统计方法有地理统计法、线性回归分析法、逻辑回归分析法和实证权重法等[17-18]. 过程模拟方法是利用成熟的污染物迁移转化模型，模拟污染物在包气带和饱水带中迁移、转化的过程，结合地下水脆弱性评价公式计算地下水脆弱性综合指数. 过程模拟法既能模拟地下水污染物的迁移转化过程，又可以预测污染物的时空分布情况，且污染物迁移时间、污染物浓度及污染面积等模拟结果均可进行量化. 由于评价过程中没有主观因素的影响，评价结果科学性和可信度较强，近年来受到越来越多的关注. 目前，常用的污染过程模拟软件和模型主要有FEFLOW、FEMWATER、GMS、Visual MODFLOW、MT3D和RT3D、HYDRUS等[19-22].

HYDRUS-2D是一款模拟饱和-非饱和区水分、热及溶质运移的软件，具有输入数据相对较少、模拟结果可靠等优点，在包气带土壤水分及溶质运移模拟研究中被广泛应用[23-25]. 因此，该研究以京津冀地区某典型污染区域为例，采用HYDRUS-2D软件，通过模拟污染物在不同包气带结构下的迁移转化过程，对地下水污染风险进行评价研究. 研究成果可为管理者合理利用地下水资源，有针对性地制定地下水保护措施提供理论支撑.

1 研究区概况

研究区位于京津冀平原区内，面积约为909.72 km2，地形由北向南倾斜，平均坡度为0.1%～0.2%. 气候类型属于典型暖温带大陆性季风气候，四季分明，春季干旱多风，夏季高温多雨，秋季天高气爽，冬季寒冷干燥.

研究区第四系沉积物广泛分布于平原和山间沟谷，沉积层的特征为北薄南厚、东薄西厚. 根据松散层的沉积规律和埋藏条件，研究区含水层由上到下可分为潜水含水层、浅层承压含水层和深层承压含水层组. 潜水含水层广泛分布于研究区，该层岩石颗粒大、渗透性强，直接接受大气降水和地表水入渗补给，为研究区农业用水的主要取水段；浅层承压含水层为多层含水层分布区，含水层埋深范围为50～120 m，其岩性均一、透水性强、富水性好，是当地水源井、自备井的主要开采层. 地下水补给主要为大气降水入渗补给、侧向径流补给、水库渗漏、河流入渗补给以及灌溉回渗补给，其中大气降水入渗补给是研究区内地下水最主要的补给来源. 天然条件下，地下水径流方向与地形地貌变化一致，即由山前向平原，由北西向南东流动. 现状条件下，由于地下水的集中开采，局部地区流向发生改变，在开采量较大的地区，形成区域性降落漏斗. 研究区地下水的排泄方式主要有两种：①自然排泄，主要是指地下水的溢出、蒸发及流向下游的侧向流出；②人工开采，包括应急水源地、区县水源地及工农业自备井开采. 天然条件下，地下水的排泄方式主要为蒸发排泄; 现状条件下，人工开采已成为地下水主要排泄方式.

2 地下水污染荷载评价

基于现场调研结果，研究区地下水污染源主要包括工业源、农业面源、居民区、垃圾填埋场和地表水体5种类型，污染源分布情况见图1. 根据各类型污染源的污染物排放方式，以硝酸盐为特征污染物计算污染荷载.

2.1 工业源

结合环统数据和现场调研资料，通过式(1)进行工业源污染荷载评价.

Qn=q×Qi×λ/Si

(1)

式中：Qn为单个工业源污染荷载值，mg/(m2·a)；q为特征污染物硝酸盐排放浓度，mg/L；Qi为污染源i的污水排放量，L/a；λ为地表入渗系数；Si为污染源i的面积，m2.

由于缺乏工业源的污水排放监测数据，参考各类型工业源排放标准中的最大允许硝酸盐排放浓度代替硝酸盐的实际排放浓度. 工业源的影响范围参照《地下水污染防治分区划分工作指南》中各类污染源缓冲区半径设定，如表1所示.

2.2 居民区

城镇居民区用水主要包括生活用水、公共服务用水、生产运营用水等，假设城镇居民区的用水全部为生活用水，其供应的水全部转化为城市污水，供水的

表1 工业源缓冲区半径的设定

损失量全部渗漏进入地下，则生活污水的入渗系数计算公式：

(2)

式中：λr为生活污水入渗系数；Q1为居民区供水量，m3/a；Q2为居民区废水排放量，m3/a；Q3为多年平均降雨量，m/a；S为居民区面积，m2；μ为降雨径流系数.

该地区2018年水资源公报显示，2018年该地区总供水量为39.3×108m3，污水排放总量为20.4×108m3，耗水量为18.9×108m3，计算出该地区生活污水入渗系数(λr)为0.77. 根据《全国水环境容量核定技术指南》，一般城市人均产污系数约为6 g/d，由此计算出居民区硝酸盐污染荷载(见表2).

表2 居民区硝酸盐污染荷载计算结果

2.3 农业面源

研究区农业面源污染物的排放浓度根据《全国水环境容量核定技术指南》确定. 《全国水环境容量核定技术指南》中规定，标准农田源强系数COD为15 t/(km2·a)、氨氮为3 t/(km2·a). 研究区主要为平原区，表层土为黏粉和粉黏，其年均降雨量为567 mm，经过研究区实际情况修正后，硝酸盐源强系数为2.64 g/(m2·a).

2.4 垃圾填埋场

垃圾填埋场的硝酸盐污染荷载，主要根据垃圾填埋场渗滤液产生量及渗滤液中硝酸盐浓度确定. 假定填埋场中填埋的垃圾水分含量为零，其渗滤液的产生量仅与年降雨量有关，根据式(3)可估算垃圾渗滤液产生量[26].

Qf=K×I×A×10-3

(3)

式中：Qf为垃圾填埋场渗滤液流量，m3/d；K为渗出系数；I为最大年降水量的日换算值，mm/d；A为垃圾填埋场面积，m2.

由于缺乏垃圾填埋场渗滤液检测数据，参考GB 16889—2008《生活垃圾填埋场污染控制标准》中现有和新建生活垃圾填埋场水污染物排放质量浓度限值，计算硝酸盐污染荷载，结果见表3.

表3 垃圾填埋场硝酸盐污染荷载计算结果

2.5 地表水体

地表水体硝酸盐污染荷载主要与河道水质、水量及底泥的渗透性有关，由于河道的面积难以确定，通过其河道的长和宽确定河道的水量. 河道硝酸盐污染荷载计算公式：

Qw=L×W×Km×C

(4)

式中：Qw为河道硝酸盐污染荷载，t/a；L为河道的长度，m；W为河道的宽度，m；Km为底泥入渗速率，参考淤泥质土的入渗速率，取值为1.5×10-9m/s；C为河道中硝酸盐浓度，mg/L.

由于研究区内的河道多为再生水排放河道，因此根据GB 18918—2002《城镇污水处理厂污染物排放标准》确定河道的硝酸盐浓度限值，计算各河道硝酸盐污染荷载.

图2 研究区硝酸盐污染荷载空间分布Fig.2 Spatial distribution of nitrate pollution load in the study area

2.6 污染荷载叠加计算

根据以上5类污染源的污染荷载量化计算结果，基于ArcGIS 10.2软件中的栅格叠加功能将不同污染源的硝酸盐污染荷载计算结果进行叠加，得到研究区内各类污染源硝酸盐污染荷载空间分布(见图2). 由图2可见，研究区硝酸盐污染荷载主要受垃圾填埋场及工业源分布的影响. 其中1号垃圾填埋场为2019年研究区重点排污单位名录中固体废物重点监管单位，垃圾渗滤液中的氨氮含量较高，造成较高的污染荷载. 2号垃圾填埋场建设时间较早，随着垃圾量的增多，出现了超负荷运营的情况，目前该垃圾填埋场已封场，但场地地下水依旧存在较高的潜在污染风险.

3 包气带污染过程模拟

地表污染物随着下渗水进入包气带运移过程中，可在包气带介质发生各种衰减作用[27-29]. 将流出包气带介质的污染物质量与进入包气带介质的污染物质量之比定义为包气带折减系数，用其来表征污染物在包气带介质中迁移转化的过程. 基于研究区钻孔资料对包气带结构进行概化分区，通过HYDRUS-2D软件模拟硝酸盐从地表到地下的过程，进而计算得到每个包气带介质分区下的折减系数.

3.1 包气带结构概化分区

包气带中硝酸盐的迁移转化过程与包气带介质的岩性、厚度、结构密切相关[30-31]，而实际环境中包气带结构复杂，为了便于计算，对其进行概化分区. 通过对研究区钻孔资料、水文地质图等资料的收集分析，将研究区概化为8个不同包气带结构分区，如图3所示. 各分区下的包气带结构概化结果见图4.

图3 研究区包气带结构概化分区Fig.3 Partition diagram of vadose zone in the study area

图4 研究区包气带结构概化结果Fig.4 Lithology generalization diagram of vadose zone structure

3.2 污染物包气带运移模拟

通过HYDRUS-2D软件，模拟硝酸盐在每种包气带结构下的迁移转化过程，并基于模拟结果计算包气带折减系数. 硝酸盐在不同包气带结构下的折减系数根据式(5)计算.

Rj=Cj/C0

(5)

式中：Rj为包气带结构分区j下的硝酸盐折减系数；Cj为硝酸盐穿过包气带到达地下水面时的浓度，mg/L；C0为硝酸盐进入包气带时的初始浓度值，mg/L.

3.2.1不同初始浓度模拟计算结果

为探究硝酸盐初始浓度(C0)是否对包气带结构分区中硝酸盐折减系数(Rj)有影响，以1区为例，在包气带结构不变的情况下，设置5种情景模式，模拟不同污染荷载条件下，硝酸盐穿过包气带到达含水层的浓度(Cj). 5种情景模式的硝酸盐初始浓度(C0)分别为5、10、15、20、50 mg/L，其初始条件、边界条件等不变. 通过模拟得到10年后硝酸盐到达含水层的浓度(Cj)，根据式(5)计算出硝酸盐折减系数，其结果见图5.

图5 1区不同硝酸盐污染荷载下垂向迁移的模拟结果Fig.5 The vertical migration simulation results of different pollution loads of nitrate in zone 1

从模拟结果来看，在同一包气带结构下，情景1～情况5进入地下水面中的硝酸盐浓度分别为3、6、9、12、30 mg/L，硝酸盐折减系数均为0.6. 模拟结果表明，同一包气带结构下，不同硝酸盐初始浓度仅影响到达地下水面的硝酸盐浓度，并不会改变垂向折减系数.

3.2.2不同包气带结构模拟计算结果

假设硝酸盐的初始浓度为20 mg/L，即不改动模拟环境中的水分边界，将溶质边界设定为20 mg/L. 通过HYDRUS-2D软件分别模拟硝酸盐在研究区各包气带结构分区下的迁移转化过程，并计算各包气带结构分区下的硝酸盐折减系数，如图6所示. 模拟结果表明，不同包气带结构分区下其硝酸盐折减系数不同，按照从大到小的顺序分别为R4

图6 不同包气带结构分区的硝酸盐折减系数Fig.6 The reduction coefficient diagram of nitrate in different vadose zone

4 地下水污染风险评价

基于研究区地表硝酸盐污染荷载计算结果(见图2)和各包气带结构分区下的硝酸盐折减系数计算结果(见图6)，通过ArcGIS 10.2软件中的栅格计算功能，将地表硝酸盐荷载量与折减系数相乘，计算硝酸盐穿过包气带介质后到达地下水面的总量，来表征研究区地下水硝酸盐污染风险，评价结果见图7. 由图7可见，地下水硝酸盐污染风险高的区域呈现零星的分布趋势，主要受垃圾填埋场和工业园区的分布影响，此外由于研究区内地表水体多为再生水排放河道，河道周边地区地下水硝酸盐污染风险也相对较高.

图7 研究区地下水硝酸盐污染风险评价结果Fig.7 Risk assessment chart of groundwater nitrate contamination in the study area

5 结论

a) 研究区硝酸盐污染荷载按从大至小的顺序分别为垃圾填埋场、工业源、农业面源、地表水体、居民区.

b) HYDRUS-2D软件模拟结果表明，包气带结构不变的情况下，污染源荷载的变化只会导致进入含水层中的污染物浓度不同，不会改变包气带硝酸盐折减系数.

c) 污染源类型、包气带介质岩性及厚度是造成研究区内地下水硝酸盐污染风险评价存在差异的主要原因.