南水北调中线受水区保定平原地下水质量演变预测研究

曹文庚，杨会峰，高媛媛，南 天，王 哲，徐素娟

（1.中国地质科学院 水文地质环境地质研究所，河北 石家庄 050061；2.自然资源部 地下水科学与工程科普基地，河北 石家庄 050061；3.水利部南水北调规划设计管理局，北京 100038；4.海河水利委员会水文局，天津 300170；5.河北省地矿局地质勘查技术中心，河北 石家庄 050081）

1 研究背景

南水北调中线受水区保定平原近40年来，由于长期、持续、大规模过度开采地下水，区域地下水水位持续下降，超采严重[1]，部分含水层被疏干或枯竭，引发了地面沉降、地裂缝、水质恶化等一系列生态与环境地质问题，严重制约当地经济社会可持续发展和生态文明建设，危及生态安全和供水安全。南水北调东中线一期工程经过10 多年的建设分别于2013年、2014年顺利通水。随着工程的通水，2018年南水北调工程开始对保定市的拒马河、白洋淀进行生态补水，地下水超采状况有所缓解，保定平原受水区供水格局已发生改变，在缓解水资源紧张局面的同时，水位恢复不当也会引起一系列的环境地质问题[2]，特别是由于水位上升，在水岩交换作用下，会引发地下水化学组份的变化，局部地区地下水存在潜在污染的可能。因此本文基于校正的地下水流模型，利用美国地调局开发的PHREEQC 软件，模拟2030年浅层地下水位上升后，水岩交换作用下保定平原地下水主要化学组分变化趋势，分析地下水质量演变机制。

目前，基于PHREEQC 的水岩交换模拟主要采用水文地球化学系统的化学平衡模拟方法[3]，国内外学者利用PHREEQC 开发出多种模型和方法进行水岩交换反应模拟，Claus 利用PHREEQC 对不同深度的地下水化学组份进行了建模，定义了一种一维的反应输运模型，该模型考虑了地下水位上升过程中部分氧化的含硫黄铁矿沉积物与地下水中所发生的水文地球化学反应模拟[4]；王广才团队建立了耦合的Monte Carlo-WATEQ4F模型对水-岩反应平衡状态进行模拟和分析[5]；Stigter 团队建立正向地球化学模型方法，模拟了地表水与地下水混合灌溉对地下水环境的影响机制[6]；汤洁团队[7]及赵江涛等[8]利用反向水文地球化学模拟方法解释了地下水化学成分的形成、分布及演化机理。通过梳理总结，目前国内外在水岩交换方面较为成熟和流行的模型包括：（1）表面吸附模型（SCMs），主要基于水相及固相的化学特性，用数学的方法刻画水岩交换作用中的离子吸附作用、竞争及协同效应，尤其擅长模拟矿物表面化学键，例如络合反应[9-10]；（2）扩散双层模型（DLL），是显示的计算扩散层组成和非静电层的表面络合的模型，可以模拟水溶液中选择性地吸附某种离子[11]。（3）电位分布-多点位络合模型（CD-MUSIC），可以描述铁氧化物及氢氧化物表面的吸附行为和离子形态分布[12-13]，还能预测环境条件改变下的离子形态的响应特征[14]，同时可以捕捉主要离子的非线性电荷效应[15]。因此，CD-MUSIC模型非常适合本次模拟的主题，故选择CD-MUSIC模型模拟地下水位恢复后的地下水化学组分的演化趋势。

2 研究区概况

2.1 自然地理概况南水北调受水区保定段主体为保定市平原区，北临京廊、西依太行、东临沧衡、南邻石家庄，行政区划上包括保定市区、满城区、清苑区、徐水区、雄安新区、涿州市、定州市、高阳县、蠡县、高碑店、定兴、望都和安国，平原区面积约1.1 万km2（图1）。该地区地势平坦，海拔多在50 m 以下，自山前向白洋淀微倾斜。地表水系属大清河流域，流经的主要河流有拒马河、易水河、漕河、界河和唐河等。

研究区地貌单元主要包括太行山山前冲洪积倾斜平原、冲积湖积微倾斜平原及平原东部的冲湖积低平原。山前冲洪积倾斜平原分布于易县、满城、唐县、望都、定州、曲阳和行唐一带，由拒马河、易水河、漕河、府河、唐河等河流冲洪积扇构成，地面高程80～30 m，地形坡降4‰～1‰，地势大致自西南向东北倾斜。冲积湖积微倾斜平原分布于定兴西南、容城以北、徐水、清苑、博野和安国一带，位于山前冲洪积倾斜平原东部。上部为近代河流冲积层或扇前洼地相堆积物，下伏冲洪积层。地势相对平坦，地面高程30～10 m，地形坡降2‰～1‰。冲湖积低平原分布于以白洋淀为中心的周边地带，主要为近代河流冲积相和湖沼相共同组成的平原，地势低平、地面高程10～7 m，地形坡降一般0.7‰～0.2‰。

