潮白河再生水补给河道对周边浅层地下水影响的数值模拟研究

姜瑞雪，韩冬梅，宋献方，李炳华

（1.中国科学院地理科学与资源研究所陆地水循环及地表过程重点实验室，北京 100101；2.中国科学院大学资源与环境学院，北京 100049；3.北京市水科学技术研究院，北京 100048）

水资源短缺已经成为制约我国经济社会可持续发展的主要瓶颈[1]。源自废水处理的再生水越来越多地作为替代水源被广泛利用，主要用于农业灌溉、景观用水、工业再利用、市政环卫用水、地下水补给等，很大程度上缓解了水资源匮乏的压力[2-4]。北京自20世纪80年代开始开展污水再利用工程[5]，目前已是我国再生水利用量最大的城市。2004年北京雨水和再生水供水2.04×108m3，占总供水量的6%[6]。2020年再生水利用量已达12×108m3，占总供水量的29.6%；其中11×108m3的再生水被用于河湖补水，占再生水利用量的92.3%[7]。

已有研究表明，再生水补给河道可促进周边地下水水位上升、储量增加[8]。经过河床沉积物和含水层介质的过滤，再生水中氮磷、总有机碳、部分重金属和部分有机物可以一定程度地被去除[3]。但长期使用再生水会造成受水区浅层地下水（埋深30 m）水化学特征逐渐与再生水特征相似[9-11]，Cl-、Na+、K+含量增加、总硬度降低[12]，部分抗生素和内分泌干扰物渗入[13-14]。潮白河流域的再生水利用工程实施已有十余年，但缺乏再生水对周边地下水的影响范围与程度的定量认识。

本文针对潮白河再生水补给的河道开展研究，运用MODFLOW 和MT3DMS 程序建立再生水补给河道周边的三维地下水流和溶质运移模型，基于10年（2007—2017）的地下水监测数据，模拟分析了浅层地下水水位及储量的变化，估算了伴随再生水补给河道，受水区地下水中氯离子（Cl-）和硝态氮（NO3-N）的负荷变化，明确了再生水补给河道对周边浅层地下水的影响范围和程度。

1 研究区概况

研究区位于北京市东北部的顺义区（图1），距离北京城区约30 km。顺义区年平均气温11.5 °C，年均降水量610 mm，降水主要集中在6—9月，多年平均水面蒸发量为1 247 mm。

为缓解潮白河供水短缺问题，恢复部分河段生态环境，北京市“引温济潮”水资源利用工程于2007年底建成通水[15]，截至2017年累计引再生水2.3×108m3入河[16]。2007—2011年底，工程一期的常年受水河道为减河和潮白河的土坝—河南村橡胶坝河段，向阳闸—土坝段每年5月和10月间歇受水。2012年工程二期实施，常年受水河道增加了河南村橡胶坝—苏庄橡胶坝河段（图1）。研究区内的潮白河为自然河道。减河为人工型河道，两侧砖砌边坡，河底为自然淤泥。

图1 研究区位置和地下水采样点分布示意图Fig.1 Location of the study area and sampling sites

研究区处于潮白河冲洪积扇的中下段，基岩之上沉积了大量的潮白河冲洪积物，表面堆积物主要是砂和亚砂土。沉积厚度自东北向西南变厚，东北向西南方向从200 m 增厚到千米以上[17]。地层主要为第四系孔隙含水层，根据南北向水文地质剖面图（图2）[9]，河道下部含水系统主要由4 个不同深度含水层构成，两两含水层间为一粉质黏土弱透水层。30 m 深和50 m深含水层自北向南从以卵砾石为主，渐变为以细砂为主。80 m 深和80～200 m 深含水层以细砂为主，局部夹有砾石。自北向南沉积物颗粒由粗变细，含水层厚度减小，层数变多，含水层间的水力联系逐渐变差。补水前区内地下水流向为由西南、南、东南向东北、北、西北流入地下水漏斗区域。

图2 研究区水文地质剖面图Fig.2 Hydrogeologic profile along P-P’

顺义是北京重要的水源地，研究区河道周边的水源地有4 处（图1）。地下水开采层位主要在埋深80 m以下。地下水除了水平方向的运动外，还存在自上而下的垂向运动。

2 地下水水流和溶质运移模型

2.1 概念模型

根据研究区水文地质条件，将地层概化为7 层，4 层含水层自上而下分别为30 m 深潜水含水层、50 m深承压含水层、80 m 深承压含水层和80～200 m 深承压含水层；两两含水层之间为一弱透水层，共3 层。模型的底部边界最深为200 m。

