典型石漠化地區(qū)巖溶水系統(tǒng)循環(huán)演化分析*

姚邦杰 劉 琦② 任 標(biāo) 王 涵

( ①同濟(jì)大學(xué)土木工程學(xué)院地下建筑與工程系 上海 200092)

( ②巖土及地下工程教育部重點(diǎn)實(shí)驗(yàn)室( 同濟(jì)大學(xué)) 上海 200092)

0 引 言

巖溶石漠化是指在巖溶環(huán)境下，受到人類不合理的社會(huì)活動(dòng)和經(jīng)濟(jì)活動(dòng)等的干擾和破壞造成的土壤嚴(yán)重侵蝕、土壤肥力和生產(chǎn)力大幅下降，并導(dǎo)致基巖大面積裸露且出現(xiàn)類似于荒漠化景觀的生態(tài)環(huán)境問題，它是土地荒漠化的主要類型之一(王世杰等，2003)。

巖溶地區(qū)因?yàn)樗畬?duì)可溶性巖石的侵蝕以及兩者之間的能量物質(zhì)交換形成了該地區(qū)特有的地形地貌形態(tài)，這種地形地貌形態(tài)可以分為地下和地表兩個(gè)子系統(tǒng)。其中管道、溶洞、地下暗河等組成了地下子系統(tǒng)，而峰叢、溶溝、溶槽、溶隙、巖溶漏斗、天窗、落水洞等是地表子系統(tǒng)的組成部分( Ford，1989; 盧耀如，1999) 。巖溶水徑流在地表發(fā)生并分配，然后在地下子系統(tǒng)中輸送和調(diào)蓄。巖溶系統(tǒng)存在孔隙、裂隙和管道3 種介質(zhì)空間，孔隙流、裂隙流和管道流并存。因此，巖溶地下水在結(jié)構(gòu)上表現(xiàn)為功能上的耗散結(jié)構(gòu)、三維的空間地域結(jié)構(gòu)、二元流場(chǎng)形態(tài)結(jié)構(gòu)還有不均一的雙重含水介質(zhì)結(jié)構(gòu)( White，1969) 。巖溶地下水一般不符合達(dá)西定律，表現(xiàn)為紊流或者混合流( Ford，1989; 石朋等，2012) ，正是由于巖溶水系統(tǒng)循環(huán)演化的特殊性，水不僅會(huì)攜帶土壤發(fā)生地表流失，而且會(huì)發(fā)生地下漏失，導(dǎo)致巖溶石漠化現(xiàn)象突出。巖溶水系統(tǒng)循環(huán)演化受到自然和人為因素控制。自然因素包括降雨量、蒸發(fā)量、年徑流量、地形地貌、地質(zhì)構(gòu)造、巖性、水文地質(zhì)條件等; 人為因素包括土地利用形式、水利水電工程、向河流排放污染水源等( Stringfield et al.，1979; Rosenberry et al.，2008) 。國內(nèi)外學(xué)者進(jìn)行了大量的研究，其研究方法主要包括:室內(nèi)實(shí)驗(yàn)法、野外調(diào)查試驗(yàn)法、理論分析法和數(shù)值模擬法( 史運(yùn)良等，1992; 盧衛(wèi)中，1995; 袁道先等，1996; 王臘春等，2000; 韓寶平等，2004; 章程等，2004; 劉延惠等，2005; 曾科，2012) 。但目前少有對(duì)典型巖溶石漠化地區(qū)的巖溶水系統(tǒng)演化問題進(jìn)行系統(tǒng)性研究，而這一循環(huán)過程恰恰是導(dǎo)致石漠化地區(qū)土壤地表流失和地下漏失方向和趨勢(shì)的前提條件和關(guān)鍵性問題。

本文以典型巖溶石漠化地區(qū)貴州貞豐—關(guān)嶺花江高原峽谷區(qū)為研究區(qū)，對(duì)該區(qū)的大氣降水、河水、泉水進(jìn)行水化學(xué)、氫氧同位素分析，研究巖溶水的補(bǔ)、徑、排特征，水化學(xué)特征的變化; 利用Modflow對(duì)研究區(qū)進(jìn)行地下水?dāng)?shù)值模擬，用PEST 法( 參數(shù)評(píng)估) 進(jìn)行模型識(shí)別與驗(yàn)證，通過構(gòu)建研究區(qū)的地下水流場(chǎng)，計(jì)算研究區(qū)的降雨入滲量、地下暗河的排泄量、地下水泄流量、水均衡情況等，以定量研究區(qū)內(nèi)的巖溶水系統(tǒng)循環(huán)演化過程，從而為進(jìn)一步研究石漠化地區(qū)土壤地表流失和地下漏失的方向和趨勢(shì)提供相關(guān)理論依據(jù)。

