泰萊盆地地下水傍河水源地準(zhǔn)保護(hù)區(qū)形態(tài)研究

楊海博，朱文峰，徐 建，趙 耘，袁金浩

（1.中國冶金地質(zhì)總局山東正元地質(zhì)勘查院，濟(jì)南250000；2.山東省第一地質(zhì)礦產(chǎn)勘查院，濟(jì)南250000）

在我國水資源的供給結(jié)構(gòu)中，地下水占據(jù)很大分量，地下水資源保護(hù)對保障用水安全具有重要意義。準(zhǔn)保護(hù)區(qū)等同于捕獲區(qū)的概念，已經(jīng)被廣泛應(yīng)用到水源地保護(hù)領(lǐng)域[1-3]。1979年，Bear 求解了均勻流場中井開采條件下區(qū)域水流被井抽取的分水線，這是早期準(zhǔn)保護(hù)區(qū)的概念基礎(chǔ)[4]。此后多井抽水、復(fù)雜區(qū)域邊界、相鄰河流或分水嶺等復(fù)雜條件下的準(zhǔn)保護(hù)區(qū)被勾勒研究[5-8]。隨著理論研究的成熟，準(zhǔn)保護(hù)的應(yīng)用研究逐漸成為國內(nèi)外熱點(diǎn)。2018年，李若怡利用粒子示蹤劃分了鄂爾多斯浩勒報(bào)吉水源地的準(zhǔn)保護(hù)區(qū)，將保護(hù)區(qū)擴(kuò)展到三維形態(tài)[9]。實(shí)際中，準(zhǔn)保護(hù)區(qū)具有更復(fù)雜的形態(tài)特征。

地下水水源地根據(jù)賦存條件可分為多種，其中傍河水源地最為普遍[10]，但傍河水源地準(zhǔn)保護(hù)區(qū)研究多停留在理論研究方面。1988年NeWsom 和Wilson 根據(jù)區(qū)域流流向的不同將排泄河流水源地準(zhǔn)保護(hù)區(qū)分成了三類：開采水來自側(cè)向徑流，開采水來自側(cè)向徑流同時(shí)河水補(bǔ)給地下水，開采水來自側(cè)向徑流和河流[11]。2011年Asadi-Aghbolaghi 等通過改變開采量也得到了以上3 種準(zhǔn)保護(hù)區(qū)[12]。此次研究將增加補(bǔ)給型河流，將準(zhǔn)保護(hù)區(qū)增加為五類：Ⅰ開采水來自側(cè)向徑流，Ⅱ開采水來自側(cè)向徑流，但受開采影響河水補(bǔ)給地下水，Ⅲ開采水來自側(cè)向徑流和河流補(bǔ)給，Ⅳ開采水來自河水，Ⅴ開采水來自側(cè)向徑流、河水以及河流另一側(cè)區(qū)域流。

目前傍河水源地準(zhǔn)保護(hù)區(qū)研究多限于理論方面，多河流研究更是處于空白。泰萊盆地地處黃河與京杭大運(yùn)河交匯處，水系發(fā)育，主要河流6 條，圍繞河流建有地下水水源地14個[13]。這獨(dú)特的地域條件為準(zhǔn)保護(hù)研究提供了基礎(chǔ)，同時(shí)對魯中地區(qū)用水安全具有重要意義。

1 研究區(qū)概況

泰萊盆地位于東經(jīng)116°20'52″～117°58'57″，北緯35°43'05″～36°34'15″，面積8 173 km2。其范圍橫跨泰安市肥城市、泰山區(qū)、新泰市的多個縣區(qū)，海拔高度范圍在37～1 545 m，最高處為“五岳獨(dú)尊”的泰山。研究區(qū)屬溫帶季風(fēng)大陸性氣候，四季分明。據(jù)多年資料統(tǒng)計(jì)，年平均氣溫在11～13 ℃之間；區(qū)內(nèi)多年平均降水量765.3 mm，日最大降雨量222.5 mm；年總蒸發(fā)量1 664.2 ～1 927.0 mm[14]。盆地內(nèi)水系發(fā)育，如圖1所示，河流主要有：大汶河、牟汶河、瀛汶河、柴汶河、康王河、匯河以及其他支流。水庫包括黃前水庫、光明水庫、上莊爐水庫和金斗山水庫等[15]。目前已建成傍河水源地14 座，均位于河流附近。由于水源地功能不同，日開采量之間存在較大差異。羊流宮里和石橫水源地日均開采量近7 萬m3，而寨里和埠陽莊水源地只有2 000～3 000 m3左右。不同水源地間開采深度和取水巖層也不相同[15]，詳細(xì)信息見表1。水源地開采會降低地下水位，導(dǎo)致地下水與河流的水力聯(lián)系更加復(fù)雜，進(jìn)而準(zhǔn)保護(hù)形態(tài)也更加豐富。

