盛 沖，許鶴華，張文濤

(1.中國科學院南海海洋研究所/中國科學院邊緣海與大洋地質(zhì)重點實驗室，廣東 廣州 510301；2.中國科學院大學，北京 100049)

南海諸島自古為中國領(lǐng)土，據(jù)統(tǒng)計，在地貌學上南海共擁有62座干出礁及49個灰沙島，其地緣政治復雜，戰(zhàn)略地位特殊，擁有豐富的海洋資源[1]。其中的永興島是西沙群島中最大的灰沙島嶼，為三沙市政治、經(jīng)濟、文化中心，島上常駐人口2 500多人。隨著近年來三沙市布局規(guī)劃的逐步開展，如何使島上的環(huán)境更加適宜人們長期居住成為了關(guān)注的重點，其中備受關(guān)注是島上淡水的補給問題。由于特殊的地質(zhì)地形結(jié)構(gòu)，島上的淡水資源奇缺，這嚴重制約了日常生活和生態(tài)環(huán)境的改善。根據(jù)1996年中國科學院南沙綜合科學考察隊對南海島礁的考察發(fā)現(xiàn)[2]，在南海某些島嶼上，部分降雨入滲補給含水層，并在補水和損失過程中維持著平衡狀態(tài)，這部分浮在底層海水之上的淡化水體，即“淡水透鏡體”(Freshwater Lens)。

在地貌未發(fā)生大規(guī)模改變前，永興島淡水透鏡體的最大厚度約13.5 m，淡水資源儲量約為147.2×104m3[3～4]。但隨著島上各項工程和規(guī)劃的逐步開展，截止到2018年初，永興島的沙丘面積增加了近50%，新增區(qū)域主要集中在島礁的北部，這一片陸域把石島完全并入了永興島。為了探究新增部分沙丘面積對島礁淡水透鏡體的影響以及之后永興島淡水資源儲量的變化，利用數(shù)值模擬進行預測顯得十分必要。

本文在結(jié)合永興島地貌模型和前人研究的基礎(chǔ)上，利用GMS軟件構(gòu)建了永興島淡水透鏡體的三維變密度耦合模型，并對與島礁淡水透鏡體形成有關(guān)的部分參數(shù)進行了敏感性分析，探討了地貌變化對永興島淡水透鏡體的影響，此外還預測了永興島淡水透鏡體再次達到穩(wěn)定狀態(tài)所需的時間及淡水資源儲量。這對永興島及其他相似島嶼淡水資源的開發(fā)和管理具有重要意義。

圖1 研究區(qū)示意圖及觀測點位置(據(jù)文[5]，有改動)Fig.1 Schematic diagram of the study area and location of the observation points (modified after[5])

1 研究區(qū)概況

永興島屬于中國西沙群島東部的宣德群島，主要由珊瑚和其他生物骨骼堆積而成，四周礁坪環(huán)繞，向海部分經(jīng)礁坪前緣以大角度進入深海盆地[5]。2009年前，永興島外形近似為橢圓平底盆地，東西長約1 950 m，南北寬約1 350 m，平均海拔高度約5 m。但隨著島上各項規(guī)劃的逐步開展，目前永興島的本島面積已由之前的2.13 km2增加到了3.2 km2(包括石島)，島嶼的形狀也發(fā)生巨大改變(圖1)。此外永興島降雨量豐富，多年平均降雨量為1 509 mm；由于島上土體的孔滲性很好且島邊緣有砂堤，降雨很少形成徑流流失，除一部分被植物截留蒸發(fā)外，大部分降至地面入滲補給島下淡水[6]。

