鄄城縣工業園區污染物下滲及遷移預測研究

畢桂真 孟令杰 韋愛菊

（山東省菏澤市水文局 菏澤 274008）

1 地下水流數值模型的建立及求解

1.1 水文地質概念模型

從含水層結構上看，研究區為多層結構區。據實測資料，園區內潛水埋深一般為1.51～4.73m，含水層厚度為14.6m，潛水層之間夾有平均厚度在5～10m的粘質砂土、粉質粘土弱透水層。此次研究把潛水作為計算目的層，通過分析鉆孔資料和水文地質調查資料，將含水層概化為非均質、各向同性含水層。所以模擬區地下水的水動力條件概化為穩定的三維流。

1.2 地下水流數學模型

根據前述的水文地質概念模型，研究區地下水運動數學模型如下：

式中：H—地下水位標高（m）；H0—第一類邊界水位標高（m）；t—時間（t）；A—模擬區四周邊界；K—滲透系數（m/d）；Z—潛水面水頭分布值（m）；x，y，z—坐標變量（m）；rw—抽水井的半徑（m）；SS—含水層的釋水率（L-1）；Qi—第 i口井的抽水流量（m3/d）；W—研究區三維流的源匯項（L/d）；μ—給水度；U—滲流研究區域；f（x，y，z）—初始水位分布（m）；（r，θ）—輔助柱坐標變量。

1.3 模型離散化及基礎資料的給定

1.3.1 空間的離散化

研究范圍是一個規則的矩形區域。將研究區垂向上剖分為4層，其中潛水含水層1層，弱透水層1層，微承壓含水層1層和隔水層1層，由于研究區地處平原，地形平坦，面積較小，地形高程和各層厚度以平面的形式輸入到模型中。

1.3.2 時間的離散化

選取2017年4月1日為模擬的起始時間，7300日后為模擬的終止時間。最大時間步長為200d，嚴格控制每次迭代的誤差。在每個應力期保持含水層補給和排泄強度不變。

1.3.3 初始條件的輸入

初始水位采用2017年4月統測的分層觀測水位，將實測水位輸入到模型中從而得到研究區數值模擬的初始流場，其中弱透水層的初始水位取上下兩個含水層水位的算術平均值。

1.3.4 邊界條件的輸入

四周邊界概化為第一類邊界條件，每月1日水位值作為水位觀測值輸入到Visual MODFLOW中。

1.3.5 水文地質參數的初值

根據2017年3月16日穩定流抽水試驗資料和水文地質勘察室內土工試驗資料，結合研究區地質、水文地質條件，將以上資料所得參數輸入到地下水流模型中。

1.3.6 源匯項的處理

潛水的補給來源包括降雨入滲補給、灌溉入滲補給和渠系滲漏補給，其中降水入滲補給占主導地位。地下水排泄方式包括側向徑流排泄、潛水蒸發排泄和人工開采。補給和排泄按補給強度和開采強度處理，人工開采按雙井開采流量計算。降水采用鄄城縣多年平均降水量609.7mm，蒸發強度采用魏樓閘水文站多年平均水面蒸發量，潛水蒸發極限埋深參考鄧集試驗站成果取4.0m。

1.4 數學模型的識別、驗證

1.4.1 模型識別

把各種水文地質資料代入模型，以迭代殘差和最大水頭改變量最小為目標。通過調整分區參數值使二者盡量小，并據此來判斷所用水文地質參數及分區是否合理。經反復調整參數，獲得了較為滿意的水文地質參數。

1.4.2 模型驗證

選用2016年（年降水量611.0mm）作為驗證年份，以周邊11#、15B#、32A#多年地下水動態監測井作為參照井進行驗證。經驗證，各監測井實測水位與計算水位差值的絕對值小于0.05m，表明各觀測孔的水位計算結果與實測結果吻合很好，驗證了所取參數的合理性。同時對比2017年4月潛水的實測水位與模型計算水位，模擬流場與實際流場的變化趨勢基本一致，可以運用到地下水水質模型中。

1.5 地下水流數學模型預測

在考慮各項補給項的條件下，使用建立的地下水流數學模型，假定園區工業取用水戶雙井穩定流開采地下水，單井出水量1200m3/d。連續抽水48h后實際觀測記錄與預測結果見表1。

2 地下水水質模型及預測

通過Visual MODFLOW中的MT3DMS模塊計算污染物質的運移情況，可以求出污染物在地下水系統中的變化規律，預測研究區污染物質在不同時刻、不同的情況下所導致的地下水污染程度。

2.1 溶質運移數學模型

根據研究區的具體條件，采用下述的溶質運移模型：

