基于MODFLOW 模型滹沱河傍河地下水源地保護區劃分

冷文鵬，陶 亞，孫若涵，夏建新

(1. 中央民族大學 生命與環境科學學院，北京 100081；2. 生態環境部環境規劃院，北京 100012)

我國飲用水源地面臨水質污染超標、突發水污染事件頻發等安全隱患，特別是具有水量豐富、水質比較穩定等特點[1]的地下水水源地被污染后不易修復，因此，要加強對地下水飲用水源地的保護。為保護地下水水質，建立地下水源地保護區是必不可少的措施[2]。但是，地下水流動比較復雜，地下水源地保護區劃分存在一定的技術困難。利用數值模擬法確定地下水水源中的污染物運移方向和距離，可為保護區劃分提供依據。鄧媛媛等[3]利用FEFLOW 軟件建立了吳忠市金積傍河型地下水水源地的三維地下水數值模型，將水源地的12 眼水井概化成“一眼、四眼、十二眼”，模擬了3 種地下水開采情景下地下水的運動，簡化復雜的邊界條件，并利用反向示蹤粒子跟蹤技術，建立了三級保護區。王涵等[4]在內蒙古呼和浩特市城市水源地運用MODFLOW 建立了該水源地的地下水流場，根據148 口抽水井單獨設置了148 個反向示蹤粒子點位，計算出反向示蹤粒子的運移跡線，確定了水源地各級保護區的界線。肖杰等[5-9]運用Visual MODFLOW 中的modflow 和modpath 模塊分別對崇州市地下水源地、渭河傍河水源地、關中盆地水源地進行了數值模擬，按照抽水井位置設置單個示蹤粒子，對其進行保護區的劃分。目前邊界條件的相關研究方向集中在單井概化，缺少對井群概化的研究和應用。邊界條件處理方法的差異可能導致模擬結果不同，按照單井概化處理的研究結果可能使保護區范圍過小，帶來污染風險。因此，通過研究兩種不同的邊界條件，分析其計算結果，對于提高保護區劃分結果的科學性和準確性具有重要意義。

本文主要以石家莊市新調整后的滹沱河應急后備地下水源地為例，設定了單井和井群兩種邊界條件，基于MODFLOW 軟件計算分析地下水流場，運用數值模擬法計算不同工況下反向示蹤粒子100 和1 000 d的運移距離，并分析入滲場使用情況下的示蹤粒子運移規律。

1 研究區概況

滹沱河應急后備地下水源地位于石家莊市新華區，西望太行山，北依滹沱河。按照供水計劃該水源地要保持供水規模10×104m3/d，布設開采井17 眼，單井開采量6 000 m3/d，井眼內直徑800 mm，利用段為第四系第Ⅱ+Ⅲ含水組，井距800 m，水源地面積約為10 km2。

石家莊市西部地處太行山中段，東部為滹沱河沖積平原。區內面積為1 521.08 km2，地勢為西北高東南低，地面坡降(1.6～2.5)‰。滹沱河流域的地下水主要賦存在第四系松散巖類孔隙中，含水層多由亞砂土、砂、卵礫石組成，粒度粗、厚度大，水動力特征為潛水、微承壓水。水源地擬取水的含水層為第Ⅱ含水組和第Ⅲ含水組。第Ⅰ含水組埋藏深度為15～20 m，該層沉積較薄，顆粒較細，巖性為粉土、細、中、粗砂及含礫石，由于當地長期取水，該層基本疏干。第Ⅱ含水組底板埋藏深度約為100 m，含水層厚度30～50 m，該層沉積厚度大，含水層顆粒較粗，且磨圓度較好；主要巖性為砂礫、卵礫石，透水性及富水性好。該層分為上、下兩段，尤以下段含水層最為豐富，單位涌水量30～40 m3/(h·m)，滲透系數一般為37～145 m/d；地下水質良好，礦化度小于0.5 g /L。第Ⅲ含水組底板埋藏深度約為220 m，自西北向東南傾斜，含水層厚度大于50 m，巖性含礫卵石、砂礫夾砂質黏土，其中砂卵石、砂礫石分選較差，該層在經濟技術開發區以西遭受了不同程度的風化，透水性和富水性均較差；開發區以東富水性較好，受地方開采井連通影響，本區水動力特征屬潛水-微承壓水，礦化度小于0.5 g/L。

