硇洲岛海水入侵数值模拟

段梅

(广东省地质局第四地质大队，湛江　524049)



硇洲岛海水入侵数值模拟

段梅

(广东省地质局第四地质大队，湛江524049)

摘要：通过对硇洲岛地质和水文地质条件的综合分析，建立了本研究区三维地下水数值模型，并获取了相关含水层参数。从而建立海水入侵数值模型，对硇洲岛在不同开采条件下的海水入侵趋势进行了模拟预测，结果表明，减少开采量能减缓硇洲岛海水入侵；在满足当地供水需求的前提下，减少浅层水开采量，在远离沿海区域的中部地区适当增加中层水开采量的开采方案优于现状开采条件下的方案。本研究也为减缓海水入侵提供了相关科学依据。

关键词：硇洲岛；海水入侵；数值模拟；模型预测

海水入侵是指由于滨海地区地下水动力条件或水文地质条件发生变化,引起海水或高矿化咸水向陆地淡水含水层运移而发生的水体侵入的过程和现象。海水入侵早在19世纪就引起了人们的注意，Badon Ghyben和Herzberg分别于1889年和1901年独立的提出了著名的计算咸淡水界面的Ghyben-Herzbeg公式，开创了海水入侵的研究序幕[1]。100多年来，海水入侵经历了从理论假设到合理概化，从理论模型、室内实验模型到数值模型等多个阶段的发展过程。Glover[2]、Cooper[3]、Henry[4]等研究了荷兰、以色列、美国等地的咸淡水混合带或过渡带，用可混溶流体的对流-弥散方程代替不可混溶的锋面表达式，沿这一方向经过多年研究形成了对流-弥散模型。Gupta 和 Yapa等[5]运用对流-弥散模型研究泰国曼谷附近的海水入侵。国内关于海水入侵的研究起步较晚，李国敏等[6]利用人工弥散加权方法建立三维有限元模型，对广西北海涠洲岛的海水入侵进行了研究。卢薇等[7]运用FEFLOW软件，利用数值法对珠江口东岸地区的海水入侵进行了数值模拟研究。梁池生等[8]采用数值模拟方法建立的管理模型对雷州半岛防海水入侵进行了定量优化研究，获得了不同区域在给定的地下水水位降深约束下的允许开采量，为雷州半岛地区防海水方案制定提供了重要依据。任雅娴[9]通过对硇洲岛地下水水质状况分析，认为海水入侵是造成地下水咸化的主要原因。罗强等[10]对湛江市硇洲岛的水文地质条件和地下水质变咸的现象进行了分析，发现当地海水入侵受超量开采和区域水文地质条件改变两个因素的共同影响，提出了及早涵养和恢复浅层水，谨慎、合理开采中、深层承压水和调整产业发展规划等防治对策。以上对硇洲岛地区海水入侵的相关研究已表明当地海水入侵现象已相当明显，并有加剧趋势。为了减轻海水入侵造成的危害，本研究从硇洲岛地质、水文地质条件和地下渗流力学的基本原理出发，选取海水中最主要的稳定常量元素氯离子为特征指标，利用数值法对硇洲岛海水入侵进行了模拟研究，预测不同开采条件下的海水入侵趋势，为了防治海水入侵及减缓其扩张趋势提供依据。

1研究区概况及水文地质条件

硇洲岛是雷州半岛近岸一个独立火山岛屿，位于湛江市东南方向，面积约56 km2。地面标高10～50 m，地形沿制高点向四周辐射缓降。岛上地表水系不发育，居民生活用水及灌溉用水、工业用水均依赖地下水。浅层火山岩孔洞裂隙水和中层承压水是岛内的主要供水水源。20世纪90年代以来，香蕉种植和高位养虾成为当地的支柱产业并不断发展，导致地下水开采强度增加，地下水水位下降，地下水水质变咸。尤其是近来来，海水入侵范围扩大的速度大大增加，已经对当地居民的生产生活造成了比较严重的影响。

