基巖島嶼地下水數值模擬發(fā)展研究現狀

王大慶， 許顥礫， 鄧正棟， 丁志斌， 倪博睿， 周澤林， 趙小蘭

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

0 引言

基巖島嶼多數為大陸島[1-2]，是大陸的延伸，其基巖多數為花崗巖、片麻巖(如我國東海、南海的近岸大陸島)； 基巖島嶼也有的是大洋島[1-2]，或由海底火山噴發(fā)形成的火山島，其基巖多數為玄武巖、安山巖(如琉球群島、硫磺列島)?；鶐r島嶼的地下水埋藏或富集方式有別于沙質島嶼[3]。沙質島嶼的淡水體以“淡水透鏡體”的形式存在[4-5]，其地下水一般經歷由咸變淡的過程，且最終形成淡水與咸水的分界面。本文所關注的基巖島嶼演變過程多數是由陸變?yōu)閸u的過程，其地下淡水資源不必經過由咸變淡的漫長過程，地下淡水以淡水蘑菇體的形式存在[6]。本文通過綜述基巖島嶼地下水的相關理論及數值模擬發(fā)展研究現狀，為今后開展基巖島嶼地下水數值模擬研究提供理論依據。

1 基巖島嶼地下水模擬相關理論

1.1 島嶼淡水的存在形式

基巖島嶼和沙質島嶼是2種不同地質巖性的島嶼，其地下淡水的存在形式也不相同(圖1)。基巖島嶼淡水呈“蘑菇體”，沙質島嶼淡水呈“透鏡體”。由于基巖島嶼原本就存在淡水，不必經歷由咸變淡的過程，加之基巖的滲透系數較小[7]，所含地下淡水不是以“透鏡體”的形式浮在咸水上，且淡水底部沒有咸淡水過渡帶。

圖1 基巖島嶼淡水“蘑菇體”(左)與沙質島嶼淡水“透鏡體”(右)

Fig.1Freshwatermushroominthebedrockislands(left)andfreshwaterlensinthesandyislands(right)

1.2 基巖島嶼地下水的均衡

1.2.1 基巖島嶼地下水的補給和排泄

島嶼地下淡水的主要補給來自大氣降雨入滲。一些近大陸島嶼的中深層或深層含水層與大陸的含水層水文地質條件一致，與大陸地下含水層有水力聯系[8]，島嶼中深層或深層的含水層地下水補給也來源于大陸地下含水層的補給。島嶼地下水其他可能的補給源還有人造水庫排放、天然水坑或采石場采石形成的水坑補給，這些補給源本質上仍然是大氣降雨。

基巖島嶼的排泄主要是地面蒸發(fā)，其次是泉水出露； 另外，地下水也會向有水力聯系的地面水坑或小湖泊進行排泄； 還有一小部分通過植物的葉面蒸騰排泄，但蒸騰量不大[9]。

1.2.2 島嶼地下水均衡方程

島嶼的地下水均衡與大陸地下水均衡相似，其中，潛水水量均衡方程為

μΔh=(X+Y+W1+Z1+R1)-

(W2+Ws+Z2+R2)，

(1)

式中：μΔh為潛水變化量，m3；X為降水入滲量，m3；Y為地表水對潛水的補給量，m3；W1，W2分別是地表水流入與流出量，m3；Z1，Z2分別為潛水的凝結水補給量和蒸發(fā)量，m3；Ws為泉涌出量，m3；R1，R2為人工補給量與排出量，m3。

承壓水水量均衡方程為

(2)

島嶼水均衡方程為

Q補-Q排=Q降雨+Q陸補+Q人工+Q地補-Q蒸發(fā)-Q泉-Q地排-Q蒸騰-Q開采，

(3)

式中：Q補-Q排為水量變化量，m3；Q降雨為降雨量，m3；Q陸補為大陸水力聯系補給量，m3；Q人工為人工水庫等排泄量，m3；Q地補為地表水坑或湖泊補給量，m3；Q蒸發(fā)為蒸發(fā)量，m3；Q泉為泉水涌出量，m3；Q地排為向地表水坑或湖泊排泄量，m3；Q蒸騰為植物蒸騰量，m3；Q開采為人工開采水量，m3。

1.3 基巖島嶼地下水的運動規(guī)律

不同巖性基巖島嶼的地下水運動理論表述形式有所不同，具體理論表述應根據島嶼具體的基巖地層巖性選擇。

1.3.1 松散地層基巖島嶼地下水運動規(guī)律