2.2 水资源状况及环境地质问题据《保定市第二次水资源评价报告》显示，保定平原多年平均地表水资源量为1.49亿m3/a，地下水资源量为15.79亿m3/a，浅层地下水可开采资源量为15.39亿m3/a。依据2016年保定市水资源公报数据，2016年保定平原总供水量22.98亿m3，其中地下水开采量达到20.34亿m3，包括浅层地下水19.37亿m3，深层地下水0.97亿m3。保定水资源整体状况匮乏，自1980年代开始大规模开采地下水，多年的地下水超采造成地下水水位持续下降，形成多个浅层地下水位降落漏斗，依据2016年6月浅层地下水位等值线（图1（b）），形成了“保定徐水区-容城浅层地下水漏斗”、“雄县地下水漏斗”以及“高阳-蠡县-清苑浅层地下水漏斗”，最大水位埋深为49.02 m，漏斗总面积达到4144.4 km3。

2.3 南水北调保定平原区外调水量包括南水北调和引黄入冀工程，南水北调计划分配给保定水量为5.0015亿m3，引黄入冀工程计划分配保定水量2.55亿m3，白洋淀生态补水1.1亿m3，2017年开始试通水。2018年9月水利部选择河北省境内的滹沱河、滏阳河、南拒马河3条典型河流的重点河段进行生态补水，2020年计划将生态补水范围扩大到14条河流。南水北调持续来水，保定平原全区皆为南水北调受水区，一定程度上替代保定平原全区现有的地下水开采，大幅度压减了地下水开采量，从而引起区域地下水流场的变化，地下水位上升明显。

图1 保定平原水文地质钻探取样点及浅层地下水流场

3 研究方法

3.1 地下水流场预测获取2030年浅层地下水位基于已校正的地下水流模型[16]，依据国务院批复的《海河流域综合规划（2012—2030）》，南水北调河北省2020年分配水量30.4亿m3，保定市年均分配江水量5.05亿m3；2030年南水北调河北省规划分配水量分别为42.3亿m3，同时考虑雄安新区建设，保定2030年南水北调分配水量预计达到8.0亿m3。在预测模型中，模拟预测期为2016.6—2030.12，共175个应力期，以月为应力期。源汇项部分，白洋淀入渗、河流、降水等采用全部平水年的数据；边界流量采用现状模型的定流量处理。依据预测结果（图1（c）（d）），2030年保定平原全区浅层地下水位呈现不同程度的上升，上升幅度最大的为原地下水漏斗区域的雄县、高阳及蠡县，水位最大上升幅度达到35 m，山前地区水位上升幅度相对较小，一般小于15 m（图1（d））。

3.2 地下水化学组分模拟点位及钻探岩心采集深度选取在全区选取32个水文地质钻探点，点位的分布基本控制全区所有地质地貌单元（图1（a）），进行地下水样品及钻孔岩心样品采集，岩心的采集深度范围与2016 至2030年浅层地下水位上升高度范围一致，见图1（d），样品采集间隔为每10 cm 采集岩心样品50 mg，然后进行混合、过筛晾晒后送往河北省地质实验测试中心进行土样的易溶盐测试，测试项目包括：Na、Ca、Mg、HCO3、Cl、SO4含量。由于每个地下水样品的地下水环境不同，钻孔易溶盐含量也存在较大的差异，因此本次模拟针对每个地下水样品点都建立独立的水文地球化学模型，分别进行32组水岩相互作用的地球化学模拟。

3.3 水文地球化学模拟方法本次模型基于美国地调局开发的PHREEQC 软件，采用的水文地球化学模拟方法是建立在平衡常数基础上的平衡常数法。通常，水文地球化学模拟方法分为正向地球化学模拟和反向地球化学模拟，前者是依据假设反应及动力学参数来预测地下水化学组分及水-岩间的质量交换；后者是依据水、土样品的化学组分数据来反向推测水岩交换反应的过程。

（1）本次研究已经掌握的数据包括：①2016年浅层地下水化学数据；②2016年水样点附近岩心土样的化学元素含量数据。

（2）模拟预测的目标为：①2030年浅层地下水化学组分；②2030年浅层地下水化学类型。

