某制藥廠廢水排放對地下水污染的數值模擬預測研究

楊 鑫,寧立波

(1.中國地質大學(武漢)地球科學學院，湖北 武漢 430074;2.中國地質大學(武漢)環境學院，湖北 武漢 430074)

制藥工業關系著國計民生，是經濟發展的支柱行業之一。但是，制藥行業屬于精細化工，制藥廢水的組成十分復雜，具有有機污染物種類多、濃度高，COD和BOD5值高，NH3-N濃度和固體懸浮物SS濃度高、色度深、毒性大等特征[1-3]。據統計，我國制藥行業廢水的排放量較大，約占全國工業企業污水排放量的2%[2]。制藥廢水在排放過程中污水管道發生滲漏可能造成的地下水污染問題一直備受關注，也是目前地下水污染中較難解決的問題之一，其中的關鍵問題在于地下水中污染物的運移特征及其范圍較難確定。數值模擬方法是解決該類問題的途徑之一，前人已做了大量的研究工作[4-6]，主要運用的數值模擬軟件包括Modflow、GMS、Tough2等[6-9]，但針對制藥廢水泄漏而產生的地下水污染范圍進行數值模擬的研究較少。新鄭市屬于鄭州市規劃的航空港區，其功能龐大，對環境質量的要求較高，但尚未進行地下水污染方面的環境影響評價工作。基于此，本文以擬進駐該區的某制藥廠為研究對象，運用GMS軟件對該制藥廠廢水排放過程中污水管道發生滲漏可能造成的地下水污染范圍進行了數值模擬研究，以為該地區地下水資源的合理開發與保護提供科學依據。

1 研究區概況

該制藥廠選址于河南省新鄭市規劃新區(見圖1)，該地區屬暖溫帶大陸性季風氣候，多年平均蒸發量為1 859.7 mm，多年平均降雨量為683.64 mm。研究區地處黃河和雙洎河沖洪積平原區，地形平坦，地表徑流遲緩，地下水埋深較淺，且包氣帶巖性大部分為粉土，局部為粉砂，結構松散，極有利于大氣降水的入滲補給。該地區地層以第四系為主，主要含水巖組為松散巖類孔隙含水巖組，其補給來源主要為大氣降水的入滲補給，地下水的水化學類型較簡單，以HCO3-Ca、HCO3-Ca·Mg型水為主。

圖1 研究區位置平面圖Fig.1 Location of the study area

2 研究區數值模型的建立

2.1 研究區水文地質條件與模型概化

水文地質概念模型是將含水層實際的邊界性質、內部結構、滲透性質、水力特征和補給排泄等條件進行概化，便于進行數學與物理模擬[4-5]。水文地質概念模型的核心要素是邊界條件、內部結構和地下水流態，通過對調查區的巖性構造、水動力場、水化學場的分析，可以確定水文地質概念模型的要素[6]。

2.1.1 研究區水文地質條件

研究區地下水屬松散巖類孔隙水，含水層由第四系上更新統和全新統的細砂、中砂、砂礫石層構成。研究區雙洎河河谷區的含水層巖性為砂、砂礫石，厚度為15～20 m，地下水水位埋深淺。雙洎河河谷區的外圍地帶和雙洎河北部支流的河谷地帶，含水層巖性主要為粉土，局部夾少量薄層粉細砂和卵礫石，厚度為20～30 m，地下水水位埋深為25 m左右；研究區內的其他地區地下水水位埋深大，含水層巖性為黃土，夾少量的砂層，富水性差。研究區的淺層含水巖組主要由黃河多次遷徒、改道、泛濫沖積而成，在含水層的顆粒粗細、埋藏分布、滲透性能等方面具有明顯的分帶性；含水層大體由西南向東北方向延伸，呈單層出現，局部地段夾有粉土或粉質黏土薄層透鏡體，厚度為20～40 m，個別地段厚度小于20 m；含水層頂板埋藏深度為5～50 m不等,巖性以粉土、粉砂為主，從上游至下游顆粒稍有變細，整體厚度約為60 m，屬潛水含水層，下伏粉質黏土、黏土隔水層。

研究區淺層地下水主要補給來源是大氣降水的入滲補給，此外雙洎河的側向滲漏也直接補給該區淺層地下水；淺層地下水的徑流方向大致與地形變化趨勢一致，總體方向由西南、西北山前向東、東南方向徑流，淺層地下水的排泄方式主要有人工開采、越流排泄和徑流排泄等，地下水埋深基本都大于4 m，蒸發排泄量很小，可以忽略。

2.1.2 研究區水文地質概念模型的建立

