基岩岛屿地下水流场数值模拟研究

许颢砾， 王大庆， 邓正栋， 丁志斌， 赵小兰， 刘志新， 许新刚， 苏合岩

(1.陆军工程大学国防工程学院，南京 210007； 2.中国矿业大学资源与地球科学学院，徐州 221116)

0 引言

岛屿的淡水资源十分宝贵，了解其地下水资源状况非常必要。美国较早开发利用了沙质岛屿的“淡水透镜体”[1-2]，并利用电法勘探确定其厚度和范围[3]。在20世纪90年代，我国海军和后勤学院研究解决西沙的供水问题时，对海岛“淡水透镜体”的有关理论开展过研究，并获得了丰硕的成果[4-5]。中国地质大学陈崇希教授和李国敏研究员运用水头、浓度相互依赖的有限元法对涠洲岛的火山碎屑孔隙含水层及玄武岩孔洞-裂隙含水层的水位、水流、水质等参数进行了公式计算和模拟分析[6]。中国科学院庞忠和研究员联合南京大学对我国东部沿海的庙岛群岛进行水文地质、物探等方面的研究，验证了基岩岛屿“淡水蘑菇体”的相关理论[7-8]。随着计算机模拟技术的发展，MODFLOW程序、Visual Modflow软件[9]、GMS软件[10]、FEFLOW软件[11]和TOUGH软件[12]等模拟系统相继出现，地下水模拟变得简单化、流程化。在模拟过程中，需要对基岩岛屿的地下水补径排及海水入侵情况进行详细调查，通过试验获得研究区地质参数，并利用软件建立模型进行地下水的模拟[13]。针对基岩岛屿的模拟，周鹏鹏等[14-15]根据有限的水井数据和降雨蒸发数据，采用Visual Modflow软件对湛江东海岛及其周边地区(面积约286 km2)进行了地下水模拟。中国地质大学温汉辉等[16]对雷州半岛(约13 225 km2)的水文、地质(主要为岩浆岩和变质岩)、地貌、气象等信息进行了较为详细的调查，利用GMS软件对该半岛进行了地下水数值模拟，并对地下水的合理开发提出了宝贵建议。上述研究是对较大尺度(几百到上万平方公里)范围的基岩岛屿地下水流场进行模拟，基于达西定律将基岩岛屿概化为孔隙介质从而得到地下水流场的大体分布。显然，这些模拟工作与实际情况有较大的出入，基岩岛屿的地下水类型多为裂隙水，仅表层的风化层为孔隙介质(除深层非基岩地下水)，且基岩岛屿的裂隙水分布不均匀。针对基岩岛屿地下水模拟的诸多难点，本文以珠海外伶仃岛为例，利用遥感技术、物探技术、水文地质方法及地下水模拟技术来研究基岩岛屿的地下水流场，改善了因地质条件复杂、水井数据等初始条件较少而难以全面对基岩岛屿地下水状况进行模拟的现状，对同类型岛屿的地下水流场模拟工作有参考意义。

1 研究区概况

以基岩岛屿珠海外伶仃岛为对象展开研究。外伶仃岛位于广东省珠海市万山海洋开发试验区的东北部, 西距万山区桂山岛17.6 n mile(1n mile=1 852 m)，离珠海市香洲区27.5 n mile，北距香港长洲6 n mile，距香港九龙尖沙咀11 n mile。其经纬度范围是114°1′11.19″～114°3′15.40″E、22°5′17.16″～22°6′57.54″N，具体位置如图1所示。

图1 研究区位置(外伶仃岛)

外伶仃岛地区岩石为燕山早期第三期花岗岩，一般呈基岩、岩柱产出，主要由黑云母花岗岩组成，属查氏分类第II类3科硅酸过饱和碱性岩石。同位素年龄为160～138 Ma，主要集中在155～140 Ma，归属于晚侏罗世。岛屿面积约为4. 33 km2，岛屿周围的海岸线长12.3 km。整个岛屿长约3 200 m，宽约2 400 m，岛屿地势呈西北中部、北部高，西北近海岸、南近海岸低，主峰外伶仃峰海拔311.8 m，岛中央与岛东南部之间存在构造形成的沟谷。岛上共有3处沙滩，海岸线多为基岩海岸，也有悬崖峭壁。岛内地下水以基岩风化裂隙水为主，也存在构造裂隙水[17]。

