截渗墙渗漏对黄水河下游地区地下水水质的影响研究❋

刘贯群，张淑琪，李 熊

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

地下水资源是中国许多城市生活、生产、生态用水的重要水源，对当地经济社会发展起着重要作用。滨海地区是人口高度集中、经济快速发展的地区，因对地下淡水的需求量过高导致过度开采，破坏了咸水与内陆淡水的动态平衡，最终出现海水侵入淡水含水层的现象[1]。目前，全球已有几十个国家和地区出现海水入侵现象，导致地下水水质恶化，引发了一系列生态环境问题，严重制约经济社会的可持续发展[2-4]。黄水河流域是龙口市工农业生产和居民生活的重要供水水源地，是该市经济发展的基础保障。20世纪70年代后期黄水河流域出现海水入侵现象，到20世纪80年代海水入侵现象加重[5]，为了控制黄水河区域海水倒灌，1995年在黄水河下游距海岸线 1.2 km处修建一座地下水库截渗墙，截渗墙具有拦截入海地下潜流和阻挡海水入侵的双重作用。截渗墙的修建有效阻止了截渗墙两侧的水力联系，同时也使内陆地区残留了一定区域的咸水体。

了解滨海含水层海水入侵的影响因素对于地下水资源的合理开发利用具有重要意义，多位学者研究发现海平面上升、气候变化、地下水开采对海水入侵程度产生影响[6-9]。针对截渗墙这一影响因素，国内外多位学者研究了截渗墙类型、高度、位置等因素变化对海水入侵的影响[10-12]，少有学者考虑截渗墙墙体出现渗漏的可能性，而水利工程施工工序较多，施工人员、外部环境及材料等因素将导致在施工过程中或建成后出现渗漏的情况[13]，海水仍可侵入淡水含水层，故有必要研究其对地下水水质的影响程度。本文运用SEAWAT模块对黄水河下游地区建立三维变密度地下水数值模拟模型，针对不同截渗墙渗漏情况设计多种预测方案，分析未来研究区地下水氯离子浓度的演变趋势。

1 研究区概况

龙口市是山东省沿海地区重要开放城市，位于胶东半岛的西北部，渤海湾南畔。龙口市南部地势较高，北部地势平坦，其地表水主要集中在南部山区，地下水主要分布在北部平原区。黄水河为龙口市境内最大的河流，发源于栖霞县北部山区，从东南向西北流经龙口市境内流入渤海，干流全长55.43 km，流域面积 1 034.47 km2。研究区域位于黄水河下游，地处龙口市东北部的平原区，地理坐标为120°30′20″E—120°38′16″E, 37°37′16″N—37°44′36″N，总面积约85 km2，具体位置如图1所示。

图1 研究区地理位置示意图

研究区属于暖温带半湿润季风型大陆性气候[14]，四季分明，年平均气温11.8 ℃，年平均蒸发量1 150～1 250 mm。根据研究区内的诸由观水文站1966—2020年降水量资料显示，年平均降水量为557 mm，降水多集中在6—9月份，约占全年总降水量的72%，7—8月份降水量最多(见图2)。

图2 诸由观水文站降水量图

2 海水入侵数值模拟研究

2.1 水文地质概念模型

黄水河下游地区主要含水层为第四系含水层，含水层厚度在20～55 m之间。潜水含水层可大致分为两层，上层含水层岩性以砾质粗砂为主，下层含水层主要由粗砂、中砂组成，局部为细砂，上层含水层渗透性比下层含水层强。研究区水文地质剖面如图3所示，剖面位置参见图1。依据含水层岩性、埋藏条件及开采条件，将研究区潜水含水层在垂直方向上划分为两层，概化为非均质各向同性的三维潜水含水层。

图3 水文地质剖面图