根據研究區的水文地質條件，構建了研究區的水文地質概念模型，見圖2。

圖2 研究區水文地質概念模型圖Fig.2 Conceptual model map of hydrogeology in the study area

(1) 模擬范圍：研究區模擬范圍為40 km2，由于該地區淺層地下水含水層的埋深約60 m，結合鉆孔資料，模擬垂直深度為60 m。

(2) 邊界條件：①水平邊界，即模擬區西與南左側邊界為定流量補給邊界，東、北與南右側邊界為定流量排泄邊界；②垂直邊界，模擬區垂向地下水補給包括大氣降水入滲補給、河流側向入滲補給、灌溉入滲補給；地下水排泄主要為側向徑流排泄與人工開采。

(3) 含水層結構特征：含水層大體由西南向東北方向或由西向東延伸，本次模型含水層為層1(粉細砂)，隔水層為層2(粉質黏土、黏土)。

研究區水文地質概念模型的水力特征參數見圖3，污染源為某藥業污水處理廠管道破裂滲漏形成，為一線形連續性污染源。

圖3 研究區水文地質概念模型的水力特征參數圖 (負值表示排泄)Fig.3 Hydraulic characteristic parameters of the study area hydrogeological conceptual model (negative values represent excretion)

2.2 數學模型的構建

有機污染物在地下水中的運移非常復雜，影響因素除對流、彌散作用以外，還存在物理、化學、微生物等作用，這些作用常常會使地下水中污染物總量減少，運移擴散速度減慢。本次模擬中假設污染物質在地下水運移過程中不與含水層介質發生反應，可以被認為是保守型污染物質。所以，本文選擇地下水流模型和溶質運移模型兩個模型，聯合求解水流方程和溶質運移方程，即可得到污染物質的空間分布。

2.2.1 地下水流模型

三維、非均質、各向異性的層流、非穩定潛層地下水流數學模型可表示為

式中：Kx、Ky、Kz分別為含水層在x、y、z方向的滲透系數(m/d)；Kxx為Kx在x方向上的投影；Kyy為Ky在y方向上的投影；Kzz為Kz在z方向上的投影;x，y，z方向為各向異性介質滲透系數的主方向;h為地下水水位(m);h0為地下水初始水位(m)；Kn為邊界上垂直于等水位線方向的滲透系數(m/d);Ω為研究區范圍；Γ2為邊界；w為降雨入滲補給量(m3/a)。

2.2.2 溶質運移模型

不考慮污染物在含水層中的吸附、交換、揮發、生物化學反應，地下水中溶質運移的數學模型可表示為

式中：α為含水層的彌散度(m)；vi，vj分別為i和j方向上的速度分量(m/d)；|v|為速度模；C為模擬污染物的濃度(mg/L)；ne為有效孔隙度；i,j為1，2，3，x1=x，x2=y，x3=z；t為模擬時間(d)；C′為源匯的污染物濃度(mg/L)；W為源匯單位面積上的通量；Dij為彌散系數(m2/d)。

2.3 模型識別與參數確定

2.3.1 模擬流場及初始條件

本次模擬以2016年5月的地下水流場作為初始流場(見圖4)，以2017年5月統測的地下水流場作為模擬流場。

圖4 模型初始地下水流場Fig.4 Initial groundwater flow field of the model

2.3.2 模擬區網格剖分

模擬區網格剖分單元格為50 m×50 m，其網格剖分見圖5。

圖5 模擬區網格剖分圖Fig.5 Mesh subdivision map of simulation area

2.3.3 模型識別與參數確定

模型的識別與驗證是整個模擬中極為重要的工作，通常需要在反復地調整參數和某些源匯項的基礎上才能達到較為理想的擬合結果[10]。本次模型識別與驗證過程采用試估-校正法，該法屬于反求參數的間接方法之一[11]。本次模擬時間為2018年1月到2022年12月，每個時間段內包括若干時間步長，時間步長由模型自動控制，嚴格控制每次迭代的誤差。模型最終識別的研究區含水層水文地質參數分區見圖6，模型識別的降雨入滲系數、含水層滲透系數和給水度見表1，灌溉入滲系數為0.10，該區地下水水位埋深普遍大于5 m，潛水蒸發排泄小可忽略不計。

圖6 研究區含水層水文地質參數分區圖Fig.6 Zoning maps of hydrogeological parameters of aquifers in the study area

編號降雨入滲系數水平方向滲透系數/(m·d-1)給水度層1層2層1層210.150.80.008 60.0470.03420.202.10.001 00.0470.03430.151.00.001 00.0470.03440.153.50.001 00.0340.03450.101.00.001 00.0340.03460.121.00.001 00.0400.03470.101.50.008 60.0470.034

