張 楠

(1.沈陽環(huán)境科學(xué)研究院，沈陽 110016；2.遼寧省城市生態(tài)重點(diǎn)實(shí)驗(yàn)室，沈陽 110167)

前 言

習(xí)近平總書記“兩山論”的提出標(biāo)志著我國生態(tài)文明建設(shè)進(jìn)入了新的階段。為貫徹落實(shí)習(xí)近平總書記關(guān)于地下水污染防治工作的重要批示精神，為保障地下水安全，讓人民群眾喝上放心水，生態(tài)環(huán)境部、自然資源部、住房和城鄉(xiāng)建設(shè)部、水利部、農(nóng)業(yè)農(nóng)村部在2019年3月聯(lián)合發(fā)布了《地下水污染防治實(shí)施方案》，該方案的發(fā)布標(biāo)志著我國地下水污染防治工作開啟了新的篇章。在2021年全國生態(tài)環(huán)境保護(hù)工作會(huì)議上的工作報(bào)告中，提出：“要推進(jìn)重點(diǎn)地區(qū)開展化工園區(qū)地下水環(huán)境狀況調(diào)查評(píng)估，持續(xù)開展地下水污染防治試點(diǎn)”。

遼寧省作為東北老工業(yè)基地，擁有省級(jí)以上工業(yè)集聚區(qū)114家，擁有化工園區(qū)21家，在這部分園區(qū)中不乏一些企業(yè)在生產(chǎn)過程中存在“跑、冒、滴、漏”等情況，這部分有毒有害物質(zhì)將會(huì)隨降水滲入到地下水中，造成園區(qū)地下水污染[1～3]。因此，為實(shí)現(xiàn)地下水資源的可持續(xù)利用，防止地下水受到污染，亟需加強(qiáng)對(duì)污染源頭的防范和治理[4-5]，同時(shí)，開展園區(qū)企業(yè)污染物滲漏模擬，摸清地下水中污染物遷移規(guī)律，對(duì)于生態(tài)環(huán)境部門制定監(jiān)督管理政策及方案具有重要意義[6-7]。

當(dāng)前國內(nèi)外地下水污染模擬軟件主要有3種，分別為Feflow、GMS和Vissual Modflow[8]。Feflow是有限元三角網(wǎng)格的計(jì)算軟件，相比其他兩種其主要優(yōu)勢(shì)在于可以任意剖分網(wǎng)格，并在局部區(qū)域加密，GMS和Vissual Modflow兩個(gè)軟件功能基本相同，均有Modflow、MT3D等經(jīng)典模塊，目前是應(yīng)用最普遍的兩個(gè)模型[9]。本文以沈陽市某高端裝備制造產(chǎn)業(yè)園為研究對(duì)象，采用GMS軟件對(duì)地下水中溶質(zhì)運(yùn)移特征進(jìn)行模擬，從而預(yù)測(cè)污染物的擴(kuò)散范圍，為地下水的污染防治提供依據(jù)。

1 研究區(qū)概況

研究區(qū)位于渾河沖擊平原，平均海拔高度36m，地勢(shì)平坦，屬于北溫帶半濕潤(rùn)季風(fēng)型大陸性氣候類型，多年平均降水量714mm，多集中在7、8月。區(qū)內(nèi)的地表水主要為渾河、細(xì)河兩條河流。

研究區(qū)含水構(gòu)造分布于渾河的高低漫灘地帶，地下水主要賦存于沖洪積的松散砂、砂礫石的孔隙之中。本區(qū)域受地貌和第四系沉積環(huán)境的控制，水文地質(zhì)條件呈現(xiàn)過渡性的特征。其上游含水層顆粒粗大，結(jié)構(gòu)簡(jiǎn)單；下游，含水層顆粒變細(xì)，結(jié)構(gòu)亦變得復(fù)雜。根據(jù)研究區(qū)內(nèi)的含水介質(zhì)、形成年代、水力特征和埋藏條件等，區(qū)內(nèi)的第四系含水層可分為潛水含水層(Q4)、淺層承壓水含水層(Q3)、深層承壓水含水層(Q1+2)。根據(jù)勘探資料，各層的埋藏條件，分布及變化規(guī)律、水位、水量、水質(zhì)有所區(qū)別，形成了三個(gè)相對(duì)獨(dú)立的含水層。包氣帶巖性為雜填土和粉質(zhì)粘土，平均厚度為7m，滲透系數(shù)2.21×10-5～5.78×10-5cm/s。

地下水的補(bǔ)給以垂向補(bǔ)給為主，側(cè)向補(bǔ)給為輔，主要接受降水入滲補(bǔ)給、稻田即渠道回滲補(bǔ)給及側(cè)向徑流補(bǔ)給；排泄方式主要有人工開采排泄、下游側(cè)向徑流排泄、河流排泄。

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

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