研究区地下水的主要补给来源包括降水入渗、再生水入渗和侧向流入。再生水补给河道后形成了约5.63 km2水面，因此将河流入渗处理为面状补给。排泄方式包括地下水开采、侧向流出和潜水面蒸发。结合模拟区的水文地质条件、地下水流系统形状信息，以及区内外的地下水水位情况，将北侧设为零流量边界，西侧和东侧设为第三类边界，南侧取第一类边界，见图3（a）。模拟面积为255 km2。潜水含水层自由水面为模型的上边界。

图3 （a）模型边界及（b）水文地质参数分区Fig.3 Boundary conditions, precipitation infiltration zones and values (a) and hydrogeological parameter zones (numbered) in the model (b)

在水流模型的基础上建立溶质运移模型，选择稳定保守的Cl-作为溶质运移因子。河道周边地下水中Cl-含量升高是再生水补给后地下水水化学变化的显著特征[10-11]，再生水补给前后Cl-浓度的变化可以反映再生水补给河道对周边地下水的影响范围和程度。根据对再生水出水和河水水质的长期监测，将河道设为第一类浓度边界，分段设置Cl-补给浓度，浓度总体为60～120 mg/L 之间。以2007年12月为初始时间，利用地下水水位和Cl-浓度的监测数据，通过插值赋值地下水的初始水位和Cl-初始浓度。

2.2 数学模型

通过对研究区水文地质条件的分析，将研究区含水层系统概化为非均质各向同性、三维非稳定的地下水水流和溶质运移模型。时间上，以2007年12月为起始时间，2027年12月为结束时间，时间步长为1 个月。2007年12月—2017年12月为模拟期，2017年12月—2027年12月为预测期。空间上，采用等间距有限差分方法离散地下水数值模型，平面剖分的网格尺寸为100 m×100 m，7 层共178 500 个网格。

再生水的入渗补给量利用水量平衡进行估算：

式中：Ql—再生水入渗补给量/（m3·a-1）；

Qrw—再生水引水量/（m3·a-1）；

Qpr—降水补给河道的水量/（m3·a-1）；

Qrc—河道蓄水量/（m3·a-1）；

Qe—水面蒸发量/（m3·a-1）；

Qt—向上游调出的水量/（m3·a-1）。

根据研究区的地表岩性，将模拟区分为三个区域，入渗系数见图3（a）。降水入渗量为：

式中：Qp—降水入渗补给量/m3；

P—降水量/mm，使用数据为密云站的月降水量[18]；

α—降水入渗系数；

F—补给区面积/km2。

水源地深层地下水的集中开采为研究区地下水主要的排泄方式，本研究着重考虑图1 中标出的4 处水源地的开采。参考2008年的地下水开采情况，开采总量为7.79×107m3/a。2014年南水进京后各大水源地减采，研究区的地下水开采量减少为4.26×107m3/a。华北平原区的潜水埋深大于3 m 时，潜水蒸发作用消失[19]，模型中的地下水蒸发极限深度取3 m。

参考前人对该地区的研究，将各含水层和弱透水层细化分区，见图3（b），对不同层不同区的水文地质参数进行赋值[20]。地下水流模拟中使用的参数包括渗透系数（K）和储水系数（S）。在Cl-运移的模拟中，使用的参数包括孔隙度（n）和弥散度（αL、αT）。

2.3 模型的识别与验证

浅层地下水（即埋深30 m 的潜水）受再生水补给的影响最为显著[9,21]，使用2007年12月至2015年12月浅层地下水水位和Cl-浓度的监测数据对模型进行参数识别，使用2016年1月到2017年12月的监测数据进行验证。水流模型率定的参数包括再生水入渗补给量、K、S。溶质运移模型率定的参数包括n、αL、αT。校准后的参数取值范围见表1。

表1 校准后的模型参数Table 1 Hydrogeological parameters used in the calibrated model

使用监测值与模拟值的均方根误差（RMSE）和决定系数（R2）评估识别的结果。图4 展示了浅层地下水水位和Cl-浓度监测值和模拟值的拟合及验证情况，监测值与模拟值变化趋势一致、误差较小。根据监测数据绘制了2017年12月浅层地下水水位等值线图，与模拟结果对比具有较好的一致性（图5）。结果显示，建立的地下水模型可以在一定程度上反映研究区真实的地下水流运动和Cl-运移情况。

图4 浅层地下水水位和Cl-浓度监测值与模拟值拟合及验证曲线Fig.4 Simulated versus observed results for groundwater tables and Cl- concentration in the shallow groundwater

