基于有限差分模拟的海水入侵滨海承压含水层水位变化研究

Study on Water Level Change of Seawater Intrusion into Coastal Confined Aquifer Based on Finite Difference Simulation

摘要: 针对滨海地区承压含水层海水入侵问题，本文基于有限差分法，建立了地下水二维稳定流数值模型。运用有限差分数值模拟软件，计算获得了海平面上升的承压含水层压力水头变化曲线，分析了由于海平面上升所引起的滨海承压含水层水位变化规律。在此基础上，结合具体算例分析了海平面上升高度、含水层渗透系数、截渗墙嵌入深度等因素对海水入侵过程的影响规律。模拟结果表明，此数值模型能准确模拟滨海地区海水入侵过程，并预测滨海承压含水层水位变化，可为滨海地区海水入侵现象的科学预测和防治保护提供支撑。

Abstract: In this paper, to address the problem of seawater intrusion into the confined aquifer in the coastal area, a two-dimensional numerical model of groundwater steady flow is established based on the finite difference method. By using the finite difference numerical simulation software, the pressure-head change curve of the confined aquifer with the rise of the sea level is obtained, and the change rule of the water level of the confined aquifer in the coastal area caused by the rise of the sea level is analyzed. On this basis, the influence of sea level rise height, permeability coefficient of aquifer and embedding depth of the cut-off wall on the process of seawater intrusion is analyzed with specific examples. The simulation results show that this numerical model can accurately simulate the process of seawater intrusion in the coastal area and predict the water level changes in the coastal confined aquifer, which can provide support for the scientific prediction of seawater intrusion in the coastal area and the prevention and protection of seawater intrusion phenomenon.

1. 引言

滨海地区是人类经济活动最活跃的地区，据统计，全球约4成的人口生活在沿海100公里范围内，滨海地区的人口平均密度约为全球平均水平的2倍。沿海含水层是全世界多达10亿人的主要饮用水来源 [1] ，而近些年越发频繁的海水入侵现象使得有限的地下水资源与居民的用水需求矛盾日益尖锐，对滨海含水层带来了致命影响 [2] [3] 。

海水入侵是因为滨海地区地下水动力条件发生改变，使得海水逐渐向陆地淡水含水层运移，而发生的水体侵入过程和现象 [4] 。目前，世界上已经有许多国家和地区发现了海水入侵现象，例如德国、意大利、法国、西班牙、葡萄牙、英国、澳大利亚、美国、中国等 [5] 。海水入侵会导致咸淡水界面逐渐向内陆含水层移动，造成滨海地区水质恶化、灌溉用水源地减少、耕地资源退化、土壤生态系统失衡、影响农业生产和工业生产等影响 [6] 。影响海水入侵的因素主要分为自然因素和人为因素，其中自然因素主要是由于海平面不断上升，而人为因素主要是由于过度开采地下水，二者均会导致滨海地带含水层的水动力条件发生改变。

随着计算机技术的发展，数值模拟逐渐被应用于地下水渗流问题的研究，作为一种地下水位变化与溶质运动趋势预测工具，这种方法可以刻画复杂的水文地质条件并耦合模拟地下水流运动与溶质运移过程，弥补了解析方法适用性的不足，常用的数值模拟法包括有限差分法、有限单元法、有限条法、有限层法等。其中，有限差分法因其便捷性和通用性较强深受广大学者喜爱，目前在地下水渗流问题中的应用最为广泛。易连兴(2007)以参数分区方法建立了单元非等参有限差分方程，能够更准确地处理地下水非均质各向异性渗流问题 [7] 。龚嘉临等(2017)探讨了MODFLOW有限差分模型和MT3DMS溶质运移模型在地下水与湖泊的相互作用研究中广泛应用的特点及应用 [8] 。杨文琳等(2017)针对地下水稳定流和溶质运移模型应用有限差分法进行求解分析，并以生活垃圾填埋场为例，模拟了该垃圾填埋场在非正常状况下运行对周边潜水含水层水质产生的影响 [9] 。吴晓莉等(2018)针对地下水稳定流和溶质运移模型应用有限差分法进行求解分析，模拟了预处理工艺区在非正常状况下运行对周边潜水层水质产生的影响 [10] 。此外，随着学者对海水入侵发生机理认识的加深，变密度海水入侵模拟模型被开发应用，如da Silva Alves等(2022)在Visual MODFLOW软件中引入SEAWAT插件，为研究区域含水层构建了一个数值模型，并模拟了目前开采量约为750 m³/小时的抽水井在20年内的海水入侵情况，获得了旱季和雨季补给量变化以及水井开采流量减少的模拟结果，有助于当地的水文地质管理 [11] 。李岱远等(2023)采用SEAWAT方法构建了海水入侵数值模型，通过情景分析法定量评估了海水入侵影响因子，为研究区域提出了合理预测未来海水入侵发展趋势的建议和地下水开发利用解决方案 [12] 。