根據(jù)地下水的流向及等水位線并結(jié)合搜集已有各類水文地質(zhì)資料，東北部邊界沿等水位線概化為流量補(bǔ)給邊界；西南部邊界概化為流量排泄邊界；西北部邊界垂直于地下水等水位線，概化為隔水邊界即零流量邊界；東南部為渾河，渾河接受評(píng)價(jià)區(qū)地下水側(cè)向滲漏補(bǔ)給，地下水通過河床滲漏補(bǔ)給地表水，因此將東南部邊界概化為河流(River)邊界，河流河床滲透系數(shù)、河水位標(biāo)高及河床標(biāo)高數(shù)據(jù)均已在野外調(diào)查以及收集資料當(dāng)中獲取。各斷面流入、流出量，根據(jù)斷面處含水層滲透系數(shù)、斷面處水力坡度和斷面面積，由Darcy定律求出。

上部的潛水含水層和淺層承壓含水層之間有較厚的連續(xù)粘土層作為隔水層，水力聯(lián)系較弱，因此可將粘土層概化為隔水底板。由于潛水含水層較承壓含水層易于污染，是最敏感的含水層，因此第四系孔隙潛水含水層作為本次的主要研究含水層，該層厚度變化不大，一般在30m左右，地下水位埋深為5～15m。地下水運(yùn)動(dòng)以水平方式為主，自東北向西南方向徑流。

2.2 地下水流數(shù)學(xué)模型的建立

研究區(qū)地下水運(yùn)動(dòng)的數(shù)學(xué)模型可概化為非均質(zhì)各向同性三維地下水非穩(wěn)定流模型。

(1)

公式(1)中，H為模擬區(qū)內(nèi)含水層水頭(L)；K為模擬區(qū)內(nèi)含水層水平滲透系數(shù)(LT-1)；μ為模擬區(qū)內(nèi)含水層的給水度；z為模擬區(qū)內(nèi)含水層底板標(biāo)高(L)；ε為降雨和渠系灌溉入滲強(qiáng)度(LT-1)；E(x，y)為在淺層降水入滲區(qū)，其值為1，非降雨入滲區(qū)，其值為0；W為模擬區(qū)內(nèi)含水層其它源匯項(xiàng)(L3T-1)；Qj為模擬區(qū)內(nèi)含水層生產(chǎn)井抽水量(L3T-1)；H0為模擬區(qū)內(nèi)含水層初始水頭(L)；L1為第二類零流量邊界；L2為第二類給定流量邊界；L3為河流邊界；Ω為模擬區(qū)滲流計(jì)算區(qū)域；q為含水層邊界單寬流量(L2T-1)；n為各邊界面的外法線方向；Qr為流地下水交換量(m3/d)；Hr為河流水位標(biāo)高(m)。

2.3 數(shù)值模擬

2.3.1 網(wǎng)格剖分

模擬區(qū)面積約為158.16km2，對(duì)整個(gè)區(qū)域模型采用矩形網(wǎng)格剖分，網(wǎng)格間距為100m，剖分為250列、250行，共剖分矩形網(wǎng)格單元62 500個(gè)，其中有效單元16 336個(gè)，無效單元46 164個(gè)，計(jì)算節(jié)點(diǎn)位于單元中心。垂直方向上剖分為1層。

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

根據(jù)區(qū)內(nèi)的地形地貌、抽水試驗(yàn)資料以及水文地質(zhì)調(diào)查資料，將研究區(qū)含水層劃分為3個(gè)水文地質(zhì)參數(shù)分區(qū)，見圖1。

圖1 水文地質(zhì)參數(shù)分區(qū)圖Fig.1 Zoning map of hydrogeological parameters

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

本次模擬以2018年10月作為模型的初始流場(chǎng)(圖2)，2019年1月作為模型識(shí)別流場(chǎng)(圖3)，2019年7月作為驗(yàn)證流場(chǎng)(圖4)。應(yīng)力期以月為單位，共劃分為10個(gè)應(yīng)力期，每個(gè)應(yīng)力期又包括若干個(gè)時(shí)間步長(zhǎng)，時(shí)間步長(zhǎng)為模型自動(dòng)控制，嚴(yán)格控制每次的迭代誤差，在同一應(yīng)力期內(nèi)地下水補(bǔ)排項(xiàng)不變。

圖2 地下水初始流場(chǎng)圖Fig.2 Initial flow field map of groundwater

圖3 地下水識(shí)別期流場(chǎng)擬合圖Fig.3 Flow field fitting map of groundwater identification period

