基于Visual Modflow 軟件的地下水污染物運(yùn)移模擬研究

陳 峰，吳 濤

（中國(guó)煤炭地質(zhì)總局 水文地質(zhì)工程地質(zhì)環(huán)境地質(zhì)勘查院，河北 邯鄲 056004）

0 引 言

水是人類(lèi)生存不能缺少的資源，隨著國(guó)民經(jīng)濟(jì)的迅速發(fā)展和人民生活水平的不斷提高，人們對(duì)水資源的需求量也不斷增長(zhǎng)。地下水是工業(yè)、農(nóng)業(yè)和人類(lèi)等生活的重要部分，根據(jù)調(diào)查，我國(guó)大中城市中地下水資源均受到一定的污染。

近年來(lái)，人們對(duì)環(huán)境保護(hù)的意識(shí)不斷增強(qiáng)，地下水環(huán)境質(zhì)量也越來(lái)越受到重視。對(duì)一些可能會(huì)造成地下水污染的企業(yè)，防范措施要求也逐漸提高，但仍不能排除特殊情況下污染物泄露對(duì)地下水造成污染。應(yīng)用科學(xué)有效的計(jì)算方法模擬污染物泄露的影響范圍及程度，對(duì)于地下水污染的防治有著極其重要的意義。目前地下水系統(tǒng)數(shù)值模擬方法主要有有限差分法（FDM）、有限單元法（FEM）、邊界元法（BEM） 和有限分析法（FAM） 等[1]，可以用來(lái)模擬地下水位下降、地下水位過(guò)量開(kāi)采和地下水污染等問(wèn)題。隨著地下水科學(xué)與計(jì)算機(jī)技術(shù)的快速發(fā)展，地下水進(jìn)行水流和水質(zhì)運(yùn)移采用數(shù)值模擬進(jìn)行預(yù)測(cè)研究也得到了廣泛的應(yīng)用[2]。此次利用Visual Modflow軟件來(lái)模擬研究區(qū)地下水污染物運(yùn)移情況，為研究區(qū)地下水污染物的防治提供科學(xué)依據(jù)[3]。

1 Visual Modflow 軟件

Visual Modflow 軟件是加拿大Waterloo 水文地質(zhì)公司在美國(guó)地質(zhì)勘探局原MODFLOW 軟件基礎(chǔ)上應(yīng)用現(xiàn)代可視化技術(shù)研發(fā)成的，該軟件是能夠建立三維地下水流動(dòng)和模擬污染物遷移模型的軟件，于1994 年8 月在國(guó)際上公開(kāi)發(fā)行[12-13]，截止目前，該軟件的三維地下水流模擬在國(guó)際上也得到認(rèn)可。地下水一旦遭到破壞，由于地下水中的污染物降解和遷移比較緩慢，地下水一旦污染水化學(xué)就很難恢復(fù)，采用Visual Modflow 軟件在地下水中建立污染物運(yùn)移模型，來(lái)模擬地下水的運(yùn)移過(guò)程。

2 研究區(qū)概況

2.1 研究區(qū)特征

研究區(qū)為河北省館陶縣工業(yè)園區(qū)，地貌特征屬平原型，地勢(shì)比較平坦，區(qū)域位于沖湖積平原區(qū)水文地質(zhì)單元。研究區(qū)地下水賦存于第四系松散巖類(lèi)孔隙含水層中，含水層巖性主要有細(xì)砂、粉砂、粉細(xì)砂。研究區(qū)范圍內(nèi)，淺層水底板一般埋深在40 ～55 m，含水層單位涌水量一般1.2 ～2.7 L/s·m，滲透系數(shù)5.12 m/d，為潛水。

該區(qū)地下水補(bǔ)給主要有大氣降水的入滲補(bǔ)給、地下水側(cè)向徑流補(bǔ)給。地下水由西南流向東北，水力坡度在1/1 000 左右。地下水的排泄方式主要為地下水徑流排泄和人工開(kāi)采等。

2.2 地下含水層劃分及特征

研究區(qū)的地下水主要賦存于第四系多層結(jié)構(gòu)的松散巖層中，以第四系地層為基礎(chǔ)，以水文地質(zhì)條件為依據(jù)，自上而下將第四系劃分4 個(gè)含水組，即Ⅰ、Ⅱ、Ⅲ、Ⅳ含水組。第Ⅰ含水組大致相當(dāng)于全新統(tǒng)（Q）4、第Ⅱ含水組大致相當(dāng)于上更新統(tǒng)（Q）3，第Ⅲ含水組大致相當(dāng)于中更新統(tǒng)（Q）2，第Ⅳ含水組大致相當(dāng)于下更新統(tǒng)（Q）1。

（1） 第Ⅰ含水組。

