巖溶區(qū)某磷石膏堆場(chǎng)地下水滲流與污染物運(yùn)移數(shù)值模擬

張科正，陳 舟，陳長(zhǎng)生，趙貴清

(1.河海大學(xué) 地球科學(xué)與工程學(xué)院，南京 210000； 2.長(zhǎng)江三峽勘測(cè)研究院有限公司，武漢 430074)

0 序言

我國(guó)西南巖溶地區(qū)的部分磷石膏堆放場(chǎng)存在選址不當(dāng)、堆場(chǎng)設(shè)計(jì)不合理、防滲措施不徹底等造成的滲漏污染問(wèn)題[1-2]，主要污染物有磷、氟、硫酸鹽等。污染物對(duì)地下水的污染及對(duì)地表水體的排放對(duì)水環(huán)境和居民生活造成了較大影響。針對(duì)此類(lèi)問(wèn)題，已有眾多學(xué)者對(duì)滲漏類(lèi)型、特征和定量分析做出了相關(guān)研究，其中，趙偉麗等通過(guò)10個(gè)工業(yè)廢渣堆場(chǎng)的巖溶含水層介質(zhì)特征及污染質(zhì)運(yùn)移特征的分析,劃出4種巖溶滲漏污染的水動(dòng)力彌散類(lèi)型，分別為管道系統(tǒng)、管道裂隙、孔隙裂隙、采空區(qū)滲漏污染類(lèi)型[3]。彭展翔等對(duì)本研究堆場(chǎng)滲漏區(qū)、滲漏類(lèi)型、滲漏途徑、滲漏主方向等滲漏特征進(jìn)行了分析，確定滲漏主方向?yàn)榘l(fā)財(cái)洞方向[4]。在定量評(píng)價(jià)方面，由于巖溶地區(qū)地質(zhì)條件的復(fù)雜性，對(duì)于巖溶滲漏污染的數(shù)值模擬與預(yù)測(cè)研究還不夠完善，現(xiàn)主要預(yù)測(cè)模型有水動(dòng)力彌散疊加模型、稀釋混合模型、基于MATLAB BP神經(jīng)網(wǎng)絡(luò)模型、一維解析解法等。趙良杰等采用MODFLOW中排水溝(Drain)和河流(River)模塊概化模擬了巖溶管道特征，對(duì)比分析并指出Drain模塊中巖溶管道僅起排水作用，不允許管道內(nèi)水流流向周?chē)畬?，而River模塊允許巖溶管道與周?chē)畬舆M(jìn)行水量交換，因而相對(duì)于Drain模塊，應(yīng)用River模塊概化模擬巖溶管道更加精確[5]?？偨Y(jié)可知，巖溶區(qū)地下水污染物遷移的數(shù)值模擬鮮有報(bào)道，本文擬以某磷石膏堆放場(chǎng)的污染滲漏為例，運(yùn)用Groundwater Modeling System(GMS)軟件搭建了該區(qū)域的三維地下水滲流場(chǎng)和污染物遷移的數(shù)值模型，以期查明污染物(磷)的分布遷移特征，為堆場(chǎng)的污染治理工作和巖溶水污染預(yù)測(cè)評(píng)價(jià)提供依據(jù)。

1 研究區(qū)基本地質(zhì)條件

2 堆場(chǎng)滲漏的數(shù)值模擬

2.1 模型建立

2.1.1 水文地質(zhì)概念模型

圖1 研究區(qū)模型邊界及類(lèi)型Figure 1 Study area model boundary and type

圖2 堆場(chǎng)—發(fā)財(cái)洞剖面圖Figure 2 Dump-Facaidong section

2)含水層結(jié)構(gòu)。研究區(qū)屬裸露型巖溶山區(qū)，碳酸鹽巖出露面積占70%以上，地下水類(lèi)型為裂隙水和巖溶水為主，區(qū)域內(nèi)的斷層、天窗使得上下含水層相互串通，成為統(tǒng)一的混合含水層。

