巖溶區(qū)某錳礦尾礦回填區(qū)地下水環(huán)境數(shù)值模擬預(yù)測

楚敬龍，林星杰，苗 雨，劉楠楠，譚海偉，劉 芳

(北京礦冶科技集團(tuán)有限公司，北京 100160)

礦業(yè)開采帶來巨大的經(jīng)濟(jì)利益，但同時(shí)也產(chǎn)生大量的礦業(yè)廢棄地，導(dǎo)致嚴(yán)重的生態(tài)破壞和環(huán)境污染[1]。這些礦業(yè)廢棄地包括尾礦庫、廢石場、采空區(qū)、堆浸場及采空回填區(qū)等，其內(nèi)所含尾砂、廢石等廢棄物常含有與礦業(yè)開采相關(guān)的重金屬等污染物，這些廢棄地若無防滲措施或防滲措施失效，在降雨淋溶或自身含水滲漏的情況下，可能對周圍地下水產(chǎn)生污染。在我國南方巖溶區(qū)存在大量采礦廢棄地，由此產(chǎn)生的巖溶滲漏污染是經(jīng)常遇到且必須解決的重大環(huán)境問題，而滲漏污染物造成的地下水及地表水污染問題，已對水環(huán)境和居民生活造成了較大影響[2-3]。

鑒于此，眾多學(xué)者對污染物在巖溶含水層中的遷移行為及巖溶水的數(shù)值模擬技術(shù)等方面做出了相關(guān)研究。趙偉麗等[4]劃分了4種巖溶滲漏的水動(dòng)力彌散類型，分別為管道型(巖溶強(qiáng)發(fā)育)、管道裂隙型(巖溶強(qiáng)發(fā)育)、裂隙孔隙型(巖溶中等—弱發(fā)育)、采空區(qū)型(類似地表徑流)，為巖溶含水層污染預(yù)測評價(jià)提供了依據(jù)；陸嘉等[5]以氯化物為預(yù)測因子，利用地下水?dāng)?shù)值模擬軟件Visual Modflow對煤層氣采出水進(jìn)入地下水后污染物的遷移路徑和污染范圍進(jìn)行模擬，結(jié)果表明污染物基本不會(huì)對周圍居民飲水健康造成影響；張科正等[3]針對巖溶區(qū)某磷石膏堆放場滲漏污染問題進(jìn)行研究，運(yùn)用地下水?dāng)?shù)值模擬軟件GMS模擬了巖溶區(qū)地下水流場及污染物(磷)在巖溶含水層的運(yùn)移情況，確定了污染物的濃度分布；胡建青等[6]結(jié)合準(zhǔn)格爾礦區(qū)水文地質(zhì)特點(diǎn)，利用地下水?dāng)?shù)值模擬方法，定量研究了深層巖溶水資源量的變化，將研究區(qū)內(nèi)的寒武、奧陶系碳酸鹽巖地下水系統(tǒng)概化為巖溶裂隙普遍發(fā)育(已構(gòu)成一細(xì)密網(wǎng)絡(luò)狀的含水結(jié)構(gòu))、地下水運(yùn)移符合達(dá)西定律的非均質(zhì)各向異性的三維非穩(wěn)定流系統(tǒng);趙良杰等[7]為了描述巖溶多重含水介質(zhì)存在非達(dá)西流問題，研究了管道流CFP模型在巖溶區(qū)數(shù)值模擬中的應(yīng)用，該方法具有一定通用性，可為實(shí)際應(yīng)用提供指導(dǎo)。

廣西覆蓋型巖溶區(qū)某堆積型錳礦賦存于第四系殘坡積層中，分布于低洼地或緩坡部位，礦體埋藏淺。主開采區(qū)包括6個(gè)堆積錳礦體，分為三期進(jìn)行露天開采，其中第一期包含①號(hào)礦體，服務(wù)年限9.6年；第二期包含②、③、④、⑤共四個(gè)礦體，服務(wù)年限12.1年；第三期包含⑥號(hào)礦體，服務(wù)年限2.3年。選廠前2.5年所產(chǎn)生的尾礦脫水后運(yùn)至CK1～CK3舊采空區(qū)回填，2.5年后所產(chǎn)生的尾礦脫水后運(yùn)至新產(chǎn)生的露天采空區(qū)進(jìn)行回填。本文以該堆積型錳礦的尾礦回填區(qū)污染滲漏為例，運(yùn)用地下水?dāng)?shù)值模擬軟件GMS，建立研究區(qū)地下水水流及溶質(zhì)運(yùn)移數(shù)值模型，預(yù)測污染物在巖溶含水層中的時(shí)空分布情況，評價(jià)尾礦回填區(qū)污染滲漏對巖溶區(qū)地下水敏感目標(biāo)的影響，以期為巖溶區(qū)的地下水污染預(yù)測及防治提供一定的研究基礎(chǔ)。

