地貌变化对永兴岛淡水透镜体影响的数值模拟

盛 冲，许鹤华，张文涛

(1.中国科学院南海海洋研究所/中国科学院边缘海与大洋地质重点实验室，广东 广州 510301；2.中国科学院大学，北京 100049)

南海诸岛自古为中国领土，据统计，在地貌学上南海共拥有62座干出礁及49个灰沙岛，其地缘政治复杂，战略地位特殊，拥有丰富的海洋资源[1]。其中的永兴岛是西沙群岛中最大的灰沙岛屿，为三沙市政治、经济、文化中心，岛上常驻人口2 500多人。随着近年来三沙市布局规划的逐步开展，如何使岛上的环境更加适宜人们长期居住成为了关注的重点，其中备受关注是岛上淡水的补给问题。由于特殊的地质地形结构，岛上的淡水资源奇缺，这严重制约了日常生活和生态环境的改善。根据1996年中国科学院南沙综合科学考察队对南海岛礁的考察发现[2]，在南海某些岛屿上，部分降雨入渗补给含水层，并在补水和损失过程中维持着平衡状态，这部分浮在底层海水之上的淡化水体，即“淡水透镜体”(Freshwater Lens)。

在地貌未发生大规模改变前，永兴岛淡水透镜体的最大厚度约13.5 m，淡水资源储量约为147.2×104m3[3～4]。但随着岛上各项工程和规划的逐步开展，截止到2018年初，永兴岛的沙丘面积增加了近50%，新增区域主要集中在岛礁的北部，这一片陆域把石岛完全并入了永兴岛。为了探究新增部分沙丘面积对岛礁淡水透镜体的影响以及之后永兴岛淡水资源储量的变化，利用数值模拟进行预测显得十分必要。

本文在结合永兴岛地貌模型和前人研究的基础上，利用GMS软件构建了永兴岛淡水透镜体的三维变密度耦合模型，并对与岛礁淡水透镜体形成有关的部分参数进行了敏感性分析，探讨了地貌变化对永兴岛淡水透镜体的影响，此外还预测了永兴岛淡水透镜体再次达到稳定状态所需的时间及淡水资源储量。这对永兴岛及其他相似岛屿淡水资源的开发和管理具有重要意义。

图1 研究区示意图及观测点位置(据文[5]，有改动)Fig.1 Schematic diagram of the study area and location of the observation points (modified after[5])

1 研究区概况

永兴岛属于中国西沙群岛东部的宣德群岛，主要由珊瑚和其他生物骨骼堆积而成，四周礁坪环绕，向海部分经礁坪前缘以大角度进入深海盆地[5]。2009年前，永兴岛外形近似为椭圆平底盆地，东西长约1 950 m，南北宽约1 350 m，平均海拔高度约5 m。但随着岛上各项规划的逐步开展，目前永兴岛的本岛面积已由之前的2.13 km2增加到了3.2 km2(包括石岛)，岛屿的形状也发生巨大改变(图1)。此外永兴岛降雨量丰富，多年平均降雨量为1 509 mm；由于岛上土体的孔渗性很好且岛边缘有砂堤，降雨很少形成径流流失，除一部分被植物截留蒸发外，大部分降至地面入渗补给岛下淡水[6]。

目前关于永兴岛已有较为详细的水文地质钻孔资料，地质勘探结果表明[3, 7～8]，原永兴岛区域地表以下0～22 m主要为全新统碳酸盐生物骨壳堆积物，即珊瑚贝壳砂层，由于未发生成岩作用，呈松散状态，孔隙率和渗透率较高，其中藏有淡水透镜体；22～169 m主要由灰白色珊瑚礁灰岩组成，该地层虽然发生了成岩作用，但由于更新世时位于海平面以上，所以岩溶体系发育，地层总的孔渗性高，海水容易渗入(图2)。根据前人在永兴岛上进行的抽水实验结果[6]，不整合面以上含水层的渗透系数为70 m/d，之下的礁灰岩地层由于溶洞体系发育，渗透系数取500 m/d[9～10]；给水度为0.14，孔隙度为0.45。弥散度由于存在一定的尺度效应，且难以通过实验准确测得，在进行数值模拟计算时一般根据实验室测值，并结合模拟结果和经验参数进行确定[11～12]，因此在进行数值模拟时，弥散度的初始值：纵向弥散度取5 m，横向弥散度取0.5 m，垂向弥散度取0.05 m。

图2 研究区C-C’线水文地质剖面(据文[5，8]，有改动)Fig.2 Hydrogeologic profile of line C-C’ in the study area (modified after [5，8])