圖1 研究區(qū)水系與水源地分布圖Fig.1 Distribution map of water system and water source area

表1 傍河水源地信息表Tab.1 Information table of riverside water sources

研究區(qū)盆地地貌典型，在盆地南部、東部和北部有大范圍泰山群火成巖和變質(zhì)巖出露，風(fēng)化程度弱，且多處于高海拔地區(qū)。河谷地區(qū)主要為坡積、殘坡積、洪坡積形成的第四系沉積物，其成分主要是片麻巖、角閃巖、灰?guī)r碎屑與泥沙形成的混合物。第四系埋深最深處位于河谷中心，約50 m 左右，由河谷向外，埋深逐漸減小。在盆地西部和河谷下伏巖層為寒武—奧陶系地層，巖性以灰?guī)r為主，其中摻雜白云巖、頁巖等。局部地區(qū)發(fā)育石炭系頁巖和古近系砂巖。盆地內(nèi)局部發(fā)育巖溶和斷層地質(zhì)構(gòu)造，其可影響局部地下水流場，但整體水流場依舊受盆地地形控制[16]。

研究區(qū)內(nèi)地下水類型齊全，分為第四系孔隙水、碳酸鹽巖類巖溶裂隙水、碳酸鹽類巖溶水、變質(zhì)巖與巖漿巖類裂隙水，如圖2 和圖3所示。第四系孔隙水主要分布在河谷兩側(cè)，厚度50 m 左右，富水性強(qiáng)，在河谷附近埋深1～7 m，年變幅1～3 m。隨著距河流越遠(yuǎn)，水位埋深變化逐漸加大。含水層介質(zhì)為由變質(zhì)巖、碳酸鹽巖類、石炭巖、砂巖、頁巖組成的散巖類?？紫端饕艽髿饨涤辍⒑恿鳚B漏、巖溶越流補(bǔ)給。地下水流運(yùn)移主要受地形控制，整體呈現(xiàn)由地勢高向地勢低洼的河谷地區(qū)徑流，最后排泄到大汶河和其支流的下游。碳酸鹽巖類巖溶水和巖溶裂隙水主要位于盆地中南部地區(qū)。隨著深度增加，巖溶裂隙的發(fā)育程度逐漸減弱。在地表出露地區(qū)，水位埋深大，富水性強(qiáng)。受氣候和地質(zhì)條件影響，大氣降雨、河流滲漏、孔隙水入滲是巖溶水和巖溶裂隙水主要補(bǔ)給源。局部范圍內(nèi)地下水運(yùn)移受巖溶和裂隙發(fā)育的影響，大范圍內(nèi)依舊受地勢作用，呈現(xiàn)由北部、東部向南部徑流的趨勢。在南部地區(qū)巖溶水發(fā)生越流，向上層第四系孔隙水或河流排泄。除此之外人工開采也是巖溶裂隙水的排泄方式之一。變質(zhì)巖與巖漿巖類裂隙水主要分布在北部泰山巖群地區(qū)，由于巖層裂隙發(fā)育和風(fēng)化程度弱，富水性極差。地下水補(bǔ)給主要來自大氣降雨，但受巖層滲透性低的影響，地下水運(yùn)移能力差[17,18]。

圖2 水文地質(zhì)簡圖Fig.2 Hydrogeological sketch map

圖3 地層巖性剖面圖Fig.3 Stratum lithology profile

2 準(zhǔn)保護(hù)區(qū)劃分

2.1 概念模型與數(shù)學(xué)模型

此次研究將泰萊盆地視作一個獨(dú)立的水文地質(zhì)單元，三維模型如圖4(b)所示。研究區(qū)地下水受地形控制由高海拔向低海拔地區(qū)徑流，盆地邊緣多為火成巖山地，隔水性強(qiáng)，因此將北部、東部、南部設(shè)定為隔水邊界。盆地西南臨近東平湖，地勢最低，地下水從西南邊界流出，受地表水控制水位穩(wěn)定，所以設(shè)定為定水位排泄邊界。地下水上邊界相對復(fù)雜，源匯項(xiàng)主要有：大氣降雨、蒸發(fā)、河水滲漏、河流排泄和水源地開采。垂向上由于第四系松散層與基巖層水循環(huán)差異大，將模型分為上第四系潛水層和下基巖承壓含水層，在基巖出露的地區(qū)設(shè)定存在很薄的虛擬層[見圖4(a)]。