综上所述，运用有限差分法对地下水渗流问题和海水入侵数值模拟研究已取得了一定进展，但对于针对由海平面上升所引起的滨海承压含水层水位变化问题的研究还较为薄弱。为了更真实地反映海水入侵过程，本文基于有限差分法，建立了地下水二维流数值模型，在模型中设置截渗墙，对由于海平面上升导致水动力条件发生改变而引起的滨海承压含水层水位变化进行了深入研究。并结合具体算例分析了海平面上升高度、含水层渗透系数、截渗墙嵌入深度等因素对海水入侵过程的影响规律。

2. 模型建立

2.1. 数学模型

为准确模拟海水入侵引起滨海含水层水位变化过程，我们建立了一个二维地下水渗流模型，如图1所示。通过将坐标原点设置在近海含水层左下角来建立坐标系，取z轴方向向上为正；在两含水层中间设置截渗墙，截渗墙下端开口长度为D，将模型划分为两个区域。一是靠近海水的近海含水层，厚度为m，计算长度为B；二是位于截渗墙右侧的内陆承压含水层，厚度为M，计算长度为C。设置近海含水层左侧边界与上部边界为定水头边界，取值为h_a；将内陆承压含水层右侧设置为定水头边界，取值为h_b，其余边界均设置为隔水边界。含水层中初始地下水位为h₀。建立数学模型所采用的基本假设如下：

1) 内陆承压含水层具有有限的横向尺寸；

2) 含水层中的流动遵循达西定律；

3) 含水层是均质，正交各向异性、等厚且水平；

4) 止水帷幕的厚度在研究中被忽略不计。

Figure 1. Schematic diagram of the computational model for seawater intrusion

图1. 海水入侵计算模型示意图

2.2. 控制方程与边界条件

近海含水层数学模型为：