2 數值模擬法劃分保護區

2.1 軟件模擬過程概要

MODFLOW 軟件由20 世紀80 年代美國地質調查局（USGS）的Mc Donald 和Harbaugh 研發，基于有限差分法原理用于孔隙介質中地下水流動的數值模擬。本研究采用MODFLOW 6.0.3 計算地下水流場，并通過MODPATH（粒子跟蹤后處理軟件程序）計算示蹤粒子的運移跡線。首先按照特定的格式整理空間展布資料（地表高程、分層數據、邊界位置）、含水層參數、初始水位、源匯項（大氣降水入滲系數、蒸發排泄系數、人工開采量）及水位動態觀測資料。隨后，根據水文地質剖面圖確定含水層數，根據以上資料輸入新建模型中的各含水層底板高程數據、初始水流場數據、含水層參數。其中，人工開采量和側向徑流補給量分別以抽水井、注水井的方式確定[10]，分別按照所需格式輸入數據。輸入11 個水文地質參數分區數值，使用PEST 模塊調整水文地質參數，使得模擬值和實測值的差值在誤差允許范圍之內，運行MODFLOW，得出符合實際條件的地下水流場數值。最后通過MODPATH 模塊設置反向示蹤粒子，計算100 和1 000 d 粒子運移的距離。

2.2 地下水流動模型

以滹沱河傍河地下水井群為中心，將研究區邊界條件、徑排補給、各區滲透系數、給水度等要素進行概化，確定研究區域范圍為30 km×45 km。石家莊市區17 口應急后備水井位于滹沱河河心島處，取水第Ⅱ、Ⅲ含水層即為潛水和微承壓水，因此本次模擬的主要含水層為微承壓水層。

MODFLOW 共有自然鄰點插值法、反距離插值法和Kriging 空間插值法等3 種插值方法。本次模擬采用Kriging 空間插值。劉光孟等[11]在幾種插值方法的比較研究中，選取疏密差異較大的均勻和不均勻離散點數據進行試驗，得出在少量的不均勻分布數據中，Kriging 空間插值法具有較好的精度。因此，基于本研究區域有限的地表數字高程資料，選擇Kriging 空間插值法計算得到模擬區域地表高程模型的效果最優。在垂向上根據含水組特點劃分為3 層，通過Kriging 空間插值對數據進行初步處理后，獲得區域內第Ⅰ含水組底板高程（38.6～118.5 m）、第Ⅱ含水組底板高程（?36.2～83.7 m）及第Ⅲ含水組底板高程（?96.3～71.6 m）。其中，第Ⅰ含水組層厚為10～49 m、平均層厚19.3 m；第Ⅱ含水組層厚為10～76 m、平均層厚43.4 m；第Ⅲ含水組層厚為10～90 m，平均層厚53.6 m。

在黃壁莊水庫西側（山區與平原銜接處）和東北側（滹沱河沖洪積扇與磁河沖洪積扇交接處）分別存在一條水文地質分區界線，將這兩條分區界線作為西側和東北側研究區邊界。西側邊界處在山區與平原銜接處，作為弱透水的定流量邊界。西北部邊界設在黃壁莊水庫大壩位置，作為定水位邊界。由歷年地下水流場圖得知：研究區的南部邊界和北部邊界等水位線形狀變化不大，且邊界與等水位線基本垂直，因此模型的南部邊界和北部邊界可作為零流量邊界處理。研究區東部邊界設在地下分水嶺處，以藁城區崗上鎮至欒城區城關鎮為界，該處水位變化穩定，作為定水位邊界。

2.3 數學方程

結合石家莊市區的水文地質條件，將研究區模型概化為非均質非穩定三維流地下水系統。

2.4 計算參數確定

滹沱河地下水源地地層多為砂、砂卵礫石層，結構松散，孔隙度大，一般厚40～50 m，隨著地下水位逐年下降，目前水源地含水層厚20～30 m，西厚東薄；區域內富水性及導水性良好，滲透系數K 為300～400 m/d，導水系數T 一般在15 000 m2/d 左右，為強富水含水砂層。由于各含水層的水文地質參數在平面和空間上存在區域差異性，因此本文結合前人水文地質研究成果[12]進行參數率定，經過模型調試和識別，13 個水文地質參數分區參數值如表1 所示。