本次模拟预测工作采用正向水文地球化学模拟，然而，正向模拟需要提供室内试验条件下的反应动力学方程以及相关模型参数。但是，由于本次工作缺乏室内试验，无法通过试验获取这些参数，因此在正向模拟之前，还需要加入反向地球化学模拟工作。依据2016年已知的地下水化学组分和岩心样品化学组分，通过反向模拟，获取了水文地球化学模拟的相关参数。然后，将参数重新代入正向地球化学模拟过程中，利用CD-MUSIC模型提供的反应方程，采用水文地球化学混合作用模型、表面络合作用模型、离子交换模型，进行正向地球化学模拟，预测得出2030年浅层地下水恢复条件下水化学组分（图2）。

图2 南水北调受水区保定段水位恢复地下水化学组分模拟技术路线

正向地球化学模拟主要基于以下4个方程式。

利用PHREEQC 进行地下水组分存在形式的计算，各种表面络合组分可有以下方程确定：

式中：Cj为第j中络合组分的活度；Kj为第j种络合组分的热力学平衡常数；αk为第k种组分的活度；Pk.j为第k种组分相对第j种络合组分的化学计量数。上式即质量作用方程，对于地下水中的主要组分，若第k种主要组分的化学分析浓度为Ak，则质量守恒方程为：

从物理意义上讲，ak和Cj必须满足如下条件：

正向模拟反应的基本步骤为：起点水溶液的化学组分+反应相→终点水溶液的化学组分+生成相。可用下式表示：

式中：p为矿物相（反应相+生成相）数目；ap为第p种矿物相进入或离开每升溶液的摩尔数；bp，k为第p种矿物相中第k种元素的化学计量数；mT.k为溶液中第k种元素的总摩尔浓度。

将32个水位点的化学组分数据以及相关热力学参数代入以上方程式进行迭代，分别建立32个水化学模拟模型，得出最终的评价结果。

3.4 模型参数选择基于本项目已掌握的数据、模拟目标及方法思路，结合保定平原浅层含水层的矿物环境特征及地下水的水化学特征，确定模型的三种参数，具体内容及选择理由为：①铁锰氧化物及氢氧化物的吸附参数，依据本区的1∶5万区域地质调查成果，此区域浅层地下水所赋存的晚更新世地层中，重矿物主要为针铁矿、磁铁矿、钛铁矿等，因此本次模型考虑针铁矿中铁锰氧化物的吸附性能，因此反应方程采用针铁矿和各离子吸附解析的平衡方程（见表1）。表1中，lgK为收集已有的化学反应的热力学常数；lgKopti为进行参数优化后的热力学常数；②表面吸附节点的参数，主要考虑表面吸附的物理吸附和微观结构特征，本次依据CD-MUSIC模型定义了针铁矿的两种吸附点位类型，第一种是只有一个铁原子的吸附点位，即Geo_uni，第二种为有三个铁原子的吸附点位，即Geo_tri，基于固相的表观形貌特征来反映吸附能力；③方解石、白云石、石膏等饱和指数。依据保定市1∶20万水文地质图及其说明书等资料，晚更新世是太行山保定段山前冲洪积扇发育最为强烈的时期，由于山区广泛分布的为蓟县系雾迷山组和铁岭组的白云岩、白云质灰岩等，受山区汇水区溶滤作用影响，平原区地下水阳离子主要为Ca、Mg 等离子，阴离子以HCO3及SO4为主，因此主要考虑方解石（CaCO3）及白云石（CaMg（CO3）2）以及石膏（CaSO4·2H2O）矿物的溶解及沉淀过程。

3.5 水文地球化学反应方程确定依据模型参数的选择，本次模拟的反应方程选取具体见表1。

3.6 模拟结果检验本次通过PHREEQC 模拟预测获取的2030年保定平原浅层地下水化学数据，需要对结果进行可靠性检查，从而保证水文地球化学模拟的质量。

本次采用阴阳离子平衡检查，天然地下水的化学成分中，阴阳离子是主要的组成部分，因此在水化学分析前进行阴阳离子的平衡检查是非常必要的。将所有阴阳离子的单位由原来的“mg/L”换算为毫克当量浓度（meq/L），通过计算水中阴阳离子的相对误差来判断预测的2030年水化学数据的可靠性。

离子平衡的检查公式为：

式中：E为相对误差，%；Mc及Ma为阳离子及阴离子的毫克当量数，meq/L。

4 结果与分析

