胡艷春，馮園園，李婷婷，石 建

(1.山東省第一地質(zhì)礦產(chǎn)勘查院，山東 濟南 250010；2.山東省富鐵礦勘查技術(shù)開發(fā)工程實驗室，山東 濟南 250010；3.臨沂市自然資源和規(guī)劃局，山東 臨沂 276301；4.山東省國土測繪院，山東 濟南 250010)

1 研究區(qū)概況

魯中南地區(qū)是山東省煤炭、煤電、煤化工礦業(yè)經(jīng)濟區(qū)，分布有章丘煤田、淄博煤田、萊蕪煤田、新汶煤田等幾大煤田[1-3]。淄博煤田賦存于淄博向斜盆地石炭、二疊系地層中，含煤地層分布廣泛，面積約為418 km2。煤層分布差異性較強，厚度變化幅度大，共有煤層20余層，其中可采或局部可采的有12層。煤礦區(qū)處在淄博向斜盆地腹部，地質(zhì)構(gòu)造極為復(fù)雜，斷層縱橫切割，煤系地層下伏有奧陶系中統(tǒng)巖溶承壓含水層，共同構(gòu)成了區(qū)域復(fù)雜的水文地質(zhì)條件[4-5]。

淄博市洪山、寨里煤礦區(qū)各個礦自從上世紀五六十年代開始開發(fā)利用，因資源枯竭分別于1990年前后、2005年陸續(xù)閉坑。隨著時間推移，礦山開采形成的礦坑排水對礦區(qū)地下水資源及地質(zhì)環(huán)境產(chǎn)生的負面影響越來越突出，已經(jīng)威脅到當(dāng)?shù)毓┧踩鹆烁鹘绲母叨戎匾?。礦山閉坑后，老坑積水因泄流通道遭到堵塞，無法及時排泄，引起地下水位持續(xù)抬升，最終超過奧灰水，因受到水壓作用，老窯水可以大量串層，進而補給奧灰水，引起奧灰?guī)r溶含水層地下水污染[6-10]。

2 邊界條件的概化

2.1 含水層結(jié)構(gòu)

研究區(qū)巖溶水是具有系統(tǒng)補給、徑流、排泄條件的水文地質(zhì)單元[11]，雖然串層區(qū)范圍不大，考慮到巖溶水系統(tǒng)水力聯(lián)系的統(tǒng)一性，將模擬計算區(qū)擴展至地下水系統(tǒng)邊界，模擬區(qū)總面積約618.75 km2。

洪山、寨里煤礦區(qū)位于城東，淄博向斜盆地腹部，地勢東高西低。復(fù)雜的含水巖組巖性組成，形成了第四系松散巖類孔隙水、石炭—二疊系碎屑巖類裂隙水、以及煤層中的老窖水及煤系地層水、奧陶系碳酸鹽巖類裂隙巖溶水。第四系孔隙水富水性較差，對分布在其下部的碎屑巖裂隙含水層補給較少；礦坑水的主要充水來源是石炭-二疊系砂頁巖裂隙承壓水，礦坑含水巖組是研究區(qū)的最主要的污染源；煤層底板與奧陶灰?guī)r之間的夾層砂頁巖相對隔水巖組；當(dāng)?shù)刂饕揽抠x存在奧陶系灰?guī)r中的巖溶水作為生活飲用水河工業(yè)生產(chǎn)用水。礦山閉坑后，礦坑水通過底板破壞帶、天窗或廢棄破裂井管、斷裂帶進入深層承壓含水層，造成巖溶地下水水質(zhì)惡化。根據(jù)對區(qū)域地層巖性以及水文地質(zhì)條件的分析，可將該地區(qū)概化為4層結(jié)構(gòu)：煤系地層風(fēng)化砂巖層孔隙裂隙含水巖組和局部地段的松散巖類孔隙水為模型的第1層；煤系地層薄層灰?guī)r及砂巖裂隙含水巖組為第2層；開采煤層底板以下的煤系地層弱透水層(因采動作用引發(fā)越流天窗的形成)作為模型第3層，奧灰水作為第4層，模擬計算時，以上各層均概化為非均質(zhì)各向異性含水層(見圖1)。

圖1 研究區(qū)地質(zhì)概念模型及邊界條件概化

2.2 邊界條件

