Visual Modflow在地下水污染物運(yùn)移模擬中的應(yīng)用

尉鵬翔

(山西省水文水資源勘測(cè)局太原分局,山西太原 030002)

近年來(lái),人們普遍注意到由于超量開(kāi)采地下水和水質(zhì)污染,地下水生態(tài)環(huán)境受到破壞,且日趨嚴(yán)重。應(yīng)用科學(xué)有效的計(jì)算方法尋找處理地下水污染問(wèn)題的方法,對(duì)于地下水的污染防治有著極其重要的作用[1]。目前地下水系統(tǒng)數(shù)值模擬方法主要有有限差分法(FDM)、有限單元法(FEM)、邊界元法(BEM)和有限分析法(FAM)等,可以用來(lái)研究海水入侵、地面沉降、地下水位降低、地下水過(guò)量開(kāi)采和地下水污染等問(wèn)題。隨著計(jì)算機(jī)技術(shù)的快速發(fā)展,數(shù)值計(jì)算方法在地下水資源分析評(píng)價(jià)中得到廣泛應(yīng)用。筆者擬借助Visual Modflow商業(yè)化軟件研究北京某污染場(chǎng)區(qū)地下水溶質(zhì)運(yùn)移問(wèn)題。

1 研究區(qū)概況

所研究污染場(chǎng)區(qū)位于北京市區(qū)東南部,地貌上屬于永定河沖洪積扇的中下部,場(chǎng)區(qū)地形基本平坦,地面以下30m深度范圍內(nèi)一般存在2個(gè)含水層,局部區(qū)域有上層滯水,單層含水層厚度一般小于10m,僅部分地區(qū)大于10m。本文只研究埋深在10～20m的含水層和埋深在25m左右的承壓含水層以及它們之間的黏土弱透水層,鑒于此,模型在空間上概化為3層:潛水含水層、弱透水層和承壓含水層。

場(chǎng)區(qū)目前潛水流向由北向南,北部邊界是側(cè)向補(bǔ)給邊界,主要接受側(cè)向補(bǔ)給,南部邊界為整個(gè)含水層系統(tǒng)的微量流出邊界,由于地下水整體流向由北向南,故模擬區(qū)東、西兩側(cè)可以流線為邊界,將它們?cè)O(shè)為隔水邊界。承壓含水層水流向由西北向東南,西北部為側(cè)向補(bǔ)給邊界,主要接受側(cè)向補(bǔ)給,東南部為側(cè)向流出邊界,含水層底部大多為透水性較差的黏土地層,對(duì)系統(tǒng)內(nèi)部影響很小,定義為隔水邊界。

場(chǎng)區(qū)含水層分布廣、厚度大,在常溫常壓下地下水運(yùn)動(dòng)符合達(dá)西定律,考慮淺、深層之間的流量交換以及軟件的特點(diǎn),地下水運(yùn)動(dòng)可概化成空間三維流;地下水系統(tǒng)的垂向運(yùn)動(dòng)主要是層間的越流,三維立體結(jié)構(gòu)模型可以很好地解決越流問(wèn)題;地下水系統(tǒng)的輸入和輸出隨時(shí)間、空間變化,故地下水為非穩(wěn)定流;參數(shù)隨空間的變化體現(xiàn)了系統(tǒng)的非均質(zhì)性,但沒(méi)有明顯的方向性,所以參數(shù)概化成各向同性。

2 地下水污染數(shù)值模擬

2.1 數(shù)學(xué)模型的建立

計(jì)算采用Visual Modflow軟件進(jìn)行模擬,在進(jìn)行地下水模擬前,首先必須對(duì)模擬區(qū)進(jìn)行矩形剖分。數(shù)值模型分潛水含水層、弱透水層和承壓含水層共剖分50×50×3個(gè)單元。此外,研究重點(diǎn)是第四系松散孔隙潛水和第一承壓含水層含水系統(tǒng)。系統(tǒng)的水均衡要素的補(bǔ)給項(xiàng)主要考慮大氣降水入滲、邊界側(cè)向補(bǔ)給等,其他形式的補(bǔ)給作用相對(duì)較弱;排泄項(xiàng)主要考慮蒸發(fā)和邊界側(cè)向流出。

在不考慮水的密度變化條件下,孔隙介質(zhì)中地下水在三維空間的流動(dòng)可用偏微分方程來(lái)表示[2-3]:

