許顥礫， 王大慶， 鄧正棟， 丁志斌， 趙小蘭， 劉志新， 許新剛， 蘇合巖

(1.陸軍工程大學(xué)國防工程學(xué)院，南京 210007； 2.中國礦業(yè)大學(xué)資源與地球科學(xué)學(xué)院，徐州 221116)

0 引言

島嶼的淡水資源十分寶貴，了解其地下水資源狀況非常必要。美國較早開發(fā)利用了沙質(zhì)島嶼的“淡水透鏡體”[1-2]，并利用電法勘探確定其厚度和范圍[3]。在20世紀(jì)90年代，我國海軍和后勤學(xué)院研究解決西沙的供水問題時，對海島“淡水透鏡體”的有關(guān)理論開展過研究，并獲得了豐碩的成果[4-5]。中國地質(zhì)大學(xué)陳崇希教授和李國敏研究員運用水頭、濃度相互依賴的有限元法對潿洲島的火山碎屑孔隙含水層及玄武巖孔洞-裂隙含水層的水位、水流、水質(zhì)等參數(shù)進(jìn)行了公式計算和模擬分析[6]。中國科學(xué)院龐忠和研究員聯(lián)合南京大學(xué)對我國東部沿海的廟島群島進(jìn)行水文地質(zhì)、物探等方面的研究，驗證了基巖島嶼“淡水蘑菇體”的相關(guān)理論[7-8]。隨著計算機(jī)模擬技術(shù)的發(fā)展，MODFLOW程序、Visual Modflow軟件[9]、GMS軟件[10]、FEFLOW軟件[11]和TOUGH軟件[12]等模擬系統(tǒng)相繼出現(xiàn)，地下水模擬變得簡單化、流程化。在模擬過程中，需要對基巖島嶼的地下水補徑排及海水入侵情況進(jìn)行詳細(xì)調(diào)查，通過試驗獲得研究區(qū)地質(zhì)參數(shù)，并利用軟件建立模型進(jìn)行地下水的模擬[13]。針對基巖島嶼的模擬，周鵬鵬等[14-15]根據(jù)有限的水井?dāng)?shù)據(jù)和降雨蒸發(fā)數(shù)據(jù)，采用Visual Modflow軟件對湛江東海島及其周邊地區(qū)(面積約286 km2)進(jìn)行了地下水模擬。中國地質(zhì)大學(xué)溫漢輝等[16]對雷州半島(約13 225 km2)的水文、地質(zhì)(主要為巖漿巖和變質(zhì)巖)、地貌、氣象等信息進(jìn)行了較為詳細(xì)的調(diào)查，利用GMS軟件對該半島進(jìn)行了地下水?dāng)?shù)值模擬，并對地下水的合理開發(fā)提出了寶貴建議。上述研究是對較大尺度(幾百到上萬平方公里)范圍的基巖島嶼地下水流場進(jìn)行模擬，基于達(dá)西定律將基巖島嶼概化為孔隙介質(zhì)從而得到地下水流場的大體分布。顯然，這些模擬工作與實際情況有較大的出入，基巖島嶼的地下水類型多為裂隙水，僅表層的風(fēng)化層為孔隙介質(zhì)(除深層非基巖地下水)，且基巖島嶼的裂隙水分布不均勻。針對基巖島嶼地下水模擬的諸多難點，本文以珠海外伶仃島為例，利用遙感技術(shù)、物探技術(shù)、水文地質(zhì)方法及地下水模擬技術(shù)來研究基巖島嶼的地下水流場，改善了因地質(zhì)條件復(fù)雜、水井?dāng)?shù)據(jù)等初始條件較少而難以全面對基巖島嶼地下水狀況進(jìn)行模擬的現(xiàn)狀，對同類型島嶼的地下水流場模擬工作有參考意義。

1 研究區(qū)概況

以基巖島嶼珠海外伶仃島為對象展開研究。外伶仃島位于廣東省珠海市萬山海洋開發(fā)試驗區(qū)的東北部, 西距萬山區(qū)桂山島17.6 n mile(1n mile=1 852 m)，離珠海市香洲區(qū)27.5 n mile，北距香港長洲6 n mile，距香港九龍尖沙咀11 n mile。其經(jīng)緯度范圍是114°1′11.19″～114°3′15.40″E、22°5′17.16″～22°6′57.54″N，具體位置如圖1所示。

圖1 研究區(qū)位置(外伶仃島)