寨里-洪山閉坑礦山地下水串層污染區(qū)域并非一獨立水文地質(zhì)單元，為了降低側(cè)向邊界對數(shù)值模型的影響，盡量以自然邊界界定計算區(qū)。研究區(qū)域巖溶水南部邊界為大土屋侵入巖脈，為隔水邊界；北部至南定鎮(zhèn)賈莊村-店子一線，為弱透水流量邊界，店子-灃水-辛莊村一線為隔水邊界；東部邊界為地下水分水嶺，是零通量邊界；西部邊界為煤系地層深埋區(qū)域，徑流滯緩，按照隔水邊界處理(見圖1)。

寨里-洪山礦山閉坑后，石炭-二疊系砂頁巖裂隙承壓水與奧灰水通過砂頁巖相對隔水巖組層中的開采遺留通道(斷層構(gòu)造的揭露、采動裂隙、采空區(qū))發(fā)生水力聯(lián)系。該區(qū)地下水運動屬于非穩(wěn)定流，滲流計算區(qū)邊界條件概化為第二類邊界條件；以各含水層區(qū)域地下水觀測孔的統(tǒng)一觀測水頭值作為水流模擬初始值。

2.3 源匯項

大氣降水和地表水入滲補給第四系松散巖類孔隙水，主要以淺部徑流、人工開采方式排泄，研究區(qū)內(nèi)埋深較大，蒸發(fā)量可以忽略；砂頁巖裂隙水補給方式主要有兩種，一種是大氣降水入滲，另一種是局部地區(qū)存在的松散巖類孔隙水補給，排泄方式主要表現(xiàn)在煤礦排水、徑流以及人工開采；巖溶水補給源在東南山區(qū)，徑流向西北方向的過程中，遇到煤系地層后地下水遇阻富集，其排泄方式多樣，既有上升泉，又包括側(cè)向徑流，人類開采也是重要排泄途徑。地下水按照類型近乎獨立分布，分屬性質(zhì)各異的補徑排系統(tǒng)，天然條件下煤系地層裂隙水與巖溶水水力聯(lián)系差，煤礦采動后局部地段砂頁巖裂隙水與巖溶水存在互補關(guān)系。研究區(qū)大氣降水概化為面狀補給源，機民井、排水礦點概化為井點。

3 地下水流數(shù)值模型

3.1 模擬軟件選用

本次研究采用Visual MODFLOW平臺，先采用實際觀測資料進行模型識別和驗證，模型識別驗證完成再進行求解，預(yù)測分析串層污染對地下水量與水質(zhì)的影響。巖溶水污染物運移模擬采用MT3DMS模塊，充分模擬地下巖溶水系統(tǒng)中各污染物質(zhì)遷移運動規(guī)律，預(yù)測煤礦閉坑后，煤系地層中賦存的受到污染地下水向深部奧灰承壓含水層的遷移轉(zhuǎn)化過程[12-14]。

3.2 水流模型的建立

充分考慮淄博煤田地下水系統(tǒng)的實際條件，分析地下水補徑排過程以及邊界條件，根據(jù)達西定律以及滲流過程連續(xù)性方程，明確了淄博煤田洪山、寨里礦區(qū)水文地質(zhì)概念模型，以及三維非穩(wěn)定流數(shù)學(xué)模型：

(1)

式中：kxx，kyy，kzz分別代表沿 x，y，z 三個方向的滲透系數(shù)，單位：m/d；h是t 時刻任一點水頭值，單位：m；W 代表匯源項，單位：1/d；Ss為任一點處的儲水率，單位1/m；t為時間，單位d；Ω為滲流區(qū)域；S2為第二類邊界條件；nx、ny、nz分別為S2的外法線沿 x、y、z 軸方向的單位矢量；q為S2上單位面積的側(cè)向補給量(m3/d)。

3.2.1 網(wǎng)格剖分

研究區(qū)含水層可分為四層，平面上模擬 300×300 m 的差分網(wǎng)格對研究區(qū)展開剖分，每層剖分為139行和103列，邊界線以外做無效處理，第1～4層有效單元格數(shù)分別為 2 908、2 888、2 963、6 875個(見圖2)。

圖2 網(wǎng)格剖分示意圖 (左為第1～3層；右圖為第4層)

3.2.2 源匯項的確定

