数值模拟法划分地下饮用水源保护区—以内蒙古呼和浩特市城市水源地为例

2018-12-13张翼龙王文中

水文地质工程地质 2018年6期

王涵，刘琦，张翼龙，王文中

(1.同济大学土木工程学院地下建筑与工程系，上海 200092；2.中国地质科学院水文地质环境地质研究所，河北石家庄 050061)

地下饮用水水源的水质直接影响着人类的健康。为了能让地下饮用水源供水井远离污染，建立地下饮用水源保护区，一方面可以通过一定的距离和时间，让污染物在到达供水井前能够被各种介质所拦截净化，使其浓度降低到目标含量；另一方面，也可以让相关部门在地下水受污染后有一定的应急处置时间，使地下饮用水源的水质达到规定的标准，满足饮用及其他生活用水要求，实现地下饮用水源地的长效发展。

划分地下水保护区的方法有很多，早期地下水保护区划分主要依据经验，随着研究的不断深入，地下饮用水源地保护区的划分方法也逐渐衍生出经验公式法、分析法、圆柱法、水文地质条件法、解析模型法和数值模拟法[1～3]。然而随着经济的发展，城市化进程的加快，土地资源越发稀缺，愈加宝贵，所以水质的好坏固然是重要的考虑因素，但同时也需要追求经济和社会价值的最大化，即应在合理的范围之内让保护区的范围更小[4]。如何更加科学、经济地划分地下饮用水源保护区成为了一个更加有现实意义的研究课题。

数值模拟法[5]适用于水文地质条件复杂，边界条件众多的地区，例如含水层较多且非均质的含水系统、非稳定流和分布多个注水井和抽水井的含水系统等[6～8]。本文主要以内蒙古呼和浩特市为例，运用数值模拟法计算得到地下饮用水源地一级和二级保护区，并根据现实情况提出合理的保护建议，以期为呼和浩特市地下水资源的健康长效发展提供参考。

1 研究区概况

本研究区为呼和浩特市城区，面积约2 286 km2，如图1所示。呼和浩特市城区饮用水源主要为地下水，目前共148口主要取水井，供水能力30.8×04m3/d，周边区域工业企业及混采井分布密集，易受污染影响。因此，科学合理地划分地下饮用水源地保护区对保障呼和浩特市的供水安全，防止地下水污染具有重要的意义。

图1 研究区地理位置Fig.1 Location of the study area

1.1 水文地质条件

呼和浩特市城区地下饮用水源地位于呼和浩特市平原区，主要分布更新统和全新统地层。地下水主要赋存于冲洪积砂砾层中，富水性强。区域主要接受大气降水和地下水径流补给，地下水流向大致为北东—南西向。研究区含水层模型见图2，研究区含水层大致分为北部大青山单一潜水含水层(约占研究区15%)和南部呼和浩特平原潜-承双层结构含水层(约占85%)。其中北部单一潜水含水层岩性以卵砾石、中粗砂为主，一般厚度为120～140 m。地下水位埋藏较深，大部分地段地下水位埋深大于60 m，水量中等，单位涌水量100～500 m3/(d·m)。地下水径流条件较好，水力坡度一般在1‰～3‰。南部双层结构含水层的上部潜水含水层厚度10～27 m，以砂砾石为主，水量丰富或极丰富。该含水层径流条件良好，地下水由东流向南西，水头高度1 018～1 080 m，水力坡度为1‰～2‰，其中北部靠近大青山地区潜水含水层疏干。下部承压含水层上覆有稳定较厚的淤泥层，承压含水层主要由湖滨相粗碎屑物组成，水量丰富，含水层厚度较大，是呼市主要地下水开采层。该含水层北东向埋藏深度较浅，一般埋深40～60 m；南西向埋藏深度逐渐增大，颗粒由粗变细，含水层顶板逐渐变厚至接近200 m。承压含水层地下水由北东流向南西，水头高度1 010～1 060 m，水力坡度为1‰～3‰。

图2 研究区水文地质结构概念图Fig.2 Conceptual diagram showing the hydrogeological structure of the study area

1.2 淤泥层分析

研究区南部大部分地区潜水含水层与承压含水层之间分布有较厚的淤泥层，该层厚度分布如图3所示。

图3 研究区淤泥层厚度等值线图Fig.3 Contour map of thickness of the silty layer in the study area

淤泥层的渗透性直接影响着两层含水层之间的连通性。对此有学者进行了详细分析，根据沉积环境和实验分析，淤泥层厚度小于70 m的区域淤泥层含砂量较高，因此从淤泥层分布边界到厚度70 m等值线之间区域淤泥层渗透性较好，易于发生潜水和承压含水层之间的越流[9]。研究区内大部分取水井分布于该区域，受污染风险较大。

