潍河下游典型卤水开采区海咸水入侵趋势研究

蒙永辉 王集宁 张丽霞 罗梅

摘要：为了掌握潍河下游卤水开采对海咸水入侵的影响，在地质、水文地质、地下水位和地下水水化学等实测资料的基础上，基于潍河下游典型卤水贮存、运动特征和开采现状，建立了潍河下游典型卤水开采区的地下水流系统数值模型，并对卤水开采区不同开采条件下海水入侵的趋势进行了预测。结果表明：在保持目前地下水和卤水资源开采状态的情况下，海咸水入侵锋线向卤水区推进，且卤水区地下水矿化度大幅度降低；若保持地下水开采、停止卤水资源开采，则卤水区地下水位逐渐恢复，卤水区地下水矿化度有降低现象；海咸水入侵锋线处于相对稳定状态，原因是其处于地下水降落漏斗中心处，南部低于锋线值和北部高于锋线值地下水均向漏斗区汇聚，使得海咸水入侵锋线处于相对稳定状态。

关键词：海咸水入侵；典型卤水开采区；数值模拟；潍河下游

中图分类号：P314.1 文献标志码：A doi：10.3969/j.issn.1000-1379.2018.01.015

地下卤水资源具有很高的经济价值，可用于军工、盐化工、电子、制药等领域[1]。

潍河下游位于莱州湾中部，具有得天独厚的地下卤水资源，是全国最大的海盐及盐化工生产基地，截至2006年年底，潍坊北部沿海地区共有卤水开发企业4000多家[2]。近年，卤水资源开采导致了一系列的生态环境问题，包括水位下降、咸淡水分界线向内陆迁移、海水入侵、海洋生态环境破坏、卤水质量浓度和有益组分含量降低等[3]，其中海水入侵的危害尤为突出，目前该地区海水入侵面积已超过4000km2，入侵陆地的速度达500～1000m/a[4]。日益加剧的海水人侵严重污染了地下水环境，使土壤次生盐渍化，对海岸区域的工农业生产和生态环境造成严重的影响。因此，研究卤水资源开采与海水入侵的关系、评估海水人侵的趋势，具有重要的理论和实际意义。国内外学者对此做了很多研究工作，较为常用的研究方法包括数值模拟、水文地球化学和同位素方法等。Oahman等[5]利用二维有限差分模型SEAWAT模拟了Gaza海水入侵状况，认为地下水开采使海水入侵状况恶化；Karahanoglu等[6]利用有限元方法模拟了土耳其Kocaeli-Darica采石场海岸含水层海水入侵状况，研究了海水入侵行为、地下水位变化规律和盐分浓度的关系；Mahlknecht等[7]利用地下水中主要组分、Sr和B同位素，评价了人类活动影响下干旱地区海岸含水层的海水入侵状况；Cimino等[8]利用地球物理和地球化学方法，评价了西西里岛北部Acquedolci海岸含水层海水入侵状况，并提出了咸淡水混合模型；我国学者林文盘等[9]探讨了莱州湾海水入侵灾害防治的进展，提出海水入侵的机理本质上是海岸带咸淡水界面在地下水位变化情况下水动力平衡被破坏；章斌等[10]根据秦皇岛洋戴河平原地下潜水和海水的水化学特征，运用数理统计和模糊数学方法研究了该地区潜水含水层的海水入侵程度；李国敏等[11]根据海水入侵研究现状，认为数值方法是现今模拟和求解海水入侵问题的最有力工具；段梅[12]综合分析硇洲岛的地质和水文地质条件，建立了海水入侵数值模型，对硇洲岛在不同開采条件下的海水入侵趋势进行了模拟预测。但是，上述研究仅局限于单一影响因素，对海平面水位、风暴潮等影响因素鲜有考虑。

本研究利用潍河下游典型卤水开采区地下水位、水质长期动态系列资料，以及海平面水位、风暴潮侵入次数和距离等数据，建立以地下水矿化度为判断因子，以地下水对流迁移为主导，并考虑分子扩散、吸附作用、密度变化因素的变密度三维地下水流模型和溶质运移模型，对未来海咸水入侵的发展趋势进行预测，旨在为研究区卤水开发利用和保护提供技术参考。

1 研究区概况

研究区位于莱州湾南岸滨海平原东部、潍河下游，面积约946km2。气候属暖温带半湿润大陆性气候区，四季分明，气候温和，多年平均气温为11.9℃，多年平均年降水量560mm，其中6-9月降水量占年降水量的55%～77%，年均水面蒸发量为1788mm。地貌类型由南至北可分为南部山前冲洪积平原、中部冲积海积平原和北部海积平原。地层主要为第四系更新统至全新统冲洪积、海积沉积层，岩性以粉细沙、中粗沙和砾沙为主。地下水主要赋存于第四系沙砾石层等含水介质中，黏性土作为相对隔水层，形成松散岩类孔隙式的多层含水结构，地下水类型主要为松散岩类孔隙水，垂向分为浅层含水岩组和深层含水岩组。浅层含水岩组主要指深度小于120m的浅层潜水、微承压水，含水层岩性从山前平原中粗沙、砾石变为滨海区的粉沙、细沙，分为浅层淡水、浅层微咸水和浅层咸水，其中：浅层淡水主要赋存于山前冲洪积扇、潍河及其古河道堆积形成的河谷、阶地含水层中，矿化度<1g/L，单井涌水量为500～3000m3/d，水化学类型为HCO3·Cl-Na；浅层微咸水分布于柳瞳以北，矿化度为1～3g/L，单井涌水量1000～3000m3 /d，水化学类型为HCO3·Cl-Na；浅层咸水分布于第四系海相地层的松散沉积物中，矿化度>3g/L，水化学类型为Cl-Na。深层含水岩组主要分布于海积、冲洪积海陆交互堆积平原，岩性以中沙和细沙为主，单井涌水量为1000～3000m3/d，矿化度为1～2g/L，水化学类型为HCO3·Cl-Na。