外伶仃岛位居北回归线以南，属亚热带季风气候，日照充足。1月份平均气温14.8 ℃，7月份平均气温27.9 ℃，3—4月份为雾季，3月份雾最多，5—7月份天气平和，8—9月份有台风，最大风速12级以上，10—12 月份天气平和。四季温差不大，冬天无严寒，夏天无酷暑。大万山气象站1973—1991年的气温资料表明，有85%的年份冬季极端低温不低于4 ℃，大多数年份夏季极端高温不超过35 ℃。外伶仃岛降雨量一年中有丰枯季之分，4—9月份为丰水期，降雨量占年降雨量的80%以上，月平均降雨量在200～350 mm之间； 11月份至次年2月份为枯水期，降雨量较少，月降雨量只有10～20 mm。全年常向风为东风，夏季风以南风为主，冬季风以东北风为主。

外伶仃岛除岛南的南湾部分水深小于10 m，其他海湾海水深都大于10 m，且在海岛的东南海域，其水深大于20 m。外伶仃岛周边区域海水pH 值为8.0～8.2； 海水温度年变化与太阳辐射量有关，每年3月开始升温，7—8月水温最高，9月开始降温，1—2月水温最低，海水温度年平均为22.7 ℃； 海水盐度1月份为33.19‰，7月份为23.59‰，平均为28.39‰。全年之中，冬季的波浪高于夏季，一般波浪不高，平均为0.9～1.9 m，受台风直接影响时，平均波高5～6 m，最大波高8 m。受强台风直接影响时，平均波高7～8 m，最大波高大于10 m。外伶仃海域的潮汐属不正规半日潮，最高潮位为3.4 m，最低潮位为0.84 m，最大潮差为2.52 m。

2 基岩岛屿地下水流场模拟

2.1 地下水水位获取方法

本文重点关注岛屿地表以下第1层含水层(基岩岛屿地下的风化裂隙水或构造裂隙水含水层)，多数情况为潜水层，有时也指上层滞水或承压水。利用地下水遥感评估法[18-19]所获得的第一层含水层的地下水遥感评估分值(S)与水井实测水量呈正相关线性关系，再考虑水的单位面积，并将井或泉位置的水量取平均(目的是将水量统一以单位投影面积来计算水位埋深)，计算公式为

(1)

式中：h为单位水位埋深，m； Δs为单位面积，m2；S为第一层含水层的地下水遥感评估值，m3； a、b为常数。

将S与h进行拟合，公式为

S=a1·eb1hΔs

，

(2)

式中：h为单位水位埋深，m； Δs为单位面积，m2；S为第一层含水层的地下水遥感评估值，m3；a1、b1为常系数。

求出系数a1、b1，可计算整个研究区第一层含水层的地下水水位埋深(H)，其公式为

H=D-h

，

(3)

式中：D为数字高程模型数据的值，m；H为第1层含水层地下水水位值，m；h为地下水埋深值，m。

将第1层含水层地下水水位H矩阵导入地下水模拟软件或程序进行地下水流场模拟的初始水位值，可以克服基岩岛屿难以全面获得第一层含水层的地下水的初始水位标高的困难。

2.2 岛屿勘察资料

珠江口外海岛的地质地貌结构简单，多为花岗岩。海岛岩层自上而下一般可分为3层，分别为地表风化层土壤(植被覆盖)和裸露岩石、半风化层以及花岗岩基岩(图2)。图2(a)是位于岛屿北侧的一处天然水池，常年有水，为泉水出露和降雨补给形成，整个岛有多处类似的天然水池(或称水坑)，图2(b)可以较好地反映地表1～2 m的松散风化层，图2(c)为半风化层，图2(d)为裸露的花岗岩石蛋，图2(e)为花岗岩基岩底层。根据外伶仃岛的有降雨补给和地下水补给形成的地表水域及以前人们使用的老井及泉水出露，已经可以说明岛屿存在地下淡水，但是容量有限。岛屿周边都是海水，地下水的唯一补给源就是雨水。对于海岛的供水系统而言，以前人口少时，主要利用水井地下水或积累雨水的水池。随着旅游业发展及当地政府部门的开发，海岛人口逐步增多，原有的供水体系已无法满足需要，于是就用船运来的水通过高压泵和管道输送到各个用水地点(如外伶仃岛)，或者是在当地建立小型水库来满足岛民的生活用水(桂山岛和东澳岛)。

(a) 天然水池 (b) 地表松散风化层

2.3 水文地质参数

对岛屿花岗岩取样测试，并根据《水文地质手册》[20]、《堤防工程手册》[21]和《地下水污染物迁移模拟》[9]所给出的水文地质参数检验值，将水平渗透系数、竖直渗透系数、给水度、储水系数和各向异性参数等水文地质参数输入地下水模拟软件或程序(表1)进行模拟计算。