1.3 水源地开发利用现状及规划

截止2015年，呼和浩特市内在用的水井共有9个公共服务业用水水源厂下辖的119眼抽水井以及29眼城市供水补压井，均为承压水抽水井，总供水能力30.8×104m3/d。本文将研究区内主要抽水井进行分类描述见图4及表1。由于呼钢水厂和城市供水补压井年抽水量较小，因此本文数值模拟部分只考虑其他110眼抽水井。

近年来，因部分开采井成井质量低，止水效果差，或未采取止水措施，以及存在一些未处理过的废弃旧井，造成承压和潜水含水层之间连通，开采中存在混采问题，且因为地下水的超采，呼市上部潜水含水层水位较下部深层含水层高，模拟时需要注意受污染的上层地下水越流污染承压水的问题。

图4 研究区内抽水井分布及承压含水层边界条件Fig.4 Distribution of the pumping wells and boundary conditions of the confined aquifer

名称抽水井数量(眼)实际单井平均供水量/(m3·d-1)一水厂161 288二水厂212 862三水厂142 911四水厂172 750五水厂141 930六水厂91 966金川水厂11686如意白塔水厂8103呼钢水厂9调节运行,实际供水量较低补压井29调节运行,实际供水量较低

2 数值模拟法保护区划分

2.1 水文地质概念模型

首先对呼和浩特市含水层的实际空间结构、渗透性质、补给排泄及边界性质等条件进行概化，以便于进行数值模拟。

呼和浩特市城区148口抽水井位于呼和浩特平原，该处北面和东面环山，其余地方地势平坦，含水层结构大致分为下部承压含水层，中部厚层淤泥质黏土层以及上部潜水含水层。该地区主要供水水源为下部承压水，因此本次模拟主要含水层对象为下部承压含水层。

模型横向长41.0 km，纵向长34.6 km，分为潜水含水层及承压含水层两层，按照500 m×600 m网格剖分。因为北部单一含水层与南部双层结构的潜水和承压水都有水力联系，为它们提供补给来源。按照地层时代，将北部单一结构含水层分为两层，上更新统和全新统为上层，中更新统和下更新统为下层。其中下层与南部双层结构中的承压含水层合并为深层含水层，上层与南部双层结构中的潜水含水层合并为潜水含水层。因北面潜水含水层部分疏干，数值模型假定淤泥层的储水能力及水平渗透系数可以忽略不计，将淤泥层概化进潜水含水层，水文地质参数即概化为两者的等效值，由不同的垂向水力传导系数反映弱透水层，即采用MODFLOW中的“准三维”方法(图2)。

模型潜水含水层顶板设置为潜水含水层的自由水面，高程为1 018～1 080 m；潜水含水层底板结合北部单一结构含水层厚度，根据淤泥层底板高度确定，高程为827～1 112 m；承压含水层底板根据第四系中更新统下段含水层厚度及地面标高确定，结合单一结构区基岩底板，高程为503～1 061 m。

依据2011年5月水位统测资料绘制模型初始地下水流场如图4所示，根据流场特征和对含水层结构的分析，将呼市工作区的北部、东部确定为定流量边界条件，为流入边界，模型的第一、二层地下水经过此边界接受山区侧向补给；将呼市工作区南西部边界确定为定水头边界，为流出边界，模型的第一、二层地下水经由此边界流出工作区，见图4。

2.2 数学模型

工作区水位和地下水补、径、排特征综合反映地下水流为非稳定流特性，水文地质参数体现出含水层的非均质特征，本文使用可考虑非稳定、非均质地下水流系统及越流效应的“准三维”数学模型[9～10]：

潜水为：

(1)

承压水为：

(2)

式中：h1——潜水含水层水位标高/m；

b——潜水含水层厚度/m；

h2——承压含水层水位标高/m；

z(x，y)——隔水底板标高/m；

M＇——弱透水层厚度/m；

Kx、Ky、Kn——x、y及边界面法线方向上的渗透系数/(m·d-1)；

Tx、Ty、Tn——x、y及边界面法线方向上的导水系数/(m2·d-1)；

K——淤泥层的垂向渗透系数/(m·d-1)；

μ——潜水含水层的重力给水度；

S——承压含水层储水系数；

ε——含水层的源汇项/(m·d-1)；

p——蒸发和降水入渗强度/(m·d-1)；

Ω——渗流区域；

Γ0——渗流区域的上边界；

Γ1——渗流区域的二类边界；

q(x，y，t)——二类边界单宽流量/(m3·d-1)，隔水边界为0，流入为正。

2.3 数值模拟

