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

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

(1. 中央民族大學(xué) 生命與環(huán)境科學(xué)學(xué)院，北京 100081；2. 生態(tài)環(huán)境部環(huán)境規(guī)劃院，北京 100012)

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

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

1 研究區(qū)概況

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

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

2 數(shù)值模擬法劃分保護區(qū)

2.1 軟件模擬過程概要

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

2.2 地下水流動模型

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

MODFLOW 共有自然鄰點插值法、反距離插值法和Kriging 空間插值法等3 種插值方法。本次模擬采用Kriging 空間插值。劉光孟等[11]在幾種插值方法的比較研究中，選取疏密差異較大的均勻和不均勻離散點數(shù)據(jù)進(jìn)行試驗，得出在少量的不均勻分布數(shù)據(jù)中，Kriging 空間插值法具有較好的精度。因此，基于本研究區(qū)域有限的地表數(shù)字高程資料，選擇Kriging 空間插值法計算得到模擬區(qū)域地表高程模型的效果最優(yōu)。在垂向上根據(jù)含水組特點劃分為3 層，通過Kriging 空間插值對數(shù)據(jù)進(jìn)行初步處理后，獲得區(qū)域內(nèi)第Ⅰ含水組底板高程（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。

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

2.3 數(shù)學(xué)方程

結(jié)合石家莊市區(qū)的水文地質(zhì)條件，將研究區(qū)模型概化為非均質(zhì)非穩(wěn)定三維流地下水系統(tǒng)。

2.4 計算參數(shù)確定

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

表1 研究區(qū)水文地質(zhì)參數(shù)分區(qū)Tab. 1 Division of hydrogeological parameters in the study area

2.5 模型識別及結(jié)果

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

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

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

3 模擬結(jié)果分析

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

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

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

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

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

工況2：將水源地的17 口抽水井概化成1 個水井群，將反向示蹤粒子設(shè)置在井群的邊界處（井群邊界即應(yīng)急后備水源地的邊界位置），計算其運移距離。經(jīng)過計算，示蹤粒子在第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 運移流線的長度相對于東南邊界較長。經(jīng)分析可知：粒子軌跡方向與地下水位落差方向基本平行，抽水井群東南邊界及該區(qū)域的水井受到井群內(nèi)部取水造成的井群內(nèi)地下水位坡度平緩的影響，其粒子運移距離相對地下水井群西北邊界及該區(qū)域水井的粒子運移距離較短。

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

4 地下水保護區(qū)范圍的確定

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

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

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

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

5 結(jié) 語

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

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

（2）考慮入滲條件，反向示蹤粒子在1 000 d 的運移距離完全在入滲場的場區(qū)范圍內(nèi)，建議當(dāng)?shù)叵嚓P(guān)主管部門加強對入滲場的保護，采取相應(yīng)措施制定水質(zhì)安全突發(fā)事件的預(yù)案，以保證應(yīng)急后備水源地安全，真正起到應(yīng)急和后備的作用。