1 研究區(qū)概況

1.1 自然地理?xiàng)l件

本研究范圍包含各礦體歷史遺留采坑在內(nèi)，西側(cè)以A河流為界，北側(cè)以B河流為界，東側(cè)沿著地表水分水嶺山脊線(涵蓋石炭系大埔組(C2pd)巖溶較發(fā)育區(qū))，南側(cè)到達(dá)C村南部一帶，為一個(gè)相對完整的水文地質(zhì)單元，總面積約14.5 km2。

區(qū)內(nèi)主體構(gòu)造線呈近南東—北西走向，山脈排列、水系發(fā)育均受構(gòu)造控制。地勢上北東側(cè)較高，西側(cè)較低為河谷，屬強(qiáng)烈剝蝕淺切割的丘陵低山地形，地貌形態(tài)為巖溶峰林、洼地、孤峰等，高程多在+150～+400 m，相對高差100～150 m，海拔最低為+122～+127 m(A河)，靠近村屯地段地形較平緩。年平均降雨量1 980 mm，A河從礦區(qū)西部自南南東向北北西流過，河床兩岸谷地寬闊平坦，多為水田或耕地。B河為A河支流，在研究區(qū)北部由東向西注入A河。

研究區(qū)內(nèi)各礦體下游分布有1個(gè)集中式飲用水源，19口分散式機(jī)/民井，5處泉水取水點(diǎn)，水源均為石炭系大埔組(C2pd)巖溶水，可作為地下水保護(hù)目標(biāo)。見圖1。

圖1 地下水保護(hù)目標(biāo)圖

1.2 水文地質(zhì)特征

研究區(qū)含水巖組分為松散巖類含水巖組、碳酸鹽巖含水巖組、基巖構(gòu)造裂隙含水巖組三類。1)松散巖類含水巖組主要分布于A河兩岸，一般分三層，上部為亞黏土，巖屑亞黏土及巖屑質(zhì)亞黏土，含少量錳礦塊；中部為堆積錳礦體，由氧化錳礦塊與亞黏土、黏土及巖屑混雜組成；下部巖性為黏土，巖屑黏土，巖屑亞黏土。該層枯季一般不含水，雨季則常具季節(jié)性的含水特性，透水性差，為相對隔水層。2)碳酸鹽巖含水巖組含裂隙溶洞水，分為中等富水性和弱富水性兩個(gè)等級。中等富水性的巖溶水賦存于泥盆系上統(tǒng)融縣組(D3r)、唐家灣組(D2t)及石炭系大埔組(C2pd)、英塘組(C1yt)地層中，巖性為灰?guī)r、白云巖，分布于礦區(qū)北東部，主要接受大氣降雨補(bǔ)給，巖溶較發(fā)育；弱富水性的巖溶水賦存于二疊系棲霞組(P2q)地層中，主要巖性為灰?guī)r，常含燧石結(jié)核及條帶，局部白云質(zhì)較多，整體巖溶弱發(fā)育。3)基巖構(gòu)造裂隙含水巖組含基巖構(gòu)造裂隙水，主要賦存于中泥盆統(tǒng)信都組(D2x)、中二疊統(tǒng)孤峰組(P2g)及上二疊統(tǒng)龍?zhí)督M(P3l)等地層中，富水性弱，對礦床充水影響弱。

圖2 水位擬合圖

研究區(qū)地下水以中等富水性的碳酸鹽巖裂隙溶洞水為主。礦區(qū)處于石炭系大埔組(C2pd)碳酸鹽巖裂隙溶洞水中等富水性的巖溶水地區(qū)，巖性為灰?guī)r、白云巖。巖溶水主要接受大氣降雨補(bǔ)給，順地形自東向西徑流，主要運(yùn)移于巖溶裂隙之中，最終排泄于A河。因上覆第四系透水性差，有相對隔水作用，因此巖溶水具承壓性質(zhì)，部分向上越流至第四系，在低洼處呈泉水排泄出地表。

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

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

