凌子涵，束龍倉，王 然

(1.河海大學水文水資源學院,江蘇 南京 210098； 2.河海大學水文水資源與水利工程科學國家重點實驗室, 江蘇 南京 210098)

我國擁有豐富的海洋資源，海島作為人類探索海洋的重要支點，對于開發(fā)海洋資源、建設海洋經濟大國具有重要意義。本文所研究的海島是西沙群島中最大的灰沙島嶼，為三沙市政治、經濟、文化中心，島上常住人口2 500多人[1]。近年來隨著三沙市對該海島的大力建設，島上現已具備良好的居住生活條件，但面臨著用水形勢嚴峻的問題，目前主要通過補給船定期補給以解決島上居民的飲用水問題[2]，但這種方式運水，水量十分有限，且成本昂貴。除此之外還有海水淡化、雨水收集等方法，但均無法滿足日趨增長的島上用水需求。

海島獨特的地質結構使得島嶼地下淡水被咸水包圍，并浮于咸水之上，被稱作“淡水透鏡體”。近年來國內外學者對淡水透鏡體的形成、演變及影響因素等進行了深入研究[3-7]。李國敏等[8]對潿洲島含水層條件進行概化，并建立數值模型，模擬了海水入侵的3D過程。Vacher[9]利用Ghyben-Herzberg公式對Bermuda島和Great Exuma島淡水透鏡體的形態(tài)進行了模擬計算。James等[10]使用SEAWAT軟件，模擬海岸侵蝕對美國佛羅里達州海島淡水透鏡體形態(tài)與厚度的影響。甄黎等[11]通過物理模型模擬了不同開采條件下倒錐的變化過程。本文在前人研究的基礎上，利用Visual MODFLOW軟件構建了所研究海島淡水透鏡體的二維模型，并在透鏡體穩(wěn)定后設置不透水面，通過改變不透水面的覆蓋范圍，模擬透鏡體形態(tài)及儲水量的變化。此外，對影響透鏡體儲水量的不同因素進行靈敏度分析，為島嶼淡水資源開采提供科學參考。

1 研究區(qū)概況

本文所研究海島為海南省三沙市永興島，朝海部分的島嶼經礁坪前緣以大角度進入深海盆地[12]，全島在北回歸線以南，緯度低，地處熱帶地區(qū)，屬赤道帶、熱帶海洋性季風氣候，島嶼年平均氣溫在26℃以上，年降雨量超過1 400 mm，其中6—11月降水量占全年85%以上。

該海島的地質結構大體上可以分為全新統和更新統，地表以下約22 m處為不整合面，不整合面以上為全新世珊瑚貝殼砂地層，淡水透鏡體主要存在于該地層；不整合面以下由巖溶十分發(fā)育的更新世石灰?guī)r構成，由于原生礁不斷受到海水侵蝕，石灰?guī)r體被孔隙通道貫通，巖溶發(fā)育，因此滲透性很強[13]。鉆孔資料及其他文獻[14-19]表明：島嶼地表以下6 m由素填土和珊瑚碎屑中砂組成，地表以下6～22 m由珊瑚碎屑礫和珊瑚碎屑礫砂組成，22 m以下為珊瑚礁灰?guī)r；島嶼介質孔隙度為0.4，給水度為0.2；含水介質縱向彌散度為5 m，橫縱比為0.1。此外，該海島淡水補給的主要來源是降雨入滲補給，排泄方式主要是垂向的蒸發(fā)排泄和側向的向海洋排泄。

2 數值模型

2.1 水文地質概念模型

本文以研究區(qū)的水文地質條件為依據，建立二維模型觀察淡水透鏡體的形態(tài)變化。選取島嶼截面長度1 980 m，模型的最大深度取到海平面以下48 m，海平面以上島嶼高程設定為5 m。將模擬區(qū)域地層整體分為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。模型頂部接受大氣降雨補給且存在蒸發(fā)，設定為流量邊界；補給濃度取雨水中氯離子質量濃度，近似為0 g/cm3；模型底部含水層與海水連通，不存在水量交換，因此可視為零流量邊界[22]。模型側向邊界在海平面以上部分質量濃度為0 g/cm3,在海平面以下的部分被海水包圍，故取海水的氯離子質量濃度，即0.019 g/cm3。模型初始水頭取0 m，概化后模型示意圖如圖1所示。