3)含水層的水力特征。結(jié)合研究區(qū)地質(zhì)資料、地下水補(bǔ)、徑、排特征及巖溶發(fā)育特征，將研究區(qū)概化為非均質(zhì)各向異性三維地下水流場(chǎng)。研究區(qū)出露泉點(diǎn)水位多年動(dòng)態(tài)變化較穩(wěn)定，所以本次模擬將建立枯水期研究區(qū)地下水穩(wěn)定流場(chǎng)。地下水溶質(zhì)主要通過(guò)對(duì)流-彌散作用隨地下水流一起運(yùn)移。

2.1.2 地下水滲流與溶質(zhì)運(yùn)移耦合數(shù)學(xué)模型

1)地下水滲流數(shù)學(xué)模型。通過(guò)概化得到的非均質(zhì)各向異性等效連續(xù)介質(zhì)模型，地下水穩(wěn)定運(yùn)動(dòng)數(shù)學(xué)模型為

(1)

式中：H—地下水水頭，m；

Kx，Ky，Kz—各向異性主滲透系數(shù)，m/d；

Γ1—模擬區(qū)域第一類(lèi)邊界；

Γ2—模擬區(qū)域第二類(lèi)邊界；

H0(x,y,z)—含水層初始水頭m；

Hг(x,y,z)—第一類(lèi)邊界條件邊界水頭，m；

q0(x,y,z)—第二類(lèi)邊界單位面積過(guò)水?dāng)嗝嫜a(bǔ)給流量，m2/d；

ε—源匯項(xiàng)強(qiáng)度(包括開(kāi)采強(qiáng)度等)，1/d；

Ω—滲流區(qū)域。

(2)地下水溶質(zhì)運(yùn)移數(shù)學(xué)模型。

(2)

式中：n—介質(zhì)孔隙度，無(wú)量綱；

C—組分的濃度，mg/L；

t—時(shí)間，d；

Dij—水動(dòng)力彌散系數(shù)張量，m2/d；

υi—地下水滲流速度張量，m/d；

qs—源和匯，1/d；

Cs—源或匯水流中組分的濃度，mg/L；

C0(x,y,z,t0)—已知濃度分布，mg/d；

Ω—模擬區(qū)域；

Γ1—通量邊界；

fi(x,y,z,t)— 邊界上已知的彌散通量函數(shù)。

2.2 參數(shù)選擇

2.2.1 滲透系數(shù)

本模擬結(jié)合鉆孔壓水試驗(yàn)，并對(duì)比了相關(guān)的工程實(shí)例，進(jìn)行了滲透系數(shù)的選取。考慮巖體各向異性，其滲透系數(shù)按照三維取值，各巖性對(duì)應(yīng)滲透系數(shù)見(jiàn)表1。

2.2.2 降雨入滲系數(shù)

模型輸入時(shí)的參數(shù)Recharge rate項(xiàng)等于降水強(qiáng)度與降雨入滲系數(shù)的積，降雨入滲系數(shù)在枯水年平均有效降雨入滲系數(shù)為0.35，參考其他文獻(xiàn)和區(qū)域的地形和巖性特征，區(qū)域的降入入滲分區(qū)及參數(shù)取值可見(jiàn)表2。

表1 研究區(qū)巖體滲透系數(shù)賦值

表2 不同巖性對(duì)應(yīng)降雨入滲系數(shù)

2.3 模型識(shí)別與校正

研究區(qū)地下水?dāng)?shù)值模型采用GMS10.2軟件進(jìn)行數(shù)值離散，采用等間距有限差分的離散方法對(duì)模擬區(qū)進(jìn)行自動(dòng)剖分，網(wǎng)格單元行×列×層為200×200×5，模擬區(qū)的三維尺度在X方向上長(zhǎng)度為8 349.66 m，Y方向上長(zhǎng)度為7 492.21 m，Z方向的長(zhǎng)度為500m，網(wǎng)格剖分及巖性三維分區(qū)見(jiàn)圖3。初始水頭設(shè)置為網(wǎng)格頂部高程。基于GMS創(chuàng)建各巖溶通道并將其屬性設(shè)置為River，通道的兩個(gè)端點(diǎn)高程設(shè)置為巖溶管道實(shí)際高程。