图5 2017年12月模拟和监测（a）浅层地下水水位等值线图和（b）氯离子浓度等值线图Fig.5 Contour map of the observed and simulated groundwater tables and Cl- concentration in December, 2017

2.4 模型结果的处理

污染负荷表示环境对污染物的负载量[22]，本文利用模型中离散的网格，根据每个单元内的地下水量和溶质浓度计算污染负荷：

式中：P—污染负荷/t；

Ci—第i单元地下水的溶质浓度/（mg·L-1）；

Vi—第i单元地下水的水量/m3，从模型中计算得到。

Cl-浓度从模型结果中获取，NO3-N 浓度通过将地下水中NO3-N 浓度的监测值在模型离散的网格中插值获取。

2.5 情景设计

为了预测不同补给排泄条件影响下，未来浅层地下水水位和水质的变化情况，通过保持其他条件不变，仅分别改变降水补给量、地下水开采量、再生水补给情况，设计了预测情景，具体见表2。

表2 模型预测情景设计Table 2 Scenario design for the groundwater model

首先设置了一个基准情景，将其作为其他情景结果的参照。设置预测期内（2018—2027）的月降水量重现模拟期（2008—2017）的月降水过程。地下水开采量和再生水补给河道水量均维持2017年的现状不变，地下水开采量为4.26×107m3/a，再生水引水量为2.89×107m3/a。

根据密云站1956—2020年的月降水数据[18]，计算出不同保证率下的典型降水年份和降水量，枯水年、平水年、丰水年的年降水量分别为454.8，574.6，710.9 mm。情景1、2、3 分别设置为预测期内所有年份的月降水量重现典型降水年的月降水过程。

为探究地下水开采量变化的影响，情景4 和情景5 分别设置为2018年开始引潮入城水源地开采量增加25%和减少25%。

2019年研究区河道周边的浅层地下水中再生水的比例已高达95%[23]，为了探究再生水的补给对周边地下水的影响范围和程度，设置情境6 为2018年开始不再使用再生水补给河道。

3 结果与分析

3.1 再生水补给河道对浅层地下水水位的影响

3.1.1 补给后浅层地下水水位变化

受深层地下水开采影响，向阳闸以北是研究区浅层地下水水位最低的区域。再生水补给河道后，向阳闸以南水位整体上有上升的趋势。由图4（b）可知，间歇受水区的浅层地下水水位（G01、G15）在补给的前2年（2007—2009年）上升了约4 m，之后在波动中下降。第7—10年（2014—2017年），基本保持稳定，仅存在季节性波动。土坝以南的常年受水区，在补给的第1年（2007—2008年）水位快速上升，G22 上升了约3.5 m，G30 上升了约2 m。第2—4年（2008—2010年），G22 水位下降约1 m 后保持稳定；G30 水位基本稳定。二期工程实施后（2012—2017年），G22 和G30 的水位均呈缓慢上升趋势。

3.1.2 浅层地下水均衡和水量变化

根据模拟结果，2007—2017年潜水含水层的水均衡项见表3，总补给量小于总排泄量，潜水含水层处于负均衡状态。降水入渗补给量占总补给量的84.6%，其次是侧向流入和再生水入渗。潜水向下层含水层的越流约占总排泄量的97.4%，大于总补给量。

表3 2007年底至2017年底潜水含水层水均衡统计Table 3 Groundwater budget in the unconfined aquifer from December, 2007 to Deccember, 2017

研究区潜水位＞第一承压水位＞第二承压水位＞第三承压水位，垂向水力梯度诱使地下水从浅层向深层越流。结合水均衡分析认为，正是深层地下水开采引起的潜水越流，导致了潜水层水量的减少。

对潜水含水层的水量进行计算（图6）。再生水补给的前4年（2007—2011年），水量逐年降低，从4.8×108m3下降到了3.9×108m3。尽管再生水补给后河道周边的浅层地下水水位上升，但向阳闸以北受深层地下水开采影响，水位仍然处于下降状态。2012年底再生水补给河道的二期工程实施，受水河道范围增加，再生水的入渗补给量增大，浅层地下水水位上升的面积增加，水量回升至4.1×108m3；之后继续下降，2014年底下降到3.76×108m3。2014年开始向阳闸以北的水源地开采量减半，深层含水层超采情况得到缓解，潜水向深层含水层的越流量从1.04×108m3/a（2008—2013年）下降到了8.81×107m3/a（2014—2017年）。加之再生水的持续入渗补给，浅层地下水的水量开始缓慢上升，到2017年底水量回升至3.85×108m3。

图6 浅层地下水的水量变化及Cl-和NO3-N 负荷变化Fig.6 Variations in the shallow groundwater storage, Cl- loads in the model area, and Cl- loads and NO3-N loads in the affected zones