式中:Kxx、Kyy、Kzz分別為滲透系數(shù)在x、y、z方向上的分量;h為水頭;W為在非均衡狀態(tài)下通過(guò)均質(zhì)、各向同性土壤介質(zhì)單位體積通量,即地下水的源和匯;Ss為孔隙介質(zhì)儲(chǔ)水率;t為時(shí)間。

2.2 定解條件的處理

污染物在地下水中運(yùn)移轉(zhuǎn)化的過(guò)程是極其復(fù)雜的,不僅有對(duì)流、彌散作用,而且可能有離子交換、化學(xué)反應(yīng)、微生物降解等化學(xué)和生物化學(xué)作用,在不同的環(huán)境下這些作用的方向和強(qiáng)度又有很大不同,本文只考慮對(duì)流和彌散作用。

地下水溶質(zhì)運(yùn)移滿足下面的對(duì)流彌散方程:

式中:c為地下水中污染物濃度;xi為沿坐標(biāo)軸各方向的距離;Dij為水力擴(kuò)散系數(shù);Vi為地下水滲流速度;qs為源和匯的單位流量;cs為源和匯的濃度;Q為含水層孔隙率;∑Rk為化學(xué)反應(yīng)項(xiàng)。

2.3 污染物質(zhì)的確定

現(xiàn)場(chǎng)調(diào)查和地下水水質(zhì)分析表明,研究區(qū)的污染情況相當(dāng)復(fù)雜,并非僅有一個(gè)明顯的污染源,而是可能存在一個(gè)主要的污染源和多個(gè)次要的污染源,可能的污染源主要來(lái)自管道滲漏和不符合標(biāo)準(zhǔn)的工業(yè)廢水直接排放。參照美國(guó)EPA標(biāo)準(zhǔn)確定的嚴(yán)重污染物,研究區(qū)水質(zhì)分析結(jié)果表明場(chǎng)區(qū)地下水的污染以單環(huán)芳烴和多環(huán)芳烴等特征污染物的污染最為突出,其中最為嚴(yán)重的是苯和萘2種污染物,利用Visual Modflow的MT3D模塊對(duì)這2種污染物進(jìn)行模擬。

2.4 模型的識(shí)別與驗(yàn)證

模型的識(shí)別與驗(yàn)證過(guò)程是整個(gè)模擬過(guò)程中的重要步驟,通常要進(jìn)行反復(fù)修改參數(shù)和調(diào)整某些源匯項(xiàng)才能達(dá)到較為理想的擬合結(jié)果。模型的識(shí)別和驗(yàn)證主要遵循以下2個(gè)原則:①模擬的地下水流場(chǎng)要與實(shí)際地下水流場(chǎng)基本一致,即要求地下水模擬等值線與實(shí)測(cè)地下水位等值線形狀相似;②識(shí)別的水文地質(zhì)參數(shù)要符合實(shí)際水文地質(zhì)條件。

根據(jù)以上原則,對(duì)模型進(jìn)行識(shí)別,先給出各個(gè)參數(shù)的范圍值,采用自動(dòng)和手動(dòng)相結(jié)合的方法,通過(guò)實(shí)際水位和計(jì)算水位的擬合分析判斷,如果計(jì)算水位與實(shí)測(cè)水位相差大,則根據(jù)參數(shù)的變化范圍和實(shí)際水位差值,再試給另一組參數(shù),直至二者擬合較好為止,通過(guò)反復(fù)調(diào)整參數(shù)和均衡要素,識(shí)別水文地質(zhì)條件,確定模型結(jié)構(gòu)、參數(shù)和均衡要素。

2.5 參數(shù)的選取

因?yàn)楸敬文M只考慮了對(duì)流和彌散的作用,所以溶質(zhì)運(yùn)移模型中系數(shù)的選取也只考慮彌散度、擴(kuò)散系數(shù)和孔隙度。彌散度是建立地下水溶質(zhì)運(yùn)移模型中最難以確定的系數(shù),一個(gè)重要原因是它的尺度效應(yīng),因此它的取值是當(dāng)今彌散理論研究的一個(gè)熱點(diǎn),如李國(guó)敏等[4]綜合世界范圍內(nèi)百余個(gè)水質(zhì)模型中所使用的縱向彌散度與空間尺度的關(guān)系,認(rèn)為水動(dòng)力彌散的尺度效應(yīng)具有分形特征,并給出了不同模型縱向彌散度尺度。綜上所述,本次溶質(zhì)運(yùn)移模擬中擬采用的相關(guān)參數(shù)為:有效孔隙度0.1～0.35;有效擴(kuò)散系數(shù)1.5×10-5m2/d;縱向彌散度10～30 m;垂向彌散度與縱向彌散度之比0.1;水平向彌散度與縱向彌散度之比0.1;降水入滲系數(shù)0.15;蒸發(fā)系數(shù)0.02。