硇洲岛位于雷琼断陷自流盆地东北部，自第三纪以来沉积了一套厚愈千米的松散-半固结沉积物，由于后期火山喷发导致地表普遍覆盖着一层数米～数十米的火山岩被。其中，由砂性土和孔洞裂隙发育的火山岩组成的含水层与相对隔水的粘性土层相互叠置，构成了多层结构含水层组。根据地下水埋藏条件、水力特征以及本岛的实际情况，岛内的地下水自上而下依次可分为第一含水岩组即浅层潜水-微承压水(埋藏深度一般≤30 m)、第二含水岩组即中层承压水(埋藏深度一般在30～200 m)和第三含水岩组即深层承压水(埋藏深度一般>200 m)三大含水层组，三层含水岩组之间各有一套厚20～30 m粘土层，为相对隔水层。

第一层地下水由火山岩孔洞裂隙水及湛江组上部微承压孔隙水共同组成，受降水补给，地下水流向由中心向四周辐射，最终以泉排泄或汇入大海或者消耗于人工开采和越流补给中层地下水；富水性中东部火山口富水性最好，向四周依次降低，主要受含水层厚度控制。第二层地下水赋存于湛江组的砂砾石中，富水性较丰富，地下水位埋深最深47 m，火山口附近一般大于27 m，岛屿边缘一般小于8 m，水位随开采量增加逐年下降。第三层地下水赋存于下洋组的砾石、砂土中，地下水由于埋深过大，供水意义极小。第二、三含水层与陆地同层含水层相连，第二、三层地下水除主要接受上层水越流补给外，还接受来自大陆方向侧向补给及通过火山口垂向补给，但其侧向补给微弱，故与第一层地下水具有相似的流向和水力特征。

硇洲岛为一个相对独立的水文地质单元，研究区源汇项主要为降雨补给及人工开采。由于地处亚热带气候区，蒸发强烈，灌溉用水一般用于蒸发，灌溉回灌与蒸发项相抵消。本研究区第一与第二含水层之间以及第二与第三含水层之间存在弱透水层粘土，火山口也是各含水层之间的很好连接通道，使三层含水层地下水之间存在水力联系。故可将研究区概化为三层结构的非稳定的地下水三维流动系统。

2海水入侵数学模型

本次海水入侵的数学模型包括水流模型和氯离子迁移模型两部分。首先要求解水流数学模型，然后将流速代入氯离子运移模型求解氯离子随时间与空间的变化规律。

在研究区，影响地下水压力分布的地表水主要是海水，由于海水的密度大于淡水的密度，因此，在研究海水入侵地下水时，若考虑了海水的密度，则计算结果更加合理。描述海水入侵的数学模型，变密度水流方程为

式中，n为有效孔隙度；q为源汇项，表示单位体积多孔介质源或者汇的流量；ej为重力方向单位矢量的第j个分量(e1=e2=0,e3=1)；ρs为海水密度；ρ为地下水密度。

描述海水入侵盐分运移的对流-弥散方程为

(1)

(2)

地下水水动力弥散参数是不确定性程度很强的参数，受很多因素影响，在实际研究中，通常将水动力弥散系数表示为

式中，u是合速度；ux、uy和uz分别是x、y和z轴方向上u的分速度；αT和αL分别是纵向弥散度和横向弥散度；Dx0、Dy0和Dz0分别是分子扩散系数。

采用有限差分法求解上述三层结构的地下水准三维非稳定流定解问题。在数值模拟计算中，我们采用波前法求解最终需要求解的线性代数方程组。

3地下水数值模型

3.1网格剖分及模拟期确定

根据硇洲岛33个钻孔及地质图、剖面图[11]建立其概化含水层三维结构模型。模型共分5层，共17 670个单元格，单元格最小100 m×110 m，最大200 m×200 m。研究区自1992年开始大量开采地下水。在此之前，地下水位基本处于天然的稳定状态。故选取模拟期为1992年1月1日至2010年1月1日。模型共18个应力期，一个应力期为1 a。每个应力期划分12个时间步长。

3.2边界条件

第一含水层被海水切割，海岸线是其定水头边界，因此将海水作为第一含水层的水头边界，由于模型中以月为最小时间步长，因此在本次模拟中忽略潮汐作用的影响。第一层地下水边界取定水头边界，水头值为0 m。第二、三含水层向海域延伸，根据“等效排泄边界”的原理，把其边界选在距海岸线5～15 km处，为定水头边界[12，13]。根据2010年野外调查中研究区边界上三口第二层地下水开采井的水位分别为-0.05 m、-0.2 m和-0.8 m，与0 m海面高程相差不过1 m。因此设第二层地下水边界为定水头边界，水头值同浅层取0 m。第三层地下水埋深一般大于200 m，供水意义小，对第一、二层地下水影响较小，故设其边界条件与第二层地下水一致。