含水組底板埋深50 ～110 m，為一套沖積洪積—沖積湖積沉積物。含水層巖性以細(xì)砂為主，含水層厚度26 ～60 m，單層厚度一般5 ～l0 m，單位涌水量一般為5 ～17 m3/h。含水組上段為潛水，下段由微承壓水變成承壓水，水位埋深26.1 ～29.4 m。館陶北縣東北局部地區(qū)下部為咸水體外，其余地區(qū)礦化度0.9 ～2.0 g/L，為淡水。

（2） 第Ⅱ含水組。

含水組底界埋深200 ～230 m，為一套沖積湖積沉積物。巖性主要為細(xì)砂和粉細(xì)砂，東南部為中砂含礫石，顆粒稍粗且較厚，含水層厚度為28.2 ～78 m，單層厚度最大可達(dá)到20 m。為承壓水，單位涌水量為6 ～l0 m3/h·m，水位埋深為29.5 ～30.1 m。主要為咸水，礦化度一般大于3 g/L。

（3） 第Ⅲ含水組。

含水組底板埋深300 ～320 m，為沖積湖積沉積物，巖性主要為細(xì)砂，厚度20 ～50 m，北薄南厚，為承壓水。基本全為淡水，礦化度小于2 g/L。

（4） 第Ⅳ含水組。

含水組巖性沖積湖積和冰水沉積物，為承壓水，巖性主要為細(xì)中砂，單位涌水量一般5 m3/h·m左右。全為淡水，礦化度小于1 g/L。

2.3 地下含動(dòng)態(tài)特征

從淺層淡水長(zhǎng)期觀測(cè)資料分析，水位變化屬降水入滲—人工開(kāi)采型，年內(nèi)水位具有明顯的季節(jié)性變化，影響水位變化的主要因素是降水量大小，每年l ～3 月份為水位平緩期，是年內(nèi)最高水位；3 ～6 月中旬為農(nóng)業(yè)開(kāi)采期，且降雨量少，水位一般下降最大，6 月中旬為年內(nèi)最低水位，7 ～9 月份為降雨期，水位逐漸回升，至下一年1 月份達(dá)到年最高水位，年變幅0.2 ～3.4 m。深層淡水較淺層淡水具有滯后效應(yīng)，滯后期l ～2 個(gè)月，年變幅2 ～4 m。

3 地下水污染物運(yùn)移模型建立

3.1 地下水流數(shù)值模型的建立

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

水文地質(zhì)概念模型是對(duì)水文地質(zhì)條件的簡(jiǎn)化，對(duì)地下水系統(tǒng)的科學(xué)概化，是把含水體的邊界性質(zhì)、地層內(nèi)部的結(jié)構(gòu)、地下水的滲透性能、地下水的水力特征、地下水的補(bǔ)給情況和排泄情況等條件概化，為了更方便建立數(shù)學(xué)和物理模擬的基本模式[4-7]。建立研究區(qū)的水文地質(zhì)概念模型是進(jìn)行預(yù)測(cè)評(píng)價(jià)的第一步。

首先對(duì)研究區(qū)水文地質(zhì)條件進(jìn)行概化處理，研究區(qū)內(nèi)水文地質(zhì)條件簡(jiǎn)單，為潛水含水層。由于研究區(qū)范圍只是水文地質(zhì)單元的一部分，模型邊界為人為劃定邊界，研究區(qū)西南邊界和東北邊界以地下水等水位線為界，按一類(lèi)邊界處理；西北邊界、東南邊界基本與地下水等水位線垂直，按隔水邊界處理；在垂向上，淺層含水層自由水面為系統(tǒng)的上邊界，通過(guò)該邊界，淺層水與系統(tǒng)外發(fā)生垂向交換，比如大氣降水的入滲補(bǔ)給；以含水層底板作為模型的下邊界。

3.1.2 數(shù)學(xué)模型的建立

通過(guò)對(duì)水文地質(zhì)概念模型的分析，依據(jù)滲流連續(xù)性方程和達(dá)西定律，建立研究區(qū)與地下水系統(tǒng)水文地質(zhì)概念模型相對(duì)應(yīng)的二維穩(wěn)定流數(shù)學(xué)模型[8-11]：

式中：H為地下水水頭，m；K為滲透系數(shù)，m/d；H0（x、y） 為開(kāi)始地下水水頭函數(shù)，m；H1（x、y）為第一類(lèi)邊界地下水水頭函數(shù)，m；q（x、y） 為含水層二類(lèi)邊界單位面積過(guò)水?dāng)嗝嫜a(bǔ)給流量函數(shù)，m2/d；ε 為源匯項(xiàng)強(qiáng)度（包括開(kāi)采強(qiáng)度等），m/d；Ω 為滲流區(qū)域；B1為水頭已知邊界，第一類(lèi)邊界；B2為流量已知邊界，第二類(lèi)邊界；n為滲流區(qū)邊界的單位外法線方向。

