許顥礫, 王大慶, 鄧正棟, 丁志斌, 趙小蘭, 劉志新, 許新剛, 蘇合巖
(1.陸軍工程大學(xué)國(guó)防工程學(xué)院,南京 210007; 2.中國(guó)礦業(yè)大學(xué)資源與地球科學(xué)學(xué)院,徐州 221116)
島嶼的淡水資源十分寶貴,了解其地下水資源狀況非常必要。美國(guó)較早開(kāi)發(fā)利用了沙質(zhì)島嶼的“淡水透鏡體”[1-2],并利用電法勘探確定其厚度和范圍[3]。在20世紀(jì)90年代,我國(guó)海軍和后勤學(xué)院研究解決西沙的供水問(wèn)題時(shí),對(duì)海島“淡水透鏡體”的有關(guān)理論開(kāi)展過(guò)研究,并獲得了豐碩的成果[4-5]。中國(guó)地質(zhì)大學(xué)陳崇希教授和李國(guó)敏研究員運(yùn)用水頭、濃度相互依賴(lài)的有限元法對(duì)潿洲島的火山碎屑孔隙含水層及玄武巖孔洞-裂隙含水層的水位、水流、水質(zhì)等參數(shù)進(jìn)行了公式計(jì)算和模擬分析[6]。中國(guó)科學(xué)院龐忠和研究員聯(lián)合南京大學(xué)對(duì)我國(guó)東部沿海的廟島群島進(jìn)行水文地質(zhì)、物探等方面的研究,驗(yàn)證了基巖島嶼“淡水蘑菇體”的相關(guān)理論[7-8]。隨著計(jì)算機(jī)模擬技術(shù)的發(fā)展,MODFLOW程序、Visual Modflow軟件[9]、GMS軟件[10]、FEFLOW軟件[11]和TOUGH軟件[12]等模擬系統(tǒng)相繼出現(xiàn),地下水模擬變得簡(jiǎn)單化、流程化。在模擬過(guò)程中,需要對(duì)基巖島嶼的地下水補(bǔ)徑排及海水入侵情況進(jìn)行詳細(xì)調(diào)查,通過(guò)試驗(yàn)獲得研究區(qū)地質(zhì)參數(shù),并利用軟件建立模型進(jìn)行地下水的模擬[13]。針對(duì)基巖島嶼的模擬,周鵬鵬等[14-15]根據(jù)有限的水井?dāng)?shù)據(jù)和降雨蒸發(fā)數(shù)據(jù),采用Visual Modflow軟件對(duì)湛江東海島及其周邊地區(qū)(面積約286 km2)進(jìn)行了地下水模擬。中國(guó)地質(zhì)大學(xué)溫漢輝等[16]對(duì)雷州半島(約13 225 km2)的水文、地質(zhì)(主要為巖漿巖和變質(zhì)巖)、地貌、氣象等信息進(jìn)行了較為詳細(xì)的調(diào)查,利用GMS軟件對(duì)該半島進(jìn)行了地下水?dāng)?shù)值模擬,并對(duì)地下水的合理開(kāi)發(fā)提出了寶貴建議。上述研究是對(duì)較大尺度(幾百到上萬(wàn)平方公里)范圍的基巖島嶼地下水流場(chǎng)進(jìn)行模擬,基于達(dá)西定律將基巖島嶼概化為孔隙介質(zhì)從而得到地下水流場(chǎng)的大體分布。顯然,這些模擬工作與實(shí)際情況有較大的出入,基巖島嶼的地下水類(lèi)型多為裂隙水,僅表層的風(fēng)化層為孔隙介質(zhì)(除深層非基巖地下水),且基巖島嶼的裂隙水分布不均勻。針對(duì)基巖島嶼地下水模擬的諸多難點(diǎn),本文以珠海外伶仃島為例,利用遙感技術(shù)、物探技術(shù)、水文地質(zhì)方法及地下水模擬技術(shù)來(lái)研究基巖島嶼的地下水流場(chǎng),改善了因地質(zhì)條件復(fù)雜、水井?dāng)?shù)據(jù)等初始條件較少而難以全面對(duì)基巖島嶼地下水狀況進(jìn)行模擬的現(xiàn)狀,對(duì)同類(lèi)型島嶼的地下水流場(chǎng)模擬工作有參考意義。
以基巖島嶼珠海外伶仃島為對(duì)象展開(kāi)研究。外伶仃島位于廣東省珠海市萬(wàn)山海洋開(kāi)發(fā)試驗(yàn)區(qū)的東北部, 西距萬(wàn)山區(qū)桂山島17.6 n mile(1n mile=1 852 m),離珠海市香洲區(qū)27.5 n mile,北距香港長(zhǎng)洲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)出,主要由黑云母花崗巖組成,屬查氏分類(lèi)第II類(lèi)3科硅酸過(guò)飽和堿性巖石。同位素年齡為160~138 Ma,主要集中在155~140 Ma,歸屬于晚侏羅世。島嶼面積約為4. 33 km2,島嶼周?chē)暮0毒€長(zhǎng)12.3 km。整個(gè)島嶼長(zhǎng)約3 200 m,寬約2 400 m,島嶼地勢(shì)呈西北中部、北部高,西北近海岸、南近海岸低,主峰外伶仃峰海拔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月份有臺(tái)風(fēng),最大風(fēng)速12級(jí)以上,10—12 月份天氣平和。四季溫差不大,冬天無(wú)嚴(yán)寒,夏天無(wú)酷暑。大萬(wàn)山氣象站1973—1991年的氣溫資料表明,有85%的年份冬季極端低溫不低于4 ℃,大多數(shù)年份夏季極端高溫不超過(guò)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; 海水溫度年變化與太陽(yáng)輻射量有關(guān),每年3月開(kāi)始升溫,7—8月水溫最高,9月開(kāi)始降溫,1—2月水溫最低,海水溫度年平均為22.7 ℃; 海水鹽度1月份為33.19‰,7月份為23.59‰,平均為28.39‰。全年之中,冬季的波浪高于夏季,一般波浪不高,平均為0.9~1.9 m,受臺(tái)風(fēng)直接影響時(shí),平均波高5~6 m,最大波高8 m。受強(qiáng)臺(tái)風(fēng)直接影響時(shí),平均波高7~8 m,最大波高大于10 m。外伶仃海域的潮汐屬不正規(guī)半日潮,最高潮位為3.4 m,最低潮位為0.84 m,最大潮差為2.52 m。
本文重點(diǎn)關(guān)注島嶼地表以下第1層含水層(基巖島嶼地下的風(fēng)化裂隙水或構(gòu)造裂隙水含水層),多數(shù)情況為潛水層,有時(shí)也指上層滯水或承壓水。利用地下水遙感評(píng)估法[18-19]所獲得的第一層含水層的地下水遙感評(píng)估分值(S)與水井實(shí)測(cè)水量呈正相關(guān)線性關(guān)系,再考慮水的單位面積,并將井或泉位置的水量取平均(目的是將水量統(tǒng)一以單位投影面積來(lái)計(jì)算水位埋深),計(jì)算公式為
(1)
式中:h為單位水位埋深,m; Δs為單位面積,m2;S為第一層含水層的地下水遙感評(píng)估值,m3; a、b為常數(shù)。
將S與h進(jìn)行擬合,公式為
S=a1·eb1hΔs
,
(2)
式中:h為單位水位埋深,m; Δs為單位面積,m2;S為第一層含水層的地下水遙感評(píng)估值,m3;a1、b1為常系數(shù)。
求出系數(shù)a1、b1,可計(jì)算整個(gè)研究區(qū)第一層含水層的地下水水位埋深(H),其公式為
H=D-h
,
(3)
式中:D為數(shù)字高程模型數(shù)據(jù)的值,m;H為第1層含水層地下水水位值,m;h為地下水埋深值,m。
將第1層含水層地下水水位H矩陣導(dǎo)入地下水模擬軟件或程序進(jìn)行地下水流場(chǎng)模擬的初始水位值,可以克服基巖島嶼難以全面獲得第一層含水層的地下水的初始水位標(biāo)高的困難。
珠江口外海島的地質(zhì)地貌結(jié)構(gòu)簡(jiǎn)單,多為花崗巖。海島巖層自上而下一般可分為3層,分別為地表風(fēng)化層土壤(植被覆蓋)和裸露巖石、半風(fēng)化層以及花崗巖基巖(圖2)。圖2(a)是位于島嶼北側(cè)的一處天然水池,常年有水,為泉水出露和降雨補(bǔ)給形成,整個(gè)島有多處類(lèi)似的天然水池(或稱(chēng)水坑),圖2(b)可以較好地反映地表1~2 m的松散風(fēng)化層,圖2(c)為半風(fēng)化層,圖2(d)為裸露的花崗巖石蛋,圖2(e)為花崗巖基巖底層。根據(jù)外伶仃島的有降雨補(bǔ)給和地下水補(bǔ)給形成的地表水域及以前人們使用的老井及泉水出露,已經(jīng)可以說(shuō)明島嶼存在地下淡水,但是容量有限。島嶼周邊都是海水,地下水的唯一補(bǔ)給源就是雨水。對(duì)于海島的供水系統(tǒng)而言,以前人口少時(shí),主要利用水井地下水或積累雨水的水池。隨著旅游業(yè)發(fā)展及當(dāng)?shù)卣块T(mén)的開(kāi)發(fā),海島人口逐步增多,原有的供水體系已無(wú)法滿(mǎn)足需要,于是就用船運(yùn)來(lái)的水通過(guò)高壓泵和管道輸送到各個(gè)用水地點(diǎn)(如外伶仃島),或者是在當(dāng)?shù)亟⑿⌒退畮?kù)來(lái)滿(mǎn)足島民的生活用水(桂山島和東澳島)。
(a) 天然水池 (b) 地表松散風(fēng)化層
對(duì)島嶼花崗巖取樣測(cè)試,并根據(jù)《水文地質(zhì)手冊(cè)》[20]、《堤防工程手冊(cè)》[21]和《地下水污染物遷移模擬》[9]所給出的水文地質(zhì)參數(shù)檢驗(yàn)值,將水平滲透系數(shù)、豎直滲透系數(shù)、給水度、儲(chǔ)水系數(shù)和各向異性參數(shù)等水文地質(zhì)參數(shù)輸入地下水模擬軟件或程序(表1)進(jìn)行模擬計(jì)算。
該地區(qū)的平均降雨量約為0.005 33 m/d(數(shù)據(jù)來(lái)自珠海氣象局),由于該基巖島嶼的特殊性,在基巖海岸邊緣是以淡水“蘑菇體”的形式存在,根據(jù)地下水模擬技術(shù)國(guó)內(nèi)外研究現(xiàn)狀[13],可以將基巖海岸線設(shè)成定為模擬的范圍邊界。
表1 各種巖層的水文地質(zhì)參數(shù)
本次采用美國(guó)GSSI探地雷達(dá)SIR-4000儀器探測(cè)島嶼地層分層情況。本文利用100 MHz天線,可以對(duì)地下8~20 m的范圍進(jìn)行探測(cè)。為佐證探地雷達(dá)的探測(cè)結(jié)果,采用法國(guó)SYSCALR2直流電法儀進(jìn)行了三極裝置直流電法探測(cè)。如圖3所示,沿島嶼西南部的一條公路進(jìn)行了約1.5 km的探地雷達(dá)作業(yè)與直流電法作業(yè)(圖3)。提取同一段水平距離220 m測(cè)線的探地雷達(dá)探測(cè)結(jié)果(圖4)和直流電法探測(cè)結(jié)果(圖5)進(jìn)行對(duì)比。圖4、圖5按沿測(cè)線的地勢(shì)起伏將探測(cè)結(jié)果繪制在地表之下,以反映地下探測(cè)深度范圍內(nèi)的地球物理分層。如圖4,根據(jù)地質(zhì)體介電常數(shù)的差異可以較清晰地看出: 厚度1~2 m的地表層與圖2所示的地表巖體風(fēng)化層對(duì)應(yīng); 厚度5~7 m是另外一個(gè)地質(zhì)層,可以對(duì)應(yīng)于半風(fēng)化層; 厚度5~7 m以下的部分可以對(duì)應(yīng)于花崗巖基巖底層。圖5可反映如下情況: 低視電阻率的藍(lán)色區(qū)域厚度為1~2 m; 中等大小視電阻率(綠色及黃色)厚度為2~5 m; 5~7 m深度之下的高視電阻率(橙色和紅色)對(duì)應(yīng)花崗巖基巖底層??梢?jiàn),探地雷達(dá)與直流電法的探測(cè)結(jié)果符合較好,與現(xiàn)場(chǎng)勘查情況對(duì)應(yīng)也很好,可以確定島嶼在地下水范圍內(nèi)的地質(zhì)分層情況為: 地表之下1~2 m為地表風(fēng)化層; 地表下5~7 m為半風(fēng)化層; 地表7 m以下為花崗巖基巖底層。
圖3 外伶仃島西南部物探測(cè)線位置
圖4 島嶼西南部的探地雷達(dá)探測(cè)220 m水平距離結(jié)果
圖5 島嶼西南部的直流電法探測(cè)220 m水平距離結(jié)果
根據(jù)研究區(qū)地下水富集性評(píng)價(jià)結(jié)果等級(jí)圖,可得出每一點(diǎn)地下水富集性遙感評(píng)估值(S),再根據(jù)平均水量直接用該點(diǎn)實(shí)際單位水位埋深h來(lái)表示,即S=a1·eb1h。S與h進(jìn)行擬合,得出系數(shù)a1、b1,再利用ENVI軟件計(jì)算整個(gè)研究區(qū),或研究區(qū)的關(guān)鍵點(diǎn)的h,將地下水的水位H矩陣作為該地區(qū)第一層含水層的初始水位數(shù)據(jù)。表2是地下水評(píng)估得分與單位水位埋深關(guān)系表。
表2 地下水評(píng)估得分與單位水位埋深關(guān)系[17]
根據(jù)井或泉觀測(cè)值的平均出水量得到單位時(shí)間單位面積的水位埋深,將h與S進(jìn)行曲線擬合(圖6),得出擬合曲線關(guān)系式為
S=0.427e0.239h,
(4)
式中:h為單位水位埋深,m;S為地下水遙感評(píng)估得分,m3; 0.427和0.239為系數(shù)。
此時(shí),擬合優(yōu)度R2=0.847 2(顯著性P<0.05),說(shuō)明井或泉中對(duì)應(yīng)的單位時(shí)間單位面積的水位埋深h與地下水富集性遙感評(píng)估得分S的擬合程度較好。
圖6 水位埋深h與評(píng)估得分S的擬合曲線
根據(jù)公式(4)反推出各個(gè)點(diǎn)的評(píng)估值S對(duì)應(yīng)h的計(jì)算關(guān)系式為
(5)
式中:h為單位水位埋深,m;S為地下水遙感評(píng)估值,m3; 0.427和0.239為系數(shù)。
根據(jù)式(3)和式(5),利用ENVI軟件將地下水遙感評(píng)估得分S值計(jì)算整個(gè)島的初始水位H值,即基巖島嶼的地下水流場(chǎng)模擬的初始水位用H值表示。
建立研究區(qū)地下水模型分為4個(gè)步驟。
(1)根據(jù)研究區(qū)的水位地質(zhì)資料背景和實(shí)際勘測(cè)島嶼資料,設(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ù)值,則補(bǔ)給為二者之間的差值,即降雨蒸發(fā)量為+0.005 33 m/d,在有水坑、小型湖泊的位置設(shè)置補(bǔ)給為0.006 03 m/d,并設(shè)置相應(yīng)井的抽水量。
(4)將DEM數(shù)據(jù)與H值的TIF導(dǎo)入GMS,并轉(zhuǎn)化為散點(diǎn)(圖7)。DEM數(shù)據(jù)的利用,相對(duì)于傳統(tǒng)的人工跑點(diǎn)用GPS點(diǎn)位而言,可以提高研究區(qū)地質(zhì)體的概化精度。
圖7 DEM的TIF文件和H的TIF文件在GMS上轉(zhuǎn)化為散點(diǎn)
在GMS軟件上進(jìn)行網(wǎng)格剖分,本模擬將該島嶼剖分為3層,分別是地表松散風(fēng)化層、第1層(含水層)和第2層(隔水層),本文主要針對(duì)地下水的流場(chǎng)分布,再利用GMS自帶的Check Simulation功能進(jìn)行自動(dòng)查錯(cuò),修正后進(jìn)行MODFLOW模塊計(jì)算。
地表松散風(fēng)化層(1~2 m)水位及流場(chǎng)模擬結(jié)果如圖8,第1層(含水層)(2~5 m)地下水水位及流場(chǎng)模擬結(jié)果如圖9。
圖8與圖9經(jīng)比較可以看出,地表松散風(fēng)化層中的水位與第1層含水層的地下水水位有明顯的差異,雖然二者的地下水流場(chǎng)的方向趨勢(shì)大致相同,但表面松散風(fēng)化層的地下水分布很不均勻,且較大范圍的風(fēng)化層處于無(wú)水狀態(tài),而地下水遍布全島,第1層含水層的地下水流場(chǎng)的水位分布較地表風(fēng)化層的地下水流場(chǎng)更明顯。從圖9反映出的第一層含水層地下水水位模擬結(jié)果來(lái)看,水位角度呈中部和東南部高、周邊低的趨勢(shì)。從島嶼第一層含水層的地下水流場(chǎng)趨勢(shì)上看,從中部向外流,從東南部向外流,周邊又有向島內(nèi)流的趨勢(shì),應(yīng)該是有海水的存在的原因。
圖8 地表松散風(fēng)化層(1~2 m)水流場(chǎng)模擬結(jié)果
圖9 第一層(含水層)(2~5 m)地下水流場(chǎng)模擬結(jié)果
本研究利用直流電法儀以對(duì)稱(chēng)四極裝置形式,在供電極距AB=30 m條件下進(jìn)行實(shí)地探測(cè)。直流電法的2條測(cè)線具體位置如圖10中的紅色測(cè)線位置,第一條測(cè)線(4個(gè)點(diǎn),約12 m長(zhǎng))位于島嶼南部水坑上方公路旁邊,起始點(diǎn)坐標(biāo)為22.096 674°N、114.038 623°E,海拔33.2 m。由圖11(a)直流電法探測(cè)剖面圖揭示的視電阻率差異可以推出,該剖面中部的地下5~13 m范圍應(yīng)該是有基巖裂隙水儲(chǔ)存,且在該地區(qū)下方約10 m處出現(xiàn)泉水出露,泉位置為22.097 272°N、114.038 793°E,海拔21.4 m。第二條測(cè)線(6個(gè)點(diǎn),約30 m長(zhǎng))位于島嶼南部水坑上方公路旁邊,其起始點(diǎn)坐標(biāo)為22.096 674°N、114.038 623°E,海拔33.2 m,結(jié)果如圖11(b)所示。
圖10 研究區(qū)南部直流電法作業(yè)測(cè)線位置
(a) 測(cè)線1
根據(jù)電阻率差異,可以推測(cè)出該剖面中部,地下5~13 m范圍應(yīng)該是有基巖裂隙水儲(chǔ)存,且在該地區(qū)下方約10 m處出現(xiàn)泉水出露。泉位置為22.095 464°N、114.038 102°E,海拔8.8 m。
將直流電法探測(cè)解譯的水位值與地下水流場(chǎng)模擬值(圖9)進(jìn)行對(duì)比,如表3所示。
表3 探測(cè)水位與模擬水位的對(duì)應(yīng)值
圖12所示為表3中同序號(hào)水位模擬值(橫坐標(biāo))與實(shí)際水位值(縱坐標(biāo))散點(diǎn)及其擬合關(guān)系曲線,可見(jiàn)兩者呈線性正相關(guān),且R2=0.872 2(P<0.05),相關(guān)性很好,說(shuō)明模擬水位與實(shí)際水位結(jié)果符合很好。因此,將基于地下水遙感評(píng)估所得出的初始水位導(dǎo)入GMS軟件并加設(shè)相應(yīng)的島嶼水文地質(zhì)參數(shù)所建立的流場(chǎng)模擬模型能夠反映該島嶼的地下水流場(chǎng)分布情況。
圖12 模擬水位與實(shí)際探測(cè)水位值的相關(guān)性分析
本文針對(duì)基巖島嶼地下水流場(chǎng)參數(shù)難以全面獲取的現(xiàn)狀,采用遙感、物探和水文地質(zhì)綜合方法與技術(shù)對(duì)外伶仃島地下水流場(chǎng)進(jìn)行了模擬研究,得到了以下認(rèn)識(shí)和結(jié)論。
(1)利用數(shù)字高程模型數(shù)據(jù)(普遍是30 m精度的數(shù)據(jù),也可利用15 m或10 m等更高精度的數(shù)據(jù)),并轉(zhuǎn)化好格式導(dǎo)入地下水模擬軟件或程序中,這樣提高了地質(zhì)體的模型概化精度。
(2)利用物探技術(shù)(探地雷達(dá)法與直流電法相互佐證)對(duì)島嶼地層進(jìn)行了探測(cè),獲取了各地層厚度及分界面等數(shù)據(jù),并結(jié)合水文地質(zhì)試驗(yàn)參數(shù)(抽水試驗(yàn)參數(shù)、不同巖性的經(jīng)驗(yàn)參數(shù))來(lái)合理地設(shè)定地下水模擬的各地層參數(shù)。
(3)基于地下水遙感的評(píng)估得分與該地區(qū)的各點(diǎn)的水位(井與泉水域范圍的探測(cè)值)進(jìn)行擬合,再利用擬合關(guān)系式推出各個(gè)點(diǎn)的初始水位,并將所得的初始水位導(dǎo)入地下水模擬軟件或程序中作為地下水的初始水位,且設(shè)定相應(yīng)的初始條件與邊界條件進(jìn)行流場(chǎng)模擬。該方法克服了基巖島嶼的初始水位難以全面獲取的困難,且模擬出的水位模擬值與探測(cè)水位值呈正相關(guān),其兩者的擬合精度也滿(mǎn)足要求。
(4)遙感、物探和水文地質(zhì)綜合方法與技術(shù)可以較好地對(duì)基巖島嶼地下水流場(chǎng)進(jìn)行模擬。但是本研究仍然是利用現(xiàn)有的地下水模擬軟件基于等效理論來(lái)模擬基巖島的裂隙水(除松散層孔隙水等以外),將來(lái)希望可以嘗試編程實(shí)現(xiàn)基于立方定律等理論的裂隙水的直接模擬。