永兴岛的北部区域主要由近年来海浪及人为因素搬运堆积的珊瑚砂砾等生物碎屑组成，面积约1.2 km2，平均厚度达6 m。为了明确该区域珊瑚砂的水理特征，对该区域的几组珊瑚砂样进行颗粒分析，具体的颗粒级配曲线见图3。该区域珊瑚砂颗粒级配不良且不连续，此外样品的有效粒径d10为40.31 μm、特征粒径d30为85.18 μm、中值粒径d50为151.3 μm和控制粒径d60为208.2 μm。计算所得不均匀系数Cu为5.615和曲率系数Cc为0.865，因此定义为细粒土质砂，其中含有珊瑚碎片。为了进一步探究该区域珊瑚砂的其他水文地质特征，还在岛上进行了双环实验和抽水实验，相关的水文地质参数见表1。弥散度需结合室内试验和经验参数确定。

图3 岛礁北部珊瑚砂的颗粒级配曲线Fig.3 Grading curve of the coral sand in the northern part of the Yongxing Island

参数取值渗透系数K/(m·d-1)5 孔隙度n0.42给水度Sy0.103弥散度α纵向弥散度3 m,横向弥散度0.3 m,垂向弥散度0.03 m贮水率Ss1×10-5

2 数值模型的构建

2.1 水文地质概念模型

通过分析研究区地貌及水文地质条件，将永兴岛模拟区域的地层概化为三层，其中第一层为地表到海平面以下1 m，主要由北部人工堆积的珊瑚砂砾和原有岛屿未固结的全新统珊瑚砂组成，四周以海岸带为边界，概化为定水头和定浓度边界，面积约3.4 km2(图4)。考虑到淡化水体向海域延伸的可能，且岛礁在礁坪前缘以大角度向海底倾斜，因此模型其他两层的四周边界向外延伸至礁坪前缘[13～16]，同样概化为定水头和定浓度边界。其中第二层主要由全新统松散珊瑚砂构成，孔渗性良好，为淡水透镜体的主要赋存层；不整合面到海平面以下45 m处概化为模型的第三层，为更新统固结的礁灰岩地层，溶洞体系发育，海水较容易渗入。由于淡水透镜体悬浮于海水之上，通过孔隙、溶洞与海水相通，潮汐会引起淡水透镜体整体涨落，但不影响透镜体内部动力学特征，因此不考虑潮汐的影响[6, 8, 17]。

图4 水文地质概念模型示意图(C-C’剖面)Fig.4 Conceptual hydrogeological model of theYongxing Island (profile along C-C’)

模型的顶部为定流量边界，由于这部分补给量受降雨强度、植被、地形等多种因素影响，难以确定，因此参考一些经验数据及相关统计资料[6, 18～19]，取年均降雨量的45%均分至每天由地表补给淡水透镜体。模型底部位于海平面以下45 m处，与海水相连且不存在水量交换，因此可以视为零流量边界。

关于模型初始条件的概化，在淡水透镜体未形成之前，岛礁海平面以下区域完全被海水充填，随着降雨不断入渗补给，逐渐在岛上形成淡水。因此模型海平面以上的Cl-初始浓度取0 g/L，海平面以下的初始浓度取海水中的Cl-浓度19 g/L；模型初始水头值取0 m。地貌变化后的岛礁淡水透镜形成初始条件：将此前形成的稳定淡水透镜体的水头值和浓度值作为原永兴岛区域的初始条件；北部新增区域海平面上、下的Cl-初始浓度分别取0 g/L和19 g/L，初始水头值取0 m。

2.2 控制方程

在构建淡水透镜体三维数值模型时，由于海水与淡水属于可混溶性液体，对于这类问题需要考虑过渡带的存在以及对流-弥散作用，因此需要2个方程进行描述：水流方程，用来描述密度不断改变的液体流动；对流-弥散方程用来描述地下水中盐分的运移[20]。基于质量守恒和达西定律确定的变密度水流方程为：

(1)

式中：K——渗透系数/(m·d-1)；

H——测压管水头/m；

c——混合流体密度/(kg·m-3)；

Ss——贮水率；

t——时间/d；

ρ0——淡水密度(参考密度)/(kg·m-3)；

qs——单位体积多空介质的源或汇强度/d-1;

η——密度耦合系数/(kg·m-3)；

ρ——混合流体密度/(kg·m-3)。

在考虑流体密度时，忽略压力和温度对流体密度的影响，将流体密度处理成溶质浓度的线性函数，表示为：

(2)

(3)

η=ε/cs

(4)

式中：ρs——海水密度/(kg·m-3)；

cs——海水密度对应的质量浓度/(kg·m-3)；

ε——密度差相对比率。

溶质运移方程为：

(5)

式中：ui——渗透流速/(m·d-1)；