因作者掌握的承压含水层静态水头资料为2011年5月31日测量数据，且2015年动态水头资料最为完善可靠，因此本次数值模拟模型的模拟期为2011年5月31日—2015年12月31日，共55个月。将整个模拟期按照对应的自然月划分为55个应力期，所有外部源汇项的强度在每个应力期中不变。模型的验证期为2015年1月—12月，共12个月，以一个月作为一个应力期。采用MODFLOW中的PCG算法进行计算[11]，计算得到的年均水均衡状况，见表2，研究区内潜水含水层向下部承压含水层的越流总量约为0.21×108m3/a。

表2年均地下水均衡表

Table2Annualgroundwaterbudget

补排项补排量/(108 m3·a-1)降雨补给0.263补给项侧向补给0.906小计1.169侧向排泄0.198排泄项地下水开采1.031小计1.229补排差-0.060

如图5所示，模拟的地下水流场基本能够反映地下水的流动规律和趋势。研究区总体水流方向为由东北、东南流往西南，地下水降落漏斗位于城区，范围与实际情况基本相符。在研究区西北部和南部，拟合相对差些。简而言之，研究区的模拟流场与实际情况总体一致，拟合较好。

图5 2015年12月31日承压水头计算值及水位观测井Fig.5 Calculation of pressure head and groundwater level observations wells on December 31, 2015

拟合的水位观测点分布情况见图5。大多数水位观测井观测值与模拟值拟合良好，两者之间相差不大，能较好地反映该点水位随时间变化趋势；研究区中部降落漏斗处和南部小部分观测点水位拟合相对较差，但随时间变化趋势一致。误差主要与获取的地质资料丰富程度、源汇项统计误差、以及控制性观测井分布情况等因素有关。此外，研究区的源汇项复杂，模拟范围大，且存在很多概化问题，因此误差不可避免。图6～图8列出了研究区的几口主要观测井水位拟合状况。

图6 OW-6观测井承压水位拟合效果图Fig.6 Fitting of groundwater levels at the OW-6 observation well

图7 OW-16观测井承压水位拟合效果图Fig.7 Fitting of groundwater levels at the OW-16 observation well

图8 OW-25观测井承压水位拟合效果图Fig.8 Fitting of groundwater levels at the OW-25 observation well

3 地下水保护区划分

根据《饮用水水源保护区划分技术规范》(HJ 338—2018)[12]，地下饮用水源一级保护区是以取水井为中心，井口周围溶质运移100 d的距离所圈定的范围；在一级保护区范围以外，井口周围粒子运移1 000 d的距离所圈定的范围为二级保护区，径流区和补给区为准保护区。

以各大自来水厂的抽水井为中心环形设置反向示踪粒子，利用Visual MODFLOW中的MODPATH功能计算出各口抽水井周围示踪粒子逆向运移100 d和1 000 d的流线，流线的最外层质点连接起来得到的范围分别就是一、二级保护区，见图9。

3.1 一级保护区划分

数值模拟法计算出的各单井一级保护区范围均不同，难以罗列，但计算得出的一级保护区形状均类似于图10所示的喇叭形，因此采用Visual MODFLOW中MODPATH的计算，各孔隙承压水水源井数值模拟法一级保护区范围计算结果见表3。

3.2 二级保护区划分

根据数值模拟法得到的示踪粒子1 000 d流线，结合现实道路、河流等划分确定了二级保护区范围见图11。

3.3 保护措施

针对呼和浩特市地下饮用水源保护区，提出相应的保护建议：

(1)对于一级保护区，取水点周围设立警示标志，区域内不得堆放垃圾、废渣，严禁大量使用农药、化肥，并对区域内存在的混采井进行安全封堵，一切工程施工需进行污染相关论证，消除一切可能导致地下水污染的因素。

图9 示踪粒子运移100 d和1 000 d模拟图Fig.9 Traces of particle movement for 100 d and 1 000 d

图10 二水厂水源井区划示意图及一级保护区范围Fig.10 Sampling of the wellfields in the second waterworks and the scope of the primary protection areas

名称地下水流向长轴/m短轴/m侧轴/m一水厂-5082 32 52 二水厂-172180 22 85 三水厂102105 23 55 四水厂-72247 25 63 五水厂19132 22 43 六水厂-151115 27 53 金川水厂-28124 20 55 如意白塔水厂-170138 21 32

(2)对于二级保护区，需保证在污染物流入一级保护区之前有充分的时间发现并消除污染，以确保水源地的安全。区域内合理开采利用地下水，停止所有混采井的使用，并建立水位水质长期监测机制，定期监测水位、水质。

(3)需严密论证潜水含水层和承压含水层的连通性问题并排查一切连通两层含水层的工程施工，严防承压含水层受到污染。