圖4 概念模型剖面圖及三維模型示意圖Fig.4 Cross-section diagram of conceptual model of groundwater flow system；schematic diagram of three-dimensional model

根據(jù)概念模型，可用數(shù)學(xué)模型描述穩(wěn)定地下水流系統(tǒng)如下，上部潛水含水層地下水流微分方程：

下部承壓含水層地下水流微分方程：

其中：

式中：h1、h2、h0、hR、hf分別為潛水含水層水位、承壓含水層水位、西南邊界水位、河水水位、地表高程；Kx、Ky分別為潛水含水層x和y方向滲透系數(shù)；T為承壓含水層導(dǎo)水系數(shù)；σ為導(dǎo)水系數(shù)；Q1i、Q2j分別為潛水水源地開采量和承壓水源地開采量；f、WR、r分別為河流分布函數(shù)、地下水與河床交換強(qiáng)度、河床滲漏系數(shù)；R、α、G分別為降雨入滲強(qiáng)度、降雨入滲系數(shù)、降雨強(qiáng)度；E、E0、Δ、Δ0分別為潛水含水層蒸發(fā)強(qiáng)度、水面蒸發(fā)強(qiáng)度、潛水位埋深、潛水極限蒸發(fā)埋深；n、S1、Γ1、S2、Γ2分別為邊界外展方向、第一層排泄邊界、第一層隔水邊界、第二層排泄邊界、第二層隔水邊界。

2.2 模型求解與校準(zhǔn)

本次研究使用GMS 軟件中的MODFLOW[19]模塊求解地下水穩(wěn)定流場，并在此基礎(chǔ)上利用MODPATH模塊采用粒子示蹤方法劃分水源地準(zhǔn)保護(hù)區(qū)。MODPATH 模塊可以記錄粒子在流場中的運(yùn)移路徑，并根據(jù)穩(wěn)定流場的模擬結(jié)果來劃定準(zhǔn)保護(hù)區(qū)。這是在半解析解粒子示蹤法的基礎(chǔ)上發(fā)展起來的。其原理是基于穩(wěn)定狀態(tài)下三維地下水流守恒的偏微分方程[20，21]：

式中：vx，vy和vz為地下水在坐標(biāo)軸3 個方向上的平均矢量流速；n為孔隙度；W為含水層單位體積內(nèi)補(bǔ)給或排泄地下水量的速率。方程（11）代表了極小體積內(nèi)的水量守恒。

數(shù)值模型求解使用PCG2 包，收斂準(zhǔn)設(shè)定為水頭變化精確到0.01，單元格直接流動殘差為0.001。模型各個區(qū)域的滲透系數(shù)利用242 個觀測孔進(jìn)行PEST 反演獲得，迭代收斂標(biāo)設(shè)定為0.001，具體參數(shù)如表2所示。

表2 模型分區(qū)參數(shù)表Tab.2 Model partition parameter table

MODFLOW 模型使用242個地下水觀測孔進(jìn)行校準(zhǔn)，其中第四系潛水觀測孔155 個，承壓水觀測孔87 個（見圖1）。觀測數(shù)據(jù)來源于《山東省地下水水源地調(diào)查評價(jià)報(bào)告》。模擬水位與實(shí)測水位相關(guān)性關(guān)系如圖5所示，潛水層標(biāo)準(zhǔn)誤差近7 m，承壓含水層標(biāo)準(zhǔn)誤差為10 m，但考慮到研究區(qū)內(nèi)包括海拔1 524 m 的泰山，海拔差異巨大，因而認(rèn)定校準(zhǔn)后的結(jié)果相對可靠，能夠用于準(zhǔn)保護(hù)區(qū)劃分。

圖5 模擬水位與實(shí)測水位對比圖Fig.5 Comparison of simulated water head and measured water head

3 結(jié)果與討論