3.3初始条件和参数分区的确定

由于硇洲岛地区缺乏初始水位及参数分区等资料，我们先建立该区的稳定流模型。利用1981年广东省地质图测绘大队绘制的1∶20万综合水文地质图中硇洲岛部分第二层地下水5 m等水位标高线[14]进行模型初步拟合，并且根据含水层的岩性、富水性的不同，分别把第一含水层和第二含水层的水文地质参数分别分为3个参数区，而第一和第二及第二和第三含水层之间的弱透水层及第三含水层不分区,即各自只有1个参数区。拟合结果见图1，水位拟合的总体趋势较好，说明采用的参数分区较为合理。利用该识别后稳定流模型模拟出水位结果即作为地下水流模型的初始水位。

图1　第二含水层水位标高等值线实测值与计算值对比图

3.4源汇项

研究区源汇项主要包括降雨入渗和人工开采。源项为降雨补给量，由降雨入渗系数和降雨量资料计算得出。汇项中浅层水开采主要用于生活和农业灌溉，开采分散，开采井密度大，成面状开采。浅层面状反补给量是根据研究区的农业灌溉用水、水产养殖和居民生活用水这三项相加，再除以研究区总面积得到的一个反补给强度。中层开采井抽水量利用模型拟合机井2008年抽水量对比得出。

3.5模型拟合及结果分析

模型的识别拟合是利用研究区内特征点、线、面、体上已知的地下水动态观测资料与模型在该点上的计算值进行对比，以此判断数值模型对地质体的仿真程度，通常要进行反复的修改参数和调整某些源汇项才能达到较为理想的拟合结果。本次利用2010年野外调查收集的20口机井的静水位标高及成井时间进行模型拟合。拟合结果见图2。

图2　第二层地下水水位标高实测值与计算值对比图(黑色为实测值)

我们把水位拟合的误差和浮动百分比计算出来，如表1所示。从表1可以看出，所有机井拟合的误差值最大为1 m，大多数在0.05 m以下。其中14口机井其拟合浮动百分比都在5%以下，拟合效果较好。个别虽然浮动大于100%，但其误差值很小，最大仅为0.62 m。通过综合分析水位拟合曲线和各种信息可知，所建立的地下水流数值模型基本达到模型精度要求，可用于预测计算。

表1第二层地下水水位标高实测值与计算值对比表

孔号成井时间实测值/m计算值/m误差/m浮动百分比/%z40022004-2.58-2.450.135.14z40052007-3-2.670.3311.03z40071991-5.3-5.250.050.86z40082005-0.050.070.12243.37z40112005-0.2-0.530.33166.53z40122002-0.8-0.800.38z40131999-3.72-3.690.030.91z40141974-7.16-7.140.020.3z40152003-2.5-2.450.051.92z40162005-2.75-2.720.031.1z40171985-10.9-9.919.14z40192009-10-9.870.131.27z40212001-3.5-3.60.12.99z40221968-7.55-7.60.050.7z40251969-12.5-12.010.493.95z40272005-4.4-4.350.051.05z40281995-5.5-5.410.091.55z40311989-2.5-2.560.062.48z40322009-3.18-3.210.031.06z40331991-0.57-1.190.62109.49

4硇洲岛海水中氯离子迁移数值模拟

氯离子迁移模型是利用水流计算结果进行溶质运移的耦合计算的。故其网格剖分、边界条件及模拟期都等同于硇洲岛地下水流动模型中的设定。

根据海水中氯离子的平均浓度，在模型中设边界的氯离子的补给浓度为20 000 mg/l，区域背景值假定为0 mg/l[15]。研究区所收集的资料中没有关于进行溶质运移计算所需要的参数资料，参考相关资料，设定纵向弥散度为20 m，横向/纵向弥散度比率为0.2，垂向/纵向弥散度比率为0.2[16-18]。