式中：C—溶解于水中的污染物濃度；n—孔隙度；xi—空間坐標；qw—源（正值）或匯（負值）的單位流量；Dij—水動力彌散系數張量；t—時間；F—固相表面的溶質濃度；C0—源匯項的濃度；Ω—空間區域；Vi—地下水滲透流速；（x，y，z）—空間位置；r—研究區邊界。

2.2 水質模型的建立

水質模型以水流模型為基礎建立，水質模型的概化與所建立的水流概念模型相符。水質模擬區范圍、含水層結構、邊界類型劃分、源匯項的概化均與水流概念模型相同，流體概化為不可壓縮的均質流體，粘度和密度均為常數。

2.2.1 模擬因子的選擇

此次研究為穩定流情況下預測研究區污染物的運移情況，考慮污染物質在含水層中的線性等溫吸附（平衡），無反應項，選擇化學需氧量（COD）為示蹤劑。

2.2.2 邊界條件、初始條件

研究區水質模型選定2017年4月1日作為初始時刻，初始時刻含水層中污染物的濃度為0。

表1 園區抽水試驗觀測記錄與預測結果對照表

2.2.3 模型參數

溶質運移模型涉及的參數中含水介質的有效孔隙度（n）由試驗所得n=0.19，其他參數的取值如下：

（1）地下水滲流速度：按照2017年3月16日實際抽水試驗結果，含水層地下水垂向滲透系數速度為1.0521×10-4cm/s，由2017年3月26日室內土工試驗，弱透水層地下水垂向滲透系數確定為2.73×10-6cm/s，隔水層為9.94×10-8cm/s。

（2）縱向彌散系數的確定：溶質運移模型需要的彌散系數是在結合抽水試驗進行的野外彌散實驗，同時取研究區的土樣，采用一維土柱彌散實驗法，進行室內土柱彌散實驗而確定的。經過計算，評價區縱向彌散系數為1.46m/d。

（3）其他參數選用經驗值。

3 預測結果

假定園區發生污染事故，污染物下滲污染當地地下水。以地面表層受到一定濃度COD污染作為初始條件。在假定的園區工業取用水戶雙井穩定流開采地下水的條件下，通過Visual MODFLOW中的MT3DMS模塊和已建立的溶質運移模型，預測不同時間段內污染物垂直和水平運移的情況。設定污染羽外圍污染物的濃度達到0.5mg/L時為該污染物的影響范圍。

3.1 不同時段污染物的垂直運移情況

以地面表層受到濃度為5000mg/L的COD污染為初始條件，經過365d（1年）的垂直運移，地下水已經受到了輕度的污染。此時，在潛水的表層有仍有少許積累外，污染物已經到達潛水含水層的中下部，中心位置的濃度達到100mg/L。經過5475d（15年）的垂直運移，污染物在垂向上已經能夠到達隔水層的下部，實現了對隔水層的滲透性運移。運移7300d（20年）后，垂向上的影響范圍沒有顯著的增加，但潛水表層的污染物濃度明顯增加，污染物在橫向和縱向上得到了積累，地下水中的污染物運移范圍進一步擴大，濃度遞增顯著，中心位置處不但濃度高達1000mg/L，而且范圍擴大了近10倍。

3.2 不同時段污染物的水平運移情況

為對比不同濃度下污染物的水平運移情況，分別對初始濃度5000mg/L、500mg/L的COD做污染物的運移模擬，對比在7300d后的運移等值線，5000mg/L濃度COD的10mg/L等值線運移了582m，污染羽（濃度為0.5mg/L）運移了680m；500mg/L濃度COD的10mg/L等值線運移了426m，污染羽運移了606m。說明污染物濃度越高，污染羽在相同時段運移距離越遠速度越快。

4 結論與建議

根據鄄城縣工業園區水文地質條件的調查，運用Visual MODFLOW模型對工業園區受到污染物污染后污染物垂直和水平運移的模擬計算，得出以下結論：

（1）當園區受到污染物污染后，隨著時間的推移，園區地下水從受到輕度污染逐步演變到重度污染。從初期污染物在垂向上到達隔水層的下部并滲透性運移過隔水層，到后期地下水中的污染物運移范圍進一步擴大，濃度遞增顯著，污染物濃度逐年線性增加，若不及時治理，污染范圍會進一步擴大。隨著抽水時間的延長，地下水埋深逐年增加，地下水位逐年降低，污染羽的運移速度在第二含水層顯著增加。

（2）園區場地基礎之下第一巖土層不能滿足天然防滲層的要求，園區內事故污染下滲會造成當地地下水的污染。鑒于園區地層不能滿足最低滲透性標準，建議采用在園區四周和底部鋪設土工膜的方法，進行人工水平防滲和垂直防滲，防止由于園區受到污染物污染對地下水水質造成影響