含水層實(shí)測(cè)流場(chǎng)和模擬流場(chǎng)的比較圖結(jié)果顯示，模擬流場(chǎng)與實(shí)測(cè)流場(chǎng)擬合較好，水位差小于0.02m的計(jì)算點(diǎn)達(dá)到80%以上，這表明模擬模型與實(shí)際地下水系統(tǒng)的水位在空間和時(shí)間上是比較吻合的，基本反映了地下水系統(tǒng)的水力特征，達(dá)到了模型精度要求，可以用于溶質(zhì)運(yùn)移預(yù)測(cè)。

圖4 地下水驗(yàn)證期流場(chǎng)擬合圖Fig.4 Flow field fitting map of groundwater verification period

識(shí)別驗(yàn)證后的水文地質(zhì)參數(shù)見表1。

表1 水文地質(zhì)參數(shù)分區(qū)表Tab.1 Hydrogeological parameter zoning table

3 溶質(zhì)運(yùn)移模型及預(yù)測(cè)結(jié)果分析

水質(zhì)模型是建立在水流模型基礎(chǔ)上，利用MT3DMS模塊進(jìn)一步來模擬預(yù)測(cè)地下水中污染質(zhì)的運(yùn)移情況。本次模擬不考慮污染物在含水層中的吸附、揮發(fā)、生物化學(xué)反應(yīng)，模型中各項(xiàng)參數(shù)予以保守性考慮。

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

地下水中溶質(zhì)運(yùn)移的數(shù)學(xué)模型可表示為：

(2)

(3)

式(2)中，C為模擬污染質(zhì)的濃度(mg/L)；ne為有效孔隙度；C′為模擬污染質(zhì)的源匯濃度(mg/L)；W為源匯單位面積上的通量；Vi為滲流速度(m/d)。

式(3)中，αijmn為含水層的彌散度；Vm，Vn為分別為m和n方向上的速度分量；|v|為速度模。

3.2 彌散度的確定

本次彌散度借鑒1996年研究區(qū)附近做的現(xiàn)場(chǎng)彌散試驗(yàn)資料[10]，彌散度確定為αL=0.19～0.47m，αT=0.09～0.23m，對(duì)比本區(qū)的具體水文地質(zhì)條件，彌散系數(shù)分區(qū)水文地質(zhì)參數(shù)分區(qū)一致，給定彌散系數(shù)值見表2。

表2 彌散系數(shù)分區(qū)表Tab.2 Partition table of dispersion coefficient

3.3 污染預(yù)測(cè)情景設(shè)定

對(duì)本裝備制造產(chǎn)業(yè)園區(qū)進(jìn)行地下水污染風(fēng)險(xiǎn)識(shí)別，得到WSCLC和HCBM是園區(qū)內(nèi)最大的風(fēng)險(xiǎn)源。本次以這兩個(gè)企業(yè)為例選擇兩種情況分別進(jìn)行預(yù)測(cè)分析。

3.3.1 連續(xù)源強(qiáng)污染物泄漏情景

情景設(shè)定：污水處理設(shè)施地面防滲層出現(xiàn)泄漏點(diǎn)導(dǎo)致防滲能力下降。

源強(qiáng)計(jì)算：假定按廢水處理量的0.01%滲漏，滲漏量按5%通過地表滲入地下水中，此情景下只對(duì)WSCLC進(jìn)行預(yù)測(cè)，WSCLC特征污染因子以氨氮為例，污染物濃度參照污水處理設(shè)施進(jìn)水水質(zhì)，計(jì)算得到氨氮的泄漏量為0.027kg/d(見表3)。

表3 各情景下污染物源強(qiáng)情況一覽表Tab.3 Pollutant source intensity under each scenario

3.3.2 瞬態(tài)源強(qiáng)污染物泄漏情景

情景設(shè)定：污水管道管線由于連接處(如法蘭、焊縫)開裂或腐蝕磨損等原因，會(huì)發(fā)生污染物泄漏。若恰好發(fā)生泄漏處的地下水防滲層斷裂或破壞，則將導(dǎo)致泄漏污染物污染地下水。

源強(qiáng)計(jì)算：設(shè)定采取的滲漏檢測(cè)發(fā)現(xiàn)及修復(fù)事故工況時(shí)間為15天；破裂泄漏孔半徑按10mm計(jì)，則污水的泄漏量為：3.14×0.012×1.5m/s×3600s/h×24h/d×15d=610.42m3。假設(shè)泄漏量全部通過地表進(jìn)入地下水，此情景下只對(duì)HCBM污水處理站進(jìn)行預(yù)測(cè)，其特征污染因子以鎳為例，污染物濃度參照相應(yīng)環(huán)評(píng)中的污染物產(chǎn)生濃度，則鎳泄漏進(jìn)入地下水中的污染物的量為：610.42m3×103L/m3×40mg/L÷1000000mg/kg=4.139kg/次(見表3)。

3.4 預(yù)測(cè)結(jié)果及分析