Dii——弥散系数/(m2·d-1)；

c*——源或汇的流体浓度/(kg·m-3)；

n——多孔介质孔隙度。

2.3 网格剖分及时间离散

模拟区域平面上的网格大小为32 m×26 m，垂向共分为30层，其中由于淡水透镜体主要存在于不整合面之上，因此对不整合面所处的位置及以上的网格进行细化，模型有效单元格共计294 171个。

模型总模拟时间24 000 d，其中初始时间步长为0.01 d，增益因子为2，最大运移步长为2 d，模型运行结果每隔365 d输出一次。

3 模型结果

3.1 永兴岛淡水透镜体的模拟结果与检验

利用GMS软件中的SEAWAT模块对模型进行求解，该模块由美国地调局(USGS)研发推出，是一款用于模拟多孔介质中三维非稳定变密度流的有限差分模型，其主要是通过同步时间步长方法，即交替使用修改过的MODFLOW和MT3DMS分别求解地下水流运动控制方程和溶质运移控制方程，实现水流和溶质运移的耦合求解[21]，目前已广泛应用于海水入侵、岛礁淡水透镜体等各种三维变密度地下水模拟问题的研究中。

模型的计算结果见图5和图6。永兴岛地貌变化前的淡水透镜体海平面以上的高度约0～0.41 m，其中最大水头值为0.41 m，正是由于这部分高出海平面的高度，保持了淡水不断地流向海洋，将海水向外推移，从而使淡水透镜体始终处于一个动态平衡状态。若以Cl-浓度0.6 g/L作为淡水透镜体的边界浓度[6]，则永兴岛淡水透镜体的最大厚度为14 m。

图5 淡水透镜体海平面以上的水头分布及3D形态Fig.5 Distribution of hydraulic heads above the sea level and the 3D form of the freshwater lens

图6 地貌变化前B-B’剖面的Cl-浓度分布Fig.6 Cl- concentration distribution of profile along B-B’ before the geomorphologic change

高密度电阻率法是以岩土体的电性差异为基础，研究在施加电场作用下，地下岩土体传导电流的变化分布规律，该方法能够很好地揭示岛上淡水透镜体的厚度与范围。Stewart详细论述了小海岛淡水透镜体电法勘探的方法原理及应用[22]，Comte利用数值法并结合电法勘探数据探讨了降雨补给对淡水透镜体的影响[13]。因此根据前人在岛上进行高密度电法勘探的研究结果[4]，结合6个观测点的淡水透镜体厚度对模型结果进行检验，具体观测点位置见图1，可以看出模拟结果基本符合永兴岛淡水透镜体的实际形态。

表2 模拟结果的识别与检验

3.2 参数的敏感性分析

模型的敏感性分析是为了对已经识别过的模型进行不确定性量化，是模型预测中不可缺少的环节。通常敏感性分析的做法是在同一时间内保证其他值不变的情况下，只改变一个参数值，观察其对模拟结果的影响。但更为有效的敏感性分析是计算敏感度，它表示一个因素的变化对模型计算结果的影响程度[23]。由于不同模型变量的量纲不同，不同参数的量纲差别也很大，为了便于不同参数间敏感度大小的比较，需要进一步无量纲化，通常表示为：

(6)

式中:βi,k——模型变量在观测点上的敏感度；

Hi——模型变量；

ak——第k个参数的取值。

其中，βi,k的绝对值越大，表明相应参数的变化对模拟结果的影响越大。考虑到研究区实际的水文地质情况，以及对淡水透镜体厚度的影响，选取渗透系数、降雨入渗补给系数和给水度进行敏感性分析。模型变量为淡水透镜体的最大厚度。由于淡水透镜体主要形成于不整合面以上未固结的珊瑚砂中，渗透系数一般为40～120 m/d，给水度为0.1～0.2[24]；降雨入渗补给系数难以确定，前人研究中海岛降雨入渗补给系数多取0.4～0.67[6, 18～19]。经过模型参数识别，以上参数的初始赋值分别为75 m/d，0.45和0.14，为了便于分析，确定各个参数在其初始值附近分别增加或减少10%、20%、30%、40%对模型进行敏感度分析，具体的敏感性分析结果见图7。从参数的敏感性分析结果可以看出：

(1)随着渗透系数的增大，淡水透镜体的厚度逐渐减小。当降雨入渗补给系数增大时，淡水透镜体的厚度也随之增大。给水度对淡水透镜体的厚度几乎没有影响。

(2)岛礁淡水透镜体的厚度对降雨入渗补给系数的变化最为敏感，其次是含水层的渗透系数，对含水层给水度的敏感度最小。