表1 研究區水文地質參數分區Tab. 1 Division of hydrogeological parameters in the study area

2.5 模型識別及結果

根據研究區的地下水位觀測資料，去除觀測序列不完整或者波動異常的水位觀測數據，選用永安村、省二院、省委黨校、西兆通村、凌透村、宋營村6 個水位監測井2016 年1 月至2017 年12 月的觀測數據進行模型率定識別，率定期內以30 d 為1 個時間步長，選取2017 年1—12 月作為模型的驗證期。本次模擬參數的率定采用試估-校正法，結合實際的水文地質條件和前人的研究成果，不斷對模型進行約束性的調試。采用置信區間、均方根、相關系數等進行誤差分析，計算值和觀測值的擬合效果見圖1。

圖1 監測井計算值與實測值的擬合Fig. 1 Fitting of calculated value and measured value of monitoring wells

從圖1 可見，監測井計算值和觀測值的最大擬合誤差是凌透村的監測點位0.934 m，最小擬合誤差是省委黨校的監測點位0.098 m，監測的誤差范圍為0.098～0.934 m，其中擬合誤差絕對值小于0.5 m 的監測井數占總監測井的83%，所有監測點位在滿足置信度為95%時的相關系數為0.999。由以上分析可得，模型的整體擬合效果頗為理想，滿足校正標準。

3 模擬結果分析

從計算流場的結果看，總體流向從西北向東南，最大流動速度小于1.0 mm/s。在滹沱河傍河地下水井群處，地下水位高程介于45～55 m，如圖2所示。

圖2 水源地區域點位及地下水流場（單位：m）Fig. 2 Point location and groundwater flow field in water source area (unit: m)

根據《飲用水水源地保護區劃分技術規范》[13]（HJ 338—2018）可知，地下飲用水源地一級保護區是以取水井為中心，溶質質點運移100 d 的距離所圈定的范圍。二級保護區是在一級保護區外，溶質質點運移1 000 d 的距離所圈定的范圍。準保護區是包括水源地的補給區和徑流區。在本研究中，溶質質點運移是利用MODPATH 粒子示蹤程序來進行計算可視化反向示蹤粒子的運移流線。

工況1：水源地的17 口水井單獨設置示蹤粒子，計算其運移距離。經過計算，示蹤粒子在第100 d 時運移距離為0.10～0.83 km，平均長度為0.54 km，其運移流線呈現長方喇叭狀；示蹤粒子在第1 000 d 時運移距離為3.0～7.4 km，平均長度為6.1 km，其運移流線同樣呈現長方喇叭狀，如圖3（a）所示。

圖3 不同工況下反向示蹤粒子的運移流線和保護區示意圖Fig. 3 Flow line of reverse tracer particles and schematic diagram of protected area under different working conditions

工況2：將水源地的17 口抽水井概化成1 個水井群，將反向示蹤粒子設置在井群的邊界處（井群邊界即應急后備水源地的邊界位置），計算其運移距離。經過計算，示蹤粒子在第100 d 時運移距離為0.12～0.83 km，平均長度為0.49 km；示蹤粒子在第1 000 d 時運移距離為3.1～7.4 km，平均長度為5.6 km，如圖3（b）所示。

由反向示蹤粒子流線圖可知：抽水井群的西北邊界100 和1 000 d 運移流線的長度相對于東南邊界較長。經分析可知：粒子軌跡方向與地下水位落差方向基本平行，抽水井群東南邊界及該區域的水井受到井群內部取水造成的井群內地下水位坡度平緩的影響，其粒子運移距離相對地下水井群西北邊界及該區域水井的粒子運移距離較短。