3.1.3 对浅层地下水水位的影响范围和程度

在前述地下水模型中去除再生水入渗补给项，其他不变，运行模型得到的结果即为无再生水补给情境下研究区地下水水位的分布情况。比较同一时间有再生水补给和无再生水补给的模拟结果，两者之间的差异即是受再生水补给影响而发生的水位变化。

如图7 所示，模拟期的第1年（2008年底）有再生水补给和无再生水补给情景模拟水位的差值为0～10.5 m，第2年（2009年底）0～12.5 m，第5年（2012年底）0～16 m，第10年（2007年底）0～19 m。随着再生水的持续补给，由再生水引起的水位上升范围逐渐扩大，水位上升的幅度也逐年增大。说明尽管再生水入渗补给在地下水补给项中并非占比最高，但正是再生水的持续补给，抬升了河道周边的地下水水位，维持水位的稳定。

图7 有再生水补给和无再生水补给条件下的浅层地下水水位差Fig.7 Difference of the groundwater table with and without recharge of reclaimed water

3.2 再生水补给河道对浅层地下水水质的影响

3.2.1 补给后浅层地下水Cl-浓度变化

由于再生水中Cl-浓度高于地下水本底值（5～75 mg/L），随着再生水的入渗，河道周边地下水中Cl-浓度也随之上升。如图4（d）所示，间歇受水区浅层地下水的Cl-浓度在补给后的前2年上升最为显著，G01、G15 的Cl-浓度从40 mg/L 左右上升到约90 mg/L。而常年受水区G22 的Cl-浓度在补给后的第1年就从20 mg/L 左右上升到约100 mg/L。之后河道周边地下水的Cl-浓度始终保持在升高后的浓度水平，但仍符合生活饮用水卫生标准（GB 5749-2006）（≤250 mg/L）[24]。

3.2.2 对浅层地下水水质的影响范围

河道周边的浅层地下水主要来自入渗的再生水和降水的混合，再生水是地下水中最主要的高浓度Cl-来源[8,23]，因此可将再生水补给后浅层地下水中Cl-浓度上升的区域视为受到再生水补给影响的区域。以每年12月底的Cl-浓度为参考，绘制浅层地下水水质受到再生水影响的范围图（图8）。

图8 浅层地下水质受再生水影响的范围变化Fig.8 Zones affected by reclaimed water for the shallow groundwater after reclaimed water restoring the river channel

可以看出，水质影响区分布于河道两侧，左岸的影响范围大于右岸。水质影响范围随着时间和补给河道的增加而逐渐扩大，从2008年底的11.7 km2逐年上升到了2017年底的26.7 km2。2009年和2012年的影响范围面积增速较快，这是由于再生水补给后的前2年，Cl-浓度上升速度最快；而2012年再生水二期工程实施，受水河道增加，水质影响范围也随之增大。河道受水范围稳定后，周边地下水中Cl-浓度逐渐稳定。2013—2017年底的水质影响范围虽逐年增大，但年增速逐年减小，说明在补给水量水质稳定的情况下，影响范围也逐渐趋于稳定。

3.2.3 补给后浅层地下水Cl-和NO3-N 负荷的变化及原因分析

通过计算浅层地下水中的Cl-负荷可以看出（图6），水质影响区内的Cl-负荷逐年升高，从2008年底的1.8×103t 上升到了2017年底的3.8×103t。再生水的补给抬升河道周边地下水水位的同时，也向地下水中输入了较高浓度的Cl-，使得河道周边的Cl-负荷升高。但相反的，模拟区的Cl-负荷从2007年底1.7×104t 逐年下降到了2017年底6.5×103t。这一方面是因为模拟区内浅层地下水向深层排泄，水量减少；另一方面是由于再生水入渗导致的Cl-浓度升高只发生在河道周边较小的范围内（即水质影响区内），而影响区外的Cl-浓度在降水和侧向径流的稀释中降低。

NO3-N 是再生水中的典型污染物[25-27]，在水质影响区范围内计算NO3-N 的负荷，结果显示（图6），除2011年外，水质影响区内的NO3-N 负荷整体呈逐年下降状态，从2008年底的29.8 t 下降到了2017年的11.9 t，其中2011年的负荷为48.8 t。表明再生水补给后周边浅层地下水的NO3-N 负荷整体下降，但在某些时刻存在高值，仍需关注。

