王 涵，劉 琦，張翼龍，王文中

(1.同濟大學(xué)土木工程學(xué)院地下建筑與工程系，上海 200092；2.中國地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所，河北 石家莊 050061)

地下飲用水水源的水質(zhì)直接影響著人類的健康。為了能讓地下飲用水源供水井遠離污染，建立地下飲用水源保護區(qū)，一方面可以通過一定的距離和時間，讓污染物在到達供水井前能夠被各種介質(zhì)所攔截凈化，使其濃度降低到目標(biāo)含量；另一方面，也可以讓相關(guān)部門在地下水受污染后有一定的應(yīng)急處置時間，使地下飲用水源的水質(zhì)達到規(guī)定的標(biāo)準(zhǔn)，滿足飲用及其他生活用水要求，實現(xiàn)地下飲用水源地的長效發(fā)展。

劃分地下水保護區(qū)的方法有很多，早期地下水保護區(qū)劃分主要依據(jù)經(jīng)驗，隨著研究的不斷深入，地下飲用水源地保護區(qū)的劃分方法也逐漸衍生出經(jīng)驗公式法、分析法、圓柱法、水文地質(zhì)條件法、解析模型法和數(shù)值模擬法[1～3]。然而隨著經(jīng)濟的發(fā)展，城市化進程的加快，土地資源越發(fā)稀缺，愈加寶貴，所以水質(zhì)的好壞固然是重要的考慮因素，但同時也需要追求經(jīng)濟和社會價值的最大化，即應(yīng)在合理的范圍之內(nèi)讓保護區(qū)的范圍更小[4]。如何更加科學(xué)、經(jīng)濟地劃分地下飲用水源保護區(qū)成為了一個更加有現(xiàn)實意義的研究課題。

數(shù)值模擬法[5]適用于水文地質(zhì)條件復(fù)雜，邊界條件眾多的地區(qū)，例如含水層較多且非均質(zhì)的含水系統(tǒng)、非穩(wěn)定流和分布多個注水井和抽水井的含水系統(tǒng)等[6～8]。本文主要以內(nèi)蒙古呼和浩特市為例，運用數(shù)值模擬法計算得到地下飲用水源地一級和二級保護區(qū)，并根據(jù)現(xiàn)實情況提出合理的保護建議，以期為呼和浩特市地下水資源的健康長效發(fā)展提供參考。

1 研究區(qū)概況

本研究區(qū)為呼和浩特市城區(qū)，面積約2 286 km2，如圖1所示。呼和浩特市城區(qū)飲用水源主要為地下水，目前共148口主要取水井，供水能力30.8×04m3/d，周邊區(qū)域工業(yè)企業(yè)及混采井分布密集，易受污染影響。因此，科學(xué)合理地劃分地下飲用水源地保護區(qū)對保障呼和浩特市的供水安全，防止地下水污染具有重要的意義。

圖1 研究區(qū)地理位置Fig.1 Location of the study area

1.1 水文地質(zhì)條件

呼和浩特市城區(qū)地下飲用水源地位于呼和浩特市平原區(qū)，主要分布更新統(tǒng)和全新統(tǒng)地層。地下水主要賦存于沖洪積砂礫層中，富水性強。區(qū)域主要接受大氣降水和地下水徑流補給，地下水流向大致為北東—南西向。研究區(qū)含水層模型見圖2，研究區(qū)含水層大致分為北部大青山單一潛水含水層(約占研究區(qū)15%)和南部呼和浩特平原潛-承雙層結(jié)構(gòu)含水層(約占85%)。其中北部單一潛水含水層巖性以卵礫石、中粗砂為主，一般厚度為120～140 m。地下水位埋藏較深，大部分地段地下水位埋深大于60 m，水量中等，單位涌水量100～500 m3/(d·m)。地下水徑流條件較好，水力坡度一般在1‰～3‰。南部雙層結(jié)構(gòu)含水層的上部潛水含水層厚度10～27 m，以砂礫石為主，水量豐富或極豐富。該含水層徑流條件良好，地下水由東流向南西，水頭高度1 018～1 080 m，水力坡度為1‰～2‰，其中北部靠近大青山地區(qū)潛水含水層疏干。下部承壓含水層上覆有穩(wěn)定較厚的淤泥層，承壓含水層主要由湖濱相粗碎屑物組成，水量豐富，含水層厚度較大，是呼市主要地下水開采層。該含水層北東向埋藏深度較淺，一般埋深40～60 m；南西向埋藏深度逐漸增大，顆粒由粗變細(xì)，含水層頂板逐漸變厚至接近200 m。承壓含水層地下水由北東流向南西，水頭高度1 010～1 060 m，水力坡度為1‰～3‰。