研究区含水层在东西两侧逐渐尖灭，东、西两侧视为隔水边界。北部边界视为定水头边界，采用渤海常年海平面平均值作为定水头边界的值，取值为0 m。研究区地下水位自南向北呈递减趋势，南部边界受到上游地下水补给，地下水侧向流入研究区域，因此将南部边界视为给定流量边界。对于溶质运移模型，东、西部边界设置为零质量通量边界。北部边界定为定浓度边界，浓度值采用研究区实测海水样氯离子值，取值为19 000 mg/L。南部边界定义为定浓度边界条件。研究区域的上边界为潜水含水层自由水面，含水层的下边界为透水性差的粘土层，组成区域性隔水底板。

2.2 数学模型

变密度过渡带模型考虑密度对水头、流速和浓度的影响，能够有效刻画复杂水文地质条件、人类活动等因素影响下的咸淡水界面运移规律[15-17]，是当前海水入侵数值模拟的重要研究方向。过渡带海水入侵的数学模型中须用两个偏微分方程来描述, 方程(1)描述了密度不断改变的混合液体(咸淡水混合物)的流动。

(1)

式中：hf为等效淡水水头值(m)；Kfx、Kfy、Kfz为沿x、y、z方向的渗透系数(m/d)；Z为计算点高程(m)；t为时间(d)；ρ为实际液体密度(mg·L-1)；ρf为等效淡水的密度(mg·L-1)；Sf为等效淡水单位贮水率(1/m)；θ为有效孔隙度；C为溶解物质的浓度(mg·L-1)；ρs为源、汇项中溶解物质的密度(mg·L-1)。qs为单位时间进入单位体积含水层源、汇项的体积(d-1);Ω为计算区；h为(x,y,z)在t时刻的水头值(m)；h0为点(x,y,z)处的初始水位；q(x,y,z)为第二类边界上单位面积补给量，q=0代表隔水边界,为第二类边界条件。

方程(2)描述了混合液体中盐分的运移。在考虑溶质运移过程，即盐分运移过程时，选用化学性质上较为稳定的氯离子作为模拟组分。在海水入侵过程中,氯离子的运移方式以对流和弥散为主，暂不考虑运移过程中可能进行的化学反应。

(2)

2.3 变密度地下水数值模型

2.3.1 网格剖分与时间离散 利用Visual MODFLOW中的SEAWAT模块建立研究区变密度地下水数值模拟模型，SEAWAT模块广泛应用于例如海水入侵这样的变密度流问题，其模拟结果具有较高的准确性与可信度[18]。根据模型区域内地层、流场、边界条件等特征，以及考虑到模型计算的精度和耗时长短，恰当地设计有限差分单元网格，经综合考虑，平面上将研究区剖分为70行、65列边长为200 m的正方形网格结构，垂向上将模型划分为2层。

整个模拟期包括两个阶段:2019年4月18日—10月31日为模型识别期，2019年11月1日—2020年9月1日为模型验证期。选择研究区内24眼地下水观测井的实地观测资料进行模型识别与验证，以研究区2019年 4 月 18 日的实测数据作为初始地下水流场和初始浓度场。初始地下水位等值线图及初始氯离子浓度等值线图如图4所示。

图4 初始地下水位等值线图和氯离子等值线图

2.3.2 水文地质参数及源汇项 水文地质参数关系到地下水模型的准确性，是进行精确科学计算地下水资源和地下水均衡的基础。水文地质参数分区及初值根据研究区水文地质条件和前人抽水试验资料确定[19-20]，针对研究成果较少的水文地质参数，则根据研究区不同岩性的经验值进行确定[21]。最终将研究区划分11个水文地质分区，具体水文地质参数分区见图5。

图5 含水层水文地质参数分区图

研究区潜水含水层主要受大气降水补给和侧向补给,近几十年研究区降水量较少，黄水河上游王屋水库无下泄流量，导致河水基本被水库截留，黄水河一直处于断流状态，不与地下水进行水量交换，故在模型建立时暂不考虑河流作用。根据研究区内黄河营、诸由观、兰高3个雨量站的降水量数据,采用泰森多边形法进行大气降水补给分区(见图6)，再根据降水入渗系数方法计算降水入渗补给量，3个雨量站的降水量数据均来自烟台水文信息网的实际统测数据。南部边界以侧向补给的形式给研究区提供补给量，采用达西公式计算其侧向径流补给量。

