基于Visual Modflow的某化工項目地下水污染預測

魏慶龍，楊金偉，靳廷朝

（河南省地礦局第一地質環境調查院，河南鄭州450045）

在地下水資源的定量分析研究中，建立能較真實反映水文地質結構與流場特點的地下水數值模擬是地下水資源評價與管理的重要研究方法[1-2]，基于Visual Modflow軟件的有限差分法可以滿足復雜地下水的數值模擬[3]。

1 模擬區水文地質概況

模擬評價區處于豫北黃河沖積平原，不是一個完整的水文地質單元，地下水可以劃分為松散巖類孔隙水、半膠結碎屑巖類裂隙水。其含水層組按埋藏條件，可進一步劃分為淺層含水層組、中深層含水層組和深層含水層組。將淺層含水層和中深層含水層之間的粘土層當作模型的隔水底板。淺層含水層底板埋深90～120m，第一個連續穩定隔水層之上的上更新統與全新統含水層，淺層水為潛水—微承壓水，含水介質為沖洪積及風積的粗、中、細、粉細砂組成，可見3層，單層厚度20～30m，總厚度為60～90m，由南向北，含水層厚度變薄。單井出水量為1800～2400m3/d，滲透系數13～20m/d，為中等富水區。

淺層地下水排泄主要以垂直蒸發為主，近些年來人工開采已成為重要的排泄方式，其次是向下游的徑流排泄，以及向中深層地下水的越流補給。大氣降水為本區淺層水的主要補給源。引黃渠道的滲漏補給也是本區淺層水的補給源。淺層水化學類型為HCO3·SO4-Na·Mg。

本區淺層地下水動態類型屬氣象-水文型。區內淺層地下水水位埋深差異較小，一般在3～7m之間，年變幅一般在0～4m左右。水位動態不全受降水的制約，而且又受引黃灌溉渠道水回滲因素影響，年內水位出現2個峰值，與降水和引黃灌溉時間相吻合。

2 地下水水動力場數值模擬

水是溶質運移的載體，地下水流場是溶質運移模擬的基礎，在溶質運移模擬前，需先建立模擬區地下水流場模型，地下水流場預測模型為modflow。

2.1 地下水流數學模型

根據模擬區水文地質特征，可將模擬區淺層孔隙水含水層概化為非均質、各向異性、三維非穩定地下水流系統，其運動的數學模型為：

式中：Ω——滲流區域；

Kxx、Kyy、Kzz——滲透系數在x、y和z方向上的分量，假定滲透系數的主軸方向與坐標軸的方向一致，m/d；

h——水位（水頭）標高，m；

ε——源匯項，m2/（d·m），ε＜0表示流出地下水系統，ε＞0表示流入地下水系統；

Ss——含水層給水性能參數，取重力給水度；

從圖2可以看出，常規治理方案一般需將脫硫裝置出口濕煙氣從50～55 ℃加熱至75～80 ℃(煙道及煙囪溫降按5 ℃考慮)，需要較大的熱量，現在電廠普遍采用除塵器前或脫硫塔前的高溫煙氣作為加熱熱源，高溫煙氣余熱無法再進行節能綜合利用，因此，治理方案應盡量減少飽和濕煙氣的加熱熱量，并考慮煙氣余熱的綜合利用效果。

t——時間，d；

h0——含水層的初始水位，m；

Γ2——流量邊界；

Γ3——淺層含水層的底板；

ˉ——邊界面的法線方向；

q(x，y，z，t)——二類邊界的單寬流量，m3/（d·m），流入為正，流出為負。

2.2 地下水流數值模擬的建立

對整個區域模型采用不等距矩形網格剖分，其中廠區進行局部加密剖分，其他區域等距離剖分，剖分為94行133列，共剖分矩形網格單元12502個，其中有效單元11118個，無效單元1384個。模擬區網格剖分見圖1，模型結構剖面圖見圖2。

圖1 模擬區網格剖分圖

2.3 模型識別

圖2 模型結構剖面示意圖

模型識別過程通常要在反復修改水文地質參數和調整某些源匯項基礎上才能達到較為理想的擬合結果。此模型的識別過程采用的方法稱為試估—校正法，屬于反求參數的間接方法之一。運行計算程序，可得到這種水文地質概念模型在給定水文地質參數和各均衡項條件下的地下水位時空分布，通過擬合同時期的地下水水位，識別水文地質參數、邊界值和其它均衡項，使建立的模型更加符合模擬區的水文地質條件。

