基于VM-MT3D的污染物質(zhì)運移數(shù)值模擬預(yù)測分析

王成文

(寧夏回族自治區(qū)水文環(huán)境地質(zhì)勘察院，寧夏 銀川 750011)

地下水作為寶貴的供水資源和重要的生態(tài)環(huán)境因子，在保證城市生活用水、工業(yè)農(nóng)業(yè)生產(chǎn)用水和生態(tài)環(huán)境平衡等方面扮演著不可替代的角色，同時也是保持“大氣降水—地表水—地下水”三水轉(zhuǎn)化、循環(huán)的主要參與者，對于保障社會經(jīng)濟發(fā)展和安全供水具有十分重要的作用。近些年，隨著城鎮(zhèn)一體化的快速發(fā)展， 地下水開發(fā)利用過程中存在著水資源供需矛盾，地下水過度被開采、地下水水位下降和地下水污染等多種問題[1-3]。

Visual Modflow是地下水水流和溶質(zhì)運移模擬軟件，廣泛應(yīng)用于地下水流場和污染溶質(zhì)運移分析中。借助Visual Modflow 軟件， 趙明明等[4]分析了玉井煤礦對地下淺含水層會產(chǎn)生較大影響，影響范圍達到了 5 000 m。張韻等[5]建立了某經(jīng)濟開發(fā)區(qū)地下水流和溶質(zhì)運移的可視化模型，并預(yù)測出Cr6+10年后運移112 m。趙娜[6]在對某倉庫進行了前期水文地質(zhì)條件的了解和資料的準(zhǔn)備，建立相應(yīng)的 GMS 數(shù)值模型，對鉛酸蓄電池中的硫酸根離子的運移特征進行了預(yù)測研究，并提出適合研究區(qū)當(dāng)?shù)氐奈廴痉雷o措施。李丹等[7]研究對預(yù)測地下水污染發(fā)展趨勢和地下水溶質(zhì)遷移機制，選擇最佳防治措施，制定相關(guān)地下水保護政策，具有重要的實際意義。

本次研究借鑒于地區(qū)開展供水水文地質(zhì)勘察資料，結(jié)合研究區(qū)的實際特征，利用Visual Modflow軟件MT3DMS模型(簡稱VM-MT3D)[8-10]進行污染物質(zhì)運移模擬預(yù)測，建立穩(wěn)定流數(shù)學(xué)模型，采用地下水流動模型選取氨氮(NH3-N)為特征污染物，預(yù)測NH3-N的運移對土壤及地下水造成的影響，為地下水資源的安全和管理提供依據(jù)。

1 研究區(qū)地下水流數(shù)值模型

1.1 研究區(qū)概況

研究區(qū)為某生活垃圾焚燒發(fā)電三期擴建廠區(qū)，所在的微水文地質(zhì)單元邊界外擴形成的面積約為20 km2的矩形區(qū)域。研究區(qū)所在地形為低山丘陵地形，含水層類型主要為火山巖類孔洞裂隙含水層，地下水的徑流條件受地形地勢的影響明顯，地下水徑流方向和地形一致，區(qū)域地下水是受地勢相對較高的地形影響，形成局部的地下水流動系統(tǒng)，整體由南部低山丘陵向西北流動，研究區(qū)基本情況見圖1。

圖1 研究區(qū)基本狀況

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

1.2.1水文地質(zhì)條件

在研究區(qū)域內(nèi)廣泛分布含水層為火山巖類孔洞裂隙含水層，均為潛水含水層，巖性主要是火山碎屑巖石、玄武巖、蜂窩狀-氣孔狀玄武巖為主，含水層直接出露地表，屬于裸露行火山巖孔洞裂隙水，研究區(qū)水文地質(zhì)剖面見圖2。

地下水主要來源為大氣降水入滲補給和側(cè)向補給，還有少量灌溉回滲補給；主要排泄方式為地下水的側(cè)向徑流排泄，其次為蒸發(fā)。研究區(qū)域周邊地下水開采利用程度較低，地下水的徑流條件受地形地勢的影響明顯，多以三維流的形式徑流，由地勢相對較高的顏春嶺低山丘陵向四周徑流，地下水徑流方向和地形一致，但水力梯度小于地形坡度，水力梯度1‰～6‰。地下水總體流動方向為南東流向北西，研究區(qū)水位等值線見圖3。

因此，本次水文地質(zhì)模型概化僅將火山巖、玄武巖類孔洞裂隙含水層概化為一層，作為模擬預(yù)測層位，厚度一般在10～75 m，其中火山巖類地下水含水層厚度為1～30 m，玄武巖含水巖層厚度30～70 m。所在區(qū)水位等值線見圖3。

圖2 研究區(qū)水文地質(zhì)剖面

圖3 研究區(qū)水位等值線

圖4 研究區(qū)模擬邊界條件

1.2.2邊界概化