圖1 水文地質模型

2.2 數學模型

由于海水與淡水可以混溶，因此所構建海島淡水透鏡體二維剖面模型存在可混溶液體間的水動力彌散問題[23]。數學模型由控制方程、初始條件以及邊界條件構成，依據上述水文地質概念模型，建立研究區(qū)地下水流與溶質運移數學模型，通過以下微分方程的定解問題來描述：

(1)

式中：K1為潛水含水層滲透系數；μ為潛水含水層給水度；H為潛水水位；B為潛水含水層底板標高；W1為源匯項；H0為初始水位；H1為研究區(qū)一類邊界點的水位；x,z為坐標；D為計算區(qū)范圍；Γ1為研究區(qū)一類邊界：t為時間。

多孔介質的溶質運移方程包括水動力彌散項、地下水對流項、源匯項以及反應項。溶質運移方程為

(2)

式中：θ為有效孔隙度；C為溶質質量濃度；qs為單位時間內進入單位體積含水層的源匯項體積；ρs為源匯項溶質質量濃度；Dij為水動力彌散系數(i,j= 1, 2, 3)；vi為流體速度；∑Rn為化學物質反應項。

2.3 模型離散化

將研究區(qū)域劃分為若干個網格單元(共30層34列)，每個網格單元內部的水文參數為常數，單元之間存在水力聯系，將分界面以上部分及側向邊界網格進行加密處理。模型模擬期為73 050 d，水頭計算初始時間步長為0.001 d。溶質運移計算初始時間步長1 d，增長因子為2，最大時間步長為10 d，模型結果以月為單位輸出。

3 模型模擬結果與分析

3.1 淡水透鏡體模擬結果

在Visual MODFLOW軟件中，根據實際逐月降雨資料與逐月蒸發(fā)資料設定模型源匯項，模擬期為100 a，選用SEAWAT程序包進行模擬，其中水流計算采用預處理共軛梯度法(PCG)來求解；對流計算采用中心差分法[24]。運行結果以氯離子600 mg/cm3質量濃度等值線[21]為淡水透鏡體下界面輪廓即咸淡水界面。

模型運行結果如圖2所示，隨著降雨入滲的不斷補給，淡水透鏡體厚度變大，在50 a后形狀達到穩(wěn)定，穩(wěn)定時透鏡體最大厚度約為15.6 m。圖3為島嶼地下水流的流速示意圖，箭頭的大小表示流速的快慢，可以看出，降雨入滲的水流由上至下流動，并向兩側排泄；透鏡體內部中心水流流速很低，由于咸淡水之間存在密度差，使得淡水透鏡體中心部分淡水受浮力較大，因此透鏡體下界面輪廓中心部分向上凸起。

圖2 淡水透鏡體形態(tài)演變過程

圖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]，蒸發(fā)量設為0，將原始島嶼模擬的水頭及質量濃度結果作為初始條件輸入模型中，模型大約經過15 a后再次達到穩(wěn)定，模擬結果如圖6所示，可以看出設置了不透水面一側的島嶼，由于降雨入滲補給的減少，淡水透鏡體向島嶼內部凹陷，淡水透鏡體的最大厚度減小至14.6 m，且淡水透鏡體厚度整體上變薄，經計算穩(wěn)定時島嶼地下淡水儲量為3 408 m3/m，相比初始時減少了20.23%。

圖6 設置不透水面后透鏡體形態(tài)

圖7為設置不透水面前后海島內部水流的流向，由圖7可知，未設置不透水面的島嶼降雨入滲補給的水流到達透鏡體形狀邊緣時，部分水流繼續(xù)向下流動，部分水流水平流動，還有一些水流有向上的速度分量；在設置不透水面后，不透水面的存在阻礙了局部的降雨入滲補給，在不透水面存在的一側，水流多為橫向流動，向下流動路徑縮短，并在島嶼側邊向海洋排泄，因此透鏡體輪廓向內凹陷，淡水儲量也相應減少。在島嶼開發(fā)過程中，隨著建設力度的增大，島嶼土地覆蓋中不透水面占比會隨之增加，因此在模型中設置了多種不透水面覆蓋范圍，以分析淡水透鏡體的形態(tài)及儲水量變化規(guī)律。結果表明，隨著不透水面覆蓋范圍的增大，存在不透水面一側的淡水透鏡體向內萎縮趨勢愈加明顯，體積不斷減小，淡水透鏡體最大厚度位置也在逐漸右移。