若基巖島嶼地層為松散地層，如鎮(zhèn)江市東海島[10]，島嶼地層主要為受風化、卸荷、地形因素影響的花崗巖、片麻巖等，這些地質體可以概化為非均質各向異性或橫向各項同性介質，用地下水運動方程(簡化的各向異性介質)或二維流運動方程描述這些地層中的地下水運動規(guī)律為

(4)

式中：K為滲透系數，m/d；H為測壓管水頭，m；Ss為貯水率，m-1；t為時間，s；x,y,z表示分布的不同方向。

式(4)為一般的各向異性介質的地下水運動方程[11]，將x=y代入式(4)中，得出橫向各項同性介質的地下水運動方程，即

(5)

式中：K為滲透系數，m/d；H為測壓管水頭，m；Ss為貯水率，m-1；t為時間，s；x,y,z分布表示不同的方向。

設Kx為水平方向的滲透系數，Kz為豎直方向的滲透系數，則二維流運動方程為

(6)

式中：K為滲透系數，m/d；H為測壓管水頭，m；Ss為貯水率，m-1；t為時間，s；x,y,z分布表示不同的方向。

其實三維流和二維流本質上是一致的(圖2)。

圖2 二維、三維滲流模型地質體模型

圖2中剖面1為地質體y-z剖面，剖面2為地質體x-z剖面，其滲透系數分別為Ky和Kz，Kx和Kz。而Ky=Kx，則說明y-z剖面與x-z剖面滲透系數一致，也就說明了3D形式的滲透系數與2D的滲透系數一致，即3D流與2D流的本質一致，僅僅是表示形式不同。

1.3.2 裂隙巖體基巖島嶼地下水的運動規(guī)律

若基巖島嶼以基巖裂隙巖體為主，則應考慮用立方定律，根據納維-斯托克斯(Navier-Stokes equation)方程簡化可得在光滑平板裂隙中水體運動的理論公式為

(7)

式中：u為平均流速，m/s；g為重力加速度，m/s2；δ是裂隙的寬度，m；v水的運動黏滯系數，m2/s；J則是水力坡度；Kf是裂隙滲透系數，m/s。

2 基巖島嶼地下水的模擬現狀

2.1 國內外基巖島嶼地下水模擬

對島嶼地下水流場模擬是為了對島嶼地下水的分布進行預測，從而指導對島嶼地下水的開采、利用及保護。1995年，李國敏等[5]運用水頭、濃度相互依賴的有限元法對潿洲島的火山碎屑孔隙含水層及玄武巖孔洞-裂隙含水層的水位、水流、水質等參數進行了計算和模擬分析，為潿洲島地下水的開采規(guī)劃與控制海水入侵提供了依據。1998年，Person等[12]利用sharp-interface模型對楠塔基島(美國)在1977—2020年間的地下水分布情況進行研究和預測，并討論了海水入侵的情況。隨著計算機軟件的發(fā)展，一些相應的地下水模擬程序及軟件相繼問世，Liu等[13]于2006年根據中國金門島的地質條件(該島東半部以花崗巖、片麻巖為主，而西半部則以紅土為主)設定不同的參數，用MODFLOW-96模擬了金門島的地下水分布，評估了地下水水位變化等級，提供了金門島地下水的利用方案。2008年Barazzuoli等[14]以意大利托斯卡納濱海地區(qū)1995—2013年的實地與實驗數據為依據，利用有限元FEFLOW數值模擬軟件建模分析了地下水自然流入海洋、河流和人工水井的情況，考慮海水入侵的影響，進行了該地區(qū)用水的預案。溫漢輝[15]于2013年在對雷州半島淺層含水層(松散巖類和火山巖空洞裂隙巖層組成)進行概化的基礎上，采用GMS數值軟件前后處理及MODFLOW模塊對該地區(qū)地下水進行了模擬研究。EI-Kadi等[16]在2014年利用GMS數值軟件對濟州島(該島是火山活動形成，玄武巖和火山碎屑巖為主)的地下水情況進行模擬，用MODFLOW模塊模擬分析地下水水位及泉水等，用SEAWATER模塊來模擬分析海水入侵情況，根據水位海拔、泉流、鹽度等指標，評估了該島地下水可持續(xù)的能力。2016年，Lathashri等[17]利用GMS數值模擬軟件的MODFLOW模塊和SEAWATER模塊來研究濱海地區(qū)的地下水。Yi等[18]采用FEMWATER模塊以變密度流體概化地下水與海水，利用GMS模擬天津的渤海灣地下水的排放和海水入侵問題。騰建標等[19]、Zhou等[20]先后對湛江的東海島進行了研究，提出該島嶼的中深層含水層和深層含水層與大陸地下水存在水力聯系，將該島嶼概化為非均質各向異性，利用Visual Modflow地下水數值模擬軟件進行模擬求解，再結合實際鉆孔資料對比研究，所建模型可以較好地反映東海島地下水水文地質情況，根據模擬結果討論了該島的地下水可持續(xù)性情況，并用于指導開采使用地下水。