圖3 網(wǎng)格剖分及巖性三維分區(qū)Figure 3 Mesh generation and lithologic 3D partitioning

在數(shù)值模擬過(guò)程中，由于地質(zhì)條件的概化、模型參數(shù)的選取等都會(huì)造成模擬地下水流場(chǎng)與實(shí)際地下水流場(chǎng)不符。因此，對(duì)模型進(jìn)行識(shí)別與校正是數(shù)值模擬的一個(gè)重要環(huán)節(jié)。模型識(shí)別與校正不能片面追求擬合效果而任意調(diào)整模型參數(shù)，甚至改變含水層的邊界條件及源(匯)項(xiàng)等原始數(shù)據(jù)，這種做法不僅容易使模型脫離實(shí)際的水文地質(zhì)條件，而且還造成模型的多解性。研究區(qū)的面積很大，鉆孔水位資料相對(duì)缺乏，但巖溶泉點(diǎn)有50個(gè)，根據(jù)泉水點(diǎn)位的出露高程插值可得研究區(qū)的初始流場(chǎng)圖，見(jiàn)圖4。此外，模型的識(shí)別利用了2016年踏勘所調(diào)查的16個(gè)巖溶泉出露水位作為監(jiān)測(cè)井，模型的初始水頭設(shè)置為網(wǎng)格頂板高程[8]。

圖4 研究區(qū)地下水初始流場(chǎng)圖Figure 4 Study area groundwater initial flow field

利用GMS軟件中的參數(shù)自動(dòng)識(shí)別模塊(PEST)與手動(dòng)調(diào)參相結(jié)合的方法進(jìn)行模型校正與參數(shù)反演，識(shí)別后模型的地下水流場(chǎng)和監(jiān)測(cè)井水位計(jì)算誤差圖可見(jiàn)圖5和圖6。

圖5 研究區(qū)地下水場(chǎng)擬合圖Figure 5 Study area groundwater field fitting diagram

圖6 各觀測(cè)井實(shí)測(cè)值與計(jì)算值的擬合效果圖Figure 6 Fitting effect diagram of measured and computed values from observation wells

由結(jié)果可知，16口觀測(cè)井的觀測(cè)水位與計(jì)算水位的平均殘差值為-0.253 7，擬合誤差的平均絕對(duì)值為0.586 4，擬合誤差小于0.8m的觀測(cè)井占總的井?dāng)?shù)量的87%以上，可見(jiàn)實(shí)測(cè)水位與計(jì)算水位之間達(dá)到了較好的擬合效果，參數(shù)校正結(jié)果比較可靠，參數(shù)的校正結(jié)果滿(mǎn)足精度要求。此外，模型River模塊反應(yīng)出了巖溶通道中水流集中排泄的特征，通道區(qū)域周?chē)牡叵滤畢R集到管道中以較快的流量進(jìn)行排泄，所建立的模型切合實(shí)際。

2.4 堆場(chǎng)滲漏預(yù)測(cè)

2.4.1 溶質(zhì)運(yùn)移模型

本研究利用GMS軟件中MT3DMS模塊進(jìn)行污染物遷移的模擬和預(yù)測(cè)。模擬參數(shù)中的孔隙度可見(jiàn)表3。

縱向彌散度由Xu和Eckstein(1995)[7]回歸分析得到公式：

aL=0.83×(logLs)2.414

=0.83×(log 2 000)2.414

=14.8m

(3)

式中：aL—縱向彌散系數(shù)；Ls—觀測(cè)尺度。橫向及垂向彌散度為縱向彌散度的十分之一。

表3 不同巖性對(duì)應(yīng)的初始和識(shí)別后孔隙度

污染滲漏工況設(shè)置如下：由文獻(xiàn)[3]可知，擺紀(jì)渣場(chǎng)的滲漏量約為300 m3/d，滲漏液中磷的濃度約為3 000 mg/L，本次模擬按照滲漏液的10%會(huì)入滲地下水的工況進(jìn)行堆場(chǎng)滲漏模擬。