1 研究區(qū)概況

研究區(qū)位于貴州省安順市關(guān)嶺縣的北盤江北岸與打邦河西岸所夾的三角地帶，東經(jīng)105°18'45″～105°52'30″、北緯25°33'45″～25°50'00″，研究區(qū)內(nèi)巖溶地貌發(fā)育，地形高差極大，地形復(fù)雜。研究區(qū)年均降雨量約為1100 mm，年均溫約18 ℃，海拔376 ～1828 m，研究區(qū)河流兩岸高差由于河流的下切作用非常大，是典型的巖溶高原峽谷地貌。該區(qū)石漠化發(fā)育程度強(qiáng)烈，石漠化面積約占研究區(qū)總面積的49%，且主要是中度和強(qiáng)度石漠化，已經(jīng)對(duì)人類的生存發(fā)展造成威脅。

該區(qū)主要出露三疊系、白堊系及第四系地層，其中三疊系地層出露面積占95%以上，碳酸鹽巖出露面積約占85%，根據(jù)巖性、巖層的水理性質(zhì)及地下水賦存狀態(tài)，將研究區(qū)地下水劃分為3 類，包括碳酸鹽巖裂隙溶洞水、碎屑巖碳酸鹽巖溶洞裂隙水和層狀巖類裂隙水。研究區(qū)內(nèi)共有主要地下河兩條，巖溶管道總長度約39 km，巖溶大泉20 個(gè)。研究區(qū)主要發(fā)育有4 組裂隙，并且主要發(fā)育1 組北西向的裂隙，巖溶發(fā)育的主控方向即主要受北西向裂隙控制，幾條主要地下暗河走向均為北西向，主要分布于研究區(qū)打邦河西岸，各含水類型及地下河的分布見圖1。

2 巖溶水化學(xué)特征及補(bǔ)、徑、排特征

2018 年9 月，在研究區(qū)進(jìn)行野外調(diào)查并采集雨水樣4 組，在采集雨水樣2 d 后采集泉水樣20 組，河水樣2 組( 采樣點(diǎn)見圖1) ，基本覆蓋研究區(qū)泉水出露點(diǎn)，采樣后進(jìn)行氫氧同位素和水化學(xué)成分測(cè)定，結(jié)果見表1。

圖1 研究區(qū)巖溶水系統(tǒng)水文地質(zhì)圖Fig. 1 Hydrogeological map of karst groundwater system in the study area

因隨著海拔的增加，氣溫相應(yīng)降低，此時(shí)由于蒸發(fā)分餾作用的影響，大氣降水中氫氧同位素組成會(huì)發(fā)生線型變化，因此研究區(qū)大氣降水中氫氧同位素的組成也具有顯著的高程效應(yīng)( 陳中笑等，2010) ，據(jù)此可推算出地下水的補(bǔ)給高程，它是判斷地下水的補(bǔ)給位置以及補(bǔ)給范圍的基礎(chǔ)，巖溶地區(qū)地下水循環(huán)更替較快，該方法更加有效( 黃荷等，2016) 。通過研究區(qū)水樣分析計(jì)算結(jié)果表明，花江雨水中δD值與海拔成正相關(guān)關(guān)系: H=－28.976 δD－998.25( 圖2) 。每個(gè)主要巖溶泉的氫氧同位素值代入研究區(qū)雨水δD 值與降雨高程的關(guān)系公式，計(jì)算出主要巖溶泉的大致補(bǔ)給高程，并根據(jù)結(jié)果標(biāo)注到地圖上( 圖3) 。大部分泉水的補(bǔ)給范圍與巖溶發(fā)育主控方向( 北西—南東) 一致，泉水的補(bǔ)給區(qū)域主要受斷層、地下管道控制，其中16 號(hào)泉水( 地下暗河) 補(bǔ)給范圍最大，從西北邊的上關(guān)鎮(zhèn)到東北部的壩坎周圍均為其匯水范圍。根據(jù)泉水補(bǔ)給范圍推測(cè)出研究區(qū)地下水流向主要為北西—南東向，靠近南部北盤江時(shí)地下水流向轉(zhuǎn)為近北—南向，局部地區(qū)如花江鎮(zhèn)東南部地下水流向因地形高差原因會(huì)有一定變化。

3 研究區(qū)地下水?dāng)?shù)值模擬

為了深入了解研究區(qū)巖溶水資源轉(zhuǎn)化量的關(guān)系，本文使用MODFLOW 軟件進(jìn)行地下水?dāng)?shù)值模擬，根據(jù)室內(nèi)外資料構(gòu)建地下水流場(chǎng)，計(jì)算研究區(qū)的水均衡情況、降雨入滲量、地下河管道的排泄量、地下水泄流量等，以定量研究區(qū)內(nèi)巖溶水轉(zhuǎn)化關(guān)系。

