呂堂安 LV Tang-an；伍靜 WU Jing；黃偉 HUANG Wei；黃華娟 HUANG Hua-juan；朱敏杰 ZHU Min-jie

（①廣西大學(xué)資源環(huán)境與材料學(xué)院，南寧 530004；②中國有色桂林礦產(chǎn)地質(zhì)研究院有限公司，桂林 541004；③廣西中科地質(zhì)勘查有限責(zé)任公司，南寧 530022）

0 引言

數(shù)值模擬技術(shù)能直觀了解污染物泄漏后運移影響范圍及濃度變化，對地下水污染的防治有著極其重要的意義，其應(yīng)用的關(guān)鍵在于對模型的概化、邊界和參數(shù)的準(zhǔn)確選取[1]。Visual MODFLOW 由于簡便的操作性和強(qiáng)大的可視化能力，已成為國際上最流行的溶質(zhì)運移模擬軟件之一[2]。國內(nèi)外學(xué)者已用其對地下水污染運移開展了廣泛的研究。國外有：Morsy Samah Mahmoud[3-5]等人對埃及西米尼亞地區(qū)、蒂魯普和卡拉季地區(qū)地下水中污染物進(jìn)行預(yù)測，為開發(fā)區(qū)建設(shè)提供防治依據(jù)；國內(nèi)有：劉莉[6]針對淄博某化工項目，預(yù)測不同工況下的地下水中污染物溶質(zhì)運移變化規(guī)律；馬國明[7]運用Visual MODFLOW 對宿州市某制藥廠地下水中氨氮污染物進(jìn)行了模擬預(yù)測，并提出了防治措施。

某改、擴(kuò)建鈦白粉廠地處紅層區(qū)，周邊人口眾多，地下水是人們生產(chǎn)生活的重要水源，對該廠地下水的防治顯得極為敏感和重要。因此，本文以該鈦白粉廠為例，在詳細(xì)了解地下水水文地質(zhì)特征基礎(chǔ)上，概化了水文地質(zhì)模型，通過識別與檢驗，利用Visual MODFLOW 軟件建立了二維數(shù)學(xué)模型，開展特征污染物數(shù)值模擬，了解了該區(qū)地下水中污染因子隨時間的衰減及運移范圍，對指導(dǎo)該地區(qū)地下水的防治和保護(hù)有著重要意義。本次研究也為同類型地區(qū)如何獲取準(zhǔn)確的水文地質(zhì)參數(shù)，建立符合實際地下水流場的數(shù)值模型提供參考。

1 研究區(qū)域概況

1.1 地質(zhì)條件 研究區(qū)位于廣西博白旺茂鎮(zhèn)，為紅層丘陵盆地地貌，依據(jù)地形及巖層特性，將其劃分為兩個次級水文單元（圖1A）。由水文地質(zhì)剖面圖（圖1B）可知，研究區(qū)內(nèi)地層為下白堊統(tǒng)大坡組（K1d），上部為黃紅色粉砂土及部分出露含礫砂巖；下部為紅色砂巖夾礫巖、紫紅色砂巖含礫砂巖，中砂－細(xì)砂，中－粗礫，磨圓度較好，分選性差，成分各異，以砂巖、石英砂巖為主，少量灰?guī)r礫石。研究區(qū)內(nèi)無大的斷裂構(gòu)造，白堊系砂礫巖層產(chǎn)狀整體上較平緩。

圖1 研究區(qū)水文地質(zhì)簡圖（A）；A-A＇ 水文地質(zhì)剖面圖（B）

1.2 水文地質(zhì)條件 研究區(qū)內(nèi)含水層主要為白堊系紅層，表層風(fēng)化帶厚度一般為10～30m，紅層風(fēng)化裂隙較為發(fā)育，局部形成風(fēng)化層裂隙水，向下為沿構(gòu)造裂隙發(fā)育的中、弱風(fēng)化帶，其脆性的砂、礫巖在向斜軸部及構(gòu)造變形強(qiáng)烈部位，形成較多的節(jié)理裂隙，存在構(gòu)造裂隙水。鑒于紅層表層風(fēng)化裂隙水含水層厚度不大，與下伏構(gòu)造裂隙水之間無隔水層，存在直接水力聯(lián)系，故本次研究將含水層均概化為單層潛水含水層，地下水類型主要為碎屑巖類構(gòu)造裂隙水，其厚度為50～150m。淺層構(gòu)造裂隙水主要接受大氣降水的補(bǔ)給；排泄主要為蒸發(fā)和側(cè)向徑流補(bǔ)給。底部脆性構(gòu)造不發(fā)育的致密砂礫巖層膠結(jié)完好，無裂隙發(fā)育地段，視為隔水層。

該區(qū)地下水枯水期埋深一般在3m～8.65m，豐水期埋深一般在1.52m～5.74m，年水位變幅0.85～3.48m，水質(zhì)類型以HCO3—Ca 型為主。根據(jù)廠區(qū)機(jī)井J 的抽水試驗結(jié)果，該區(qū)含水層滲透系數(shù)為1.663×10-4cm/s，透水性中等，給水度μ 取0.085。鑒于本次研究區(qū)的紅層構(gòu)造裂隙含水介質(zhì)，其持水度較小，因此，有效孔隙度n 近似等于給水度，n=μ=0.085，水力坡度值I 為20.8‰，實際流速u 為0.035m/d。