通過粒子示蹤方法得到的傍河地下水水源地準(zhǔn)保護(hù)區(qū)分布如圖6所示，結(jié)合表1 信息可知，受開采量、河流、開采深度的影響，準(zhǔn)保護(hù)區(qū)的范圍各不相同。其中S13準(zhǔn)保護(hù)區(qū)范圍最大，其水源地具有開采量大、開采深度大、距離河流遠(yuǎn)的特點(diǎn)。S7、S8、S9、S10、S12、S14 準(zhǔn)保護(hù)區(qū)面積小，其都具有距離河流近的特點(diǎn)，同時(shí)S7、S8 開采量小，S9、S12 開采深度淺。由上可推測，開采量越大，開采深度越深，準(zhǔn)保護(hù)區(qū)面積越大。當(dāng)水源地距離河流近時(shí)，可以接受充足的水源補(bǔ)給，所以河流具有限制準(zhǔn)保護(hù)區(qū)范圍的能力。水源地S1、S3 形成的準(zhǔn)保護(hù)區(qū)不接觸河流，表明其開采的地下水均來自大氣降水補(bǔ)給，在水質(zhì)監(jiān)測時(shí)可以重點(diǎn)監(jiān)測保護(hù)區(qū)淺層水；S7、S8、S12、S14 水源地準(zhǔn)保護(hù)區(qū)只與河流的一側(cè)接觸，表明部分河段滲漏的河水經(jīng)過地下含水層被水源地開采，在水質(zhì)監(jiān)測時(shí)需對接觸段河水進(jìn)行監(jiān)測；S2、S4、S5、S6、S9、S10、S11、S13水源地準(zhǔn)保護(hù)區(qū)穿越了河流，能夠抽取河流另一側(cè)的地下水，接觸段河流滲漏的河水也會被水源地開采。通常情況下，河流往往會被認(rèn)定為水文地質(zhì)單元的邊界，但由于地下水開采，河流另一側(cè)的區(qū)域流也能被水源地抽取，所以河流另一側(cè)的地下水水質(zhì)也需要監(jiān)測。

圖6 傍河地下水水源地準(zhǔn)保護(hù)區(qū)分布圖Fig.6 Distribution map of quasi-protected zones of groundwater sources near rivers

傍河水源地開采使地下水與河流間的水力聯(lián)系發(fā)生變化，因此形成的準(zhǔn)保護(hù)區(qū)形態(tài)也更加復(fù)雜。在NeWsom和Wilson研究的基礎(chǔ)上，將傍河水源地分成了五類：Ⅰ開采水來自側(cè)向徑流，Ⅱ開采水來自側(cè)向徑流，但受開采影響河水補(bǔ)給地下水，Ⅲ開采水來自側(cè)向徑流和河流補(bǔ)給，Ⅳ開采水來自河水，Ⅴ開采水來自側(cè)向徑流、河水以及河流另一側(cè)區(qū)域流。分類在實(shí)際應(yīng)用時(shí)，會出現(xiàn)側(cè)向徑流的水來自河流的上游入滲，導(dǎo)致無法對保護(hù)區(qū)進(jìn)行分類的情況。這時(shí)需要根據(jù)當(dāng)?shù)貧夂驐l件和污染物運(yùn)移衰減程度對水源地周邊范圍進(jìn)行限制，當(dāng)滲漏河段處于限定范圍以外時(shí)，認(rèn)為開采的水來自側(cè)向徑流，當(dāng)滲漏河段位于限定范圍內(nèi)時(shí)，認(rèn)定抽取的水來自河流。如圖7所示，水源地S1、S3、S7、S8、S9 抽取的水來自附近區(qū)域的側(cè)向徑流，屬于第Ⅰ類；S5 抽取的水來自河流和周邊的側(cè)向徑流，屬于第Ⅲ類，同時(shí)在S5 附近出現(xiàn)了第Ⅱ類具有特征，這是因?yàn)橄薅ǚ秶笤斐傻?；S12、S14 水源地準(zhǔn)保護(hù)區(qū)只與河流的一側(cè)接觸但距離河流很近，因此劃歸為第Ⅳ類；S2、S4、S6、S10、S11、S13 水源地準(zhǔn)保護(hù)區(qū)穿越了河流，歸為第Ⅴ類。

圖7 水源地附近流線分布圖Fig.7 Streamline distribution near the water source

