不透水面影響下海島淡水透鏡體的形成與演化

凌子涵，束龍倉，王 然

(1.河海大學(xué)水文水資源學(xué)院,江蘇 南京 210098； 2.河海大學(xué)水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 江蘇 南京 210098)

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

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

1 研究區(qū)概況

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

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

2 數(shù)值模型

2.1 水文地質(zhì)概念模型

本文以研究區(qū)的水文地質(zhì)條件為依據(jù)，建立二維模型觀察淡水透鏡體的形態(tài)變化。選取島嶼截面長度1 980 m，模型的最大深度取到海平面以下48 m，海平面以上島嶼高程設(shè)定為5 m。將模擬區(qū)域地層整體分為3層：第一層為地表至地表以下6 m，該層主要是填土及珊瑚碎屑；第二層為地表以下6～22 m，為全新世地層；地表以下22～54 m為更新世地層。參照該海島的地理特性及文獻(xiàn)[20]，3層地層結(jié)構(gòu)滲透系數(shù)k的初始值分別取60 m/d、150 m/d、1 000 m/d。

海島淡水透鏡體的補(bǔ)給來源為降雨補(bǔ)給。Ghassemi等[21]根據(jù)統(tǒng)計(jì)資料對Home島淡水透鏡體進(jìn)行三維數(shù)值模擬時(shí)，降雨入滲補(bǔ)給系數(shù)取0.44；周從直等[22]對永興島淡水透鏡體進(jìn)行三維數(shù)值建模時(shí)，取多年平均降雨的40%作為補(bǔ)給量。本文結(jié)合之前學(xué)者的經(jīng)驗(yàn)，考慮該海島的水文氣象條件，降雨入滲補(bǔ)給系數(shù)取0.4。模型其他參數(shù)設(shè)置如下：淡水密度為1 g/cm3，海水密度為1.025 g/cm3，海水質(zhì)量濃度為0.019 g/cm3，彈性貯水率為10-5m-1。模型頂部接受大氣降雨補(bǔ)給且存在蒸發(fā)，設(shè)定為流量邊界；補(bǔ)給濃度取雨水中氯離子質(zhì)量濃度，近似為0 g/cm3；模型底部含水層與海水連通，不存在水量交換，因此可視為零流量邊界[22]。模型側(cè)向邊界在海平面以上部分質(zhì)量濃度為0 g/cm3,在海平面以下的部分被海水包圍，故取海水的氯離子質(zhì)量濃度，即0.019 g/cm3。模型初始水頭取0 m，概化后模型示意圖如圖1所示。

圖1 水文地質(zhì)模型

2.2 數(shù)學(xué)模型

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

(1)

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

多孔介質(zhì)的溶質(zhì)運(yùn)移方程包括水動(dòng)力彌散項(xiàng)、地下水對流項(xiàng)、源匯項(xiàng)以及反應(yīng)項(xiàng)。溶質(zhì)運(yùn)移方程為

(2)

式中：θ為有效孔隙度；C為溶質(zhì)質(zhì)量濃度；qs為單位時(shí)間內(nèi)進(jìn)入單位體積含水層的源匯項(xiàng)體積；ρs為源匯項(xiàng)溶質(zhì)質(zhì)量濃度；Dij為水動(dòng)力彌散系數(shù)(i,j= 1, 2, 3)；vi為流體速度；∑Rn為化學(xué)物質(zhì)反應(yīng)項(xiàng)。

2.3 模型離散化

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

3 模型模擬結(jié)果與分析

3.1 淡水透鏡體模擬結(jié)果

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

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

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

圖3 島嶼地下水水流流速

該島嶼淡水透鏡體最大厚度隨時(shí)間的變化如圖4所示，可以看出透鏡體形成初期(0～10 a)厚度增加速度逐漸加快，到中期(10～40 a)速度逐漸下降，40 a后透鏡體厚度變化十分微小直至不變化。如圖5所示，透鏡體潛水面呈向上凸起的形狀，中間厚、兩側(cè)薄，中間水頭最大超過0.33 m，邊緣水頭最小為0 m，水頭向四周遞減形成水力坡降，可以保證淡水從透鏡體最厚的地方不斷排泄入海。

圖4 淡水透鏡體最大厚度隨時(shí)間變化

圖5 淡水透鏡體潛水面形狀

3.2 不透水面對透鏡體的影響

島嶼城市化過程中土地利用率不斷加大，海水入侵現(xiàn)象更加明顯，島嶼淡水儲(chǔ)量會(huì)隨之減少。為探究不透水面對淡水透鏡體的影響，在模型0～594 m部分將降雨入滲補(bǔ)給系數(shù)設(shè)置為0.01，將該部分界面作為不透水面[25]，蒸發(fā)量設(shè)為0，將原始島嶼模擬的水頭及質(zhì)量濃度結(jié)果作為初始條件輸入模型中，模型大約經(jīng)過15 a后再次達(dá)到穩(wěn)定，模擬結(jié)果如圖6所示，可以看出設(shè)置了不透水面一側(cè)的島嶼，由于降雨入滲補(bǔ)給的減少，淡水透鏡體向島嶼內(nèi)部凹陷，淡水透鏡體的最大厚度減小至14.6 m，且淡水透鏡體厚度整體上變薄，經(jīng)計(jì)算穩(wěn)定時(shí)島嶼地下淡水儲(chǔ)量為3 408 m3/m，相比初始時(shí)減少了20.23%。

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

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

圖7 海島內(nèi)部地下水流向

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

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

3.3 模型靈敏度分析

本文采用擾動(dòng)分析法中的Morris法作為靈敏度計(jì)算方法，表示參數(shù)變化對模擬結(jié)果影響程度的靈敏度S的計(jì)算公式[26]為

(3)

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

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

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

圖9 參數(shù)的靈敏度分析

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

4 結(jié) 論

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

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

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