污染物遷移模型的識(shí)別和校正的數(shù)據(jù)利用了2016年發(fā)財(cái)洞處污染物(磷)監(jiān)測(cè)數(shù)據(jù)。觀測(cè)值與計(jì)算值的誤差圖可見(jiàn)圖7?？紫抖茸R(shí)別后的參數(shù)值可見(jiàn)表4，彌散度識(shí)別后的數(shù)值為12.5 m。

圖7 磷濃度計(jì)算值與觀測(cè)值比較Figure 7 Comparison of phosphorus measured and observed concentration values

如圖7所示，發(fā)財(cái)洞處一年內(nèi)磷濃度的計(jì)算值與觀測(cè)值較為一致，數(shù)據(jù)間相關(guān)系數(shù)為0.997，可認(rèn)為模型較為真實(shí)地反映出地下水中污染物的運(yùn)移情況。

2.4.2 滲漏預(yù)測(cè)

為預(yù)測(cè)持續(xù)滲漏情況下污染物的濃度變化，利用識(shí)別后的模型預(yù)測(cè)了堆場(chǎng)連續(xù)滲漏后1 a、5 a、10 a和20 a后污染物在含水層的運(yùn)移情況，具體可見(jiàn)圖8。此外，為了對(duì)比含水層污染物遷移的垂直分布和三維分布情況，以污染物遷移10 a后為例展示了相關(guān)圖件，具體可見(jiàn)圖9。

由預(yù)測(cè)可知，污染物持續(xù)滲漏十年后發(fā)財(cái)洞處總磷濃度可達(dá)400 mg/L左右，磷石膏堆場(chǎng)～發(fā)財(cái)洞方向?yàn)槎褕?chǎng)的主要滲漏方向，老落凼巖溶管道與擺郎發(fā)財(cái)洞巖溶管道為主要滲漏通道，滲漏污染物隨地下水水流方向河流處排泄， 擺郎發(fā)財(cái)洞巖溶管道水系統(tǒng)為主導(dǎo)方向，形成堆場(chǎng)至發(fā)財(cái)洞污染物濃度梯次遞減形分布的規(guī)律。

圖8 不同時(shí)段污染物遷移模擬圖Figure 8 Pollutant migration simulated diagram of different time intervals

圖9 10 a后污染物運(yùn)移預(yù)測(cè)三維及剖面圖Figure 9 Predicted pollutant migration 3D map and section after 10 years

3 結(jié)語(yǔ)

為了定量研究巖溶區(qū)某磷石膏堆放場(chǎng)地下水滲流和污染物滲漏污染問(wèn)題，本文利用地下水?dāng)?shù)值模擬軟件GMS建立了巖溶區(qū)地下水流場(chǎng)及污染物(磷)在巖溶含水層的運(yùn)移情況，在模型的識(shí)別和校正后預(yù)測(cè)了1～20 a后污染物的濃度分布。

由結(jié)果可知：River模塊能較好的模擬巖溶區(qū)的巖溶管道滲流。堆場(chǎng)至發(fā)財(cái)洞方向?yàn)檠芯繀^(qū)的滲漏主方向，老落凼巖溶管道與擺郎發(fā)財(cái)洞巖溶管道為主要滲漏通道，形成堆場(chǎng)至發(fā)財(cái)洞污染物濃度梯次遞減分布的規(guī)律。預(yù)測(cè)持續(xù)滲漏情況下，10 a后發(fā)財(cái)洞處總磷濃度可達(dá)400mg/L左右。

由于基礎(chǔ)資料不夠充足且?guī)r溶區(qū)的部分特征未能在GMS模型中體現(xiàn)，數(shù)值模擬的結(jié)果還存在很大的完善空間；溶質(zhì)運(yùn)移模擬中，本文只考慮了污染物運(yùn)移過(guò)程中的對(duì)流、彌散作用，并未考慮吸附作用及不同管道特征的影響，這也是相關(guān)工作需進(jìn)一步深入之處。

參考文獻(xiàn):

[1]段先前,韋俊發(fā),丁堅(jiān)平.貴州某磷石膏堆場(chǎng)滲漏污染評(píng)價(jià)[J].資源環(huán)境與工程, 2008, 22(2):218-221.