因研究區(qū)面積較大，水文地質(zhì)條件復(fù)雜，采用統(tǒng)一水文地質(zhì)參數(shù)計算誤差將十分巨大，因此根據(jù)研究區(qū)特征，將其劃分為四個分區(qū)。第I亞區(qū)，位于洪山礦區(qū)羅村鎮(zhèn)附近，是本次數(shù)值模擬的重點區(qū)域。亞區(qū)II，位于研究區(qū)北部，區(qū)內(nèi)有漫泗河斷層。亞區(qū)III，位于研究區(qū)中西部，地下水埋深大，也是地下水排泄區(qū)。亞區(qū) IV，位于東南西部，為地下水分水嶺，區(qū)內(nèi)巖溶地下水大部分由該區(qū)補給。年均降水補給量在研究區(qū)內(nèi)平均分布。

3.2.3 模型識別

利用間接法，采用1989-2002 年的歷史序列資料對模型參數(shù)進行調(diào)試，根據(jù)各觀測孔的模擬值和觀測值，反復(fù)地調(diào)整參數(shù)(包括滲透系數(shù)、儲水系數(shù)、給水度、孔隙度等)和調(diào)整某些源匯項進行擬合，從而實現(xiàn)對水文地質(zhì)參數(shù)和源匯項的反演。根據(jù)地質(zhì)、水文地質(zhì)條件等，對研究區(qū)初步分區(qū)賦值，運行程序，得出給定條件下的地下水位隨時間和空間的分布，并借助對同時期流場和長觀孔歷時曲線的擬合，對邊界值以及水文地質(zhì)參數(shù)等進行識別。經(jīng)過多次調(diào)整、調(diào)參擬合，得到1～4層含水層的水文地質(zhì)參數(shù)及分區(qū)圖(見圖3)。

圖3 地下水流動參數(shù)分區(qū)示意圖

表1 水文地質(zhì)參數(shù)分區(qū)表

3.2.4 模型檢驗

為驗證模型可靠性，利用2003年5月5日至2012年10月5日的實測資料進行驗證，劃分20個應(yīng)力期，計算周期10個水文年。經(jīng)運行，擬合程度較好，研究區(qū)各觀測結(jié)點計算水位與觀測水位過程線趨勢相同，水位擬合絕對誤差較小，模擬流場與實際流場基本擬合，總體流動方向相同，可利用該模型進行模擬與預(yù)測(圖4和圖5)。

圖4 2003-2012年110號 圖5 2012年巖溶水平面孔水位擬合過程線 流場擬合圖

4 溶質(zhì)運移模型

4.1 水動力彌散方程的數(shù)學(xué)模型

本次以水流模型為基礎(chǔ)，進一步豐富建立溶質(zhì)運移模型，模型概化完全繼承自水流模型，流體統(tǒng)一概化為均值流體，且該流體不可壓縮，其粘度以及密度保持不變。

閉坑礦山老窯水相當(dāng)于一個污染源，即石炭—二疊系砂頁巖裂隙承壓水中的溶質(zhì)污染物通過水力聯(lián)系區(qū)不斷進入奧灰水含水層。這一過程可概化為溶質(zhì)在飽和模型中的運移，其溶質(zhì)運移方程為：

式中：Ck為溶質(zhì)k在地下水中的濃度，(m/L3)；n為孔隙度；Dij為水力彌散系數(shù)張量，(L3/s)；∑Rn為反應(yīng)項，[m/(L3·s)]；Cks為物質(zhì) k 源匯項的濃度，(m/L3)；qs為含水層中流體單位體積的過流量，(1/s)；vi為滲流速率，(L/s)。Co為源匯項的濃度，(m/L3)；Γ1研究區(qū)邊界；C1(x，y，z)為第一類邊界上的已知濃度，(m/L3)。

4.2 彌散度的確定

水質(zhì)模型中涉及的參數(shù)包括彌散度、有效孔隙度等。參數(shù)根據(jù)各觀測孔的模擬值和觀測值，結(jié)合實際水文地質(zhì)條件，經(jīng)過反復(fù)調(diào)試擬合，由模型反演得到(見圖6)。

圖6 彌散參數(shù)分區(qū)示意圖

模型反演結(jié)果表明，煤系地層彌散度取值范圍為5～10.3 m，奧灰含水層取值范圍為1.24～2.4 m，串層點彌散度最大取值30 m。根據(jù)與研究區(qū)具有相似水文地質(zhì)條件的堠皋地區(qū)和大武水源地勘查結(jié)果，奧陶系碳酸鹽類裂隙巖溶含水巖組彌散度取值范圍為1.4～3.5 m[15-18]，補給區(qū)地下水流更新快，巖溶發(fā)育程度相對較好，水力坡度相對較大，水流速度相對較快，彌散度相對較大，如位于補給區(qū)的湖田地區(qū)彌散度為2.6 m。本次模擬反演結(jié)果基本呈現(xiàn)該規(guī)律，模擬反演結(jié)果與臨區(qū)數(shù)值基本一致。

