不透水面影响下海岛淡水透镜体的形成与演化

2022-01-15凌子涵束龙仓

水利水电科技进展 2022年1期

凌子涵，束龙仓，王然

(1.河海大学水文水资源学院,江苏南京 210098； 2.河海大学水文水资源与水利工程科学国家重点实验室, 江苏南京 210098)

我国拥有丰富的海洋资源，海岛作为人类探索海洋的重要支点，对于开发海洋资源、建设海洋经济大国具有重要意义。本文所研究的海岛是西沙群岛中最大的灰沙岛屿，为三沙市政治、经济、文化中心，岛上常住人口2 500多人[1]。近年来随着三沙市对该海岛的大力建设，岛上现已具备良好的居住生活条件，但面临着用水形势严峻的问题，目前主要通过补给船定期补给以解决岛上居民的饮用水问题[2]，但这种方式运水，水量十分有限，且成本昂贵。除此之外还有海水淡化、雨水收集等方法，但均无法满足日趋增长的岛上用水需求。

海岛独特的地质结构使得岛屿地下淡水被咸水包围，并浮于咸水之上，被称作“淡水透镜体”。近年来国内外学者对淡水透镜体的形成、演变及影响因素等进行了深入研究[3-7]。李国敏等[8]对涠洲岛含水层条件进行概化，并建立数值模型，模拟了海水入侵的3D过程。Vacher[9]利用Ghyben-Herzberg公式对Bermuda岛和Great Exuma岛淡水透镜体的形态进行了模拟计算。James等[10]使用SEAWAT软件，模拟海岸侵蚀对美国佛罗里达州海岛淡水透镜体形态与厚度的影响。甄黎等[11]通过物理模型模拟了不同开采条件下倒锥的变化过程。本文在前人研究的基础上，利用Visual MODFLOW软件构建了所研究海岛淡水透镜体的二维模型，并在透镜体稳定后设置不透水面，通过改变不透水面的覆盖范围，模拟透镜体形态及储水量的变化。此外，对影响透镜体储水量的不同因素进行灵敏度分析，为岛屿淡水资源开采提供科学参考。

1 研究区概况

本文所研究海岛为海南省三沙市永兴岛，朝海部分的岛屿经礁坪前缘以大角度进入深海盆地[12]，全岛在北回归线以南，纬度低，地处热带地区，属赤道带、热带海洋性季风气候，岛屿年平均气温在26℃以上，年降雨量超过1 400 mm，其中6—11月降水量占全年85%以上。

该海岛的地质结构大体上可以分为全新统和更新统，地表以下约22 m处为不整合面，不整合面以上为全新世珊瑚贝壳砂地层，淡水透镜体主要存在于该地层；不整合面以下由岩溶十分发育的更新世石灰岩构成，由于原生礁不断受到海水侵蚀，石灰岩体被孔隙通道贯通，岩溶发育，因此渗透性很强[13]。钻孔资料及其他文献[14-19]表明：岛屿地表以下6 m由素填土和珊瑚碎屑中砂组成，地表以下6～22 m由珊瑚碎屑砾和珊瑚碎屑砾砂组成，22 m以下为珊瑚礁灰岩；岛屿介质孔隙度为0.4，给水度为0.2；含水介质纵向弥散度为5 m，横纵比为0.1。此外，该海岛淡水补给的主要来源是降雨入渗补给，排泄方式主要是垂向的蒸发排泄和侧向的向海洋排泄。

2 数值模型

2.1 水文地质概念模型

本文以研究区的水文地质条件为依据，建立二维模型观察淡水透镜体的形态变化。选取岛屿截面长度1 980 m，模型的最大深度取到海平面以下48 m，海平面以上岛屿高程设定为5 m。将模拟区域地层整体分为3层：第一层为地表至地表以下6 m，该层主要是填土及珊瑚碎屑；第二层为地表以下6～22 m，为全新世地层；地表以下22～54 m为更新世地层。参照该海岛的地理特性及文献[20]，3层地层结构渗透系数k的初始值分别取60 m/d、150 m/d、1 000 m/d。