4.1 反向水文地球化学模拟结果保定平原从1970年代开始大规模开采地下水，地下水位不断下降，本次选择某一个地下水样品点，例如以SH0538C 点为例，此点位于保定市博野县内，2016年水位标高为-9.45 m，预测2030年水位标高恢复至20 m，水位恢复接近30 m。本文通过钻探采集SH0538C 水样点附近岩心，选择水位恢复段岩心样品后，进行过筛、混合，对混合样品测试易溶盐组分，获取水位恢复段的岩土易溶盐组分特征，由于易溶盐是采用水溶法，即用4倍的纯水进行溶解，并测试溶解出来主要离子含量。所以，易溶盐组分特征实际上是岩土的水溶态。

表1 主要离子表面络合反应

基于以上条件，可以得到以下假设：水位恢复段的岩土水溶态含量基本等于自20世纪70年代水位下降过程中，表面络合到岩土颗粒上的吸附量。基于此假设，尝试建立反向地球化学模型，以岩土易溶盐含量作为模拟的终点，反向模拟强开采条件下水位下降过程中水岩交换反应，对研究方法中3个参数进行调整，包括：表面吸附点位参数、铁氧化物和氢氧化物的吸附参数以及方解石、石膏的饱和指数。

假设验证：通过调整参数，对Ca、Mg、Na、HCO3、SO4、Cl 共6种离子进行拟合，拟合度结果见表2。从表2可以看出，整体拟合度较好，Ca、Na、HCO3、SO4拟合度最优（图3及图4），Cl 拟合度最差，其原因是Cl 的化学性质较为稳定，在水岩交换过程中不发生任何交换反应，溶解度较高，因此造成Cl含量实测值与预测值拟合度较低，Cl 的拟合度基本可以满足本次模拟的精度。

表2 水位恢复段岩土易溶盐含量指标与反向水文地球化学模拟计算的表面吸附量拟合度评价指标

图3 易溶盐指标含量实测值与反向水文地球化学模拟计算的表面吸附量对比

4.2 正向水文地球化学预测结果基于反向地球化学模拟得出的表面吸附节点参数、铁锰氧化物吸附参数及方解石石膏的饱和指数，代入正向地球化学模型中，进行正向地球化学预测工作，基于方程（1）—（3），导入反向地球化学模型模拟得出表面络合参数、铁锰吸附参数以及矿物的饱和指数，在给定参数后就可计算出地下水中各种组分存在形式的浓度（活度）。建立水位恢复上升的一维流水-岩作用模型，基于PHREEQC 本身自带表面吸附模型和离子交换模型，对2016年每个地下水样品点进行正向模拟，分别输入各点水位上升前水化学数据，分别确定各自表面吸附过程、离子交换过程及矿物的溶解沉淀过程。最终获取保定平原全区2030年浅层地下水化学组分（表3）。

采用研究方法中3.6 的模拟检验公式，对预测出的2030年浅层地下水化学组分进行检验，32组水化学数据电荷平衡误差皆小于10%，本次预测结果较为可靠。

5 讨论

依据模拟预测结果，2030年南水北调受水区保定平原域浅层地下水位恢复明显，区内地下水化学组分也随之发生变化。为进一步分析南水北调供水对保定平原浅层地下水环境的影响，对保定平原2016年及2030年的地下水化学规律进行对比研究，分析水化学类型及影响地下水质量的主要指标的变化趋势。

5.1 水化学类型变化趋势依据2016年及2030年保定平原浅层地下水化学数据，分别制作了2016年浅层地下水化学类型分区图以及2030年浅层地下水化学类型分区图，对比二者可以看出：至2030年，南水北调供水条件下，由于水位大面积抬升，全区浅层地下水化学类型发生了较大的变化。

表3 2030年保定平原浅层地下水化学组分预测 （单位：mg/L）

总体上，浅层地下水逐渐向盐化方向发展，HCO3-SO4型水代替HCO3型水成为保定平原主要的阴离子水化学类型，平原东部的地下水排泄区出现了Cl型水；Na-Ca型水成为平原中部主要的阳离子水化学类型（图4）。

分带特征变化明显，由山前至平原中部，2016年浅层地下水化学类型由山前的Ca-Mg-HCO3型水，至平原中部逐渐演化为Na-Ca-HCO3-Cl型水；2030年由山前至平原中部，山前冲洪积平原的Ca-Mg-HCO3型水面积大幅下降，由2016年6749 km2下降至3351 km2，平原中部Na-Mg-HCO3-SO4型水面积由2016年的921 km2扩大至4109 km2（表4及图4）。大部分地区，Na和Cl 是地下水化学组分中增长比例最大的离子，造成此问题的原因是由于水位恢复段岩心存在的大量Na+及Cl-，岩盐的溶解是造程平原东部大部分地区Cl元素含量上升的主要原因，高蠡清浅层地下水漏斗是Cl含量增幅最大的区域，经过预测，至2030年，在平原东部地下水漏斗中心，行政区位置属于高阳县东北部，形成了Na-Cl型水，面积达到121.99 km2（表4）。