由于未收集到数据暂未拟合，但在本模型参数标准下，对氯离子运移进行研究还是有一定的指导意义的。模拟结果与实际情况是基本吻合的。模拟结果表明，以氯离子浓度为200 mg/l作为判断海水入侵的标准，硇洲岛第一含水层沿海1 km以内地区几乎全部受到了海水入侵的影响。第二含水层在西部和北部沿海也已逐步受到浅层水咸化的影响。该结果与陈苑[19]2011年对硇洲岛地下水咸化与海水入侵调查得出硇洲岛内地势较为低缓的西南、西北与东北部沿海区域已经出现海水入侵，海水最大纵深渗透距离达 1 000 m 左右，入侵区内 30～80 m 深度范围的地下含水层水质普遍发生不同程度咸化的结论是相符的，这表明本次建立的模拟模型是合理的。

5海水入侵预测

5.1现状开采条件下2020年海水入侵预测

保持开采量及边界条件不变，降雨补给量取1992年～2009年平均值0.004 0 m/d，模拟时段为2010年～2020年，共10个应力期，每个应力期12个时间步长。模拟结果表明浅层含水层受海水入侵影响的范围在这10 a有所扩展，平面扩展速率约20～40 m/a。由于第一含水层的地形优势，除沿海地区海水向淡水渗流，其余大部分地区地下水流向仍保持由高地势地区向低地势地区渗流。所以第一含水层的海水入侵速率相对来说有所降低。

第二含水层受海水入侵影响的范围急剧增加，几乎占了全区面积的三分之一。2009年底第二含水层仅西部和北部沿海地区受到海水入侵的影响，而经过10 a的发展，几乎所有沿海地区全部发生了海水入侵。第二含水层是硇洲岛地区的主要供水层位。假如维持现状开采，仅仅10 a第二含水层就将面临着十分严重的海水入侵这一地质环境问题。

5.2优化开采条件下2020年海水入侵预测

本次分别预测了以下两种开采方案下，2020年氯离子的浓度分布特征：方案1,第一含水层开采量减小而第二含水层开采量不变；方案2，第一含水层开采量减小但第二含水层开采量在远离沿海区域的中部地区增加以保证总开采量稳定。

方案1模拟结果与维持现开采条件下的模拟结果进行对比，第一含水层水位有所上升，海水入侵的范围减小；第二含水层仍不可避免的受到了海水入侵的影响，扩展范围比现状开采条件下略有缩小。

方案2模拟结果表明增加第二层地下水开采量没有导致第一含水层海水入侵的进一步发展，第二含水层受海水入侵的影响范围相对于方案1略有增加，证明了第二含水层的海水入侵程度在一段时间内会无可避免的加重，必须采取除减小开采量以外更有效的措施来阻止海水入侵程度的发展。

为了更好对比上诉三种开采方案，选取z4022井所处单元格在这三种开采条件下中该单元格中氯离子浓度随时间变化的特征，来具体分析三种开采条件下对中层含水层海水入侵程度的影响。如图3所示。

图3不同开采条件下z4022所处单元格氯离子浓度随时间变化曲线

从图3中可以看出，虽然方案2中增加了中层含水层的开采量，但z4022所处单元格中氯离子浓度却低于维持现状开采条件下的氯离子浓度。这说明了在不减少开采量的情况下，合理的开采布局是可以在一定程度上抑制海水入侵的发展的。

6结论

海水入侵对硇洲岛的社会和经济发展会带来不利影响，因此本文根据硇洲岛水文地质的特点和钻孔资料，通过水文地质条件概化建立硇洲岛的三维稳定流模型及瞬变流模型到氯离子的迁移模型，运用收集到的地下水位资料对模型进行识别和校验，运用识别后的模型模拟预测不同开采条件下海水入侵趋势。结果表明，减少开采量能减缓海水入侵；在满足当地的供水需求的前提下，减小浅层水开采量，在远离沿海区域的中部地区适当增加中层水开采量的开采方案优于现状开采条件下的方案。

由于硇洲岛专门的海水入侵勘探工作量较少，缺乏长观资料和边界资料及试验数据，现有的机井观测点分布不均，本研究仅为初步模拟研究。