图6 降水分区图

研究区地下水的主要排泄方式包括城市生活、工业、农业等人工开采；研究区地下水埋深4～18 m，均大于3 m，因此不考虑蒸发排泄。根据烟台市水资源公报和统计年鉴，确定识别期和验证期的地下水开采总量，农业开采量按照山东水科院提供的农业开采月比例系数分配到每月，工业开采和生活开采量按月平均分配。

2.4 模型识别与验证

采用试估-校正法对模型进行校正，反复调整水文地质参数以达到较为理想的拟合效果，识别验证后水文地质参数值见表1。图7(a)和(b)分别为地下水位和氯离子浓度值的拟合情况，识别与验证期水位拟合绝对误差小于0.5 m的分别达到82%和78%，均高于70%；识别与验证期氯离子浓度绝对误差小于15%所占比分别为90%和87%，均高于50%，满足精度要求，说明模型已达到较为理想的拟合效果。图7(c)和(d)分别为研究区从北到南三眼观测井识别期地下水位和氯离子浓度计算值与实测值的拟合情况，图7(e)和(f)为验证期拟合情况，可以看出模型计算结果与实测值的趋势一致，说明通过识别和验证，得到了能够比较真实反映研究区地下水流场和氯离子运移实际情况的数值模型。

图7 模型识别期与验证期拟合情况

表1 各区水文地质参数值

3 截渗墙渗漏对地下水的影响研究

3.1 预测方案构建

对1966—2020年黄水河流域降水量进行降水频率分析，得出研究区实际降水量变化规律为丰-枯-平-枯-枯；丰水年(P=25%)降水量为710.7 mm，平水年(P=50%)降水量为517.5 mm，枯水年(P=75%)降水量为400 mm。以2020年9月1日为基准年，预测未来5年地下水氯离子浓度变化趋势,根据降水量年内变化规律，将一年划分为4个应力期(9—10月、11—4月、5—6月、7—8月)，预测期共有20个应力期。

在丰-枯-平-枯-枯水年条件下，探究地下水库截渗墙不同渗漏位置对研究区地下水水质的影响，具体预测方案如表2所示。首先，考虑截渗墙不同方位渗漏对研究区地下水的影响，将截渗墙边界按地理位置和水文地质条件，划分西部、中部和东部分区，并在3个分区内各选取1处渗漏点，每一处裂缝宽度取为1 m，对渗漏位置所在的单元格进行加密处理，将边长为200 m的单元格剖分为宽度为1 m的单元格；截渗墙用Visual MODFLOW中的水平阻滞(HFB)程序包来模拟，考虑渗漏时移除渗漏位置处的墙体。其次，选取截渗墙西部渗漏的情况，分析截渗墙位于上层含水层和下层含水层的部分分别渗漏时对研究区地下水水质的影响，渗漏点的宽度为1 m且上部和下部渗漏点的平面位置相同。

表2 海水入侵预测方案

3.2 预测结果分析

以氯离子高于250 mg/L的面积作为海水入侵总面积，表3显示了各预测方案下研究区海水入侵程度，S0代表现状年(2020年9月1日)海水入侵程度。从表3和图8中可以看出，在丰-枯-平-枯-枯水年条件下，考虑地下水库截渗墙渗漏的情况时，研究区海水入侵范围扩大。截渗墙西部和中部出现渗漏时的海水入侵范围增加幅度相近，与现状年相比分别增加9.52%和9.11%。当截渗墙东部位置的墙体出现渗漏时，对研究区地下水的影响最大，氯离子浓度高于250 mg/L面积高达23.38 km2，与现状年相比增加了3.49 km2，增长了17.58%。这是由于现状年时研究区东北部地下水为淡水(<250 mg/L)，当海水侵入时，地下水逐渐咸化，导致研究区高于250 mg/L的面积相比现状年增加幅度最大。

