（河北省地矿局秦皇岛矿产水文工程地质大队，河北 秦皇岛066001）

海水入侵是现代社会所特有的环境问题，严重破坏了我国部分沿海地区的环境，阻碍了当地工农业和城市经济的发展。目前海水入侵问题已经引起了国际社会的广泛关注[1]。秦皇岛是我国北方重要的能源输出港口、著名的旅游城市。随着社会经济建设地快速发展和城市化进程地加快，地下水开采强度不断加大，打破了地下水的天然平衡，导致海水入侵，使地下水环境恶化，从而引起区域环境的破坏和生态系统的失衡，海水入侵地质灾害频发[2]。

Visual Modflow是由加拿大Waterloo水文地质公司在原Modflow 软件的基础上应用可视化技术开发研制的，其独特的可视化菜单结构允许操作者简单快捷地定义模拟区范围、选择单位、方便地分配模型的属性和边界条件[3]。这个软件包主要包括Modflow（水流评价）、Modpath（平面和剖面流线示踪分析）、MT3D（溶质运移评价）和Zone Budget（水量均衡计算）四大模块。

1 海水入侵剖面选取

本研究选取秦皇岛市海水入侵氯离子剖面—HQ6剖面作为示例剖面（见图1）。

图1 HQ6模拟区域示意图

本区咸淡水界面存在一个较宽的过渡带，滨海含水层中的海水入侵问题是可混液体间的水动力弥散问题[4]。鉴于本区海水入侵的特征和剖面观测数据的限制，本研究建立了一个考虑密度差异的三维对流－弥散海水入侵模型[5]。

2 数学模型建立

该模型的数学模型由地下水三维水流运移模型和污染质运移模型耦合而成[6]。Modflow是一个三维有限差分地下水流动模型，方程如下：

∂/∂x（Kxx∂h/∂x）+∂/∂y（Kyy∂h/∂y）+∂/∂z（Kzz∂h/∂z）-W=Ss∂h/∂t

其中：Kxx，Kyy和Kzz为渗透系数在x，y和z方向上的分量。假定渗透系数的主轴方向与坐标轴方向一致，量纲为（LT-1）；h为水头（L）；W为单位体积流量（T-1），用以代表流进汇或来自源的水量；Ss为空隙介质的贮水率（L-1）；t为时间（T）。

边界条件：

第一类边界（水头边界）：

h—Γ1=h1（x，y，z）（x，y，z）∈Γ1

第二类边界（流量边界）：

K∂h/∂n—Γ2=q（x，y，z）（x，y，z）∈Γ2

第三类边界（混合边界）：

