肖 斌，許 模，曾 科，王 梅

(1.成都理工大學地質(zhì)災害防治與地質(zhì)環(huán)境保護國家重點實驗室，四川成都610059;2.四川省地質(zhì)環(huán)境監(jiān)測總站，四川成都610084)

0 前言

巖溶介質(zhì)在空間上的多尺度性決定了巖溶介質(zhì)中地下水流態(tài)的多重性，巖溶介質(zhì)的模擬概化是水文地質(zhì)模擬上的一大難題。目前，對于多孔介質(zhì)的模擬已經(jīng)較為成熟，因此，在模擬巖溶介質(zhì)時，國內(nèi)外學者多采用等效連續(xù)多孔介質(zhì)模型[1]。如Georgios Panagopoulps采用等效多孔介質(zhì)模型模擬了希臘Trifilia地區(qū)巖溶含水層[2];毛邦燕等人運用等效多孔介質(zhì)模型對不同型的巖溶介質(zhì)進行分區(qū)概化[3]，都取得了較好的效果。但這種方法只適用于模擬范圍大，巖溶發(fā)育程度較弱的含水層。對于含水層中存在對地下水具有強烈控制作用的大尺度(1～104 cm)巖溶管道[4]時，這種方法并不能對其進行有效的刻畫。那么，能否在現(xiàn)有成熟的等效連續(xù)多孔介質(zhì)模型的理論基礎之上，找到一種恰當?shù)膶r溶介質(zhì)管道概化模擬的方法呢?筆者在本文中嘗試利用三維有限差分地下水流動模型(Modflow)，尋求一種對于巖溶管道的合理概化，從而使得對強巖溶發(fā)育區(qū)地下水流的模擬更加合理、精確。

Modflow是由美國地質(zhì)調(diào)查局的McDonald和Harbaugh開發(fā)出來的，基于連續(xù)多孔介質(zhì)理論的地下水流模擬軟件，其最顯著的一個特點就是軟件的模塊化結構，即包括一個主程序以及若干相對獨立的子模塊。模塊化的設計使得模型中各參數(shù)與邊界條件的輸入和修改相對獨立，計算過程又能使模型各個單獨的模塊緊密的聯(lián)系起來[5]。這種模塊化的設計為我們尋求巖溶管道的概化方法提供了可能性。

1 Modflow中現(xiàn)有可用模塊

目前，在Modflow的各個模塊中，可以考慮用于概化模擬巖溶管道的有:單元滲流計算子模塊(BCF)、排水溝子模塊(DRN)、河流子模塊(RIV)和溪流子模塊(STR)。

1)單元滲流計算子模塊(BCF)

在Modflow中，單元滲流計算子模塊用于計算相鄰單元之間的水力傳導系數(shù)以及計算單元之間的地下水滲流量，也可以用于計算含水層由于貯水量的變化所吸收或釋放的水量[5]。在實際應用中，是通過滲透系數(shù)(K)和貯水系數(shù)(S)的設置來實現(xiàn)的。在模擬強巖溶介質(zhì)含水層時，可溶巖巖體中巖溶管道發(fā)育規(guī)模較小的部分，可根據(jù)實際鉆孔抽水試驗獲得的數(shù)據(jù)進行賦值;而根據(jù)陳崇希的管道K值表式(1)，可將其中的大型巖溶管道部分賦值為由該式折算得出的滲透系數(shù)，貯水系數(shù)按100%孔隙率計算。

式中:K為巖溶管道滲透系數(shù)(LT－1);d為巖溶管道內(nèi)徑(L);γ為流體的重率(FL－3);μ為流體的動力粘度(FTL－2);n為空隙率(無量綱)[6，7，8]。

2)排水溝子模塊(DRN)

在Modflow中，排水溝子模塊原設計目的是模擬農(nóng)用排水溝渠的排水效果。排水溝從含水層中排泄地下水，排水量正比于含水層的水頭與排水溝的高程之差，其中要求含水層水頭高于排水溝的高程，而當含水層水頭低于該排水溝高程時，就不起排水作用[5]。圖1表示某計算單元的橫截面圖，假定排水溝中僅部分有水，因而排水溝的水頭大約等于其半高處的高程di，j，k。模型所計算的計算單元水頭(hi，j，k)是計算單元的的水頭平均值。排水溝的水頭di，j，k只局限于溝渠之內(nèi)而不是將其延伸至整個計算單元。在排水溝與含水層水頭hi，j，k之間的區(qū)域，垂向面上存在輻射流或半輻射流，其特點一般為，越接近溝渠，水頭梯度越大。在地下水流向排水溝過程中的水頭損失是形成水頭差hi，j，k－ di，j，k的原因。

圖1 排水溝計算單元(i，j，k)的匯聚水流之水頭損失橫剖面示意圖