DUAN X Q,WEI J F, et al.Evaluation on seepage pollution of a phosphogypsum yard in GuiZhou Province[J].Resources Environment And Engineering, 2008, 22(2):218-221.(in Chinese)

[2]鄭文成,賈洪彪,郭明,等.平壩磷石膏堆場(chǎng)滲漏分析與污染防治對(duì)策[J].中國(guó)巖溶, 2010, 29(1):75-80.

ZHENG W C, JIA H B, GUO M, et al. Analysis and countermeasures of pollution prevention and control of Pingba phosphogypsum yard leakage[J].Carsologica Sinica,2010,29(1):75-80.(in Chinese)

[3]趙偉麗,褚學(xué)偉,董毓,等.巖溶含水層滲漏污染彌散類(lèi)型分析——以貴州廢渣堆場(chǎng)污染為例[J].地下水, 2011, 33(2):6-7.

ZHAO W L,CHU X W, el al. Analysis on leakage pollution dispersion type in karst aquifer[J].Groundwater, 2011,33(2):6-8.(in Chinese)

[4]彭展翔, 褚學(xué)偉.擺紀(jì)磷石膏堆場(chǎng)滲漏分析[J].地下水, 2012, 34(5).

PENG Z X, ZHU X W. Analysis on the leakage of the phosphorus gypsum storage yard in Baiji[J].Groundwater,2012,34(5).(in Chinese)

[5]趙良杰,夏日元,楊楊,等.基于MODFLOW的巖溶管道水流模擬方法探討與應(yīng)用[J].中國(guó)巖溶, 2017, 36(3):346-351.

ZHAO L J,XIA R Y, el al. Discussion and application of simulation methods for karst conduit flow based on MODFLOW[J].Carsolgica Sinica,2017,36(3):346-351.(in Chinese)

[6]張景輝, 孫曙光. 地下水非穩(wěn)定運(yùn)動(dòng)的數(shù)學(xué)模型及應(yīng)用[J].上海國(guó)土資源, 1981(4):43-52.

ZHANG J H,SUN S G. The mathematical model of unsteady groundwater movement and application[J].Shanghai geological,1981, 41-50.(in Chinese)

[7]M Xu,Y & Eckstein. Use of weighted least-squares method in evaluation of the relationship between dispersivity and field scale[J].Groundwater, 1995, 33(6):905-908.

[8]陳舟,趙貴清,王志光,等.巖溶區(qū)某磷石膏堆放場(chǎng)滲漏特征分析[J].水文地質(zhì)工程地質(zhì),2017,44(2):144-150.

CHEN Z,ZHAO G Q, el al. Leakage Characteristics of a Phosphorus Gypsum Storage Site in Karst Area[J].Hydrogeology And Engineering Geology,2017,44(2):144-150.(in Chinese)

[9]祝曉彬.地下水模擬系統(tǒng)(GMS)軟件[J].水文地質(zhì)工程地質(zhì), 2003, 30(5):53-55.

ZHU X L. Groundwater modeling system(GMS) software[J].Hydrogeology and engineering geology,2003,5:53-55.(in Chinese)

[10]譚文清,孫春,胡婧敏,等.GMS在地下水污染質(zhì)運(yùn)移數(shù)值模擬預(yù)測(cè)中的應(yīng)用[J].東北水利水電, 2008, 26(5):54-55.

TAN W Q,SUN C, el al. Application of GMS in simulation of pollutants migration for groundwater[J].Water Resources And Hydropower of Northeast China,2008,5:54-55.(in Chinese)

[11]駱祖江,于丹丹,劉艷,等.江蘇長(zhǎng)強(qiáng)鋼鐵廠地下水環(huán)境影響評(píng)價(jià)[J].中國(guó)煤炭地質(zhì), 2014, 26(5):34-37.

LUO Z J,YU D D, el al. Jiangsu Changqiang Steel Mill Groundwater Environment Impact Assessment[J].Coal Geology of China,2014,26(5):34-37.(in Chinese)