q（x，y，z）—Γ3= k′（h-h0）/(B′）（x，y，z）∈Γ3

式中Γ1、Γ2、Γ3分别表示一、二、三类边界。

初始条件：

h|t=0=h0（x，y，z，t）Kx、Ky、Kz分别为沿X、Y、Z坐标轴方向的水力传导率；h（x，y，z，t）为水头；K为渗透系数。

污染质运移模型（MT3DMS）的基本方程为：（∂（nC））/∂t=∂/∂x[n（Dxx∂C/∂x+Dxy∂C/∂y+Dxz∂C/∂z-CVx）]+∂/∂y[n（Dxy∂C/∂x+Dyy∂C/∂y+Dyz∂C/∂z-CVy）]+∂/∂z[n（Dxz∂C/∂x+Dyz∂C/∂y+Dzz∂C/∂z-CVz）]+GCe-WC，（x，y，z）∈Ω，t＞0

C（x，y，z，t）—t=0=C0（x，y，z），（x，y，z）∈Ω，t=0

C（x，y，z，t）—Г1=C1（x，y，z，t），（x，y，z）∈Г1，t＞0

D·C/（∂n→）—Г2=f2（x，y，z，t），（x，y，z）∈Г2，t＞0

式中：n为有效孔隙度；C（x，y，z，t）为污染溶质浓度；Ce为污染源浓度；D为弥散系数张量，用矩阵表示为：

矩阵中各个分量分别为：

V为地下水渗透速率，V=（Vx2+Vy2+Vz2）Vx、Vy、Vz为V在x、y、z轴上的投影Vx=-1/nKx∂h/∂x Vy=-1/nKy∂h/∂y Vz=-1/nKz∂h/∂z

3 海水入侵剖面模型建立

3.1 概念模型

实例区南部近海，北部为洋河，东西向剖面不做研究，可设定为人为边界。自上而下可分为8个含水层。其中，亚粘土透水性较差，为弱透水层；而砂卵石层、砂砾石层、细砂层为透水层；下伏基岩为隔水底板。

根据Visual Modflow 用户手册的格式和操作步骤录入数据，将模拟区平面剖分为10行、95列，垂向分为5层，共剖分8 550个六面体单元。本次建模选择的长度单位为m，时间单位为d，渗透系数单位为m/s，抽水率单位为m/d，垂直补给单位为mm/a，质量单位为kg，浓度单位为mg/L。在软件中生成地形、浅层含水层底板、弱透水层底板和深层含水层底板数据数字高程图形，导入渗透系数等参数。

3.2 模型边界条件

3.2.1 地下水流场边界

模拟区南部近海，北部为洋河，故可以设定为水头边界；下伏基岩为隔水边界；剖面的东西向不考虑，可设定为人为边界。

3.2.2 浓度边界

剖面北向南向海逐渐接近，故将剖面北侧（即模型左侧）设定为淡水定浓度，模型南侧（即模型右侧）设定为近海水定浓度。

3.3 时间离散

本模型总模拟时段为10 a（2008年5月1日至2018年5月1日），考虑到实际情况，不同的剖面划分不同的应力期，起始运移步长设定为0.001 d，递增系数为1.1，最大运移步长设定为200 d。

3.4 模型初始条件

3.4.1 初始流场

采用2008年5月常观井和统测井实际监测数据，结合给定的边界条件作为已知水位值的结点，通过Kriging插值方法生成模型的初始水位场。

3.4.2 初始浓度场

由于氯离子是地下水系统中稳定的离子，不易被含水层的孔隙介质所吸附，而且其浓度与矿化度呈线性关系，因此将其作为评价研究区海水入侵标志。以剖面计算区内2008年5月统测井和常观井的实测浓度值和模型中第一类浓度边界（北部边界和海岸边界）上的结点的浓度边界值作为已知浓度值的结点，通过二次Akima插值法生成模型的初始浓度场。

3.5 模型的验证

本模型选取2008年5月至2012年5月枯水期（5月份取样）作为模型的检验期，选取浓度观测井的实测值与模型计算值进行对比，以验证模型的精度。观测数据和模型运算数据比对见图2。

从观测井和模型运算的数据比对可以看出，模型的运算数据都比较平滑，这主要是由于模型建立对各种实际情况的概化引起的，但是就整体趋势而言，模型和实际观测数据比较切合。

图2 Hq6-6 观测数据和模型运算数据比对曲线

4 海水入侵监测剖面预测

模型充分考虑了研究区的实际情况。模型的模拟再现了2008年5月至2012年5月海水入侵的全过程。由于该段时间内的工况条件没有发生改变，模型运算的结果也与实际情况比较拟合，该过程也完整地反映了Cl-浓度场的变化过程。模型运行初期的浓度变化较大，但是总体趋势与实际观测数据相符，截止到模型运行1 825 d，模型已经完全平衡，为了尽可能准确地预测不同工况条件下的海水入侵变化趋势，按照不同的工况条件给予预测（图3）。

图3 2 800 m3/d 工况条件下1 825 d 氯离子变化

以往的预测针对降水概率等进行不同工况的划分，从模型的输入源进行调整，从而改变了地下水的均衡条件。本次模拟采用改变抽水井抽水量的方式，从排泄项方面来改变地下水的平衡，从而模拟不同工况条件下的海水入侵状况。

模型0～1 825 d的稳定抽水最终确定为2 800 m3/d（图3），分两种工况条件以1 825 d作为预测的起点，分别以2 400 m3/d、3 100 m3/d，进行模拟预测2015年5月、2018年5月、2020年5月的入侵情况（见图4、图5）。

图4 2 400 m3/d 工况条件下氯离子浓度变化

图5 3 100 m3/d 工况条件下氯离子浓度变化

5 模拟结果分析

由模拟结果可以看出，海水入侵剖面在抽水量2 800 m3/d工况条件下，模型运行1 825 d基本区域稳定。以1 825 d为预测起点，在2 800 m3/d工况下，海水明显呈海退现象，而在3 100 m3/d工况下，海水入侵继续发生。

因此秦皇岛海水入侵灾害最主要是由不合理的人类活动（超采地下水，沿海工程）引起的，地下水补给量的减少则对海水入侵起到加速作用。因此，适当地减少开采量，防止地下水位降低，可以有效地缓解海水入侵，实现地下水资源的可持续开发利用以支持该区经济社会的可持续发展。

6 结论

（1）应用Visural Modflow软件建立海水入侵剖面模型可行，该软件可有效地将海水入侵可视化，以直观反映海水入侵问题，所建模型和预警系统可信度高，为进一步的研究提供了可靠依据。

（2）本区的地下水位下降是由于大量开采地下水所造成的，地下水位的快速下降引起海水入侵，使地下水中Cl-含量急剧增加。从模型的预测结果看，适当地减少开采量，防止地下水位降低，可以有效地缓解海水入侵，实现地下水资源的可持续开发利用以支持该区经济社会的可持续发展。保护好秦皇岛地下水资源，防治海水入侵，必须从本区的具体地质，水文地质条件，地下水开采现状及水资源状况等实际情况出发，确定有效的保护方法，才能防止海水入侵的发展。

[1]高学平，杨建华，涂向阳，等.帷幕防治海水入侵的数值模拟研究[J].海洋环境科学，2006，25（4）：55-58.

[2]孙娟，杨燕雄.秦皇岛市海水入侵特征[J].中国环境管理干部学院学报，2007，17（2）：51-53，68.

[3]薛禹群，谢春红.地下水数值模拟[M].北京：科学出版社，2007.

[4]艾康洪，陈崇希.漫尾岛咸淡水界面运移剖面二维流水质模型模拟研究[J].勘探科学技术，1994（6）：3-9.

[5]成建梅，陈崇希，吉孟瑞，等.山东烟台夹河中、下游地区海水入侵三维水质数值模拟研究[J].地质前缘，2001，8（1）：179-184.

[6]李俊亭.地下水流数值模拟[M].北京：地质出版社，1989.