表2 溶質(zhì)運移參數(shù)表

4.3 參數(shù)識別與驗證

選取1989-5-5日至2002-9-30日作為模擬識別時段，每年分為豐、枯兩個時段，共27個應(yīng)力期，每個應(yīng)力期分為10個時間步長。利用1989-2002年長期觀測數(shù)據(jù)，可以建立研究區(qū)地下水污染物運移模型，并通過MT3D 軟件運算求解，對比研究區(qū)內(nèi)現(xiàn)有監(jiān)測井2002年的濃度值。受到水質(zhì)資料限制，濃度邊界條件及初始濃度刻畫不夠清晰。模擬水質(zhì)與監(jiān)測水質(zhì)變化趨勢基本一致(見圖7)，過程線在一定程度上能夠反映污染物的遷移轉(zhuǎn)化規(guī)律。

圖7 1989-2002年寨里鄉(xiāng)觀測孔TDS以及 SO42-擬合曲線(計算值，實測值)

4.4 模擬結(jié)果分析

受礦區(qū)串層影響，礦坑老窖水隨水流進入串層通道，對奧灰水造成污染。礦區(qū)串層污染的典型污染物是SO42-離子，因此以SO42-離子為例進行污染預(yù)測。根據(jù)本次工作形成的數(shù)值模型，在不采取應(yīng)對措施的前提下，以2013年污染物濃度場作為模型初始濃度條件，對研究區(qū)進行10、20、25、30年模擬預(yù)測(見圖8)。

圖8 未來30年寨里監(jiān)測點SO42-濃度歷時曲線

受串層因素影響，寨里村監(jiān)測點2 720 d之前SO42-濃度迅速增加，由70.2 mg/L增加至96.76 mg/L；隨后水頭差逐漸減小，串層補給作用強度減弱，濃度增幅逐漸變緩，9 218 d以后奧灰水濃度升高趨勢減緩，此時水質(zhì)污染范圍擴大至下游10.55 km處。礦坑老窖水在漫長的時間里，持續(xù)影響奧灰水水質(zhì)(見圖9)。

5 污染防治措施

針對研究區(qū)實地情況，采取注漿封堵的方案可以阻斷地下水污染路徑。在已建立的數(shù)值模型的基礎(chǔ)上，將污染串層部位作注漿區(qū)，設(shè)定滲透系數(shù)1×10-6cm/s。以SO42-為指示因子進行污染治理效果預(yù)測。

自模擬初期2013年起算，采取治理措施10年后，預(yù)測洪山礦區(qū)奧灰水SO42-超標范圍依然約為0.225 km2，與不進行治理奧灰水SO42-預(yù)測超標范圍2.25 km2相差一個數(shù)量級(見圖10)。

圖10 洪山礦區(qū)2023年預(yù)測SO42-濃度平面分布對比圖(左側(cè)為不治理,右側(cè)為治理)

模擬結(jié)果顯示，采取治理措施30年后，全區(qū)SO42-濃度均降低到250 mg/L以下，即煤礦串層引起的水質(zhì)污染得到逐步緩解，全區(qū)域巖溶水中SO42-離子濃度基本達到飲用水水質(zhì)標準，硬度、TDS等參數(shù)模擬結(jié)果也具有相同趨勢(見圖11)。

圖11 研究區(qū)30年后奧灰水SO42-濃度平面等值線(左側(cè)為不治理,右側(cè)為治理后)

6 結(jié)語

(1)礦坑老窖水污染奧灰水的持續(xù)時間漫長，如果不進行治理，其污染擴散范圍將繼續(xù)擴大，2020年之前SO42-濃度增加幅度較大，隨后水位差逐漸變小，串層補給作用減小，后期濃度增幅逐漸變緩，2038年以后，奧灰水濃度升高趨勢減緩，奧灰水水質(zhì)污染范圍將擴大至下游10.55 km處。

(2)對串層污染區(qū)采取注漿封堵的修復(fù)工程措施，在注漿層滲透性1×10-6cm/s方案下，在一定程度上可控制該區(qū)的礦坑積水、老窖水對奧灰水的影響，使奧灰水水質(zhì)得到改善。治理30 a后，煤礦串層引起的水質(zhì)污染可逐步緩解，全區(qū)域巖溶水中SO42-等組分含量達到飲用水水質(zhì)標準。