该地区的平均降雨量约为0.005 33 m/d(数据来自珠海气象局)，由于该基岩岛屿的特殊性，在基岩海岸边缘是以淡水“蘑菇体”的形式存在，根据地下水模拟技术国内外研究现状[13]，可以将基岩海岸线设成定为模拟的范围边界。

表1 各种岩层的水文地质参数

2.4 岛屿地质体分层情况

本次采用美国GSSI探地雷达SIR-4000仪器探测岛屿地层分层情况。本文利用100 MHz天线，可以对地下8～20 m的范围进行探测。为佐证探地雷达的探测结果，采用法国SYSCALR2直流电法仪进行了三极装置直流电法探测。如图3所示，沿岛屿西南部的一条公路进行了约1.5 km的探地雷达作业与直流电法作业(图3)。提取同一段水平距离220 m测线的探地雷达探测结果(图4)和直流电法探测结果(图5)进行对比。图4、图5按沿测线的地势起伏将探测结果绘制在地表之下，以反映地下探测深度范围内的地球物理分层。如图4，根据地质体介电常数的差异可以较清晰地看出： 厚度1～2 m的地表层与图2所示的地表岩体风化层对应； 厚度5～7 m是另外一个地质层，可以对应于半风化层； 厚度5～7 m以下的部分可以对应于花岗岩基岩底层。图5可反映如下情况： 低视电阻率的蓝色区域厚度为1～2 m； 中等大小视电阻率(绿色及黄色)厚度为2～5 m； 5～7 m深度之下的高视电阻率(橙色和红色)对应花岗岩基岩底层。可见，探地雷达与直流电法的探测结果符合较好，与现场勘查情况对应也很好，可以确定岛屿在地下水范围内的地质分层情况为： 地表之下1～2 m为地表风化层； 地表下5～7 m为半风化层； 地表7 m以下为花岗岩基岩底层。

图3 外伶仃岛西南部物探测线位置

图4 岛屿西南部的探地雷达探测220 m水平距离结果

图5 岛屿西南部的直流电法探测220 m水平距离结果

2.5 地下水水位分布参数

根据研究区地下水富集性评价结果等级图，可得出每一点地下水富集性遥感评估值(S)，再根据平均水量直接用该点实际单位水位埋深h来表示，即S=a1·eb1h。S与h进行拟合，得出系数a1、b1，再利用ENVI软件计算整个研究区，或研究区的关键点的h，将地下水的水位H矩阵作为该地区第一层含水层的初始水位数据。表2是地下水评估得分与单位水位埋深关系表。

表2 地下水评估得分与单位水位埋深关系[17]

根据井或泉观测值的平均出水量得到单位时间单位面积的水位埋深，将h与S进行曲线拟合(图6)，得出拟合曲线关系式为

S=0.427e0.239h，

(4)

式中：h为单位水位埋深，m；S为地下水遥感评估得分，m3； 0.427和0.239为系数。

此时，拟合优度R2=0.847 2(显著性P<0.05)，说明井或泉中对应的单位时间单位面积的水位埋深h与地下水富集性遥感评估得分S的拟合程度较好。

图6 水位埋深h与评估得分S的拟合曲线

根据公式(4)反推出各个点的评估值S对应h的计算关系式为

(5)

式中：h为单位水位埋深，m；S为地下水遥感评估值，m3； 0.427和0.239为系数。

根据式(3)和式(5)，利用ENVI软件将地下水遥感评估得分S值计算整个岛的初始水位H值，即基岩岛屿的地下水流场模拟的初始水位用H值表示。

3 研究区地下水模型建模

3.1 模型建立步骤

建立研究区地下水模型分为4个步骤。

(1)根据研究区的水位地质资料背景和实际勘测岛屿资料，设置如表1的水文地质参数。

(2)导入底图，再导入该地区的数字高程模型(DEM)数据的TIF文件，并绘制相应的CAD底图，方便该地区的地表水域的区域、岛屿区域、分块分区及边界域。

(3)将该地区的降雨量设置为正值，蒸发量设置为负值，则补给为二者之间的差值，即降雨蒸发量为+0.005 33 m/d，在有水坑、小型湖泊的位置设置补给为0.006 03 m/d，并设置相应井的抽水量。

(4)将DEM数据与H值的TIF导入GMS，并转化为散点(图7)。DEM数据的利用，相对于传统的人工跑点用GPS点位而言，可以提高研究区地质体的概化精度。

图7 DEM的TIF文件和H的TIF文件在GMS上转化为散点

3.2 网格剖分与模型计算