根據所掌握的資料，以2013年底流場作為初始流場（見圖3），識別期選為2014年1月～2015年6月，2015年7～12月為驗證期，應力期以月為單位，共劃分為24個應力期。地下水位擬合效果見圖4，可見地下水位模擬水位和實測水位誤差均小于0.5m，所建立的數學模型達到了模型精度要求，可以利用該模型對研究區進行地下水污染情景預報。

3 地下水污染模擬預測

本次地下水污染模擬過程未考慮污染物在含水層中的吸附、揮發、生物化學反應，模型中各項參數予以保守性考慮。溶質運移模型采用MT3DMS進行求解。

3.1 溶質運移數學模型

地下水中溶質運移的數學模型為：

式中：ρb——介質密度，mg/(d·m)3；

θ——介質孔隙度，無量綱；

C——組分的濃度，mg/L；

t——時間，d；

x、y、z——空間位置坐標，m；

Dij——水動力彌散系數張量，m2/d；

Vi——地下水滲流速度張量，m/d；

Cs——組分的濃度，mg/L。

圖3 初始流場

3.2 地下水污染預測情景設定

通過對廠區平面布置、污水收集排放系統等綜合分析，確定規劃區污水處理站為此次模擬泄露點。根據廠區規劃特點，特征污染因子選為COD，地下水污染主要預測非正常工況防滲不正常情景

而非正常工況下，構筑物的最大滲水量取100倍的正常工況的滲水量進行計算。由于廠區有較完備的污染事故監測預警系統，考慮不正常工況持續的時間為2d，非正常工況排除后，污水池滲水量按2L/（m2·d）進行計算。

3.3 場區生產對地下水水質的影響評價

根據非正常工況防滲不正常情景下污染模擬預測結果，依據《地下水質量標準》（GB/T14848-1993）中CODMn的Ⅰ類和Ⅲ類水的要求，特征污染因子在1、5、10、20、30年對地下水的影響范圍、超標范圍見表1和圖5。第1年CODMn濃度最大值分別為12.478mg/L；第30年CODMn濃度最大值分別為82.0544mg/L。第30年CODMn影響范圍已到達廠區邊界。此外，第27d，CODMn濃度超過地下水Ⅰ類標準，其濃度最大值為1.0288 mg/L，如果以廠區范圍作為污染控制范圍、以CODMn作為控制性污染因子，則非正常工況排除的時間不超過2d。

圖4 地下水位擬合效果圖

4 結論

根據研究水文地質資料和地下水觀測資料，識別和驗證了區域地下水運動水流運動模型，數值模擬模型具有較高的模擬精度。

根據廠區規劃特點，特征污染因子選為COD。地下水污染預測情景設定為非正常工況，構筑物的最大滲水量取200L/（m2·d），不正常工況持續的時間為2d，非正常工況排除后，污水池滲水量按2L/（m2·d）進行計算。

若以《地下水質量標準》（GB/T14848-1993）中Ⅰ類和Ⅲ類水CODMn濃度標準值分別作為地下污染影響范圍和超標范圍，則非正常工況下，第1年CODMn濃度最大值為12.478mg/L；第30年CODMn濃度最大值為82.0544mg/L。第30年CODMn影響范圍已到達廠區邊界。此外，第27d，CODMn濃度超過地下水Ⅰ類標準，其濃度最大值為1.0288mg/L，如果以廠區范圍作為污染控制范圍、以CODMn作為控制性污染因子，則非正常工況排除的時間不超過2d。

圖5 非正常工況下CODMn影響范圍

表1 非正常工況下污染物對地下水的影響范圍

[1]薛禹群，謝春紅.地下水數值模擬[M].北京:科學出版社，2007.

[2]錢家忠，吳劍鋒，朱學愚，等.地下水資源評價與管理數學模型的研究進展[J].科學通報，2011,46（2）:99-103.

[3]張洪霞，宋文.地下水數值模擬的研究現狀與展望[J].水利科技與經濟，2007,13（11）:794-796.

[4]郭瑞.興平市地下水動態特征及數值模擬研究[D].西北農林科技大學，2008:5-8.

[5]徐娟花.地下水數值模擬現狀[J].工程地質與計算機應用，2012(4):10-14.

[6]陳正國.數值模擬在地下水資源評價中的應用研究[J].黑龍江水利科技，2012,34(12):69-74.

[7]孫從軍.地下水數值模擬的研究與應用進展[J].環境工程，2013,31(5):9-17.

[8]孫道林.典型垃圾堆放場對地下水污染模擬及控制技術研究[D].合肥工業大學，2008:23-35.

[9]劉娟，楊國勇，蔣兆華，等.FEFLOW在某市巖溶地下水四氯化碳污染中的應用[J].地下水，2011,33(4):87-91.