目前關(guān)于永興島已有較為詳細的水文地質(zhì)鉆孔資料，地質(zhì)勘探結(jié)果表明[3, 7～8]，原永興島區(qū)域地表以下0～22 m主要為全新統(tǒng)碳酸鹽生物骨殼堆積物，即珊瑚貝殼砂層，由于未發(fā)生成巖作用，呈松散狀態(tài)，孔隙率和滲透率較高，其中藏有淡水透鏡體；22～169 m主要由灰白色珊瑚礁灰?guī)r組成，該地層雖然發(fā)生了成巖作用，但由于更新世時位于海平面以上，所以巖溶體系發(fā)育，地層總的孔滲性高，海水容易滲入(圖2)。根據(jù)前人在永興島上進行的抽水實驗結(jié)果[6]，不整合面以上含水層的滲透系數(shù)為70 m/d，之下的礁灰?guī)r地層由于溶洞體系發(fā)育，滲透系數(shù)取500 m/d[9～10]；給水度為0.14，孔隙度為0.45。彌散度由于存在一定的尺度效應，且難以通過實驗準確測得，在進行數(shù)值模擬計算時一般根據(jù)實驗室測值，并結(jié)合模擬結(jié)果和經(jīng)驗參數(shù)進行確定[11～12]，因此在進行數(shù)值模擬時，彌散度的初始值：縱向彌散度取5 m，橫向彌散度取0.5 m，垂向彌散度取0.05 m。

圖2 研究區(qū)C-C’線水文地質(zhì)剖面(據(jù)文[5，8]，有改動)Fig.2 Hydrogeologic profile of line C-C’ in the study area (modified after [5，8])

永興島的北部區(qū)域主要由近年來海浪及人為因素搬運堆積的珊瑚砂礫等生物碎屑組成，面積約1.2 km2，平均厚度達6 m。為了明確該區(qū)域珊瑚砂的水理特征，對該區(qū)域的幾組珊瑚砂樣進行顆粒分析，具體的顆粒級配曲線見圖3。該區(qū)域珊瑚砂顆粒級配不良且不連續(xù)，此外樣品的有效粒徑d10為40.31 μm、特征粒徑d30為85.18 μm、中值粒徑d50為151.3 μm和控制粒徑d60為208.2 μm。計算所得不均勻系數(shù)Cu為5.615和曲率系數(shù)Cc為0.865，因此定義為細粒土質(zhì)砂，其中含有珊瑚碎片。為了進一步探究該區(qū)域珊瑚砂的其他水文地質(zhì)特征，還在島上進行了雙環(huán)實驗和抽水實驗，相關(guān)的水文地質(zhì)參數(shù)見表1。彌散度需結(jié)合室內(nèi)試驗和經(jīng)驗參數(shù)確定。

圖3 島礁北部珊瑚砂的顆粒級配曲線Fig.3 Grading curve of the coral sand in the northern part of the Yongxing Island

參數(shù)取值滲透系數(shù)K/(m·d-1)5 孔隙度n0.42給水度Sy0.103彌散度α縱向彌散度3 m,橫向彌散度0.3 m,垂向彌散度0.03 m貯水率Ss1×10-5

2 數(shù)值模型的構(gòu)建

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

通過分析研究區(qū)地貌及水文地質(zhì)條件，將永興島模擬區(qū)域的地層概化為三層，其中第一層為地表到海平面以下1 m，主要由北部人工堆積的珊瑚砂礫和原有島嶼未固結(jié)的全新統(tǒng)珊瑚砂組成，四周以海岸帶為邊界，概化為定水頭和定濃度邊界，面積約3.4 km2(圖4)?？紤]到淡化水體向海域延伸的可能，且島礁在礁坪前緣以大角度向海底傾斜，因此模型其他兩層的四周邊界向外延伸至礁坪前緣[13～16]，同樣概化為定水頭和定濃度邊界。其中第二層主要由全新統(tǒng)松散珊瑚砂構(gòu)成，孔滲性良好，為淡水透鏡體的主要賦存層；不整合面到海平面以下45 m處概化為模型的第三層，為更新統(tǒng)固結(jié)的礁灰?guī)r地層，溶洞體系發(fā)育，海水較容易滲入。由于淡水透鏡體懸浮于海水之上，通過孔隙、溶洞與海水相通，潮汐會引起淡水透鏡體整體漲落，但不影響透鏡體內(nèi)部動力學特征，因此不考慮潮汐的影響[6, 8, 17]。

圖4 水文地質(zhì)概念模型示意圖(C-C’剖面)Fig.4 Conceptual hydrogeological model of theYongxing Island (profile along C-C’)

模型的頂部為定流量邊界，由于這部分補給量受降雨強度、植被、地形等多種因素影響，難以確定，因此參考一些經(jīng)驗數(shù)據(jù)及相關(guān)統(tǒng)計資料[6, 18～19]，取年均降雨量的45%均分至每天由地表補給淡水透鏡體。模型底部位于海平面以下45 m處，與海水相連且不存在水量交換，因此可以視為零流量邊界。