圖2 研究區(qū)水文地質(zhì)結(jié)構(gòu)概念圖Fig.2 Conceptual diagram showing the hydrogeological structure of the study area

1.2 淤泥層分析

研究區(qū)南部大部分地區(qū)潛水含水層與承壓含水層之間分布有較厚的淤泥層，該層厚度分布如圖3所示。

圖3 研究區(qū)淤泥層厚度等值線圖Fig.3 Contour map of thickness of the silty layer in the study area

淤泥層的滲透性直接影響著兩層含水層之間的連通性。對此有學(xué)者進行了詳細(xì)分析，根據(jù)沉積環(huán)境和實驗分析，淤泥層厚度小于70 m的區(qū)域淤泥層含砂量較高，因此從淤泥層分布邊界到厚度70 m等值線之間區(qū)域淤泥層滲透性較好，易于發(fā)生潛水和承壓含水層之間的越流[9]。研究區(qū)內(nèi)大部分取水井分布于該區(qū)域，受污染風(fēng)險較大。

1.3 水源地開發(fā)利用現(xiàn)狀及規(guī)劃

截止2015年，呼和浩特市內(nèi)在用的水井共有9個公共服務(wù)業(yè)用水水源廠下轄的119眼抽水井以及29眼城市供水補壓井，均為承壓水抽水井，總供水能力30.8×104m3/d。本文將研究區(qū)內(nèi)主要抽水井進行分類描述見圖4及表1。由于呼鋼水廠和城市供水補壓井年抽水量較小，因此本文數(shù)值模擬部分只考慮其他110眼抽水井。

近年來，因部分開采井成井質(zhì)量低，止水效果差，或未采取止水措施，以及存在一些未處理過的廢棄舊井，造成承壓和潛水含水層之間連通，開采中存在混采問題，且因為地下水的超采，呼市上部潛水含水層水位較下部深層含水層高，模擬時需要注意受污染的上層地下水越流污染承壓水的問題。

圖4 研究區(qū)內(nèi)抽水井分布及承壓含水層邊界條件Fig.4 Distribution of the pumping wells and boundary conditions of the confined aquifer

名稱抽水井?dāng)?shù)量(眼)實際單井平均供水量/(m3·d-1)一水廠161 288二水廠212 862三水廠142 911四水廠172 750五水廠141 930六水廠91 966金川水廠11686如意白塔水廠8103呼鋼水廠9調(diào)節(jié)運行,實際供水量較低補壓井29調(diào)節(jié)運行,實際供水量較低

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

2.1 水文地質(zhì)概念模型

首先對呼和浩特市含水層的實際空間結(jié)構(gòu)、滲透性質(zhì)、補給排泄及邊界性質(zhì)等條件進行概化，以便于進行數(shù)值模擬。

呼和浩特市城區(qū)148口抽水井位于呼和浩特平原，該處北面和東面環(huán)山，其余地方地勢平坦，含水層結(jié)構(gòu)大致分為下部承壓含水層，中部厚層淤泥質(zhì)黏土層以及上部潛水含水層。該地區(qū)主要供水水源為下部承壓水，因此本次模擬主要含水層對象為下部承壓含水層。

模型橫向長41.0 km，縱向長34.6 km，分為潛水含水層及承壓含水層兩層，按照500 m×600 m網(wǎng)格剖分。因為北部單一含水層與南部雙層結(jié)構(gòu)的潛水和承壓水都有水力聯(lián)系，為它們提供補給來源。按照地層時代，將北部單一結(jié)構(gòu)含水層分為兩層，上更新統(tǒng)和全新統(tǒng)為上層，中更新統(tǒng)和下更新統(tǒng)為下層。其中下層與南部雙層結(jié)構(gòu)中的承壓含水層合并為深層含水層，上層與南部雙層結(jié)構(gòu)中的潛水含水層合并為潛水含水層。因北面潛水含水層部分疏干，數(shù)值模型假定淤泥層的儲水能力及水平滲透系數(shù)可以忽略不計，將淤泥層概化進潛水含水層，水文地質(zhì)參數(shù)即概化為兩者的等效值，由不同的垂向水力傳導(dǎo)系數(shù)反映弱透水層，即采用MODFLOW中的“準(zhǔn)三維”方法(圖2)。

模型潛水含水層頂板設(shè)置為潛水含水層的自由水面，高程為1 018～1 080 m；潛水含水層底板結(jié)合北部單一結(jié)構(gòu)含水層厚度，根據(jù)淤泥層底板高度確定，高程為827～1 112 m；承壓含水層底板根據(jù)第四系中更新統(tǒng)下段含水層厚度及地面標(biāo)高確定，結(jié)合單一結(jié)構(gòu)區(qū)基巖底板，高程為503～1 061 m。