2 地下水?dāng)?shù)值模型建立

2.1 水文地質(zhì)概念模型 根據(jù)地質(zhì)與水文地質(zhì)條件調(diào)查，將其概化為一個潛水含水層，淺層地下水運動以水平方向為主；地下水流向基本與下隔水層面平行，垂向分量忽略不計，符合二維流特征，因此概化為平面二維流[8]。本次模擬將區(qū)域地下水含水系統(tǒng)概化為非均質(zhì)、各向同性二維無壓非穩(wěn)定流。研究區(qū)含水層厚度為50～150m，西側(cè)以老屋水庫為地下水排泄邊界，視為定水頭邊界；東側(cè)及西北側(cè)以地下水分水嶺為邊界，視為二類邊界；東北側(cè)以金甲河為北側(cè)排泄邊界，視為定水頭邊界；上邊界為降水補(bǔ)給、蒸發(fā)排泄邊界，下邊界以基巖作為相對隔水邊界。

2.2 地下水流數(shù)學(xué)模型 在實際情況中，由于缺少三維實測數(shù)據(jù)，在必要時須做一些假設(shè)，這可能造成三維模型的結(jié)果失真[9]。在實際工作時，當(dāng)含水層與弱透水層或隔水層的貯水系數(shù)比值大于10，而弱透水層或隔水層越流系數(shù)小于1·10-6m/d 時可用二維模型[10]。本次研究建立各向同性二維非穩(wěn)定流數(shù)學(xué)模型，即式（1）[11]：

式中：μd為重力給水度，l/m；Kx、Ky分別為x、y 方向上的滲透系數(shù)，m/d；t 為時間，d；h 為潛水面高度，m；z 為含水層隔水底板高度，m；ε 為含水層匯源項，m3/d；Ω—滲流區(qū)域；H0為滲流區(qū)Ω 內(nèi)各點的初始水頭，m；Γ1為滲流區(qū)域Ω 的第一類邊界；H1為第一類邊界Γ1上的已知水頭函數(shù)，m；Γ2為Ω 的二類邊界；K 為邊界Γ2法線方向的滲透系數(shù)，m/d；n 為邊界Γ2上某點處的外法線方向；q 為Γ2上的單位面積流量函數(shù)，m3/d。

本次模擬網(wǎng)格剖分精度為20m×20m，在分水嶺及污染物滲漏處加密網(wǎng)格為10m×10m，模擬面積約為2.17km2。研究區(qū)主要補(bǔ)給來源為大氣降雨，老屋水庫及金甲河為主要排泄邊界。降雨蒸發(fā)參數(shù)根據(jù)《玉林市地下水利用與保護(hù)規(guī)劃（2016-2030）》，多年平均逐月降雨量為1583mm，多年平均蒸發(fā)量為1506mm，蒸發(fā)臨界深度取5m。

2.3 數(shù)值模型動態(tài)識別與檢驗 本次模擬以枯水期水位作為初始水位，按照差分法得到含水層初始流場，枯水期至豐水期為模型識別期。根據(jù)水文地質(zhì)試驗設(shè)置含水層滲透系數(shù)Kx和Ky初始值為1.663×10-4cm/s，采用試估—校正法對模型進(jìn)行檢驗，獲得與實際流場相符程度高的模擬流場。將枯、豐水期實測點位數(shù)據(jù)與模擬水位進(jìn)行檢驗，置信區(qū)間相關(guān)系數(shù)分別為0.995 和0.972，表明模擬流場精度較高，適用于本次模擬研究。最終確定研究區(qū)水文地質(zhì)單元Ⅰ和Ⅱ的溶質(zhì)模擬滲透系數(shù)Kx、Ky取值為1.039×10-4cm/s。

2.4 溶質(zhì)運移數(shù)學(xué)模型 本文從保守角度考慮，本次溶質(zhì)模擬計算的模型識別和計算沒有考慮污染物在含水層中的吸附、揮發(fā)、生物化學(xué)反應(yīng)，考慮最大值[12]。地下水以水平運動為主，垂向彌散作用忽略不計。在此前提下，本次調(diào)查研究建立地下水二維溶質(zhì)運移模型，二維水動力彌散方程的數(shù)學(xué)模型見式（2）[11]：

式中：C 為模擬污染物質(zhì)的濃度，mg/L；θ 為地層介質(zhì)的孔隙度，無量綱；t 為時間，d；xi為沿直角坐標(biāo)系軸向的距離，m；Dij為水動力彌散系數(shù)，m2/d；vi為孔隙水平均實際流速，m/d；qs為單位體積含水層流量，代表源和匯，m3/d；Cs為源或匯水流中組分的濃度，mg/L。

彌散系數(shù)結(jié)合水文地質(zhì)條件和參考《水文地質(zhì)手冊》[13]，最終選確定為8m2/d。