礦區(qū)堆積型錳礦賦存于第四系殘坡積層中，分布于山前低洼地或緩坡部位，礦體埋藏淺，其中①號(hào)、②號(hào)、③號(hào)、④號(hào)、⑤號(hào)、⑥號(hào)礦體為凹陷開采，并作為尾礦的回填區(qū)，主要影響第四系下伏石炭系大埔組(C2pd)巖溶水，因此考慮堆積型錳礦的凹陷開采形成的尾礦回填區(qū)對礦區(qū)巖溶地下水的影響。

模擬范圍：北邊界為B河，概化為定水頭邊界；東邊界為石炭系大埔組(C2pd)巖溶較發(fā)育區(qū)與二疊系棲霞組(P2q)巖溶弱發(fā)育區(qū)地層分界線，概化為弱流量邊界；南邊界為垂直于地下水等水位線，概化為零流量邊界；西邊界為A河，概化為定水頭邊界；上邊界為第四系與石炭系大埔組的交界面，接受大氣降水補(bǔ)給(以東側(cè)裸露山區(qū)的大氣降水補(bǔ)給為主)，并越流補(bǔ)給上層第四系相對隔水層，同時(shí)以地下水側(cè)向徑流及人工開采的形式進(jìn)行地下水排泄；下邊界為深部巖溶弱發(fā)育的相對隔水層。

模擬范圍之內(nèi)的地下水主要在石炭系大埔組(C2pd)碳酸鹽巖類裂隙溶洞水含水層中流動(dòng)。根據(jù)礦區(qū)水文地質(zhì)條件、水文地質(zhì)平面圖、剖面圖及鉆孔柱狀圖，綜合考慮，將模擬范圍內(nèi)的含水層概化為一層，即石炭系大埔組(C2pd)碳酸鹽巖類裂隙溶洞水承壓含水層，巖性主要為灰?guī)r、白云巖，平均發(fā)育厚度在第四系相對隔水層下約20～30 m。

在巖溶相對中等—弱發(fā)育區(qū)，巖溶介質(zhì)溶孔、溶蝕裂隙發(fā)育相對來說比較穩(wěn)定，含水層中水力聯(lián)系相對比較緊密，巖溶介質(zhì)的發(fā)育具有相對均勻性，巖溶滲漏可認(rèn)為符合均勻介質(zhì)的運(yùn)移規(guī)律，濃度水動(dòng)力彌散曲線與理論曲線相符合[4]。因此研究區(qū)地下水系統(tǒng)可概化為非均質(zhì)各向異性、空間單層結(jié)構(gòu)穩(wěn)定地下水流系統(tǒng)。

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

1)地下水流數(shù)學(xué)模型

基于2.1中的水文地質(zhì)概念模型，可以用以下數(shù)學(xué)模型描述地下水流系統(tǒng)：

式中：Ω—地下水滲流區(qū)域；H—地下水水頭，m；S1—模型的第一類邊界；S2—模型的第二類邊界；kxx、kyy、kzz分別表示x、y、z主方向的滲透系數(shù)，m/d；W—源匯項(xiàng)，包括降水入滲補(bǔ)給等，m3/d；H1(x,y,z)—第一類邊界已知地下水水頭函數(shù)，m；q(x,y,z)—第二類邊界單位面積流量函數(shù)，m3/d；n—邊界S2上的外法線方向[8]。

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

式中：C—地下水中組分的溶解相濃度，mg/L；θ—地層介質(zhì)的孔隙度，無量綱；t—時(shí)間，d；xi—沿直角坐標(biāo)系軸向的距離，m；Dij—水動(dòng)力彌散系數(shù)張量，m2/d；vi—地下水平均實(shí)際流速，m/d；qs—單位體積含水層流量，代表源和匯，m3/d；Cs—源或匯水流中組分的濃度，mg/L；∑Rn—化學(xué)反應(yīng)項(xiàng)，mg/(m3·d)；Ω—溶質(zhì)運(yùn)移的區(qū)域；C0—組分的初始濃度，mg/L；Г2—Neumann邊界。