關(guān)于模型初始條件的概化，在淡水透鏡體未形成之前，島礁海平面以下區(qū)域完全被海水充填，隨著降雨不斷入滲補給，逐漸在島上形成淡水。因此模型海平面以上的Cl-初始濃度取0 g/L，海平面以下的初始濃度取海水中的Cl-濃度19 g/L；模型初始水頭值取0 m。地貌變化后的島礁淡水透鏡形成初始條件：將此前形成的穩(wěn)定淡水透鏡體的水頭值和濃度值作為原永興島區(qū)域的初始條件；北部新增區(qū)域海平面上、下的Cl-初始濃度分別取0 g/L和19 g/L，初始水頭值取0 m。

2.2 控制方程

在構(gòu)建淡水透鏡體三維數(shù)值模型時，由于海水與淡水屬于可混溶性液體，對于這類問題需要考慮過渡帶的存在以及對流-彌散作用，因此需要2個方程進行描述：水流方程，用來描述密度不斷改變的液體流動；對流-彌散方程用來描述地下水中鹽分的運移[20]。基于質(zhì)量守恒和達西定律確定的變密度水流方程為：

(1)

式中：K——滲透系數(shù)/(m·d-1)；

H——測壓管水頭/m；

c——混合流體密度/(kg·m-3)；

Ss——貯水率；

t——時間/d；

ρ0——淡水密度(參考密度)/(kg·m-3)；

qs——單位體積多空介質(zhì)的源或匯強度/d-1;

η——密度耦合系數(shù)/(kg·m-3)；

ρ——混合流體密度/(kg·m-3)。

在考慮流體密度時，忽略壓力和溫度對流體密度的影響，將流體密度處理成溶質(zhì)濃度的線性函數(shù)，表示為：

(2)

(3)

η=ε/cs

(4)

式中：ρs——海水密度/(kg·m-3)；

cs——海水密度對應的質(zhì)量濃度/(kg·m-3)；

ε——密度差相對比率。

溶質(zhì)運移方程為：

(5)

式中：ui——滲透流速/(m·d-1)；

Dii——彌散系數(shù)/(m2·d-1)；

c*——源或匯的流體濃度/(kg·m-3)；

n——多孔介質(zhì)孔隙度。

2.3 網(wǎng)格剖分及時間離散

模擬區(qū)域平面上的網(wǎng)格大小為32 m×26 m，垂向共分為30層，其中由于淡水透鏡體主要存在于不整合面之上，因此對不整合面所處的位置及以上的網(wǎng)格進行細化，模型有效單元格共計294 171個。

模型總模擬時間24 000 d，其中初始時間步長為0.01 d，增益因子為2，最大運移步長為2 d，模型運行結(jié)果每隔365 d輸出一次。

3 模型結(jié)果

3.1 永興島淡水透鏡體的模擬結(jié)果與檢驗

利用GMS軟件中的SEAWAT模塊對模型進行求解，該模塊由美國地調(diào)局(USGS)研發(fā)推出，是一款用于模擬多孔介質(zhì)中三維非穩(wěn)定變密度流的有限差分模型，其主要是通過同步時間步長方法，即交替使用修改過的MODFLOW和MT3DMS分別求解地下水流運動控制方程和溶質(zhì)運移控制方程，實現(xiàn)水流和溶質(zhì)運移的耦合求解[21]，目前已廣泛應用于海水入侵、島礁淡水透鏡體等各種三維變密度地下水模擬問題的研究中。

模型的計算結(jié)果見圖5和圖6。永興島地貌變化前的淡水透鏡體海平面以上的高度約0～0.41 m，其中最大水頭值為0.41 m，正是由于這部分高出海平面的高度，保持了淡水不斷地流向海洋，將海水向外推移，從而使淡水透鏡體始終處于一個動態(tài)平衡狀態(tài)。若以Cl-濃度0.6 g/L作為淡水透鏡體的邊界濃度[6]，則永興島淡水透鏡體的最大厚度為14 m。

圖5 淡水透鏡體海平面以上的水頭分布及3D形態(tài)Fig.5 Distribution of hydraulic heads above the sea level and the 3D form of the freshwater lens

圖6 地貌變化前B-B’剖面的Cl-濃度分布Fig.6 Cl- concentration distribution of profile along B-B’ before the geomorphologic change