依據(jù)2011年5月水位統(tǒng)測資料繪制模型初始地下水流場如圖4所示，根據(jù)流場特征和對含水層結(jié)構(gòu)的分析，將呼市工作區(qū)的北部、東部確定為定流量邊界條件，為流入邊界，模型的第一、二層地下水經(jīng)過此邊界接受山區(qū)側(cè)向補給；將呼市工作區(qū)南西部邊界確定為定水頭邊界，為流出邊界，模型的第一、二層地下水經(jīng)由此邊界流出工作區(qū)，見圖4。

2.2 數(shù)學(xué)模型

工作區(qū)水位和地下水補、徑、排特征綜合反映地下水流為非穩(wěn)定流特性，水文地質(zhì)參數(shù)體現(xiàn)出含水層的非均質(zhì)特征，本文使用可考慮非穩(wěn)定、非均質(zhì)地下水流系統(tǒng)及越流效應(yīng)的“準(zhǔn)三維”數(shù)學(xué)模型[9～10]：

潛水為：

(1)

承壓水為：

(2)

式中：h1——潛水含水層水位標(biāo)高/m；

b——潛水含水層厚度/m；

h2——承壓含水層水位標(biāo)高/m；

z(x，y)——隔水底板標(biāo)高/m；

M＇——弱透水層厚度/m；

Kx、Ky、Kn——x、y及邊界面法線方向上的滲透系數(shù)/(m·d-1)；

Tx、Ty、Tn——x、y及邊界面法線方向上的導(dǎo)水系數(shù)/(m2·d-1)；

K——淤泥層的垂向滲透系數(shù)/(m·d-1)；

μ——潛水含水層的重力給水度；

S——承壓含水層儲水系數(shù)；

ε——含水層的源匯項/(m·d-1)；

p——蒸發(fā)和降水入滲強度/(m·d-1)；

Ω——滲流區(qū)域；

Γ0——滲流區(qū)域的上邊界；

Γ1——滲流區(qū)域的二類邊界；

q(x，y，t)——二類邊界單寬流量/(m3·d-1)，隔水邊界為0，流入為正。

2.3 數(shù)值模擬

因作者掌握的承壓含水層靜態(tài)水頭資料為2011年5月31日測量數(shù)據(jù)，且2015年動態(tài)水頭資料最為完善可靠，因此本次數(shù)值模擬模型的模擬期為2011年5月31日—2015年12月31日，共55個月。將整個模擬期按照對應(yīng)的自然月劃分為55個應(yīng)力期，所有外部源匯項的強度在每個應(yīng)力期中不變。模型的驗證期為2015年1月—12月，共12個月，以一個月作為一個應(yīng)力期。采用MODFLOW中的PCG算法進行計算[11]，計算得到的年均水均衡狀況，見表2，研究區(qū)內(nèi)潛水含水層向下部承壓含水層的越流總量約為0.21×108m3/a。

表2年均地下水均衡表

Table2Annualgroundwaterbudget

補排項補排量/(108 m3·a-1)降雨補給0.263補給項側(cè)向補給0.906小計1.169側(cè)向排泄0.198排泄項地下水開采1.031小計1.229補排差-0.060

如圖5所示，模擬的地下水流場基本能夠反映地下水的流動規(guī)律和趨勢。研究區(qū)總體水流方向為由東北、東南流往西南，地下水降落漏斗位于城區(qū)，范圍與實際情況基本相符。在研究區(qū)西北部和南部，擬合相對差些。簡而言之，研究區(qū)的模擬流場與實際情況總體一致，擬合較好。

圖5 2015年12月31日承壓水頭計算值及水位觀測井Fig.5 Calculation of pressure head and groundwater level observations wells on December 31, 2015

擬合的水位觀測點分布情況見圖5。大多數(shù)水位觀測井觀測值與模擬值擬合良好，兩者之間相差不大，能較好地反映該點水位隨時間變化趨勢；研究區(qū)中部降落漏斗處和南部小部分觀測點水位擬合相對較差，但隨時間變化趨勢一致。誤差主要與獲取的地質(zhì)資料豐富程度、源匯項統(tǒng)計誤差、以及控制性觀測井分布情況等因素有關(guān)。此外，研究區(qū)的源匯項復(fù)雜，模擬范圍大，且存在很多概化問題，因此誤差不可避免。圖6～圖8列出了研究區(qū)的幾口主要觀測井水位擬合狀況。

圖6 OW-6觀測井承壓水位擬合效果圖Fig.6 Fitting of groundwater levels at the OW-6 observation well

圖7 OW-16觀測井承壓水位擬合效果圖Fig.7 Fitting of groundwater levels at the OW-16 observation well