外伶仃島地區(qū)巖石為燕山早期第三期花崗巖，一般呈基巖、巖柱產(chǎn)出，主要由黑云母花崗巖組成，屬查氏分類第II類3科硅酸過飽和堿性巖石。同位素年齡為160～138 Ma，主要集中在155～140 Ma，歸屬于晚侏羅世。島嶼面積約為4. 33 km2，島嶼周圍的海岸線長12.3 km。整個島嶼長約3 200 m，寬約2 400 m，島嶼地勢呈西北中部、北部高，西北近海岸、南近海岸低，主峰外伶仃峰海拔311.8 m，島中央與島東南部之間存在構(gòu)造形成的溝谷。島上共有3處沙灘，海岸線多為基巖海岸，也有懸崖峭壁。島內(nèi)地下水以基巖風(fēng)化裂隙水為主，也存在構(gòu)造裂隙水[17]。

外伶仃島位居北回歸線以南，屬亞熱帶季風(fēng)氣候，日照充足。1月份平均氣溫14.8 ℃，7月份平均氣溫27.9 ℃，3—4月份為霧季，3月份霧最多，5—7月份天氣平和，8—9月份有臺風(fēng)，最大風(fēng)速12級以上，10—12 月份天氣平和。四季溫差不大，冬天無嚴(yán)寒，夏天無酷暑。大萬山氣象站1973—1991年的氣溫資料表明，有85%的年份冬季極端低溫不低于4 ℃，大多數(shù)年份夏季極端高溫不超過35 ℃。外伶仃島降雨量一年中有豐枯季之分，4—9月份為豐水期，降雨量占年降雨量的80%以上，月平均降雨量在200～350 mm之間； 11月份至次年2月份為枯水期，降雨量較少，月降雨量只有10～20 mm。全年常向風(fēng)為東風(fēng)，夏季風(fēng)以南風(fēng)為主，冬季風(fēng)以東北風(fēng)為主。

外伶仃島除島南的南灣部分水深小于10 m，其他海灣海水深都大于10 m，且在海島的東南海域，其水深大于20 m。外伶仃島周邊區(qū)域海水pH 值為8.0～8.2； 海水溫度年變化與太陽輻射量有關(guān)，每年3月開始升溫，7—8月水溫最高，9月開始降溫，1—2月水溫最低，海水溫度年平均為22.7 ℃； 海水鹽度1月份為33.19‰，7月份為23.59‰，平均為28.39‰。全年之中，冬季的波浪高于夏季，一般波浪不高，平均為0.9～1.9 m，受臺風(fēng)直接影響時，平均波高5～6 m，最大波高8 m。受強(qiáng)臺風(fēng)直接影響時，平均波高7～8 m，最大波高大于10 m。外伶仃海域的潮汐屬不正規(guī)半日潮，最高潮位為3.4 m，最低潮位為0.84 m，最大潮差為2.52 m。

2 基巖島嶼地下水流場模擬

2.1 地下水水位獲取方法

本文重點關(guān)注島嶼地表以下第1層含水層(基巖島嶼地下的風(fēng)化裂隙水或構(gòu)造裂隙水含水層)，多數(shù)情況為潛水層，有時也指上層滯水或承壓水。利用地下水遙感評估法[18-19]所獲得的第一層含水層的地下水遙感評估分值(S)與水井實測水量呈正相關(guān)線性關(guān)系，再考慮水的單位面積，并將井或泉位置的水量取平均(目的是將水量統(tǒng)一以單位投影面積來計算水位埋深)，計算公式為

(1)

式中：h為單位水位埋深，m； Δs為單位面積，m2；S為第一層含水層的地下水遙感評估值，m3； a、b為常數(shù)。

將S與h進(jìn)行擬合，公式為

S=a1·eb1hΔs

，

(2)

式中：h為單位水位埋深，m； Δs為單位面積，m2；S為第一層含水層的地下水遙感評估值，m3；a1、b1為常系數(shù)。

求出系數(shù)a1、b1，可計算整個研究區(qū)第一層含水層的地下水水位埋深(H)，其公式為

H=D-h

，

(3)

式中：D為數(shù)字高程模型數(shù)據(jù)的值，m；H為第1層含水層地下水水位值，m；h為地下水埋深值，m。

將第1層含水層地下水水位H矩陣導(dǎo)入地下水模擬軟件或程序進(jìn)行地下水流場模擬的初始水位值，可以克服基巖島嶼難以全面獲得第一層含水層的地下水的初始水位標(biāo)高的困難。

2.2 島嶼勘察資料