結(jié)合研究區(qū)水文地質(zhì)特征，由于東側(cè)邊界和西側(cè)邊界垂直于等水位線，設(shè)置為零通量邊界；南側(cè)邊界為地下水的補給邊界，北側(cè)邊界為地下水的排泄邊界；在垂向上，潛水含水層自由水面作為模型上邊界，通過該邊界潛水與系統(tǒng)外發(fā)生垂向上的水量交換，如大氣降水入滲補給、蒸發(fā)排泄，模擬區(qū)邊界條件見圖4。

1.3 地下水?dāng)?shù)學(xué)模型

a)水流模型。根據(jù)水文地質(zhì)概念模型，地下水流概化為非均質(zhì)水平各向同性、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)，可用如下方程的定解問題來描述[11-12]：

b)溶質(zhì)模型。三維溶質(zhì)運移可用以下微分方程的定解問題描述[12-14]：

x、y、z∈Ω,t≥0

(2)

污染物在地下水中的運移轉(zhuǎn)化過程是極其復(fù)雜的[15]，本次模擬研究不考慮污染物運移過程中的吸附、化學(xué)反應(yīng)和生物降解等作用，只考慮對流彌散作用對污染物運移的影響，彌散系數(shù)與孔隙的平均流速呈線性關(guān)系，其比值為彌散度，在模型中流速是自動計算的，溶質(zhì)運移模型需要給定縱向彌散度，橫向彌散度為縱向彌散度的1/10。本次評價縱向彌散度根據(jù)前人的研究成果和一些類似水文地質(zhì)條件的模擬結(jié)果確定[16]：縱向彌散度取10 m，橫向彌散度取1 m。溶質(zhì)運移參數(shù)見表1。

表1 NH3-N溶質(zhì)運移模型參數(shù)

1.4 數(shù)值模擬相關(guān)參數(shù)

1.4.1水文地質(zhì)參數(shù)

研究區(qū)含水層主要以氣孔狀、火山碎屑巖石為主，水文地質(zhì)條件較為均勻，滲透系數(shù)空間變化較小，將滲透系數(shù)和給水度等含水層參數(shù)分區(qū)值以面狀形式導(dǎo)入模型中，滲透系數(shù)為0.09～0.14 m/d，給水度為0.05～0.15，水文地質(zhì)參數(shù)見表2。

表2 不同歷時時間NH3-N對地下水的影響情況

1.4.2源匯項概算

研究區(qū)源匯項比較簡單，地下水的主要補給來源為大氣降水入滲補給，主要排泄途徑為蒸發(fā)排泄、少量人工開采及地下水側(cè)向徑流流出。

a)補給項。模擬區(qū)源匯項比較簡單，補給主要為大氣降水補給，大氣降水量取多年平均降水量1 700 mm，當(dāng)?shù)亟涤耆霛B系數(shù)0.09～0.13，因此，研究區(qū)大氣降水入滲補給地下水量根據(jù)式計算[17]可得，大氣降水入滲速率為0.000 42～0.000 61 m/d。

V=α·P·10-3/365

(3)

式中V——降雨入滲速率，m/d；α——降雨入滲系數(shù)；P——降雨量，mm/a。

b)排泄項。地下水排泄主要為側(cè)向流出和蒸發(fā)。研究區(qū)及周邊區(qū)域地下水水位埋藏深度較淺，地下水水位一般在1.93～7.74 m，當(dāng)?shù)卣舭l(fā)極限深度為3.5 m，因此在水位埋藏深度較淺的區(qū)域，地下水會以蒸發(fā)的方式排泄，蒸發(fā)是項目所在區(qū)域地下水的排泄方式之一。蒸發(fā)量取當(dāng)?shù)囟嗄昶骄舭l(fā)量1 524.3 mm，蒸發(fā)極限深度取3.5 m；側(cè)向流出量根據(jù)達西定律計算得到，地下水年變幅源匯項見表3。

表3 地下水源匯項概算變化量 單位：億m3/a

1.5 模型識別與檢驗

通過模型的識別和驗證，使模型達到所需的精度的情況下進行模擬預(yù)測。通過模擬的水位等值線和水文地質(zhì)勘查過程中6個水文地質(zhì)鉆孔的實測水位和模型的模擬水位進行模型的識別驗證。研究區(qū)域地下水基本處于無開采狀態(tài)，地下水水位年際動態(tài)變化幅度較小，水位較穩(wěn)定。因此，本次模擬采用非穩(wěn)定流模型進行水位的識別，模型地下水流場見圖5，可知模擬區(qū)水位等值線和實際水位等值線相比，水位值變化不大，水位等值線的形狀基本一致，水力梯度基本不發(fā)生變化，污染物質(zhì)隨地下水流向徑流，研究區(qū)水位可以反映出真實地下水中NH3-N濃度隨水流的歷時變化過程。

圖5 模擬地下水流場滲流方向

在水位點的識別驗證中，研究區(qū)6個水位地質(zhì)鉆孔的模擬預(yù)測水位和實際水位情況見表4，水位點模擬水位和真實水位的擬合情況見圖6。

通過分析可知，6個水文地質(zhì)鉆孔模擬水位水位和實際水位差均在0～1.2 m，水位誤差較小，滿足模擬預(yù)測所需精度要求，可以用做地下水溶質(zhì)運移模擬預(yù)測。