在Asadi-Aghbolaghi 的研究中，水源地準(zhǔn)保護(hù)區(qū)的類別隨著抽水量的增大而發(fā)生變化，如圖8(a)～圖8(c)所示。Q為抽水量，Q1、Q2為兩個關(guān)鍵值，當(dāng)Q＜Q1時(shí)，準(zhǔn)保護(hù)區(qū)為第Ⅰ類；當(dāng)Q1＜Q＜Q2時(shí)，準(zhǔn)保護(hù)區(qū)屬于第Ⅱ類；當(dāng)Q2＜Q時(shí)，準(zhǔn)保護(hù)區(qū)為第Ⅲ類?，F(xiàn)開采條件下S8 準(zhǔn)保護(hù)區(qū)為第Ⅰ類，在模型中增大S8 的開采量，流線變化如圖8(d)～圖8(f)所示。當(dāng)S8 開采量增大時(shí)會成為第Ⅱ類[圖8(d)]和第Ⅲ類[圖8(e)]，當(dāng)?shù)冖箢悳?zhǔn)保護(hù)區(qū)水源地開采量繼續(xù)增大時(shí)會變成第Ⅴ類[圖8(f)]準(zhǔn)保護(hù)區(qū)，這里存在一個關(guān)鍵值Q3。第Ⅳ類準(zhǔn)保護(hù)區(qū)出現(xiàn)于補(bǔ)給型河流附近，當(dāng)開采量增大到Q3＜Q時(shí)，第Ⅳ類會和第Ⅲ類一樣轉(zhuǎn)變?yōu)榈冖躅悳?zhǔn)保護(hù)區(qū)。為了分析關(guān)鍵值Q3的影響因素，在現(xiàn)有參數(shù)基礎(chǔ)上簡化模型，進(jìn)行敏感性分析。如圖9所示，Q3隨河流滲漏強(qiáng)度增大而增大，當(dāng)入滲強(qiáng)度大于0.1 m3/d 時(shí)，增大幅度急劇上升，當(dāng)水源地距離河流越近時(shí)，Q3越小。同時(shí)Q3還受區(qū)域流水力梯度和開采深度的影響，區(qū)域流水力梯度與Q3成正相關(guān)，而開采深度與Q3呈負(fù)相關(guān)。Asadi-Aghbolaghi 和此次研究都將河流概化成一條線，但在實(shí)際中河流具有一定寬度，所以Q1、Q2、Q3會從一個關(guān)鍵值轉(zhuǎn)化成一個關(guān)鍵區(qū)間。泰萊盆地14個水源地保護(hù)區(qū)分類情況與信息總結(jié)在表3中。

圖8 S8水源地附近流線分布對比圖Fig.8 Comparison of streamline distribution near the S8 water source

圖9 關(guān)鍵值Q3敏感性分析圖Fig.9 Key value Q3 sensitivity analysis diagram

表3 水源地準(zhǔn)保護(hù)區(qū)信息表Tab.3 Information table of water source quasi-protected zones

4 結(jié)論

此次研究基于NeWsom 和Wilson 對準(zhǔn)保護(hù)區(qū)分區(qū)理論研究，增加河流類型，并劃分了五類準(zhǔn)保護(hù)區(qū)?；谔┤R盆地多天然水系和傍河水源地特征，通過粒子示蹤劃分各個水源地準(zhǔn)保護(hù)區(qū)形態(tài)，并驗(yàn)證了NeWsom和Wilson理論研究具有實(shí)際支撐。傍河水源地準(zhǔn)保護(hù)區(qū)面積受到開采量、開采深度和水源地到河流距離影響，最后根據(jù)水源地附近流線分析得出開采水的來源區(qū)域以及附近河水對水質(zhì)的影響。根據(jù)Asadi-Aghbolaghi 關(guān)于抽水井開采量改變準(zhǔn)保護(hù)類型的理論研究，結(jié)合泰萊盆地水源地實(shí)際劃分結(jié)果和對S8進(jìn)行增大開采量模擬，完善了水源地開采量變化對準(zhǔn)保護(hù)區(qū)類型的影響，同時(shí)對關(guān)鍵值Q3進(jìn)行敏感性分析，發(fā)現(xiàn)水源地到河流距離、河流滲漏強(qiáng)度以及區(qū)域流水力梯度與Q3成正相關(guān)，開采深度與Q3呈負(fù)相關(guān)。在驗(yàn)證前人理論研究的同時(shí)，結(jié)合泰萊盆地實(shí)際多水源地準(zhǔn)保護(hù)區(qū)特征，完善總結(jié)了研究內(nèi)容，為水源地污染防治提供科學(xué)依據(jù)。