特别感谢张莉在本文写作过程中提出宝贵建议和给以的帮助！

参考文献

[1]蔡祖煌，马凤山．海水入侵的基本理论及其在入侵发展预测中的应用[J].中国地质灾害与防治学报，1996,7(3)：1-9.

[2]Glover R E．The pattern of fresh water flow in a coastal aquifer. J Geophys Res, 1959, 64(4):457-459.

[3]Cooper H H, F A Kohout, H R Henry, and R E Glover．Sea water in coastal aquifers.U.S.Geol. Surv. Water Supply Paper, 1964, 1614-C:70-84.

[4]Henry H R．Effect of dispersion on salt encroachment in coastal aquifers.U.S.Geol. Surv., Water Supply Paper. 1964, 1613-C: 70-84.

[5]Oude Essink G H P, Boekelman R H．Problems with large-scale modeling of salt water intrusion in 3D. 14thSalt Water Intrusion Meeting, 1996.

[6]李国敏，陈崇希，沈照理，等．涠洲岛海水入侵模拟[J].水文地质工程地质，1995，(5)：1-5.

[7]卢薇，朱照宇，刘卫平.基于FEFLOW的海水入侵数值模拟[J].地下水，2010，32(3)：19-21.

[8]梁池生,周训,陈明佑．雷州半岛防海水入侵最优地下水开采方案探讨[J].现代地质，2002,16(4)：429-434.

[9]任雅娴.硇洲岛地下水咸化成因及变化趋势分析[J].广东水利水电，2013，(5)：23-27.

[10]罗强，苏肇汉.硇洲岛地下水海水入侵的特征分析及治理对策[J]，地质灾害与环境保护，2007，18(2)：28-32.

[11]广东省地质矿产局.湛江综合区域地质调查报告(1:50000)[R].1989.

[12]陈崇希,林敏，舒本媛，等．滨海承压含水层地下水的等效排泄边界——以北海禾塘水源地为例[J]．水文地质工程地质，1990，17(4):2-4．

[13]成建梅,黄丹红,胡进武.海水入侵模拟理论与方法研究进展[J].水资源保护,2004，20(2):3-8.

[14]广东省地质局．雷州半岛区域水文地质普查报告(1:200000)[R].1981.

[15]付旭.泉州市典型区域海水入侵实验与数值模拟研究[D].华侨大学,2013.

[16]孙讷正.地下水水质的数学模拟(四)水动力弥散方程的数值解法[J].水文地质工程地质，1982,4(1):49-55.

[17]邹立芝,潘俊,杨昌兵，等.含水层水力参数的尺度效应研究现状[J].长春地质学院学报,1994,24(1):66-69.

[18]张宇峰，张雪英，徐炎华，等.土壤中水动力弥散系数的研究进展[J].环境污染治理技术与设备，2003,4(7):8-12.

[19]陈苑.硇洲岛地下水水质咸化和海水入侵主要成因及对策[J].大众科技，2011，(4)：114-115.

NUMERICAL SIMULATION OF SEAWATER INTRUSION IN NAOZHOU ISLAND

DUAN Mei

(The Fourth Geology ofGuangdong Geological Bureau ,Zhanjiang Guangdong524049,China)

Abstract:Through the comprehensive analysis of geological and hydrogeological conditions of the NaoZhou island, we have established the three-dimensional numerical model of groundwater and obtained the aquifer parameters in the study area. The trend of seawater intrusion of the NaoZhou island under different mining conditions have been simulated and predicted through the seawater intrusion numerical model. The results showed that reducing production could slow down the seawater intrusion of the NaoZhou island; and on the premise of meet the demand of the local water supply, the solution scheme that reduces shallow water production, and increases the middle water production appropriately which is away from the central part of the coastal areas, is superior to the present situation of mining conditions. This study also provided related scientific basis for slowing down seawater intrusion.

Key words：Naozhou Island； seawater intrusion; numerical simulation； model prediction

作者简介：段梅(1985-)女，工程师，环境工程专业，主要从事水文、工程和环境地质研究工作。E-mail:41548316@qq.com

中图分类号：P736;P641

文献标识码：A

收稿日期：2015-08-20改回日期：2015-11-17

文章编号：1006-4362(2016)01-0108-05