表4 研究區(qū)觀測點模擬水位與實際水位對比

圖6 水位觀測點擬合

2 模型預(yù)測結(jié)果與分析

2.1 模擬預(yù)測設(shè)定

根據(jù)擴建工程滲漏事故發(fā)生后污染羽的變化情況和服務(wù)年限劃分模擬預(yù)測周期，利用Visual Modflow軟件MT3D模型進行數(shù)值模擬。在網(wǎng)格剖分中，大區(qū)域上將模擬區(qū)剖分為10 m×10 m的網(wǎng)格，為了詳細刻畫項目區(qū)地下水流動特征、模擬預(yù)測事故工況下污染物對區(qū)域地下水環(huán)境的影響情況，模型共分為415行×411列，其中有效單元格共119 414個，網(wǎng)格剖分見圖4。

根據(jù)研究區(qū)地下水流場特征、含水層的分布情況和水文地質(zhì)單元的邊界情況，確定本次模擬預(yù)測的范圍為包含擴建三期工程單元及下游區(qū)域的面積約12 km2的區(qū)域，模擬預(yù)測單元為焚燒廠垃圾池滲液處理站，選擇污染物質(zhì)因子為NH3-N，設(shè)定污染源強濃度為20 000 mg/L，模擬預(yù)測的情景為垃圾池的垃圾滲濾液連續(xù)滲漏90天，每天滲漏1.0 m3的情況，該情景設(shè)計是基于滲濾液處理站每天處理的垃圾能力和容載能力而定，并參照曾經(jīng)滲漏檢測的污染物質(zhì)濃度值。

2.2 污染物運移預(yù)測結(jié)果分析

在模型識別驗證的基礎(chǔ)上，將污染源輸入模型，預(yù)測時間分別為滲漏事故后100天、1 000天、10年、20年、25年(擴建工程服務(wù)期滿)和30年(服務(wù)期滿后5年)共6個時段，模擬NH3-N污染羽的分布情況，滲漏事故發(fā)生后NH3-N污染羽變化情況見圖7。

a)100天

b)1 000天

c)10年

d)20年

e)25年

f)30年

結(jié)合表5可知，滲漏事故發(fā)生后，受滲漏的高濃度的滲濾液的影響NH3-N污染羽影響距離隨時間不斷增長，污染羽范圍隨時間逐漸增大，滲漏事故發(fā)生30年后，NH3-N最遠影響距離為496.2 m；污染羽范圍為36 800 m2；污染羽范圍內(nèi)NH3-N濃度隨時間逐漸減小，滲漏事故發(fā)生30年后，研究區(qū)NH3-N最大濃度為2.19 mg/L，高于GB/T 14848—2017《地下水質(zhì)量標(biāo)準(zhǔn)》[18]中NH3-N的Ⅲ類標(biāo)準(zhǔn)限值。

表5 不同歷時點NH3-N對地下水的影響情況

滲濾液處理站滲漏事故發(fā)生處和滲濾液處理站下游廠界處NH3-N濃度歷時曲線見圖8、9，可知滲漏事故發(fā)生后，滲濾液處理站滲漏處NH3-N隨時間逐漸減小，滲漏事故發(fā)生7.4年后NH3-N降至地下水Ⅲ類NH3-N 0.5 mg/L標(biāo)準(zhǔn)限值；滲濾液處理站下游廠界處NH3-N隨時間先增大后逐漸減小，滲濾液滲漏事故發(fā)生后1.1～14.8年，NH3-N的污染濃度高于地下水質(zhì)量標(biāo)準(zhǔn)限值，滲漏事故發(fā)生30年后，滲濾液處理站下游廠界處NH3-N降為0.002 3 mg/L。

圖8 處理站滲漏處地下水中NH3-N濃度歷時曲線

圖9 處理站下游廠界處地下水中NH3-N濃度歷時曲線

3 主要結(jié)論

a)利用Visual Modflow軟件MT3D模型進行地下水NH3-N溶質(zhì)運移數(shù)值模擬，研究區(qū)滲漏處理站的NH3-N污染羽影響距離隨時間不斷增長，范圍不斷增大；模擬預(yù)測滲漏事故發(fā)生30年后，NH3-N最大濃度為2.19 mg/L，最遠影響距離為496.2 m，污染羽范圍為36 800 m2。

b)滲濾液處理站滲漏處NH3-N隨時間逐漸減小，滲漏事故發(fā)生7.4年后NH3-N降至地下水Ⅲ類標(biāo)準(zhǔn)；處理站下游廠界處NH3-N隨時間先增大后逐漸減小，滲濾液滲漏事故發(fā)生后1.1～14.8年，NH3-N的污染濃度高于地下水質(zhì)量標(biāo)準(zhǔn)限值，直到30年后下游廠界處NH3-N降為0.002 3 mg/L。

c)通過對地下水水位的驗證識別，發(fā)現(xiàn)此次地下水溶質(zhì)模擬預(yù)測可靠，可以為更復(fù)雜的地下水溶質(zhì)運移研究提供一定參考。