高密度電阻率法是以巖土體的電性差異為基礎(chǔ)，研究在施加電場作用下，地下巖土體傳導電流的變化分布規(guī)律，該方法能夠很好地揭示島上淡水透鏡體的厚度與范圍。Stewart詳細論述了小海島淡水透鏡體電法勘探的方法原理及應用[22]，Comte利用數(shù)值法并結(jié)合電法勘探數(shù)據(jù)探討了降雨補給對淡水透鏡體的影響[13]。因此根據(jù)前人在島上進行高密度電法勘探的研究結(jié)果[4]，結(jié)合6個觀測點的淡水透鏡體厚度對模型結(jié)果進行檢驗，具體觀測點位置見圖1，可以看出模擬結(jié)果基本符合永興島淡水透鏡體的實際形態(tài)。

表2 模擬結(jié)果的識別與檢驗

3.2 參數(shù)的敏感性分析

模型的敏感性分析是為了對已經(jīng)識別過的模型進行不確定性量化，是模型預測中不可缺少的環(huán)節(jié)。通常敏感性分析的做法是在同一時間內(nèi)保證其他值不變的情況下，只改變一個參數(shù)值，觀察其對模擬結(jié)果的影響。但更為有效的敏感性分析是計算敏感度，它表示一個因素的變化對模型計算結(jié)果的影響程度[23]。由于不同模型變量的量綱不同，不同參數(shù)的量綱差別也很大，為了便于不同參數(shù)間敏感度大小的比較，需要進一步無量綱化，通常表示為：

(6)

式中:βi,k——模型變量在觀測點上的敏感度；

Hi——模型變量；

ak——第k個參數(shù)的取值。

其中，βi,k的絕對值越大，表明相應參數(shù)的變化對模擬結(jié)果的影響越大?？紤]到研究區(qū)實際的水文地質(zhì)情況，以及對淡水透鏡體厚度的影響，選取滲透系數(shù)、降雨入滲補給系數(shù)和給水度進行敏感性分析。模型變量為淡水透鏡體的最大厚度。由于淡水透鏡體主要形成于不整合面以上未固結(jié)的珊瑚砂中，滲透系數(shù)一般為40～120 m/d，給水度為0.1～0.2[24]；降雨入滲補給系數(shù)難以確定，前人研究中海島降雨入滲補給系數(shù)多取0.4～0.67[6, 18～19]。經(jīng)過模型參數(shù)識別，以上參數(shù)的初始賦值分別為75 m/d，0.45和0.14，為了便于分析，確定各個參數(shù)在其初始值附近分別增加或減少10%、20%、30%、40%對模型進行敏感度分析，具體的敏感性分析結(jié)果見圖7。從參數(shù)的敏感性分析結(jié)果可以看出：

(1)隨著滲透系數(shù)的增大，淡水透鏡體的厚度逐漸減小。當降雨入滲補給系數(shù)增大時，淡水透鏡體的厚度也隨之增大。給水度對淡水透鏡體的厚度幾乎沒有影響。

(2)島礁淡水透鏡體的厚度對降雨入滲補給系數(shù)的變化最為敏感，其次是含水層的滲透系數(shù)，對含水層給水度的敏感度最小。

(3)各個參數(shù)在其初始賦值附近對稱性增加或減少時，其敏感度大小呈現(xiàn)出對稱性的變化規(guī)律，即隨著參數(shù)的變化幅度增大時，敏感度系數(shù)也隨之增大。

圖7 參數(shù)的敏感性分析Fig.7 Sensitivity analysis of parameters(a)、(b) and (c) The maximum thickness of the freshwater lens varyingwith each parameter；(d) Sensitivity comparison of each parameter

圖8 沙丘面積增加后永興島的淡水透鏡體形成過程(a)沙丘面積增加前的淡水透鏡體形態(tài)；(b)、(c)、(d)、(e)、(f)分別為永興島沙丘面積增加后1年、5年、10年、15年、25年的淡水透鏡體形態(tài)Fig.8 Formation process of the freshwater lens in the Yongxing Island after increase in the dune area(a)state of the freshwater lens before the increase in the dune area; (b)、(c)、(d)、(e) and (f) State of the freshwater lens in 1, 5, 10, 15 and 25 years after the increase in the dune area in the Yongxing Island

