截渗工程条件下滨海地下水源地地下水补给源解析

黄俞铭， 阎 妮，2，3， 郑西来，2，3， 刘贯群，2，3

(1. 中国海洋大学环境科学与工程学院， 山东 青岛 266100；2. 中国海洋大学山东省海洋环境地质工程重点实验室， 山东 青岛 266100；3. 中国海洋大学海洋环境与生态教育部重点实验室， 山东 青岛 266100)

地下水是宝贵的淡水资源，为全世界大多数的城市提供生活、农业和工业用水。在沿海地区，地下水需求量的增加引起地下水过量开采，从而导致地面沉降、海水入侵等环境地质问题，海水入侵引起的地下水水质恶化已经成为了全球性的水文地质问题[1]。为了防治海水入侵，在滨海地区修建了各类防治海水入侵的工程措施。滨海含水层本身就具有地层结构复杂、水利条件脆弱等特点[2]，在防治工程措辞的影响下，滨海含水层变成了一个更为复杂的地质体，地下水来源解析也因此变得十分困难。

国内外研究者已经对滨海含水层地下水来源展开了多方面的研究。地下水来源解析方法主要分为水化学法和模型计算法。水化学方法包括统计分析法[3]、离子比值法[4-5]、Piper图解法[6]和吉布斯图解法[7]等。水化学法已经被广泛应用于不同海水入侵区地下水来源分析的研究，但是水化学法作为一种定性分析方法，无法定量解析出地下水的来源，特别是化学成分相似的补给源，仅计算1～2个保守物种来确定混合比例。同时研究者发现，对于同一地下水样品，当选择不同的保守组分进行混合比计算时，计算出的混合比可能会有很大的不同[8-9]。

基于这种情况，为了对地下水来源进行定量解析，研究者建立了多种结合多元统计分析、水化学分析和数学计算在内的模型计算方法。最小二乘法首先被广泛应用于地下水来源解析的研究[8]，其主要目的是为了克服分析误差，将样品点投影到不同端元的混合线上。

但是最小二乘法也仅仅使用了2～3种指标进行来源解析，缺少对样品的整体认识。在最小二乘法的基础上，Hooper等[10]结合主成分分析(PCA)提出了端元混合分析(End-Member Mixing Analysis, EMMA)，该方法综合选取地下水样品中的保守和非保守组分同时作为参数，期望用最少的端元数量来解释特定区域的地下水来源，该方法被广泛应用于地下水、地表水等不同水体以及沉积物、大气污染物等研究领域[11]。但是，端元混合分析法不能计算混合比，因此，需要结合其他方法。Carrera等[12]基于最大似然法提出了混合计算(MIX calculations)程序，该方法应用最大似然法计算计算混合比和端元组成的方法。该方法承认端元组分的不确定性，通过非线性优化考虑所有混合物，选择合适的选择示踪剂来计算潜在的补给源和混合比例。验证表明，相比水化学法和其他模型计算法，混合计算程序的计算结果更接近真实情况。但是模型计算法简化了地下水流动过程和其中可能存在的水文地球化学反应，假定了地下水流动的连续性，并在此基础上建立溶质和水的质量平衡计算，可能无法准确反映研究区复杂的水文地质条件。因此，对特定研究区进行长期水文地质监测，为模型计算提供依据是地下水源解析过程中必不可少的步骤。

基于以上分析，本文以黄水河滨海地下水源地为研究对象，通过现场的地下水水质的系统监测，在水化学解析方法的基础上，采用端元混合分析确定地下水潜在的补给来源，利用混合计算程序计算不同地下水来源的混合比，确定研究区地下水补给来源和混合比，为该区域水资源开发利用和地下水污染治理提供科学依据和技术支持。

1 研究区的环境背景

1.1 自然地理概况

研究区位于中国山东省烟台市龙口市东北部，坐标北纬37°37′—37°43′，东经120°30′—120°37′(见图1)。北部地区的海拔不到1 m，南部地区达到25 m,属温带季风气候区。根据1978—2018年的气象记录资料，研究区年平均降雨量大约680.5 mm，多年平均蒸发强度为1 932 mm。研究区内降雨在年内呈现出分配不均的特点，每年6—9月为汛期，降水量占全年降水量的70%～77%。年平均气温约11.2～12.4 ℃，1月份平均气温最低，约-4.3—-7.0 ℃，历史最高低温为-17.9 ℃；7月份平均气温最高，约25.5～27.8 ℃，历史最高气温为38.9 ℃[13]。