圖8 OW-25觀測井承壓水位擬合效果圖Fig.8 Fitting of groundwater levels at the OW-25 observation well

3 地下水保護區(qū)劃分

根據(jù)《飲用水水源保護區(qū)劃分技術(shù)規(guī)范》(HJ 338—2018)[12]，地下飲用水源一級保護區(qū)是以取水井為中心，井口周圍溶質(zhì)運移100 d的距離所圈定的范圍；在一級保護區(qū)范圍以外，井口周圍粒子運移1 000 d的距離所圈定的范圍為二級保護區(qū)，徑流區(qū)和補給區(qū)為準(zhǔn)保護區(qū)。

以各大自來水廠的抽水井為中心環(huán)形設(shè)置反向示蹤粒子，利用Visual MODFLOW中的MODPATH功能計算出各口抽水井周圍示蹤粒子逆向運移100 d和1 000 d的流線，流線的最外層質(zhì)點連接起來得到的范圍分別就是一、二級保護區(qū)，見圖9。

3.1 一級保護區(qū)劃分

數(shù)值模擬法計算出的各單井一級保護區(qū)范圍均不同，難以羅列，但計算得出的一級保護區(qū)形狀均類似于圖10所示的喇叭形，因此采用Visual MODFLOW中MODPATH的計算，各孔隙承壓水水源井?dāng)?shù)值模擬法一級保護區(qū)范圍計算結(jié)果見表3。

3.2 二級保護區(qū)劃分

根據(jù)數(shù)值模擬法得到的示蹤粒子1 000 d流線，結(jié)合現(xiàn)實道路、河流等劃分確定了二級保護區(qū)范圍見圖11。

3.3 保護措施

針對呼和浩特市地下飲用水源保護區(qū)，提出相應(yīng)的保護建議：

(1)對于一級保護區(qū)，取水點周圍設(shè)立警示標(biāo)志，區(qū)域內(nèi)不得堆放垃圾、廢渣，嚴(yán)禁大量使用農(nóng)藥、化肥，并對區(qū)域內(nèi)存在的混采井進行安全封堵，一切工程施工需進行污染相關(guān)論證，消除一切可能導(dǎo)致地下水污染的因素。

圖9 示蹤粒子運移100 d和1 000 d模擬圖Fig.9 Traces of particle movement for 100 d and 1 000 d

圖10 二水廠水源井區(qū)劃示意圖及一級保護區(qū)范圍Fig.10 Sampling of the wellfields in the second waterworks and the scope of the primary protection areas

名稱地下水流向長軸/m短軸/m側(cè)軸/m一水廠-5082 32 52 二水廠-172180 22 85 三水廠102105 23 55 四水廠-72247 25 63 五水廠19132 22 43 六水廠-151115 27 53 金川水廠-28124 20 55 如意白塔水廠-170138 21 32

(2)對于二級保護區(qū)，需保證在污染物流入一級保護區(qū)之前有充分的時間發(fā)現(xiàn)并消除污染，以確保水源地的安全。區(qū)域內(nèi)合理開采利用地下水，停止所有混采井的使用，并建立水位水質(zhì)長期監(jiān)測機制，定期監(jiān)測水位、水質(zhì)。

(3)需嚴(yán)密論證潛水含水層和承壓含水層的連通性問題并排查一切連通兩層含水層的工程施工，嚴(yán)防承壓含水層受到污染。

圖11 二級保護區(qū)范圍Fig.11 Secondary protected area

4 結(jié)論

本文以內(nèi)蒙古呼和浩特市為例，采用數(shù)值模擬法對主要供水井的一、二級保護區(qū)范圍進行求解，得出以下結(jié)論：

(1)采用VisualMODFLOW數(shù)值模擬軟件，建立了研究區(qū)的地下水滲流場，并對各主要抽水井MODPATH進行粒子反向示蹤模擬，根據(jù)井口周圍示蹤粒子運移100 d和1 000 d的流線范圍劃分一、二級保護區(qū)。根據(jù)水源地的特點，從社會經(jīng)濟和工程技術(shù)等方面提出了合理的地下水保護建議及保護區(qū)管理措施。

(2)本文利用“準(zhǔn)三維”方法簡化大流域復(fù)雜含水層環(huán)境，實現(xiàn)了地下水流場較為準(zhǔn)確的模擬。但由于源匯項統(tǒng)計誤差、地層資料的精度以及控制點分布不均等因素，對于細(xì)部地區(qū)的地下水模擬仍有較大誤差。建議在進行一級保護區(qū)劃分時將研究區(qū)劃分為多個小流域分別進行數(shù)值模擬，結(jié)果將更加準(zhǔn)確可靠。