3.1.3 網(wǎng)格剖分

對(duì)研究區(qū)進(jìn)行網(wǎng)絡(luò)剖分，網(wǎng)格剖分在水平方向上用正交網(wǎng)格剖分法，將模型剖分成80×80 的單元格，在重點(diǎn)研究區(qū)域進(jìn)行了細(xì)化剖分，在原剖分基礎(chǔ)上進(jìn)行了2 倍加密，如圖1 所示。

圖1 模型剖分及邊界條件示意Fig.1 Model division and boundary conditions

3.1.4 污染物質(zhì)的確定

根據(jù)現(xiàn)場(chǎng)調(diào)查分析得出，研究區(qū)的污染不是一個(gè)明顯的污染源，而是有存在多個(gè)次要或者一個(gè)主要污染源，分析得出評(píng)價(jià)區(qū)內(nèi)污染物COD 負(fù)荷比73.56%、氨氮負(fù)荷比26.44%，COD 為主要污染物。此次利用Visual Modflow 軟件對(duì)這兩種污染物進(jìn)行數(shù)字模擬。

3.1.5 水文地質(zhì)參數(shù)的選取

為了更加準(zhǔn)確地表達(dá)出該區(qū)水文地質(zhì)情況，模型中參數(shù)的確定主要依據(jù)研究區(qū)已有的水文地質(zhì)勘查成果，并結(jié)合研究區(qū)的區(qū)域水文地質(zhì)資料和多年的經(jīng)驗(yàn)值，確定了含水層參數(shù)值，見(jiàn)表1。

表1 水文地質(zhì)模型參數(shù)Table 1 Parameters of hydrogeological model

3.1.6 源匯項(xiàng)處理

（1） 大氣降雨入滲補(bǔ)給量。

降雨入滲補(bǔ)給量是指大氣降雨滲入到土壤并在重力作用下滲透補(bǔ)給地下水的水量。年降雨入滲補(bǔ)給量的計(jì)算公式為：

式中：Pr為年降雨入滲補(bǔ)給量；P為年降雨量；a為年降雨入滲補(bǔ)給系數(shù)；F為計(jì)算區(qū)面積。

研究區(qū)多年平均降雨量為445 mm/a，降雨入滲系數(shù)參照研究區(qū)內(nèi)水文地質(zhì)資料，取值為0.18。

（2） 側(cè)向補(bǔ)給。

側(cè)向補(bǔ)給量用達(dá)西公式計(jì)算，公式如下：

式中：Q為側(cè)向補(bǔ)給量，m3/d；K為滲透系數(shù)，m/d；D為剖面寬度，m；M為含水層厚度，m；I為垂直于剖面的水力坡度。在模型中，根據(jù)補(bǔ)給邊界的以上參數(shù)自動(dòng)計(jì)算邊界的流入量。

（3） 側(cè)向排泄量。

側(cè)向排泄量用達(dá)西公式計(jì)算，公式如下：

式中：Q為側(cè)向排泄量，m3/d；K為滲透系數(shù)，m/d；D為剖面寬度，m；M為含水層厚度，m；I為垂直于剖面的水力坡度。

（4） 蒸發(fā)。

蒸發(fā)是指潛水（其中埋深小于4 m 時(shí)） 在毛細(xì)管力的作用下由下向上運(yùn)動(dòng)，最終以參加陸面蒸散發(fā)形式散逸到大氣中的水分損失量。此次研究區(qū)內(nèi)淺層地下水的埋深均超過(guò)了4 m，因此此次地下水蒸發(fā)量按零計(jì)。

3.1.7 模型的識(shí)別及驗(yàn)證

該階段是模擬過(guò)程中的重要部分，先對(duì)模型進(jìn)行識(shí)別，給出每個(gè)參數(shù)的范圍值，采用軟件自動(dòng)進(jìn)行模擬，期間也經(jīng)過(guò)手動(dòng)進(jìn)行調(diào)整，如果實(shí)際與計(jì)算水位相差較大，然后根據(jù)變化情況在進(jìn)行調(diào)試，直到兩者擬合較好為止。

根據(jù)研究區(qū)地下水分布規(guī)律、地下水動(dòng)態(tài)特征以及收集的地下水動(dòng)態(tài)資料，將識(shí)別期定為2018年5 月。地下水流場(chǎng)擬合如圖2 所示。由圖2 可知，在識(shí)別期實(shí)際地下水流場(chǎng)與模擬的地下水流場(chǎng)基本一致，可以認(rèn)為實(shí)際值與模擬值擬合情況較好。通過(guò)以上分析，此次在研究區(qū)建立的模擬地下水模型基本符合研究區(qū)的水文地質(zhì)條件。