图1 研究区位置与采样点分布图Fig.1 Location of study area and distribution of sampling points

黄水河为研究区的主要河流，发源于栖霞市，经龙口境内流入渤海，年平均流量约4.77 m3/s，平均径流量约1.51 m3/s，流域面积1 066 km2。为充分利用黄水河径流，在最上游修建了王屋水库，王屋水库为大(二)型水库，总库容1.21亿m3。同时在河流中下游修建了拦河闸和地下坝等工程，组成了黄水河下游地下水库，总库容5 329万m3，同时具有储存地下水与防止海水入侵的作用[14]。

1.2 水文地质条件概况

1.2.1 地下水赋存条件 研究区含水层为第四纪孔隙潜水含水层，厚度在0～70 m之间，含水层为冲洪积细砂、砾质粗砂、中粗砂和卵砾石[12]。研究区含水层大致分为四个部分：① alQ4砾质粗砂，为高渗透层，渗透系数为237.39 m/d；② mQ4中粗砂；③ pl-alQ3粗砂，分布于黄水河Ⅱ级阶地，渗透系数约为126.5 m/d；④ pl-alQ2中粗砂，局部夹杂卵砾石，该层的渗透系数为28.09～130.26 m/d(见图2)[14]。

图2 研究区I-I’水文地质剖面图Fig.2 I-I’ hydrogeological profile of the study area

1.2.2 地下水的补给、径流和排泄 研究区地下水补给来源主要为大气降水、地下水侧向补给和地表水，其中最主要的补给来源为大气降水和地表水。地下水排泄途径主要是人为开采，由于地下水埋深均大于3 m，因此蒸发较弱。在天然条件下，地下水面形态与地形基本一致，即南高北低。地下水水力坡度较小，为1～3‰，水流缓慢。

1.2.3 地下水动态变化规律 研究区地下水受人类活动和降雨量影响较大，因此在年内呈现出明显的季节变化规律。在春季，地下水水位显著下降，5—6月达到地下水水位最低值；6—9月进入汛期，地下水水位开始上升，8—9月达到地下水水位最高值。地下水水位的年际变化呈现出明显的分配不均的现象，在枯水年或连续枯水年，地下水水位大幅度持续下降；丰水年或连续丰水年地下水水位又持续上升。

1.2.4 地下水水化学特征

1.2.4.1地下水样品采集 分别在2019年4和9月以及2020年9月对黄水河地下水源地进行3次采样，每次采集样品21个，包括19个地下水样品、1个上游地表水样品和1个海水样品，并在2020年8月采集大气降水样品1个。采样点选择的主要依据包括研究区水文地质条件、土地利用类型、海水入侵程度等。分别选取滨海区和农业活动区内具有代表性的采样点。采样井主要为工业用井、农业用井和已发表研究成果[15]的监测井。样品收集方法以及室内水质监测方法参照《中华人民共和国地质矿产行业标准 地下水质检验方法》(DZ/T 0064.1～0064.90—93)。对得到的离子组分进行电荷平衡误差(Charge balance errors, C.B.E.)计算，保证所有样品的电荷平衡误差小于5%。在现场用便携式YSI型多参数水质测定仪原位测定地下水的温度、TDS、pH、EC、DO和ORP指标。

1.2.4.2 地下水水化学特征 表1列出了研究区所有地下水样品点的物理指标和化学组成情况，包括每个指标的最大值、最小值、标准偏差和平均值。研究区多数地下水样品没有表现出明显受到海水入侵的地下咸水的离子组成特征，水样总体的TDS值较低，且γ(Na+)的平均值低于γ(Ca2+)，与1994年[15]相比，大部分地下水样品中Cl-浓度显著降低。接近海边的地下水样品(1、2、3、4和19号)中Cl-浓度较高，体现出了受海水入侵影响的特征。因此根据采样点地理位置、水化学特征和土地利用类型，将研究区分为滨海区和农业活动区，其中滨海区包括的地下水采样点为1、2、3、4和19号，其余采样点位于农业活动区内(见图1)。

表1 研究区地下水水化学分析结果统计Table 1 The statistical table of groundwater hydrochemical monitoring in the study area

2 研究区地下水来源解析

2.1 水化学分析方法