表3 不同预测方案下海水入侵程度对比

图8 模型运行5年后截渗墙不同渗漏位置氯离子浓度变化图

14号井靠近截渗墙的西部地区，距离西部、中部、东部渗漏点分别约190、3 000和5 800 m。图9(a)为14号地下水井不同渗漏位置氯离子浓度变化图，可以看出14号井氯离子浓度值受截渗墙西部渗漏的影响最大，受截渗墙中部和东部渗漏的影响较小，预计2025年9月1日的Cl-浓度比中部、东部渗漏时的Cl-浓度高1 100 mg/L左右。17号井靠近截渗墙中部地区，距离西部、中部、东部渗漏点分别约2 430、1 080和2 800 m。

图9 不同位置单独渗漏时各地下水井Cl-浓度变化图

从图9(b)中可以看出17号井受截渗墙中部渗漏的影响较大，截渗墙西部渗漏和东部渗漏对其影响相近，预计5年后中部渗漏时Cl-浓度比西部、东部渗漏时大约高80 mg/L。靠近截渗墙东部的3号井距离西部、中部、东部渗漏点分别为5 150、1 350和3 080 m，模型运行5年后，当东部出现渗漏时，氯离子浓度与西部和中部渗漏相比仅增加30 mg/L左右(见图9(c))。因此，可以得出当截渗墙出现渗漏时，距离截渗墙越近的地下水井氯离子浓度受其渗漏的影响越大，且靠近一处渗漏位置的地下水井，受其他两处渗漏时的影响较小。

表3和图10反映了截渗墙上部和下部分别渗漏时对研究区地下水水质的影响，可以看出，模型运行5年后，与截渗墙下部渗漏相比，当上部出现渗漏时海水入侵总面积更大，地下水咸化程度更严重；上部渗漏时 14号井的氯离子浓度比下部渗漏高300 mg/L左右，这是由于上部含水层渗透性较强，当截渗墙出现渗漏，海水涌入时，向内陆入侵的速度更快，导致内陆地区地下水咸化程度更重。

图10 截渗墙位于不同含水层位置渗漏时14号井中Cl-浓度变化图

4 结论

基于SEAWAT模块建立变密度地下水流模型和溶质运移模型，利用研究区域水文地质资料及地下水动态监测数据，通过反复调整水文地质参数，对模型进行了识别和验证，得到了可以较好的反映研究区地下水流场和氯离子运移过程的地下水数值模型。为了探究地下水库截渗墙渗漏对研究区地下水的影响，在丰-枯-平-枯-枯水年，针对截渗墙不同位置出现渗漏设计多种预测方案，预测结果如下：

(1)当截渗墙不同位置出现渗漏时，研究区地下水氯离子浓度与现状年相比均升高。以地下水氯离子浓度高于250 mg/L的面积作为海水入侵总面积，西部渗漏和中部渗漏时，海水入侵面积增速相近。当截渗墙东部出现渗漏时，对研究区地下水的影响最大，海水入侵面积与现状年相比增长了17.58%。

(2)当截渗墙出现渗漏时，距离截渗墙越近的地下水井受其渗漏的影响越大，且靠近一处渗漏位置的地下水井受其他两处渗漏时的影响较小。

(3)针对截渗墙位于不同含水层部分的渗漏情况进行预测分析，得出截渗墙上部出现渗漏时地下水咸化程度更严重。

(4)截渗墙渗漏将会加重黄水河下游地区海水入侵程度，影响研究区地下水资源的开发利用，因此应高度重视黄水河地下水库截渗墙渗漏的探测工作，可适当采取防渗加固措施，以确保截渗墙发挥良好的截渗效果。此外，应减少研究区地下水开采量，抬高截渗墙南侧区域的地下水位，缩小截渗墙南北两侧的水力梯度，以减少截渗墙出现渗漏时所造成的影响。模拟及预测结果可为未来地下水管理提供理论依据，对经济社会的可持续发展具有重要的现实意义。