2.3 數(shù)值模型識(shí)別與驗(yàn)證

在建立水文地質(zhì)概念模型、數(shù)學(xué)模型的基礎(chǔ)上，運(yùn)用基于有限差分法的GMS軟件建立了地下水流數(shù)值模型。應(yīng)用MODFLOW模塊，將模擬區(qū)剖分為50 m×50 m的矩形網(wǎng)格，模型剖分共計(jì)12 960個(gè)計(jì)算單元。源匯項(xiàng)包括降雨入滲、地下水側(cè)向徑流、地下水與地表水的相互補(bǔ)排、人工開采等，均由相應(yīng)軟件包進(jìn)行刻畫，然后分配到相應(yīng)單元格。

模型的識(shí)別與檢驗(yàn)是整個(gè)模擬工作中極為重要的一步，是一個(gè)不斷調(diào)節(jié)水文地質(zhì)參數(shù)，使模型結(jié)果盡可能與實(shí)際調(diào)查資料相吻合的過程。根據(jù)研究區(qū)專項(xiàng)水文地質(zhì)勘查提供的滲透系數(shù)及其他水文地質(zhì)參數(shù)的建議值，運(yùn)行數(shù)值模型，并根據(jù)抽水試驗(yàn)數(shù)據(jù)及實(shí)際水文地質(zhì)條件進(jìn)行反復(fù)調(diào)參，得到最終水文地質(zhì)參數(shù)分布情況，并確定地下水流場。流場擬合情況見圖2所示。由圖2可知，對模擬區(qū)調(diào)查的22個(gè)水位觀測點(diǎn)進(jìn)行擬合，誤差小于1 m的數(shù)據(jù)為21個(gè)，占比為95.5%，誤差在1～2 m之間數(shù)據(jù)為1個(gè)(最西北側(cè))，占比為4.5%，反應(yīng)了水位擬合點(diǎn)處的水位計(jì)算值和實(shí)際觀測值匹配較好，所建立的地下水流數(shù)值模型基本達(dá)到精度要求。

3 預(yù)測結(jié)果及分析

3.1 預(yù)測情景

本次預(yù)測情景針對礦區(qū)①～⑥號(hào)礦體回填區(qū)尾礦淋濾液滲漏進(jìn)行模擬預(yù)測。本研究區(qū)的主要含水層為第四系下伏石炭系大埔組碳酸鹽巖含水巖組，而錳礦開采位于第四系地層，形成的采空區(qū)回填尾礦，在降雨淋溶條件下，尾礦淋濾液會(huì)滲透至下伏的巖溶含水層。尾礦回填是逐步開展的，即先回填已有舊采空區(qū)(CK1～CK3)，同時(shí)進(jìn)行礦體開采，舊采空區(qū)填滿后，再回填不斷形成的新采空區(qū)，因而尾礦的降雨淋濾液是逐步產(chǎn)生并隨地下水往下游遷移的。根據(jù)實(shí)際情況，將各個(gè)尾礦回填區(qū)污染源概化為補(bǔ)給濃度連續(xù)面源污染。

第一期開采共9.6年：0～2.5年，CK1～CK3等舊采坑全部填滿尾礦；2.5～9.6年，①號(hào)礦體全部開采完并填滿尾礦。第二期開采共12.1年：9.6～10.3年，②號(hào)礦體全部開采完并填滿尾礦；10.3～11.8年，③號(hào)礦體全部開采完并填滿尾礦；11.8～20.7年，④號(hào)礦體全部開采完并填滿尾礦；20.7～21.7年，⑤號(hào)礦體全部開采完并填滿尾礦。第三期開采共2.3年：21.7～24年，⑥號(hào)礦體全部開采完并填滿尾礦。24～30年，為服務(wù)期滿后。以此作為時(shí)間序列，對本情景進(jìn)行污染物遷移模擬預(yù)測。