2.1.1 主要阳离子与氯离子的比值 在受到海水入侵的地下水中，通常TDS主要受到海水的影响，因此主要阳离子与氯离子的比值通常接近海水[16]。主要阳离子与氯离子的比值与对应样品的氯离子浓度相结合，可以作为判断地下水中离子的海洋来源的重要依据。图3为研究区不同水体主要阳离子与氯离子的关系图。在Na+与Cl-的比值图以及K+与Cl-的比值图中，低Cl-浓度的样品中更接近海水中的比值，这与通常的海水入侵的研究结果不同。低Cl-浓度的地下水的Na+和K+与Cl-的比值接近海水的主要原因可能是岩石的溶滤风化作用导致岩盐溶解；而滨海区地下水中Na+与Cl-的比值低与海水的原因是截渗工程使海水入侵程度降低，滨海区地下水由咸水转变为微咸水，Na+和Cl-的浓度降低，同时发生了Na+与Ca2+、Mg2+之间的阳离子交换作用。

图3 研究区不同地下水中阳离子与氯离子关系图Fig.3 Relationships between main cations vs. chloride of groundwater samples in the study area

Ca2+与Cl-的比值以及Mg2+与Cl-的比值随着Cl-浓度的增加而降低，但是都远高于海水中的比值。说明在Cl-浓度较高的地下水中，Cl-与Ca2+、Mg2+的比例更接近海水，但是均高于海水。说明海水入侵一定程度上影响了地下水离子组成，但是研究区地下水总体表现出低TDS的特征。同时滨海区地下水中Mg2+与Cl-的比值比Ca2+与Cl-的比值更接近海水中的比值，主要是由于Ca2+为非保守离子，在地下水形成过程中，参与了多种矿物的沉淀——溶解过程和阳离子交换反应。

Ca2+与Cl-的比值以及Mg2+与Cl-的比值随着Cl-浓度的增加而降低，但是都远高于海水中的比值。说明在Cl-浓度较高的地下水中，Cl-与Ca2+、Mg2+的比例更接近海水，但是均高于海水。说明海水入侵一定程度上影响了地下水离子组成，但是研究区地下水总体仍表现出低TDS的特征。同时滨海区地下水中Mg2+与Cl-的比值比Ca2+与Cl-的比值更接近海水中的比值，主要是由于Ca2+为非保守离子，在地下水形成过程中，参与了多种矿物的沉淀-溶解过程和阳离子交换反应。

图4 研究区不同地下水吉布斯图Fig.4 Gibbs maps of groundwater samples in study area

2.1.3 优化Piper图解法 为了识别海水入侵的影响，Gopinath等[3]使用7种不同的分类标准，将传统的Piper图法进行了优化。图5为使用改良Piper三线图对研究区地下水化学类型进行分析。结果表明，滨海区的地下水样品均位于海水入侵水区域，农业活动区的地下水样品主要分布在轻微海水入侵水区域，但是也有Cl-浓度较低的地下水样品落在海水入侵水区，这些地下水样品中Na+和Cl-离子的毫克当量百分比较高，可能是受到岩盐溶解的影响，而非海水入侵造成的。

图5 研究区地下水样品优化Piper图Fig.5 Optimized Piper diagram of groundwater samples in the study area

2.2 耦合模型计算法

2.2.1 模型的原理和主要计算方法 本研究针对黄水河滨海地下水源地，使用Tubau等[11]提出的一种结合使用端元混合分析和混合计算来计算混合比的方法对研究区进行地下水来源解析。该方法包括4个主要步骤：①研究区水文地质条件调查与水质分析；②选择要在分析中使用的指标；③ 确定潜在的地下水补给源；④使用混合计算程序计算混合比。其中步骤②和步骤③需要在端元混合分析中进行，步骤④在混合计算程序中进行。

2.2.1.1 端元混合分析(EMMA) 端元混合分析提供了一种利用水化学数据的特征向量和残差分析来估计端元的数量的方法，它的目标是确定解释该研究区地下水补给来源的所需的端元的最小数目。在端元混合分析过程中，主要遵守的原则和步骤为：①对所有样品进行主成分分析，该步骤主要是为了消除特征值最大的指标；②对于分析后特征向量的选取，遵循“1规则”[10,19]，即保留所有特征值大于或等于1的特征向量，数量为n；③最终选取的端元的数量为n+1；④最后将每个样品点的特征值投影到特征向量所定义的空间中。根据投影图对端元进行初步选择，并对落在投影定义区域之外的点进行分析。