表1 研究區(qū)泉水、地表水、雨水同位素及水化學(xué)分析數(shù)據(jù)Table 1 Spring，surface water，rainwater isotope and water chemical analysis data of the study area

圖2 研究區(qū)雨水δD 值的高程效應(yīng)Fig. 2 Elevation effect of rainfall δD value in the study area

3.1 地下水?dāng)?shù)值模擬

按照研究區(qū)的水文地質(zhì)條件圈定研究區(qū)的邊界，并將研究區(qū)的所有匯水面積圈定為模型范圍，最終模型橫向長31.7 km，縱向長38.3 km。其中大面積發(fā)育的碳酸鹽巖為模型的主要含水介質(zhì)。研究區(qū)的邊界主要分為河流和斷層兩類，北部由兩條大型壓性斷層構(gòu)成邊界，南部由北盤江和打邦河兩條大型河流構(gòu)成研究區(qū)邊界。模型的邊界條件主要分為4 類，第1 類為大型河流構(gòu)成的河流邊界，分為兩段，模型南面為北盤江，東面為北盤江支流打邦河，大部分地下水補(bǔ)給入這兩條河流，因此將這兩條河流概化為模型的河流邊界和流出邊界; 第2 類為地下暗河構(gòu)成的排水溝邊界，模型邊界內(nèi)主要發(fā)育有許凹地下河，運(yùn)用MODFLOW 中的“排水溝”子程序來模擬該條地下河，為流出邊界; 第3 類為模型西南邊海拔最高處的定流量邊界; 第4 類為補(bǔ)給邊界。在穩(wěn)定流計(jì)算時(shí)該地區(qū)的降雨補(bǔ)給量設(shè)定為400 mm·a－1。非穩(wěn)定流計(jì)算時(shí)，降雨補(bǔ)給量按氣象站實(shí)測(cè)值的0.33 倍給出。模擬范圍、研究區(qū)邊界及模型邊界條件如圖4 所示。

圖3 研究區(qū)主要巖溶泉推測(cè)補(bǔ)給范圍Fig. 3 Speculation range of major karst springs in the study area

表2 水文地質(zhì)參數(shù)初值Table 2 Initial values of hydrogeological parameters

依據(jù)上述水文地質(zhì)概念模型，可將研究區(qū)的裂隙巖溶地下水流問題概化為非均質(zhì)各向異性介質(zhì)中的三維非穩(wěn)定流問題，研究區(qū)在垂向上分為6 層，含水層均為潛水含水層，水文地質(zhì)參數(shù)初值如表2 所示。對(duì)模擬區(qū)進(jìn)行網(wǎng)格剖分。模型在垂向上除了原有的分區(qū)以外不進(jìn)行更細(xì)致的劃分。模型的網(wǎng)格剖分總體上采用154×159 m2的形式。剖分后活動(dòng)單元格數(shù)約為24 500 個(gè)，模擬面積1214.1 km2。以研究區(qū)內(nèi)主要巖溶泉的高程值作為該點(diǎn)地下水位的擬合值，采用MODFLOW 中的PEST 算法進(jìn)行參數(shù)反演。

3.2 數(shù)值模擬結(jié)果

3.2.1 地下水位變化

研究區(qū)地下水流場(chǎng)模擬結(jié)果如圖5 所示，地下水流向大致呈北西—南東向，分別匯入南部北盤江和東部打邦河，地下水流向與該地區(qū)巖溶發(fā)育主控方向一致，較好反映了研究區(qū)的地下水流場(chǎng)。OW-4觀測(cè)點(diǎn)的水位變化情況( 圖6) ，該區(qū)地下水位年變化幅度較大，季節(jié)性強(qiáng)，3 ～5 月為地下水枯水期，8～10 月為地下水豐水期，地下水位變化情況相對(duì)降雨量變化情況有一定的滯后性。

3.2.2 地下暗河流量

圖4 模型范圍、研究區(qū)邊界及模型邊界條件Fig. 4 Model range，study area boundaries，and model boundary conditions

模擬得到地下暗河流量變化情況( 圖7) ，地下暗河流量變化與降雨量變化情況存在一定的滯后性，地下暗河流量隨雨季來臨而上升速度加快，模擬地下暗河流量的回落速度在降雨量減小時(shí)較慢。結(jié)合模擬地下河出口的總排泄流量與實(shí)際觀測(cè)流量的擬合效果分析，模擬雨季最大流量小于實(shí)測(cè)雨季最大流量，而模擬枯季最小流量大于實(shí)測(cè)枯季最小流量，可見模擬流量變化情況比實(shí)際情況稍緩，但誤差不大，模擬過程和實(shí)際條件下地下水系統(tǒng)排泄量的變化特征基本吻合。