2 数据来源和研究方法

2.1 数据来源

研究所用数据来自山东省地质环境监测总站和潍坊市矿产资源管理中心，包括气象、水文、地形地貌、地质、水文地质、地下水位、地下水水化学和地下水开采量等数据。目前，研究区卤水资源年开采量为4800万m3。

2.2 水文地质概念模型

2.2.1 地下水系统边界条件的概化

研究区处于潍河下游，南部以潍河冲洪积扇顶部为界，北部到海岸线，东部以胶莱河为界，西部以潍河与白浪河冲积地层的交接处为界。南部、东部和西部边界均为流量边界；北部海岸线边界为给定水位边界，即按照渤海湾海平面的观测，并考虑海水潮位变化，结合实际观测的海水潮位变化给定水位边界。垂向上，由于很多区域浅部卤水含水层（潜水层）地下水位已低于含水层底板，含水层被疏干，因此将浅部卤水含水层与第一承压（微承压）卤水含水层概化为一层，下部各承压卤水含水层概化为一层，咸水层以下深层淡水概化为一层。按照地质剖面地层结构，从北部卤水区向南将整个含水层结构概化为3层含水层和2层相对隔水层（弱透水层），把地下水系统概化为非均质各向异性三维非稳定流。

2.2.2 模型初始值确定

（1）初始流场和初始浓度场。将2015年12月的地下水流场作为模型的初始流场（水位等值线），见图1。把地下水矿化度作为海咸水入侵的鉴别因子，浅部地下水初始浓度场：南部以2015年年底地下水矿化度实测值为初始浓度值；北部卤水区根据2007年卤水勘察的卤水波美度值，依据卤水波美度与浓度的转化关系（表1）和研究区卤水主要阴阳离子比例分布，将卤水波美度转化为浓度，近似作为卤水矿化度值，并将其作为第一含水层和第二含水层的矿化度初值，形成第一、第二含水层模拟的初始浓度场（图2）。依据南部深层水界线上的矿化度情况，设定深层含水层（第三含水层）地下水矿化度均为2000mg/L。

（2）源汇项处理。源汇项主要为大气降水入渗、潜水蒸发和地下水人工开采。2016年研究区地下水资源开采量为5 177.5万m3，卤水资源开采量为4800万m3，地下水和卤水资源开采区均位于研究区北部的漏斗区，在模型中均以面状开采形式体现。

2.3 数学模型

根据上述概化的水文地质概念模型，建立研究区地下水概念模型和地下水溶质运移模型，其数学模型（控制方程）和定解条件为式中：H为地下水位（相对于淡水）；Kij为渗透系数张量（i，j=1，2，3）；η为密度耦合系数，η=ε/Cs，ε为密度差率，ε=（ρs-ρ0）/ρ0，Cs为与流体最大密度ρs对应的浓度，ρ0为参考密度（淡水密度）；ρ为混合溶液的密度；ρ*为现状条件下混合液密度；C为溶液浓度；ej为重力方向单位矢量第j个分量；AR为贮水率；φ为孔隙率；q为单位体积孔隙介质源（或汇）流量；xi、xj为笛卡儿坐标（i，j=1，2，3）；H0为地下水初始水位；HB为边界Г1上的给定地下水位；qB2为边界Г2上补排强度；t为时间；W'为潜水面边界Г3上的补排强度；H*为潜水面边界Г3上各点地下水位；He为潜水面边界Г3上各点的高程；ni为边界上外法线单位矢量；D为水动力弥散系数张量；μi为地下水实际流速在xi方向的分量；C*为抽出或注入液体的浓度；C0为初始浓度；CB为边界Г1上的浓度；C'为降水入渗的浓度。

2.4 模型识别及结果

利用2015—2016年的地下水动态资料对模型进行调试与识别。把2015年年底实测的地下水流场作为初始流场，对典型观测孔实测的地下水位动态曲线和2016年年底的地下水流场、地下水浓度（矿化度）场进行拟合。图3、图4、图5分别为长观孔地下水位、地下水流场和地下水浓度场拟合结果。经模型识别及检验，将全区分为5个降雨入渗系数分区，将含水层和弱透水层分别分为12、8个水文地质参数分区，见表2～表4和图6～图8。由于V区为晒盐池，且已作防渗处理，因此该区大气降水入渗系数为0。考虑到第四系含水层中夹有黏性土以及垂向上的水流沉积作用，模型中第四系水平渗透系数取垂向渗透系数的10倍。