3.3 地貌變化后永興島淡水透鏡體的模擬結(jié)果

根據(jù)模型的計算結(jié)果，地貌變化后的永興島淡水透鏡體隨時間的變化見圖8，可以看出在之前5年時間內(nèi)，隨著島上降雨的不斷入滲補給，永興島這部分新增區(qū)域地下水體的離子濃度不斷得到淡化，但這一階段還沒有淡水體的產(chǎn)生，可以稱之為淡水透鏡體形成的準備階段；之后隨著時間的不斷推移，包括新增區(qū)域在內(nèi)，整個永興島淡水透鏡體的厚度及儲量都得到了增加。大約在25年后能夠再次達到穩(wěn)定狀態(tài)。

再次穩(wěn)定時的永興島淡水透鏡體在形態(tài)上明顯不對稱，即永興島西南區(qū)域的淡水透鏡體厚度較大，最大深度達到了17.5 m，而永興島東北區(qū)域淡水透鏡體的最大厚度僅為8.5 m(圖8)。初步分析認為，該區(qū)域的沙丘面積較小，且平均寬度僅為685 m，接受降雨入滲補給的面積有限，根據(jù)前面的敏感性分析可知，降雨的入滲補給對島礁淡水透鏡體厚度的影響最大；過窄的區(qū)域?qū)挾葧е陆涤耆霛B補給的淡水很快流入海洋，故而使得再次穩(wěn)定時的永興島淡水透鏡體產(chǎn)生了明顯的不對稱現(xiàn)象。

穩(wěn)定時淡水透鏡體的水頭分布見圖9。淡水透鏡體中央處最大水頭值為0.49 m，根據(jù)淡水透鏡體的貯量計算公式[6]：

(7)

式中:V——貯水量/m3；

μ——含水層的給水度；

S——平面上每個差分網(wǎng)格的面積/m2；

hf,i,j——每個網(wǎng)格的淡水厚度/m。

圖9 25年后永興島海平面以上的水頭分布Fig.9 Distribution of hydraulic head abov the sea level in theYongxing Island 25 years later

計算出25年后的永興島淡水透鏡體的貯量約為328×104m3。雖然島嶼面積僅增加了50%，但與地貌變化前的淡水儲量相比，沙丘面積增加后的淡水儲量增加了近一倍，這是因為隨著島礁面積的增大，不僅新增區(qū)域的海平面以下形成了淡水，原先的島礁淡水透鏡體厚度也隨之增加。

4 結(jié)論

(1)永興島北部新增的這部分區(qū)域主要由人為因素搬運堆積的珊瑚砂礫等生物碎屑組成，具有明顯的顆粒級配不良且不連續(xù)的特征，結(jié)構(gòu)疏松且孔滲性良好，有利于島礁淡水透鏡體的形成。

(2)結(jié)合永興島的地貌地層條件，建立了永興島淡水透鏡體的三維數(shù)值模型，計算得出地貌變化前的永興島淡水透鏡體的最大厚度達14 m，其中海平面以上的高度為0～0.41 m，將模擬結(jié)果與島上高密度電阻率法的觀測數(shù)據(jù)進行對比擬合，模型能夠很好地反映淡水透鏡體的特征。

(3)對影響島礁淡水透鏡體形成的主要參數(shù)，如滲透系數(shù)，降雨的入滲補給系數(shù)，給水度進行了敏感性分析，可以看出，降雨入滲補給系數(shù)對島礁淡水透鏡體形成的影響最大，其次是含水層的滲透系數(shù)，而含水層的給水度對島礁淡水透鏡體的形成影響最小。

(4)在保持現(xiàn)有島礁面積不變的情況下，永興島大約在25年后再次形成穩(wěn)定的淡水透鏡體，此時淡水透鏡體的最大厚度為17.5 m，其中島礁東北部區(qū)域淡化水體的最大厚度為8.5 m；計算所得的淡水資源儲量增加了近一倍，約328×104m3。

(5)淡水透鏡體再次穩(wěn)定時，其東北區(qū)域的厚度較小，這主要是由于該區(qū)域的沙丘面積和寬度較小；進一步可以認為若沙丘面積和寬度過小，如面積增加前的永興島東北區(qū)域，降雨入滲補給的淡水會很快流入海洋，難以形成淡水透鏡體。