$\{\begin{array}{l}{K}_x\frac{{\partial }^2{h}_1}{\partial x^2}+K_z\frac{{\partial }^2{h}_1}{\partial z^2}=0\left(0\le x\le B\right)\\ \frac{\partial {h^{\prime }}_1}{\partial z}|\begin{array}{l}\\ {}_{z=0}\end{array}=0\\ {h^{\prime }}_1|{}_{z=m}={h^{\prime }}_1|{}_{x=0}=0\\ \frac{\partial {h^{\prime }}_1}{\partial x}|\begin{array}{l}\\ {}_{x=B}\end{array}=\{\begin{array}{l}q\left(z\right),0\le z\le D\\ 0,D\le z\le m\end{array}\end{array}$ (1)

式中：h₁为近海含水层地下水水头高度，L； ${h^{\prime }}_1$ 为近海含水层地下水水头变化高度，计算公式为 ${h^{\prime }}_1=h_1-h_a$ ，L；K_x、K_z分别为含水层水平和垂直渗透系数，LT⁻¹；B为近海含水层宽度，L；m为近海含水层厚度，L；q(z)为两含水层公共边界流量函数；D为截渗墙下开口长度，L。

内陆承压含水层数学模型为：

$\{\begin{array}{l}{K}_x\frac{{\partial }^2{h}_2}{\partial x^2}+K_z\frac{{\partial }^2{h}_2}{\partial z^2}=0\left(B\le x\le L\right)\\ \frac{\partial {h^{\prime }}_2}{\partial z}|\begin{array}{l}\\ {}_{z=0}\end{array}=\frac{\partial {h^{\prime }}_2}{\partial z}|\begin{array}{l}\\ {}_{z=M}\end{array}=0\\ {h^{\prime }}_2|{}_{x=L}=0\\ \frac{\partial {h^{\prime }}_2}{\partial x}|\begin{array}{l}\\ {}_{x=B}\end{array}=\{\begin{array}{l}q\left(z\right),0\le z\le D\\ 0,D\le z\le M\end{array}\end{array}$ (2)

式中：h₂为内陆承压含水层地下水水头高度，L； ${h^{\prime }}_2$ 为内陆承压含水层地下水水头变化高度，计算公式为 ${h^{\prime }}_2=h_2-h_b$ ，L；M为内陆承压含水层厚度，L。

两含水层之间的连续性条件为：

${h^{\prime }}_1|{}_{x=B}={h^{\prime }}_2|{}_{x=B}=0$ (3)

2.3. 计算参数选取与网格划分

海水入侵引起滨海含水层水位变化模型计算参数的选取如表1所示。模型计算区域为250 × 40 m，划分采用矩形网格，总计网格数目为9000，每个网格单元尺寸为1 m²，如图2所示。从图中可以看出，左边近海含水层区域的左侧边界与上部边界、内陆承压含水层区域的右侧边界为橙色网格，这是设置了定水头边界的效果，图中未设置橙色网格的边界默认为隔水边界。网格图中的蓝色区域，即为截渗墙设置的区域，截渗墙的水平和垂直渗透系数我们设置为1 × 10⁻¹⁰ m/d，以此来模拟截渗墙在海水入侵过程中对海水的阻隔作用。

Table 1. Table of model calculation parameters

表1. 模型计算参数表

Figure 2. Grid division diagram for the finite difference model

图2. 有限差分模型网格剖分图

3. 海平面上升引起内陆含水层水位变化数值模拟结果分析

3.1. 不同含水层深度处的地下水位变化特征

经过上述操作将模型划分且输入对应计算参数后，运算可得MODFLOW压力水头计算结果，我们选取了z轴方向上5组数据进行研究，运用origin软件将数据整理绘制得到图3，观察下图可以发现，z = 10 m、z = 15 m处含水层内压力水头上升幅度较含水层顶板(z = 40 m)和靠近含水层顶板(z = 30 m)的大很多，这是因为海水是从近海含水层一侧流入内陆含水层中，在海水附近的地下水渗流路径最短，地下水可以在很短时间内得到释放，从而导致靠近海水一侧的水位变化幅度明显，而靠近内陆含水层顶部区域的地下水渗流路径较远，因而水位变化幅度明显较为平缓。同时，从图中可以看出，当x = 50 m时，由于截渗墙的作用，z = 10 m、z = 15 m处含水层内压力水头分布曲线不仅发生偏折，而且在很小的范围内，水位降深的变化幅度较大，反映出截渗墙对海水入侵的阻隔作用。此外，随着水平距离x的增加，不同深度处的压力水头逐渐减小并趋于一致，说明当水平距离增大到一定程度时，不同深度处的压力水头不再发生变化，地下水流趋于一维流。

Figure 3. Pressure-head distribution curve along x-axis at different depths of the aquifer

图3. 不同含水层深度沿x轴压力水头分布曲线图

3.2. 海平面上升高度对海水入侵发展的影响

近十年来，全世界范围内海平面年平均上升速率已经达到4.5毫米，且未来全球海平面加速上升趋势仍将继续。海平面的不断上升，严重威胁着滨海城市和岛屿国家的生存。同时，海平面上升会使滨海地区地下水动力条件发生改变，对海水入侵有重要的影响。就目前形势来看，研究在沿海含水层中海平面上升对海水入侵的影响和怎样采取积极措施对海水入侵进行防治有十分重要的意义。为此，本文建立了二维地下水稳定流模型，通过对近海含水层边界设置不同的定水头数值，以此来模拟海平面升高的情况，为方便观测，本文下述图表中的数据均取自于z = 10 m的含水层。海平面上升引起内陆含水层水位变化结果如图4所示，从图中可以看出，随着海平面高度的不断上升，含水层内压力水头上升幅度也在不断增加。

Figure 4. Head distribution curve under changes in sea level rise heights

图4. 海平面上升高度变化下的水头分布曲线图

3.3. 渗透系数对海水入侵发展的影响

为了研究渗透系数对海水入侵发展的影响，对验证算例中渗透系数取不同的值，模型将K_x与K_z设置相同的数值，记为k，取值情况如图5右上角所示，其余参数同第二节算例。使用GMS对不同渗透系数下的压力水头进行计算，结果见图5。可以看出，当渗透系数由0.1 m/d增大至10 m/d时，x = 100 m处的含水层的压力水头减小了约2.4 m，说明渗透系数增加会使得内陆含水层水位变化幅度减小；并且随着渗透系数增大到一定数值，含水层内水头变化逐渐减小，地下水流逐渐趋于稳定。

Figure 5. Head distribution curve under changes in permeability coefficient

图5. 渗透系数变化下的水头分布曲线图

3.4. 截渗墙嵌入深度对海水入侵发展的影响

图6给出了滨海含水层截渗墙不同埋设深度(l_c = 5 m, 10 m, 15 m)时所得的等势线图。由图可知，由于截渗墙所处位置水平和垂直渗透系数较小，设置为1.0 × 10⁻¹⁰ m/d，等势线出现急剧变化，体现出截渗墙对海水入侵的阻隔作用，对减缓海水入侵过程可以产生有利影响；截渗墙两端等势线呈现拱门的形状，越靠近截渗墙，等势线分布越密集，右侧内陆含水层等势线趋于垂直，且间距逐渐增大；同时，随着截渗墙插入深度的增加，两个含水层之间的等势线分布也更加密集，计算结果均符合地下水水力特性，验证了计算模型的正确性。

Figure 6. Influence of buried depth of cut-off wall l_c = 5 m, 10 m, and 15 m on groundwater level distribution

图6. 截渗墙埋设深度l_c = 15 m、10 m、5 m对地下水位分布的影响

4. 结论

为了更好地开发利用滨海含水层地下水资源，掌握地下水运动规律，并运用这些规律指导实践，为滨海地区海水入侵现象的科学预测和防治保护提供支撑，本文基于有限差分法，针对海平面上升引起的滨海含水层水位变化的问题，建立了承压含水层二维稳定流地下水数值模型，运用有限差分数值模拟软件，计算获得了海平面上升所引起的承压含水层压力水头变化曲线。并结合具体算例分析了截渗墙嵌入深度、含水层渗透系数等因素对海水入侵过程的影响规律，具体结论如下：

1) 越靠近海水的含水层水位变化幅度越明显，并且随着深度的增加，压力水头数值越大，含水层内水位也越高；同时，随着内陆承压含水层水平距离的增大，不同深度处的压力水头逐渐减小并趋于一致，说明当水平距离增大到一定程度时，不同深度处的压力水头不再发生变化，地下水流趋于一维流。

2) 随着海平面高度的上升，海水从近海含水层一侧不断流入内陆承压含水层，使得内陆承压含水层水位上升，将会加快海水入侵速率。

3) 渗透系数增加会使得内陆含水层水位变化幅度减小，并且随着渗透系数增大到10 m/d后，含水层内水头变化逐渐减小，地下水流逐渐趋于稳定。

4) 截渗墙所处位置的含水层内压力水头分布曲线不仅发生偏折，而且在很小的范围内，水位降深的变化幅度较大，等势线分布也越密集，说明在滨海承压含水层设置截渗墙可以对海水入侵过程起到明显的阻隔作用，对减缓海水入侵过程可以产生有利影响。