圖2 地下水流場(chǎng)擬合圖Fig.2 Groundwater flowfield fitting

3.2 地下水溶質(zhì)運(yùn)移數(shù)值模型建立

（1） 數(shù)學(xué)模型的建立。

溶質(zhì)運(yùn)移的二維水動(dòng)力彌散方程的數(shù)學(xué)模型如下[8-10]：

式中：右端前二項(xiàng)為彌散項(xiàng)，后二項(xiàng)為對(duì)流項(xiàng)，最后一項(xiàng)為由于化學(xué)反應(yīng)或吸附解析所產(chǎn)生的溶質(zhì)的增量；D為彌散系數(shù)；μx，μy為x、y方向的實(shí)際水流速度；c為溶質(zhì)濃度：ML-3；Ω 為溶質(zhì)滲流的區(qū)域，量綱：L2；c0為初始濃度：ML-3。

（2） 彌散度的確定。

研究污染物在地下水中遷移轉(zhuǎn)化規(guī)律的最重要參數(shù)之一是彌散度，反映滲流系統(tǒng)彌散特征的一個(gè)綜合參數(shù)是彌散系數(shù)D，在忽略其分子的擴(kuò)散時(shí)，它是介質(zhì)彌散度僅和孔隙流速V的函數(shù)。在地下水溶質(zhì)運(yùn)移方程中，表征含水層介質(zhì)彌散特征的參數(shù)是水動(dòng)力彌散系數(shù)，它可表示為：

式中：αL，αT分別為縱向和橫向孔隙尺度彌散度，是僅與介質(zhì)特性有關(guān)的參數(shù)。結(jié)合區(qū)域水文地質(zhì)條件特征，確定潛水含水層縱向彌散度取值0.2 m。

（3） 模擬時(shí)段設(shè)定。

此次預(yù)測(cè)模擬工作以5 500 d 為模擬時(shí)間，研究污染物濃度時(shí)空變化過(guò)程與規(guī)律，更好的評(píng)價(jià)廠區(qū)污染物泄露后對(duì)地下水環(huán)境可能造成的直接影響和間接危害提供依據(jù)。

4 地下水污染物運(yùn)移模擬預(yù)測(cè)

模擬預(yù)測(cè)根據(jù)污染風(fēng)險(xiǎn)分析的情景設(shè)計(jì)，在選定優(yōu)先控制污染物的基礎(chǔ)上，分別對(duì)地下水污染物在不同時(shí)段的運(yùn)移距離和影響范圍進(jìn)行模擬預(yù)測(cè)。

根據(jù)對(duì)研究區(qū)模擬計(jì)算地下水的水質(zhì)情況和污染源的分布及類(lèi)型，選取對(duì)地下水環(huán)境質(zhì)量影響較大的污染物作為污染物溶質(zhì)運(yùn)移對(duì)象，通過(guò)標(biāo)準(zhǔn)指數(shù)法排序，確定了本次污染物，研究區(qū)污染物收集池中未200 mg/L。非正常工況下，假設(shè)廠區(qū)污水收集池防滲層發(fā)生破損，則確定廠區(qū)污水收集池為模擬泄漏點(diǎn)。研究區(qū)地下水現(xiàn)狀監(jiān)測(cè)指數(shù)0.8 mg/L，以此數(shù)據(jù)為邊界進(jìn)行計(jì)算。污染物發(fā)生泄露后，運(yùn)移模擬情況如圖3 所示。預(yù)測(cè)結(jié)果表明，在5 000 d 內(nèi)污染物羽狀污染范圍最大為0.2 km2，最大運(yùn)移距離436 m，在5 000 d時(shí)污染物濃度達(dá)到邊界值，對(duì)下游村莊不會(huì)造成影響。

圖3 污染物運(yùn)移模擬圖Fig.3 Simulation of pollutant transport

5 結(jié) 論

此次模擬計(jì)算是依據(jù)研究區(qū)內(nèi)地下水的水化學(xué)現(xiàn)狀、污染源的分布及類(lèi)型，選擇對(duì)地下水環(huán)境質(zhì)量影響較大的污染組分；在同樣的濃度和同樣體積的污水注入含水層的條件下，如果污染物指數(shù)含量不超標(biāo)，而其余污染物則更不會(huì)超標(biāo)。模擬結(jié)果表明，污染物在地下水運(yùn)移過(guò)程中，污染范圍由最初的近圓形逐漸變成長(zhǎng)軸方向與水流方向一致的近橢圓形。

污染物在地下水中隨水流運(yùn)移，隨著時(shí)間的變化，污染羽中心向下游運(yùn)移，運(yùn)移范圍逐步擴(kuò)大，但濃度降低，5 000 d 時(shí)污染物濃度接近臨界值，最大運(yùn)移436 m，對(duì)下游的村莊沒(méi)有影響。