3 模擬結(jié)果分析

根據(jù)目前該研究區(qū)污染物濃度,對(duì)苯和萘2種污染物分別進(jìn)行預(yù)測(cè),預(yù)測(cè)結(jié)果見(jiàn)圖1、圖2。

圖1 苯的污染羽遷移模擬結(jié)果

圖2 萘的污染羽遷移模擬結(jié)果

根據(jù)模型預(yù)測(cè)結(jié)果,由圖1(a)、(b)可見(jiàn),監(jiān)測(cè)點(diǎn)GW1苯的污染羽比監(jiān)測(cè)點(diǎn)GW2的污染羽大,運(yùn)移到第10年的時(shí)候,前者沿水流方向總的污染距離約150m,后者沿水流方向總的污染距離約80m,這是因?yàn)楸O(jiān)測(cè)點(diǎn)GW1的現(xiàn)狀污染程度比GW2的現(xiàn)狀污染程度大。從圖1(c)、(d)可見(jiàn),苯的污染羽在垂向上的遷移是從第一層慢慢向下面含水層滲透,在第6年的時(shí)候穿透中間的隔水層開(kāi)始滲入到承壓含水層。

從圖2(a)、(b)可見(jiàn),萘的污染羽隨時(shí)間逐漸增大,但是相對(duì)于同一監(jiān)測(cè)點(diǎn)GW1中萘的污染羽比苯的污染羽要小得多,運(yùn)移到第10a的時(shí)候萘沿水流方向的污染距離約100m,小于苯的150m,這是因?yàn)檩恋默F(xiàn)狀污染程度比苯的污染程度小。從圖2(c)、(d)萘的污染羽在垂向上的遷移狀況可見(jiàn),萘的污染羽也是從第一層慢慢向下遷移,但是由于現(xiàn)狀污染程度比較小,運(yùn)移10年時(shí)沒(méi)有穿透隔水層。

4 結(jié) 論

a.污染物在地下水遷移過(guò)程中,由于機(jī)械彌散和對(duì)流的共同作用,污染范圍的形狀由最初近似圓形逐漸變?yōu)殚L(zhǎng)軸方向與水流向一致的近似橢圓形。

b.污染物在含水層中隨著地下水而遷移,在遷移過(guò)程中隨著時(shí)間變化,苯和萘的污染羽中心都在隨著地下水的流向遷移,污染范圍逐漸增大。

c.溶質(zhì)運(yùn)移模擬顯示,研究區(qū)最主要的污染源可能是位于GW1點(diǎn)周?chē)膮^(qū)域,另外可能還存在多處小的污染源。

d.污染場(chǎng)地的污染源不完全是點(diǎn)源,而本次模擬時(shí)按照點(diǎn)源(有限面源)的方式作了簡(jiǎn)化處理,降低了污染風(fēng)險(xiǎn),建議盡早對(duì)該污染場(chǎng)地的地下水進(jìn)行治理,以防污染物通過(guò)運(yùn)移進(jìn)入到承壓含水層,造成更大的危害。

[1] 王克三.地下水污染及其監(jiān)測(cè)治理問(wèn)題[J].水文地質(zhì)工程地質(zhì),1998(5):46-47.

[2] 郭衛(wèi)星,盧國(guó)平.模塊化三維有限差分地下水流動(dòng)模型:美國(guó)地質(zhì)調(diào)查局水資源調(diào)查技術(shù)叢書(shū)之六[M].南京:南京大學(xué)出版社,1999.

[3] MC DONALD M G,HARBAUGH A W.A modular threedimensional finite-difference ground-water flow model:U S Geological Survey Techniques of Water-Resources Investigations[M].Reston,Virginia:USGS,1988.

[4] 李國(guó)敏,陳崇希.空隙介質(zhì)水動(dòng)力彌散尺度效應(yīng)的分形特征及彌散度初步估計(jì)[J].地球科學(xué)-中國(guó)地質(zhì)大學(xué)學(xué)報(bào),1995,20(4):405-409.