此外，通過入滲場補給水源時，計算過程同工況1，各單井設置反向示蹤粒子，計算其運移距離。滹沱河在南水北調干渠與滹沱河相交斷面上游建設長約5.5 km、寬約0.6 km 的入滲場，以補給地下水，入滲場的補給水量為2×108m3/a。滹沱河傍河擬建水源井分布在石家莊地下水庫范圍內，假設區域自黃壁莊水庫方向沒有地下水量補充，而是由滹沱河內入滲場作為主要的地下水補給來源。在此情景下滹沱河傍河地下水井的粒子運移流線將發生改變，經過反向示蹤粒子法計算，各擬建抽水井1 000 d 運移距離為2.4～3.5 km，平均長度為3.0 km。在滹沱河入滲場補給水量為2×108m3/a 的影響下，地下水井群局部地下水流場及水力坡度將會發生劇烈變化，此種情境下17 口擬建抽水井的補水來源主要是滹沱河入滲場，相應的反向示蹤粒子的流線方向也將變更為入滲場水源補給水源地的方向。此時，抽水井群1 000 d 的粒子流線范圍在地下水井群南部邊界與滹沱河入滲場北部邊界之間。同理，反向示蹤粒子在100 d 的運移流線也是在入滲場范圍內，如圖3（c）所示。

4 地下水保護區范圍的確定

在將地下水井按照設定的不同工況條件概化時，通過反向示蹤粒子可得到污染物的遷移距離，見表2。

表2 不同工況下的模擬計算結果Tab. 2 Simulation results under three working conditions

在工況1 中，當抽水井作為基準點單獨設置示蹤粒子進行計算時，100 d 的粒子運移流線長度范圍為0.10～0.83 km。又因各抽水井彼此相距0.80 km，部分0.10 km 的跡線范圍可能導致相鄰水井的一級保護區之間會有閑置區域。如果閑置區域存在人為活動，那么會增加水源地受到污染的風險，沒有起到相應的保護作用。在工況2 中，由于示蹤粒子的運移是隨著地下水流場的流動而進行的，運移流線的不規則分布和長度的差異可以直觀地表現出抽水時水源地對附近地下水流場的影響；同時，示蹤粒子反向運移的軌跡幾何形狀較工況1 更加規則，且遷移距離相對較短，能夠有效保證土地資源的利用。在入滲場補給水源地的條件下，滹沱河傍河地下水井的粒子運移流線將發生改變，100 和1 000 d 的運移流線都在入滲場范圍內。綜上所述，選取工況2 的模擬情景，并將特殊補給情況下的入滲場區作為一級保護區的補充范圍，可以做到科學合理劃分保護區，滿足水源保護的要求。

因此，根據工況2 模擬計算結果，將水源地一級保護區的邊界范圍設為：西側上游從水源地邊界外推0.49 km，東側下游外推0.22 km。北側外推到河南岸大堤外坡腳，南側從邊界外推0.12 km；同時，考慮特殊補水情況，故入滲場整個區域也要作為一級保護區的范圍。污染物在1 000 d 內運移到水源地的最長距離為7.4 km，而滹沱河地下水源地上游邊界距離黃壁莊水庫大壩只有15.9 km，因此，二級保護區從一級保護區上游邊界直接延伸到黃壁莊水庫壩下。兩側及下游二級保護區范圍主要參考數值模擬計算結果確定，上游外推到黃壁水庫壩下，下游以一級保護區邊界外推1.0 km；北部邊界考慮河堤的物理邊界作用，從河流北岸外推50 m；南部邊界從井群邊外推3.1 km。水源地保護區示意圖如圖3（d）所示。

5 結 語

以石家莊市滹沱河應急后備水源地為例，采用數值模擬法設置了2 種工況，對水源地污染物運移軌跡進行模擬計算，進而劃分一級、二級保護區。主要結論如下：

（1）利用MODFLOW 和MODPATH 等程序，建立了研究區的水文地質模型，模擬出地下水流場，每口抽水井單獨設置反向示蹤粒子和17 口水井概化成1 個水井群兩種不同的工況來計算反向示蹤粒子的流線，通過比較2 種工況下100 和1 000 d 粒子反向示蹤運移流線可知，將水源地概化成井群時反向示蹤粒子運移軌跡更加精準，運移距離相對較短，因此選擇工況2 并考慮入滲條件，更利于水源地的環境保護。

（2）考慮入滲條件，反向示蹤粒子在1 000 d 的運移距離完全在入滲場的場區范圍內，建議當地相關主管部門加強對入滲場的保護，采取相應措施制定水質安全突發事件的預案，以保證應急后備水源地安全，真正起到應急和后備的作用。