珠江口外海島的地質(zhì)地貌結(jié)構(gòu)簡單，多為花崗巖。海島巖層自上而下一般可分為3層，分別為地表風(fēng)化層土壤(植被覆蓋)和裸露巖石、半風(fēng)化層以及花崗巖基巖(圖2)。圖2(a)是位于島嶼北側(cè)的一處天然水池，常年有水，為泉水出露和降雨補給形成，整個島有多處類似的天然水池(或稱水坑)，圖2(b)可以較好地反映地表1～2 m的松散風(fēng)化層，圖2(c)為半風(fēng)化層，圖2(d)為裸露的花崗巖石蛋，圖2(e)為花崗巖基巖底層。根據(jù)外伶仃島的有降雨補給和地下水補給形成的地表水域及以前人們使用的老井及泉水出露，已經(jīng)可以說明島嶼存在地下淡水，但是容量有限。島嶼周邊都是海水，地下水的唯一補給源就是雨水。對于海島的供水系統(tǒng)而言，以前人口少時，主要利用水井地下水或積累雨水的水池。隨著旅游業(yè)發(fā)展及當(dāng)?shù)卣块T的開發(fā)，海島人口逐步增多，原有的供水體系已無法滿足需要，于是就用船運來的水通過高壓泵和管道輸送到各個用水地點(如外伶仃島)，或者是在當(dāng)?shù)亟⑿⌒退畮靵頋M足島民的生活用水(桂山島和東澳島)。

(a) 天然水池 (b) 地表松散風(fēng)化層

2.3 水文地質(zhì)參數(shù)

對島嶼花崗巖取樣測試，并根據(jù)《水文地質(zhì)手冊》[20]、《堤防工程手冊》[21]和《地下水污染物遷移模擬》[9]所給出的水文地質(zhì)參數(shù)檢驗值，將水平滲透系數(shù)、豎直滲透系數(shù)、給水度、儲水系數(shù)和各向異性參數(shù)等水文地質(zhì)參數(shù)輸入地下水模擬軟件或程序(表1)進(jìn)行模擬計算。

該地區(qū)的平均降雨量約為0.005 33 m/d(數(shù)據(jù)來自珠海氣象局)，由于該基巖島嶼的特殊性，在基巖海岸邊緣是以淡水“蘑菇體”的形式存在，根據(jù)地下水模擬技術(shù)國內(nèi)外研究現(xiàn)狀[13]，可以將基巖海岸線設(shè)成定為模擬的范圍邊界。

表1 各種巖層的水文地質(zhì)參數(shù)

2.4 島嶼地質(zhì)體分層情況

本次采用美國GSSI探地雷達(dá)SIR-4000儀器探測島嶼地層分層情況。本文利用100 MHz天線，可以對地下8～20 m的范圍進(jìn)行探測。為佐證探地雷達(dá)的探測結(jié)果，采用法國SYSCALR2直流電法儀進(jìn)行了三極裝置直流電法探測。如圖3所示，沿島嶼西南部的一條公路進(jìn)行了約1.5 km的探地雷達(dá)作業(yè)與直流電法作業(yè)(圖3)。提取同一段水平距離220 m測線的探地雷達(dá)探測結(jié)果(圖4)和直流電法探測結(jié)果(圖5)進(jìn)行對比。圖4、圖5按沿測線的地勢起伏將探測結(jié)果繪制在地表之下，以反映地下探測深度范圍內(nèi)的地球物理分層。如圖4，根據(jù)地質(zhì)體介電常數(shù)的差異可以較清晰地看出： 厚度1～2 m的地表層與圖2所示的地表巖體風(fēng)化層對應(yīng)； 厚度5～7 m是另外一個地質(zhì)層，可以對應(yīng)于半風(fēng)化層； 厚度5～7 m以下的部分可以對應(yīng)于花崗巖基巖底層。圖5可反映如下情況： 低視電阻率的藍(lán)色區(qū)域厚度為1～2 m； 中等大小視電阻率(綠色及黃色)厚度為2～5 m； 5～7 m深度之下的高視電阻率(橙色和紅色)對應(yīng)花崗巖基巖底層?？梢姡降乩走_(dá)與直流電法的探測結(jié)果符合較好，與現(xiàn)場勘查情況對應(yīng)也很好，可以確定島嶼在地下水范圍內(nèi)的地質(zhì)分層情況為： 地表之下1～2 m為地表風(fēng)化層； 地表下5～7 m為半風(fēng)化層； 地表7 m以下為花崗巖基巖底層。

圖3 外伶仃島西南部物探測線位置

圖4 島嶼西南部的探地雷達(dá)探測220 m水平距離結(jié)果

圖5 島嶼西南部的直流電法探測220 m水平距離結(jié)果

2.5 地下水水位分布參數(shù)