2.2.1.2 混合计算(MIX calculations)程序 混合计算是Carrera等[12]提出的一种对地质年代学线性拟合方法的扩展，其主要步骤是：①初始化，即假设端元浓度已知，定义初始混合比；②通过极大似然函数重新评估端元和地下水样品的组分；③使用地下水样品组分的浓度的期望值，利用极大似然函数获得最终混合比；④重复步骤②和步骤③直到收敛。针对端元组分和混合比构建的极大似然函数分别为：

(1)

式中：yp是第p个样品中所有指标的向量；Ap是它们的协方差矩阵；F是选择的端元的指标的ns×ne维矩阵；δp是混合比组成的向量。极大似然函数的具体求解方法可以参考Carrera等[12]。

图6 EMMA特征向量组分相对贡献系数图Fig.6 Relative contribution coefficient of eigenvector components of EMMA

根据端元混合分析的特征向量投影图(见图7)，选取了截渗墙外咸水(19号)、大气降水(R1)以及农业灌溉地表水(11号)作为端元。

图7 研究区地下水样品特征向量投影图Fig.7 The projection of the eigenvector of groundwater samples in the study area

2.2.3 不同地下水补给来源的混合比

2.2.3.1 混合计算输入数据 端元混合分析得到的3个端元 (19号、R1和11号)分别定义为:EM1、EM2、EM3。将2019年9月采集的样品中选定组分的浓度作为端元浓度进行混合计算，由于混合计算承认端元浓度的不确定性，计算对每个选定的端元都赋予了标准偏差，标准偏差根据3次采样的结果进行计算得到。端元的初始浓度及输入标准差见表2。

表2 输入端元初始浓度及标准差Table 2 Initial parameters of MIX program

2.2.3.2 混合比计算结果 3个端元对每组地下水的贡献如表3所示。对于研究区整体而言，EM3对研究区的贡献最大，混合比平均值为51.17%；EM1和EM2的混合比分别为10.22%和38.60%。

表3 研究区地下水混合比Table 3 Mixing ratios of groundwater samples in the study area %

截渗墙下游的地下咸水对滨海区地下水的影响十分明显，平均混合比为34.03%。截渗墙下游地下咸水对3号样品点影响最为明显，混合比为59.40%。大气降水和灌溉地表水对滨海区地下水的平均混合比分别为29.70%和36.27%。表明滨海地区地下水虽然受到海水入侵的影响，但是海水入侵程度已经逐渐减弱，大气降水和地表水的贡献增加。

对于农业活动区地下水，贡献率最高的是灌溉地表水，平均混合比为55.75%。截渗墙下游的地下咸水对其的影响十分微弱，平均混合比仅为2.91%，其中混合比最高的是8号样品，混合比为11.30%。同时发现大气降水混合大于50%的样品点，除5号样品外，全部分布在黄水河河道两侧，说明大气降水通过进入河道后入渗补给地下水。

3 结论

本研究利用水化学和数学计算方法来探究黄水河流域地下水源地特殊水文地质条件下的地下水补给源问题。主要结论为:

(1)端元混合分析与混合计算耦合的数学方法在判断特定研究区的地下水补给来源和计算每个补给来源的混合比方面，可以得到传统的源解析方法无法得到的定量计算结果，并具有更高的准确性，具有良好的应用价值。

(2)水化学初步判断可以得出，研究区地下水盐分形成过程中海水入侵、矿物溶解—沉淀和阳离子交换均发挥了不同程度的作用。海水入侵影响了滨海区地下水的离子组分；矿物溶解——沉淀主要是岩盐的溶解是农业区地下水中TDS的主要影响因素；阳离子交换则主要发生在滨海区，随着海水入侵过程的加剧，影响了地下水中Ca2+和Na+的浓度。

(3)根据采样点地理位置、水化学组分和土地利用类型，可以将研究区分为滨海区与农业活动区，其中截渗墙下游的地下咸水、大气降水和灌溉的地表水对滨海区地下水的平均混合比分别为34.03%、29.70%和36.27%，对农业活动区地下水的平均混合比分别为6.69%、49.03%和44.28%。对于研究区整体而言，3个补给源的平均混合比分别为10.22%、38.60%和51.17%。