(3)各个参数在其初始赋值附近对称性增加或减少时，其敏感度大小呈现出对称性的变化规律，即随着参数的变化幅度增大时，敏感度系数也随之增大。

图7 参数的敏感性分析Fig.7 Sensitivity analysis of parameters(a)、(b) and (c) The maximum thickness of the freshwater lens varyingwith each parameter；(d) Sensitivity comparison of each parameter

图8 沙丘面积增加后永兴岛的淡水透镜体形成过程(a)沙丘面积增加前的淡水透镜体形态；(b)、(c)、(d)、(e)、(f)分别为永兴岛沙丘面积增加后1年、5年、10年、15年、25年的淡水透镜体形态Fig.8 Formation process of the freshwater lens in the Yongxing Island after increase in the dune area(a)state of the freshwater lens before the increase in the dune area; (b)、(c)、(d)、(e) and (f) State of the freshwater lens in 1, 5, 10, 15 and 25 years after the increase in the dune area in the Yongxing Island

3.3 地貌变化后永兴岛淡水透镜体的模拟结果

根据模型的计算结果，地貌变化后的永兴岛淡水透镜体随时间的变化见图8，可以看出在之前5年时间内，随着岛上降雨的不断入渗补给，永兴岛这部分新增区域地下水体的离子浓度不断得到淡化，但这一阶段还没有淡水体的产生，可以称之为淡水透镜体形成的准备阶段；之后随着时间的不断推移，包括新增区域在内，整个永兴岛淡水透镜体的厚度及储量都得到了增加。大约在25年后能够再次达到稳定状态。

再次稳定时的永兴岛淡水透镜体在形态上明显不对称，即永兴岛西南区域的淡水透镜体厚度较大，最大深度达到了17.5 m，而永兴岛东北区域淡水透镜体的最大厚度仅为8.5 m(图8)。初步分析认为，该区域的沙丘面积较小，且平均宽度仅为685 m，接受降雨入渗补给的面积有限，根据前面的敏感性分析可知，降雨的入渗补给对岛礁淡水透镜体厚度的影响最大；过窄的区域宽度会导致降雨入渗补给的淡水很快流入海洋，故而使得再次稳定时的永兴岛淡水透镜体产生了明显的不对称现象。

稳定时淡水透镜体的水头分布见图9。淡水透镜体中央处最大水头值为0.49 m，根据淡水透镜体的贮量计算公式[6]：

(7)

式中:V——贮水量/m3；

μ——含水层的给水度；

S——平面上每个差分网格的面积/m2；

hf,i,j——每个网格的淡水厚度/m。

图9 25年后永兴岛海平面以上的水头分布Fig.9 Distribution of hydraulic head abov the sea level in theYongxing Island 25 years later

计算出25年后的永兴岛淡水透镜体的贮量约为328×104m3。虽然岛屿面积仅增加了50%，但与地貌变化前的淡水储量相比，沙丘面积增加后的淡水储量增加了近一倍，这是因为随着岛礁面积的增大，不仅新增区域的海平面以下形成了淡水，原先的岛礁淡水透镜体厚度也随之增加。

4 结论

(1)永兴岛北部新增的这部分区域主要由人为因素搬运堆积的珊瑚砂砾等生物碎屑组成，具有明显的颗粒级配不良且不连续的特征，结构疏松且孔渗性良好，有利于岛礁淡水透镜体的形成。

(2)结合永兴岛的地貌地层条件，建立了永兴岛淡水透镜体的三维数值模型，计算得出地貌变化前的永兴岛淡水透镜体的最大厚度达14 m，其中海平面以上的高度为0～0.41 m，将模拟结果与岛上高密度电阻率法的观测数据进行对比拟合，模型能够很好地反映淡水透镜体的特征。

(3)对影响岛礁淡水透镜体形成的主要参数，如渗透系数，降雨的入渗补给系数，给水度进行了敏感性分析，可以看出，降雨入渗补给系数对岛礁淡水透镜体形成的影响最大，其次是含水层的渗透系数，而含水层的给水度对岛礁淡水透镜体的形成影响最小。

(4)在保持现有岛礁面积不变的情况下，永兴岛大约在25年后再次形成稳定的淡水透镜体，此时淡水透镜体的最大厚度为17.5 m，其中岛礁东北部区域淡化水体的最大厚度为8.5 m；计算所得的淡水资源储量增加了近一倍，约328×104m3。

(5)淡水透镜体再次稳定时，其东北区域的厚度较小，这主要是由于该区域的沙丘面积和宽度较小；进一步可以认为若沙丘面积和宽度过小，如面积增加前的永兴岛东北区域，降雨入渗补给的淡水会很快流入海洋，难以形成淡水透镜体。