在預(yù)測因子選擇方面，唐文杰等[9]通過對錳礦區(qū)土壤重金屬全量及其形態(tài)分布的分析，認(rèn)為在錳礦區(qū)Cd與Mn這2種重金屬可利用態(tài)與潛在可利用態(tài)較高，是錳礦區(qū)污染土壤修復(fù)的主要修復(fù)對象。陳春強(qiáng)等[1]對廣西荔浦、平樂、來賓3個(gè)錳礦恢復(fù)區(qū)的土壤和農(nóng)作物進(jìn)行調(diào)查和重金屬含量分析，得出3個(gè)錳礦恢復(fù)區(qū)土壤的主要重金屬污染因子是Cd與Mn。根據(jù)現(xiàn)有老采坑礦坑水濃度監(jiān)測值可知，錳的最大濃度值為0.18 mg/L，超標(biāo)0.8倍，其余監(jiān)測因子均滿足標(biāo)準(zhǔn)要求，因此本次預(yù)測選取Mn作為預(yù)測因子。

3.2 結(jié)果及分析

各個(gè)污染源(CK1～CK3、①～⑥號(hào)礦體回填區(qū))按不同的濃度時(shí)間序列輸入模型，模擬期為30年，得到Mn擴(kuò)散預(yù)測結(jié)果。圖3給出了滲漏后9.6年、21.7年、24年、30年，Mn在巖溶含水層中的遷移范圍。按照《地下水質(zhì)量標(biāo)準(zhǔn)》(GB/T14848-2017)中Ⅲ類標(biāo)準(zhǔn)，該污染范圍的外邊界即為0.1 mg/L等濃度線。表1為Mn在巖溶含水層中的擴(kuò)散范圍及影響的保護(hù)目標(biāo)[10]。

由預(yù)測結(jié)果可知：

1)在CK1～CK3、①～⑥號(hào)礦體采空區(qū)逐步進(jìn)行尾礦回填的過程中，滲濾液逐步產(chǎn)生，并隨地下水由東向西遷移，在30年的模擬期內(nèi)，污染羽前鋒在礦體外最大遷移距離為300 m。其中④號(hào)礦體的污染物遷移范圍在24年時(shí)首先到達(dá)A河。①～③號(hào)礦體的污染物，由于其下游有S43、S43-1、S44、S46、S53取水井的存在，而被限制在一定范圍內(nèi)，沒有繼續(xù)往下游擴(kuò)散。⑤號(hào)礦體及其周圍未出現(xiàn)超標(biāo)污染羽。⑥號(hào)礦體污染物未擴(kuò)散到集中式水源保護(hù)區(qū)，鑒于集中式水源保護(hù)區(qū)的敏感性，應(yīng)對集中式水源保護(hù)區(qū)進(jìn)行地下水水質(zhì)監(jiān)測，并采取相應(yīng)預(yù)防措施，防止集中式水源地受到污染。

2)在污染物遷移擴(kuò)散過程中，受影響的保護(hù)目標(biāo)為：①號(hào)礦體下游：S44(泉水)、S46(泉水)、S53(機(jī)井)；②號(hào)礦體下游：S43(民井)、S43-1(機(jī)井)；④號(hào)礦體下游：S37(泉水)、S40(泉水)；⑥號(hào)礦體下游：S30(民井)、S32(民井)，其取水層位均為巖溶水含水層。一旦這些村屯的飲用水井及泉水受到影響，應(yīng)立即采取措施，啟動(dòng)應(yīng)急預(yù)案，確保村民取水不受影響。

圖3 不同時(shí)間的Mn遷移擴(kuò)散圖

4 結(jié)論

1)本文以Mn為預(yù)測因子，應(yīng)用地下水模擬軟件GMS對研究區(qū)進(jìn)行了水文地質(zhì)概化，建立了研究區(qū)地下水?dāng)?shù)值模型，預(yù)測了錳礦尾礦回填區(qū)污染物滲漏對地下水環(huán)境的影響。結(jié)果表明，污染物隨地下水由東向西往A河遷移，模擬期內(nèi)，有9個(gè)分散式居民飲用水井及泉水受到影響，但不會(huì)影響到⑥號(hào)礦體北側(cè)的集中式水源保護(hù)區(qū)。

2)本研究區(qū)處于覆蓋型巖溶區(qū)，下伏巖溶含水層富水性中等，巖溶較發(fā)育，溶孔、溶蝕裂隙發(fā)育相對均勻，水力聯(lián)系相對密切，可概化為均勻介質(zhì)含水層進(jìn)行模擬預(yù)測，為巖溶區(qū)的地下水污染預(yù)測提供一定的借鑒。