在GMS软件上进行网格剖分，本模拟将该岛屿剖分为3层，分别是地表松散风化层、第1层(含水层)和第2层(隔水层)，本文主要针对地下水的流场分布，再利用GMS自带的Check Simulation功能进行自动查错，修正后进行MODFLOW模块计算。

3.3 地下水流场模拟结果分析

地表松散风化层(1～2 m)水位及流场模拟结果如图8，第1层(含水层)(2～5 m)地下水水位及流场模拟结果如图9。

图8与图9经比较可以看出，地表松散风化层中的水位与第1层含水层的地下水水位有明显的差异，虽然二者的地下水流场的方向趋势大致相同，但表面松散风化层的地下水分布很不均匀，且较大范围的风化层处于无水状态，而地下水遍布全岛，第1层含水层的地下水流场的水位分布较地表风化层的地下水流场更明显。从图9反映出的第一层含水层地下水水位模拟结果来看，水位角度呈中部和东南部高、周边低的趋势。从岛屿第一层含水层的地下水流场趋势上看，从中部向外流，从东南部向外流，周边又有向岛内流的趋势，应该是有海水的存在的原因。

图8 地表松散风化层(1～2 m)水流场模拟结果

图9 第一层(含水层)(2～5 m)地下水流场模拟结果

4 结果验证与分析

4.1 直流电法剖面实测数据

本研究利用直流电法仪以对称四极装置形式，在供电极距AB=30 m条件下进行实地探测。直流电法的2条测线具体位置如图10中的红色测线位置，第一条测线(4个点，约12 m长)位于岛屿南部水坑上方公路旁边，起始点坐标为22.096 674°N、114.038 623°E，海拔33.2 m。由图11(a)直流电法探测剖面图揭示的视电阻率差异可以推出，该剖面中部的地下5～13 m范围应该是有基岩裂隙水储存，且在该地区下方约10 m处出现泉水出露，泉位置为22.097 272°N、114.038 793°E，海拔21.4 m。第二条测线(6个点，约30 m长)位于岛屿南部水坑上方公路旁边，其起始点坐标为22.096 674°N、114.038 623°E，海拔33.2 m，结果如图11(b)所示。

图10 研究区南部直流电法作业测线位置

(a) 测线1

根据电阻率差异，可以推测出该剖面中部，地下5～13 m范围应该是有基岩裂隙水储存，且在该地区下方约10 m处出现泉水出露。泉位置为22.095 464°N、114.038 102°E，海拔8.8 m。

将直流电法探测解译的水位值与地下水流场模拟值(图9)进行对比，如表3所示。

表3 探测水位与模拟水位的对应值

4.2 相关性分析

图12所示为表3中同序号水位模拟值(横坐标)与实际水位值(纵坐标)散点及其拟合关系曲线，可见两者呈线性正相关，且R2=0.872 2(P<0.05)，相关性很好，说明模拟水位与实际水位结果符合很好。因此，将基于地下水遥感评估所得出的初始水位导入GMS软件并加设相应的岛屿水文地质参数所建立的流场模拟模型能够反映该岛屿的地下水流场分布情况。

图12 模拟水位与实际探测水位值的相关性分析

5 结论

本文针对基岩岛屿地下水流场参数难以全面获取的现状，采用遥感、物探和水文地质综合方法与技术对外伶仃岛地下水流场进行了模拟研究，得到了以下认识和结论。

(1)利用数字高程模型数据(普遍是30 m精度的数据，也可利用15 m或10 m等更高精度的数据)，并转化好格式导入地下水模拟软件或程序中，这样提高了地质体的模型概化精度。

(2)利用物探技术(探地雷达法与直流电法相互佐证)对岛屿地层进行了探测，获取了各地层厚度及分界面等数据，并结合水文地质试验参数(抽水试验参数、不同岩性的经验参数)来合理地设定地下水模拟的各地层参数。

(3)基于地下水遥感的评估得分与该地区的各点的水位(井与泉水域范围的探测值)进行拟合，再利用拟合关系式推出各个点的初始水位，并将所得的初始水位导入地下水模拟软件或程序中作为地下水的初始水位，且设定相应的初始条件与边界条件进行流场模拟。该方法克服了基岩岛屿的初始水位难以全面获取的困难，且模拟出的水位模拟值与探测水位值呈正相关，其两者的拟合精度也满足要求。

(4)遥感、物探和水文地质综合方法与技术可以较好地对基岩岛屿地下水流场进行模拟。但是本研究仍然是利用现有的地下水模拟软件基于等效理论来模拟基岩岛的裂隙水(除松散层孔隙水等以外)，将来希望可以尝试编程实现基于立方定律等理论的裂隙水的直接模拟。