基巖島嶼地下水模擬在一些集成化的程序及軟件未出現前，以解析計算為主流，其可視化效果不好，且對模擬者的編程水平要求很高。隨著MODFLOW、Visual Modflow、GMS和FEFLOW等軟件相繼出現，地下水模擬變得簡單化、流程化，模擬過程中，需要對基巖島嶼的地下水補徑排及海水入侵情況進行詳細調查，再通過調查試驗獲得研究區(qū)地質參數，然后利用軟件建立模型即可進行地下水模擬。對于火山碎屑巖、玄武巖孔隙介質、風化程度較高的花崗巖或片麻巖及其他松散巖類的基巖島嶼的概化是沒有問題的，符合達西定律和地下水運動方程或二維流運動方程(式(4)～(6))。但對于一些以裂隙基巖為主的基巖島嶼，例如花崗巖、片麻巖，其地下水的運動應主要符合立方定律(式(7))，若仍然要概化為孔隙模型，類似于等效多孔連續(xù)介質模型[21-22]，則需要對局部地區(qū)的巖性、給水度、滲透系數等參數進行相應修改，如：將滲透系數設大以達到類似于裂隙流的目的，隨后還必須進行精度驗證和模型校正。若直接考慮為裂隙模型，則需要用離散裂隙網絡滲流模型[23-24]來模擬裂隙水流運動。

若基巖裂隙為主的島嶼淺部地層風化程度高，已達到松散巖類的程度，其中，深部地層仍然是未風化的基巖裂隙，則應該考慮雙層介質或多層介質。如果用孔隙模型模擬，可用Visual Modflow、GMS和FEFLOW等軟件進行模擬，則所設定的分層參數應該符合各個地層的性質； 如果用裂隙模型模擬，則可用TOUGH2軟件或基于GMS的三維TOUGH2來模擬裂隙[25]。目前，美國勞倫斯伯克利實驗室推出了TOUGH3軟件，是TOUGH2的升級版，增加了編碼修改功能，實現了更大范圍的模擬。MODFLOW-SURFACT也是一個可以模擬雙重介質的程序模塊，可以模擬裂隙井的模塊[26]。CONNETFLOW是一款針對裂隙的有限元地下水模擬軟件，其中NAMMU模塊針對多孔介質模型，NAPSAC模塊針對離散裂隙網絡滲流模型，并且該軟件可以應用于雙重介質。

2.2 基巖島嶼海水入侵模擬

1987年，Voss等[27]用SUTRA軟件(二維有限元軟件)模擬夏威夷Dahu島南部的海水入侵。Boxton等[28]于1992年利用MODFLOW地下水流動程序與粒子追蹤程序模擬了紐約長島的海水入侵情況。陳耀登等[29]利用MODFLOW軟件，借鑒湖泊地下水模擬的處理方式，將海岸邊界設為具有周期性和高滲透性的定水頭邊界也能在一定程度上反映地下水的變化情況。盧薇等[30]利用FEFLOW軟件建立珠江口東岸地區(qū)海水入侵三維溶質模型,通過設定合理的濃度場、流場參數模擬淡水與海水，主要考慮了Clˉ濃度。陳開榮等[31]利用Visual Modflow軟件的SEAWATER模塊模擬了大連海洋大學新校區(qū)變密度條件下的海水入侵。SEAWATER模塊是基于MODFLOW與MT3DMS的多孔介質三維非穩(wěn)定變密度地下水流及溶質運移模擬模型，已被用于各種變密度地下水模擬及海水入侵的模擬問題。Lu等[32]利用FEFLOW軟件針對中國深圳濱海地區(qū)建立3D變密度模型研究海水入侵，并研究了潮汐的影響?？梢?，在模擬海水入侵問題上，SUTRA軟件是較早的針對海水入侵模擬的二維軟件，三維軟件Visual MODFLOW軟件與GMS軟件可以考慮Clˉ濃度(或其他礦化度，Brˉ、鈉離子吸附比，咸化系數[33]也能反映海水與淡水的差異，也可當作判定指標),利用其MT3DMS模塊來模擬海水。近些年推出的SEAWATER模塊專門針對變密度，且基于MODFLOW的水流運動模型與MT3DMS的溶質運移模型，利用軟件的模塊的相互結合可以解決變水頭、變密度、變濃度(即水流場、密度場、濃度場的耦合)的問題，實現島嶼淡水與海水間的模擬，實現島嶼海水入侵模擬。