氮素的转化主要受反硝化作用和硝化作用的影响。再生水中NO3-N 浓度峰值为20.2 mg/L，而浅层地下水NO3-N 浓度峰值为14.7 mg/L，多数时间为接近检测限水平（0.03 mg/L）。NO3-N 高值多出现在土坝以北区域。总体上，再生水入渗穿过河流-含水层界面易于形成反硝化热点，NO3-N 在地层中发生反硝化反应，被极大程度的去除[28-29]。地下水中的溶解氧、总有机碳和溶解性有机碳低于河水，也反映了反硝化作用的发生[11,30]。在对地下水的监测中也观察到，再生水补给后，浅层地下水中的NO3-N 浓度降低[10]。加之土坝以南地层的黏土层较厚，更有利于氮素的去除，G22监测井附近NO3-N 的衰减率高达99.6%[26]。因此土坝以南浅层地下水中NO3-N 浓度降低，NO3-N 负荷也降低。而土坝以北区域，砾石层厚、渗透性好，易于保持氧化环境，对NO3-N 的去除能力相对较弱[17]。且土坝以北是间歇性受水河道，仅在5月、10月补水，淹水期和干燥期交替出现，NH4-N 在淹水期易被沉积物吸附，在干燥期发生硝化反应转化为NO3-N[31]。当反应条件合适时，土坝以南间歇性受水河道附近更易出现高浓度的NO3-N（≤14.7 mg/L），导致某些时刻NO3-N负荷的增加。

3.3 模型预测

其他补排项不变，在不同降水量情景下，浅层地下水水位丰水年＞基线情景＞平水年＞枯水年，见图9（a）（b）。基准情景的浅层地下水水位的年际变化最小，近似为平衡状态。图9（b）可见，2027年底丰水年情景下的浅层地下水水位分别比平水年和枯水年情景高1.8，3.3 m。

图9（c）（d）展示了不同地下水开采量情景下的浅层地下水水位。结果显示，减少地下水开采可以恢复浅层地下水水位，离开采井越近受影响越明显。2027年底减采情景下的G01 水位相比增采情景高2.1 m，而在G22 监测井这一差距几乎可忽略。

保持降水量和地下水开采量与基准情景相同，仅改变再生水补给条件的情景下，图9（e）（f）展示了无再生水补给后周边浅层地下水水位的变化。结果显示，与基准情景相比，无再生水补给后周边浅层地下水水位快速下降，逐渐达到相对稳定的状态。2027年底G01 和G22 监测井基准情景和无再生水补给的水位差分别达到了2.5，4.5 m。进一步说明了正是再生水的持续补给维持了浅层地下水水位的稳定。

图9 不同情境下的浅层地下水水位变化预测结果Fig.9 Changes of the groundwater tables under different scenarios

浅层地下水中Cl-浓度的预测结果在不同降水情景下（情景1～3）的差异小于10 mg/L。不同的地下水开采量情景下（情景4 和5），Cl-浓度的预测结果几乎没有差异。由于再生水是模型中高浓度Cl-的唯一来源，无再生水补给的情景下（情景6），浅层地下水中的Cl-浓度迅速被降水和侧向径流稀释降低。

4 结论

（1）再生水补给后的2年（2007—2009年），周边浅层地下水水位上升约3～4 m，之后保持稳定。深层地下水的开采导致向阳闸以北区域浅层地下水水位持续下降，2007年底到2017年底浅层地下水水量处于负均衡状态。

（2）地下水的压采措施和再生水的持续补给可促进浅层地下水水位和水量的恢复。引潮入城水源地减采25%可使得周边浅层地下水水位在10年中（2018—2027年）抬升约1 m。若停止使用再生水补给河道，周边的浅层地下水水位将在10年中（2018—2027年）下降约5 m。

（3）受再生水入渗的影响，河道周边浅层地下水中的Cl-浓度在补给后的2年从补给前的5～75 mg/L增加到了50～130 mg/L，之后保持稳定。浅层地下水水质受再生水影响的范围从2008年底的11.7 km2增长至2017年底的26.7 km2。水质影响区内的Cl-负荷逐年增加，而NO3-N 负荷总体呈下降趋势，是由于NO3-N 在入渗过程中不断衰减。

利用再生水补给河道可以有效抬升河道周边浅层地下水的水位，增加水量，维持水位水量稳定；对水质的影响仅存在于在河道周边有限范围内，但长期使用再生水可能带来的污染隐患不容忽视，尤其是新型污染物（如内分泌干扰物、微塑料等）。未来应继续对再生水补给河道周边的地下水进行长期监测，关注地下水水质的动态变化。建议使用数值模拟方法优化再生水补给方案，增强再生水利用地区地下水的保护和管理，实现水资源的可持续利用。