研究區不同巖性介質彌散度的初值依據前人研究成果及經驗值確定，最終確定研究區不同巖性介質的平均縱、橫向彌散度，見表2。

表2 研究區不同巖性介質彌散度經驗參數表

3 模擬結果與分析

污染源位置為遂成藥業污水處理廠，假定污水管道出現長10 m、寬2 cm的裂縫，管道天然基礎層的滲透系數取值為0.047 m/d，則滲漏量約為10×0.02×0.047×1 000=9.4 kg/d，管道在無防滲設置的情況下發生滲漏，即為連續的線性污染源。根據污水中主要污染物監測指標，考慮污染物與巖土之間的反應和自身降解性，選取CODCr為地下水污染預測因子，根據遂成藥業污水處理廠進水水質狀況，進水中CODCr濃度為419 mg/L。

本文運用GMS軟件對擬建的遂成藥業污水處理廠廢水排放過程中污水管道發生滲漏可能造成的地下水污染范圍進行了數值模擬預測研究，模型運行時間起始于2018年1月1日，模擬預測了100 d、730 d 和1 825 d(2022年12月31日)三個時間段內研究區地下水中污染物CODCr的最大運移距離和污染暈的運移分布情況，其預測結果見表3和圖7。

污染組分的擴散具有彌散特征，初期是以污染源為中心向周圍彌散，之后隨地下水的流向污染暈的主軸方向與其趨向一致(向南運移)，所以，污染組分CODCr的檢出范圍以整個擴散范圍進行統計，最大運移距離以污染源(見圖7中顏色最深的中心點)為起點向東南方向的運移距離為準。預測結果表明：該制藥廠廢水排放過程中管道滲漏發生100 d后，地下水中污染物CODCr的檢出范圍為148 520 m2，超標范圍為0 m2，地下水中污染物CODCr的最大運移距離為238 m；滲漏發生730 d后，地下水中污染物CODCr的檢出范圍為165 500 m2,超標范圍為0 m2，地下水中污染物CODCr的最大運移距離為250 m；滲漏發生1 825 d后，地下水中污染物CODCr的檢出范圍224 000 m2，超標范圍為0 m2，地下水中污染物CODCr的最大運移距離為362 m(見表3)。受該地區西側與西南側兩個地下水降落漏斗的影響，該污染物CODcr滲漏后傾向于西南側運移且運移速度緩慢，至規劃年2022年末距東側地表水體距離較遠，對周邊水體的影響較小。

表3 研究區地下水中污染物CODCr的預測結果

圖7 研究區地下水中污染物CODCr的污染暈運移分布 預測圖Fig.7 Prediction map of CODCr pollution halo transport in groundwater in the study area

由圖7可見，研究區地下水中污染物CODCr的污染暈隨地下水的流向而擴散，其擴散方向與地下水的流向一致；地下水中污染物CODCr的運移距離隨時間的增加而增加，但整體上看在預測期內地下水中污染物CODCr的運移距離仍處在該制藥廠的建設區范圍內，不會對區域地下水造成明顯的污染。

4 結論與建議

(1) 制藥廠建設所排放的廢水含有復雜且大量的有機污染物，可能會對廠區范圍內的地下水造成污染，所以必須在制藥廠建設前對廢水排放可能造成的地下水污染范圍進行預測研究。

(2) 本文運用GMS軟件對擬建的遂成藥業污水處理廠廢水排放過程中污水管道發生滲漏可能造成的地下水污染范圍進行了數值模擬研究，模型運行時間起始于2018年1月1日，模擬預測了100 d、730 d和1 825 d(2022年12月31日)三個時間段內研究區地下水中污染物CODCr的最大運移距離。模擬結果表明：滲漏發生100 d后，地下水中污染物CODCr的最大運移距離為238 m；滲漏發生730 d后，地下水中污染物CODCr的最大運移距離為250 m；滲漏發生1 825 d后，地下水中污染物CODCr的最大運移距離為362 m。

(3) 研究區地下水中污染物CODCr的污染暈隨地下水的流向而擴散，其擴散方向與地下水的流向一致；地下水中污染物CODCr的運移距離隨時間的增加而增加，但整體上看，在預測期內不會對區域地下水造成明顯的污染。但是，鑒于該制藥廠污染物的極大危害性，建議在廠區范圍內及其周邊設置監測井，定期檢測地下水中污染組分的變化，一旦發現地下水受到污染，應立即采取相應的措施進行治理。