圖7 海島內部地下水流向

圖8展示了不同不透水面覆蓋范圍占比下島嶼單寬地下淡水儲量的變化，隨著不透水面覆蓋范圍占比的不斷加大，地下淡水儲量逐漸減少，且減少速度逐漸加快，當不透水面覆蓋范圍占比超過40%后，二者逐漸轉變?yōu)榫€性關系。同時可以看出，當不透水面覆蓋范圍占比超過島嶼總長度的50%時，地下淡水儲量減少至原來的50%左右；當不透水面覆蓋范圍占比為60%時，地下淡水儲量已經不足原儲量的40%。

圖8 島嶼地下淡水儲量隨不透水面覆蓋范圍占比變化

3.3 模型靈敏度分析

本文采用擾動分析法中的Morris法作為靈敏度計算方法，表示參數變化對模擬結果影響程度的靈敏度S的計算公式[26]為

(3)

式中：yi為參數變化后模型輸出值；y0為模型原始輸出值；P0為參數初始值；Pi為參數變化后的值；Δi為Pi相對于P0的變化幅度。

影響海島淡水透鏡體形態(tài)及體積的因素有很多，可以分為自然因素和人為因素，自然因素包括降雨入滲補給、蒸發(fā)和島嶼介質滲透性等，人為因素主要是不透水面對透鏡體的影響。為了研究該海島地下淡水儲量在不同影響因素變化下的敏感程度，本文對4項主要影響因素進行靈敏度分析，采用局部靈敏度分析方法，研究單一影響因素變化時對淡水透鏡體體積(文中指地下淡水儲量)的影響。依次選用不透水面覆蓋范圍、含水層給水度、降雨入滲補給系數和滲透系數進行靈敏度計算。

以上4種參數初始值設置如下：不透水面覆蓋范圍594 m；含水層給水度0.2；降雨入滲補給系數0.4；3層地層結構滲透系數分別為60 m/d、150 m/d和1 000 m/d。為了便于計算，參數變化幅度取±10%、±20%和±30%，靈敏度計算及分析結果見圖9。

圖9 參數的靈敏度分析

從靈敏度的計算結果可以得到：隨著不透水面覆蓋范圍或滲透系數的增大，海島地下淡水儲量隨之減少；當降雨入滲補給系數或給水度增大時，海島地下淡水儲量也隨之增加。島嶼地下淡水儲量對于給水度的變化最為敏感，其次是降雨入滲補給系數和島嶼介質滲透系數，再次是不透水面覆蓋范圍。降雨入滲補給系數和滲透系數的靈敏度在初始值對稱增減時呈對稱分布，且隨變幅的增大靈敏度也隨之增大。給水度和不透水面覆蓋范圍靈敏度隨著其自身的增大而增大。

4 結 論

a.對研究島嶼進行水文地質概化，建立數值模型，模型運行期0～10 a透鏡體厚度增長較快，10～40 a增長速度放緩，經過大約50 a達到穩(wěn)定，模擬結果顯示島嶼地下淡水透鏡體最大厚度約為15.6 m，透鏡體潛水面形狀中間厚、兩邊薄。

b.在島嶼長度0～594 m部分設置不透水面，模型經過15a后重新達到穩(wěn)定，淡水透鏡體整體上變薄，最大厚度減小至14.6 m，設置了不透水面一側的透鏡體下界面向內凹陷，島嶼地下淡水儲量為3 408 m3/m，相比初始島嶼減少了20.23%。

c.靈敏度分析結果表明，島嶼地下淡水儲量對于給水度最為靈敏，其次是降雨入滲補給系數、含水層滲透系數和不透水面覆蓋范圍，其中不透水面覆蓋范圍和滲透系數與島嶼地下淡水儲量呈負相關關系，降雨入滲補給系數和給水度與島嶼地下淡水儲量呈正相關關系。