3.2.3 水均衡

模擬計(jì)算出的研究區(qū)年均地下水均衡情況( 表3) ，本次模擬期中降雨補(bǔ)給量是主要的補(bǔ)給來源; 側(cè)向補(bǔ)給量是次要補(bǔ)給來源。河流側(cè)向排泄量占排泄量的60.21%，雖然占比較大，但唯一一條地下河的管道排泄量占總排泄量的24.92%，符合巖溶地區(qū)的實(shí)際情況。

將2018 年每個(gè)月研究區(qū)的水均衡情況歸納如圖8 所示。在雨季，大量增加的地下水無法迅速通過地下暗河或河流排泄，其中的大部分被地下含水層儲(chǔ)存，而地下水向地下暗河及河流的排泄量迅速增大。在旱季，因該地區(qū)幾乎沒有大氣降水的補(bǔ)給，地下含水層中儲(chǔ)存的地下水開始發(fā)揮補(bǔ)給功能，而地下水向地下暗河及河流的排泄量迅速減小。一年之中，側(cè)向排泄量和側(cè)向補(bǔ)給量變化較小，地下暗河的排泄量變化幅度顯著大于河流排泄量的變化幅度。

根據(jù)模擬計(jì)算結(jié)果看出，研究區(qū)33%的大氣降水通過表層巖溶裂隙帶的下滲作用轉(zhuǎn)化為地下水，其中24.9%的地下水通過地下暗河排泄、60.2%的地下水通過河流排泄作用進(jìn)入北盤江和打邦河成為地表水，其余的14.9%的地下水側(cè)向排泄成為區(qū)域外的地下水。

4 結(jié) 論

通過野外調(diào)查、水樣分析及數(shù)值模擬等方法，對(duì)貴州貞豐—關(guān)嶺花江高原峽谷石漠化治理示范區(qū)巖溶水系統(tǒng)的補(bǔ)、徑、排特征，水化學(xué)特征及水資源量等進(jìn)行分析，得到如下結(jié)論:

圖5 OW-4 觀測(cè)點(diǎn)地下水位模擬值Fig. 5 Groundwater level simulation value of OW-4 observation point

圖6 模擬地下水流場(chǎng)和地下水位擬合點(diǎn)Fig. 6 Simulated groundwater flow field and groundwater level fitting point

表3 年均地下水均衡表Table 3 Annual average groundwater balance table

(1) 研究區(qū)內(nèi)大氣降水的氫氧同位素組成具有顯著的高程效應(yīng)，據(jù)此計(jì)算出了研究區(qū)主要巖溶泉的補(bǔ)給高程。綜合水化學(xué)分析推斷出的補(bǔ)給區(qū)地層，以及氫氧同位素分析計(jì)算出的補(bǔ)給高程，得到了研究區(qū)主要巖溶泉的補(bǔ)給范圍。得出研究區(qū)河水上下游的水化學(xué)、氫氧同位素變化與河流兩岸的海拔、巖性的變化高度一致，可見地下水與河水之間聯(lián)系密切。

圖7 地下暗河流量變化模擬情況Fig. 7 Simulation of underground dark river flow changes

圖8 研究區(qū)2018 年水均衡變化情況Fig. 8 Water balance changes in the study area in 2018

( 2) 建立了關(guān)嶺—貞豐花江地區(qū)MODFLOW 數(shù)值模型，根據(jù)實(shí)際情況將含水巖層分為6 層，模擬過程和實(shí)際條件下地下水系統(tǒng)排泄量的變化特征基本吻合，根據(jù)模擬結(jié)果得到，研究區(qū)33%的大氣降水通過表層巖溶裂隙帶的下滲作用轉(zhuǎn)化為地下水，其中24.9%的地下水通過地下暗河排泄、60.2%的地下水通過河流排泄作用進(jìn)入北盤江和打邦河成為地表水，其余的14.9%的地下水側(cè)向排泄成為區(qū)域外的地下水。其結(jié)果符合巖溶地區(qū)的實(shí)際情況，但由于巖溶地區(qū)地下水流情況十分復(fù)雜，且地層資料精度較差、地下水位觀測(cè)數(shù)據(jù)缺乏，地下水的模擬仍有一定誤差。

分析結(jié)果基本反映了實(shí)際巖溶水系統(tǒng)排泄量變化及地下水流向特征，初步闡明了該區(qū)域巖溶水系統(tǒng)循環(huán)演化特征，為巖溶石漠化地區(qū)的土壤地表流失和地下漏失的方向和趨勢(shì)研究及防治提供理論依據(jù)。