3 模擬預(yù)測

根據(jù)地質(zhì)和水文地質(zhì)調(diào)查，滲漏點位于鈦白粉廠廠區(qū)西部，屬于水文地質(zhì)單元Ⅰ，地下水流向為北東向西徑流，污染物滲漏后將向西側(cè)老屋水庫一帶運移，因此本次模擬重點關(guān)注水文地質(zhì)單元Ⅰ。

3.1 污染源設(shè)定 由工藝分析可知：鈦白粉廠產(chǎn)生的酸性廢水中主要污染物為COD、氨氮和氟化物，模擬滲漏區(qū)為鈦白粉污水處理車間和硫鐵礦制酸污水處理站。根據(jù)《環(huán)境影響評價技術(shù)導(dǎo)則地下水環(huán)境》（HJ 160），當(dāng)系統(tǒng)老化等原因?qū)е抡{(diào)節(jié)池滲漏≥2L/m3·d 時視為非正常狀況，結(jié)合工程驗收監(jiān)測數(shù)據(jù)，計算得出鈦白粉污水處理車間污水滲漏量為864L/d；鐵礦制酸污水處理站污水滲漏量為96L/d，預(yù)測污染物持續(xù)滲漏60 天可發(fā)現(xiàn)。

3.2 污染物遷移模擬預(yù)測結(jié)果 ①COD 預(yù)測評價。滲漏液初始濃度為417mg/L，持續(xù)滲漏60 天時調(diào)節(jié)池底部地下水COD 濃度達(dá)到最大值4.80mg/L；第100d 時COD中心濃度最大為3.9mg/L，滲漏點下游0～12m 范圍內(nèi)COD濃度高于地下水Ⅲ類標(biāo)準(zhǔn)，溶質(zhì)運移范圍0～25m，未出廠界（見圖2A）；滲漏第1000d 時，污染物中心濃度最大為1.00mg/L，污染物運移至下游70m（見圖2B）；在滲漏第7300d 時污染物COD 已運移至廠界外，到達(dá)西側(cè)老屋水庫一帶，中心濃度最大為0.05mg/L，濃度低于地下水Ⅲ類標(biāo)準(zhǔn)（見圖2C），擴(kuò)散距離為75～269m，位于下游的敏感目標(biāo)地古嶺民井距研究區(qū)420m，不受到研究區(qū)滲漏影響。滲漏后COD 平均運移速度為13.45m/a。

圖2 COD 運移路徑及濃度變化圖；100d（A）；1000d（B）；7300d（C）

②氨氮預(yù)測評價。滲漏液初始濃度為33.75mg/L，持續(xù)滲漏60 天時調(diào)節(jié)池底部地下水氨氮濃度約0.38mg/L；滲漏第100d 時氨氮中心濃度最大為0.30mg/L，污染物運移至下游0～6m 范圍內(nèi)，未出廠界（圖3A）；滲漏第365d 時調(diào)節(jié)池底部污染物氨氮濃度為0.21mg/L，污染物運移至下游0～41m 范圍內(nèi)，未出廠界（圖3B）；滲漏第730d 時調(diào)節(jié)池底部污染物氨氮濃度為0.05mg/L，污染物運移至下游0～58m 范圍內(nèi)（圖3C），遠(yuǎn)低于地下水Ⅲ類標(biāo)準(zhǔn)（0.5mg/L）。滲漏后氨氮平均運移速度為2.9m/a。

圖3 氨氮運移途徑及濃度變化圖；100d（A）；365d（B）；730d（C）

③氟化物預(yù)測評價。滲漏液初始濃度為195.25mg/L，持續(xù)滲漏60 天時調(diào)節(jié)池底部地下水氟化物濃度約2.4mg/L；滲漏第100d 氟化物中心濃度最大為1.90mg/L，污染物運移至下游0～30m 范圍內(nèi)，未出廠界（圖4A）；滲漏第365d 污染物中心濃度為1.10mg/L，污染物運移至下游0～41m 范圍內(nèi)，未出廠界（圖4B），滲漏第1000d 污染物中心濃度為0.50mg/L，污染物運移至下游0～56m 范圍內(nèi)（圖4C），此時濃度已經(jīng)低于地下水Ⅲ類標(biāo)準(zhǔn)（1.00mg/L）。滲漏后氟化物平均運移速度為2.8m/a。

圖4 氟化物運移途徑及濃度變化圖；100d（A）；365d（B）；1000d（C）

4 結(jié)論

本次研究表明在符合一定水文地質(zhì)條件下的紅層地區(qū)，采用Visual MODFLOW 能建立與實際流場相符程度較高的二維流模型。在此基礎(chǔ)上進(jìn)行模擬預(yù)測，表明該區(qū)地下水在鈦白粉廠發(fā)生滲漏后對下游紅層地下水有一定范圍影響，但影響未擴(kuò)散至地谷嶺村屯。企業(yè)應(yīng)重點對污水站下游和廠區(qū)、地下水下游邊界監(jiān)測水點進(jìn)行監(jiān)測，隨時掌握不同時期地下水的水質(zhì)水位動態(tài)資料。