图4 2016年及2030年水位恢复条件下保定平原浅层地下水水化学类型分区对比

表4 保定平原2016年及2030年浅层地下水水化学类型面积对比

5.2 TDS变化趋势分析前人多项成果已经表明，地下水超采以及地下水的快速下降，会造成地下水TDS含量的增加。依据本次模型，水位回升条件控制下，2030年全区浅层地下水TDS含量变化表现出明显的水文地质分带性。依据2016年及2030年浅层地下水TDS含量制作的TDS含量变差图（图5）可以看出，自山前向平原东部，TDS含量由增加逐渐变为减少，保定平原北部永定河冲洪积平原、平原中部的易水河冲洪积平原、平原南部的唐河-沙河冲洪积扇，都表现为TDS 上升的趋势，Cl及SO4含量上升是造成TDS 增大的主要原因；保定市区、安新、雄县地区TDS含量表现为下降趋势，其分布范围与白洋淀附近的高水位区基本一致，TDS含量下降最大处分别位于保定市区及白洋淀周边地区。

5.3 硬度变化趋势分析硬度是研究地下水水质的重要指标之一，也是区域上影响地下水质量分级的主控因素。虽然目前已有大量的文献对硬度与地下水位、地下水开采的关系进行了分析论证，通常认为水位下降、地下水超采是造成地下水硬度升高的主要原因。随着南水北调供水、地下水开采减少，水位逐渐恢复，包气带厚度减小，部分包气带转化为饱和带，地下水环境由氧化环境逐渐变为还原环境，地下水入渗途径缩短，大大减少了白云石等矿物溶解带来的Ca2+及Mg2+含量。依据2016年至2030年浅层地下水硬度变差图（图5）可以看出，全区大部分地区，浅层地下水硬度下降明显。至2030年，硬度下降幅度最大的区域主要为雄安新区及其周边地区，硬度下降幅度超过200 mg/L，雄安新区及其附近水质有逐渐好转的趋势；而清苑南部以及保定市北部徐水容城一带，浅层地下水硬度有小幅度上升。

图5 2016—2030年保定平原浅层地下水TDS及硬度含量变幅

6 结论

针南水北调受水区保定平原开展地下水位回升条件下的水质演化预测研究，基于水岩相互作用过程中的表面吸附作用、铁锰氧化物及氢氧化物的吸附作用、离子的竞争与协同吸附效应，以及不同矿物的溶解平衡，分析得出：

（1）总体上，南水北调受水区保定平原浅层地下水逐渐向盐化方向发展，HCO3-SO4型水代替HCO3型水成为保定平原主要的阴离子水化学类型，平原东部高阳-蠡县地区出现了Cl型水，Na-Ca型水成为保定平原中东部地区主要的阳离子水化学类型，岩盐的溶解是造成平原东部大部分地区Cl元素含量上升的主要原因，高阳-蠡县浅层地下水漏斗区是Cl含量增幅最大的区域。

（2）至2030年浅层地下水位整体回升后，保定平原的东西部地下水质量演化趋势存在较大的差异。西部的山前平原浅层地下水TDS含量增大明显，其中拒马河冲洪积扇、易水河冲洪积扇以及唐河-沙河冲洪积扇增幅达到100 mg/L 以上，水位上升后地层中岩盐及石膏的矿物溶解是造成TDS 上升的主要原因。中东部平原的保定市区、安新、雄县、高阳等地区TDS含量表现为下降趋势，其分布范围与浅层地下水降落漏斗基本一致；至2030年，全区浅层地下水硬度下降明显，硬度降幅最大的区域主要为地下水漏斗中心的雄县、高阳及蠡县地区，硬度下降幅度超过200 mg/L，阳离子交换作用是造成Ca、Mg含量下降的主要原因。南水北调供水及压采后，水位恢复对地下水漏斗区域水质改善具有一定的积极作用。

（3）地表水的混合主要发生在山前的几条补水河道中，这种混合作用对局部河段影响带地下水质量的影响是比较大，但是由于补水期较短、补水量较小，且未见长期补水规划，对区域地下水质量影响相对较小。因此，本次水文地球化学模拟工作主要考虑的是“南水北调来水后，本地的地下水压采产生的水位抬升”，未能考虑地人工回补的地表水和地下水的混合作用，这是本次工作的缺憾，笔者下一步工作会开展相关研究，在模型中加入地表水混合作用。