海岛淡水透镜体的补给来源为降雨补给。Ghassemi等[21]根据统计资料对Home岛淡水透镜体进行三维数值模拟时，降雨入渗补给系数取0.44；周从直等[22]对永兴岛淡水透镜体进行三维数值建模时，取多年平均降雨的40%作为补给量。本文结合之前学者的经验，考虑该海岛的水文气象条件，降雨入渗补给系数取0.4。模型其他参数设置如下：淡水密度为1 g/cm3，海水密度为1.025 g/cm3，海水质量浓度为0.019 g/cm3，弹性贮水率为10-5m-1。模型顶部接受大气降雨补给且存在蒸发，设定为流量边界；补给浓度取雨水中氯离子质量浓度，近似为0 g/cm3；模型底部含水层与海水连通，不存在水量交换，因此可视为零流量边界[22]。模型侧向边界在海平面以上部分质量浓度为0 g/cm3,在海平面以下的部分被海水包围，故取海水的氯离子质量浓度，即0.019 g/cm3。模型初始水头取0 m，概化后模型示意图如图1所示。

图1 水文地质模型

2.2 数学模型

由于海水与淡水可以混溶，因此所构建海岛淡水透镜体二维剖面模型存在可混溶液体间的水动力弥散问题[23]。数学模型由控制方程、初始条件以及边界条件构成，依据上述水文地质概念模型，建立研究区地下水流与溶质运移数学模型，通过以下微分方程的定解问题来描述：

(1)

式中：K1为潜水含水层渗透系数；μ为潜水含水层给水度；H为潜水水位；B为潜水含水层底板标高；W1为源汇项；H0为初始水位；H1为研究区一类边界点的水位；x,z为坐标；D为计算区范围；Γ1为研究区一类边界：t为时间。

多孔介质的溶质运移方程包括水动力弥散项、地下水对流项、源汇项以及反应项。溶质运移方程为

(2)

式中：θ为有效孔隙度；C为溶质质量浓度；qs为单位时间内进入单位体积含水层的源汇项体积；ρs为源汇项溶质质量浓度；Dij为水动力弥散系数(i,j= 1, 2, 3)；vi为流体速度；∑Rn为化学物质反应项。

2.3 模型离散化

将研究区域划分为若干个网格单元(共30层34列)，每个网格单元内部的水文参数为常数，单元之间存在水力联系，将分界面以上部分及侧向边界网格进行加密处理。模型模拟期为73 050 d，水头计算初始时间步长为0.001 d。溶质运移计算初始时间步长1 d，增长因子为2，最大时间步长为10 d，模型结果以月为单位输出。

3 模型模拟结果与分析

3.1 淡水透镜体模拟结果

在Visual MODFLOW软件中，根据实际逐月降雨资料与逐月蒸发资料设定模型源汇项，模拟期为100 a，选用SEAWAT程序包进行模拟，其中水流计算采用预处理共轭梯度法(PCG)来求解；对流计算采用中心差分法[24]。运行结果以氯离子600 mg/cm3质量浓度等值线[21]为淡水透镜体下界面轮廓即咸淡水界面。

模型运行结果如图2所示，随着降雨入渗的不断补给，淡水透镜体厚度变大，在50 a后形状达到稳定，稳定时透镜体最大厚度约为15.6 m。图3为岛屿地下水流的流速示意图，箭头的大小表示流速的快慢，可以看出，降雨入渗的水流由上至下流动，并向两侧排泄；透镜体内部中心水流流速很低，由于咸淡水之间存在密度差，使得淡水透镜体中心部分淡水受浮力较大，因此透镜体下界面轮廓中心部分向上凸起。

图2 淡水透镜体形态演变过程

图3 岛屿地下水水流流速

该岛屿淡水透镜体最大厚度随时间的变化如图4所示，可以看出透镜体形成初期(0～10 a)厚度增加速度逐渐加快，到中期(10～40 a)速度逐渐下降，40 a后透镜体厚度变化十分微小直至不变化。如图5所示，透镜体潜水面呈向上凸起的形状，中间厚、两侧薄，中间水头最大超过0.33 m，边缘水头最小为0 m，水头向四周递减形成水力坡降，可以保证淡水从透镜体最厚的地方不断排泄入海。

图4 淡水透镜体最大厚度随时间变化

图5 淡水透镜体潜水面形状

3.2 不透水面对透镜体的影响

