王大慶, 許顥礫, 鄧正棟, 丁志斌, 倪博睿, 周澤林, 趙小蘭
(1.陸軍工程大學(xué)國防工程學(xué)院,南京 210007; 2.中國礦業(yè)大學(xué)資源與地球科學(xué)學(xué)院,徐州 221116)
基巖島嶼多數(shù)為大陸島[1-2],是大陸的延伸,其基巖多數(shù)為花崗巖、片麻巖(如我國東海、南海的近岸大陸島); 基巖島嶼也有的是大洋島[1-2],或由海底火山噴發(fā)形成的火山島,其基巖多數(shù)為玄武巖、安山巖(如琉球群島、硫磺列島)。基巖島嶼的地下水埋藏或富集方式有別于沙質(zhì)島嶼[3]。沙質(zhì)島嶼的淡水體以“淡水透鏡體”的形式存在[4-5],其地下水一般經(jīng)歷由咸變淡的過程,且最終形成淡水與咸水的分界面。本文所關(guān)注的基巖島嶼演變過程多數(shù)是由陸變?yōu)閸u的過程,其地下淡水資源不必經(jīng)過由咸變淡的漫長過程,地下淡水以淡水蘑菇體的形式存在[6]。本文通過綜述基巖島嶼地下水的相關(guān)理論及數(shù)值模擬發(fā)展研究現(xiàn)狀,為今后開展基巖島嶼地下水?dāng)?shù)值模擬研究提供理論依據(jù)。
基巖島嶼和沙質(zhì)島嶼是2種不同地質(zhì)巖性的島嶼,其地下淡水的存在形式也不相同(圖1)?;鶐r島嶼淡水呈“蘑菇體”,沙質(zhì)島嶼淡水呈“透鏡體”。由于基巖島嶼原本就存在淡水,不必經(jīng)歷由咸變淡的過程,加之基巖的滲透系數(shù)較小[7],所含地下淡水不是以“透鏡體”的形式浮在咸水上,且淡水底部沒有咸淡水過渡帶。
圖1 基巖島嶼淡水“蘑菇體”(左)與沙質(zhì)島嶼淡水“透鏡體”(右)
Fig.1Freshwatermushroominthebedrockislands(left)andfreshwaterlensinthesandyislands(right)
1.2.1 基巖島嶼地下水的補(bǔ)給和排泄
島嶼地下淡水的主要補(bǔ)給來自大氣降雨入滲。一些近大陸島嶼的中深層或深層含水層與大陸的含水層水文地質(zhì)條件一致,與大陸地下含水層有水力聯(lián)系[8],島嶼中深層或深層的含水層地下水補(bǔ)給也來源于大陸地下含水層的補(bǔ)給。島嶼地下水其他可能的補(bǔ)給源還有人造水庫排放、天然水坑或采石場采石形成的水坑補(bǔ)給,這些補(bǔ)給源本質(zhì)上仍然是大氣降雨。
基巖島嶼的排泄主要是地面蒸發(fā),其次是泉水出露; 另外,地下水也會(huì)向有水力聯(lián)系的地面水坑或小湖泊進(jìn)行排泄; 還有一小部分通過植物的葉面蒸騰排泄,但蒸騰量不大[9]。
1.2.2 島嶼地下水均衡方程
島嶼的地下水均衡與大陸地下水均衡相似,其中,潛水水量均衡方程為
μΔh=(X+Y+W1+Z1+R1)-
(W2+Ws+Z2+R2),
(1)
式中:μΔh為潛水變化量,m3;X為降水入滲量,m3;Y為地表水對潛水的補(bǔ)給量,m3;W1,W2分別是地表水流入與流出量,m3;Z1,Z2分別為潛水的凝結(jié)水補(bǔ)給量和蒸發(fā)量,m3;Ws為泉涌出量,m3;R1,R2為人工補(bǔ)給量與排出量,m3。
承壓水水量均衡方程為
(2)
島嶼水均衡方程為
Q補(bǔ)-Q排=Q降雨+Q陸補(bǔ)+Q人工+Q地補(bǔ)-Q蒸發(fā)-Q泉-Q地排-Q蒸騰-Q開采,
(3)
式中:Q補(bǔ)-Q排為水量變化量,m3;Q降雨為降雨量,m3;Q陸補(bǔ)為大陸水力聯(lián)系補(bǔ)給量,m3;Q人工為人工水庫等排泄量,m3;Q地補(bǔ)為地表水坑或湖泊補(bǔ)給量,m3;Q蒸發(fā)為蒸發(fā)量,m3;Q泉為泉水涌出量,m3;Q地排為向地表水坑或湖泊排泄量,m3;Q蒸騰為植物蒸騰量,m3;Q開采為人工開采水量,m3。
不同巖性基巖島嶼的地下水運(yùn)動(dòng)理論表述形式有所不同,具體理論表述應(yīng)根據(jù)島嶼具體的基巖地層巖性選擇。
1.3.1 松散地層基巖島嶼地下水運(yùn)動(dòng)規(guī)律
若基巖島嶼地層為松散地層,如鎮(zhèn)江市東海島[10],島嶼地層主要為受風(fēng)化、卸荷、地形因素影響的花崗巖、片麻巖等,這些地質(zhì)體可以概化為非均質(zhì)各向異性或橫向各項(xiàng)同性介質(zhì),用地下水運(yùn)動(dòng)方程(簡化的各向異性介質(zhì))或二維流運(yùn)動(dòng)方程描述這些地層中的地下水運(yùn)動(dòng)規(guī)律為
(4)
式中:K為滲透系數(shù),m/d;H為測壓管水頭,m;Ss為貯水率,m-1;t為時(shí)間,s;x,y,z表示分布的不同方向。
式(4)為一般的各向異性介質(zhì)的地下水運(yùn)動(dòng)方程[11],將x=y代入式(4)中,得出橫向各項(xiàng)同性介質(zhì)的地下水運(yùn)動(dòng)方程,即
(5)
式中:K為滲透系數(shù),m/d;H為測壓管水頭,m;Ss為貯水率,m-1;t為時(shí)間,s;x,y,z分布表示不同的方向。
設(shè)Kx為水平方向的滲透系數(shù),Kz為豎直方向的滲透系數(shù),則二維流運(yùn)動(dòng)方程為
(6)
式中:K為滲透系數(shù),m/d;H為測壓管水頭,m;Ss為貯水率,m-1;t為時(shí)間,s;x,y,z分布表示不同的方向。
其實(shí)三維流和二維流本質(zhì)上是一致的(圖2)。
圖2 二維、三維滲流模型地質(zhì)體模型
圖2中剖面1為地質(zhì)體y-z剖面,剖面2為地質(zhì)體x-z剖面,其滲透系數(shù)分別為Ky和Kz,Kx和Kz。而Ky=Kx,則說明y-z剖面與x-z剖面滲透系數(shù)一致,也就說明了3D形式的滲透系數(shù)與2D的滲透系數(shù)一致,即3D流與2D流的本質(zhì)一致,僅僅是表示形式不同。
1.3.2 裂隙巖體基巖島嶼地下水的運(yùn)動(dòng)規(guī)律
若基巖島嶼以基巖裂隙巖體為主,則應(yīng)考慮用立方定律,根據(jù)納維-斯托克斯(Navier-Stokes equation)方程簡化可得在光滑平板裂隙中水體運(yùn)動(dòng)的理論公式為
(7)
式中:u為平均流速,m/s;g為重力加速度,m/s2;δ是裂隙的寬度,m;v水的運(yùn)動(dòng)黏滯系數(shù),m2/s;J則是水力坡度;Kf是裂隙滲透系數(shù),m/s。
對島嶼地下水流場模擬是為了對島嶼地下水的分布進(jìn)行預(yù)測,從而指導(dǎo)對島嶼地下水的開采、利用及保護(hù)。1995年,李國敏等[5]運(yùn)用水頭、濃度相互依賴的有限元法對潿洲島的火山碎屑孔隙含水層及玄武巖孔洞-裂隙含水層的水位、水流、水質(zhì)等參數(shù)進(jìn)行了計(jì)算和模擬分析,為潿洲島地下水的開采規(guī)劃與控制海水入侵提供了依據(jù)。1998年,Person等[12]利用sharp-interface模型對楠塔基島(美國)在1977—2020年間的地下水分布情況進(jìn)行研究和預(yù)測,并討論了海水入侵的情況。隨著計(jì)算機(jī)軟件的發(fā)展,一些相應(yīng)的地下水模擬程序及軟件相繼問世,Liu等[13]于2006年根據(jù)中國金門島的地質(zhì)條件(該島東半部以花崗巖、片麻巖為主,而西半部則以紅土為主)設(shè)定不同的參數(shù),用MODFLOW-96模擬了金門島的地下水分布,評估了地下水水位變化等級,提供了金門島地下水的利用方案。2008年Barazzuoli等[14]以意大利托斯卡納濱海地區(qū)1995—2013年的實(shí)地與實(shí)驗(yàn)數(shù)據(jù)為依據(jù),利用有限元FEFLOW數(shù)值模擬軟件建模分析了地下水自然流入海洋、河流和人工水井的情況,考慮海水入侵的影響,進(jìn)行了該地區(qū)用水的預(yù)案。溫漢輝[15]于2013年在對雷州半島淺層含水層(松散巖類和火山巖空洞裂隙巖層組成)進(jìn)行概化的基礎(chǔ)上,采用GMS數(shù)值軟件前后處理及MODFLOW模塊對該地區(qū)地下水進(jìn)行了模擬研究。EI-Kadi等[16]在2014年利用GMS數(shù)值軟件對濟(jì)州島(該島是火山活動(dòng)形成,玄武巖和火山碎屑巖為主)的地下水情況進(jìn)行模擬,用MODFLOW模塊模擬分析地下水水位及泉水等,用SEAWATER模塊來模擬分析海水入侵情況,根據(jù)水位海拔、泉流、鹽度等指標(biāo),評估了該島地下水可持續(xù)的能力。2016年,Lathashri等[17]利用GMS數(shù)值模擬軟件的MODFLOW模塊和SEAWATER模塊來研究濱海地區(qū)的地下水。Yi等[18]采用FEMWATER模塊以變密度流體概化地下水與海水,利用GMS模擬天津的渤海灣地下水的排放和海水入侵問題。騰建標(biāo)等[19]、Zhou等[20]先后對湛江的東海島進(jìn)行了研究,提出該島嶼的中深層含水層和深層含水層與大陸地下水存在水力聯(lián)系,將該島嶼概化為非均質(zhì)各向異性,利用Visual Modflow地下水?dāng)?shù)值模擬軟件進(jìn)行模擬求解,再結(jié)合實(shí)際鉆孔資料對比研究,所建模型可以較好地反映東海島地下水水文地質(zhì)情況,根據(jù)模擬結(jié)果討論了該島的地下水可持續(xù)性情況,并用于指導(dǎo)開采使用地下水。
基巖島嶼地下水模擬在一些集成化的程序及軟件未出現(xiàn)前,以解析計(jì)算為主流,其可視化效果不好,且對模擬者的編程水平要求很高。隨著MODFLOW、Visual Modflow、GMS和FEFLOW等軟件相繼出現(xiàn),地下水模擬變得簡單化、流程化,模擬過程中,需要對基巖島嶼的地下水補(bǔ)徑排及海水入侵情況進(jìn)行詳細(xì)調(diào)查,再通過調(diào)查試驗(yàn)獲得研究區(qū)地質(zhì)參數(shù),然后利用軟件建立模型即可進(jìn)行地下水模擬。對于火山碎屑巖、玄武巖孔隙介質(zhì)、風(fēng)化程度較高的花崗巖或片麻巖及其他松散巖類的基巖島嶼的概化是沒有問題的,符合達(dá)西定律和地下水運(yùn)動(dòng)方程或二維流運(yùn)動(dòng)方程(式(4)~(6))。但對于一些以裂隙基巖為主的基巖島嶼,例如花崗巖、片麻巖,其地下水的運(yùn)動(dòng)應(yīng)主要符合立方定律(式(7)),若仍然要概化為孔隙模型,類似于等效多孔連續(xù)介質(zhì)模型[21-22],則需要對局部地區(qū)的巖性、給水度、滲透系數(shù)等參數(shù)進(jìn)行相應(yīng)修改,如:將滲透系數(shù)設(shè)大以達(dá)到類似于裂隙流的目的,隨后還必須進(jìn)行精度驗(yàn)證和模型校正。若直接考慮為裂隙模型,則需要用離散裂隙網(wǎng)絡(luò)滲流模型[23-24]來模擬裂隙水流運(yùn)動(dòng)。
若基巖裂隙為主的島嶼淺部地層風(fēng)化程度高,已達(dá)到松散巖類的程度,其中,深部地層仍然是未風(fēng)化的基巖裂隙,則應(yīng)該考慮雙層介質(zhì)或多層介質(zhì)。如果用孔隙模型模擬,可用Visual Modflow、GMS和FEFLOW等軟件進(jìn)行模擬,則所設(shè)定的分層參數(shù)應(yīng)該符合各個(gè)地層的性質(zhì); 如果用裂隙模型模擬,則可用TOUGH2軟件或基于GMS的三維TOUGH2來模擬裂隙[25]。目前,美國勞倫斯伯克利實(shí)驗(yàn)室推出了TOUGH3軟件,是TOUGH2的升級版,增加了編碼修改功能,實(shí)現(xiàn)了更大范圍的模擬。MODFLOW-SURFACT也是一個(gè)可以模擬雙重介質(zhì)的程序模塊,可以模擬裂隙井的模塊[26]。CONNETFLOW是一款針對裂隙的有限元地下水模擬軟件,其中NAMMU模塊針對多孔介質(zhì)模型,NAPSAC模塊針對離散裂隙網(wǎng)絡(luò)滲流模型,并且該軟件可以應(yīng)用于雙重介質(zhì)。
1987年,Voss等[27]用SUTRA軟件(二維有限元軟件)模擬夏威夷Dahu島南部的海水入侵。Boxton等[28]于1992年利用MODFLOW地下水流動(dòng)程序與粒子追蹤程序模擬了紐約長島的海水入侵情況。陳耀登等[29]利用MODFLOW軟件,借鑒湖泊地下水模擬的處理方式,將海岸邊界設(shè)為具有周期性和高滲透性的定水頭邊界也能在一定程度上反映地下水的變化情況。盧薇等[30]利用FEFLOW軟件建立珠江口東岸地區(qū)海水入侵三維溶質(zhì)模型,通過設(shè)定合理的濃度場、流場參數(shù)模擬淡水與海水,主要考慮了Clˉ濃度。陳開榮等[31]利用Visual Modflow軟件的SEAWATER模塊模擬了大連海洋大學(xué)新校區(qū)變密度條件下的海水入侵。SEAWATER模塊是基于MODFLOW與MT3DMS的多孔介質(zhì)三維非穩(wěn)定變密度地下水流及溶質(zhì)運(yùn)移模擬模型,已被用于各種變密度地下水模擬及海水入侵的模擬問題。Lu等[32]利用FEFLOW軟件針對中國深圳濱海地區(qū)建立3D變密度模型研究海水入侵,并研究了潮汐的影響。可見,在模擬海水入侵問題上,SUTRA軟件是較早的針對海水入侵模擬的二維軟件,三維軟件Visual MODFLOW軟件與GMS軟件可以考慮Clˉ濃度(或其他礦化度,Brˉ、鈉離子吸附比,咸化系數(shù)[33]也能反映海水與淡水的差異,也可當(dāng)作判定指標(biāo)),利用其MT3DMS模塊來模擬海水。近些年推出的SEAWATER模塊專門針對變密度,且基于MODFLOW的水流運(yùn)動(dòng)模型與MT3DMS的溶質(zhì)運(yùn)移模型,利用軟件的模塊的相互結(jié)合可以解決變水頭、變密度、變濃度(即水流場、密度場、濃度場的耦合)的問題,實(shí)現(xiàn)島嶼淡水與海水間的模擬,實(shí)現(xiàn)島嶼海水入侵模擬。
海水模擬的一大關(guān)鍵是如何設(shè)定海水邊界。對于基巖島嶼而言,可以直接將海岸線作為含水系統(tǒng)的海水邊界[34]。為了解基巖島嶼淡水的存在形式,龐忠和[35]利用物探方法對廟島群島的南長島北部進(jìn)行了電阻率測深與電剖面測定,根據(jù)海水、淡水的電阻率不同,再結(jié)合鉆孔資料,可以看出海水、淡水邊界幾乎為該島嶼的海岸線。對于基巖島嶼將海水邊界設(shè)定在島嶼海岸線也符合圖1(左)所示的淡水“蘑菇體”的基巖島嶼淡水的貯存形式。若要更為精準(zhǔn)地考慮定邊界問題,則需要利用物探方法結(jié)合鉆孔數(shù)據(jù)來確定。
綜上,有關(guān)基巖島嶼海水入侵的研究實(shí)例可見,在基巖島嶼海水入侵?jǐn)?shù)值模擬時(shí),應(yīng)該將基巖海岸線設(shè)置為隔水邊界; 而基巖島嶼的海岸出現(xiàn)砂質(zhì)地質(zhì)體部分(如沙灘)時(shí),在數(shù)值模擬時(shí)應(yīng)將這類海岸設(shè)置為定水頭邊界,此時(shí),與沙質(zhì)島嶼關(guān)于海水入侵的邊界模擬相似。
由于人工開發(fā)島嶼資源,島嶼地下水可能受到以下污染:菜地、廁所、餐廳等人們的生活污水污染,采石場的礦渣、廢石等殘?jiān)廴荆l(fā)電廠機(jī)器運(yùn)轉(zhuǎn)的廢油及污水排放污染,垃圾廠的垃圾溶濾滲漏溶液的污染,污水處理廠的污水排放污染,以及海水入侵的影響[36]等等。目前,多數(shù)研究者仍然延用鄭春苗等[37-38]開發(fā)的MT3D、MT3DMS程序,如:Curtis等[39]于2006年利用MODFLOW模擬水流場,用MT3DMS程序模擬鈾(VI)的反應(yīng)運(yùn)移; 劉娟等[40]利用FEFLOW模擬地下水中四氯化碳污染羽的濃度場; 長安大學(xué)高小文等[41]利用MT3DMS溶質(zhì)運(yùn)移模塊模擬尾礦庫周邊環(huán)境污染,為環(huán)境治理提供參考依據(jù)。可見,目前污染物運(yùn)移模擬依舊利用三大熱門的地下水模擬軟件,主要還是依托于MT3DMS溶質(zhì)運(yùn)移模塊。
(1)基巖島嶼多數(shù)為大陸島,其地下水以淡水“蘑菇體”的形式存在,無需經(jīng)歷由咸變淡的過程,且無咸淡水過渡帶。
(2)基巖島嶼地下水的補(bǔ)給主要是降雨,其次是島嶼中、深層含水層與大陸的水力聯(lián)系補(bǔ)給。
(3)基巖島嶼的地質(zhì)模型可概化為3種:孔隙型、裂隙型和孔隙-裂隙型?;鹕綅u、松散巖類或風(fēng)化程度較高的基巖島嶼,可概化為孔隙型,其滲流符合非均質(zhì)各向異性或橫向各項(xiàng)同性的地下水運(yùn)動(dòng)方程,可直接用Visual Modflow、GMS和FEFLOW等軟件模擬; 島嶼地質(zhì)體以裂隙基巖為主,應(yīng)概化為裂隙型,其地下水運(yùn)動(dòng)規(guī)律符合立方定律,可采用MODFLOW-SURFACT模塊、TOUGH軟件、CONNETFLOW軟件進(jìn)行模擬; 對于孔隙-裂隙型模型,則需基于詳細(xì)的地質(zhì)體建模,結(jié)合TOUGH軟件或SURFACT程序模塊等軟件進(jìn)行耦合模擬。
(4)對于基巖島嶼的海水邊界設(shè)定,應(yīng)根據(jù)島嶼海岸線及物探與鉆探資料圈定出海水邊界。
(1)基巖島嶼地下水?dāng)?shù)值模擬需做到精細(xì)化,要將地質(zhì)體的各種參數(shù)刻畫細(xì)致,如每種巖層的給水度、滲透系數(shù)等。還應(yīng)細(xì)致刻畫,包括湖泊、河流等的地表水,對于降雨補(bǔ)給、蒸發(fā)排泄、潮汐影響以及氣候引起的水量變化也需仔細(xì)確定。此外,還應(yīng)注意與大陸的中、深層含水層的水力聯(lián)系的補(bǔ)給項(xiàng)。
(2)基巖島嶼多數(shù)為孔隙與裂隙并存巖體,其地下水賦存地質(zhì)模型需考慮雙重或多重介質(zhì)耦合模型。目前,TOUGH軟件、MODFLOW-SURFACT程序模塊或CONNETFLOW軟件可實(shí)現(xiàn)雙重或多重介質(zhì)的模擬,但該方面的文獻(xiàn)較少,仍然是一個(gè)研究難點(diǎn),希望廣大研究者可以研發(fā)出更完美的模擬軟件。
(3)對于基巖島嶼地下水模擬的驗(yàn)證也是重中之重,數(shù)值模擬結(jié)果應(yīng)與島嶼水量均衡,且與實(shí)測水量、水位相一致,模擬是一個(gè)不斷修正與完善的過程,不能一蹴而就。