海水模擬的一大關鍵是如何設定海水邊界。對于基巖島嶼而言，可以直接將海岸線作為含水系統(tǒng)的海水邊界[34]。為了解基巖島嶼淡水的存在形式，龐忠和[35]利用物探方法對廟島群島的南長島北部進行了電阻率測深與電剖面測定，根據海水、淡水的電阻率不同，再結合鉆孔資料，可以看出海水、淡水邊界幾乎為該島嶼的海岸線。對于基巖島嶼將海水邊界設定在島嶼海岸線也符合圖1(左)所示的淡水“蘑菇體”的基巖島嶼淡水的貯存形式。若要更為精準地考慮定邊界問題，則需要利用物探方法結合鉆孔數據來確定。

綜上，有關基巖島嶼海水入侵的研究實例可見，在基巖島嶼海水入侵數值模擬時，應該將基巖海岸線設置為隔水邊界； 而基巖島嶼的海岸出現砂質地質體部分(如沙灘)時，在數值模擬時應將這類海岸設置為定水頭邊界，此時，與沙質島嶼關于海水入侵的邊界模擬相似。

2.3 基巖島嶼地下水污染物模擬

由于人工開發(fā)島嶼資源，島嶼地下水可能受到以下污染：菜地、廁所、餐廳等人們的生活污水污染，采石場的礦渣、廢石等殘渣污染，發(fā)電廠機器運轉的廢油及污水排放污染，垃圾廠的垃圾溶濾滲漏溶液的污染，污水處理廠的污水排放污染，以及海水入侵的影響[36]等等。目前，多數研究者仍然延用鄭春苗等[37-38]開發(fā)的MT3D、MT3DMS程序,如：Curtis等[39]于2006年利用MODFLOW模擬水流場，用MT3DMS程序模擬鈾(VI)的反應運移； 劉娟等[40]利用FEFLOW模擬地下水中四氯化碳污染羽的濃度場; 長安大學高小文等[41]利用MT3DMS溶質運移模塊模擬尾礦庫周邊環(huán)境污染，為環(huán)境治理提供參考依據。可見，目前污染物運移模擬依舊利用三大熱門的地下水模擬軟件，主要還是依托于MT3DMS溶質運移模塊。

3 認識與建議

3.1 認識

(1)基巖島嶼多數為大陸島，其地下水以淡水“蘑菇體”的形式存在，無需經歷由咸變淡的過程，且無咸淡水過渡帶。

(2)基巖島嶼地下水的補給主要是降雨，其次是島嶼中、深層含水層與大陸的水力聯系補給。

(3)基巖島嶼的地質模型可概化為3種：孔隙型、裂隙型和孔隙-裂隙型?；鹕綅u、松散巖類或風化程度較高的基巖島嶼，可概化為孔隙型，其滲流符合非均質各向異性或橫向各項同性的地下水運動方程，可直接用Visual Modflow、GMS和FEFLOW等軟件模擬； 島嶼地質體以裂隙基巖為主，應概化為裂隙型，其地下水運動規(guī)律符合立方定律，可采用MODFLOW-SURFACT模塊、TOUGH軟件、CONNETFLOW軟件進行模擬； 對于孔隙-裂隙型模型，則需基于詳細的地質體建模，結合TOUGH軟件或SURFACT程序模塊等軟件進行耦合模擬。

(4)對于基巖島嶼的海水邊界設定，應根據島嶼海岸線及物探與鉆探資料圈定出海水邊界。

3.2 建議

(1)基巖島嶼地下水數值模擬需做到精細化，要將地質體的各種參數刻畫細致，如每種巖層的給水度、滲透系數等。還應細致刻畫，包括湖泊、河流等的地表水，對于降雨補給、蒸發(fā)排泄、潮汐影響以及氣候引起的水量變化也需仔細確定。此外，還應注意與大陸的中、深層含水層的水力聯系的補給項。

(2)基巖島嶼多數為孔隙與裂隙并存巖體，其地下水賦存地質模型需考慮雙重或多重介質耦合模型。目前，TOUGH軟件、MODFLOW-SURFACT程序模塊或CONNETFLOW軟件可實現雙重或多重介質的模擬，但該方面的文獻較少，仍然是一個研究難點，希望廣大研究者可以研發(fā)出更完美的模擬軟件。

(3)對于基巖島嶼地下水模擬的驗證也是重中之重，數值模擬結果應與島嶼水量均衡，且與實測水量、水位相一致，模擬是一個不斷修正與完善的過程，不能一蹴而就。