根據(jù)研究區(qū)地下水富集性評價結(jié)果等級圖，可得出每一點地下水富集性遙感評估值(S)，再根據(jù)平均水量直接用該點實際單位水位埋深h來表示，即S=a1·eb1h。S與h進(jìn)行擬合，得出系數(shù)a1、b1，再利用ENVI軟件計算整個研究區(qū)，或研究區(qū)的關(guān)鍵點的h，將地下水的水位H矩陣作為該地區(qū)第一層含水層的初始水位數(shù)據(jù)。表2是地下水評估得分與單位水位埋深關(guān)系表。

表2 地下水評估得分與單位水位埋深關(guān)系[17]

根據(jù)井或泉觀測值的平均出水量得到單位時間單位面積的水位埋深，將h與S進(jìn)行曲線擬合(圖6)，得出擬合曲線關(guān)系式為

S=0.427e0.239h，

(4)

式中：h為單位水位埋深，m；S為地下水遙感評估得分，m3； 0.427和0.239為系數(shù)。

此時，擬合優(yōu)度R2=0.847 2(顯著性P<0.05)，說明井或泉中對應(yīng)的單位時間單位面積的水位埋深h與地下水富集性遙感評估得分S的擬合程度較好。

圖6 水位埋深h與評估得分S的擬合曲線

根據(jù)公式(4)反推出各個點的評估值S對應(yīng)h的計算關(guān)系式為

(5)

式中：h為單位水位埋深，m；S為地下水遙感評估值，m3； 0.427和0.239為系數(shù)。

根據(jù)式(3)和式(5)，利用ENVI軟件將地下水遙感評估得分S值計算整個島的初始水位H值，即基巖島嶼的地下水流場模擬的初始水位用H值表示。

3 研究區(qū)地下水模型建模

3.1 模型建立步驟

建立研究區(qū)地下水模型分為4個步驟。

(1)根據(jù)研究區(qū)的水位地質(zhì)資料背景和實際勘測島嶼資料，設(shè)置如表1的水文地質(zhì)參數(shù)。

(2)導(dǎo)入底圖，再導(dǎo)入該地區(qū)的數(shù)字高程模型(DEM)數(shù)據(jù)的TIF文件，并繪制相應(yīng)的CAD底圖，方便該地區(qū)的地表水域的區(qū)域、島嶼區(qū)域、分塊分區(qū)及邊界域。

(3)將該地區(qū)的降雨量設(shè)置為正值，蒸發(fā)量設(shè)置為負(fù)值，則補給為二者之間的差值，即降雨蒸發(fā)量為+0.005 33 m/d，在有水坑、小型湖泊的位置設(shè)置補給為0.006 03 m/d，并設(shè)置相應(yīng)井的抽水量。

(4)將DEM數(shù)據(jù)與H值的TIF導(dǎo)入GMS，并轉(zhuǎn)化為散點(圖7)。DEM數(shù)據(jù)的利用，相對于傳統(tǒng)的人工跑點用GPS點位而言，可以提高研究區(qū)地質(zhì)體的概化精度。

圖7 DEM的TIF文件和H的TIF文件在GMS上轉(zhuǎn)化為散點

3.2 網(wǎng)格剖分與模型計算

在GMS軟件上進(jìn)行網(wǎng)格剖分，本模擬將該島嶼剖分為3層，分別是地表松散風(fēng)化層、第1層(含水層)和第2層(隔水層)，本文主要針對地下水的流場分布，再利用GMS自帶的Check Simulation功能進(jìn)行自動查錯，修正后進(jìn)行MODFLOW模塊計算。

3.3 地下水流場模擬結(jié)果分析

地表松散風(fēng)化層(1～2 m)水位及流場模擬結(jié)果如圖8，第1層(含水層)(2～5 m)地下水水位及流場模擬結(jié)果如圖9。

圖8與圖9經(jīng)比較可以看出，地表松散風(fēng)化層中的水位與第1層含水層的地下水水位有明顯的差異，雖然二者的地下水流場的方向趨勢大致相同，但表面松散風(fēng)化層的地下水分布很不均勻，且較大范圍的風(fēng)化層處于無水狀態(tài)，而地下水遍布全島，第1層含水層的地下水流場的水位分布較地表風(fēng)化層的地下水流場更明顯。從圖9反映出的第一層含水層地下水水位模擬結(jié)果來看，水位角度呈中部和東南部高、周邊低的趨勢。從島嶼第一層含水層的地下水流場趨勢上看，從中部向外流，從東南部向外流，周邊又有向島內(nèi)流的趨勢，應(yīng)該是有海水的存在的原因。