从典型观测孔模拟水位和实测水位的拟合曲线可看出，水位模拟误差较小，所建模型能够反映研究区地下水流场的动态变化，可以预测不同开采条件下地下水的动态变化。

2.5 预测情景

海咸水入侵预报设置了两种情景：一是保持目前地下水开采和卤水开采状态，预测未来5、10、15a咸淡水界面变化情况；二是保持目前地下水开采状态，卤水开采在5a后停止，预测未来10、15a海咸水入侵情况。预报过程中，保持研究区水文地质参数不变，大气降水量按照2001-2015年的实际降水序列加载到模型中。在上述条件下，利用验证后的模型进行海咸水入侵预测。

3 结果与分析

预测情景一和情景二的前5a是相同的，均保持目前地下水资源和卤水资源的开采量，预测结果见图9；图10，图11分别为按情景一和情景二预测的10、15a海咸水入侵情况（图中地下水位单位为m，矿化度单位为mg/L）。

从图9可以看出，若保持目前地下水开采和卤水资源开发状态，则卤水区水位降速大于地下水资源开采区的水位降速，海咸水入侵锋线向卤水区推进。

从图10可以看出，按照情景一保持目前地下水开采和卤水资源开发状态，未来10、15a卤水区水位降速仍大于地下水资源开采区的水位降速，海咸水入侵锋线向卤水区推进，且卤水区地下水矿化度大幅度降低。

从图11可以看出，按照情景二保持目前地下水开采和卤水资源开发5a后，停止卤水资源的开发，仍保持地下水开采，则地下水资源开采区水位不断下降、卤水区地下水位逐渐恢复，由于补给水源矿化度相对较低，因此卤水区地下水矿化度有所降低，但是之前卤水区矿化度太高，所以其变化不太明显。海咸水入侵锋线与前5a基本一致，主要原因是海咸水入侵锋线处于地下水降落漏斗中心处，南部低于锋线值和北部高于锋线值均向漏斗区汇聚，使得海咸水入侵锋线处于相对稳定状态。

4 结论

（1）潍河下游典型卤水开采区保持目前地下水和卤水资源开采状态时，卤水区水位降速大于地下水资源开采区的水位降速，海咸水入侵鋒线向卤水区推进，且卤水区地下水矿化度大幅度降低；若保持地下水开采、停止卤水资源的开采，则卤水区地下水位逐渐恢复、地下水矿化度有所降低。

（2）研究区海咸水入侵锋线处于相对稳定状态，主要原因是其处于地下水降落漏斗中心处，南部低于锋线值和北部高于锋线值均向漏斗区汇聚，使得其处于相对稳定状态。

参考文献：

[1]邹祖光，张东生，谭志容.山东省地下卤水资源及开发利用现状分析[J].地质调查与研究，2008，31（3）：214-221.

[2]林存菊，姚英强，付娟.黄河三角洲高效生态经济区卤水资源开采潜力评价[J].山东国土资源，2014，30（9）；48-52.

[3]陈广泉，徐兴永，彭昌盛，等.莱州湾地区海水入侵灾害风险评价研究[J].自然灾害学报，2010，19（2）：103-112.

[4]苗青，陈广泉，刘文全，等.莱州湾地区海水入侵灾害演化过程及成因[J].海岸工程，2013，32（2）：69-78.

[5]OAHMAN K，LARARI A.Evaluation and Numerical Modelingof Seawater Intrusion in the Gaza Aquifer（Palestine）[J].Hydrogeology Journal，2006，14（5）：713-728.

[6]KARAI4ANOGLU N，DOYURAN V.Finite Element Simula-tion of Seawater Intrusion into a Quarry-Site Coastal Aquifer，Kocaeli-Darica，Turkey[J].Environmental Geology，2003，44（4）：456-466.

[7]MAHLKNECHTJ，MERCHAN D，ROSNER M，et al.As-sessing Seawater Intrusion in an Arid Coastal Aquifer UnderHigh Anthropogenic Influence Using Major Constituents，Srand B Isotopes in Groundwater[J].Science of the Total En-vironment，2017，587：282-295.

[8]CIMINO A，COSENTINO C，OIENI A，et al.A Geophysicaland Geochemical Approach for Seawater Intrusion Assessmentin the Acquedolci Coastal Aquifer（Northern Sicily）[J].Envi-ronmental Geology，2008，55（7）：1473-1482.

[9]林文盘，尹泽生.莱州湾海水入侵灾害防治研究进展[J].中国减灾，1992，2（1）：31-34.

[10]章斌，宋献方，韩冬梅，等.运用數理统计和模糊数学评价秦皇岛洋戴河平原的海水入侵程度[J].地理科学，2013，33（3）：342-348.

[11]李国敏，陈崇希.海水入侵研究现状与展望[J].地学前缘，1996，3（1）：161-168.

[12]段梅.硇洲岛海水入侵数值模拟[J].地质灾害与环境保护，2016，27（1）：108-112.