本次模擬預(yù)測(cè)根據(jù)污染風(fēng)險(xiǎn)分析的情景設(shè)計(jì)，在選定優(yōu)先控制污染物的基礎(chǔ)上，分別對(duì)地下水污染物在不同時(shí)段的運(yùn)移距離和超標(biāo)范圍進(jìn)行模擬預(yù)測(cè)。氨氮和鎳的超標(biāo)范圍均參照《地下水質(zhì)量標(biāo)準(zhǔn)》(GB/T14848-2017)中Ⅲ類水的要求，見表3。

3.4.1 連續(xù)源強(qiáng)下污水泄漏引起氨氮污染

1年、5年、10年、20年后氨氮濃度擴(kuò)散情況見圖5。氨氮在1年后最大運(yùn)移距離為208m，超標(biāo)范圍為56 821m2；5年后最大運(yùn)移距離為235m，超標(biāo)范圍為84 463m2；10年后最大運(yùn)移距離為416m，超標(biāo)范圍為131 352m2；20年后最大運(yùn)移距離為550m，超標(biāo)范圍為201 347m2。

圖5 氨氮濃度擴(kuò)散圖Fig.5 Diffusion of ammonia concentration

3.4.2 瞬態(tài)源強(qiáng)下污水泄漏引起的鎳污染

1年、5年、10年、20年后鎳濃度擴(kuò)散情況見圖6。鎳在1年后最大運(yùn)移距離為314m，超標(biāo)范圍為79 174m2；5年后最大運(yùn)移距離為425m，超標(biāo)范圍為138 617m2；10年后最大運(yùn)移距離為546m，超標(biāo)范圍為204 479m2；20年后最大運(yùn)移距離為710m，超標(biāo)范圍為234 393m2。

圖6 鎳濃度擴(kuò)散圖Fig.6 Diffusion of nickel concentration

4 結(jié)果與討論

本文以某裝備制造產(chǎn)業(yè)園內(nèi)的兩個(gè)企業(yè)為例，采用GMS軟件進(jìn)行數(shù)值模擬的方法，對(duì)兩種情況(連續(xù)源強(qiáng)和瞬態(tài)源強(qiáng))可能造成的地下水環(huán)境影響做出預(yù)測(cè)分析，得出以下結(jié)論。

4.1 連續(xù)源強(qiáng)污染物泄漏情景下，氨氮在1年后最大運(yùn)移距離為208m，超標(biāo)范圍為56 821m2；5年后最大運(yùn)移距離為235m，超標(biāo)范圍為84 463m2；10年后最大運(yùn)移距離為416m，超標(biāo)范圍為131 352m2；20年后最大運(yùn)移距離為550m，超標(biāo)范圍為201 347m2。污染物泄漏量持續(xù)增加，但只要地面防滲措施不出現(xiàn)問題，對(duì)地下水造成污染的風(fēng)險(xiǎn)較低，如果疊加出現(xiàn)防滲層破損情況，規(guī)劃區(qū)域內(nèi)的企業(yè)一旦發(fā)生污染物泄漏事故，隨時(shí)間延續(xù)，會(huì)對(duì)地下水下游方向的地下水水質(zhì)產(chǎn)生不同程度的污染影響，污染風(fēng)險(xiǎn)較大。

4.2 瞬態(tài)源強(qiáng)污染物泄漏情景下，鎳在1年后最大運(yùn)移距離為314m，超標(biāo)范圍為79 174m2；5年后最大運(yùn)移距離為425m，超標(biāo)范圍為138 617m2；10年后最大運(yùn)移距離為546m，超標(biāo)范圍為204 479m2；20年后最大運(yùn)移距離為710m，超標(biāo)范圍為234 393m2。污染物在一定時(shí)間內(nèi)對(duì)地下水環(huán)境造成不同程度的污染影響，隨著時(shí)間的增加和水動(dòng)力作用，污染物擴(kuò)散范圍雖然增大，但是濃度大幅降低。

5 結(jié) 論

5.1 地下水?dāng)?shù)值模擬方法可以直觀清晰地展現(xiàn)地下水中污染物的運(yùn)移情況，得出污染物對(duì)地下水環(huán)境敏感目標(biāo)的影響程度，進(jìn)而為污染類項(xiàng)目地下水配套相應(yīng)污染防治措施提供技術(shù)支持。

5.2 園區(qū)企業(yè)在重點(diǎn)區(qū)域采取嚴(yán)格的防滲措施可以有效降低污染物污染地下水環(huán)境，同時(shí)園區(qū)應(yīng)建立地下水水質(zhì)監(jiān)測(cè)網(wǎng)絡(luò)，實(shí)現(xiàn)地下水水質(zhì)實(shí)時(shí)監(jiān)測(cè)預(yù)警，從而最大程度降低企業(yè)對(duì)于地下水環(huán)境的影響。