根據(jù)排水溝子模塊的上述特性，可以將排水溝子模塊應用于巖溶含水介質(zhì)向巖溶管道排水的計算[9]，用以反映管道和含水介質(zhì)間的水力交換。模型中，排水量正比于含水層水頭與暗河管道高程之差，當含水層水頭低于排水溝高程時，無排水效果[5]:

式中:QDi，j，k為排水溝排水量，CDi，j，k為等效水力傳導系數(shù)，用于描述排水溝和計算單元(i，j，k)間的總水頭損失。在巖溶介質(zhì)的暗河管道附近巖體導水能力強，可以采用較大的等效水力傳導系數(shù)反映巖溶管道快速排泄地下水的特點。

3)河流子模塊(RIV)

在Modflow中，河流與溪流子模塊是用來模擬地表水體對地下水流的影響，包括河流、溪流向地下水系統(tǒng)提供水，以及地下水系統(tǒng)向河流、溪流排泄地下水，具體是補給還是排泄取決于河流單元與相鄰單元地下水之間的水力梯度。河流子模塊在模型中模擬其與含水層計算單元(i，j，k)之間的水力聯(lián)系時由下面一套方程來實現(xiàn)[5]:

式中:QR為河流與含水層之間的流量，水流由河流流向含水層時取正值;HR為河流的水位;RBOT為河床基底處高程;hi，j，k為計算單元的水頭值;CR為河流－含水層互相聯(lián)接的水力傳導系數(shù)(由河流長度、寬度、河床底積物厚度、河床底積物滲透系數(shù)計算得出)。由式3－a、3－b可以看出，該函數(shù)和排水溝與含水層間流量的函數(shù)相似，不同之處在于，河流不但能反映出周邊地下水向管道中匯流的特征，亦能反映出管道內(nèi)水流向周邊含水層的補給。

4)溪流子模塊(STR)

在Modflow中，溪流子模塊的作用是計算溪流中的總流量和模擬地表水體和地下水的相互作用。溪流子模塊計算地表水和地下水交換過程的函數(shù)與河流子模塊并無二處。不同的是，在設置溪流時，溪流被分成節(jié)和段，每一段對應于單個的單元網(wǎng)格，而節(jié)是由以順流順序連接在一起的一組網(wǎng)格單元組成，通過指定每一節(jié)中第一段的溪流流入量來計算溪流流量，計算每段中鄰近下游溪流段的流量等于上游溪流段流入量加上(減去)上游溪流段含水層流入(流出)的滲流量。相比于河流子模塊，溪流子模塊在輸入時要求輸入溪流流量，這就對巖溶管道模擬時添加了一個控制條件，模擬時不但減少了模型的校正過程，同時也提高了模擬精度。

2 模塊實際應用效果檢驗

以一簡易理想模型對各模塊實際應用效果進行驗證說明。該驗證模型范圍為一長3 000m，寬2 000m的區(qū)域，在X=800處為地表分水嶺，分水嶺右側(cè)Y=果行驗證處為一由西向東發(fā)育的巖溶管道，模型右側(cè)為一高程為200的河流排泄基準面，同時在模型X=2000，Y=1000處設置一地下水位觀測點，以檢驗暗河水位的動態(tài)變化。模型運行時間為一年，假定其氣候類型為北半球亞熱帶季風氣候，即夏季多雨冬季少雨。

圖2 單元滲流計算子模塊模擬水位等值線

使用單元計算子模塊概化模擬巖溶管道的模擬結果(見圖2)可以看出，通過為巖溶管道區(qū)域賦上折算的水文地質(zhì)參數(shù)后，模型反應出了巖溶管道周圍地下水的運動規(guī)律:在管道發(fā)育區(qū)域，管道周圍地下水水頭較水流通暢的管道內(nèi)的水頭高，表現(xiàn)出管道圍巖的地下水向管道中匯集。一般說來，大型的巖溶管道往往是巖溶地下水的控水介質(zhì)，它控制著巖溶水的運動，并使之成為具獨特形狀的地下水位曲面[4]，但這一特性在模型中并未表現(xiàn)。同時，從監(jiān)測鉆孔的水位動態(tài)變化(見圖6)也可以看出，模型中暗河水位也沒有呈現(xiàn)出現(xiàn)實中隨干枯兩季大幅變化的性質(zhì)，且模型中暗河水位的變化與大氣降雨也呈現(xiàn)出明顯的滯后性。這是因為使用折算含水參數(shù)概化巖溶管道時，并沒有改變其為孔隙含水介質(zhì)的本質(zhì)，這種方法可以簡單的模擬巖溶管道的儲水和輸水功能，但不能有效的模擬巖溶管道的控水功能。同時，此方法的前提條件是溶蝕管道處于飽水帶中，對于補給源自地表水的伏流型管道并不能進行模擬。因此，這種運用折算參數(shù)模擬巖溶管道的方法只適用于埋深大、環(huán)境相對封閉、處于飽水狀態(tài)的巖溶管道，此時管道主要起到儲水和輸水的作用，而不適用于模擬埋深淺、動態(tài)變化大、處于包氣帶或季節(jié)交替帶的巖溶管道。