圖8 地表松散風(fēng)化層(1～2 m)水流場模擬結(jié)果

圖9 第一層(含水層)(2～5 m)地下水流場模擬結(jié)果

4 結(jié)果驗證與分析

4.1 直流電法剖面實測數(shù)據(jù)

本研究利用直流電法儀以對稱四極裝置形式，在供電極距AB=30 m條件下進(jìn)行實地探測。直流電法的2條測線具體位置如圖10中的紅色測線位置，第一條測線(4個點，約12 m長)位于島嶼南部水坑上方公路旁邊，起始點坐標(biāo)為22.096 674°N、114.038 623°E，海拔33.2 m。由圖11(a)直流電法探測剖面圖揭示的視電阻率差異可以推出，該剖面中部的地下5～13 m范圍應(yīng)該是有基巖裂隙水儲存，且在該地區(qū)下方約10 m處出現(xiàn)泉水出露，泉位置為22.097 272°N、114.038 793°E，海拔21.4 m。第二條測線(6個點，約30 m長)位于島嶼南部水坑上方公路旁邊，其起始點坐標(biāo)為22.096 674°N、114.038 623°E，海拔33.2 m，結(jié)果如圖11(b)所示。

圖10 研究區(qū)南部直流電法作業(yè)測線位置

(a) 測線1

根據(jù)電阻率差異，可以推測出該剖面中部，地下5～13 m范圍應(yīng)該是有基巖裂隙水儲存，且在該地區(qū)下方約10 m處出現(xiàn)泉水出露。泉位置為22.095 464°N、114.038 102°E，海拔8.8 m。

將直流電法探測解譯的水位值與地下水流場模擬值(圖9)進(jìn)行對比，如表3所示。

表3 探測水位與模擬水位的對應(yīng)值

4.2 相關(guān)性分析

圖12所示為表3中同序號水位模擬值(橫坐標(biāo))與實際水位值(縱坐標(biāo))散點及其擬合關(guān)系曲線，可見兩者呈線性正相關(guān)，且R2=0.872 2(P<0.05)，相關(guān)性很好，說明模擬水位與實際水位結(jié)果符合很好。因此，將基于地下水遙感評估所得出的初始水位導(dǎo)入GMS軟件并加設(shè)相應(yīng)的島嶼水文地質(zhì)參數(shù)所建立的流場模擬模型能夠反映該島嶼的地下水流場分布情況。

圖12 模擬水位與實際探測水位值的相關(guān)性分析

5 結(jié)論

本文針對基巖島嶼地下水流場參數(shù)難以全面獲取的現(xiàn)狀，采用遙感、物探和水文地質(zhì)綜合方法與技術(shù)對外伶仃島地下水流場進(jìn)行了模擬研究，得到了以下認(rèn)識和結(jié)論。

(1)利用數(shù)字高程模型數(shù)據(jù)(普遍是30 m精度的數(shù)據(jù)，也可利用15 m或10 m等更高精度的數(shù)據(jù))，并轉(zhuǎn)化好格式導(dǎo)入地下水模擬軟件或程序中，這樣提高了地質(zhì)體的模型概化精度。

(2)利用物探技術(shù)(探地雷達(dá)法與直流電法相互佐證)對島嶼地層進(jìn)行了探測，獲取了各地層厚度及分界面等數(shù)據(jù)，并結(jié)合水文地質(zhì)試驗參數(shù)(抽水試驗參數(shù)、不同巖性的經(jīng)驗參數(shù))來合理地設(shè)定地下水模擬的各地層參數(shù)。

(3)基于地下水遙感的評估得分與該地區(qū)的各點的水位(井與泉水域范圍的探測值)進(jìn)行擬合，再利用擬合關(guān)系式推出各個點的初始水位，并將所得的初始水位導(dǎo)入地下水模擬軟件或程序中作為地下水的初始水位，且設(shè)定相應(yīng)的初始條件與邊界條件進(jìn)行流場模擬。該方法克服了基巖島嶼的初始水位難以全面獲取的困難，且模擬出的水位模擬值與探測水位值呈正相關(guān)，其兩者的擬合精度也滿足要求。

(4)遙感、物探和水文地質(zhì)綜合方法與技術(shù)可以較好地對基巖島嶼地下水流場進(jìn)行模擬。但是本研究仍然是利用現(xiàn)有的地下水模擬軟件基于等效理論來模擬基巖島的裂隙水(除松散層孔隙水等以外)，將來希望可以嘗試編程實現(xiàn)基于立方定律等理論的裂隙水的直接模擬。