岛屿城市化过程中土地利用率不断加大，海水入侵现象更加明显，岛屿淡水储量会随之减少。为探究不透水面对淡水透镜体的影响，在模型0～594 m部分将降雨入渗补给系数设置为0.01，将该部分界面作为不透水面[25]，蒸发量设为0，将原始岛屿模拟的水头及质量浓度结果作为初始条件输入模型中，模型大约经过15 a后再次达到稳定，模拟结果如图6所示，可以看出设置了不透水面一侧的岛屿，由于降雨入渗补给的减少，淡水透镜体向岛屿内部凹陷，淡水透镜体的最大厚度减小至14.6 m，且淡水透镜体厚度整体上变薄，经计算稳定时岛屿地下淡水储量为3 408 m3/m，相比初始时减少了20.23%。

图6 设置不透水面后透镜体形态

图7为设置不透水面前后海岛内部水流的流向，由图7可知，未设置不透水面的岛屿降雨入渗补给的水流到达透镜体形状边缘时，部分水流继续向下流动，部分水流水平流动，还有一些水流有向上的速度分量；在设置不透水面后，不透水面的存在阻碍了局部的降雨入渗补给，在不透水面存在的一侧，水流多为横向流动，向下流动路径缩短，并在岛屿侧边向海洋排泄，因此透镜体轮廓向内凹陷，淡水储量也相应减少。在岛屿开发过程中，随着建设力度的增大，岛屿土地覆盖中不透水面占比会随之增加，因此在模型中设置了多种不透水面覆盖范围，以分析淡水透镜体的形态及储水量变化规律。结果表明，随着不透水面覆盖范围的增大，存在不透水面一侧的淡水透镜体向内萎缩趋势愈加明显，体积不断减小，淡水透镜体最大厚度位置也在逐渐右移。

图7 海岛内部地下水流向

图8展示了不同不透水面覆盖范围占比下岛屿单宽地下淡水储量的变化，随着不透水面覆盖范围占比的不断加大，地下淡水储量逐渐减少，且减少速度逐渐加快，当不透水面覆盖范围占比超过40%后，二者逐渐转变为线性关系。同时可以看出，当不透水面覆盖范围占比超过岛屿总长度的50%时，地下淡水储量减少至原来的50%左右；当不透水面覆盖范围占比为60%时，地下淡水储量已经不足原储量的40%。

图8 岛屿地下淡水储量随不透水面覆盖范围占比变化

3.3 模型灵敏度分析

本文采用扰动分析法中的Morris法作为灵敏度计算方法，表示参数变化对模拟结果影响程度的灵敏度S的计算公式[26]为

(3)

式中：yi为参数变化后模型输出值；y0为模型原始输出值；P0为参数初始值；Pi为参数变化后的值；Δi为Pi相对于P0的变化幅度。

影响海岛淡水透镜体形态及体积的因素有很多，可以分为自然因素和人为因素，自然因素包括降雨入渗补给、蒸发和岛屿介质渗透性等，人为因素主要是不透水面对透镜体的影响。为了研究该海岛地下淡水储量在不同影响因素变化下的敏感程度，本文对4项主要影响因素进行灵敏度分析，采用局部灵敏度分析方法，研究单一影响因素变化时对淡水透镜体体积(文中指地下淡水储量)的影响。依次选用不透水面覆盖范围、含水层给水度、降雨入渗补给系数和渗透系数进行灵敏度计算。

以上4种参数初始值设置如下：不透水面覆盖范围594 m；含水层给水度0.2；降雨入渗补给系数0.4；3层地层结构渗透系数分别为60 m/d、150 m/d和1 000 m/d。为了便于计算，参数变化幅度取±10%、±20%和±30%，灵敏度计算及分析结果见图9。

图9 参数的灵敏度分析

从灵敏度的计算结果可以得到：随着不透水面覆盖范围或渗透系数的增大，海岛地下淡水储量随之减少；当降雨入渗补给系数或给水度增大时，海岛地下淡水储量也随之增加。岛屿地下淡水储量对于给水度的变化最为敏感，其次是降雨入渗补给系数和岛屿介质渗透系数，再次是不透水面覆盖范围。降雨入渗补给系数和渗透系数的灵敏度在初始值对称增减时呈对称分布，且随变幅的增大灵敏度也随之增大。给水度和不透水面覆盖范围灵敏度随着其自身的增大而增大。