圖3 排水溝子模塊模擬水位等值線

圖4 河流子模塊模擬水位等值線

如圖3所示，使用排水溝子模塊模擬巖溶管道可以很好的反應出巖溶管道介質(zhì)的控水作用:管道兩側(cè)地下水向管道匯集，地下水位形成一定的漏斗形態(tài)。同時，監(jiān)測鉆孔顯示的模型中暗河水位的動態(tài)特征也符合自然狀態(tài)中暗河水量隨干枯兩季的起伏變化。但是，如果巖溶管道發(fā)育于包氣帶，管道中的地下水補給源自地表水時，模擬中排水溝水頭高于含水層水頭時，排水溝并不會像天然狀態(tài)下那樣表現(xiàn)出暗河管道內(nèi)的水對地下水的補給。因此，通過適當?shù)馁x值將排水溝子模塊用于模擬常年處于周邊地下水位以下的、開放式的巖溶管道是合適的，但并不適用于處于地下水位以上的伏流型巖溶管道的模擬。

如圖4所示，河流子模塊在模擬巖溶管道時體現(xiàn)出了巖溶管道在開放系統(tǒng)中的控水作用，同時也能模擬出管道內(nèi)地下水和含水介質(zhì)中地下水的相互交換。監(jiān)測鉆孔也顯示(圖6)，模型中鉆孔水位隨干枯兩季劇烈變化。因此，河流子模塊能很好的用于模擬各類大型巖溶管道。由公式(3)知，在缺乏全面的地質(zhì)資料時，使用河流－含水層互相連接的水力傳導系數(shù)(CR)來表示一個三維流過程是一種經(jīng)驗做法，在模型調(diào)試時需要經(jīng)過進行大量的校正。

圖5 溪流子模塊模擬水位等值線

圖6 不同模塊模擬下監(jiān)測鉆孔水位動態(tài)變化特征

使用溪流子模塊概化巖溶管道的模擬結果(圖5)同河流子模塊模擬的結果十分相似。溪流子模塊不但模擬出了巖溶管道應有的特性，同時，有了流入端流量數(shù)據(jù)的控制，溪流在模擬有伏流入口的巖溶管道時更加精準。

表1 各子模塊模擬巖溶管道適宜性對比

3 結論

管道－非管道介質(zhì)的地下水流動模擬問題是地下水流動模擬領域難度最大的問題。筆者在現(xiàn)今已十分成熟的等效連續(xù)多孔介質(zhì)模型的基礎之上，將復雜抽象的管道條件用相對具體且成熟的模型模塊進行概化，同時保證其在水力聯(lián)系方面的統(tǒng)一性。通過模型實踐證明，在Modflow中使用河流子模塊和溪流子模塊概化能最好地模擬巖溶管道周邊地下水的流動和動態(tài)變化;此二者在實踐中需要較多參數(shù)的支持以及較多的模型矯正。這種概化巖溶管道的方法在保證巖溶區(qū)地下水滲流場較為精確模擬的前提下，大大簡化了模型的模擬概化過程，為實際應用提供一種思路與可能。由于巖溶含水介質(zhì)自身的特殊性和復雜性，該方法在實際應用過程中仍待改進。

[1]劉國，毛邦燕，許模等.復雜巖溶含水介質(zhì)概化初探[J].物探化探計算技術，2007，29(5):439 －442.

[2]Georgios Panagopoulos.Application of MODFLOW for simulating groundwater flow in the Trifilia karst aquifer，Greece[J].Environ Earth Sci，2012，67(7):1877 －1889.

[3]毛邦燕.復雜巖溶介質(zhì)礦井涌水量的三維數(shù)值模擬研究[D].四川成都:成都理工大學，2005.

[4]陳雨孫，邊際.巖溶水的介質(zhì)和運動[J].中國巖溶，1988，7(3):229 －234.

[5]Michael G.McDonald，Arlen W.Harbaugh.A modular three－dimensional finite－difference ground－water flow model:U.S.Geological Survey Techniques of Water－ Resources Investigations[M].Washington:UNITED STATES GOVERNMENT PRINTING OFFICE，1988.

[6]陳崇希.巖溶管道－裂隙－孔隙三重空隙介質(zhì)地下水模型及模擬方法研究[J].地球科學，1995，20(4):361 －366.

[7]陳崇希，胡立堂.滲流－管流耦合模型及其應用綜述[J].水文地質(zhì)工程地質(zhì)，2008，35(3):70 －75.

[8]陳建梅，陳崇希.廣西北山巖溶管道－裂隙－孔隙地下水數(shù)值模擬初探[J].水文地質(zhì)與工程地質(zhì)，1998，25(4):50－54.

[9]John J.Quinn，David Tomasko，James A.Kuiper.Modeling complex flow in a karst aquifer[J].Sedimentary Geology.2006，184(3):343 －353.